ОБЗОР
Количественный и безошибочный анализ данных массированного секвенирования с использованием молекулярного баркодирования
1 Российский национальный исследовательский медицинский университет имени Н. И. Пирогова, Москва, Россия
2 Группа структурной организации Т-клеточного иммунитета, отдел геномики адаптивного иммунитета, Федеральное государственное бюджетное учреждение науки Институт биоорганической химии им. академиков М.М. Шемякина и Ю.А. Овчинникова Российской академии наук, Москва
Для корреспонденции: Дмитрий Михайлович Чудаков
117997, Москва, ул. Островитянова, д. 1; ur.liam@mdvokaduhc
Финансирование: работа поддержана Российским научным фондом, грант №14-35-00105.
ОШИБКИ ПЦР И СЕКВЕНИРОВАНИЯ.
МЕТОДИКИ ИХ КОРРЕКЦИИ
Возможность надежно идентифицировать редкие варианты генов в сложных образцах делает высокопроизводительное секвенирование (High Throughput Sequencing, HTS) эффективным подходом в онкодиагностике [1], пренатальной диагностике [2], анализе неоднородности и изменчивости опухолей [3], в исследованиях бактериальных [4] и вирусных [5] инфекций и микробиомов [6], а также в эволюционных исследованиях [7] и исследованиях репертуаров иммунных рецепторов [8–13]. Однако одной из ключевых проблем в анализе данных HTS является накопление ошибок полимеразной цепной реакции (ПЦР) и секвенирования. Некоторые направления практического применения технологии, например ранняя онкодиагностика, требуют надежной детекции редких мутаций, присутствующих в образце в концентрации менее 1 % и даже в концентрации порядка 0,1 % [13–18]. В таких случаях редкие подварианты последовательностей генов часто невозможно отличить от подвариантов, являющихся результатом ошибок секвенирования и предшествующей ПЦР-амплификации. Самый простой способ устранения ошибок основан на том, что редко встречающиеся подварианты рассматриваются как ошибочные производные одной и той же часто встречающейся последовательности и по этому признаку исключаются из последующего анализа (frequency-based filtration) [19–21]. Фильтрация минорных подвариантов может опираться на более или менее достоверную модель накопления ошибок ПЦР, но в целом произвольно регулируется выбором соотношения часто и редко встречающихся последовательностей в массиве данных. Такой способ фильтрации хотя и устраняет большую часть накопленных ошибок, влечет за собой также потерю значимой части реального разнообразия редко представленных подвариантов нуклеотидных последовательностей в образце. В результате наиболее часто используемые методы HTS позволяют достоверно обнаружить мутации, присутствующие в образце только в значительной концентрации порядка 5 % [22–24] До недавнего времени не существовало метода, с помощью которого можно было бы проводить безошибочное глубокое секвенирование сложных библиотек генов при сохранении реального разнообразия гомологичных вариантов последовательностей. Прорыв в области количественного и безошибочного массированного секвенирования произошел с внедрением так называемого уникального молекулярного баркодирования [25, 26]. При этом подходе каждая исходная анализируемая молекула ДНК или кДНК маркируется уникальной нуклеотидной последовательностью, и эта информация затем используется в ходе программного анализа выходных данных массированного секвенирования (рисунок).
Опишем возможный вариант анализа молекулярно-баркодированных данных на примере разработанного нами подхода, получившего название Molecular Identifier Groups-based Error Correction (MIGEC) [10]. Он базируется на двухстадийном биоинформатическом анализе. Первая стадия относительно проста и основана на идее «безопасного секвенирования» (“Safe-Seqs”), предложенной в оригинальной работе проф. Vogelstein и соавторов [25]. Прочитанные последовательности, несущие один и тот же уникальный молекулярный идентификатор (баркод), объединяются в одну группу (кластер) — Molecular Identifier Group (MIG) (рисунокБ). Наличие идентичного молекулярного идентификатора показывает, что данные прочитанные последовательности были наработаны с одной и той же стартовой молекулы ДНК или кДНК.
Соответственно, кластеризация данных секвенирования по уникальному идентификатору позволяет установить исходную нуклеотидную последовательность по доминирующей последовательности в каждой группе (рисунокВ).
Таким образом можно устранить многочисленные ошибки, накопленные в процессе амплификации, и исправить практически все ошибочно установленные в процессе секвенирования нуклеотиды. Однако на практике применения такого прямолинейного алгоритма оказывается недостаточно, для того чтобы добиться безошибочного глубокого анализа библиотеки интереса. Дело в том, что ошибки ПЦР, произошедшие на ранних стадиях амплификации, могут в определенной доле случаев образовывать доминирующую нуклеотидную последовательность в MIG вследствие стохастической природы ПЦР-амплификации, и результирующий ошибочный вариант последовательности интерпретируется как реально существующее разнообразие (рисунок, Г). Такие события составляют значимый процент анализируемой информации при глубоком HTS-анализе библиотек и не могут быть безопасно (для реального разнообразия образца) отсечены фильтрацией редких подвариантов.
Любопытно, что такие ошибочные события, во всяком случае наиболее представленные из них, не устраняются также независимыми повторностями амплификации и секвенирования образца, так как частота определенных ошибок ПЦР в определенном контексте ДНК оказывается хорошо воспроизводима [10]. Для того чтобы идентифицировать и устранить подобные ошибочные варианты, мы ввели вторую стадию анализа данных (рисунок, Д). Она основывается на том, что высокочастотные (в каждом конкретном контексте ДНК) ошибки ПЦР носят повторяющийся характер, что позволяет отличить их от реального разнообразия. Такие ошибки «выдают» себя тем, что встречаются в качестве минорного подварианта в большом числе MIG и, соответственно, могут быть идентифицированы на основе относительной частоты встречаемости варианта последовательности в виде «мажора» или «минора» в MIG [10].
Двухстадийный алгоритм позволяет отфильтровать ошибочные варианты последовательностей с высокой точностью, сохраняя при этом естественное разнообразие библиотеки, и предоставляет возможность проводить глубокий безошибочный анализ сложных библиотек, как было показано нами на примере анализа данных секвенирования библиотек иммуноглобулинов и Т-клеточных рецепторов [10], а также на примере мультиплексного генетического анализа для онкодиагностики (наши неопубликованные данные). На примере анализа гомологичных контрольных последовательностей с различной представленностью мы показали, что анализ MIGEC устраняет практически всё искусственное разнообразие, которое содержится в данных секвенирования. Соотношение сигнал-шум, посчитанное как соотношение числа прочтений наиболее крупной контрольной последовательности к самому представленному ошибочному подварианту с одной и двумя нуклеотидными заменами, возросло с 1 000:1 и 20 000:1 для стандартно процессированных данных до 12 000:1 и 60 000:1 соответственно при обработке MIGEC. В то же время контрольные низкочастотные подварианты последовательности с заменами одного и двух нуклеотидов были сохранены.
КОЛИЧЕСТВЕННЫЙ АНАЛИЗ ДАННЫХ HTS
Независимо от выбранной технологии создания библиотек — на основе геномной ДНК либо на основе кДНК — невозможно обеспечить равную эффективность реакций на всех этапах пробоподготовки и секвенирования. Любая из стадий подготовки библиотеки генов к HTS-анализу (забор образца клеток, выделение ДНК или РНК, ПЦР-амплификация), как и само секвенирование, может приводить и непременно приводит к потерям и искажению информации о реальном количестве молекул и соотношении их подвариантов в исходном образце. Использование молекулярного баркодирования при анализе данных высокопроизводительного секвенирования позволяет контролировать число реально анализируемых стартовых молекул, успешно прошедших все стадии амплификации и секвенирования. В дальнейшем сравнительном анализе образцов также становится возможным оперировать не числом прочтений секвенирования, а количеством стартовых молекул. Таким образом, этот подход принципиально важен для понимания качества проведенного анализа и нормированного сравнения многих библиотек, в том числе полученных в разное время и разными лабораториями.
Так, в эксперименте, где с исходного образца, содержащего 1 000 гомологичных молекул, получают 106 прочтений секвенирования, число уникальных вариантов последовательности может варьироваться от 1 до 1 000 в зависимости от разнообразия этого образца. Однако детекция, например, 30 уникальных вариантов последовательностей среди полученного миллиона прочтений не дает однозначной информации о реальном составе образца. Действительно, эти 30 вариантов могут составлять разнообразие всей 1 000 исходных молекул или только лишь 30 молекул, успешно вошедших в амплификацию. В то же время обнаружение, например, 100 уникальных вариантов молекулярных баркодов позволяет с уверенностью говорить, что были проанализированы последовательности именно 100 стартовых молекул. Очевидно, что разрешающая способность такого эксперимента не может быть повышена за счет увеличения покрытия секвенирования, так как с увеличением числа прочтений число реально анализируемых молекул все равно останется равным 100. С помощью молекулярного баркодирования можно не только корректировать ошибки секвенирования и ПЦР, но и контролировать реальную узость «горлышка бутылки» для конкретного эксперимента, а также эффективно нормировать образцы для сравнительного анализа [9].
Становится возможным проводить точное сравнение двух и более библиотек генов, даже если они были получены с разного количества клеток и секвенированы с разной глубиной покрытия. Для этого, например, можно в ходе анализа использовать из каждого образца равное число случайно выбранных последовательностей, помеченных разными уникальными молекулярными баркодами. Поскольку каждая прочтенная последовательность с уникальным баркодом соответствует отдельной молекуле ДНК либо кДНК, такая нормализация данных секвенирования позволяет кардинально снизить уровень искажения количественной информации, накопленной в ходе амплификации и секвенирования кДНК-библиотеки.
В результате снижается разброс данных для независимо полученных реплик и повышается относительное детектируемое разнообразие, так как каждое новое прочтение соответствует новой молекуле ДНК или кДНК. Применение молекулярного баркодирования также позволяет снизить количественные искажения относительной представленности вариантов последовательностей внутри каждой конкретной библиотеки [12], так как в ходе анализа элиминируются искажения, связанные с неравной эффективностью работы праймеров при мультиплексной ПЦР, стохастической природой ПЦР, предпочтений ПЦР и секвенирования относительно различных матриц [9, 25–27].
ПОДВОДНЫЕ КАМНИ МОЛЕКУЛЯРНОГО БАРКОДИРОВАНИЯ
Несмотря на все преимущества применения молекулярного баркодирования, эта технология имеет ограничения, которые необходимо учитывать при проведении качественного глубокого анализа данных HTS. В частности, ПЦР-ошибки внутри самой последовательности уникального идентификатора могут вносить существенную погрешность в определение количества стартовых молекул кДНК или ДНК. По нашему опыту молекулярный баркод длиною в 12 случайных нуклеотидов после амплификации и секвенирования 104 раз обычно производит до 10–20 искусственных подвариантов баркода, и в совокупности они могут быть представлены 100–200 прочтениями секвенирования. В итоге после группировки прочтений по уникальным идентификаторам мы получаем от 11 до 21 стартовой молекулы кДНК, в то время как на самом деле она была только одна. Ещё большее число искусственных подвариантов молекулярных баркодов было обнаружено в модельной системе с последовательностями уникальных идентификаторов в 16 нуклеотидов [28]. Однако существует простой подход, позволя ющий достаточно эффективно отфильтровывать такие искусственные подварианты молекулярных баркодов. Как правило, при биоинформатическом анализе заметный процент данных составляют баркоды, прочитанные в ходе секвенирования один или два раза и по последовательности отличающиеся от более представленных «родительских» вариантов всего на 1 нуклеотид. Такая ситуация типична и воспроизводилась в серии проанализированных нами данных для различных экспериментов, в которых стартовый образец прошел 27–35 раундов амплификации. Подавляющее большинство вариантов баркодов с низким покрытием секвенирования представляет собой искусственное разнообразие, возникающее из-за ошибок ПЦР на поздних раундах амплификации. Устранять такое искусственное разнообразие библиотек можно с помощью фильтрации баркодов по количеству полученных прочтений. То есть в дальнейшем анализе используются последовательности уникальных молекулярных баркодов, секвенированных не менее определенного числа раз. Величина оптимального порога, выраженного в минимальном количестве прочтений на каждый уникальный баркод, может варьироваться в зависимости от размера стартовой библиотеки и достигнутого покрытия секвенирования. Также возможно проводить frequency-based фильтрацию по принципу наличия более представленных «родительских» вариантов молекулярных баркодов.
Другим важным моментом является естественное повторение случайно синтезированных уникальных идентификаторов («коллизий») [12, 28–29]. Так, при использовании в качестве молекулярных баркодов последовательностей длиною в 12 случайных нуклеотидов их теоретическое разнообразие составляет порядка 1,7Ч107 уникальных вариантов. Стоит отметить, что реальное наблюдаемое разнообразие всегда несколько меньше теоретического: по нашим оценкам, порядка 1,4Ч107 для последовательности из 12 случайных нуклеотидов. Такое снижение происходит вследствие неравномерного синтеза праймеров и, предположительно, наличия некоторых нуклеотидных предпочтений в ходе амплификации. Вероятность того, что в образце со случайно синтезированными, например, 10 000 молекулярными баркодами, найдутся два идентичных или различающихся лишь на 1 нуклеотид варианта достаточно мала. Однако очевидно, что при глубоком секвенировании, когда стартовое число молекул достигает нескольких миллионов, коллизии 12-тинуклеотидных баркодов случаются гораздо чаще. В таком случае бывает сложно отличить естественные коллизии от искусственного разнообразия молекулярных баркодов, вызванного накоплением в них ошибок ПЦР и секвенирования. Например, для 1 миллиона анализируемых молекул можно ожидать более 30 000 пар случайно синтезированных 12-тинуклеотидных вариантов молекулярных баркодов с идентичной последовательностью. Таким образом, при анализе данных глубокого секвенирования с использованием молекулярного баркодирования необходимо учитывать и теоретическое, и реальное разнообразие используемых уникальных идентификаторов. Предпочтительно использовать большее разнообразие вариантов (то есть большую длину случайно синтезируемой последовательности в составе используемого адаптера) при глубоком анализе больших библиотек.
МАССИРОВАНОЕ СЕКВЕНИРОВАНИЕ В ИССЛЕДОВАНИЯХ АДАПТИВНОГО ИММУНИТЕТА
В иммунологии молекулярное баркодирование открывает новые возможности для анализа сложных репертуаров вариабельных фрагментов антител и Т-клеточных рецепторов, разнообразие которых в организме каждого человека может составлять сотни миллионов вариантов.
Разнообразие Т-клеточных рецепторов и антител внутри одного образца также может быть чрезвычайно велико, а высокогомологичные варианты могут присутствовать в различных пропорциях и быть практически неотличимыми от накопленных ошибок при HTS-анализе. Молекулярное баркодирование делает возможным точный анализ и сопоставление информации о репертуарах антител и Т-клеточных рецепторов для различных субпопуляций лимфоцитов, органов и тканей здоровых и больных индивидуумов и позволяет надежно отслеживать изменения в разнообразии иммунных репертуаров и судьбу клональных популяций лимфоцитов с течением времени или после проведенной терапии. С использованием метода молекулярного баркодирования нами был проведен нормированный сравнительный анализ разнообразия репертуаров бета-цепей Т-клеточных рецепторов периферической крови здоровых людей различного возраста [9].
Было показано, что наблюдаемое разнообразие Т-клеточных рецепторов практически линейно убывает в течение жизни. Одновременно с возрастом происходит активное заполнение гомеостатического пространства высокопредставленными клонами Т-клеток при значительном падении процентного содержания единично представленных клонотипов. Учитывая относительно стабильное общее количество Т-лимфоцитов (оно лишь незначительно снижается с возрастом) разрастание субпопуляций эффекторных Т-лимфоцитов и Т-лимфоцитов памяти неизбежно приводит к снижению относительного количества наивных Т-лимфоцитов. Это приводит к снижению наблюдаемого и экстраполируемого разнообразия вариантов Т-клеточных рецепторов и вероятности инициации эффективного иммунного ответа против новых патогенов и онкологических заболеваний. Молекулярное баркодирование также повышает надежность анализа репертуаров рецепторов иммунных клеток для малых чисел исследуемых лимфоцитов. Это касается ситуаций, когда работа проводится с малыми популяциями сортированных или культивируемых лимфоцитов либо с образцами ткани, содержащими небольшое количество лимфоцитов.
При глубоком секвенировании библиотек с малым числом стартовых молекул искусственное разнообразие накопленных ошибок амплификации начинает существенно превалировать над реальным, а численные значения представленности вариантов последовательностей искажаются вследствие стохастической природы ПЦР. Использование молекулярного баркодирования эффективно устраняет как ложные подварианты последовательностей, так и накопленные количественные искажения, при этом полностью сохраняя нативную информацию о разнообразии репертуаров иммунных клеток в образце [12].
Исследование индивидуальных репертуаров антител и Т-клеточных рецепторов методом HTS в последнее время все чаще применяется в медицине: при количественной оценке минимальной остаточной болезни при терапии лимфопролиферативных заболеваний [30–34]; при отслеживании результатов аутологичной трансплантации гемопоэтических клеток крови [35–37]; при количественном определении наличия и клональности проникающих в солидную опухоль лимфоцитов (Tumor infiltrating lymphocytes, TILs) [38]; при отслеживании изменений периферического T-клеточного репертуара, вызванных иммунотерапией рака [39]; при поиске противораковых Т-клеточных рецепторов среди проникающих в солидную опухоль лимфоцитов [40–42] и др.
ЗАКЛЮЧЕНИЕ
Технология молекулярного баркодирования представляет собой мощный инструмент для нормированного безошибочного HTS-анализа. При достаточном уровне покрытия секвенирования использование уникальных молекулярных идентификаторов позволяет применять методы эффективной коррекции ошибок ПЦР и секвенирования, сохраняя при этом естественное разнообразие исследуемого образца. Последнее исключительно важно для анализа библиотек, включающих высоко гомологичные варианты последовательностей. Применение молекулярного баркодирования также делает возможным точный контроль глубины анализа данных высокопроизводительного секвенирования — в единицах анализируемых молекул исходного образца, а также проведение корректного сравнения данных по секвенированию сложных библиотек генов в условиях различного количества клеток/молекул на старте, различного качества ДНК, РНК или кДНК и различной глубины секвенирования. В биологии молекулярное баркодирование активно используется для решения разнообразных задач, таких как анализ геномов, транскриптомов [25, 26, 43] и микробиомов [6], исследование точности работы полимераз [25], оценка уровня ошибок в ходе транскрипции [44], синтез праймеров [25] или собственно секвенирование [29], а также анализ разнообразия репертуаров рецепторов иммунных клеток [8–11].В медицине распространение молекулярного бакродирования в HTS-анализе позволит достичь принципиально нового уровня надежности и чувствительности, что откроет новые возможности для ранней онкодиагностики и пренатальной диагностики по ДНК плазмы крови, для анализа гетерогенности и изменчивости опухолей, инфекционных агентов и микробиомов, достоверного анализа репертуаров рецепторов иммунных клеток, который все чаще находит применение в медицинской практике.