Использование графовых моделей для биоинформатического анализа гипервариабельных биологических последовательностей тема диссертации и автореферата по ВАК РФ 05.13.18, кандидат наук Сафонова, Яна Юрьевна

• Сафонова, Яна Юрьевна
• кандидат наук
• 2016

Оглавление диссертации кандидат наук Сафонова, Яна Юрьевна

Оглавление

Стр.

Введение

Глава 1. Общие аспекты биоинформатики

1.1 Появление технологий секвенирования

1.2 Эпоха секвенирования нового поколения (СНП)

1.3 Омики

1.4 Секвенирование ДНК и биоинформатика

1.5 Формулировка вычислительных задач, связанных с анализом ДНК

1.6 Сборка геномов

1.6.1 OLC подход

1.6.2 Подход на основе графа де Брюйна

1.6.3 Области применения различных подходов для решения задачи сборки генома

1.7 Иммуногеномика

1.7.1 Адаптивные иммунные репертуары

1.7.2 Иммунопротеогеномика

1.8 Настоящая работа

1.8.1 Сборка высокополиморфных диплоидных геномов

1.8.2 Построение адаптивных иммунных репертуаров

Глава 2. Сборка высокополиморфных диплоидных геномов

2.1 Сборка геномов с помощью графа де Брюйна

2.2 Высокополиморфные диплоидные геномы

2.3 Сборка высокополиморфных диплоидных геномов

2.4 Алгоритм dipSPAdes

2.4.1 Общая идея алгоритма dipSPAdes

2.4.2 Маскирование полиморфизмов в графе де Брюйна

2.4.3 Конструирование консенсусных и двойных контигов

2.5 Анализ эффективности работы dipSPAdes

2.5.1 Симулирование реалистичного референсного диплоидного

генома с помощью S. commune

Стр.

2.5.2 Сравнение dipSPAdes с другими алгоритмами сборки

2.5.3 Валидация диплоидной сборки с помощью гапломов S. commune

2.5.4 Биологические применения dipSPAdes

Глава 3. Задачи вычислительной иммуногеномики

3.1 Адаптивная иммунная система и репертуар антител

3.2 Секвенирование адаптивных иммунных репертуаров

3.3 Проблема построения репертуаров и программный комплекс

IgRepertoireConstructor

3.3.1 Представление репертуара антител

3.4 Алгоритм IgRepertoireConstructor

3.4.1 Идея алгоритма IgRepertoireConstructor

3.4.2 Использование графа Хэмминга из прочтений для построения репертуара антител

3.4.3 Иммунопротеогеномный поиск

3.4.4 Сравнение геномных и протеомных кратностей

3.5 IgSimulator

Глава 4. Дальнейшее развитие предложенных задач

4.1 Сборка высокополиморфных диплоидных геномов

4.2 Вычислительная иммуногеномика

4.2.1 Оценка качества алгоритмов построения репертуара антител

4.2.2 Эволюционный анализ репертуара антител

4.2.3 Анализ эффекта многих цепей

4.3 Связь сборки высокополиморфных диплоидных геномов с задачами вычислительной иммуногеномики для популяционного анализа локуса иммуноглобулинов

Заключение

Список сокращений и условных обозначений

Словарь терминов

Стр.

Список литературы

Список рисунков

Список таблиц

Введение диссертации (часть автореферата) на тему «Использование графовых моделей для биоинформатического анализа гипервариабельных биологических последовательностей»

Введение

Актуальность темы. Появление технологий секвенирования открыло новые горизонты биологии и перенесло анализ биологических данных из экспериментальной лаборатории в область вычислительного анализа. Последовавшее накопление биологических данных колоссального объема привело к развитию биоинформатики, огромной междисциплинарной науки, объединяющей компьютерные науки, математику, статистику и разработку программного обеспечения для решения разнообразных биологических задач. Одним из важнейших объектов биоинформатических исследований являются вариации в биологических последовательностях: нуклеотидных последовательностях ДНК и РНК, а также аминокислотных последовательностях белков. Особенно трудным является анализ гипервариабельных последовательностей, например, высокополиморфных диплоидных геномов или адаптивных иммунных репертуаров.

В данной работе описываются решения двух проблем, возникающих при проведении анализа гипервариабельных биологических последовательностей: проблема сборки высокополиморфных диплоидных геномов и проблема построения адаптивных иммунных репертуаров. В обоих случаях в связи с чрезвычайно высоким уровнем разнообразия в биологических последовательностях стандартные подходы не в состоянии предоставить универсальное решение. При этом обе задачи преследуют одну общую цель: наиболее точно сохранить естественное разнообразие последовательностей и отличить их от артефактов подготовки биологического образца. Для решения этих задач были разработаны программные комплексы: dipSPAdes, IgRepertoireConstructor и IgSimulator. В результате вычислительных экспериментов было показана их эффективность. Предложенные решения оказались первыми доступными инструментами для поставленных задач.

Целью работы являются

1. разработка и реализация программного комплекса для сборки высокополиморфных диплоидных геномов, тестирование разработанного алгоритма на реальных и симулированных данных;

2. разработка и реализация программного комплекса для по строения репертуара антител с помощью данных иммуносеквенирования и его валида-ция с помощью иммунопротеогеномного анализа;

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

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

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

Все программные комплексы были реализованы на языках программирования C++ и Python.

Практическая ценность и научная новизна. Результатами выполненной работы стали следующие программные комплексы:

- dipSPAdes: инструмент для сборки высокополиморфных диплоидных геномов (доступен как часть геномного сборщика SPAdes: bioinf.

spbau.ru/ru/spades).

- IgRepertoireConstructor: алгоритм для построения адаптивных иммунных репертуаров из данных иммуносеквенирования (yana- safonova. github.io/ig_repertoire_constructor).

- IgSimulator: универсальный инструмент для симуляции иммунных репертуаров и данных иммуносеквенирования (yana-safonova.github. io/ig_simulator).

Предложенные алгоритмы стали первыми доступными инструментами для решения поставленных задач.

Апробация работы. Основные результаты работы неоднократно докладывались семинарах лаборатории «Центр алгоритмической лаборатории» (СПбГУ), а также на следующих конференциях и встречах:

1. Y-Tools, a toolkit for construction and analysis of adaptive immune repertoires. Устный доклад на конференции RepSeq 2016. Гаага, Нидерланды, 2016.

2. Big Data challenges in analysis of adaptive immune repertoires. Устный доклад на семинаре в UCSD, Сан-Диего, США. 2016

3. New algorithmic challenges in analysis of adaptive immune repertoires. Устный доклад на семинаре в компании Genentech. Сан-Франциско, США. 2016

4. A novel toolkit for analysis of adaptive immune repertoires using immunosequencing and mass-spectra data. Постерный доклад на конференции RECOMB 2016. Санта-Моника, США.

5. IgTools: a toolkit for analysis of antibody repertoire from immunosequencing and mass spectra data. Постерный доклад на конференции Immunogenomics 2015. Хантсвилл, США.

6. IgRepertoireConstructor: a novel algorithm for antibody repertoire construction and immunoproteogenomics analysis. Устный доклад на конференции ISMB/ECCB 2015. Дублин, Ирландия.

7. IgTools: a toolkit for construction of antibody repertoire using immunosequencing data. Постерный доклад на конференции ISMB/ECCB 2015. Дублин, Ирландия.

8. dipSPAdes: assembler for highly polymorphic diploid genomes. Устный доклад на конференции RECOMB 2014. Питтсбург, США.

9. New Frontiers of Genome Assembly with SPAdes 3.0. Постерный доклад на конференции JGI User Meeting 2014, Уолнат-Крик, США.

10. dipSPAdes: assembler for highly polymorphic diploid genomes. Устный доклад на семинаре в Joint Genome Institute. Уолнат-Крик, США, 2014.

11. NGS assembly of highly polymorphic diploid genomes. Постерный доклад на конференции RECOMB 2013. Пекин, Китай.

Публикации. Основные результаты по теме диссертации изложены в 4 публикациях [1—4], изданных в журналах, индексируемых базами данных SCOPUS и Web of Science.

Личный вклад. В статье [1] автор играл главную роль в разработке, реализации и тестировании алгоритма dipSPAdes. Статья [2] была написана в соавторстве с лабораторией эволюционной геномики МГУ, в ней автор принимал участие в вычислительных экспериментах, связанных со сборкой популяции S. commune, на основе которых был выполнен популяционный анализ статьи. В статье [3] автор предложил и реализовал алгоритм IgRepertoireConstructor для построения репертуара антител из данных иммуносеквенирования с помощью Хэмминг графа. Автор также принимал участие в экспериментах и обсуждениях, связанных с им-мунопротеогеномным поиском. В статье [4] автор выполнил разработку и реализацию алгоритма IgSimulator.

Объем и структура работы. Диссертация состоит из введения, четырёх глав и заключения. Полный объём диссертации составляет 71 страницу, включая 14 рисунков и 4 таблицы. Список литературы содержит 90 наименований.

Первая глава описывает основные объекты изучения биоинформатики и показывает актуальность выбранной темы. Далее формулируются задачи сборки высокополиморфных диплоидных геномов и построения адаптивных иммунных репертуаров.

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

Третья глава посвящена задачам вычислительной иммуногеномики. В начале главы дается описание адаптивных иммунных репертуаров — основных компонентов адаптивной иммунной системы — и приводится несколько формулировок вычислительных задач, связанных с их построением. Далее приводятся описания алгоритма IgRepertoireConstructor для построения репертуара антител с помощью кластеризации графа Хэмминга и процедуры иммунопротеогеномного поиска.

Четвертая глава содержит обсуждение полученных результатов и описывает направления дальнейших исследований.

Глава 1. Общие аспекты биоинформатики

1.1 Появление технологий секвенирования

Появление технологий секвенирования открыло биологии новые горизонты, что в свою очередь привело к развитию ряда новых областей современной биологической науки. В 1990 году был начат проект «Геном человека». Это был первый в своем роде международный проект, целью которого было просеквени-ровать, собрать и провести анализ всей ДНК человека. В результате приложенных усилий была собрана почти вся последовательность генома, а его аннотация позволила биологам понять природу ряда комплексных заболеваний с помощью ге-нотипирования вирусов и бактерий, а также идентифицировать онкологические мутации пациентов, подбирать индивидуальные средства лечения и с большей точностью предсказывать результаты его действия. Секвенирование ДНК показало себя как незаменимый инструмент геномного анализа и стало основополагающей технологией для многих биологических исследований: от персонализированной медицины до изучения эволюции. Более того, были просеквенированы и проанализированы геномы модельных организмов и наиболее распространенных возбудителей заболеваний, например, бактерия Mycoplasma capricolum (Bork и др., 1995 год [5]), пекарские дрожжи Saccharomyces cerevisiae (Goffeau и др., 1995 год [6]), кишечная палочка Escherichia coli (Blattner и др., 1997 год [7]), нематода Caenorhabditis elegans (C. elegans Sequencing Consortium, 1998 год [8]), дрозофила обыкновенная Drosophila melanogaster (Ketchum и др., [9] и Myers и др. [10], 2000 год), мышь домовая Mus musculus (Waterson и др., 2002 год [11]) и многие другие. Появление технологий секвенирования ДНК в корне изменило методологию современной биологии и стало ключевым элементом таких областей науки как: персональная медицина, судебно-медицинские прикладные науки, биотопливо, агрокультура, животноводство, биопереработка, биоархеологоия, антропология и эволюция.

1.2 Эпоха секвенирования нового поколения (СНП)

Первые технологии секвенирования, использовавшиеся, в том числе и в проекте «Геном человека», были основаны на методе Сэнгера, который производит длинные (до 900 нуклеотидов) и очень точные прочтения. К сожалению, этот подход был очень дорогостоящим, и практически вся работа проводилась вручную. За один запуск анализировалось « 380 проб ДНК, на что уходило несколько часов и тратилось около 2,400 долларов США. Всего на секвенирование генома человека потребовалось почти 3 миллиарда долларов США. По причине трудоемкости и высокой стоимости метод Сэнгера плохо подходил для работы с большими проектами. В середине 90-х годов начался период быстрого развития технологий секвенирования, так называемая эпоха секвенирования нового поколения (СНП). С появлением методов СНП стоимость процесса секвенирования значительно снизилась, открыв дорогу масштабному высокопропускному секве-нированию. С помощью этого способа стало возможно секвенировать большие объемы ДНК или РНК проб по сравнительно низкой цене. Семейства СНП технологий включают следующие платформы: SOLiD (секвенирование на основе лиги-рования), SOLEXA/Шumma (секвенирование путём синтеза), 454 (пиросеквени-рование), ^п^гш^ (ионное полупроводниковое секвенирование). По сравнению с методом Сэнгера все эти технологии производят короткие прочтения по низкой цене и с высокой скоростью (Таблица 1). Разнообразие существующих СНП технологий и их доступность значительно повлияли на скорость развития геномики (науки, посвященной изучению генома) и биоинформатических подходов для обработки данных секвенирования.

1.3 Омики

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

Таблица 1 — Сравнение параметров секвенирования с помощью метода Сэнгера и самых популярных технологий СНП.

RL (bp) — длина прочтения в bp, Rs / run — количество прочтений на запуск, Cost / 1 Mbp — стоимость 1 Mbp, RT — продолжительность работы.

Методы

RL (bp) Rs / run

ER

Cost /1 Mbp

RT

SOLiD

50+50

SOLEXA / < 300 Illumina

454 700

IonTorrent < 400

Sanger < 900

1.2 — 1.4 Gbp больше ошибок в А-Т бо- $ 0.13 1 — 2 неде-

гатых регионах. В сред- ли нем < 0.1 %

< 6 Gbp замены (0.1 %) $ 0.15 1 день

1 МЬр вставки, удаления (1 %) $ 10 24 часа

< 80 МЬр вставки, удаления (1 %) $ 1 2 часа Ж/А замены (0.001 %) $ 2400 < 3 часов

кул РНК), сравнительная геномика (сравнительные исследования структуры и функций геномов разных организмов), иммуногеномика (изучение компонентов иммунной системы). Стремительная эволюция СНП технологий приводит к постоянному развитию методов анализа биологических последовательностей. Так, например, технология Illumina MiSeq, выпущенная в 2013 г., производит прочтения, чья длина позволяет осуществить полноразмерное сканирование адаптивных иммунных репертуаров, представляющих собой набор гипервариабельных последовательностей циркулирующих антител и T-клеточных рецепторов. Развитие иммуногеномики с помощью данных секвенирования Illumina MiSeq стало ключевым этапом развития персонализированной медицины, оценки эффективности вакцин и других средств лечения аутоиммунных, онкологических и инфекционных заболеваний. Однако, даже такие современные технологии, как Illumina MiSeq, не могут идеально решить проблему восстановления генома высокополиморфных диплоидных организмов, чья частота геномных мутаций соизмерима с уровнем вариабельности антител и Т-клеточных рецепторов. Таким образом, результатом развития биотехнологий является постоянная потребность в создании новых алгоритмов для их обработки.

1.4 Секвенирование ДНК и биоинформатика

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

1.5 Формулировка вычислительных задач, связанных с анализом ДНК

Молекулы ДНК и РНК состоят из простых единиц, называющихся нуклео-тидами. Каждый нуклеотид молекулы ДНК однозначно характеризуется одним из следующих нуклеотидных оснований: аденин (А), цитозин (С), гуанин (С) и тимин (Т). В случае молекул РНК тимин заменяется урацилом (и). Таким образом, молекулы ДНК и РНК можно представить как цепочки, состоящие из букв {А, С, С, Т} (не уменьшая общности, мы будем считать Т равноценной и). Целью секвенирования молекул ДНК/РНК является восстановление последовательности нуклеотидных оснований в них. В результате этого процесса создается набор так называемых прочтений — коротких фрагментов исходной молекулы. Прочтения также представляют собой цепочку символов в алфавите {А,С,С,Т}. В идеале прочтения являются точной копией участка исходной молекулы. К сожалению,

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

1.6 Сборка геномов

Прочтения являются очень короткими подстроками генома, поэтому часто их длины не хватает для решения биоинформатических задач. Например, для аннотации геномной последовательности нужны фрагменты генома, превосходящие по длине последовательности генов, которые могут достигать до нескольких десятков тысяч нуклеотидов. В результате, возникает задача склеивания геномных прочтений в более длинные фрагменты генома, которые могут быть использованы для дальнейшего анализа. Процесс такого восстановления генома называется сборкой генома. Чтобы восстановить цепочку из коротких прочтений, требуется найти продолжение для каждого из них и склеить перекрывающиеся прочтения в единую последовательность. За счет дупликации генов и присутствия транс-позонов, в структуре генома встречается множество повторов (Рисунок 1.1а). Это значительно усложняет процесс сборки геномов. Обычно результатом сборки геномов является не единая последовательность, а набор неперекрывающихся фрагментов генома (контигов). Алгоритмы, созданные для сборки геномов, направлены на вычисление длинных и точных контишв. Существует несколько устоявшихся подходов решения проблемы сборки геномов: OLC подход (от английского overlap-layout-consensus), основанный на построении графа перекрытий, и подход, основанный на построении графа де Брюйна.

1.6.1 OLC подход

Изначально OLC подходы были разработаны для сборки длинных прочтений, созданных по методике Сэнгера. OLC алгоритмы строят графы перекрытия, в которых прочтения соответствуют вершинам. Две смежные в графе вершины означают перекрывающиеся прочтения. Поиск неразветвленных путей в построенном графе перекрытий позволяет выделить контиги (Рисунок 1.1 б). Подход OLC реализован в нескольких программных продуктах: PHRAP [12], GAP [13], TIGR [14], ARACHNE [15] и CELERA [10]. В худшем случае, граф перекрытий может быть практически полным (то есть, содержать O(n2) ребер, где n — количество исходных прочтений), что делает нецелесообразным применение этого подхода к данным высокопропускных технологий СНП.

1.6.2 Подход на основе графа де Брюйна

Для сборки геномов из коротких прочтений, созданных с помощью технологий СНП, было предложено использовать графы де Брюйна. Вершинами графа де Брюйна являются возможные k-меры (строки длины k), выделенные из исходных прочтений. Два k-мера, соединяются в графе ребром, если они являются префиксом и суффиксом k + 1-мера, также представленного в исходных прочтениях (Рисунок 1.1в). Построение контигов соответствует поиску путей в графе де Брюйна, консистентных исходным прочтениям. Идея использования графа де Брюйна была предложена в 2001 году [16] и до сих пор широко используется в современных программных комплексах для сборки генома: Velvet [17], ABySS [18], IDBA-UD [19], SOAPdenovo [20], SPAdes [21].

R

В

R

С

R

D

т>

-

---^ с

б) в)

а) Геном состоит из четырех фрагментов А, В, С и И, разделенных повтором Я. Прочтения, соответствующие геному показаны под ним. Некоторые прочтения покрывают повтор Я, в то время как другие покрывают фрагменты А, В, С и И. б) Граф перекрытий, построенный из прочтений, показанных на Рисунке а). В результате поиска неперекрывающихся однозначных путей можно выделить пять контигов: А, В, С, И и Я. в) Граф де Брюйна, построенный из прочтений, показанных на Рисунке а), иллюстрирует геном (длина А>мера не превышает длину повтора Я и фрагментов А, В, С и £>).

Рисунок 1.1 — Основные подходы к решению задачи сборки геномов.

1.6.3 Области применения различных подходов для решения задачи сборки

генома

Оба подхода — OLC и де Брюйн граф — могут быть использованы для сборки данных различных технологий секвенирования и разных проб (единичные клетки, метагеномы, диплоидные геномы). При этом только подход с использованием графов де Брюйна позволяет наглядно отобразить структуру геномных повторов (Рисунок 1.1 в). Благодаря этому свойству, подход, использующий граф де Брюйна, легко адаптируется для решения других задач биоинформатики: сравнительный анализ геномов [22], V(D)J классификация антител и Т-клеточных рецепторов [23], de novo классификация повторов [24] и сборка транскриптомов [25].

1.7 Иммуногеномика

1.7.1 Адаптивные иммунные репертуары

Объектами изучения иммуногеномики являются наборы гипервариабельных компонентов адаптивной иммунной системы: циркулирующих антител и Т-клеточных рецепторов (вместе называемых адаптивными иммунными репер-туарами). Адаптивные иммунные репертуары представимы в виде множества последовательностей (нуклеотидных строк, каждая из которых характеризует уникальное антитело или Т-клеточный рецептор) с кратностями (число, отражающее сколько раз данное антитело/Т-клеточный рецептор встретился в репертуаре). Адаптивные иммунные репертуары отражают текущее состояние организма и поэтому являются важным объектом в биомедицинских исследованиях. К примеру, мониторинг изменения репертуаров антител после вакцинации или иммунотерапии является стандартным шагом при разработке вакцин и лекарств. Однако, несмотря на важную роль в биологии и медицине, анализ репертуаров антител и Т-клеточных рецепторов остается сравнительно малоизученной задачей в иммуногеномике. Вплоть до 2009 года вычислительный анализ антител проводился с помощью подходов, основанных на протеомике [26], без использования технологий секвенирования ДНК/РНК. Weinstein и др. [27] первыми продемонстрировали возможность использования секвенирования для анализа репертуаров антител, тем самым положив начало эпохе СНП в иммуногеномике. Несмотря на то, что вслед за этим исследованием последовало множество других [28-32], до 2012 года не было предпринято ни единой попытки объединить технологии СНП и масс-спектрометрии (далее МС) для анализа репертуара антител. Отсутствие этой интеграции (далее называемой иммунопротеогеномикой) и общедоступных инструментов для предобработки данных иммуносеквенирования было преградой на пути к развитию новых биомедицинских подходов, которые могли бы работать c большими и сложными репертуарами антител и Т-клеточных рецепторов.

1.7.2 Иммунопротеогеномика

Cheung и др. [33] первыми предложили новый иммунопротеогеномный способ выделения циркулирующих моноклональных антител (антител, произведенных клонами одной клетки) из сыворотки крови, позволяющий эффективно восстанавливать репертуар. Несмотря на то, что секвенирование моноклональных антител вскоре стало рутинным процессом [26; 34; 35], возможность секвениро-вания нескольких антител из смешанных проб (поликлональные антитела) стало значительным научным прорывом с обширным потенциалом применения в сфере биомедицины. Стало очевидно [33], что для того, чтобы иметь возможность определять как антитела взаимодействуют с определенными потенциально вредоносными агентами или антигенами (см. также [36-40]), необходимо провести анализ, использующий комбинацию СНП и МС технологий. В частности, Cheung и др. [33] показали, что наиболее представленные в репертуаре антитела, выделенные только с помощью СНП, могут иметь слабую представленность в данных МС, и тем самым иметь второстепенное значение с точки зрения биомедицины. Таким образом, иммунопротеогеномика имеет ключевую роль в развитии новых технологий для анализа антител. Несмотря на это, до 2015 года не существовало ни одного иммунопротеогеномного программного продукта. Репертуар антител оказался наиболее удобным типом базы данных для проведения масс-спектрометрического поиска. При этом, создание репертуаров антител также оказалось непростой задачей: гены антител не закодированы в зародышевой линии напрямую, а диверсифицированы соматическими рекомбинациями и мутациями [37]. Быстрое развитие технологий секвенирования сделало возможным полноразмерное сканирование адаптивных иммунных репертуаров и выявило новые задачи иммуноинформатики (см. последние обзоры: Georgiou и др. [38]; Robinson [41]; Yaari и Kleinstein [42]; Grieff и др. [43]). Высокая точность адаптивных иммунных репертуаров позволяет не только выполнять иммунопротеогеном-ный анализ, но и решать такие иммунологические вопросы как анализ клональ-ных линий [44;45], статистический анализ рекомбинаций и вторичной диверсификации [46; 47], анализ формирования иммунного ответа [32; 48], популяционный анализ локусов иммуноглобулина и TCR [49] (Рисунок 1.2).

УН1 УН2 У'Н.1 \'Н4 УН5 УН6

Эволюционный анализ

Л

Статистический анализ рекомбинаций

□ X 6

В клетки

Репертуар антител

Иммунопротеогеномика

Анализ эффекта множественных цепей

Динамика развития иммунного ответа

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

Рисунок 1.2 — Различные вычислительные задачи иммуногеномики.

1.8 Настоящая работа

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

1.8.1 Сборка высокополиморфных диплоидных геномов

Сборка высокополиморфных диплоидных геномов является сложной вычислительной задачей. Попытки сборки геномов с высоким уровнем полиморфизма с использованием стандартных подходов приводят к получению сильно фрагментированных контигов. В случаях со средним уровнем полиморфизма применение существующих алгоритмов сборки генома также приводит к неудовлетворительным результатам: принимая естественное разнообразие за ошибки в геномных прочтениях, они агрессивно удаляют значительную часть разнообразия аллелей. Таким образом, целью первой части данной работы стала разработка алгоритма для сборки высоко полиморфных геномов, способного строить длинные контиги и сохранять естественные вариации, и создание программного комплекса на его основе. Решением этой задачи стала разработка и реализация алгоритма dipSPAdes [1]. Для оценки качества работы dipSPAdes был предложен новый подход, использующий комбинацию высокополиморфных геномов гриба Schizophyllum commune. Анализ структуры генома S. commune позволил в дальнейшем выполнить высококачественную сборку популяции 24 грибов, которая послужила основой популяционного анализа [2].

Список литературы диссертационного исследования кандидат наук Сафонова, Яна Юрьевна, 2016 год

Список литературы

1. Safonova Y., Bankevich A., Pevzner P.A. dipSPAdes: Assembler for highly polymorphic diploid genomes // J Comp Biol. — 2015. — Vol. 22, no. 6. — Pp. 528-45.

2. Extraordinary genetic diversity in a wood decay mushroom / Maria A. Baranova, Maria D. Logacheva, Aleksey A. Penin et al. // Molecular Biology and Evolution.

— 2015. — Vol. 32, no. 10. — Pp. 2775-83.

3. IgRepertoireConstructor: A novel algorithm for antibody repertoire construction and immunoproteogenomics analysis / Y. Safonova, S. Bonissone, E. Kurpilyansky et al. // Bioinformatics. — 2015. — Vol. 31, no. 12. — Pp. i53-61.

4. Safonova Y., Lapidus A., Lill J. IgSimulator: A versatile immunosequencing simulator// Bioinformatics. — 2015. — Vol. 31, no. 19. — Pp. 3213-5.

5. Exploring the Mycoplasma capricolum genome: a minimal cell reveals its physiology / P. Bork, C. Ouzounis, G. Casari et al. // Molecular Microbiology. — 1995.

— Vol. 16, no. 5. — Pp. 955-67.

6. Life with 6000 Genes / A. Goffeau, B. G. Barrell, H. Bussey et al. // Science. — 1995. — Vol. 274, no. 546-6.

7. The complete genome sequence of Escherichia coli K-12 / F. R. Blattner, G. 3rd Plunkett, C. A. Bloch et al. // Science. — 1997. — Vol. 5331, no. 277.

— Pp. 1453-62.

8. C. elegans Sequencing Consortium. Genome sequence of the nematode C. elegans: a platform for investigating biology // Science. — 1998. — Vol. 5396, no. 282. — Pp. 2012-8.

9. The genome sequence of Drosophila melanogaster / K.A. Ketchum, R.A. Hoskins, X. Wang, et al // Science. — 2000. — Vol. 287, no. 5461. — Pp. 2185-95.

10. A whole-genome assembly of Drosophila / E. W. Myers, G. G. Sutton, A. L. Delch-er et al. // Science. — 2000. — Vol. 287, no. March. — Pp. 2196-205.

11. Initial sequencing and comparative analysis of the mouse genome / R. H. Waterston, K. Lindblad-Toh, E. Birney, et al // Nature. — 2002. — Vol. 420, no. 6915. — Pp. 520-62.

12. Green P. Documentation for PHRAP. — 1994. — URL: http://www.genome. washington.edu/UWGC/analysis-tools/phrap.htm.

13. Bonfield J. K., Smith K. F., Staden R. A new DNA sequence assembly program // Nucleic acids research. — 1995. — Vol. 23, no. 24. — Pp. 4992-99.

14. TIGR assembler: a new tool for assembling large shotgun sequencing projects / G. G. Sutton, O. White, M. D. Adams, A. R. Kerlavage // Genome Science and Technology. — 1995. — Vol. 1, no. 1. — Pp. 9-19.

15. Batzoglou S. ARACHNE: A whole-genome shotgun assembler // Genome Research. — 2002. — Vol. 12, no. 1. — Pp. 177-89.

16. Pevzner P. A., Tang H., Waterman M. S. An Eulerian path approach to DNA fragment assembly // PNAS. — 2001. — Vol. 98, no. 17. — Pp. 9748-53.

17. Zerbino D. R., Birney E. Velvet: Algorithms for de novo short read assembly using de Bruijn graphs // Genome Research. — 2008. — Vol. 18, no. 5. — Pp. 821-9.

18. ABySS: A parallel assembler for short read sequence data ABySS : A parallel assembler for short read sequence data / J. T. Simpson, K. Wong, S. D. Jackman et al. // Genome Res. — 2009. — Vol. 19. — Pp. 1117-23.

19. IDBA-UD: A de novo assembler for single-cell and metagenomic sequencing data with highly uneven depth / Y. Peng, H. C. M. Leung, S. M. Yiu, F. Y. L. Chin // Bioinformatics. — 2012. — Vol. 28, no. 11. — Pp. 1420-8.

20. SOAPdenovo2: An empirically improved memory-efficient short-read de novo assembler / R. Luo, B. Liu, Y. Xie et al. // GigaScience. — 2012. — Vol. 1, no. 1. — P. 18.

21. SPAdes: A New Genome Assembly Algorithm and Its Applications to Single-Cell Sequencing / A. Bankevich, S. Nurk, D. Antipov et al. // Journal of Computational Biology. — 2012. — Vol. 19. — Pp. 455-77.

22. Pham S. K., Pevzner P. A. DRIMM-Synteny: Decomposing genomes into evolutionary conserved segments // Bioinformatics. — 2010. — Vol. 26, no. 20. — Pp. 2509-16.

23. Bonissone S., Pevzner P. A. Immunoglobulin classification using the colored antibody graph//LectNotedComputSC, RECOMB 2015. — 2015. — Pp. 44-59.

24. Pevzner P. A., Tang H., Tesler G. De novo repeat classification and fragment assembly: from de Bruijn to A-Bruijn graphs // Genome Res. — 2004. — Vol. 14.

— Pp. 1786-96.

25. Trinity: Reconstructing a full-length transcriptome without a genome from RNA-Seq data / M. G. Grabherr, B. J. Haas, M. Yassour et al. // Nature Biotech. — 2013.

— Vol. 29, no. 7. — Pp. 644-52.

26. Automated de novo protein sequencing of monoclonal antibodies / N. Bandeira, V. Pham, P. Pevzner et al. // Nat Biotechnol. — 2008. — Vol. 26, no. 12. — Pp. 1336-8.

27. High-throughput sequencing of the zebrafish antibody repertoire / J.A. Weinstein, N. Jiang, R.A. White et al. // Science. — 2009. — Vol. 324. — Pp. 807-10.

28. Determinism and stochasticity during maturation of the zebrafish antibody repertoire / N. Jiang, J.A. Weinstein, L. Penland et al. // P Natl Acad Sci USA. — 2011.

— Vol. 108, no. 13. — Pp. 5348-53.

29. High-resolution description of antibody heavy-chain repertoires in humans / R. Ar-naout, W. Lee, P. Cahill et al. // PloS one. — 2011. — Vol. 6, no. 8. — P. e22365.

30. Lineage structure of the human antibody repertoire in response to influenza vaccination / N. Jiang, J. He, J. A. Weinstein et al. // Sci TranslMed. — 2013. — Vol. 5, no. 171. —P. 171ra19.

31. Genetic measurement of memory B-cell recall using antibody repertoire sequencing / C. Vollmers, R. V. Sit, J. A. Weinstein et al. // P Natl Acad Sci USA. — 2013.

— Vol. 110, no. 33. — Pp. 13463-8.

32. High-resolution antibody dynamics of vaccine-induced immune responses / U. Laserson, F. Vigneault, D. Gadala-Maria et al. // Proc Natl Acad Sci USA. — 2014. — Vol. 111. — Pp. 4928-33.

33. A proteomics approach for the identification and cloning of monoclonal antibodies from serum / W. C. Cheung, S. A. Beausoleil, X. Zhang et al. // Nat Biotechnol. — 2012. — Vol. 30, no. 5. — Pp. 447-52.

34. Automated protein (re)sequencing with MS/MS and a homologous database yields almost full coverage and accuracy / X. Liu, Y. Han, D. Yuen, B. Ma // Bioinformat-ics. — 2009. — Vol. 25, no. 17. — Pp. 2174-80.

35. Resurrection of a clinical antibody: Template proteogenomic de novo proteomic sequencing and reverse engineering of an anti-lymphotoxin-a antibody / N.E. Castellana, K. McCutcheon, V.C. Pham et al. // Proteomics. — 2011. — Vol. 11, no. 3.

— Pp. 395-405.

36. Proteomics-directed cloning of circulating antiviral human monoclonal antibodies / S. Sato, S.A. Beausoleil, L. Popova et al. // Nat Biotechnol. — 2012. — Vol. 30, no. 11. —Pp. 1039-43.

37. Molecular deconvolution of the monoclonal antibodies that comprise the polyclonal serum response / Y. Wine, D. R. Boutz, J. J. Lavinder et al. // P Natl Acad Sci USA.

— 2013. — Vol. 110, no. 8. — Pp. 2993-8.

38. The promise and challenge of high-throughput sequencing of the antibody repertoire / G. Georgiou, G. C. Ippolito, J. Beausang et al. // Nat Biotechnol. — 2014.

— Vol. 32, no. 2. — Pp. 158-68.

39. Identification and characterization of the constituent human serum antibodies elicited by vaccination / J. J. Lavinder, Y. Wine, C. Giesecke et al. // P Natl Acad Sci USA. — 2014. — Vol. 111, no. 6. — Pp. 2259-64.

40. Predicting immunogenic tumour mutations by combining mass spectrometry and exome sequencing / M. Yadav, S. Jhunjhunwala, Q. T. Phung et al. // Nature. — 2014. — Vol. 515, no. 7528. — Pp. 572-6.

41. Robinson W. H. Sequencing the functional antibody repertoire-diagnostic and therapeutic discovery // Nat Rev Rheumatol. — 2015. — Vol. 11, no. 3. — Pp. 171-82.

42. Yaari G., Kleinstein S. H. Practical guidelines for B-cell receptor repertoire sequencing analysis // Genome Med. — 2015. — Vol. 7, no. 121. — Pp. 1-14.

43. Bioinformatic and statistical analysis of adaptive immune repertoires / V. Greiff, E. Miho, U. Menzel, S. T. Reddy // Trends Immunol. — 2015. — Vol. 36, no. 11.

— Pp. 738-49.

44. IgTree©: creating immunoglobulin variable region gene lineage trees / M. Barak, N. S. Zuckerman, H. Edelman et al. // J Immunol Methods. — 2008. — Vol. 338, no. 1-2. — Pp. 67-74.

45. Change-O: a toolkit for analyzing large-scale B cell immunoglobulin repertoire sequencing data / N. T. Gupta, J. A. Vander Heiden, M. Uduman et al. // Bioinfor-matics. — 2015. — Vol. 31, no. 20. — Pp. 3356-8.

46. Statistical inference of the generation probability of T-cell receptors from sequence repertoires / A. Murugan, T. Mora, a. M. Walczak, C. G. Callan // PNAS. — 2012.

— Vol. 109, no. 40. — Pp. 16161-6.

47. Inferring processes underlying B-cell repertoire diversity / Y. Elhanati, Z. Sethna, Q. Marcou et al. // Philos Trans R Soc LondB Biol Sci. — 2015. — Vol. 1676, no. 370.— P. 20140243.

48. Next generation sequencing for TCR repertoire profiling: platform-specific features and correction algorithms / D.A. Bolotin, I.Z. Mamedov, O.V. Britanova et al. // Eur J Immunol. — 2012. — Vol. 42, no. 11. — Pp. 3073-83.

49. Automated analysis of high-throughput B-cell sequencing data reveals a high frequency of novel immunoglobulin V gene segment alleles / D. Gadala-Maria, G. Yaari, M. Uduman, S. H. Kleinstein // PNAS. — 2015. — Vol. 112, no. 8. — Pp. E862-70.

50. CompeauP. E. C., PevznerP. A., Tesler G. Howto apply de Bruijn graphs to genome assembly // Nature Biotech. — 2011. — Vol. 29, no. 11. — Pp. 987-91.

51. The genome sequence of the malaria mosquito Anopheles gambiae / R. A. Holt, S. Broder, G. M. Subramanian et al. // Science (New York, NY). — 2002. — Vol. 298, no. 5591. — Pp. 129-49.

52. Whole-genome shotgun assembly and analysis of the genome of Fugu rubripes / S. Aparicio, J. Chapman, E. Stupka et al. // Science. — 2002. — Vol. 297, no. 5585.— Pp. 1301-10.

53. The draft genome of Ciona intestinalis: insights into chordate and vertebrate origins / P. Dehal, Y. Satou, R. K. Campbell et al. // Science. — 2002. — Vol. 298, no. 5601. —Pp. 2157-67.

54. The diploid genome sequence of Candida albicans. / T. Jones, N. A. Federspiel, H. Chibana et al. // Proceedings of the National Academy of Sciences of the United States of America. — 2004. — Vol. 101, no. 19. — Pp. 7329-34.

55. Noncoding regulatory sequences of Ciona exhibit strong correspondence between evolutionary constraint and functional importance / D. S. Johnson, B. Davidson, C. D. Brown et al. // Genome Research. — 2004. — Vol. 14, no. 12. — Pp. 244856.

56. Assembly of polymorphic genomes: algorithms and application to Ciona savignyi / J. P. Vinson, D. B. Jaffe, K. O'Neill et al. // Genome Research. — 2005. — Vol. 15, no. 8.— Pp. 1127-35.

57. Genome sequence of the model mushroom Schizophyllum commune. / R. A. Ohm, J. F. de Jong, L. G. Lugones et al. // Nature biotechnology. — 2010. — Vol. 28, no. 9. — Pp. 957-63.

58. HaploMerger: reconstructing allelic relationships for polymorphic diploid genome assemblies. / S. Huang, Z. Chen, G. Huang et al. // Genome Res. — 2012. — Vol. 22. — Pp. 1581-8.

59. Mauve: Multiple alignment of conserved genomic sequence with rearrangements / A.C.E. Darling, B. Mau, F.R. Blattner et al. // Genome Res. — 2004. — Vol. 14.

— Pp. 1394-403.

60. Hozumi N., Tonegawa S. Evidence for somatic rearrangement of immunoglobulin genes coding for variable and constant regions // ProcNatl AcadSci USA. — 1976.

— Vol. 73. — Pp. 3628-32.

61. The joining of V and J gene segments creates antibody diversity / M. Weigert, R. Perry, D. Kelley et al. // Nature. — 1980. — Vol. 283. — Pp. 497-9.

62. Passenger transgenes reveal intrinsic specificity of the antibody hypermutation mechanism: clustering, polarity, and specific hot spots / A.G. Betz, C. Rada, R. Pan-nell et al. // Proc Natl Acad Sci USA. — 1993. — Vol. 90. — Pp. 2385-8.

63. IgBLAST: an immunoglobulin variable domain sequence analysis tool / Jian Ye, Ning Ma, Thomas L Madden, James M Ostell // Nucleic Acids Res. — 2013. — Vol. 41. — Pp. W34-40.

64. iHMMune-align: hidden Markov model-based alignment and identification of germline genes in rearranged immunoglobulin gene sequences / B. A. Gaëta, H. R. Malming, K. J. L. Jackson et al. // Bioinformatics. — 2007. — Vol. 23, no. 13. — Pp. 1580-7.

65. Comprehensive assessment of T-cell receptor beta-chain diversity in alphabeta T cells / H.S. Robins, P.V. Campregher, S.K. Srivastava et al. // Blood. — 2009. — Vol. 114, no. 19. — Pp. 4099-107.

66. Profiling the T-cell receptor beta-chain repertoire by massively parallel sequencing / J.D. Freeman, R.L. Warren, J.R. Webb et al. // Genome Res. — 2009. — Vol. 19, no. 10. — Pp. 1817-24.

67. Overlap and effective size of the human CD8+ T cell receptor repertoire / H.S. Robins, S.K. Srivastava, P.V. Campregher et al. // Sci Transl Med. — 2010.

— Vol. 2, no. 47. — Pp. 47-64.

68. Exhaustive T-cell repertoire sequencing of human peripheral blood samples reveals signatures of antigen selection and a directly measured repertoire size of at least 1 million clonotypes / R.L. Warren, J.D. Freeman, T. Zeng et al. // Genome Res. — 2011. — Vol. 21. — Pp. 790-7.

69. MiXCR: software for comprehensive adaptive immunity profiling / D. A. Bolotin, S. Poslavsky, I. Mitrophanov et al. // Nat Meth. — 2015. — Vol. 12, no. 5. — Pp. 380-1.

70. IMSEQ — a fast and error aware approach to immunogenetic sequence analysis / L. Kuchenbecker, M. Nienen, J. Hecht et al. // Bioinformatics. — 2015. — Vol. 31, no. 18. — Pp. 1-9.

71. Burnet F. M. A Modification of Jerne's Theory of Antibody Production using the Concept of Clonal Selection// CA Cancer J Clin. — 1976. — Pp. 119-21.

72. Error correction of high-throughput sequencing datasets with non-uniform coverage / P. Medvedev, E. Scott, B. Kakaradov, P. Pevzner // Bioinformatics. — 2011.

— Vol. 27, no. 13. — Pp. i137-41.

73. Nikolenko S.I., Korobeynikov A., Alekseyev M.A. BayesHammer: Bayesian clustering for error correction in single-cell sequencing // BMC Genomics. — 2013. — Vol. 14. — P. S7.

74. Garey M. R., Johnson D. S. Computers and intractability, a guide to the theory of NP-completeness. — W.H. Freeman and Company, New York, 1979.

75. Ben-Dor A., Shamir R., Yakhini Z. Clustering gene expression patterns // J Comp Biol. — 1999. — Vol. 6, no. 3-4. — Pp. 281-97.

76. Karypis G., Kumar V. A Fast and High Quality Multilevel Scheme for Partitioning Irregular Graphs // SIAM Journal on Scientific Computing. — 1999. — Vol. 20, no. 1. — P. 359-92.

77. Galinier P., Habib M., Paul C. Chordal graphs and their clique graphs // In WG '95. — Springer-Verlag, 1995. — Pp. 358-71.

78. Target-decoy approach and false discovery rate: when things may go wrong / N. Gupta, N. Bandeira, U. Keich, P. A. Pevzner // J Am Soc Mass Spectr. — 2011. — Vol. 22, no. 7. — Pp. 1111-20.

79. Kim S., Pevzner P. A. MS-GF+ makes progress towards a universal database search tool for proteomics // Nat Commun. — 2014. — Vol. 5. — P. 5277.

80. Measurement and clinical monitoring of human lymphocyte clonality by massively parallel VDJ pyrosequencing / S. D. Boyd, E. L. Marshall, J. D. Merker et al. // SciTransl Med. — 2009. — Vol. 1. — Pp. 12-23.

81. Towards error-free profiling of immune repertoires / M. Shugay, O.V. Britanova, E.M. Merzlyak et al. // Nat Methods. — 2014. — Vol. 11. — Pp. 653-5.

82. pIRS: Profile-based Illumina pair-end reads simulator / X. Hu, J. Yuan, Y. Shi et al. // Bioinformatics. — 2012. — Vol. 28, no. 11. — Pp. 1533-5.

83. ART: a next-generation sequencing read simulator / W. Huang, L. Li, J. R. Myers, G. T. Marth // Bioinformatics. — 2012. — Vol. 28, no. 4. — Pp. 593-94.

84. Grinder: A versatile amplicon and shotgun sequence simulator / F. E. Angly, D. Willner, F. Rohwer et al. // Nucleic Acids Research. — 2012. — Vol. 40, no. 12. — Pp. 1-8.

85. Detecting heterozygosity in shotgun genome assemblies: Lessons from obligately outcrossing nematodes / A. Barriere, S.P. Yang, E. Pekarek et al. // Genome Research. — 2009. — Vol. 19. — Pp. 470-80.

86. Price A. L., Eskin E., Pevzner P. A. Whole-genome analysis of Alu repeat elements reveals complex evolutionary history // Genome Res. — 2004. — Vol. 14, no. 11.

— Pp. 2245-52.

87. Cordaux R., Batzer M. A. The impact of retrotransposons on human genome evolution// Nat Rev Genet. — 2009. — Vol. 10, no. 10. — Pp. 691-703.

88. Next-generation sequencing and protein mass spectrometry for the comprehensive analysis of human cellular and serum antibody repertoires / J. J. Lavinder, A. P. Horton, G. Georgiou, G. C. Ippolito // Current Opinion in Chemical Biology.

— 2015. — Vol. 24. — Pp. 112-20.

89. Ultra-high-throughput sequencing of the immune receptor repertoire from millions of lymphocytes / J. R. McDaniel, B. J. DeKosky, H. Tanno et al. // Nat Protocols.

— 2016. — Vol. 11. — P. 429-42.

90. IMGT®, the international ImMunoGeneTics information system® / M.-P. Lefranc, V. Giudicelli, C. Ginestoux et al. // Nucleic Acids Res. — 2009. — Vol. 37, no. suppl 1. —Pp. D1006-12.

Список рисунков

1.1 Основные подходы к решению задачи сборки геномов....................15

1.2 Различные вычислительные задачи иммуногеномики......................18

2.1 Результаты сборки высокополиморфных диплоидных геномов.....23

2.2 Преимущество использования гаплоконтигов для сборки высокополиморфных диплоидных геномов................24

2.3 Алгоритм маскирования полиморфизмов в графе де Брюйна, построенного по диплоидному геному. ..................26

3.1 Структура антитела..............................33

3.2 Различные способы представления репертуара антител..........35

3.3 Процессы секвенирования и построения репертуара антител......37

3.4 Граф Хэмминга, построенный из данных иммуносеквенирования репертуара антител, и его декомпозиция на плотные подграфы. ........40

3.5 Построение репертуара антител на основе декомпозиции графа Хэмминга на плотные подграфы. ..........................................41

3.6 Результаты иммунопротеогеномного поиска................43

3.7 Схема работы инструмента IgSimulator...................46

4.1 Эволюционное развитие репертуара антител...............50

4.2 Образование новых сегментов в локусе иммуноглобулинов в результате неравный кроссинговера. ......................................52

Список таблиц

1 Сравнение параметров секвенирования с помощью метода Сэнгера и самых популярных технологий СНП........................................11

2 Совместная сборка гапломов двух особей Schizophyllum commune. ... 28

3 Сборка диплоидного генома Candida albicans...............28

4 Результаты отдельных сборок гапломов двух особей S. commune из Таблицы 2...................................29

SAINT PETERSBURG STATE UNIVERSITY

Printed as a manuscript

Safonova Yana Yuryevna

APPLICATION OF GRAPH MODELS TO BIOINFORMATICS ANALYSIS OF HIGHLY VARIABLE BIOLOGICAL

SEQUENCES

Specialty 05.13.18 — «Mathematical modeling, numerical methods and complexes of programs»

PhD thesis

Scientific advisor: PhD, professor Pevzner Pavel Arkadyevich

Saint Petersburg — 2016

Table of contents

P.

Introduction........................................................................4

Chapter 1. Introduction in bioinformatics......................................8

1.1 The emergence of sequencing technologies................................8

1.2 The Next Generation Sequencing (NGS) era..............................8

1.3 Omics......................................................................10

1.4 DNA sequencing and bioinformatics......................................10

1.5 Computational formulation of problems related to DNA analysis .... 11

1.6 Genome assembly..........................................................11

1.6.1 OLC approach ....................................................12

1.6.2 De Bruijn graph approach ........................................12

1.6.3 Applications of different assembly approaches....................13

1.7 Immunogenomics..........................................................14

1.7.1 Adaptive immune repertoires......................................14

1.7.2 Immunoproteogenomics ..........................................14

1.8 Present work................................................................15

1.8.1 Assembly of highly polymorphic diploid genomes................15

1.8.2 Constructing adaptive immune repertoires ........................17

Chapter 2. Assembly of highly polymorphic diploid genomes..................18

2.1 Genome assembly using de Bruijn graphs................................18

2.2 Highly polymorphic diploid genomes......................................18

2.3 Assembly of highly polymorphic diploid genomes........................19

2.4 dipSPAdes algorithm......................................................19

2.4.1 Motivation for dipSPAdes development..........................19

2.4.2 Algorithm for masking polymorphisms in de Bruijn graph ... 21

2.4.3 Construcing consensus and double contigs ........................22

2.5 dipSPAdes benchmarking ..................................................23

2.5.1 Simulation of highly realistic diploid genome using S.

commune genome.........................23

2.5.2 Benchmarking dipSPAdes against existing assembly tools ... 23

P.

2.5.3 Validation of diploid assembles using S. commune haplomes . . 25

2.5.4 Biological applications of dipSPAdes ..............25

Chapter 3. Problems of computational immunogenomics...........27

3.1 Adaptive immune system and antibody repertoire............27

3.2 Sequencing adaptive immune repertoires.................29

3.3 Repertoire construction problem and IgRepertoireConstructor tool ... 31 3.3.1 Representation of full-length antibody repertoire ........31

3.4 IgRepertoireConstructor algorithm ........................................32

3.4.1 Motivation for IgRepertoireConstructor development......32

3.4.2 Construction of antibody repertoire using Hamming graph ... 33

3.4.3 Immunoproteogenomics search ....................................34

3.4.4 Comparison of genomics and proteomics multiplicties .....37

3.5 IgSimulator ................................................................37

Chapter 4. Further development........................41

4.1 Assembly of highly polymorphic diploid genomes............41

4.2 Computational immunogenomics ..........................................41

4.2.1 Quality assessment of antibody repertoire construction algorithms ..........................................................42

4.2.2 Evolutionary analysis of antibody repertoires..........42

4.2.3 Analysis of multi-chain effect ..................43

4.3 Application of assembly of highly polymorphic diploid genomes and immunogenomics analysis to population studies of immunoglobulin loci 44

Conclusions.....................................46

List of abbverations and acronyms........................47

Glossary ......................................48

References..........................................................................53

Table of figures......................................................................62

Table of tables ......................................................................63

Introduction

Motivation. Emergence of DNA sequencing opened new horizons of modern biology and allowed one to transpose the biological analysis from wet lab experiments to the analysis of digital data. The following accumulation of biological data led to the active development of bioinformatics, a huge interdisciplinary field of science that covers a wide range of biological problems and uses computer science, mathematics, statistics and software engineering to solve them. Many bioinformatics studies are focused on analysis of natural diversity of biological sequences: DNA, RNA, proteins. High level of variability (typical for highly polymorphic diploid genomes or adaptive immune repertoires) turns such analysis into extremely challenging computational problem.

In this work, computational challenges of two bioinformatics problems — assembly of highly mutated genomes and analysis of adaptive immune repertoires — were addressed. In both of these cases the standard approaches are unable to construct yield adequate results due to an extreme level of diversity. Ultimately, both problems have the same goal: distinguishing between natural variations and sequencing errors. In order to solve these problems, dipSPAdes, IgRepertoireConstructor and IgSimula-tor tools were developed. Computational experiments proved efficiency of proposed tools. Proposed solutions are the first publicly available and efficient tools for solving considered problems.

The goals of this work cover

1. design and development of an algorithm for assembly of highly polymorphic diploid genome; its benchmarking on simulated and real datasets;

2. design and development of an algorithm for construction of adaptive immune repertoires using immunosequencing data and its immunoproteogenomics validation;

3. design and implementation of a versatile tool for simulation of adaptive immune repertoires; its application to benchmarking of repertoire construction algorithms.

Methods. For assembly of highly polymorphic diploid genomes, dipSPAdes, a first de Bruijn graph-based algorithm was proposed. In contrast to conventional assemblers, dipSPAdes takes advantage of diploidy and high level of polymorphism. For

dipSPAdes testing, highly realistic simulation of diploid genome using Schizophyllum commune was developed.

For constructing adaptive immune repertoires using immunosequencing data, IgRepertoireConstructor algorithm was proposed. IgRepertoireConstructor constructs Hamming graph and uses its clustering for repertoire construction. For repertoire calculation using Hamming graph, procedure of dense subgraph finding using minimum triangulation of graph was designed. Repertoires generated by IgRepertoireConstructor were used for immunoproteogenomics search, results of which outperformed results of recent studies. For simulation of adaptive immune repertoires and IgRepertoireConstructor testing, IgSimulator tool was designed and implemented.

All tools were implemented using C++ and Python languages.

The practical value and the novelty. This work has resulted in the following

tools:

- dipSPAdes: an assembler for highly polymorphic diploid genomes (available as a part ofSPAdes genome assembler: bioinf.spbau.ru/ru/spades).

- IgRepertoireConstructor: an algorithm for error correction of immunosequencing reads and immunoproteogenomics analysis (yana-safonova. github.io/ig_repertoire_constructor).

- IgSimulator: versatile immunosequencing simulator (yana-safonova. github.io/ig_simulator).

Proposed algorithms became a first publicly available solutions for the selected problems.

Work presentation. Results of this work were presented on seminars of the «Center for Algorithmic Biotechnology» (SPbSU) and the following conferences and meetings:

1. Y-Tools, a toolkit for construction and analysis of adaptive immune repertoires. Talk at RepSeq 2016. The Hague, Netherlands. September, 2016.

2. Big Data challenges in analysis of adaptive immune repertoires. Talk at the UCSD. San Diego, USA. April, 2016

3. New algorithmic challenges in analysis of adaptive immune repertoires. Talk at the Genentech. San Francisco, USA. April, 2016

4. A novel toolkit for analysis of adaptive immune repertoires using im-munosquencing and mass-spectra data. Poster on RECOMB 2016. Santa Monica, USA. April, 2016.

5. IgTools: a toolkit for analysis of antibody repertoire from im- munosequenc-ing and mass spectra data. Poster on Immunogenomics 2015. Huntsville, USA. September, 2016

6. IgRepertoireConstructor: a novel algorithm for antibody repertoire construction and immunoproteogenomics analysis. Talk at ISMB/ECCB 2015. Dublin, Ireland. July, 2015.

7. IgTools: a toolkit for construction of antibody repertoire using immunose-quencing data. Poster on ISMB/ECCB 2015. Dublin, Ireland. July, 2015.

8. dipSPAdes: assembler for highly polymorphic diploid genomes. Talk at RE-COMB 2014. Pittsburgh, USA. April, 2014.

9. New Frontiers of Genome Assembly with SPAdes 3.0. Poster on the JGI User Meeting 2014, Walnut Creek, USA. March, 2014.

10. dipSPAdes: assembler for highly polymorphic diploid genomes. Talk at the Joint Genome Institute. Walnut Creek, USA, 2014. March, 2014.

11. NGS assembly of highly polymorphic diploid genomes. Poster on RECOMB 2013. Beijing, China. April, 2013.

Publications. Results of this work were described in 4 papers [1-4] that were published in journals indexed in SCOPUS and Web of Science databases.

Personal contribution. In the paper [1], the author plays the main role in the design, implementation and benchmarking of dipSPAdes algorithm. The paper [2] is joint with the laboratoty of Evolutionary genomics at the Moscow State University. The author participated in genome assembly of S. commune population. In the paper [3], the author proposed and implemented IgRepertoireConstructor algorithm The author also participated in immunoproteogenomics experiments and discussions. In the paper [4], the author designed and implemented IgSimulator tool.

Structure of the work. The dissertation includes introduction, 4 chapters, and conclusions. The text of the dissertation is comprised of 63 pages, including 14 figures and 4 tables. References include 90 citations.

In the first chapter, the general aspects of bioinformatics and motivation of this work were described. Further, problems of highly polymorphic genome assembly and construction of adaptive immune repertoires were formulated.

In the second chapter, the complexity of highly polymorphic genome assembly was shown. Also, dipSPAdes algorithm was proposed. At the end of the chapter, results of computational experiments and dipSPAdes benchmarking were provided.

The third chapter is dedicated to problems of computational immunogenomics. At the start of the chapter, the adaptive immune system and its main components were described. Also, several formulations of repertoire construction problem were provided. Further, IgRepertoireConstructor algorithm and immunoproteogenomics search results were described. At the end of the chapter, description of IgSimulator tool was given.

In the fourth chapter, the proposed algorithms and their further development were discussed.

Chapter 1. Introduction in bioinformatics 1.1 The emergence of sequencing technologies

The appearance of DNA sequencing technologies opened up new horizons for biological research and contributes to the development of a number of new studies. The Human Genome Project (HGP) that was started in 1990 was the first world-wide attempt to sequence, reconstruct and analyze the entirety of the human DNA sequence. As a result of the HGP project, a nearly complete nucleotide sequence of the human genome was assembled and annotated. This data allowed biologists to understand complex human diseases related to genotyping of specific viruses and bacteria, identify cancer-specific mutations, assign personal treatment and predict its effects more accurately. Thus, DNA sequencing has proved itself to be a powerful method of genome analysis and became the starting point for numerous biological studies, ranging from personal medicine to evolution research. Since that time the genomes of many model organisms and important human pathogens have been sequenced and analyzed: bacteria Mycoplasma capricolum by Bork, et al. in 1995 [5]; Baker's yeast Saccharomyces cerevisiae by Goffeau, et al. in 1995 [6]; colon bacillus Escherichia coli by Blattner, et al in 1997 [7]; nematode Caenorhabditis elegans by C. elegans Sequencing Consortium in 1998 [8]; fruit fly Drosophila melanogaster by Ketchum, et al. [9] and by Myers, et al. [10] in 2000; house mouse Mus musculus by Waterson, et al. in 2002 [11]. One could say that the emergence of DNA sequencing technologies has revolutionized biological research and contributed to the development of such areas of study as: personal medicine; forensic applied sciences; alternative sources of energy, such as biofuels; agronomy; animal husbandry; bioprocessing; bioarcheology; anthropology; and evolution studies.

1.2 The Next Generation Sequencing (NGS) era

HGP and other early studies relied on the Sanger method, which produces long (up to 900 nt) and highly accurate reads (fragments of DNA). Unfortunately, Sanger

Table 1 — Comparison of the characteristics of the Sanger sequencing method and those of the most popular NGS technologies on the market.

RL (bp) — read length in bp, Rs / run — number of reads per run, Cost / 1 Mbp — cost per 1 Mbp, RT — running time.

Methods

RL (bp) Rs / run

ER

Cost /1 Mbp

RT

SOLiD

SOLEXA / II-

lumina

454

IonTorrent Sanger

50 + 50 1.2 - 1.4 Gbp A-T bias (< 0.1 %) $ 0.13 1 - 2 w

< 300 < 6 Gbp substitutions (0.1 %) $ 0.15 1 d

700 1 Mbp indels (1 %) $ 10 24 h

< 400 < 80 Mbp indels (1 %) $ 1 2 h

< 900 N/A substitutions (0.001 %) $ 2400 < 3 h

sequencing was almost completely manual and extremely expensive. A single Sanger run is capable of producing « 380 DNA samples in several hours within and costs up to $ 2400. Particularly, the total cost of sequencing the human genome in HGP, for example, amounted to nearly $ 3 billion. This made the usage of the Sanger method impractical for large scale sequencing projects.

In the mid 1990s, the period of rapid development of sequencing technologies was started. This period is called the Next Generation Sequencing (NGS) era. The emergence of NGS technologies significantly reduced the costs associated with the sequencing process and made massive parallel sequencing possible. As a result, it became possible to sequence large DNA and RNA samples at very low cost. There are several families of NGS technologies, including: SOLiD (sequencing by ligation), SOLEX-A/Illumina (sequencing by synthesis), 454 (pyrosequencing), IonTorrent (Ion semiconductor sequencing). Compared to the Sanger method, all of these technologies produce reads of very short length, have significantly lower costs and a much higher throughput (Table 1). The variety of NGS technologies and their availability significantly contribute to the continued development of genomics, a discipline that uses a combination of genome analysis and bioinformatics methods to processing sequencing data.

1.3 Omics

The emergence of NGS technologies significantly enhanced genome analysis and allowed scientists to solve a number of problems related to the composition of the genome. This in turn gave a huge boost to the development of various -omics fields: genomics (study of genomes), metagenomics (study of complex communities and environmental samples), transcriptomics (study of RNA molecules, their structure and functions), comparative genomics (study of genome structure and function across various biological species), and immunogenomics (study of immune system components). As the NGS technologies continue to evolve, new ways of analyzing biological sequences are becoming available. E.g., Illumina MiSeq technology, released in 2013, is able to perform full-length scanning of adaptive immune repertoires (set of circulating antibodies and T cell receptors), making deep and accurate analysis of the immune system possible. The ability to solve immunogenomics problems opens new possibilities for the development of personalized medicine, estimation of vaccine efficiency and design of drugs for autoimmune disorders, cancer and infection diseases. However, even modern sequencing technologies, such as Illumina MiSeq, are unable to ideally solve the problem of reconstruction of highly polymorphic diploid genome, whose mutation rate is compatible to diversity of adaptive immune repertoires. As a result, evolving sequencing technologies leads to increasing need of more elaborate and sophisticated algorithmic approaches for processing new biological data.

1.4 DNA sequencing and bioinformatics

During the HGP, two important computational problems were formulated: to reconstruct the raw human DNA sequence (genome assembly) and identify the genes encoded within (genome annotation). Both of them can be computationally formalized using observations on input data (sequencing reads or raw DNA sequence). Despite the fact that no one can annotate a project as well as a human expert biologist, dealing with such a vast amount of data would be extremely time consuming and require a lot of memory. Computer programs allow one to significantly increase the rate of data throughput, proving themselves to be a more adequate solution for the purpose of

genome sequencing projects. The development and design of algorithms for automatic genome assembly and genome annotation was the start of the active development of bioinformatics, a modern field of science, whose aim is to come up with new methods and pieces of software to gain a better understanding of biological data. Nowadays, bioinformatics has become a huge interdisciplinary field of study that covers a wide range of biological problems and uses computer science, mathematics, statistics and software engineering to solve them.

1.5 Computational formulation of problems related to DNA analysis

DNA and RNA molecules are composed of simple units called nucleotides. Each nucleotide of the DNA molecule is composed of the following nucleobases: adenine (A), cytosine (C), guanine (G) and thymine (T). In RNA, thymine is replaced with uracil (U). Thus, we can consider that a DNA/RNA molecule can be represented as a sequence of {A, C,G,T} letters (we consider that T is equivalent to U without loss of generality). The sequencing DNA/RNA molecule is a process during which the nu-cleotide sequence is calculated. The product of the sequencing process is a collection of reads, i.e., short fragments of the initial molecule. The reads themselves are also a sequence of {A, C, G, T} letters. In an ideal situation, the reads end up being exact substrings of the initial molecule. Unfortunately, sequencing methods are not sufficiently precise and introduce sequencing errors. For example, the Illumina sequencing method is known to be prone to random substitutions (average error rate is 1 %). We also consider that start positions of reads are distributed uniformly. DNA sequencing allows one to transpose the biological analysis from wet lab experiments to the analysis of digital data.

1.6 Genome assembly

Since sequencing reads represent very short fragments of the genome, in order to be able to proceed with further analysis (gene annotation, comparative analysis, population analysis) they need to be assembled into longer fragments. The process of such

genome reconstruction is called genome assembly. In order to be able to reconstruct a random string from short reads, one needs to find the extension of each of the reads and glue the overlapping portions into a single sequence. However, genome sequence is highly repetitive due to the duplication of genes and the presence of transposable elements (Figure 1.1a). This greatly complicates genome assembly. Typically, genome assembly does not yield a single string, but a set of contiguous fragments of the genome (contigs). Genome assembly algorithms attempt to compute long, non-overlapping and accurate contigs. There are several conventional approaches for assembling a genome: the overlap-layout-consensus approach (OLC) and the de Bruijn graph approach.

1.6.1 OLC approach

The OLC-based approaches were originally developed for the purpose of assembling Sanger reads. The OLC algorithms calculate the overlap graph, in which all vertices represent reads. Every two vertices that are adjacent in the graph correspond to a pair of overlapping reads. Contigs are then found by locating the non-branching paths in the constructed overlap graph (Figure 1.1b). The overlap-layout-consensus approach is used in such tools as: PHRAP [12], GAP [13], TIGR [14], ARACHNE [15], and CELERA [10]. In the worst case scenario, the number of edges in the overlap graph is quadratic (i.e., graph contains O(n2) edges, where n is a number of original reads) making this approach impractical for high throughput NGS reads.

1.6.2 De Bruijn graph approach

To assemble genomes from short and high throughput NGS reads algorithms using de Bruijn graph were proposed. Vertices of the de Bruijn graph are all possible k-mers (strings of length k) extracted from input reads. Two k-mers are adjacent on the de Bruijn graph if they are a prefix and a suffix of the k + 1-mer presented in the input reads (Figure 1.1c). Computation of contigs is equivalent to finding paths that are consistent to input reads in the constructed de Bruijn graph. The de Bruijn graph

R

B

R

C

R

D

D

-

---^ c

b) c)

a) Genome consists of four unique fragments A, B, C, and D separated by a repeat R. Sequencing reads corresponding to genome are shown under the genome. Some reads cover copies of repeat R, some reads cover unique regions A, B, C, and D. b) Overlap graph constructed from sequencing reads shown in Fig. a), c) De Bruijn graph corresponding to genome shown in Fig. a). Note that A>mer length does not exceed a length of the repeat R and length of unique regions A, B, C, and D.

Figure 1.1 — Genome assembly approaches.

approach was proposed in 2001 [16] and is widely used in such genome assemblers as: Velvet [17], ABySS [18], IDBA-UD [19], SOAPdenovo [20], SPAdes [21].

1.6.3 Applications of different assembly approaches

Both the OLC approach and the de Bruijn graph approaches can be applied to assembly of various types of sequencing technologies (e.g., paired-end reads, mate pairs, Pacbio reads) and sample types (e.g., single-cell, metagenomes, diploid genomes). However, the unique feature of the de Bruijn graph is its ability to comprehensively represent the repeat structure of the genome (Figure 1.1c). This makes the de Bruijn graph a powerful tool that can be used for the purpose of comparative analysis of genomes [22], V(D)J labeling of immunoglobulins and TCRs [23], de novo repeat classification [24], and transcriptome assembly [25].

1.7 Immunogenomics

1.7.1 Adaptive immune repertoires

The subjects of immunogenomics studies are collections of circulating antibodies and T-cell receptors (or so-called adaptive immune repertoires). An adaptive immune repertoire can be represented by a set of sequences (each of them is characterized a unique antibody or T-cell receptor) with corresponding multiplicities (multiplicity shows the number of occurrences of the antibody/T-cell receptor in the repertoire). Adaptive immune repertoires reflect the state of an organism that make them important subjects of various biomedical studies. For example, monitoring of immune repertoire dynamics is a standard step of drug development pipeline. However, despite the importance in biology and medicine, the analysis of antibody and T-cell receptor repertoires remains poorly studied problem in immunogenomics. Until 2009, the computational analysis of antibodies had been performed using proteomics techniques [26] and did not rely on the DNA sequencing technologies. Weinstein, et al. [27] were the first to demonstrate the power of DNA sequencing for the purpose of analyzing antibody repertoires and to open an "NGS era" in antibody analysis. While many other immunose-quencing studies have quickly followed in their footsteps [28-32], until 2012 there were no attempts to integrate NGS and MS approaches for antibody analysis. The absence of this integration (immunoproteogenomics) was a bottleneck, preventing the emergence of new approaches for analyzing large-scale and complex repertoires of antibodies and T-cell receptors.

1.7.2 Immunoproteogenomics

Cheung, et al. [33] pioneered a new immunoproteogenomics approach for identifying circulating monoclonal antibodies from serum that enables high-throughput antibody development. While sequencing purified monoclonal antibodies has now become routine [26;34;35], sequencing multiple antibodies from a complex sample (poly-clonal antibodies) represents a breakthrough with great biomedical potential. The im-

portant conclusion in [33] is that antibody analysis should combine NGS and MS to infer antibodies interacting with a specific antigen, i.e., foreign and potentially harmful agent (see also [36-40]). In particular, Cheung, et al. [33] showed that the most well-represented transcripts in the antibody repertoire (revealed by NGS alone) may not be the most biomedically relevant. Thus, immunoproteogenomics is the key ingredient of the emerging new technology for antibody analysis. However, until 2015 no publicly available immunoproteogenomics software was available. An antibody repertoire (rather than a set of all DNA reads as in previous immunoproteogenomics studies) represents a sensible choice of a database for the follow up MS/MS searches. However, the construction of an antibody repertoire is a difficult problem since antibody genes are not directly encoded in the germline, but are diversified by somatic recombination and mutations [37]. Rapid development of sequencing technologies enabled deep full-length scanning of adaptive immune repertoires and posed new immunoinfor-matics challenges (see recent reviews by Georgiou, et al. [38]; Robinson [41]; Yaari and Kleinstein [42]; Grieff, et al. [43]). High accuracy of immune repertoires allows one to perform immunoproteogenomics analysis and solve such immunological problems as: analysis of clonal lineages [44; 45], statistical analysis of recombination events and secondary diversification [46; 47], analysis of immune response development [32; 48], population analysis of immunoglobulin and TCR loci [49] (Figure 1.2).

1.8 Present work

This work describes solutions of two bioinformatics problems: assembly of highly polymorphic diploid genomes and construction of adaptive immune repertoires using immunosequencing data. Both of these problems have various biological and biomedical applications.

1.8.1 Assembly of highly polymorphic diploid genomes

Assembly of highly polymorphic diploid genomes is a complicated computational problem. Existing genome assembly approaches typically report extremely frag-

/A

VIII VI 12 VH3 VH4 VHS VH6

Clonal lineages

h

B cells

0 100 Statistical inference of recombination events

Antibody repertoire

* \

Immunoproteogenomics

Multi-chain effect analysis

Immune response dynamics

Antibody repertoires are starting points for solving for various immunogenomics problems: analysis of clonal lineages and analysis of somatic hypermutations, analysis of recombination events, population analysis of adaptive immune locus and detection of novel alleles, analysis of immune response development and treatment monitoring, and immunoproteogenomics analysis.

Figure 1.2 — Problems of computational immunogenomics.

mented contigs in case of high polymorphism rate. Also, in case of medium polymorphism rate, existing assemblers are very aggressive and remove a significant portion of the allelic diversity, mistaking natural variations for sequencing errors. Thus, the goal of the first part of this work was to develop a genome assembly approach that is able to compute long contigs, while retaining genome diversity. The dipSPAdes algorithm was designed and implemented [1] as a solution for the first problem. In order to benchmark dipSPAdes, a novel approach using a combination of highly polymorphic genomes of Schizophyllum commune fungi was proposed. The analysis of the structure of 24 S. commune genomes allowed me to assemble them with high accuracy. The resulting assemblies were then used as a base for population analysis in [2].

1.8.2 Constructing adaptive immune repertoires

Construction of adaptive immune repertoires has very similar goal: to compute a correct repertoire sequences and retain their diversity. Achieving this goal requires accurate investigation of the specific features of the datasets and distinguishing between natural variation and sequencing errors. In [3], the IgRepertoireConstructor tool for construction of antibody repertoires and immunoproteogenomics analysis was proposed. To validate produced repertoires, we capitalize on the fact that sequencing and MS data have a different nature of errors. Computational experiments on polyclonal repertoires showed that the proposed method improves upon the results of a recent immunoproteogenomics study on monoclonal repertoires [33] that present much more simple case. We also released IgSimulator [4], a versatile immunosequencing simulator. IgSimulator tool generates antibody repertoires along with immunosequencing libraries, thus providing immunosequencing datasets with known references. Computational experiments demonstrated that IgRepertoireConstructor reports highly accurate repertoires for both simulated and real data and thus improves upon the state-of-the-art methodologies in the immunogenomics field. Published tools raised interest of leading research groups and biotechnology companies and allowed our lab to establish new productive collaborations.

Chapter 2. Assembly of highly polymorphic diploid genomes

2.1 Genome assembly using de Bruijn graphs

Let DB (Genome, k) be the de Bruijn graph [50] of a genome Genome and its reverse complement Genome', where vertices and edges correspond to k-mers and k + 1-mers, respectively. Each chromosome in Genome and Genome' corresponds to a path in this graph; a set of these paths represents the genome traversal of the graph. In this work, we will work with condensed de Bruijn graphs [21], where each edge represents non-branching path in the de Bruijn graph. An edge in the condensed de Bruijn graph is assigned a length (in k + 1-mers) and the length of a path is the sum of its edge lengths. Let DB (Reads, k) be the condensed de Bruijn graph constructed from a set of reads Reads from Genome and their reverse complements. For simplicity, we first consider an ideal scenario with a full coverage of Genome by error-free reads. In this case, the graphs DB (Reads, k) and DB (Genome, k) coincide. In practice, sequencing reads are error-prone (Table 1). As a result, all existing genome assemblers that use the de Bruijn graph are able to handle sequencing errors in reads.

2.2 Highly polymorphic diploid genomes

A diploid genome Genome contains a double set of chromosomes and thus can be viewed as a combination of two similar genomes Genome\ U Genome2. Genome\ and Genome2 are also called haplomes. Typically, the differences between the hap-lomes (polymorphisms) are represented as a collection of SNPs and short indels. Given a pairwise alignment, we use percent identity (percent of matches among all columns of the alignment) to measure the degree of similarity between haplomes. Correspondingly, divergence is 100 minus percent identity. Divergence of human genome does not exceed 0.01 %. In this work, we address a challenging problem of assembling haplomes that differ from each other by 0.4 — 10 %. We call such genomes highly polymorphic (HP). A high rate of polymorphism rate is typical for almost all diploid genomes (e.g., sea squirts, fungus, plants, insects) except for mammalian genomes. As a result, as-

sembly of HP diploid genomes presents a special interest for genome analysis of model organisms [51-53] and pathogens [54], evolutionary and population analysis [55; 56], and search for new antibiotics [57].

2.3 Assembly of highly polymorphic diploid genomes

The standard assembly approaches fail to reconstruct individual haplomes in such HP genomes. Assembly of a diploid genome can result in two types of contigs: hap-locontigs (contigs representing subsequences of individual haplomes) and consensus contigs (contigs representing a consensus of both haplomes for the orthologous regions) (Figure 2.1). Consensus contigs do not adequately represent haplomes, but are rather a mosaic of segments from both haplomes. Thus, the alleles present in a consensus contig of each of the polymorphic sites are somewhat arbitrarily chosen from one of the haplomes. In practice, since some regions of the HP genome are less polymorphic than others, HP genome assembly using conventional assemblers generates a mixture of haplocontigs and consensus contigs while assembling HP genomes. Pairs of haplo-contigs representing both haplomes for the same genomic region are defined as double contigs (Figure 2.1).

2.4 dipSPAdes algorithm

2.4.1 Motivation for dipSPAdes development

We propose dipSPAdes, a new algorithm that takes advantage of the de Bruijn graph constructed by SPAdes assembler [21] to generate both consensus contigs and haplocontigs (Figure 2.1), as a new tool for assembling HP diploid genomes. dipSPAdes uses the de Bruijn graph to mask polymorphisms in contigs and to produce a more comprehensive representation of the genome by both consensus contigs and haplocontigs. First, dipSPAdes constructs a de Bruijn graph from reads corresponding to the HP diploid genome Genome = Genome\ U Genome2 and compute haplocontigs using the

Two approaches for assembling highly polymorphic (HP) diploid genomes (two haplomes are shown on top as solid and dotted segments): conventional assemblers (black arrows) and dipSPAdes (green arrows). Conventional assemblers generate haplocontigs from both haplomes that are shown in red and blue (colors of haplocontigs are unknown in practice). dipSPAdes uses the de Bruijn graph to generate consensus contigs by combining and extending haplocontigs. Afterward, dipSPAdes restores allelic relations by using alignment of haplocontigs against consensus contigs.

Figure 2.1 — Assembly of highly polymorphic diploid genomes.

SPAdes algorithm. Typically, the constructed haplocontigs are highly fragmented and do not adequately represent the results of genome assembly. However, these haplocontigs have a highly important common feature: the endpoints of the haplocontigs from different haplomes do not match since the breakpoints of the haplocontigs are often located at different positions on different haplomes (Figure 2.2). As a result, the areas of overlaps between these haplocontigs allow one to assemble them into longer sequences using the OLC approach. In practice, this approach will not work because the overlapping haplocontigs that came from different haplomes may be highly diverged. This prevents one from detecting overlaps between two haplocontigs. To address this problem, dipSPAdes uses a polymorphism masking algorithm that essentially suppresses the differences between overlapping haplocontigs.

Haplocontig 1

Haplome 1

Haplome 2

Haplocontig 2

Polymorphic regions of haplomes of the diploid genome are shown in red and blue. Fragment of the de Bruijn graph corresponding to them is shown is the center. Red and blue regions may contain repeats unique for a haplome due to high divergence of polymorphic regions. As a result, constructed haplo-contigs may break at the different positions of haplomes and overlap. For example, the blue haplocontig breaks before the red haplocontig on the example.

Figure 2.2 — dipSPAdes takes advantage of using haplocontigs for assembly of HP diploid genomes.

Analysis of the de Bruijn graph corresponding to HP diploid genome revealed that most polymorphisms look like a pair of alternative paths in the graph that share a start and end vertices (Figure 2.3). Unfortunately, such structures may also correspond to diverged genome repeats. To distinguish polymorphic regions from diverged repeats dipSPAdes uses information about similarity of nucleotide sequences corresponding to paths. We expect that divergence does not exceed polymorphism rate of the genome. Thus, to mask polymorphisms in the de Bruijn graph one needs to find all pairs of alternative paths and glue them pairwise. Polymorphism masking algorithm changes structure of the de Bruijn graph and significantly simplifies it.

2.4.2 Algorithm for masking polymorphisms in de Bruijn graph

A E

C)

Blue and red paths represent polymorphic regions corresponding to homologous positions on haplomes. a) shows structure of the Bruijn graph corresponding to HP diploid genome. Red and blue polymorphic regions correspond to a pair of alternative paths in the graph, b) shows intermediate step of the polymorphism masking algorithm that searches for masking direction and projects one of alternative paths to other one. In our case, polymorphism masking is directed from red to blue, and the red path will be masked. As a result, algorithm creates projections of intermediate vertices of the red path on the blue path, c) After polymorphism masking, red path is removed from the graph, and all edges incident to its vertices (A and E on our example) will be rerouted to created projections (B' and D') on the blue path.

Figure 2.3 — Algorithm for masking polymorphisms in the de Bruijn graph corresponding to HP

diploid genome.

2.4.3 Construcing consensus and double contigs

Haplocontigs constructed at the first step of the dipSPAdes algorithm ideally correspond to paths in the graph before polymorphism masking. So, to mask polymorphisms in haplocontigs one needs to find their paths in the de Bruijn graph with masked polymorphisms. We refer to the resulting contigs as masked haplocontigs and acknowledge that polymorphism masking of this type may produce a consensus sequence that belongs to neither Genomei nor Genomeo. After that dipSPAdes searches for overlaps in masked haplocontigs and extends them, thus, improving the quality of assembly. As a result, the algorithm reconstructs double contigs in both haplomes Genomei and Genomeo.

2.5 dipSPAdes benchmarking

Each attempt to sequence the HP diploid genomes runs into the same difficult problem of needing to figure out how to evaluate accuracy of the resulting assemblies. This question often remains unsolved (this applies both to previous studies conducted in the pre-NGS era and recent studies) since there is no gold standard for checking the validity of HP assemblies. To provide a first comprehensive benchmarking of HP assemblies, we took advantage of a unique dataset generated in the course of a recent massive effort to sequence 24 genomes of Schizophyllum commune conducted in Dr. Alexey Kondrashov laboratory at the Moscow State University.

2.5.1 Simulation of highly realistic diploid genome using S. commune genome

S. commune is a model organism (wood-degrading mushroom) whose genome is ideally suited for benchmarking HP genome assemblers. Genome of S. commune is diploid, however, its haplome can be extracted using sequencing mycelia from a single meiospore. The unique feature of the widely distributed S. commune is that the haplome of any two different organisms differ by 7 — 12 % even if they were collected on the same continent (and up to 25 % if collected on different continents). Thus, combining reads from two S. commune haplomes perfectly models an HP genome, yet allows one to do that, which so far has been not possible to achieve — test the quality of results produced by assembly algorithms for diploid genomes.

2.5.2 Benchmarking dipSPAdes against existing assembly tools

We have selected several tools for benchmarking purposes: HaploMerger [58], SPAdes [21], and Velvet [17]. dipSPAdes and HaploMerger generate consensus contigs, while SPAdes and Velvet report haplocontigs. In addition, we benchmarked SPAdes*, a modified version of SPAdes (that also includes aggressive polymorphism collapsing in graph, but, in contrast to dipSPAdes, does not include masking in haplocontigs). We

introduced SPAdes* in our benchmarking to illustrate the advantages of dipSPAdes as compared to standard assemblers run in the aggressive polymorphism collapsing mode. While dipSPAdes and HaploMerger can be applied to any set of haplocontigs, in this study we applied them to SPAdes haplocontigs. To generate results for SPAdes, we turned off the procedure of sequencing error removal and removed all edges with low coverage instead. This allows us to avoid errors of sample preparation. We expect Velvet and SPAdes to produce assemblies with total length close to double length of a haplome, and SPAdes*, dipSPAdes, and HaploMerger to produce assemblies with a total length close to the that of a haplome. Benchmarking on both simulated using two S. commune haplomes (Table 2) and real fungi (Table 3) datasets (with polymorphism rate varying from 0.4 to 10 percent) demonstrated that dipSPAdes significantly improves assemblies of HP genomes.

Table 2 — Joint assembly of two Schizophyllum commune individuals.

Average polymorphism rate is 10 %, estimated haplome size is 38.9 Mbp. HaploMerger failed to produce results on these haplocontigs since it typically requires haplocontigs with N50 exceeding tens of kbp. Columns «ETL» and «CTL» refer to expected total length and computed total length (Mbp), respectively.

Velvet SPAdes HaploMerger SPAdes* dipSPAdes

ETL (Mbp) 77.8 77.8 38.9 38.9 38.9

CTL (Mbp) 39.15 60.33 N/A 45.91 38.36

# contigs 34,406 26,820 N/A 5721 3764

Largest contig 37,580 44,596 N/A 231,443 239,371

N50 1219 3598 N/A 24,931 27,245

N75 761 1694 N/A 8477 11,330

Table 3 — Assembly of diploid genome of Candida albicans.

Average polymorphism rate is 0.4 %, estimated haplome size is 14.5 Mbp

Velvet SPAdes HaploMerger SPAdes* dipSPAdes

ETL (Mbp) 29.0 29.0 14.5 14.5 14.5

CTL (Mbp) 11.28 17.37 2.84 14.85 13.16

# contigs 6731 4007 337 1540 1119

Largest contig 34,870 112,388 92,126 116,985 116,985

N50 2276 8788 23,529 25,691 28,039

N75 1155 3300 8115 10,639 12,809

2.5.3 Validation of diploid assembles using S. commune haplomes

Haplomes used in Table 4 can be used for validation of diploid assemblies constructed using dipSPAdes algorithm. For the benchmarking purposes, S. commune individuals were assembled separately using SPAdes and Velvet genome assemblers (Table 4). Assembly of a single S. commune haplome is much more simple problem compared to assembly of diploid genome simulated from two haplomes. As a result, contigs generated by SPAdes h Velvet (Table 4) are longer than contigs computed by dipSPAdes (Table 2). Table 4 also showed that SPAdes outperformed Velvet, so SPAdes contigs were selected for further analysis. Contigs from separate S. commune assemblies were used as a reference for quality assessment of contigs correponding to simulated diploid genome. To find potential assembly errors, dipSPAdes contigs were aligned against SPAdes using Mauve tool [59] (Mauve is able to deal with high polymorphism rate). Analysis of computed alignments revealed that 150 dipSPAdes contigs (all contigs of length at least 50 kbp) are substrings of SPAdes contigs except for one case (probable assembly error). Proposed validation illustrates high accuracy of dipSPAdes algorithm.

Table 4 — Separate assembly of two S. commune individuals used for simulation in Table 2.

Individuali Individual

Velvet SPAdes Velvet SPAdes

Total length (Mbp) 35.37 36.47 35.37 36.53

# contigs 3886 1893 4472 1804

Largest contig 186,960 524,070 162,349 780,644

N50 23,509 98,988 16,744 104,280

N75 11,217 44,666 8656 46,746

2.5.4 Biological applications of dipSPAdes

The ability to understand the complexity of S. commune genome structure allows us to modify the assembly procedure to accurately assemble a single individual genome. More specifically, we modified the de Bruijn graph cleaning and contig extraction algorithms to accurately process large paralogs (duplicated genes) typical for S. commune.

An accurate assembly of a single individual genome enabled large-scale analysis of S. commune population [2] (12 from the United States and 12 from European Russia). Results presented in [2] showed that S. commune genomes have an extreme level of diversity caused by an unusually high effective population size and a very high mutation rate. This exceptionally high level of nucleotide diversity also leads to an extreme amino acid diversity of protein-coding genes. The combination of the high variability of S. commune with a small genome size and ease of cultivation makes it a promising model organism for population, quantitative, and evolutionary genetics.

Chapter 3. Problems of computational immunogenomics

3.1 Adaptive immune system and antibody repertoire

The adaptive immune system is a subsystem of the overall immune system that is composed of highly specialized cells and processes meant to eliminate or prevent pathogen growth. After an initial response adaptive immunity creates an immunological memory to the specific pathogen, which that leads to a fast response to future encounters with that pathogen. Unlike the innate immune system, the adaptive immune system is highly specific to the particular pathogens. Adaptive immunity can also provide long-lasting protection: someone who has recovered from measles is protected against this disease for the rest of his or her life. In other cases, however, it does not provide a lifetime worth of protection. An example of this case would be: for example, chickenpox. The adaptive system response destroys the invading pathogens and any toxic molecules they produce. Sometimes the adaptive system is unable to accurately identify foreign molecules, as a result we get such afflictions as hay fever, asthma or other allergies. Antigens are any substances that elicit an adaptive immune response. The cells that carry out the adaptive immune response are white blood cells known as lymphocytes. The two main broad classes, antibody responses and cell mediated response, are also provided by two different lymphocytes (B lymphocytes/B cells and T lymphocytes/T cells, respectively). In antibody responses, B cells are activated to secrete antibodies, proteins that are also known as immunoglobulins (antibody structure is shown in Figure 3.1). Antibodies travel through the bloodstream and bind to the foreign antigens causing them to deactivate. Further we will describe all the algorithms in terms of antibodies since they present the most complicated case due to the presence of somatic hypermutagenesis. However, all proposed algorithms can be extended for repertoires of T cell receptors (TCRs), antigen recognizing products of T cells. The set of all antibodies in the human organism is called an antibody repertoire. The diversity of the antibody repertoire is achieved via mechanisms for antibody generation and diversification: V(D)J recombination, heavy and light chains pairing and somatic hypermutagenesis [60-62]. As a result of clonal selection, antibodies have an extremely uneven distribution of multiplicities. For example, several highly abundant antibodies are typical for repertoires that were taken during the immune response.

An antibody or immunoglobulin is a protein produced by B lymphocyte or B cell. Antibody consists of two identical heavy chains (shown in light blue) and two identical light chains (shown in yellow). Immunoglobulin genes encoding the chains are not presented in germline genome, but are the result of the so-called V(D)J recombination of germline DNA. During the V(D)J recombination of heavy chains, the V (variable), D (diversity) and J (join) gene segments (V and J for light chain) are randomly chosen from the corresponding repetitive families and glued into a single DNA sequence. V(D)J recombination form primary diversity of the antibody repertoire. The resulting V(D)J sequence is joined with a constant region and thus forms a new immunoglobulin gene. The V(D)J region of the immunoglobulin gene is also called a variable region. Lengths of variable regions for heavy and light chains are & 400 nt and & 350 nt, respectively. Variable regions contain three antigen binding sites or complementarity-determining regions shown in orange and red. The complementarity-determining region 3 (CDR3) is the most variable part of immunoglobulin gene since it covers the area where V, D, and J gene segments join. After successfully binding with an antigen, the corresponding B cell undergoes clonal expansion and somatic hypermutagenesis. The last process randomly introduces somatic hypermutations (SHMs) into DNA sequence of immunoglobulin genes and thus somehow changes antibody affinity (SHMs are shown by green circles).

Figure 3.1 — Antibody structure.

3.2 Sequencing adaptive immune repertoires

Length of variable region of an antibody (the most interesting part of an antibody for bioinformatics analysis) varies from 350 to 400 nucleotides (Figure 3.1). Complexity of antibody structure is a result of V(D)J recombination that does not allow one to apply genome assembly algorithms for reconstruction of full-length antibody sequences in an antibody repertoire. Thus, in all immunogenomics studies, sequencing technologies allowing unambiguously link each read with antibody were used. However, until recently, there were few attempts to construct full-length antibody repertoires (full-length classification) due to limitation of available sequencing technologies. For example, 454 technology used in early studies produces reads covering entire variable region of antibodies, but has extremely high error rate (Table 1). This makes 454 technology inefficient for immunogenomics studies. Moreover, total number of reads produced in a single 454 run is too small for comprehensive repertoire analysis. As a result, 454 reads were mainly used for V(D)J classification: finding the closest V, D, and J segments from the database which combination produce variable sequences for a given antibody (Figure 3.2a). Later released Illumina technology had high accuracy, but produced short reads (75 —100 nt) that were not able to cover variable regions of antibodies. Illumina short reads were used for CDR3 classification. CDR3 classification is focused on short region of variable region of antibody: CDR3 or third complementarity determining region (Figure 3.1). CDR3 classification is the more detailed representation of antibody repertoire compared to V(D)J classification. However, CDR3 classification does not allow to distinguish between antibodies that share CDR3s, but differ in other positions of the variable region (Figure 3.2b). Construction of full-length repertoires in turn provides an ideal representation of variable regions of antibodies (Figure 3.2c) and differs from the well-studied VDJ classification [23; 47; 63; 64] and CRD3 classification [65-68] problems. In fact, VDJ classification, CDR3 classification and full-length classification are three different clustering problems with increasing granularity of partition into clusters (from very rough to ideal) and different biological applications (Figure 3.2d). Emergence of MiSeq Illumina sequencing machine in 2013 opened horizons of adaptive immune repertoires investigation using highly accurate 250 x 2 Illumina reads. Availability of Illumina MiSeq reads raised interest of bioinformaticians to the problem and in 2015 three new tools for solving the issue of repertoire construction were released: IgRepertoireConstructor [3], MiXCR [69] and IMSEQ [70].

cdr3

-b

•b

a)

b)

-B

C)

d)

a) V(D)J classification, b) CDR3 classification and c) full-length classification present three types of decomposition of immunosequencing reads with increasing granularity d). V(D)J classification provides a very rough decomposition and allows one to reveal general properties of the repertoire. CDR3 classification gives an insight into a number of non-related antibodies, but glues antibodies that differ by SHMs in V segments only. Full-length classification derives all unique antibodies from immunosequencing reads and thus provide a comprehensive representation of the entire repertoire.

Figure 3.2 — Various representations of an antibody repertoire.

3.3 Repertoire construction problem and IgRepertoireConstructor tool

If we view an antibody as a center of a cluster formed by reads derived from this antibody, then construction of a repertoire corresponds to a difficult clustering problem with many closely located centers so that the radius of a cluster may exceed the distance from one cluster to another one. Since the standard clustering techniques (like &-means clustering) are not applicable to such problems, we have designed IgRepertoireCon-structor, a novel algorithm for constructing antibody repertoires.

3.3.1 Representation of full-length antibody repertoire