Часть 4. Инструменты интерпретации и анализ данных секвенирования в Excel

Часть 4. Инструменты интерпретации и анализ данных секвенирования в Excel

Интерпретация данных секвенирования – самый сложный этап анализа, который лучше доверить специалистам медицинской лаборатории. Поэтому дальнейшая информация дана скорее в ознакомительных целях.

Excel поможет нам в «домашних» условиях приблизить анализ к многофакторному. Фильтры Excel позволят комбинировать оценки различных инструментов интерпретации, маркируя значимые варианты. А расширенные автофильтры позволят использовать генетические панели для различных заболеваний и, тем самым, учитывать фенотип.

Ключевые инструменты интерпретации: прогноз повреждения гена snpEff, клиническая значимость от Clinvar, частота аллеля (AF), оценка консервативности участка.

Тем, кто хочет познакомиться с тактикой клинической интерпретации данных, я рекомендую прочитать эту статью. Также можно найти официальное руководство  по интерпретации на русском языке.

Суть статей мне видится в том, что каждый из инструментов оценки значимости варианта имеет свои слабые стороны. Поэтому их последовательное применение (один за другим) может привести к отбрасыванию потенциально патогенных вариантов. Но существует и противоположная проблема – слишком большое число вариантов с «неопределенной клинической значимостью» (VUS), которые остаются после оценки любым инструментом.

Чтобы разрешить это противоречие нужен многофакторный подход. Варианты при помощи различных инструментов клинической интерпретации расставляют по приоритетам. Приоритет – это вероятность повреждения вариантом функции гена (и быть причиной заболевания). При этом особое внимание необходимо уделять семейной истории и фенотипу.

Инструменты интерпретации

Для начала давайте рассмотрим основные инструменты интерпретации, которые доступны при аннотировании VCF-файла с программой snpEff. А также ограничения этих инструментов.

Прогноз повреждения гена snpEff

Первая аннотация, которую делает программа snpEff на основе собственной базы данных – это прогноз сохранности функции белка в зависимости от расположения варианта. В главе о поломках генов мы уже рассмотрели, как делается такой прогноз, поэтому не будем здесь повторяться.

Вы, наверное, обратили внимание, что в большей части полей ANN (или EFF) указывается более одной аннотации. Гены часто имеют несколько транскриптов (подробнее тут). Поэтому один и тот же вариант может по разному влиять на транскрипцию нескольких генов. Например, вариант может быть повреждающим для одного гена и нейтральным для другого гена. Поэтому snpEff сообщает о влиянии варианта на каждый транскрипт.

Но предсказанное повреждение не означает патогенность. Практика показывает, что компьютерные модели snpEff, предсказывающие повреждение гена даже со степенью HIGH и с потерей функции (LOF), не позволяют утверждать, что вариант был причиной заболевания. Связано это не только с рецессивностью большинства заболеваний. Но также и с тем, что оценки повреждений уточняются со временем.

Расположение и структура многих генов пока не до конца ясны. Экзонные координаты генов, стартовый и стоп-кодоны иногда пересматриваются. Более того, неизвестно даже приблизительное число генов, что хорошо видно по расхождениям в ключевых базах данных. База Ensembl содержит 26 998 кодирующих белок генов, а GenBank – 21 104. Вариант может лежать в кодирующем экзоне в модели гена одного базы данных, но находиться в интроне или даже межгенной области в другой модели.

Клиническая значимость от Clinvar

ClinVar – это общедоступная база данных клинических отчетов о взаимосвязи между SNP и заболеваниями. ClinVar – очень удобный для работы, но ненадежный инструмент.

Уровень уверенности в клинической значимости варианта зависит от числа и качества подтверждающих данных. Если данные противоречивы, то обычно это так и указывается. Эта информация собирается из опубликованных медицинских статей, триалов и исследований.

Из базы данных ClinVar будет извлечено более 10 параметров. Это уникальный номер варианта, название ассоциированных нарушений (одного или нескольких), клиническая значимость и прочее

Какая бывает клиническая значимость в БД Clinvar

Клиническая значимость включает такие варианты, как “pathogenic” и “likely pathogenic” – патогенный и вероятно патогенный, на которые нужно обращать особое внимание. Впрочем, большая часть вариантов будет промаркирована как “benign” и “likely benign” – доброкачественный вариант (слабого действия). Также возможны варианты: drug response – связан с переносимостью лекарств; association – вариант связан с предрасположенностями из базы GWAS; risk factor – не вызывающий болезнь, но увеличивающий риск; affects – не заболевание, а нарушение, например, непереносимость лактозы; protective – снижающий вероятность заболевания; uncertain significance – неопределенная значимость; conflicting – противоречивые данные; other – вариант, у которого не нашлось подходящего термина из перечисленных; not provided – вариант был упомянут в связи с нарушением, но не описано клиническое значение; ‘-‘ – вариант описан в сочетании с другим аллелем и значимость нужно уточнить.

Впрочем, даже пометка “pathogenic” при гомозиготном варианте не означает, что он приведет к заболеванию. Хотя бы потому, что в базе Clinvar есть множество «патогенных» вариантов с высокой частотой аллеля, которые вряд ли можно считать опасными. Их скорее можно отнести к предрасположенностям.

Частота аллеля (AF)

Один из наиболее мощных инструментов генетического анализа одновременно и один из наиболее простых. Понятно, что разрушительные мутации редко встречаются именно из-за своей разрушительности. Поэтому именно в них стоит искать причину при подозрении редких менделевских наследственных заболеваний.

Но также понятно, что редкий не обязательно значит патогенный, поскольку белок может сохранять свою функцию. Поэтому частота аллеля – очень важный, но несамостоятельный инструмент.

В процессе аннотирования в VCF-файл будет добавлено сразу 18 столбцов с префиксом AF. Это частоты аллелей из баз данных dbSNP, Clinvar, dbNSFP.

Наиболее полная база данных, которая содержит частоты аллелей – это gnomAD. Поэтому хорошей практикой было бы использовать ее для аннотирования VCF-файла. Однако огромный размер самой базы данных и громоздкие аннотации могут стать проблемами при работе на домашнем компьютере.

Консервативность участка

Еще один прогноз повреждений основан на оценке консервативности участка. 10 лет назад, когда еще не было накоплено достаточной статистики по частоте аллелей, этот инструмент был одним из главных.

Степень консервативности устанавливается при сравнении белков человека с подобными белками других живых организмов. Говоря научным языком, было устанавливается соответствие между гомологичными белками в филогенетических исследованиях. Наиболее важные части ДНК в ходе эволюции оказались наименее подвержены изменениям у широкого разнообразия животных. И если несовпадения в белке (и на соответствующем участке ДНК) допускаются чаще, то вариант менее консервативен и менее разрушителен.

При аннотировании используются несколько инструментов оценки консервативности. Основные из них  – SIFT (Sorts Intolerant From Tolerant) и Polyphen2. Оценка сохраняется в виде чисел, а в зависимости от них, также в виде буквенных обозначений: D (Damaging) – вероятно повреждающий вариант, P (Probably) – возможно повреждающий вариант и B (benign) – доброкачественный вариант (слабого действия).

Что означают конкретные оценки, можно посмотреть в этой таблице.

Естественно, у такого инструмента будут свои ограничения. Некоторые человеческие белки неконсервативны (то есть отличаются от других позвоночных), поэтому связанные с ними варианты не определяются как вредные. Тем не менее, повреждение в них может приводить к болезням. Как и наоборот, несоответствия далеко не всегда будут патогенными.

Аннотации dbSNP

97% вариантов отфильтрованного экзома нашлись в базе данных dbSNP, то есть были снабжены уникальными RefSNP номерами (rs). А многие оставшиеся без аннотаций записи оказались инделями или смешанными вариантами. Они тоже заслуживают внимания, если нарушают функцию белка или попадают в область подозрительных генов.

Записи были аннотированы из базы dbSNP по различным признакам в 46 столбцах.

Помимо самих RefSNP ссылок, на мой взгляд, могут быть полезны пометки о том, что вариант упоминается в каких-либо медицинских источниках.

Например, в медицинской библиотеке Pubmed (столбец PMC), в справочнике менделевских заболеваний OMIM (столбец OM), полногеномных исследований GWAS и PAGE (столбец MT) или реакции на лекарства PharmGKB (столбец TPA). Важным может оказаться поле MUT (мутация), то есть, низкочастотный вариант, который цитируется в авторитетных источниках. Остальные поля, которые во многом дублируют другие базы данных. Их описание можно найти здесь.

Также очень важно, что база dbSNP содержит многие частоты аллелей (столбцы CAF), которых нет в других базах данных, которые мы будем использовать.

Подготовка файла Excel

Мы получили тысячи записей, каждую из которых вряд ли захочется проверять вручную. Но файл extracted.txt можно открыть в Excel и использовать его инструменты фильтрации для анализа.

Если вы сразу хотите увидеть, как может выглядеть подготовленный для анализа файл с формулами Excel, то можете скачать этот пример. Естественно, в вашем файле будет не две, а тысячи строк.

Кроме текста, процесс подготовки файла для анализа я описал в этом видео.

При импорте файла не забудьте указать, что десятичные знаки разделены точками. После вставки выровняем ширину столбцов и добавим автофильтр.

После столбца HET (гетерозиготность) вставим 7 новых столбцов с именами Flag, AF max, AF av, Pred, RS link, G1 link, G2 link.

Заполнять их мы будем так:

1. Первый столбецмы назовем Flag. Будем ставить ноль для потенциально важных вариантов, которые нужно будет рассмотреть. И туда же будем вносить цифровые пометки для просмотренных вариантов, которые можно будет отфильтровывать.

2. В столбцы AF max и AF av поместим информацию о максимальной и средней частотах аллелей, собранные из баз данных. Искать частоты мы будем при помощи формул для сгруппированных вместе 21-й записи из различных источников (от CAF[1], до dbNSFP_1000Gp3_AF[0]).

Формула для максимальной частоты аллеля:

=МАКС(DV2:EP2)

Формула для средней частоты аллеля

=СРЗНАЧ(DV2:EP2)

3. В столбец Pred сведем оценки консервативности и степени повреждения участка по различным моделям. Эти данные содержатся в 8 столбцах. Формально, мы учитываем записи, содержащие символы (D, P, A, H, M) в столбцах dbNSFP_MetaSVM_pred, dbNSFP_Polyphen2_HDIV_pred, dbNSFP_MutationTaster_pred, dbNSFP_MutationAssessor_pred, dbNSFP_Polyphen2_HVAR_pred, dbNSFP_SIFT_pred, dbNSFP_LRT_pred, dbNSFP_PROVEAN_pred.

=ИЛИ(ЕЧИСЛО((ПОИСКПОЗ("*D*";AG2:AN2;0)));ЕЧИСЛО((ПОИСКПОЗ("*P*";AG2:AN2;0)));ЕЧИСЛО((ПОИСКПОЗ("*A*";AG2:AN2;0)));ЕЧИСЛО((ПОИСКПОЗ("*H*";AG2:AN2;0)));ЕЧИСЛО((ПОИСКПОЗ("*M*";AG2:AN2;0))))

4. Далее сформируем столбец RS link с гиперссылками на сайт NCBI, чтобы быстро просматривать информацию о вариантах.

Второй (G1 link) и третий (G2 link)  столбцы сформируем с именами генов на GeneCards.

Гиперссылка будет составлена из пар сцепленных строк. Это адрес сайта и имя гена или RefSNP номер варианта в базе данных dbSNP, который мы возьмем из столбца RS.

Формула для rs записей выглядит так:

=ГИПЕРССЫЛКА(СЦЕПИТЬ("https://www.ncbi.nlm.nih.gov/search/all/?term=rs";BZ2);BZ2) 

А для генов так:

=ГИПЕРССЫЛКА(СЦЕПИТЬ("https://www.genecards.org/cgi-bin/carddisp.pl?gene=";Y2);Y2)
=ГИПЕРССЫЛКА(СЦЕПИТЬ("https://www.genecards.org/cgi-bin/carddisp.pl?gene=";Z2);Z2)

Таким же образом можно сформировать гиперссылки и на другие интересующие вас базы данных.

Все нулевые и некорректные строки можно удалить.

Анализ данных в Excel

1. Отфильтруем предельную частоту аллеля

Будем исходить из того, что мы ищем неизвестное редкое наследственное нарушение. Для серьезных менделевских заболеваний высокая частота аллелей в популяции может исключать патогенность варианта.

Какую именно следует устанавливать максимальную частоту аллеля, однозначных рекомендаций нет. Она рассчитывается исходя из частоты наиболее распространенного из интересующих наследственных заболеваний в популяции. Но в нескольких источниках [здесь и здесь] считается, что вариант будет скорее всего доброкачественным, если частота аллеля будет превышать 0.5% (0.005) для аутосомно-рецессивных нарушений. Для аутосомно-доминантных заболеваний условия жестче – частота аллеля не должна превышать 0.01% (0.0001).

По этим критериям мы можем отсортировать варианты в столбце AF av.

1) Устанавливаем на столбце HOM флаг true и отфильтровываем варианты AF av меньше 0,005.

2) Устанавливаем на столбце HOM флаг false и отфильтровываем варианты AF av меньше 0,0001.

В процессе заполняем Flag нулевыми значениями. Так мы будем знать, что вариант нуждается в проверке.

2. Проверим варианты по Clinvar

Прежде чем проводить детальный анализ можно провести «экспресс-тест». Суть его в проверке редких патогенных вариантов по клиническим источникам.

Во-первых, проверим нет ли в столбце MUT записей с пометкой true. Если таких записей нет, то проверим, не описаны ли найденные варианты в Clinvar.

Далее отфильтровываем подозрительные варианты по данным Clinvar (столбец CLNSIG) – pathogenic, likely pathogenic, uncertain significance, conflicting, not provided или ‘-‘.

Устанавливаем на столбец Flag фильтр 0.

Если в стобце CLNDISDB есть упоминание OMIM или Orphanet, то проверяем вариант по номеру записи в этих базах данных (просто ищем по номеру в поисковике). После проверки помечаем вариант в столбце Flag единицей.

Анализ патогенности варианта рассматривается в этой статье.

Понятно, что считать «патогенный» вариант подозрительным можно в том случае, если описание соответствует проявлениям нарушения. Особое внимание нужно обратить на раздел Molecular Genetics.

Ключевые моменты:

а) характерный тип наследования заболевания (рецессивный или доминантный);

б) возраст, в котором у человека проявляются признаки;

в) распространенность в популяции (регион, этническая группа);

г) пенетрантность и степень повреждения гена для доминантного наследования;

д) типы вариантов (SNP, повторы, крупные изменения), которые приводили к заболеванию.

Естественно, наиболее подозрительные варианты стоит перепроверить в других источниках и найти его напрямую в прочитанной последовательности ДНК в BAM-файле.

3. Проверяем варианты с потерей функции гена

Если экспресс-метод не дал ясных результатов, можно переходить к более широкому анализу. Рассмотрим редкие варианты, которые связаны с потерей функции гена или сильно влияют на его функцию.

Снимем фильтр с Clinvar и отфильтруем записи по столбцу LOF[*].PERC, с условием «не пустой». Таким критериям соответствует немного строк.

Аналогичным образом проверяем варианты с прогнозируемым повреждением функции гена степени HIGH.

Сначала снимем фильтр со столбца с потерей функции. Затем в столбце IMPACT отфильтровываем записи, содержащие слово HIGH.

Устанавливаем фильтр с нулевыми значениями на столбец Flag. То есть проверяем редкие варианты с сильным влиянием на функцию гена.

Сначала проверяем ссылки на странице NCBI (RS link), а если там пусто – проверяем заболевания, ассоциированные с дефектами гена (по ссылкам G1 link и G2 link).

4. Проверяем консервативность участка для варианта.

Проверяем консервативность участка и оценку повреждений на основе моделей.

Снимем все фильтры.

Установим флаг на пустые значения, чтобы не проверять варианты повторно.

Отфильтровываем значения «ИСТИНА» в столбце Pred. Это означает, что вариант затронул консервативный участок или вариант может привести к нарушению функции белка по какой-то модели.

Понимая, что высокочастотные варианты не могут быть патогенными, отфильтруем значения средней частоты аллеля более 0.5%, а также варианты с неизвестной частотой аллеля. Помечаем их нулями и затем проверяем.

В первую очередь стоит проверить варианты по базе Clinvar. Затем по базам NCBI и GeneCards.

5. Фильтруем варианты по панели генов.

Если и такой поиск не выявит ключевые причины заболевания, можно использовать клинические панели генов или панели SNP, связанные с фенотипом. Панели генов, ассоциированных с эпилепсией, аутизмом, ДЦП и многими другими нарушениями, несложно найти на сайтах генетических лабораторий. Например, здесь. Важно, чтобы имена генов в панели совпадали с именами в листе Excel.

Для фильтрации по панели можно воспользоваться расширенным автофильтром Excel. Для этого нужно будет составить столбец из интересующих нас генов на новом листе.

Также целесообразно рассмотреть по панелям сложные гетерозиготные варианты, то есть возможное сочетание двух различных рецессивных мутаций в одном гене.

6. Перебираем «бедную породу»

Проверяем, не осталось ли низкочастотных вариантов, с уровнем влияния на функцию гена MODERATE.

Кроме того, многие варианты не были аннотированы по частоте аллеля. Если у варианта есть rs-идентификатор, то вероятнее всего его можно найти онлайн по номеру из dbSNP. И, вероятно, там будет частота аллеля.

Также без аннотаций могут оказаться многие индели или смешанные варианты, которые ведут к потери функции гена, согласно оценке snpEff. Естественно, они могут иметь значение только в том случае, если будет прослеживаться связь нарушения ребенка с описанием функции гена.

7. Анализ предрасположенностей

Наконец, имеет смысл изучить распространенные варианты, с высокой частотой аллеля, но связанные с теми предрасположенностями и факторами риска, которые могли повлиять на развитие плода или ребенка в раннем возрасте.

Проверяем распространенные варианты по базам Clinvar и GWAS. Например, связанные с рисками беременности (если она протекала с осложнениями) или развития. Неврология может быть связана с предрасположенностями к иммунным и гормональным нарушениями, инфекциям или другим болезням.

Дисклеймер: Эта статья очень поверхностно раскрывает тему анализа NGS-данных. Она не может заменить учебники по генетике, мануалы и Pubmed. Ее цель – помочь на старте человеку без специального образования. Я не генетик (моя область – анализ данных). Поэтому прошу извинить меня за возможные неточности. Пожалуйста, напишите, если их заметите!

Содержание:

Часть 1. Коротко о секвенировании
1.1. Когда делают секвенирование?
1.2. Что ожидать от результатов?
1.3. Что лучше экзом или геном?
Часть 2. Немного теории: чтение ДНК, аллели, поломки генов
2.1. Что и как секвенируют?
2.2. Аллели и наследственные менделевские заболевания
2.3. Что может поломаться в гене?
Часть 3. Обработка файлов секвенирования от А до Я
3.1. Выравнивание данных в Galaxy: от FASTQ к VCF-файлу
3.2. Что мы будем делать с VCF-файлами?
3.3. Аннотирование VCF-файла c программой snpEff

Часть 4. Инструменты интерпретации и анализ данных секвенирования в Excel
4.1. Инструменты интерпретации
4.1.1. Прогноз повреждения гена snpEff
4.1.2. Клиническая значимость от Clinvar
4.1.3. Частота аллеля (AF)
4.1.4. Консервативность участка
4.1.5. Аннотации dbSNP

4.2. Подготовка файла Excel
4.3. Анализ данных в Excel