Исследование вариаций числа копий (CNV) в геноме у женщин на основании данных секвенирования внеклеточной ДНК плазмы крови

Размещено на http://www.allbest.ru/

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

ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ

«НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТЕТ

«ВЫСШАЯ ШКОЛА ЭКОНОМИКИ»

ФАКУЛЬТЕТ КОМПЬЮТЕРНЫХ НАУК

МАГИСТЕРСКАЯ ДИССЕРТАЦИЯ

Исследование вариаций числа копий (CNV) в геноме у женщин на основании данных секвенирования внеклеточной ДНК плазмы крови

по направлению подготовки 01.04.02 Прикладная математика и информатика образовательная программа «Анализ данных в биологии и медицине»

Кононкова Анна Дмитриевна

Научный руководитель к.б.н. Е. Шубина

Москва 2018



Содержание

Введение

Глава 1. Современные представления о вариации числа копий в геноме человека

1.1 Повторы в геноме как основные источники CNV. Механизмы формирования и классификация CNV

1.2 Методы детекции CNV

1.3 Методы детекции CNV на основе данных NGS

1.4 Программы, основанные на «глубине ридов» , для определения CNV при малом покрытии

1.5 Базы данных CNV

Глава 2. Материалы и методы

2.1 Образцы для поиска CNV

2.2 Тестирование алгоритмов поиска CNV

2.3 Сравнение результатов работы DWAQseq, cn.MOPS и QDNAseq на большом количестве образцов

2.4 Генерация случайных локусов для пермутационного анализа

2.5 Аннотация полученных результатов

Глава 3. Результаты

3.1 Тестирование различных инструментов для определения CNV на малом покрытии

3.2 Анализ CNV в российской популяции у женщин на основе данных неинвазивной пренатальной диагностики

3.3 Анализ результатов QDNAseq с учетом X - хромосомы

3.4 Поиск CNV с помощью программ DWAC-seq и cn.MOPS

3.5 Каталог всех детектированных CNV

Заключение

Литература



Введение

геном рид локус детектированный

Для генома человека характерна стабильная, эволюционно обусловленная структура; ее индивидуальные различия составляют в среднем менее 0.01%. Генетическим разнообразием внутри вида мы обязаны небольшим изменениям в геноме на уровне нуклеотидных замен, вставок, выпадений либо более крупномасштабным событиям, таким как вариации числа копий (CNV). Вариация числа копий в современных терминах - это изменение стандартного количества копий участка генома размером более 50 п.о. CNV могут приходиться как на гены, так и на некодирующие последовательности, что в свою очередь приводит к изменению дозы гена, меняет его экспрессию на регуляторном уровне или никак не проявляется. Некоторые СNV имеют клинические проявления, оказываются причиной патологии, тогда как другие не влияют на фенотип их носителя. Результаты нескольких исследований обнаруживают связь между CNV в определенных участках генома и заболеваниями ЦНС, предрасположенностью к психическим заболеваниям, аутизму. С внедрением в клиническую практику методов молекулярного кариотипирования количество данных по CNV значительно возросло. В существующих базах данных, например DGV, содержится полезная информация о локализации CNV, их частоте, затронутых генах, что позволяет сделать правильные выводы по результатам исследований индивидуальных генетических особенностей пациента. Однако в международных базах данных CNV российская популяция представлена лишь небольшими группами в составе крупномасштабных исследований [29, 49]. Сравнение результатов анализа CNV у представителя российской популяции с доступными на сегодняшний день данными может оказаться некорректным, так как частоты вариаций могут иметь отличия между популяциями. Например, показано, что результаты широкогеномных ассоциативных исследований (в особенности из-за паттернов неравновесного сцепления), меняются в зависимости от исследуемой популяции. [38] Более рационально было бы использовать в клинической практике в качестве референсных данные по CNV с учетом популяционных отличий, либо нужны качественные исследования, указывающие на их отсутствие. В настоящий момент запущен проект “Российские геномы”, отсеквенировано более 500 геномов, часть из которых уже оцифрована, однако создание базы данных CNV по-прежнему рассматривается в перспективе.

В лаборатории молекулярно-генетических методов ФГБУ НМИЦ АГП им. Кулакова в рамках проведения неинвазивного пренатального ДНК-скрининга (НИПС) анеуплоидий было отсеквенировано с ультранизким покрытием (менее 0.4х) порядка 1000 геномов беременных женщин.[48] Доля плодовой ДНК составляет порядка 4 - 20%. Оставшаяся доля ДНК приходится на материнскую и обычно никак не используется в дальнейшем. Существуют отдельные исследования влияния CNV у матери на получение ложноположительных результатов НИПС.[42] Также в литературе освещен клинический случай обнаружения в ходе НИПС хромосомных перестроек у матери вследствие множественной миеломы на ранней стадии.[51] Однако в основном информация, полученная в ходе НИПС, исчерпывается исследованием хромосомного статуса плода (возможна также детекция CNV размером более 10 Мб при достаточной доле плодовой ДНК -15 % [52]).В то же время к настоящему моменту появилось большое количество разнообразных биоинформатические инструментов, которые позволяют исследовать CNV на основе данных NGS. Например, в третьей фазе проекта 1000 геномов информация по структурным вариациям была получена именно этим способом.[43] Некоторые методы подходят для генома с низким покрытием. Таким образом, учитывая масштаб применения НИПС во многих странах, эти данные могли бы послужить источником ценной информации о CNV в популяции у женщин и быть использованными в качестве референсных для сравнения в клинике в частности и популяционных исследованиях в общем.

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

Основные задачи:

- проанализировать существующие биоинформатические инструменты для поиска CNV в геноме с экстремально низким покрытием

- определить параметры и критерии для отбора CNV для каталога из полученных данных

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



Глава 1. Современные представления о вариации числа копий в геноме человека

Первые упоминания о вариации числа копий восходят к началу 1900-х годов, когда с помощью микроскопа было обнаружено, что мужскому фенотипу отвечает Y хромосома, а женскому X. [1] В 1930х впервые была показана трисомия по 21 хромосоме, CNV можно было изучать на уровне анеуплоидий и они прежде всего были связаны с патологией. [2] Только в 2004 году было проведено масштабное систематическое исследование на предмет наличия CNV у здоровых индивидов. [3]

В настоящее время известно, что CNV - один из основных источников полиморфизма в геноме человека. Для большей части генов мы унаследовали по одной копии от каждого из родителей, однако для некоторых фрагментов генома количество копий может отличаться от стандартного диплоидного набора. Так, большая часть бета-дефензинов может присутствовать в геноме в количестве копий от 2 до 7. Другой яркий пример - полипептид B17 семейства УДФ-глюкуронилтрансферазы 2. Восточно-азиатские популяции являются носителями гомозиготной делеции, что приводит к нулевому количеству копий. [2] Также для высокочастотных CNV показано функциональное обогащение генами, связанными с ответом на лекарственные препараты, сенсорное восприятие. Ранние исследования показывали наличие корреляции между CNV и генной экспрессией, однако позже выяснилось, что для 20% случаев корреляция оказалась обратной. [4] Таким образом, влияние CNV на функционировние генома неоднородно, нужно также учитывать возможность дозовой компенсации.

В современных терминах CNV - вариабельный участок генома размером более 50 п.о., может быть унаследованы или образоваться de novo (темп образования 10-6 - 10-8, что значительно выше по сравнению с однонуклеотидными вариациями) [5]; CNV, присутствующие более чем в 1% популяции, называются полиморфизмом числа копий (CNP). Считается, что на долю CNV приходится 4.8-9.5% генома, в геноме насчитывают порядка 100 генов, отсутствие которых никак не влияет на фенотип. [6] Около 25000 полиморфных CNV обнаружено у здоровых индивидов, среди них порядка 1000 с длиной, превышающей 50 кб. CNV также существуют в геномах других млекопитающих. Сравнительный анализ геномов человека и шимпанзе показывает, что некоторые CNV находятся под отбором; характерные для человека участки обогащены генами, связанными с клеточной пролиферацией и ответом на воспаление, что указывает на адаптивное значение CNV. [7]. Количество CNV у одного человека варьируется в диапазоне от 5 до 50 и может быть разным по популяциям. [8] Вообще подобные различия могут быть использованы, например, для исследования миграционных путей в популяционной генетике. Так, в работе 2011 года Chen et al показали, что кластеризация по СNV неплохо соотносится с количеством популяций, исследуемых в работе. [9] Для крупных CNV такая стратификация не обнаружена. Крупные CNV (более 500 кб) редко бывают общими для популяции, но каждый индивид имеет в своем геноме одну и более. CNV длиной более 1 Мб составляют 1-2%. Для них характерно большее количество дупликаций по сравнению с делециями. Противоположная тенденция наблюдается для малых CNV. [10] Около 66%, 22% и 7% приходится на длину 100-250 кб, 250-500 кб и 500кб-1Мб.

1.1 Повторы в геноме как основные источники CNV. Механизмы формирования и классификация CNV

Все CNV условно делятся на 2 больших класса: нерекуррентные и рекуррентные. Точки разрывов и длина рекуррентных CNV схожи у разных индивидов; около 20-40% полиморфных и множество de novo патогенных CNV имеют рекуррентные границы. CNV локализуются в геноме не случайным образом, их позиции хорошо коррелируют с повторяющимися последовательностями в геноме и в меньшей степени с теломерными и центромерными участками. [16] Повторы в геноме также можно условно разделить на 2 категории - низкокопийные повторы (сегментные дупликации длиной более 1 кб) с высокой идентичностью (обычно для удобства рассматривают более 95%) и высококопийные, покрывающие около 44% генома. Наиболее распространенный вид высококопийных повторов - диспергированные повторы, значительную часть которых составляют ретротранспозоны, представленные короткими диспергированными повторами (SINEs), длинными диспергированными повторами (LINEs) и эндогенными ретровирусами человека. SINE - короткие последовательности ДНК (100-400 п.о. в длину), составляют около 11% генома. Самый распространенный SINE повтор - это Alu - элемент (около 300 п.о.), он занимает важное место в эволюции приматов. Самый частый среди LINE-повторов - L1-элемент (6-8 кб), покрывает около 17%. [12] Также среди повторяющихся структурных элементов генома встречаются варьирующие по числу тандемные повторы - VNTR - вид минисателлитных ДНК, которые содержат 10-15 нуклеотидов в повторяющейся последовательности, и собственные цепи (SC) размером 150-1000 п.н. [13] Сегментные дупликация - это не что иное, как зафиксировавшиеся внутри вида CNV; они могут содержать псевдогены, эндогенные ретровирусные последовательности, гены. [14]

Риc.1 Типы CNV по классификации точек разрывов. [15]

Рекуррентные CNV в общем случае образуются вследствие внутри- и межхромосомной неаллельной гомологичной рекомбинации в мейозе (NAHR), обычно обусловленной ошибочным выравниванием фланкирующих сегментных дупликаций длиной 10-100 кб (28% CNV) или транспозонов длиной 300-500 п.о. [16]

Уровень NAHR положительно коррелирует с длиной низкокопийных повторов, но обратно пропорционален расстоянию между ними. NAHR в 7 раз чаще встречается в больших CNV по сравнению с VNTR, которые характерны для малых. [17] Известно, что для неаллельной гомологичной рекомбинации нужен гомологичный участок длиной не менее 300 п.о. Относительно короткие Alu - элементы также часто встречаются в точках разрывов рекуррентных CNV, хотя в большей степени они способствуют возникновению сегментных дупликаций, особенно древних.[14] Кроме того показано, что 83% генома человека склонно к LINE-LINE рекомбинации, и наиболее частый элемент LINE - L1 - может также приводить к NAHR. [12] Механизм NAHR порождает больше делеций, чем дупликаций. Это обусловлено тем, что в этом случае способов образования делеций больше, чем дупликаций. [16]

Риc.2 Основные механизмы образования CNV. [18]

Риc. 3 Механизмы образование CNV вследствие репарации разрывов из-за остановки вилки репликации. [19]

Для нерекуррентных событий характерны сложные структурные перестройки. Несколько делеций могут чередоваться с дупликациями, трипликациями, инсерциями и стандартным количеством копий, соединенных между собой небольшими гомологичными участками. [20]

Нерекуррентные CNV не связаны с сегментными дупликациями, хотя часто встречаются вблизи участков с повышенной плотностью SD, что говорит скорее о сходстве механизма их появления в эволюционном контексте. Точки разрывов обычно относятся к участкам микрогомологии (2-15 п.о.), что мало для NAHR. [21] Положение нерекуррентных CNV хорошо коррелирует с VNTR. Источником возникновения может быть негомологичное или обусловленное микрогомологией спаривание концов (NHEJ и MMEJ соответственно). У однояйцевых близнецов находят различные CNV, т.е. их образование имеет место не только в мейозе. Действительно, некоторые исследования показывают, что формирование нерекуррентных CNV связано скорее с митозом. [16] Например, нерекуррентные CNV, не ассоциированные с сегментными дупликациями и связанные с интеллектуальными нарушениями, имеют отцовское происхождение. Короткие участки гомологии, особенно инвертированные повторы, могут приводить к образованию вторичных структур ДНК, остановке вилки репликации и переключению матрицы (механизм FoSTeS). [22] Обусловленная микрогомологией репликация после разрыва - MMBIR - также дает вклад в образование нерекуррентных CNV. FoSTeS и MMBIR могут быть причиной CNV длиной от 100 п.о. до нескольких мегабаз. [22] Эксперимены показывают, что ингибирование репликации индуцирует появление de novo CNV, которые схожи с нерекуррентными по размеру, распределению и структуре концов. [16]

Также наблюдается корреляция между временем репликации и образованием CNV: делеции связаны с участками поздней репликации, а дупликации с участками ранней. [23] Показано, что малая скорость репликативной вилки и частые остановки также могут быть связаны с CNV. Наблюдается тенденция в пространственной колокализации граничных точек CNV и их одновременной репликации. [16] Анализ Hi-C данных говорит о том, что 7,5% граничных точек CNV совпадают с границами топологически ассоциированных доменов (TAДов), а 30.7% лежат в непосредственной близости (120 кб) от них, что существенной больше по сравнению со случайно разбросанными точками разрывов. [24]

1.2 Методы детекции CNV

Существует несколько методов детекции CNV. Кариотипирование имеет ограниченное разрешение и определяет структурные перестройки размером не менее 5 мб. Методом флуоресцентной гибридизации in situ (FISH) можно опередлить CNV размером 100 кб - 1 мб. Меченные флюорофором/биотином специфические последовательности (зонды) гибридизуются с таргетными участками ДНК, после чего с помощью флуоресцентного микроскопа можно определить делеции, дупликации и другие структурные перестройки. FISH-метод позволяет исследовать клетки на любой стадии клеточного цикла, в пробе может присутствовать генетический материал одной клетки, а не усредненная по клеткам смесь ДНК. Недостатком является то, что анализ проводится только для фрагментов ДНК с заранее заготовленными зондами.

Рис. 4 Общая схема реализации FISH-метода и aCGH.

Матричная сравнительная геномная гибридизация - array CGH - относительно новая генетическая технология. Матричная CGH основана на сравнении генома пациента с некоторым контрольным геномом и позволяет провести скрининг всего генома. Для проведения анализа контрольная и анализируемая ДНК разрезается на малые фрагменты, которые гибридизуются с олигонуклеотидами, закрепленными в определенном месте микрочипа. Благодаря различной флюоресцентной окраске исследуемой и контрольной ДНК по цвету продукта гибридизации можно определить делеции/дупликации. Принцип работы однонуклеотидных микрочипов (SNP-array) схож с устройством CGH-array, но окрашивается и гибридизуется только тестируемый образец: на микрочипах имеются зонды длиной около 25 нуклеотидов с возможными вариантами аллелей, лучше всего гибридизуются те образцы, которые полностью комплементарны зонду. Количество полностью гибридизовавшихся фрагментов ДНК непосредственно влияет на яркость сигнала. Протяженный участок, в котором детектируется один/ни одного аллеля, отражает наличие делеции. В случае с дупликацией наблюдается противоположная тенденция. SNP-array позволяют детектировать случаи однородительской дисомии, но исследовать можно не все участки генома, а только те, в которых присутствуют SNP (т.е. имеется зависимость от распределения SNP). Разрешение для CGH-array и SNP-array от 10кб до 1 мб, в среднем 50 кб. [25, 26, 27] Эти методы до сих пор широко применяются в медицине в основном за счет хорошо отработанной методики и относительно небольшой стоимости, однако обладают рядом недостатков: сложности в определении новых и редких мутаций, низкое разрешение, зашумление данных гибридизации. Также для определения CNV в конкретных участках генома проводят количественную ПЦР в реальном времени: готовят праймеры для интересующего локуса, проводят реакцию для тестового и референсного участка с добавлением интеркалирующего/флуоресцирующего красителя, на каждом цикле ПЦР флуоресценция удваивается. Таким образом, по достижении некоторого порогового значения по номеру проведенного цикла можно узнать первоначальное количество копий таргетного участка. Например, этим способом определили вариации числа копий SMN1 и SMN2 при спинальной мышечной атрофии. [28]

1.3 Методы детекции CNV на основе данных NGS

В 2008 году вышла первая статья о секвенировании генома человека методом NGS. [29] Массовое параллельное секвенирование открыло широкие возможности для медицинской диагностики в масштабах целого генома с разрешением до одного нуклеотида. Определение CNV на основе данных NGS позволило решить сразу несколько проблем, присущих микрочипам: улучшило разрешение и покрытие, позволило уточнить границы, исследовать новые, редкие CNV и сложные структурные вариации, например, инверсий. [30] Повсеместное использование NGS в сфере геномных исследований как в медицине, так и в популяционных исследованиях, послужило стимулом для создания множества инструментов для определения CNV. Инструменты отличаются методическим подходом и разным уровнем требований к качеству NGS данных, в первую очередь, к покрытию.

Методы определения CNV по данным NGS основаны на разделении ридов (Split Read (SD)), исследовании парных ридов (Read-Pair (RP)) и глубины ридов (Read Depth (RD)), сборке генома (Assembly (AS)). Также имеет место комплексный подход: некоторые инструменты сочетают несколько методов.

Read-Pair

Метод разработан для поиска CNV на основе данных NGS с парно-концевыми чтениями. В основе метода лежит сравнение среднего размера вставки для некоторой пары ридов по сравнению с ожидаемой для референсного генома. Ожидается, что в парно-концевом секвенировании фрагменты ДНК распределены определенным образом относительно размера вставки. Если для закартированных парных ридов наблюдается несоответствие размера вставки ожидаемому среднему, для данного участка генома констатируется факт CNV. Read-pair метод не чувствителен к малым делециям/инсерциям и не подходит для определения CNV в участках низкой сложности, богатых сегментными дупликациями. Инструменты, разработанные на основе данного метода: PEMer, Hydra, Ulysses и BreakDancer. RP предполгает использование информации о парно-концевых чтениях, таким образом ограничивая выбор платформы секвенирования.

Split Read

Использует информацию о закартированных парно-концевых чтениях: если с хорошим качеством картируется на геном только один из ридов, то второй, незакартированный, рассматривается как потенциальный источник точек разрыва при CNV. Картирование ридов вокруг точек разрывов структурных вариаций позволяет четко определить границы. Данный метод плохо годится для определения широкомасштабных CNV. Однако некоторые инструменты обходят и это ограничение, применяя выравнивание Нидлмана-Вуншa, например, алгоритм Prism. Другие инструменты, основанные на данном методе: Pindel, Gustaf, SVseq2.

Read Depth

Данный метод основан на предположении о том, что между глубиной покрытия некоторого участка генома и его количеством копий имеется корреляция. Для применения этого метода в зависимости от модификации может потребоваться только исследуемый геном либо в сочетании с контрольным/средним по популяции геномом. В отличие от RP и SR методов RD позволяет определить точное число копий, тогда как 2 предыдущих метода дают информацию о границах и типе CNV. Алгоритм предполагает выравнивание ридов на референсный геном, расчет количества ридов для заданного окна, нормализацию с учетом повторяющихся элементов и различия в GC-составе и затем сегментацию генома на участки со схожим значением скорректированного числа ридов. После этого проверяют статистическую значимость полученных результатов и делают окончательные выводы о наличии и типе CNV. Разработанные инструменты: CNV-seq, BIC-seq, Cn.MOPS, CNVnator, ERDS, RDXplorer, ReadDepth, SedSeq и другие.

Assembly

Теоретически все формы структурных вариаций можно определить с помощью пересборки генома и последующего сравнения с референсным. Однако этот метод требует серьезных вычислительных ресурсов и позволяет определять только гомозиготные структурные вариации. [31] Magnolia - один из примеров реализации данного подхода.



Рис. 5 Методы определения CNV на основе данных NGS.

Также существуют комбинированные подходы, например, SVDetect - один из первых инструментов, который сочетает RP и RD методы. В последние годы именно комбинированные алгоритмы пользуются наибольшей популярностью, поскольку позволяют скомпенсировать недостатки отдельно взятых методов. [32]

Выбор алгоритма определения CNV может быть ограничен как запросом исследователя (например, нас интересует распределение только крупных CNV в нескольких популяциях), так и техническими возможностями: вычислительными ресурсами, покрытием, типом полученных чтений (одноконцевые/парные). Поэтому для практического применения особенно интересен обзор инструментов, позволяющих определять CNV в условиях ограниченных технических возможностей. Метод подсчета ридов (RD) позволяет определять крупные CNV (более 50 кб) для одноконцевых чтений даже с небольшим покрытием.



1.4 Программы, основанные на «глубине ридов» , для определения CNV при малом покрытии

QDNAseq (Ilari Scheinin, 2014).

Первоначально алгоритм был разработан для анализа зафиксированных формалином в парафине (FFPE) образцов опухолевых тканей. ДНК, полученная из таких образцов, сильно деградирована, и обычно после секвенирования покрытие не велико (0.3х-0.4х). Входные данные - одноконцевые чтения. В отличие от других алгоритмов, в QDNAseq применяется одновременная коррекция по GC-составу и картированию, а также фильтрация сложных участков, которые вносят шум в анализ, но при этом не дают значимой информации (яркий пример - прецентромерные области). На первом этапе геном разбивают на неперекрывающиеся окна фиксированного размера. Разрешение метода задается величиной окна, минимальное значение которой составляет 1 кб и зависит от покрытия. Затем удаляется 12893 бина (для разрешения 15 кб), нуклеотидный состав которых не определен для референсного генома (т.е. там стоит N). В результате получается 179187 аутосомных бинов. Далее определяется количество ридов в каждом бине и рассчитывается их медиана как функция GC-состава и картируемости. Важно, что в данном случае эти 2 параметра рассматриваются как зависимые величины. Получается некоторая поверхность с изобарами, которая затем сглаживается с помощью локальной регрессии (LOESS). Рассчитывается логарифм (log2-трансформация) отношения исходного количества ридов в окне к корректированному с помощью LOESS значению медианы количества ридов, посчитанной для окон со схожим GC-составом и картируемостью. Подобная коррекция дает хороший результат, но также порождает некоторое количество артефактов: в прецентромерных и теломерных участках наблюдаются отклонения от нормального log2-отношения. Поэтому авторы составили так называемый «черный список» участков генома, порождающих артефакты. В основном это участки с большим количеством повторов. В проекте ENCODE также содержится подобная информация, однако ее недостаточно для получения корректного результата, и при разработке QDNAseq был составлен список сложных участков генома на основе данных ENCODE и 1000Genome.

Рис. 6 Коррекция числа ридов в QDNAseq: A и C - профиль числа копий для некорректированного и корректированного числа ридов, B - медиана числа ридов на бин как функция GC-состава и картируемости, C - результат применения LOESS - сглаживания.

Для 38 образцов из проекта 1000Genomes было рассчитано среднее значение разности количества ридов и медианы. Участки генома с большим значением разности были включены в «черный список». Таким образом, для бинов размером 15 кб вместе с неаннотированными участками генома из рассмотрения было исключено 954 протяженных участка, что позволило лучше отфильтровать данные для дальнейшего анализа. В результатах всегда присутствует шум, генерируемый несколькими источниками на разных этапах, начиная с обработки образца и заканчивая картированием чтений на геном; особенно важно было оценить влияние глубины секвенирования на точность полученных CNV. Поскольку разброс для профиля вариации числа копий некоторого образца является величиной аддитивной, т. е. учитывает различные источника шума суммированием соответствующих дисперсий, шум, вносимый статистикой количества чтений на бин (1/N, где N - среднее количество ридов на бин), легко учесть в совокупности остальных составляющих.

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

Для исследования взаимосвязи общей дисперсии профиля CNV с глубиной секвенирования была построена зависимость разброса результатов для некоторого образца от статистики по количеству ридов на бин, полученной на разных покрытиях путем случайной выборки соответствующего количества ридов. Дисперсия для образца рассчитывалась по значению разности количества ридов между соседними бинами, при этом случаи, выходящие за пределы 0.1% на распределении, не учитывались при оценке дисперсии. Полученная функциональная зависимость имеет наклон 1.026 и пересечение с осью ординат в значении 0.00107, что хорошо соотносится с теоретически рассчитанными значениями (1. и 0. соответственно). Авторы также проанализировали вклад в шум процедуры приготовления библиотеки и секвенирования и пришли к выводу, что эти этапы дают малый вклад в дисперсию. Таким образом основной вклад в разброс значений профиля CNV дает статистика по количеству чтений в бине, соответственно, выбирая некоторый уровень шума, пользователь может оптимальным образом определить размер бина (разрешение). По оценкам авторов, количество ридов на бин, равное 30, и соответствующее величине окна в 15 кб для генома с покрытием 0.3-0.4, вполне подходит для анализа CNV.

QDNAseq представляет собой R-пакет, доступный в Bioconductor, принимает на вход файлы в формате BAM, и относительно быстро их обрабатывает. На выходе пользователь получает корректированное, нормированное и log2-трансформированное значение количества ридов на бин. Дальнейшая сегментация и определение типа CNV происходит с помощью встроенного пакета CGHcall, который изначально разработан для аннотации данных с микрочипов. Авторы тестировали работу алгоритма на ПК с параметром тактовой частоты 2.3 GHz и процессором Intel Core i5 CPU. Для сравнения были получены результаты для того же образца и покрытия с помощью инструмента FREEC и на платформе Agilent array CGH: величина разброса для QDNAseq оказалась не больше, чем у FREEC и меньше, чем у Agilent array CGH, и в целом профиль менее зашумленный. Применение CGHcall на заключительном этапе включает сегментацию с помощью алгоритма DNAcopy и определение принадлежности сегмента одному из 6 классов: норма, гетерозиготная и гомозиготная делеции и дупликации, амплификация. Сначала алгоритмом EM делают оценку значения вероятности для каждого класса: среднее значение log2-отношения моделируется как смесь 6 нормальных распределений. Затем по значениям постериорных вероятностей проводят классификацию по правилу Байеса .

DWAQseq (D.M. van Beek , 2012)

Авторы позиционируют данный алгоритм как хорошую альтернативу микрочипам CGH c такой же стоимостью, но лучшим разрешением и другими преимуществами NGS- детекции CNV. На первом этапе проводится коррекция ошибок, связанных с картированием ридов на повторы в геноме, и удаление ПЦР-дубликатов. Далее проводится анализ закартированных ридов, где основной упор делается на фильтрацию ложно-положительных результатов. Алгоритм основан на сравнении тестируемого генома с контрольным. В качестве контрольного использовался гибридный искусственный геном, составленный из большого числа настоящих геномов; авторы использовали коммерческий вариант, а также гибридный геном из 27 образцов, полученных из крови беременных женщин с долей плодовой ДНК порядка 5%. Результаты сравнивали для алгоритмов CNVnator [Abyzov et al., 2011], RDXplorer [Yoon et al., 2009], DWAC-seq [Koval and Guryev, 2011], CNV-seq [Xie and Tammi, 2009]) и микрочипа aCGH Agilent 180k oligo microarray. Наиболее сложным оказался подбор параметров величины окна, анализ результатов с контрольным геномом разного покрытия, выбор значения границы (трэшхолда) и фильтрация. Авторы подчеркивают, что результаты не устойчивы к пропорциональному изменению величины окна при изменении покрытия (величина окна измеряется в количестве ридов, а не в парах оснований), что выборка для конструирования гибридного генома должна быть максимально разнородной и большой (не менее 100 геномов), что полученный профиль отношения тестового количества ридов к контрольному важно визуализировать, чтобы отличить истинные CNV от артефактов в связи с наличием высокочастотной CNV в популяции. Фильтрация на завершающем этапе проводится следующим образом: если тестовый и контрольный образец содержали пустые бины (без ридов), но их было не более 20%, то предполагаемые участки генома рассматривались как CNV. Если тестовый или контрольный образец содержали пустые бины, в отличие от остальных, и при этом покрытие в контрольном образце было не менее 50%, определение CNV также считалось корректным. Бины объединялись в случае близкого значения log2-отношения.



Рис.7 Пример визуализации истинно-положительного CNV при тестировании DWAC-seq. Верхняя панель - log2 отношение числа ридов тестового и контрольного геномов, вторая сверху - число ридов в образце (после нормализации), третья панель - число ридов в контрольном образце, четвертая панель - визуализация делеции после сегментации.

Данный алгоритм был протестирован на нескольких разрешениях. Разрешение менее 0.23х дает слишком много ложноположительных CNV. Увеличение покрытия контрольного образца в три раза (до 4,5х) показало худшие результаты по сравнению с исходным. Наименьшая величина окна, при которой количество истинно-положительных максимально, а ложноположительных не велико, составляет 100 ридов на бин. Также было обнаружено, что не все CNV, определенные с помощью aCGH, являются истинно-положительными, и результаты детекции CNV с помощью микрочипов не стоит рассматривать в качестве золотого стандарта. Сравнение нескольких инструментов показало в первую очередь значительное отличие в количестве CNV. У DWAC-seq и RDXplorer оно было меньше, чем у CNVnator и CNV-seq. При этом DWAC-seq выдал большее число истинно-положительных результатов по сравнению с RDXplorer (48% и 25%). В качестве тестовых были взяты 4 образца с разным покрытием (от 0.4 до 3.6), и на них также было получено разное число CNV. При определении CNV в клинической практике важным является обнаружение патогенных вариантов. С этой задачей лучше всего справился CNVnator, однако из-за огромного количества ложных результатов авторы сочли DWAC-seq более предпочтительным. [34]

CNVnator (Abyzov A, 2011)

Популярный инструмент, который часто используется для сравнения при разработке новых алгоритмов. В основе CNVnator разделение генома на неперекрывающиеся окна фиксированной ширины, подсчет количества ридов после предварительной GC- коррекции, сегментация полученного профиля на участки предполагаемых CNV и отбор последних в соответствии со статистической значимостью. RD-сигнал для каждого бина рассчитывался по формуле

.

Где - это RD сигнала для данного участка длиной L, - среднее для гауссовского приближения RD сигнала для бинов данного размера, с - масштабирующий фактор, равный 2 для всех хромосом, кроме X и Y у мужчин (с=1). Такая нормализация не устойчива к выбросам, связанным с повторами в геноме и плохой картируемостью ридов, на что указывают сами авторы. RD для участка рассчитывается суммированием для бинов. Разделение на сегменты производили методом сдвига среднего, изначально разработанного для компьютерного зрения и нашедшего применение в анализе данных микрочипов. Сначала, сравнивая RD - сигнал данного и соседнего участка, определяли направление вектора сдвига среднего. Вектор был направлен сторону бина с наиболее близким RD - сигналом. Затем граница сегмента определялась в том случае, если 2 соседних вектора имели противоположное направление, но при этом не указывали друг на друга.

Рис.8 Метод сдвига среднего.

Сложным этапом для данной задачи является проблема неуникального картирования ридов. В некоторых алгоритмах фрагменты генома с неуникально картируемыми ридами просто не рассматриваются. Авторы СNVnator выбрали иную стратегию, они учитывают все участки генома: неуникально картируемые риды оставляют и случайным образом выбирают место на геноме из возможных. Так геном оказывается покрыт относительно равномерно.

Также авторы тестировали алгоритм для определения наилучших параметров: оптимальная ширина бина для покрытия 4-60х составила 500 п.о. С одной стороны, авторы рекомендуют инструмент как для больших, так и для малых покрытий. С другой, сравнение с aCGH показали на малом покрытии несколько худшие результаты. Для 90% граничных точек разрешение составило 200п.о. Анализ полного генома занял несколько часов на 2.5-GHz Intel Core 2 Duo CPU. [35]

GROM-RD (Grigoriev A, 2015 )

Некоторые участки генома представлены повторяющимися элементами, количество ридов, закартированных туда, может не в полной мере отражать покрытие генома. Авторы предложили маскировать подобные сложные участки и проводить расчет для 2 вариантов (маскированного генома и обычного), чтобы затем иметь возможность сравнивать результаты. GC-коррекция проводится в расчете на одно основание, где j - основание, которое может повлиять на RD основания i, - вес основания j, равный сумме средней длины ридов с разными координатами начала, которые перекрывают основания i и j, равно 1, если основание G/C и 0 для других. Такая коррекция обусловлена корреляцией GC-состава с эффективностью ПЦР всего амплифицируемого фрагмента ДНК, а не рассматриваемого в данный момент участка фиксированной длины (величины бина). В случае стандартной поправки не учитывается GC-состав последовательности, лежащий вне бина, но, возможно, влияющей на количество ридов для данного окна. Стандартная процедура нормализации количества ридов по нуклеотидному составу, применяемая в большинстве алгоритмов, тем не менее дает неравномерное распределение дисперсии количества ридов по GC-составу. Поэтому авторы предлагают квантильную нормализацию. Также GROM-RD - это первый инструмент, учитывающий смещение RD в связи АТ-повторами длиной более 10 п.о. Для получения более точных границ CNV применяется метод скользящего окна с переменной шириной. Инструмент показывает значительное преимущество по стандартным параметрам оценки точности при сравнении с CNVnator и RDXplorer на геномах с малым и большим покрытием. Расчет одного генома занимает порядка 2-3 часов на процессоре (Intel Xeon E31270, 3.4 GHz) c 16Gb оперативной памяти. [36]

CN.MOPS (G Klambauer, 2012)

Подготовительный этап, как и для других алгоритмов, заключается в картировании ридов. Допускается несовпадение в двух позициях, если невозможно выбрать одну наилучшую позицию, рид картируется случайным образом на одну из возможных. Далее для подсчета числа ридов геном делится на неперекрывающиеся сегменты. По умолчанию длина сегментов одинакова, но это не принципиально, так как для каждого сегмента алгоритм строит свою модель, и на следующих этапах сегментация основана на ожидаемом количестве копий, а не на длине сегмента. Модель предполагает пуассоновское распределение количества ридов после нормализации данных. Новизна данного метода состоит в моделировании вариации числа ридов, для того чтобы отделять разброс вследствие технических шумов от разброса в связи с наличием CNV. Модель предполагает линейную зависимость между средним количеством ридов в сегменте и CNV. В отличии от других методов, cn.MOPS дает на выходе количество копий с доверительным интервалом. Алгоритм предполагает использование множества образцов (как минимум 6).



Рис. 9 Общая схема поиска CNV на основе данных NGS. Отличительные особенности алгоритма cn.MOPS.

RDXplorer (Yoon et al., 2009) основан на event-wise testing (EWT) алгоритме, оценивает CNV в непересекающихся интервалах в индивидуальном геноме. EWT алгоритм быстро находит участки с высоким/низким RD, удовлетворяющие статистическим критериям, и затем кластеризует их. [38]

Авторы XCAVATOR (Magi A, 2017) показали, что RD данные лучше моделируются отрицательным биномиальным распределением, чем Пуассоновским. Количество ридов было нормировано через медиану, что минимизирует влияние GC-состава и картируемости. Для сегментации генома была использована модель сдвига уровня (shifting level models (SLM)). Этот метод основан на гауссовском распределении, поэтому сигнал должен быть нормально распределен (нормализованные данные нужно log-трансформировать). [39]

Таким образом, все Read-depth алгоритмы можно условно разделить на 2 класса: предполагающие наличие контрольного образца либо не требующие его. Каждый метод предполагает сегментацию генома, некоторую GC-коррекцию, модель статистической оценки количества ридов. Основные трудности возникают на этапе GC- коррекции, определении источника вариабельности количества ридов и обработки повторов в геноме. Каждый алгоритм предлагает свои методы решения.

Для поиска CNV в геноме с экстремально низким покрытием также существуют коммерческие инструменты, например, Golden helix VarSeq. По словам разработчиков, данный инструмент эффективно находит CNV размером более 1 мб в геноме с покрытием 0.02.

1.5 Базы данных CNV

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

dbVar - база данных, поддерживаемая NCBI и содержащая информацию о патогенных и нейтральных CNV более 50 п.о. в геноме человека и других видов со ссылкой на исследование; синхронизована с DGVa. Формат данных: исследование, координаты, тип CNV. [40]

DGVa основана Европейским Институтом Биоинформатики (EIB), также, как и dbVar, предполагает прямое размещение информации исследователями. Подтверждение для размещения дается после проверки в соответствующей области научной литературы. Эта база данных подходит скорее для размещения экспериментальных данных конкретных исследований в свободном доступе, чем для быстрого получения информации об интересующей CNV. [40]

DGV поддерживается Канадским центром прикладной геномики и содержит информацию о непатогенных CNV более 50 п.о. Имеет удобный интерфейс с подробной информацией по СNV. Получает данные из dbVar и DGVa и проверяет их. В настоящий момент является самой обширной курируемой контрольной базой данных, включает 72 исследования, насчитывает 16,126,754 варианта. CNV, подтвержденные более 2 исследованиями для двух различных образцов, включены в золотой стандарт DGV, который регулярно обновляется с учетом версии сборки референсного генома. [41]

DECIPHER разработана в Институте Сенгера как база данных патогенных субмикроскопических структурных вариаций (в настоящий момент насчитывается более 29000 патогенных CNV связанных более чем с 75000 патологическими фенотипами). Версия референсного генома - GRCh37. Имеет удобный браузер: по координатам можно получить полную информацию об имеющихся вариантах в контрольной группе, о специфике фенотипа при наличии патологии, и наоборот.

ISCA (International Standards for Cytogenomic Arrays) содержит патогенные, нейтральные, предположительно нейтральные/патогенные и не определенные CNV, подтвержденные комитетом ISCA (клинической лабораторией). Дата последнего обновления - январь 2013 года.

ECARUCA содержит клиническую и цитогенетическую информацию по редким хромосомным нарушениям, включая микроделеции и микродупликации. В настоящее время проводится подготовка к обновлению базы данных как с технической точки зрения, так и с содержательной. К лету 2018 года планируется дополнить базу данных 3000 случаями, таким образом итоговое их количество будет около 8000. Представленная информация предназначена для врачей, для доступа необходима регистрация.

Некоторые базы специализируются по определенным заболеваниям.

SFARI - база данных по CNV при расстройствах аутистического спектра. Использует подход системной биологии: предполагаемые гены-кандидаты оцениваются согласно правилам, разработанным консультирующими специалистами в области аутизма, и затем классифицируются по нескольким категориям. Информация в базе данных разделена на несколько модулей:

- гены человека (human gene) - содержит аннотированный список генов, исследуемых когда-либо в контексте аутизма

- CNV - содержит информацию о моногенных и мультигенных делециях и дупликациях и их связь с аутизмом.

- модели на животных (animal models) - здесь представлена информация о линиях генно-модифицированных мышей и крыс, используемых в качестве моделей в исследованиях по данному вопросу

- белковое взаимодействие - включает все известные взаимодействия белков генов, вовлеченных в аутизм (в том числе с нуклеиновыми кислотами).

База данных постоянно обновляется, информация о новых CNV берется из научных статей. Каждая ассоциированная с аутизмом CNV считается валидированной, если она подтверждена независимо другой методикой. Обладает удобным современным интерфейсом. Также существует курируемая база данных Autism Chromosome Rearrangement Database.

SZGR2.0 (Schizophrenia gene resource 2) помимо точечных мутаций, дифференциальной экспрессии и статуса метилирования включает также информацию о CNV (12 локусов), ассоциированных с шизофренией. Последняя версия вышла в июне 2016 года и поддерживается Научным центром здоровья при Техасском университете в Хьюстоне.

Также некоторые веб-ресурсы (OMIM), не специализирующиеся на данных по CNV, содержат информацию о связи отдельных участков генома с возможным патологическим фенотипам. Т. е. подавая на вход координаты CNV, можно получить информацию об ассоциированных заболеваниях.

Поскольку статистика по CNV в норме может иметь различия по популяциям, а наиболее обширная информация представлена для населения Европы и Америки, появилось несколько этноспецифичных баз данных:

- CNV control database - содержит информацию с микрочипов affymetrix 6.0 о CNV в геномах 160 японцев. Поддерживает версию генома hg19.

- Korean Genomic Variant Database (KGVDB) - контрольная база данных. Данные aCGH и SNP-микрочипов получены в ходе исследования 4694 геномов корейцев. Поддерживает версию генома hg18.

- Thai CNV Database - контрольная база данных, созданная по результатам исследования на SNP-микрочипах 3017 представителей Таиланда. Поддерживает версию генома hg18.

Кроме того, существуют программы (AnnTools, ANNOVAR), позволяющие легко аннотировать интересующий геномный локус. Например, ANNOVAR - программный пакет в открытом доступе, дает информацию о генах, о гомологах белков у различных видов, пересечении с сегментными дупликациями и др. Инструмент постоянно обновляется. Существует также коммерческая реализация подобных инструментов: VarSeq1.4.7 от компании Golden Helix.



Глава 2. Материалы и методы

2.1 Образцы для поиска CNV

Для составления каталога CNV в российской популяции были использованы данные, полученные в ходе неинвазивного пренатального ДНК-скрининга анеуплоидий (НИПС) в лаборатории молекулярно-генетических методов ФГБУ НМИЦ АГП им. В.И. Кулакова. [48] Это отсеквенированная на приборе Ion Proton или Ion S5 внеклеточная ДНК, полученная из плазмы крови беременных женщин, с долей фетальной ДНК 4-20%. Полученные в результате секвенирования риды были картированы на референсный геном человека GRCh37 с помощью сервера, управляющего прибором (Torrent Server, TMAP v 5.4). Также для последующей проверки устойчивости результатов к картированию на разные варианты референсного генома риды нескольких образцов были выровнены на геном GRCh38. Была проведена фильтрация ПЦР дубликатов (Torrent Server, FilterDuplicates v 5.4). Bam-файл для каждого образца содержал около 6 млн. ридов длиной 150-180 п.о., т. е. покрытие составило примерно 0,35х. Вся работа, предшествующая получению bam-файлов, была проведена ранее и выходит за рамки данного исследования.

2.2 Тестирование алгоритмов поиска CNV

Для выбора оптимальных инструментов, подходящих для работы с малым покрытием, было протестировано несколько программ: CNVnator, cn.MOPS, GROM-RD, XCAVATOR, QDNAseq, DWAC-seq. Все они доступны для свободного пользования. Тестирование проводилось в 2 этапа. На первом этапе для расчета были использованы образцы с известной делецией (хромосома 13, 48 Mб) и дупликацией (хромосома 9, 45 Мб) с количеством чтений порядка 1,5 млн и 1,3 млн, величина бина задавалась равной 30, 100 и 1000 кб для каждого инструмента, за исключением DWAQ-seq. Для cn.MOPS поиск CNV проводился на 14 геномах, 2 из которых содержали вышеуказанные CNV, количество чтений в остальных составляло около 5-6 млн. На втором этапе все программы были запущены на образце с локальным повышением в покрытии на участке 64500000-69000000 хромосомы 4 с общим количеством ридо 5,6 млн. Тестирование проводилось на величине бина 100 кб, для DWAQ-seq этот параметр не задается, но для него было выбрано покрытие референсного генома 1,2х. Расчет для cn.MOPS проводился на 12 геномах со схожим количеством ридов.

Рис. 10 Изменение в покрытии для образца, использованного на втором этапе тестрования.

CNVnator, ХCAVATOR, GROM-RD: имеют интерфейс командной строки. Контрольный образец не требовался.

CNVnator

Входные данные: образец в формате bam-файла с индексом и последовательность референсного генома по хромосомам в формате fasta.

Для получения данных о CNV необходимо было провести последовательное выполнение нескольких команд, суть которых состоит в следующем:

- извлечение ридов из bam-файла

- построение распределения ридов

- статистический анализ

- сегментация

- определение CNV

Данные на выходе: тип CNV, координаты, размер, нормированное число ридов, e-val1, рассчитанное из t-test статистики, e-val2 - вероятность того, что число ридов находится в «хвосте» гауссовского распределения.

XCAVATOR

Входные данные: bam-файл с индексом, BigWig-файл, референсный геном в формате fasta c индексом.

Для конвертации bam-файла в BigWig-формат были использованы команды genomeCoverageBed из пакета bedtools - перевод bam-файла в BedGraph, сортировка полученного файла, bedGraphToBigWig из пакета kentUtils - перевод BedGraph в BigWig-файл. BigWig - это бинарное представление формата wiggle, который содержит дополнительную информацию/статистику для образца, например, GC-состав, и обычно используется для визуализации данных в IGV.

Алгоритм работает в несколько этапов:

- сбор информации по GC-составу и картируемости для окон заданного размера

- подсчет числа ридов и нормализация

- сегментация методом сдвига среднего и классификация сегментов.

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

cn.MOPS

Программа представляет собой R-пакет, доступный в Bioconductor. Согласно описанию, требуется не менее 6 геномов (bam-файлы с индексом) для эффективной работы. Расчет проводился для аутосом и Х-хромосомы. На первом этапе выполняется команда getReadCountsFromBAM, которая подсчитывает количество ридов в bam-файлах, далее выполняется команда, запускающая основной алгоритм поиска CNV для каждого образца. Инструмент использует информацию одновременно со всех образцов для определения CNV в каждом из них. Сn.MOPS выдает результаты в формате GRanges, число копий для каждого сегмента при этом можно представить в целочисленном виде: CN0, CN1 и т.д. Полученные результаты легко переводятся в DataFrame (команда as.data.frame()).

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

DWAC-seq

Программа DWAQ-seq написана на языке perl и запускается из командной строки. Для работы с данным инструментом был составлен контрольный геном из 128 образцов. Для объединения геномов применялась команда merge из программы samtools, c последующей сортировкой, индексацией и уменьшением количества ридов с помощью команды samtools view -h -s x.y, где x обозначает значение для запуска генератора случайных чисел (выбирается случайным образом функцией RANDOM, нужен для генерации случайной выборки ридов) , у - количество ридов в финишном геноме как доля от первоначального (например 25 - соответствует 25% от исходного bam-файла) Слияние геномов проводилось рекурсивно по 2 или 4 генома с последующей сортировкой, индексацией и уменьшением покрытия. Данная процедура повторялась несколько раз, поскольку ограничения по объему дисковой памяти и количеству временных файлов не позволяют объединять, сортировать и индексировать большое количество геномов за одну итерацию. Всего было получено 2 гибридных генома с покрытием около 1.2х и 2.4x.

...