Участки генома, кодирующие белки, составляют лишь 1 % всего человеческого генома, но именно в них сосредоточены 85 % мутаций, определяющих возникновнеие и развитие различных заболеваний [1]. В связи с этим экзомное секвенирование, а также секвенирование с использованием специально разработанных панелей для обогащения, фокусирующихся на тех участках экзонов, в которых могут быть обнаружены значимые мутации, нашли наибольшее применение в клинической практике [2].

Один из важнейших вопросов клинического использования экзомного секвенирования — выбор адекватного числа покрытий (количество раз, которое был отсеквенирован каждый нуклеотид; обычно обозначается как 10x, 20x, 50x и т. д.) [3]. Именно оно позволяет выявить возможные ошибки считывания нуклеотидов на машине и определить истинные позиции в геноме. При этом мы сталкиваемся с двумя разнонаправленными факторами, влияющими на определение количества покрытий. Первый фактор — это время и стоимость секвенирования, которые увеличиваются с ростом числа покрытий. Второй фактор — статистический: какое минимальное количество покрытий позволяет свести до уровня допустимой погрешности ошибку при выявлении мутаций, причем по этому фактору не существует единого мнения. С чем это связано?

Используя технологию коротких прочтений Illumina, Bentley и соавт. в 2008 г. определили, что почти все единичные нуклеотидные варианты (SNV) в гомозиготе обнаруживаются при покрытии 15x, тогда как для обнаружения SNV в гетерозиготе необходимо покрытие 33x [4]. Основываясь на полученных данных, в большинстве последующих работ, связанных с обнаружением SNV, авторы использовали значение покрытия 33x в качестве стандартного [5, 6]. В 2011 г. Ajay и соавт. опубликовали статью, в которой показали, что для определения 95 % SNV, а также коротких инсерций и делеций необходимо устанавливать покрытие 50x. Однако последующие эксперименты с использованием усовершенствованных реагентов и программного обеспечения для обработки данных позволили авторам получить такой же результат, снизив значение покрытия до 35x [7]. В 2014 г. вышла статья Fang и соавт., в которой было найдено, что для обнаружения 95 % вставок и делеций необходимо устанавливать покрытие 60x [8].

Разброс в числе покрытий, представленный выше, показывает, что говорить об универсальности этого значения в настоящее время становится все сложнее, поскольку количество прочтений одного и  ого же участка для обнаружения вариантов напрямую зависит от качества прочтения этого участка. На качество прочтения влияет не только технология секвенирования, используемые реагенты, но также и подготовка образца. Например, трудности амплификации GC-богатых участков при проведении полимеразной цепной реакции (ПЦР) приводят к ухудшению качества прочтения и, как следствие, к необходимости увеличения числа покрытий. В настоящее время существует химия для ПЦР, позволяющая улучшить качество реакции и тем самым дальнейшее качество прочтения при секвенировании. В 2013 г. Meynert и соавт. определили, что для обнаружения 95 % SNP требуется покрытие от 20x до 46x в зависимости от используемой химии [9]. Этими же авторами в 2014 г. было обнаружено, что для детекции 95 % SNP достаточной глубиной покрытия является 14 ридов [10]. Кроме того, Li и соавт. показали, что глубина покрытия также зависит и от количества индивидуальных образцов, используемых при секвенировании [11]. Так, при обнаружении мутации, встречающейся с частотой менее 0,2 %, секвенирование 3 000 индивидуальных образцов с покрытием 4x дает аналогичный результат, что и секвенирование менее 2 000 индивидуальных образцов при покрытии 30x. Таким образом, мы видим, что факторов больше, чем кажется на первый взгляд, и количество покрытий может быть эффективно уменьшено за счет ряда деталей при проведении исследования и исходя из его конкретных целей.

В представленной работе показано, что при использовании панели для обогащения Genotek01 минимальное покрытие, достаточное для корректного определения гетерозигот и SNV, соответствует 12x, при этом расхождение результатов секвенирования и результатов валидирования по методу Сэнгера составляет 0,5 %.

МАТЕРИАЛЫ И МЕТОДЫ

Выделение ДНК, приготовление и секвенирование ДНК-библиотек

ДНК выделялась из цельной венозной крови, полученной от пациентов, страдающих наследственными заболеваниями, с помощью набора QIAmp DNA Mini Kit (Qiagen, Германия). Контроль качества геномной ДНК проводили с помощью агарозного гель-электрофореза. Важным показателем было отсутствие деградации ДНК и загрязнения РНК. Измерение концентрации выделенной ДНК проводили на приборе Qubit 3.0 Fluorometer (Life Technologies, США). ДНК-библиотеки готовили с помощью набора NEBNext Ultra DNA Library Prep Kit for Illumina (New England Biolabs, США) с использованием адаптеров для секвенирования на платформе Illumina, согласно протоколу производителя. Двойное баркодирование выполнялось с помощью ПЦР с набором NEBNext Multiplex Oligos for Illumina (Dual Index Primers Set 1) той же фирмы. Контроль качества полученных библиотек фрагментов ДНК был проведен с помощью Agilent 2100 Bioanalyzer (Agilent Technologies, США). Для таргетного обогащения кодирующих регионов генома использовался набор MYbaits (MYcroarray, США). Секвенирование проводилось на геномном анализаторе HiSeq 2500 System (Illumina, США) методом парных прочтений длиной 100 нуклеотидов с использованием наборов HiSeq Rapid PE Cluster Kit v2 и HiSeq Rapid SBS Kit v2 (Illumina) по протоколу производителя.

Секвенирование по Сэнгеру

Для проведения процедуры мечения ампликонов флуоресцентыми метками использовали наборы BigDye Terminator v3.1 Cycle Sequencing Kit (Thermo Fisher Scientific, США) по протоколу производителя. Секвенирование по Сэнгеру проводилось на генетическом анализаторе ABI PRISM 3500 Genetic Analyzer (Applied Biosystems, США) по протоколу производителя.

Биоинформатический анализ

Полученные чтения выравнивали на референсный геном hg19 с помощью программы BWA. Дедупликацию ридов выполняли командой rmdup в SAMtools, коллинг (variant calling) осуществляли с помощью пакета GATK (Genome Analysis ToolKit). Было обнаружено 89 мутаций: 10 гомо- и гемизигот, 79 гетерозигот; 80 точечных мутаций (SNP) и 9 коротких инсерций и делеций (инделов). Были также определены генотипы в регионах в 200 нуклеотидов влево и вправо от мутации. Все позиции из регионов, проанализированных в ходе анализа секвенирования, были валидированы с использованием секвенирования по Сэнгеру, которое считается золотым стандартом для детекции коротких мутаций. Хроматограммы были обработаны единообразно. Коллинг мутаций из хроматограмм был осуществлен с помощью собственного ПО, разработанного «Генотек» на основе BioPython и пакетов R: sangerseqR, seqinr, Biostrings и Rsubread. Нами было проведено сравнение генотипов, полученных методом высокопроизводительного секвенирования (NGS), с результатами секвенирования по Сэнгеру и рассчитаны чувствительность и специфичность.

РЕЗУЛЬТАТЫ ИССЛЕДОВАНИЯ

Валидация мутаций секвенированием по Сэнгеру

Секвенированием по Сэнгеру не были подтверждены 15 из 89 обнаруженных с помощью коллинга мутаций, т. е. они или имели другой генотип, чем был определен секвенированием, или отсутствовали. При этом 8 из 15 неподтвержденных мутаций имели гетерозиготный генотип, в то время как валидация по Сэнгеру определила наличие гомозиготной мутации. Стоит отметить, что в данном случае гетерозигота по данным NGS была поддержана только одним ридом с референсным аллелем (на рис. 1 кластер мутаций в правом нижнем углу графика).

Симулирование различного покрытия

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

Оценка качества секвенирования возможна с помощью метрики Phred quality score (Q score), выдаваемой секвенатором для каждого нуклеотида [12]. Однако данное качество является лишь оценкой точности секвенирования и не дает надежных значений для имеющихся данных. Мы проверяли, совпадает ли каждый имеющийся рид, перекрывающий такую позицию, по генотипу с референсным, и если находили отличие, то считали это ошибкой коллинга.

Были проанализированы 372 443 нуклеотида. Из них 276 были нуклеотидами, отличающимися от референсного, а остальные совпадали с ним. Таким образом, частота ошибок составила 0,0741 %, что эквивалентно примерно Q31 по Phred quality score.

Для 69 позиций, в которых имелись подтвержденные гетерозиготные мутации, мы оценили распределение доли ридов, поддерживающих альтернативный аллель (рис. 2).

Используя полученное распределение ридов из гетерозиготы и частоту ошибок коллинга, была рассчитана частота возникновения комбинаций для различных покрытий от 2x до 50x для различного числа ридов от 0 до максимально возможного, поддерживающих альтернативный аллель. Полученные частоты применили для расчета частоты двух типов ошибок: определение гетерозиготной позиции как референсной гомозиготы; и определение гетерозиготной позиции как альтернативной гомозиготы при различных порогах отсечки референсных и альтернативных гомозигот.

Для отсечки референсных гомозигот варьировали фиксированное число ридов (от 2 до 10). Для отсечки альтернативных гомозигот варьировали долю ридов, поддерживающих альтернативный аллель, между 70, 80 и 90 % (таблица).

Установлено, что для коротких мутаций (SNP и инделов) точность метода достигает 99,7 %, чувствительность — 98 % при покрытии в точке 12x. Меньшее покрытие приводит к значительному снижению чувствительности в сигмоидном характере и поэтому не может быть рекомендовано.

Для планирования эксперимента в лаборатории важно знать, с каким средним покрытием нужно секвенировать образец, чтобы таргетные регионы были покрыты на 12x. Для расчета этой величины мы построили корреляцию между средним покрытием образца и долей таргетного региона, покрытого 12 ридами (рис. 3).

Было обнаружено, что для покрытия таргетных регионов не менее чем на 12x на 90 % и более необходимо минимум 40x среднее покрытие дедуплицированными ридами.

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Секвенированием по Сэнгеру не удалось подтвердить 15 из 89 мутаций, при этом 8 из 15 неподтвержденных мутаций имели гомозиготный, а не гетерозиготный генотип. Наличие данных мутаций определяется моделью ошибок GATK, программой, позволяющей получить набор вариантов в исследуемом геноме, которая во время коллинга по-разному трактует одиночные риды с нереференсным и референсным аллелями. Дело в том, что в программе GATK используется модель «уверенности в референсе» (reference confidence model) в сочетании с когортным анализом [13, 14]. Поэтому в случае получения одиночного рида, совпадающего с нереференсным аллелем, GATK считает варианты нуклеотидов в данном риде ошибкой секвенирования и не учитывает их при расчете генотипа. А в случае вариантов, найденных в одиночном риде, совпадающем с референсным аллелем, она считает ошибку маловероятной и выдает гетерозиготный (а не гомозиготный) генотип. Кроме того, большинство мутаций, которые не были подтверждены секвенированием по Сэнгеру, имели низкое покрытие (≤ 10x). Полученные результаты подтверждают, что для корректного коллинга мутаций необходимо секвенирование с достаточно большим покрытием, чтобы одиночные ошибки секвенирования не искажали результаты [15].

Необходимое покрытие в точке является вероятностной величиной и может быть рассчитано с достоверной точностью. Нами было показано, что частота ошибок в результатах, получаемых с использованием HiSeq 2500 System, соответствует погрешности прибора, и рассчитано минимально необходимое покрытие в точке для практического применения — 12x. Эта величина меньше в сравнении с результатами Bentley и соавт. [4]. Данное улучшение может быть вызвано меньшим числом ошибок за счет использования усовершенствованного прибора и новых реагентов для секвенирования. Более точные современные методы биоинформатики также позволяют увереннее фильтровать ошибки секвенирования без потери чувствительности.

ВЫВОДЫ

Нами было показано, что для достижения не менее чем 90 % покрытия таргетных регионов на 12x или более, необходимо получить количество данных, эквивалентное среднему покрытию таргета 40x после дедупликации. Данная величина может зависеть от используемого набора и протокола для обогащения, типа и длины ридов. К тому же в зависимости от протокола создания библиотек и получения исходного биоматериала обнаруживается разная степень дупликации при одинаковом объеме данных, что также должно быть учтено при расчете необходимого числа нуклеотидов на образец.