Разработка вычислительных методов для анализа и визуализации эукариотических геномных сборок тема диссертации и автореферата по ВАК РФ 03.01.09, кандидат наук Михеенко Алла Александровна

• Михеенко Алла Александровна
• кандидат наук
• 2021

Оглавление диссертации кандидат наук Михеенко Алла Александровна

Введение

Глава 1. Оценка качества эукариотических геномных сборок

1.1 Введение

1.1.1 Секвенирование генома

1.1.2 Сборка генома

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

1.1.4 Оценка качества геномных сборок: подходы на основе референсного генома

1.1.5 Проблемы анализа больших эукариотических геномов

1.2 Методы

1.2.1 Базовые метрики QUAST

1.2.2 Выравнивание сборок на референсный геном

1.2.3 Выбор наилучшего набора выравниваний

1.2.4 Обработка мобильных генетических элементов

1.2.5 Показатели качества на основе &-меров

1.2.6 Концепция теоретически оптимальной сборки

1.3 Результаты

1.3.1 Производительность QUAST-LG

1.3.2 Сравнение сборок

Глава 2. Анализ и оценка качества сборок длинных тандемных

повторов

2.1 Введение

2.1.1 Классификация тандемных повторов

2.1.2 Центромеры и их структурная организация

2.1.3 Трудности в сборке сверхдлинных тандемных повторов

2.1.4 Проблемы в оценке сборок сверхдлинных тандемных повторов

2.2 Методы

Стр.

2.2.1 Модуль выравнивания ридов

2.2.2 Модуль полировки

2.2.3 Модуль оценки качества

2.3 Результаты

2.3.1 Результаты TandemMapper на симулированных данных

2.3.2 Результаты TandemQUAST на симулированных данных

2.3.3 Анализ новых сборок центромерных регионов хромосом человека

2.4 Заключение

Глава 3. Новые подходы к визуализации данных в геномных

исследованиях

3.1 Введение

3.2 Icarus: новый инструмент визуализации, предназначенный для оценки качества сборки

3.2.1 Визуализация выравниваний сборки на референсный геном

3.2.2 Определение подобия между сборками

3.2.3 Визуализация структуры сборок

3.3 Применения программы Icarus для анализа и оценки геномных сборок

3.3.1 Набор данных S.cerevisiae

3.3.2 Набор данных C.elegans

3.4 Assembly Graph Browser: программное обеспечение для визуализации графов сборки

3.4.1 Упрощение графов сборки

3.4.2 Варианты представления графов сборки

3.5 Применение AGB для анализа графов сборки

3.5.1 Набор данных S.cerevisiae

3.5.2 Набор данных C.elegans

3.5.3 Набор данных H.sapiens

Заключение

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

Стр.

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

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

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

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

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

Введение

Актуальность работы. Первый проект сборки генома человека был завершен почти двадцать лет назад [1]. Human Genome Project был международным исследовательским проектом, который потребовал объединения усилий нескольких исследовательских групп по всему миру, 13 лет и 3 миллиарда долларов. Всего несколько лет спустя, с появлением технологий секвенирования нового поколения (англ. Next Generation Sequencing, NGS) в 2007 году, количество времени и денег, необходимых для секвенирования человеческого генома, резко сократилось. Ученые смогли секвенировать и анализировать сотни, а затем и тысячи различных геномов. Однако ограничения технологий NGS серьезно препятствовали созданию полных безошибочных сборок больших эукариотических геномов.

Затем появление и быстрое распространение технологий секвенирования на основе длинных ридов от компаний Oxford Nanopore Technologies и Pacific Biosciences привели к значительному увеличению количества впервые секвени-рованных эукариотических геномов. Большинство алгоритмов, ранее разработанных для данных NGS, не были достаточно эффективны для работы с новыми данными. Наряду с новыми геномными сборщиками потребовались новые программы для анализа и валидации создаваемых сборок.

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

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

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

Основные задачи данной работы:

- разработка и программная реализация высокоэффективных алгоритмов анализа и оценки качества сборок больших эукариотических геномов;

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

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

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

- QUAST-LG: инструмент для оценки качества сборок больших эукариотических геномов(http://cab.spbu.ru/software/quast-lg/) и Icarus (http://cab.spbu.ru/software/icarus/): модуль визуализации для пакета QUAST (зарегистрирован в Роспатенте, регистрационный номер 2016619067);

- AGB: инструмент для интерактивной визуализации и анализа графов сборки (http://cab.spbu.ru/software/agb/).

- TandemTools: программный пакет для анализа длинных тандемных повторов (http://cab.spbu.ru/software/tandemtools/).

Разработанные вычислительные методы широко применяются специалистами в области сборки генома. Icarus входит в пакет QUAST начиная с версии 4.0.0, QUAST-LG входит в пакет QUAST начиная с версии 5.0.0 и имеет почти 50 000 загрузок1. Представленные в программах Icarus и AGB подходы к визуализации доказали свою полезность для разработчиков геномных сборщиков [2],

1 Общее количество загрузок всех версий пакета, начиная с версии 5.0.0 на https://sourceforge. net/p/quast (дата обращения 31 июля 2020 года) и https://anaconda.org/bioconda/quast (дата обращения 31 июля 2020 года).

а также для исследователей в различных проектах по сборке генома [6-9]. Стоит отметить, что TandemTools является на данный момент единственным специализированным методом анализа и обнаружения ошибок в сборках высокоповторя-ющихся областей [10], и используется для валидации и обнаружения ошибок в центромерных областях новых сборок генома человека, создаваемых консорциумом Telomere-to-Telomere (рукопись находится в стадии подготовки).

Методы. Исследования и разработки, проведенные в рамках этой работы, основывались на использовании:

- алгоритмов обработки строк (в частности, алгоритмов выравнивания последовательностей);

- теории графов;

- методов эффективного управления вычислительными ресурсами, в частности многопоточного программирования.

Все инструменты были реализованы с использованием языков программирования Python и JavaScript.

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

1. TandemTools: mapping long reads and assessing/improving assembly quality in extra-long tandem repeats. ISMB 2020, виртуальная конференция. Июль 2020 года.

2. Versatile genome assembly evaluation with QUAST-LG. ISMB 2018, Чикаго, США. Июль 2018 года.

3. Analysis and visualization of segmental duplications in mammalian genomes. BiATA 2018. Санкт-Петербург, Россия. Июль 2018 года.

4. Icarus: visualizer for de novo assembly evaluation. ECCB 2016, Гаага, Нидерланды. Сентябрь 2016 года.

Основные научные результаты, выносимые на защиту:

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

2. Разработаны и реализованы интерактивные инструменты для анализа и визуализации сборок генома (Icarus) и графов сборки (AGB).

3. Разработан программный пакет TandemTools для анализа длинных тан-демных повторов. Пакет включает в себя инструмент TandemMapper для выравнивания длинных ридов, подверженных ошибкам, на сборки длинных тандемных повторов, и инструмент TandemQUAST для анализа и оценки качества сборок.

Публикации. Основные результаты данной работы описаны в четырех статьях [10-13], опубликованных в журналах, индексируемых в базах данных Web of Science Core Collection и Scopus (все журналы в настоящее время ранжируются как Q1).

Личный вклад. Автору диссертационного исследования принадлежит ведущая роль в исследованиях, разработке, программной реализации и написании текста статей в работах, посвященных Icarus [11], AGB [13] и TandemTools [10]. В статье QUAST-LG [12] автор играла главную роль в разработке и реализации алгоритмов, тестировании программы на различных наборах данных, а также участвовала в подготовке текста статьи.

Объем и структура работы. Диссертация включает введение, три главы и заключение. Полный объем диссертации составляет 108 страниц, включая 31 рисунок и 10 таблиц. Список литературы содержит 124 наименования.

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

Вторая глава описывает новые алгоритмы обработки длинных тандемных повторов. Обсуждаются недавние достижения и сохраняющиеся проблемы в области сборки и анализа длинных тандемных повторов. Описан пакет TandemTools для выравнивания длинных ридов на сборки длинных тандемных повторов и анализа и улучшения их качества. Данная глава частично адаптирована из [10].

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

струмент AGB для визуализации графов сборок. Данная глава основана на материалах работ [11] и [13].

Глава 1. Оценка качества эукариотических геномных сборок

1.1 Введение 1.1.1 Секвенирование генома

Секвенирование ДНК - это процесс определения точной последовательности нуклеотидов в молекуле ДНК. Первая технология секвенирования, использованная для получения полной последовательности генома бактериофага ФХ174, была основана на методе, предложенном Фредериком Сэнгером в 1977 году [14]. Эта технология обладает высокой точностью, но вместе с тем чрезвычайно высокой стоимостью, что ограничивает сферу ее применения. В настоящее время секвенирование по методу Сэнгера используется в основном для секвенирования небольших фрагментов генома, отдельных генов или для проверки результатов современных технологий секвенирования.

Следующим этапом в развитии методов секвенирования стали технологии высокопроизводительного секвенирования или секвенирования нового поколения (NGS). Наиболее популярные технологии NGS включали Illumina/Solexa [15] (секвенирование синтезом), Ion Torrent [16] (полупроводниковое секвенирование), Roche 454 [17] (пиросеквенирование). Их главным преимуществом было значительное снижение стоимости и продолжительности процесса секвенирования. Появление технологий NGS привело к взрывному росту числа вновь секвениро-ванных геномов. Однако до недавнего времени рутинная генерация высококачественных геномных сборок организмов с относительно большим размером генома была практически невозможна из-за ограничений технологий NGS. Технологии NGS генерируют относительно короткие геномные фрагменты, называемые ридами. Наиболее часто используемая технология NGS, Illumina, производит риды не длиннее 350 нуклеотидов. Хотя Illumina все еще широко используется в се-квенировании для клинических и исследовательских целей, короткая длина генерируемых последовательностей приводит к разрывам в сборках, так как богатые повторами области, составляющие до 15% генома, обычно не могут быть разрешены с помощью коротких ридов и остаются плохо представленными [18].

Новая эра в развитии технологий секвенирования началась с появления технологий секвенирования длиннымиридами или секвенирования третьего поколения. Эти технологии представлены методом одномолекулярного секвенирования в реальном времени от компании Pacific Biosciences [19] и секвенирования с помощью нанопор от Oxford Nanopore Technologies [20]. В отличие от технологий NGS, секвенирование длинными ридами позволяет прочитывать фрагменты генома длиной до сотен тысяч нуклеотидов. Основным недостатком этих технологий по сравнению с NGS является их высокая частота ошибок (15-30% в зависимости от платформы). Однако недавно компанией Pacific Biosciences был сделан важный шаг в развитии этих технологий. Технология высоконадежного секвенирования одной молекулы (англ. high-fidelity sequencing, HiFi) генерирует значительно более длинные (до 10-20 кб) риды [21], чем технология Illumina, сохраняя при этом чрезвычайно низкую частоту ошибок (до 0,1%).

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

1.1.2 Сборка генома

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

Для решения задачи сборки генома были разработаны различные алгоритмические подходы. Первые сборщики генома, предназначенные для работы с результатами секвенирования по методу Сэнгера, такие как Celera Assembler [22], ARACHNE [23], были основаны на подходе перекрытия-компоновки-консенсуса (англ. overlap-layout-consensus, OLC). Этот подход состоит из следующих шагов:

попарное выравнивание всех ридов, обнаружение перекрытий между всеми парами ридов и построение графа перекрытий. В графе перекрытий вершины представляют собой риды, и две вершины соединяются ребром, если соответствующие риды перекрываются. Стоит отметить, что обнаружение перекрывающихся областей является трудоемкой задачей даже при использовании дополнительных эвристик [22]. Поэтому подход OLC нецелесообразно использовать с данными NGS, так как они содержат огромное количество коротких последовательностей. Тем не менее, метод OLC вновь обрел популярность в эпоху секвенирования длинными ридами. Современные сборщики на основе длинных ридов, такие как Canu [24], FALCON [25] и другие, широко адаптировали парадигму OLC. Новые алгоритмы улучшили поиск перекрытий с учетом высокой частоты ошибок в длинных ридах и существенно ускорили попарное сравнение всех ридов [26].

Альтернативный подход к сборке с использованием графа де Брюйна (DBG) был впервые предложен Певзнером в 2001 году [27]. Этот подход основан на построении DBG из всех k-меров (последовательностей длины k) из ридов. Вершины графа представлены k-мерами, и две вершины соединены ребром, если они являются префиксом и суффиксом k + 1-мера, представленного во входном наборе ридов. Этот алгоритм был использован для создания многих популярных сборщиков на основе коротких ридов, таких как Velvet [28], ABySS [29], SOAPdenovo [30], и SPAdes [31]. Кроме того, подход на основе DBG также был реализован в гибридных сборщиках [32], которые используют одновременно короткие и длинные риды, и даже в сборщиках, использующих только длинные риды [33].

1.1.3 Оценка качества геномных сборок: подходы, не использующие

референсный геном

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

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

Одним из наиболее важных критериев при последующем анализе геномных сборок является то, насколько хорошо собраны отдельные гены. Здесь особенно важно учитывать различия между прокариотическими и эукариотическими организмами. Основное различие между эукариотическими и прокариотическими генами состоит в том, что большинство эукариотических генов имеют прерывистую структуру и состоят из кодирующих областей - экзонов и некодирующих вставок - интронов. Длина экзонов обычно колеблется от 100 до 600 п.н., а интронов -от сотен до многих тысяч нуклеотидов. Интроны могут составлять до 75% длины гена. У прокариот несколько генов могут быть организованы на одном опероне: они расположены рядом и управляют ферментами, которые осуществляют последовательные или близкие реакции синтеза. Существуют различные инструменты для предсказания генов в последовательности ДНК, которые учитывают вышеупомянутые особенности [34-36].

Альтернативный подход заключается в поиске почти универсальных высококонсервативных генов, то есть набора генов, которые, как ожидается, находятся в одной копии (т.е., не происходит их дублирования или потери) в любом секве-нированном геноме из соответствующей филогенетической группы. Инструмент оценки BUSCO [37] использует свои собственные базы данных консервативных ортологичных генов для оценки относительной полноты вновь секвенированных геномов.

Другая большая группа подходов основана на анализе ридов, которые использовались для сборки генома. Алгоритм, реализованный в программе KAT [38], сравнивает распределение &-меров в сборке и в ридах. Программный продукт REAPR [39] идентифицирует ошибки сборки, используя информацию о выравнивании рида на сборку генома. В идеале покрытие генома ридами должно быть равномерным, поэтому резкое изменения глубины покрытия может указывать на ошибку. Данный подход сильно зависит от качества выравнивания и не подходит для оценки сборок с неоднородным покрытием (например, полученных при секвенировании одной клетки).

Кроме того, существует группа инструментов, которые оценивают качество сборки как некоторое число. Например, ALE [40] и LAP [41] выравнивают риды на сборку и предсказывают вероятность правильности сборки. Инструмент

PDR [42] нацелен на объединение трех основных характеристик генома - полноты (какая часть генома собрана), правильности (сколько ошибок содержит сборка) и непрерывности (сколько фрагментов состоит сборка и какой они длины) -в один балл. Хотя этот подход может быть полезен для быстрой оценки качества, он не дает подробной информации о количестве и расположении ошибок сборки. Кроме того, в большинстве проектов по сборке генома нелегко добиться высокого качества с точки зрения как правильности, так и непрерывности. Более консервативный подход приводит к более фрагментированной сборке, в которой каждый фрагмент содержит меньше ошибок. С другой стороны, менее фрагментирован-ная сборка с большей вероятностью будет содержать ошибки сборки. Детальный анализ качества сборки может быть полезен для выбора стратегии сборки в каждом конкретном случае в зависимости от цели геномного проекта.

1.1.4 Оценка качества геномных сборок: подходы на основе референсного

генома

В случае, когда референсный геном для изучаемого организма существует, может быть очень полезно сравнить черновую сборку генома с референсным геномом. Одним из наиболее популярных программных продуктов, использующих референсный подход для оценки качества сборки, является QUAST [43]. С момента своего выпуска в 2012 году QUAST стал широко использоваться в научном сообществе в качестве стандартного метода оценки качества сборки. Конвейер QUAST состоит из выравнивания сборок генома на референсный геном и вычисления различных метрик качества, оценивающих полноту, непрерывность и правильность сборки. Основные метрики QUAST описаны ниже в Разделе 1.2.1.

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

мобильных генетических элементов (МГЭ) существенно усложняют оценку качества сборки.

1.1.5 Проблемы анализа больших эукариотических геномов

Большинство алгоритмов для оценки качества геномной сборки, в том числе QUAST [43] и подходы на основе коротких ридов, такие как REAPR [39], ALE [40] и LAP [41], были разработаны в первую очередь для небольших геномов, собранных из данных NGS, и с трудом могут быть применимы к анализу новых сборок эукариотических геномов на основе длинных ридов из-за технических и алгоритмических ограничений.

Во-первых, большой размер генома (например, длина генома человека превышает длину генома E. coli на три порядка) требует реализации новых высокоэффективных алгоритмов. Во-вторых, эукариотические геномы содержат огромное количество МГЭ, которые приводят к появлению множества различий между референсным геномом и геномом изучаемого организма. Таким образом, чрезвычайно важно различать истинные ошибки сборки и различные вариации, в том числе вызванные МГЭ. Более того, качество референсного генома обычно недостижимо в черновых сборках генома, поскольку высококачественные рефе-ренсные последовательности обычно являются результатом совместных усилий нескольких исследовательских групп и комбинации различных технологий. Таким образом, полезной может быть попытка более реалистичной оценки потенциально достижимого качества сборки с учетом используемого набора ридов. Наконец, большие эукариотические геномы содержат много высокоповторяющихся областей, что следует учитывать при выравнивании сборок на референсный геном и обработке результатов.

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

1.2 Методы

QUAST-LG является значительным расширением программы QUAST [43], разработанным для эффективной и всесторонней оценки геномных сборок любой длины. Алгоритмы программы QUAST были значительно оптимизированы для быстрого анализа больших геномов. Конвейер QUAST-LG включает в себя следующие этапы: (1) расчет базовых метрик качества, (2) выравнивание сборок на референсный геном и поиск ошибок сборки, (3) расчет метрик на основе уникальных &-меров, (4) аннотация сборки, (5) создание отчетов. Все шаги, кроме (1) и (5), являются необязательными и зависят от предоставленных данных и выбора пользователя.

1.2.1 Базовые метрики QUAST

Базовые метрики качества были подробно описаны в [43] и могут быть рассчитаны без использования референсного генома.

- # контигов: общее количество контигов в сборке.

- Total length: общее количество нуклеотидных оснований в сборке.

- GC (%): общее количество нуклеотидов G и C в сборке, деленное на общую длину сборки.

- Nx метрики (x может варьироваться от 0 до 100, наиболее широко используемые значения - N50, N75 и N90). Nx - это длина контига, для которой все контиги этой или большей длины составляют не менее x% длины сборки.

- Lx метрики показывают общее количество контигов с длиной, большей или равной Nx.

1.2.2 Выравнивание сборок на референсный геном

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

Чтобы точно оценить полноту и правильность сборки, QUAST-LG выравнивает сборки на высококачественный референсный геном. Полнота сборки рассчитывается как отношение числа позиций в референсе, покрытых по меньшей мере одним контигом сборки, к общей длине референсного генома. Правильность сборки обычно характеризуется количеством ошибок сборки, называемых также мисассемблами. [43] определяет точку разрыва сборки как позицию в контиге, удовлетворяющую одному из следующих критериев:

- последовательности с левой и правой стороны от данной точки выравниваются на разные цепи ДНК (инверсия) или на разные хромосомы (транслокация);

- расстояние между выравниваниями левой и правой последовательностей в референсе превышает некоторый порог ExtThreshold или выравнивания левой и правой последовательностей в референсе перекрываются более чем на ExtThreshold (перемещение). По умолчанию ExtThreshold равен 1 кб в случае прокариотических геномов и 7 кб в случае эукариотических геномов. Обоснование выбора этого значения описано в Разделе 1.2.4. Перемещения, транслокации и инверсии относятся к значительным ми-сассемблам;

- расстояние между выравниваниями левой и правой последовательностей в референсе больше, чем LocalThreshold (значение по умолчанию = 85 п.н., может быть задано пользователем), но не превышает ExtThreshold (локальный мисассембл). Считается, что локальные мисассемблы менее негативно влияют на качество сборки, чем значительные, поэтому их количество считается отдельно.

Ошибки сборки используются при расчете некоторых дополнительных метрик качества, наиболее важными из которых являются NGAx и LGAx. Эти метрики аналогичны Шх и Lx, однако они основаны на длине референсного генома, а не на

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

Используя информацию о выравниваниях, QUAST также рассчитывает Реконструированную долю генома (%) (англ. genome fraction) как общее количество выровненных нуклеотидов в референсе, разделенное на длину референсного генома. Нуклеотид в референсе считается выровненным, если хотя бы один контиг имеет по меньшей мере одно выравнивание, проходящее через этот нуклеотид.

В случае больших геномов выравнивание сборки на референсный геном обычно является наиболее трудоемким шагом в конвейере QUAST-LG. В отличие от оригинального QUAST, QUAST-LG использует быстрый и эффективный инструмент выравнивания minimap2 [44]. Minimap2 основан на типичной процедуре seed-chain-align. Он собирает набор специальных k-меров (которые будут использоваться в качестве семян), называемых минимайзеры. Минимайзеры обычно выбираются следующим образом [45]: дана входная последовательность S и параметры k и w, выбирается самый маленький (в соответствии с предварительно определенным порядком) k-мер в каждом окне из w последовательных k-меров в S (такое окно содержит w + k—1 нуклеотидов). Использование минимайзеров основывается на следующем наблюдении: если две последовательности имеют общую подпоследовательность определенной длины, гарантируется, что они имеют общий минимайзер. Минимайзеры представляют собой лишь небольшую часть семян, поэтому они значительно ускоряют вычисление соответствия строк. Для каждой последовательности S minimap2 использует её минимайзеры в качестве семян, находит точные совпадения с минимайзерами референса и определяет наборы последовательных соответствий как цепочки.

Список литературы диссертационного исследования кандидат наук Михеенко Алла Александровна, 2021 год

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

1. Consortium International Human Genome Sequencing. Initial sequencing and analysis of the human genome, volume = 409, year = 2001 // Nature. — Pp. 860921.

2. Nurk S. at al. HiCanu: accurate assembly of segmental duplications, satellites, and allelic variants from high-fidelity long reads, bioRxiv. — URL: https://www. biorxiv.org/content/10.1101/2020.03.14.992248v3.

3. Bzikadze A., Pevzner P.A. Automated assembly of centromeres from ultra-long error-prone reads.

4. Miga K.H. at al. Telomere-to-telomere assembly of a complete human X chromosome. — URL: https://www.nature.com/articles/s41586-020-2547-7.

5. Miga K.H. Centromeric Satellite DNAs: Hidden Sequence Variation in the Human Population// Genes. — Vol. 10. — P. 352.

6. Paskova V. at al. Characterization of NDM-Encoding Plasmids From Enterobac-teriaceae Recovered From Czech Hospitals // Front. Microbiol. — Vol. 9. — Pp. 1-12.

7. Islam T. at al. Analysis of Subtelomeric REXTAL Assemblies Using QUAST //

IEEE/ACM Transactions on Computational Biology and Bioinformatics.

8. d'Avila Levy C.M. at al. First Draft Genome of the Trypanosomatid Herpetomonas muscarum ingenoplastis through MinlON Oxford Nanopore Technology and Illumina Sequencing // Trop. Med. Infect. Dis. — Vol. 5. — P. 25.

9. Mostajo N.F. at al. A comprehensive annotation and differential expression analysis of short and long non-coding RNAs in 16 bat genomes // NAR Genomics and Bioinformatics. — Vol. 2.

10. TandemTools: mapping long reads and assessing/improving assembly quality in extra-long tandem repeats / A. Mikheenko, A. Bzikadze, A. Gurevich at al. // Bioinformatics. — Vol. 36. — P. i75-i83.

11. Icarus: visualizer for de novo assembly evaluation / A. Mikheenko, G. Valin,

A. Prjibelski et al. // Bioinformatics. — 2016. — Nov. — Vol. 32, no. 21. — Pp. 3321-3323.

12. Versatile genome assembly evaluation with QUAST-LG / A. Mikheenko, A. Prjibelski, V. Saveliev et al. // Bioinformatics. — 2018. — June.

13. Mikheenko A., Kolmogorov M. Assembly Graph Browser: interactive visualization of assembly graphs // Bioinformatics. — Vol. 35. — P. 3476-3478.

14. Nucleotide sequence of bacteriophage phi X174 DNA / F. Sanger, G. M. Air,

B. G. Barrell et al. // Nature. — 1977. — Feb. — Vol. 265, no. 5596. — Pp. 687695.

15. Bennett S. Solexa Ltd // Pharmacogenomics. — 2004. — Jun. — Vol. 5, no. 4.

— Pp. 433-438.

16. Pennisi E. Genomics. Semiconductors inspire new sequencing technologies // Science. — 2010. — Mar. — Vol. 327, no. 5970. — P. 1190.

17. Marguiles M., et al. Genome sequencing in microfabricated high-density picolitre reactors // Nature. — 2005. — no. 437. — P. 376-380.

18. Consortium 1000 Genomes Project, et al. A global reference for human genetic variation//Nature. — 2015. — Vol. 526. — Pp. 68-74.

19. Real-time DNA sequencing from single polymerase molecules / J. Eid, A. Fehr, J. Gray et al. // Science. — 2009. — Jan. — Vol. 323, no. 5910. — Pp. 133-138.

20. Schneider G. F., Dekker C. DNA sequencing with nanopores // Nat. Biotechnol.

— 2012. — Apr. — Vol. 30, no. 4. — Pp. 326-328.

21. Wenger A.M., Peluso P. Rowell W.J., Chang P.C. Accurate circular consensus long-read sequencing improves variant detection and assembly of a human genome // Nat. Biotechnology. — Vol. 37. — Pp. 1155-1162.

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

23. ARACHNE: a whole-genome shotgun assembler / S. Batzoglou, D. B. Jaffe, K. Stanley et al. // Genome Res. — 2002. — Jan. — Vol. 12, no. 1. — Pp. 177-189.

24. Koren S. et al. Canu: scalable and accurate long-read assembly via adaptive k-mer weighting and repeat separation // Genome Res. — 2017. — 05. — Vol. 27, no. 5.

— Pp. 722-736.

25. Phased diploid genome assembly with single-molecule real-time sequencing / C. S. Chin, P. Peluso, F. J. Sedlazeck et al. // Nat. Methods. — 2016. — Dec.

— Vol. 13, no. 12. — Pp. 1050-1054.

26. Ruan J., Li H. Fast and accurate long-read assembly with wtdbg2 // Nat. Methods.

— Vol. 17. — P. 155-158 2020.

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

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

29. Simpson J.T. et al. ABySS: a parallel assembler for short read sequence data // Genome Res. — 2009. — Jun. — Vol. 19, no. 6. — Pp. 1117-1123.

30. Li R. et al. De novo assembly of human genomes with massively parallel short read sequencing // Genome Res. — 2010. — Feb. — Vol. 20, no. 2. — Pp. 265272.

31. 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.

32. hybridSPAdes: an algorithm for hybrid assembly of short and long reads / D. Antipov, A. Korobeynikov, J. S. McLean, P. A. Pevzner // Bioinformatics. — 2016.

— Apr. — Vol. 32, no. 7. — Pp. 1009-1015.

33. Assembly of long error-prone reads using de Bruijn graphs / Y. Lin, J. Yuan, M. Kolmogorov et al. // Proc. Natl. Acad. Sci. U.S.A. — 2016. — Dec. — Vol. 113, no. 52. — Pp. E8396-E8405.

34. Gene prediction in novel fungal genomes using an ab initio algorithm with unsupervised training. / Vardges Ter-Hovhannisyan, Alexandre Lomsadze, Yury O. Chernoff, Mark Borodovsky // Genome research. — 2008. — Dec.. —

Vol. 18, no. 12. —Pp. 1979-1990. —URL: http://dx.doi.org/10.1101/gr.081612. 108.

35. Zhu W. et al. Ab initio gene identification in metagenomic sequences // Nucleic Acids Res. — 2010. — Vol. 38, no. 12. — P. e132.

36. Tang S., Lomsadze A., Borodovsky M. Identification of protein coding regions in RNA transcripts // Nucleic Acids Research. — 2015. —Apr. — Vol. 43. — P. e78.

37. Simao F. et al. BUSCO: assessing genome assembly and annotation completeness with single-copy orthologs // Bioinformatics. — 2015. — Jun.

38. Mapleson D. et al. KAT: a K-mer analysis toolkit to quality control NGS datasets and genome assemblies // Bioinformatics. — 2017. — 02. — Vol. 33, no. 4. — Pp. 574-576.

39. HuntM. et al. REAPR: a universal tool for genome assembly evaluation // Genome Biol. — 2013. — Vol. 14, no. 5. — P. R47.

40. Clark S. et al. ALE: a generic assembly likelihood evaluation framework for assessing the accuracy of genome and metagenome assemblies // Bioinformatics. — 2013. — Vol. 29, no. 4. — Pp. 435-443.

41. Ghodsi M. et al. De novo likelihood-based measures for comparing genome assemblies // BMC research notes. — 2013. — Vol. 6, no. 1. — P. 334.

42. Xie Luyu, Limsoon Wong. PDR: a new genome assembly evaluation metric based on genetics concerns // Bioinformatics. — 2020.

43. Gurevich A. et al. QUAST: quality assessment tool for genome assemblies // Bioinformatics. — 2013. — Vol. 29, no. 8. — Pp. 1072-1075.

44. Li H. Minimap2: versatile pairwise alignment for nucleotide sequences // Bioinformatics. — Vol. 34. — P. 3094-3100.

45. Reducing storage requirements for biological sequence comparison / Michael Roberts, Wayne Hayes, Brian R. Hunt at al. // Bioinformatics. — Vol. 20. — P. 3363-3369.

46. Myers Gene, Miller Webb. Chaining Multiple-alignment Fragments in Sub-quadratic Time // Proceedings of the Sixth Annual ACM-SIAM Symposium on Discrete Algorithms. — SODA '95. — Philadelphia, PA, USA: Society for Industrial and Applied Mathematics, 1995. — Pp. 38-47. — URL: http://dl.acm.org/citation.cfm?id=313651.313661.

47. Kurtz S. et al. Versatile and open software for comparing large genomes // Genome Biol. — 2004. — Vol. 5, no. 2. — P. R12.

48. rnaQUAST: a quality assessment tool for de novo transcriptome assemblies / E. Bushmanova, D. Antipov, A. Lapidus et al. // Bioinformatics. — 2016. — Jul. — Vol. 32, no. 14. — Pp. 2210-2212.

49. Wessler S. R. Eukaryotic Transposable Elements and Genome Evolution Special Feature // Proc Natl Acad Sci USA. — 2006. — Vol. 103, no. 47. — P. 1760017601.

50. Genomic evolution of the long terminal repeat retrotransposons in hemiascomyce-tous yeasts / C. Neuveglise, H. Feldmann, E. Bon et al. // Genome Res. — 2002.

— Jun. — Vol. 12, no. 6. — Pp. 930-943.

51. Collins J. J., Anderson P. The Tc5 family of transposable elements in Caenorhab-ditis elegans // Genetics. — 1994. — Jul. — Vol. 137, no. 3. — Pp. 771-781.

52. Kaminker J. S. et al. The transposable elements of the Drosophila melanogaster euchromatin: a genomics perspective // Genome Biol. — 2002. — Vol. 3, no. 12.

— P. RESEARCH0084.

53. Prak E. T., Kazazian H. H. Mobile elements and the human genome // Nat Rev Genet. — 2000. — Vol. 1. — Pp. 134-144.

54. Myers J. S. et al. A comprehensive analysis of recently integrated human Ta L1 elements // Am J Hum Genet. — 2002. — Vol. 71. — Pp. 312-326.

55. Putnam N. H. et al. Chromosome-scale shotgun assembly using an in vitro method for long-range linkage // Genome Res. — 2016. — Mar. — Vol. 26, no. 3. — Pp. 342-350.

56. Meraculous2: fast accurate short-read assembly of large polymorphic genomes / J. A. Chapman, I. Y. Ho, E. Goltsman, D. S. Rokhsar // ArXiv e-prints. — 2016. — Aug..

57. Kokot M., Dlugosz M., Deorowicz S. KMC 3: counting and manipulating k-mer statistics //Bioinformatics. — 2017. — Sep. — Vol. 33, no. 17. — Pp. 2759-2761.

58. Bresler G., Bresler M., Tse D. Optimal assembly for high throughput shotgun sequencing // BMC Bioinformatics. — 2013. — Vol. 14 Suppl 5. — P. S18.

59. Li H. Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM // arXiv:1708.01492. — 2013.

60. Girgis H. Z. Red: an intelligent, rapid, accurate tool for detecting repeats de-novo on the genomic scale // BMC Bioinformatics. — 2015. — Jul. — Vol. 16. — P. 227.

61. Miller J. R., Koren S., Sutton G. Assembly algorithms for next-generation sequencing data // Genomics. — 2010. — Jun. — Vol. 95, no. 6. — Pp. 315-327.

62. Jackman S. D. et al. ABySS 2.0: resource-efficient assembly of large genomes using aBloom filter// Genome Res. — 2017. — 05. — Vol. 27, no. 5. — Pp. 768777.

63. Assembly of long, error-prone reads using repeat graphs / M. Kolmogorov, J. Yuan, Y. Lin, P.A. Pevzner // Nat. Biotechnology. — Vol. 37. — P. 540-546.

64. The MaSuRCA genome assembler / Aleksey V. Zimin, Guillaume Marcais, Daniela Puiu et al. // Bioinformatics. — 2013. — Vol. 29, no. 21. — Pp. 2669-2677. — URL: http://bioinformatics.oxfordjournals.org/content/29/21/ 2669.abstract.

65. Gnerre S. et al. High-quality draft assemblies of mammalian genomes from massively parallel sequence data // Proc Natl Acad Sci USA. — 2011. — Vol. 108, no. 4.— Pp. 1513-1518.

66. Chapman J. A. et al. Meraculous: de novo genome assembly with short paired-end reads // PLoS ONE. — 2011. — Vol. 6, no. 8. — P. e23501.

67. KajitaniR. et al. Efficient de novo assembly of highly heterozygous genomes from whole-genome shotgun short reads // Genome Res. — 2014. — Aug. — Vol. 24, no. 8. — Pp. 1384-1395.

68. 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.

69. Li H. Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences // Bioinformatics. — 2016. — Jul. — Vol. 32, no. 14. — Pp. 21032110.

70. Sahlin K. et al. BESST-efficient scaffolding of large fragmented assemblies // BMC Bioinformatics. — 2014. — Aug. — Vol. 15. — P. 281.

71. Assembling short reads from jumping libraries with large insert sizes /1. Vasi-linetc, A. D. Prjibelski, A. Gurevich et al. // Bioinformatics. — 2015. — Oct. — Vol. 31, no. 20. — Pp. 3262-3268.

72. BoetzerM., Pirovano W. SSPACE-LongRead: scaffolding bacterial draft genomes using long read sequence information // BMC Bioinformatics. — 2014. — Jun.

— Vol. 15.— P. 211.

73. Krzywinski M. et al. Circos: an information aesthetic for comparative genomics // Genome Res. — 2009. — Sep. — Vol. 19, no. 9. — Pp. 1639-1645.

74. Zhongming Z. et al. Investigating single nucleotide polymorphism (SNP) density in the human genome and its implications for molecular evolution // Gene. — 2003. — no. 312. — Pp. 207-213.

75. Bacolla A. at al. Abundance and length of simple repeats in vertebrate genomes are determined by their structural properties // Genome Res. — Vol. 18. — P. 15451553.

76. Yunis J.J., Yasmineh W.G. Heterochromatin, satellite DNA, and cell function. Structural DNA of eukaryotes may support and protect genes and aid in speci-ation // Science. — Vol. 174. — P. 1200-1209.

77. McFarland K.N. at al. SMRT Sequencing of Long Tandem Nucleotide Repeats in SCA10 Reveals Unique Insight of Repeat Expansion Structure // PLoS One. — Vol. 10.— P. 0135906.

78. Giunta S., Funabiki H. Integrity of the human centromere DNA repeats is protected by CENP-A, CENP-C, and CENP-T // Proc. Natl. Acad. Sci. USA. — Vol. 114. — P. 1928-1933.

79. Song J.H.T., Lowe C.B., Kingsley D.M. Characterization of a Human-Specific Tandem Repeat Associated with Bipolar Disorder and Schizophrenia // Am. J. Human Genetics. — Vol. 103. — P. 421-430.

80. Black E.M., Giunta S. Repetitive Fragile Sites: Centromere Satellite DNA As a Source of Genome Instability in Human Diseases // Genes. — Vol. 9. — P. 615.

81. Willems T., Gymrek M., Highnam G. The landscape of human STR variation // Genome Res. — Vol. 24. — P. 1894-1904.

82. Gymrek M. at al. Abundant contribution of short tandem repeats to gene expression variation in humans // Nat. Genetics. — Vol. 48. — P. 22-29.

83. Saini S. at al. Reference haplotype panel for genome-wide imputation of short tandem repeats // Nat. Commun. — Vol. 9. — P. 4397.

84. Manuelidis L., Wu J.C. Homology between human and simian repeated DNA // Nature. — Vol. 276. — P. 92-94.

85. Willard H.F., Waye J.S. Chromosome-specific subsets of human alpha satellite DNA: analysis of sequence divergence within and between chromosomal subsets and evidence for an ancestral pentameric repeat // J. Mol. Evol. — Vol. 25. — P. 207-214.

86. Willard H.F., Waye J.S. Hierarchical order in chromosome-specific human alpha satellite DNA // Trends in Genetics. — Vol. 3. — P. 192-198.

87. Hayden K.E. Human centromere genomics: now it's personal // Chrom Res. — Vol. 20.— P. 621-633.

88. Zimin A.V. at al. Hybrid assembly of the large and highly repetitive genome of Aegilops tauschii, a progenitor of bread wheat, with the MaSuRCA mega-reads algorithm // Genome Res. — Vol. 27. — P. 787-792.

89. Chin C.S. at al. Phased diploid genome assembly with single-molecule real-time sequencing //Nat. Methods. — Vol. 13. — P. 1050-1054.

90. Li H. Minimap and miniasm: fast mapping and de novo assembly for noisy long sequences // Bioinformatics. — Vol. 32. — P. 2103-2110.

91. Jain M. at al. Linear assembly of a human centromere on the Y chromosome // Nat. Biotechnology. — Vol. 36. — P. 321-323.

92. Genomic Characterization of Large Heterochromatic Gaps in the Human Genome Assembly / N. Altemose, K.H. Miga, M. Maggioni, H.F. Willard // PLoS Comput Biol. — Vol. 10. — P. e1003628.

93. Lomsadze A. et al. Gene identification in novel eukaryotic genomes by self-training algorithm // Nucleic Acids Res. — 2005. — Vol. 33, no. 20. — Pp. 64946506.

94. Marcais G. et al. MUMmer4: A fast and versatile genome alignment system // PLoS Comput. Biol. — 2018. — Jan. — Vol. 14, no. 1. — P. e1005944.

95. Fast Approximate Algorithm for Mapping Long Reads to Large Reference Databases / C. Jain, A. Dilthey, S. Koren, A. Phillippy // J. Comput. Biol. — Vol. 25. — P. 766-779.

96. Weighted minimizer sampling improves long read mapping / C. Jain, A. Rhie, H. Zhang at al. // bioRxiv.

97. Li H. Minimap2: versatile pairwise alignment for nucleotide sequences // Bioinformatics. — Vol. 34. — P. 3094-3100.

98. Chin C.S. at al. Nonhybrid, finished microbial genome assemblies from long-read SMRT sequencing data // Nat. Methods. — Vol. 10. — P. 563-569.

99. Loman N.J. at al. A complete bacterial genome assembled de novo using only nanopore sequencing data // Nat. Methods. — Vol. 12. — P. 733-735.

100. Lin Y. at al. Assembly of long error-prone reads using de Bruijn graphs // Proc. Natl. Acad. Sci. USA. — Vol. 113. — P. 8396-8405.

101. Fast and accurate de novo genome assembly from long uncorrected reads / R. Vaser, I. Sovic, N. Nagarajan, M. Sikic // Genome Res. — Vol. 27. — P. 737-746.

102. Dvorkina T., Bzikadze A.V., A. Pevzner P. The String Decomposition Problem and its Applications to Centromere Assembly, Bioinformatics. — in press).

103. HaafT., Willard H.F. Orangutan alpha-satellite monomers are closely related to the human consensus sequence // Mamm. Genome. — Vol. 9. — P. 440-447.

104. Hall S.E., Kettler G., Preuss D. Centromere satellites from Arabidopsis populations: maintenance of conserved and variable domains // Genome Res. — Vol. 13.

— P. 195-205.

105. NanoSim: nanopore sequence read simulator based on statistical characterization / C. Yang, J. Chu, R.L. Warren, I. Birol // Gigascience. — Vol. 6. — P. 1-6.

106. Robinson J. et al. Integrative genomics viewer // Nat. Biotechnol. — 2011. — Jan. — Vol. 29, no. 1. — Pp. 24-26.

107. Wick R. R et al. Bandage: interactive visualization of de novo genome assemblies // Bioinformatics. — 2015. — Oct. — Vol. 31, no. 20. — Pp. 3350-3352.

108. Krumsiek J., Arnold R., T. Rattei. Gepard: a rapid and sensitive tool for creating dotplots on genome scale // Bioinformatics. — Vol. 23. — P. 1026-1028.

109. SpeirM. etal. The UCSC Genome Browser database: 2016 update // Nucleic Acids Res. — 2016. — Jan. — Vol. 44, no. D1. — Pp. D717-725.

110. Buels R. at al. JBrowse: a dynamic web platform for genome visualization and analysis// Genome Biol. — Vol. 17.

111. Lee E. et al. Web Apollo: a web-based genomic annotation editing platform // Genome Biol. — 2013. — Vol. 14, no. 8. — P. R93.

112. Gordon D. et al. Consed: a graphical tool for sequence finishing // Genome Res.

— 1998. — Mar. — Vol. 8, no. 3. — Pp. 195-202.

113. Simpson J. T., Durbin R. Efficient construction of an assembly string graph using the FM-index // Bioinformatics. — 2010. — Jun. — Vol. 26, no. 12. — Pp. i367-373.

114. Schmid Michael et al. Pushing the limits of de novo genome assembly for complex prokaryotic genomes harboring very long, near identical repeats // Nucleic Acids Research. — 2018. — Vol. 46, no. 17. — Pp. 8953-8965.

115. Pu L. et al. Detection and analysis of ancient segmental duplications in mammalian genomes // Genome Res. — 2018. — 06. — Vol. 28, no. 6. — Pp. 901— 909.

116. Nielsen C. B. etal. ABySS-Explorer: Visualizing Genome Sequence Assemblies // IEEE Transactions on Visualization and Computer Graphics. — 2009. — Nov.

— Vol. 15, no. 6. — Pp. 881-888.

117. Boisvert Sébastien, Laviolette François, Corbeil Jacques. Ray: simultaneous assembly of reads from a mix of high-throughput sequencing technologies // Journal of computational biology. — 2010. — Vol. 17, no. 11. — Pp. 1519-1533.

118. Kunyavskaya Olga, Prjibelski Andrey D. SGTK: a toolkit for visualization and assessment of scaffold graphs // Bioinformatics. — 2018.

119. Nielsen C. et al. Visualizing genomes: techniques and challenges // Nat. Methods.

— 2010. — Mar. — Vol. 7, no. 3 Suppl. — Pp. S5-S15.

120. McIlwain S. J. et al. Genome Sequence and Analysis of a Stress-Tolerant, Wild-Derived Strain of Saccharomyces cerevisiae Used in Biofuels Research // G3: GENES, GENOMES, GENETICS. — 2016. — Vol. 6, no. 6. — Pp. 1757-1766.

121. Istace B. et al. de novo assembly and population genomic survey of natural yeast isolates with the Oxford Nanopore MinlON sequencer // GigaScience. — 2017.

— Vol. 6, no. 2.

122. Jaffe DB et al. The FASTG format specification (v1.00) 2012. — 2012. — URL: http://fastg.sourceforge.net/FASTG_Spec_v1.00.pdf.

123. Gansner Emden R., North Stephen C. An open graph visualization system and its applications to software engineering // Softw., Pract. Exper. — 2000. — Vol. 30, no. 11. — Pp. 1203-1233. — URL: https://doi.org/10.1002/1097-024X(200009) 30:11<1203::AID-SPE338>3.0.C0;2-N.

124. Karypis George, Kumar Vipin. Multilevel Algorithms for Multi-Constraint Graph Partitioning // Proceedings of the 1998 ACM/IEEE conference on Supercomputing. — IEEE Computer Society, 1998. — P. 28.

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

1.1 Поиск надежных выравниваний.......................22

1.2 Определение ложных ошибок сборки, вызванных МГЭ.........24

1.3 Количество ошибок сборки и возможных МГЭ при разных значениях

X........................................25

1.4 Конструкция теоретически оптимальной сборки.............28

1.5 Средство просмотра выравниваний в круговой форме для набора данных Human^p...............................35

2.1 Распределение частот уникальных и надежных 19-меров в сборке центромеры X.................................43

2.2 Пример надежных &-меров a and b в сборке A и риде R.........44

2.3 Пример сегмента цепи и расширенного сегмента цепи в сборке.....48

2.4 Различные типы уникальных надежных &-меров.............49

2.5 Графики покрытия и разрывов сборки...................56

2.6 Распределение различных типов уникальных надежных &-меров. ... 57

2.7 График метрики, основанной на мономерах................57

2.8 График распределения длин мономеров..................58

2.9 Метрика точек разрыва для сборок T2T4, T2T4polish, T2T7, centroFlye

и centroFlyePoHsh...............................61

2.10 Распределение различных типов уникальных надежных &-меров в сборках T2T4 и T2T4polish..........................62

2.11 Распределение длин мономеров.......................63

2.12 Графики парных выравниваний.......................65

3.1 Различные методы визуализации геномных сборок............67

3.2 Конвейер Icarus................................70

3.3 Окно просмотра выравниваний контигов....................................71

3.4 Исследование условия подобия.......................73

3.5 Окно просмотра размеров контигов..........................................75

3.6 Часть окна просмотра выравнивания контигов для сборок S.cerevisiae. 76

3.7 Главное меню Icarus для набора данных C.elegans............77

3.8 Часть окна просмотра выравнивания контигов для набора данных C.elegans, хромосома III...........................78

3.9 Пример вывода AGB............................. 80

3.10 Пример сжатого графа сборки, набор данных S.cerevisiae........84

3.11 Визуализация сборки C.elegans с помощью AGB, режим, основанный

на референсе.................................86

3.12 Визуализация сборки C.elegans с помощью AGB, режим, основанный

на контигах.................................. 87

3.13 Подграф большого графа сборки, набор данных H.sapiens........88

3.14 Визуализация сборки H.sapiens с помощью AGB, режим с фокусом

на повторах.................................. 89

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

1 Описание наборов данных, использованных для тестирования............31

2 Сборщики, использованные в исследовании................32

3 Производительность инструментов QUAST и QUAST-LG........32

4 Сравнение сборок двух протестированных наборов данных.......33

5 Тестирование инструментов TandemMapper, minimap2, и Winnowmap

на сборке DXZ1 del..............................55

6 Тестирование инструментов TandemMapper, minimap2, и Winnowmap

на сборке D6ZIsimulated...........................56

7 Распределение различных типов уникальных надежных &-меров в сборках T2T4, T2T4polis^, T2T7, centroFlye.................62

8 Количество разных типов юнитов в сборках T2T4, T2T7 и centroFlye

и наборе ридов................................64

9 Производительность Icarus на двух наборах данных...........75

10 Результаты поиска подобия между сборками в анализе набора

данных C.elegans...............................78

SAINT PETERSBURG STATE UNIVERSITY

Printed as a manuscript

Mikheenko Alla Alexandrovna

DEVELOPMENT OF COMPUTATIONAL METHODS FOR ANALYSIS AND VISUALIZATION OF EUKARYOTIC

GENOME ASSEMBLIES

03.01.09 — «Mathematical biology, bioinformatics»

Dissertation is submitted for the degree of Candidate of Physico-mathematical

Sciences

Translation from Russian.

Scientific advisor:

Candidate of Physico-mathematical Sciences, Professor

Pevzner Pavel Arkadyevich

Saint Petersburg — 2020

Table of contents

P.

Introduction....................................112

Chapter 1. Quality assessment of eukaryotic genome assemblies.......116

1.1 Introduction................................116

1.1.1 Genome sequencing.......................116

1.1.2 Genome assembly........................117

1.1.3 Quality assessment of genome assemblies: reference-free approaches............................118

1.1.4 Quality assessment of genome assemblies: reference-based approaches............................120

1.1.5 Challenges in analysis of large eukaryotic genomes.......120

1.2 Methods..................................121

1.2.1 Basic QUAST metrics......................121

1.2.2 Aligning assemblies to the reference genome..........122

1.2.3 Best set of alignments selection.................124

1.2.4 Handling transposable elements.................127

1.2.5 ^-mer-based quality metrics...................130

1.2.6 Upper bound assembly concept.................131

1.3 Results ................................... 133

1.3.1 QUAST-LG performance.....................135

1.3.2 Assemblies comparison ..................... 135

Chapter 2. Analysis and quality evaluation of long tandem repeat assemblies 140

2.1 Introduction................................140

2.1.1 Classification of tandem repeats.................140

2.1.2 Centromeres and their structural organization..........140

2.1.3 Challenges in assembly of ETRs.................141

2.1.4 Challenges in evaluation of ETR assemblies...........142

2.2 Methods..................................143

2.2.1 Read mapping module......................143

P.

2.2.2 Polishing module.........................147

2.2.3 Quality assessment module....................148

2.3 Results...................................154

2.3.1 TandemMapper results on simulated data............155

2.3.2 TandemQUAST results on simulated data............156

2.3.3 Analysis of novel human centromere assemblies ........159

2.4 Conclusion ................................164

Chapter 3. New approaches to data visualization in genome studies.....165

3.1 Introduction................................165

3.2 Icarus: a novel visualization software for assembly quality evaluation . 167

3.2.1 Visualization of assembly alignments to the reference genome . 168

3.2.2 Identifying similarity between assemblies............170

3.2.3 Visualization of contig structure of assemblies .........171

3.3 Icarus applications to assembly analysis and evaluation.........172

3.3.1 S.cerevisiae dataset........................173

3.3.2 C.elegans dataset.........................174

3.4 Assembly Graph Browser: software for visualization of assembly graphs 176

3.4.1 Simplification of assembly graphs................178

3.4.2 Modes of assembly graph representations............179

3.5 AGB applications to analysis of assembly graphs............180

3.5.1 S.cerevisiae dataset........................180

3.5.2 C.elegans dataset.........................180

Conclusion.....................................186

List of abbreviations and acronyms .......................188

Glossary ......................................189

References.....................................191

Table of figures...................................203

Table of tables ................................... 205

Introduction

Motivation. The first draft of the human genome assembly was completed nearly twenty years ago [1]. The Human Genome Project was an international research project that required joined efforts of multiple research groups over the world, 13 years, and $3 billion. Just a few years later, with the arrival of Next Generation Sequencing (NGS) technologies in 2007, the amount of time and money needed to sequence a human genome decreased dramatically. Scientists could sequence and analyze hundreds and then thousands of various genomes. However, the limitations of NGS technologies mostly prevented generation of complete error-free assemblies of large eukaryotic genomes.

Then, emergence and rapid spread of long-read sequencing technologies by Pacific Biosciences and Oxford Nanopore Technologies resulted in a highly increased amount of newly sequenced eukaryotic genomes. Most algorithms previously developed for NGS data were not efficient enough to work with long-read data. Alongside with novel genome assemblers, new programs for analysis and validation of generated assemblies were required.

One of the biggest challenges in a full reconstruction of genome sequences is still represented by repeat-rich regions that comprise a significant part of many eukaryotic genomes, including the human genome. Some of the largest repetitive regions in the human genome are represented by centromeric and pericentromeric regions. Only recently, novel genome assemblers capable of producing high-quality assemblies of repeat-rich regions, including centromeres, have been developed [2;3]. However, these assemblies cannot be validated using previously developed computational methods due to their extremely complex repetitive structure and great variability of repetitive regions among the individuals that do not allow to use a reference genome as a ground truth [4; 5].

It is particularly important in the analysis of large complex datasets to complement various quality metrics with data visualization. Efficient visualization techniques are proven to be helpful not only in presenting the known information but in revealing hidden patterns in the data. All existing visualization tools have their scope of application and certain limitations. Novel approaches to the visualization of various kinds of genomic data, including assembly graphs and genome assemblies, can be of great use

not only in various biological studies but also for the developers of assembly software to validate and improve their algorithms.

The main goals of this work include:

- design and implementation of highly effective algorithms for analysis and quality assessment of large eukaryotic genomes;

- development of computational methods for processing assemblies of long tandem repeats including accurate read mapping to repeat-rich regions, detecting assembly errors, and quality assessment of such assemblies;

- design and implementation of visualization methods for various kinds of genomic data, including genome assemblies and assembly graphs.

Practical value and novelty. This work has resulted in the following novel software tools:

- QUAST-LG: a novel tool for quality assessment of large eukaryotic assemblies (http://cab.spbu.ru/software/quast-lg/) and Icarus (http://cab.spbu.ru/software/icarus/): a visualization module for QUAST package (registered with Rospatent, the registration number is 2016619067);

- AGB: a tool for interactive visualization and analysis of assembly graphs (http://cab.spbu.ru/software/agb/).

- TandemTools: a software package for analysis of long tandem repeats (http: //cab.spbu.ru/software/tandemtools/).

The developed computational methods are widely accepted in the genome assembly field. Icarus is included in the QUAST package starting from version 4.0.0, QUAST-LG is included in the QUAST package starting from version 5.0.0 and has almost 50,000 downloads1. Approaches to visualization implemented in Icarus and AGB tools are proven to be useful for the developers of assembly software [2] and for researchers in various genome assembly projects [6-9]. It is worth to note that TandemTools is the only specialized method for analysis and error detection in assemblies of highly repetitive regions at the moment [10], and it is used for validation and error detection in centromeric regions of novel human genome assemblies generating by the Telomere-to-Telomere consortium (the manuscript is in preparation).

Methods. The research and development carried within this work were based on:

1The total download number of all package versions started from 5.0.0 on https://sourceforge.net/ p/quast (accessed on July 31, 2020) and https://anaconda.org/bioconda/quast (accessed on July 31,2020).

- string processing algorithms (especially sequence alignment algorithms);

- graph theory;

- techniques for efficient management of computational resources, in particular, multithreaded programming.

All tools were implemented using Python and JavaScript languages.

Work presentation. Results of this work were presented on seminars of the Center for Algorithmic Biotechnology (St. Petersburg State University) and the following international conferences either as an oral presentation or as a poster:

1. TandemTools: mapping long reads and assessing/improving assembly quality in extra-long tandem repeats. ISMB 2020, virtual conference. July 2020.

2. Versatile genome assembly evaluation with QUAST-LG. ISMB 2018, Chicago, USA. July 2018.

3. Analysis and visualization of segmental duplications in mammalian genomes. BiATA 2018. St. Petersburg, Russia. July 2018.

4. Icarus: visualizer for de novo assembly evaluation. ECCB 2016, Hague, Netherlands. September 2016.

Main results of the dissertation:

1. The QUAST-LG tool for quality assessment of large eukaryotic genome assemblies was developed. The tool is capable of evaluating various quality metrics of assemblies taking into account a complex structure of eukaryotic genomes. Efficient algorithms implemented in QUAST-LG allow to perform analysis of human genome assembly in a few hours.

2. The interactive tools for analysis and visualization of genome assemblies (Icarus) and assembly graphs (AGB) were designed and implemented.

3. The TandemTools software package for analysis of long tandem repeats was developed. The package includes TandemMapper tool for mapping long error-prone reads to assemblies of long tandem repeats and TandemQUAST tool for assembly analysis and validation.

Publications. Main results of this work are described in four papers [10-13] published in journals indexed in the Web of Science Core Collection and Scopus databases (all journals currently ranked as Q1).

Personal contribution. The author was the primary author responsible for the research, design, software implementation, and text writing in the Icarus [11], AGB [13], and TandemTools [10] papers. In the QUAST-LG [12] paper the author played the main role in designing and implementation of the algorithms, maintained the software, per-

formed testing, and benchmarking on various datasets, contributed to the preparation of the manuscript.

Structure of the work. The dissertation includes an introduction, three chapters, and a conclusion. The text of the dissertation is comprised of 97 pages, including 31 figures and 10 tables. References include 124 citations.

The first chapter contains the description of the genome assembly problem and various approaches to quality assessment of genome assemblies. The distinctive features of eukaryotic genomes in comparison to prokaryotes are discussed. The QUAST-LG tool developed for analysis and quality evaluation of large eukaryotic assemblies is described in detail. This chapter is in part adapted from [12].

The second chapter is devoted to novel algorithms for the analysis of extra-long tandem repeats (ETRs). Recent progress and remaining challenges in assembly and analysis of long tandem repeats are discussed. TandemTools package for mapping long reads to ETR assemblies and analyzing and improving their quality is described. This chapter is in part adapted from [10].

The third chapter is devoted to data visualization in genomic studies. The existing approaches, their fields of application, and limitations are discussed. The Icarus tool for visualization of genome assemblies and the AGB tool for visualization of assembly graphs are presented. This chapter is in part adapted from [11] and [13].

Chapter 1. Quality assessment of eukaryotic genome assemblies

1.1 Introduction 1.1.1 Genome sequencing

DNA sequencing is the process of determining the exact sequence of nucleotides within a DNA molecule. The first sequencing technology used to obtain a complete genome sequence of the QX174 bacteriophage was based on the method proposed by Frederick Sanger in 1977 [14]. This technology has high accuracy, but an extremely high cost that limited the scope of its application. Currently, Sanger sequencing is used mostly for sequencing small fragments of the genome, individual genes, or to validate the results of modern sequencing technologies.

The next stage in the development of sequencing methods was the technology of high-throughput sequencing or next-generation sequencing (NGS). The most widely used NGS technologies include Illumina/Solexa [15] (sequencing by synthesis), Ion Torrent [16] (semiconductor sequencing), Roche 454 [17] (pyrosequencing). Their main advantage was a significant reduction in the cost and duration of the sequencing process. The emergence of NGS technologies led to explosive growth in the number of newly sequenced genomes. However, until recently, the routine generation of high-quality genome assemblies of organisms with relatively large genome size was practically impossible due to the limitations of NGS technologies. NGS generates relatively short genomic fragments called reads. The most commonly used NGS technology, Illumina, produces reads not longer than 350 nucleotides. Although Illumina is still widely used in clinical and research sequencing, a short length of generated sequences results in gapped genome assemblies where highly repeat-rich regions that comprise up to 15% of the genome usually cannot be resolved using short reads and remain poorly represented [18].

The next era in the development of sequencing technologies began with the emergence of the long-read sequencing or third-generation sequencing technologies. These technologies are presented by the method of single-molecule sequencing in real-time from Pacific Biosciences [19] and nanopore sequencing from Oxford Nanopore Tech-

nologies [20]. Unlike NGS technologies, long-read sequencing allows to read genome fragments up to hundreds of thousands of nucleotides in length. The main disadvantage of these technologies in comparison with NGS is their high error rate (15% - 30% depending on the platform). However, an important step in the development of these technologies has been recently made by Pacific Biosciences. The Single Molecule High-Fidelity reads technology [21] (HiFi reads) generates significantly longer reads than Illumina technology (10-20 kbp length) while maintaining an extremely low error rate (up to 0.1%).

The rapid spread of long-read sequencing technologies led to a great increase in the amount of newly sequenced eukaryotic genomes. Progress in sequencing technologies set out algorithmic challenges in both genome assembly and downstream analysis.

1.1.2 Genome assembly

Genome assembly is a process of reconstruction of a genomic sequence from many shorter sequencing reads. Software products aiming to solve this problem are called genome assemblers or simply assemblers. Contiguous genome fragments reconstructed in the process of assembly are called contigs. Some assemblers may also order contigs separated by coverage gaps and join neighboring contigs into longer fragments called scaffolds.

Different algorithmic approaches have been developed to solve the problem of genome assembly. The first genome assemblers designed to work with Sanger reads, such as the Celera Assembler [22], ARACHNE [23], were based on the Overlap-Layout-Consensus (OLC) approach. This approach consists of the following steps: pairwise alignment of all reads, detection of overlaps between all pairs of reads, and constructing of an overlap graph. In the overlap graph, vertices represent reads and two vertices are connected with an edge if corresponding reads overlap. It is worth to note that the detection of overlapping regions is a time-consuming task even with additional heuristics [22]. Therefore, the OLC approach is impractical to use with the NGS data, since it contains a huge amount of short sequences. Nevertheless, the OLC method regained popularity within the era of long-read sequencing. Modern long-read assemblers, such as Canu [24], FALCON [25], and others, widely adopted the OLC

paradigm. New algorithms have improved an overlap detection to handle extreme error rates of long reads and substantially accelerated the all-vs-all read comparison [26].

An alternative assembly approach using a de Bruijn graph (DBG) was first proposed by Pevzner in 2001 [27]. This approach is based on the construction of a DBG from all k-mers from the reads (sequences of length k). The vertices of the graph are represented by k-mers, and two vertices are connected with an edge if they are a prefix and a suffix of the k + 1-mer presented in the input set of reads. This algorithm has been used to create many popular short-read genome assemblers such as Velvet [28], ABySS [29], SOAPdenovo [30], and SPAdes [31]. Moreover, the DBG-based approach has also been implemented in hybrid assemblers [32] that use simultaneously short reads and long reads and even assemblers using long-read data only [33].

1.1.3 Quality assessment of genome assemblies: reference-free approaches

A lot of different approaches to the evaluation of assembly quality were developed. Basically, they can be divided into two groups depending on using a known reference genome for the evaluation - the genome sequence of the studied (or closely related) organism that has been previously reconstructed with a high degree of accuracy. In genome assembly projects, the finished reference sequence is often absent and analysis is based on other sources of information. Basic characteristics of an assembly quality such as the number of contigs, the largest contig, the total assembly length, etc can be easily calculated, but they do not provide an estimation of assembly correctness and completeness and are insufficient for a comparison between different assemblies of the same genome.

One of the most important criteria in the downstream analysis of genome assemblies is how well individual genes are assembled. Here it is particularly important to take into account the differences between prokaryotic and eukaryotic organisms. The main difference between eukaryotic genes and prokaryotic genes is that most eukary-otic genes have a discontinuous structure and consist of coding regions - exons and non-coding insertions - introns. The exon length usually ranges from 100 to 600 bp, and of introns, from hundreds to many thousands of nucleotides. Introns can comprise up to 75% of a gene length. In prokaryotes, several genes can be organized on one operon: they are located nearby and control enzymes that carry out sequential or close

synthesis reactions. Various tools for predicting genes in the DNA sequence were developed to take into account the aforementioned features [34-36].

An alternative approach is to search only for near-universal highly conserved genes, i.e., the set of genes that are expected to be found in single-copy (without duplications or losses occurring) in any newly sequenced genome from the appropriate phylogenetic clade. The BUSCO assessment tool [37] uses its own databases of conservative orthologous genes to assess the relative completeness of newly sequenced genomes.

Another large group of approaches is based on the analysis of reads that were used to assemble the genome. The algorithm implemented in the KAT [38] program compares the distribution of &-mers in the assembly and the distribution of k-mers in reads. The REAPR [39] software product identifies assembly errors using information about read mapping against a genome assembly. Ideally, read coverage of the genome should be uniform, so a sharp drop or peak in coverage depth could indicate an error. This approach is highly dependent on the quality of the mapping and is not suitable for evaluating assemblies with non-uniform coverage (for example, obtained by single-cell sequencing).

In addition, there is a group of tools attempting to estimate an assembly quality in the form of a certain number. For example, ALE [40] and LAP [41] map reads to assemblies and predict the likelihood of assembly correctness. The PDR [42] tool aims to combine three main characteristics of the genome — completeness (what fraction of the genome is assembled), correctness (how many errors the assembly contains) and continuity (how many fragments the assembly consists of and how long they are) — into a single score. While this approach can be useful for a quick quality assessment, it does not provide any detailed information about the number and location of assembly errors. Moreover, in most genome assembly projects it is not easy to achieve high quality in terms of both correctness and continuity. A more conservative approach results in a more fragmented assembly, while each fragment contains fewer errors. On the other hand, a less fragmented assembly is more likely to contain assembly errors. A detailed analysis of the assembly quality can be useful to choose the best assembly strategy in each case depending on the purpose of the genomic project.

1.1.4 Quality assessment of genome assemblies: reference-based approaches

In the case when a reference genome for the studied organism exists, one may greatly benefit from the comparison of a draft genome assembly with the reference genome. One of the most popular software products using a reference-based approach for assembly quality evaluation is the QUAST tool [43]. Since its release in 2012, QUAST became widely used in the scientific community as a standard method of assembly quality evaluation. The QUAST pipeline consists of aligning genome assemblies to the reference genome and calculating various quality metrics representing assembly completeness, continuity, and correctness. Basic QUAST metrics are described below in Section 1.2.1.

However, the reference-based approach should be used with caution if a reference genome differs from the studied organism. If the reference genome matches exactly the dataset being collected, differences can be caused by errors in the assembly algorithm or sequencing errors such as chimeric reads. Otherwise, the differences can also be caused by true structural variations such as insertions, deletions, relocations, and so on. Such variations, as well as the presence of many transposable elements, significantly complicate the assessment of assembly quality.

1.1.5 Challenges in analysis of large eukaryotic genomes

Most of the algorithms for assembly quality evaluation, including QUAST [43] and short-read-based approaches such as REAPR [39], ALE [40], and LAP [41], were designed mostly for small bacterial genomes generated from NGS data and are hardly applicable to the analysis of novel eukaryotic genome assemblies based on long-read data due to technical and algorithmic limitations.

First, the large genome size (for example, the length of the human genome exceeds the length of E.coli genome by three orders of magnitude) requires the implementation of novel highly effective algorithms. Second, eukaryotic genomes contain a huge amount of transposable elements (TEs) that cause a lot of differences between the reference genome and the genome of the studied organism. Thus, it is extremely important to distinguish between true assembly errors and various discrepancies including those

caused by TEs. Moreover, the quality of a reference genome is usually unachievable in the draft genomes, because high-quality reference genome sequences are usually a result of joined efforts from multiple research groups and a combination of different technologies. So, one may benefit from a more realistic estimation of assembly quality with respect to the used set of reads. Finally, large eukaryotic genomes contain a lot of highly repetitive regions that should be taken into account when aligning assemblies to the reference genome and processing the results.

In the scope of this work, a novel evaluation tool QUAST-LG was developed to address the aforementioned concerns. This chapter describes in detail major algorithm modifications and improvements and demonstrates QUAST-LG performance on various datasets.

1.2 Methods

QUAST-LG is a significant extension of the QUAST tool [43] designed for efficient and versatile evaluation of genome assemblies of any length. The QUAST algorithms were greatly optimized for the fast analysis of large genomes. The pipeline of QUAST-LG includes the following steps: (1) calculation of basic quality metrics, (2) alignment of assemblies to a reference genome and detecting assembly errors, (3) calculation of metrics based on unique k-mers, (4) assembly annotation, (5) reports generation. All steps except for (1) and (5) are optional and depend on the provided data and the user choice.

1.2.1 Basic QUAST metrics

Basic reference-free metrics were described in details in [43] and can be calculated without a reference genome.

- # contigs: The total number of contigs in the assembly.

- Total length: The total number of nucleotide bases in the assembly.

- GC (%) is the total number of G and C nucleotides in the assembly, divided by the total length of the assembly.

- Nx metrics (x can range from 0 to 100, most widely used values are N50, N75, and N90). Nx is the contig length for which all contigs of this or greater length make up at least x% of the length of the assembly.

- Lx metrics show the total number of contigs with the length greater or equal to Nx.

1.2.2 Aligning assemblies to the reference genome

As was previously discussed, the three main characteristics of the assembly quality are completeness, correctness, and continuity. Basic continuity statistics, such as the N50 metric, can be assessed even in the absence of a reference genome.

To accurately measure completeness and correctness of the assembly, QUAST-LG aligns assemblies to a high-quality reference genome. The assembly completeness is calculated as the ratio between the number of the reference positions covered by at least one assembly contig and the total length of the reference genome. The assembly correctness is usually characterized by the number of large assembly errors called mis-assemblies. [43] defines a misassembly breakpoint as a position in an assembled contig that satisfies one of the following criteria:

- the flanking sequences align to opposite DNA strands (inversion) or to different chromosomes (translocation);

- the distance between the alignments of the left and right flanking sequences in the reference exceeds some threshold ExtThreshold or the alignments of the left and right flanking sequences in the reference overlap by over ExtThreshold (relocation). By default, ExtThreshold equals to 1 kbp in a case of prokaryotic genomes and 7 kbp in a case of eukaryotic genomes. The reason for selecting this threshold is described in Section 1.2.4. Relocations, inversions and translocations are called extensive misassemblies;

- the distance between the alignments of the left and right flanking sequences in the reference is larger than LocalThreshold (the default value = 85 bp, can be specified by a user) but does not exceed ExtThreshold (local misassembly). Local misassemblies are considered to be less harmful than extensive misas-semblies and reported separately.

Misassemblies are used to calculate some additional quality metrics, the most important of which are NGAx/LGAx. These metrics are similar to Nx/Lx, however, they are based on the length of the reference genome instead of the assembly length, and the contigs are first broken at each extensive misassembly position. The fragments not aligned to the reference genome are ignored.

Using information about alignments, QUAST also calculates Genome fraction (%) as the total number of aligned bases in the reference, divided by the length of the reference sequence. A base in the reference genome is considered as aligned if at least one contig has at least one alignment to this base.

In a case of large genomes, aligning the assembly against the reference genome is usually the most time-consuming step in the QUAST-LG pipeline. Unlike original QUAST, QUAST-LG uses the fast and efficient aligner tool minimap2 [44]. Minimap2 is based on a typical seed-chain-align procedure. It collects a set of special k-mers (to be used as seeds) called minimizers. Minimizers are usually selected as follows [45]: given an input sequence S and parameters k and w, select the smallest k-mer (according to a predefined ordering) in each window of w consecutive k-mers in S (such a window contains w + k-1 bases). Using the minimizer technique is based on the observation that if two sequences share a substring of a specified length, then they are guaranteed to have a matching minimizer. Minimizers represent only a small fraction of seeds so they dramatically speed up string-matching computations. For each query sequence, minimap2 takes query minimizers as seeds, finds exact matches to the reference, and identifies sets of colinear matches as chains.

Minimap2 offers a wide range of options greatly affecting the procedure of minimizer selection and chaining. In the default mode, QUAST-LG runs minimap2 with the parameters aiming at the maximal accuracy which is suitable for small genomes. However, these parameters are not optimal for large and complex genomes in terms of memory and time consumption. Moreover, using default options, minimap2 generates long but suboptimal alignments, allowing long gaps and stretches of mismatches in the alignment. Taking these drawbacks into account, QUAST-LG runs minimap2 in a way to produce more accurate, yet more fragmented alignments and, in addition, QUAST-LG implements its own post-processing procedure (described in details in Section 1.2.3) aimed to find an optimal set of alignments for each contig. After the rigorous benchmarking, the following important running parameters for minimap2 were chosen:

---no-long-join option to prevent extending an alignment through a long gap;

---min-occ-floor parameter is set to 200. This parameter forces minimap2 to

always use k-mers occurring --min-occ-floor times or less. Basically, a key property of the minimizer algorithm is that if two sequences share a substring, they are guaranteed to have a matching minimizer. However, minimap2 discards the most frequently occurring minimizers, so this guarantee is lost and mapping accuracy is reduced in repetitive regions. Increasing the maximal frequency improves minimap2' performance, although at the cost of increasing running time. An alternative approach to mapping long reads to highly repetitive regions is described in Chapter 3.

---mask-level parameter is set to 1 in the case of a metagenomic assembly and

0.9 otherwise. This parameter affects which chains will be marked as secondary and discarded. In metagenomic assembly, the same parts of an assembly tend to map simultaneously to several reference genomes, especially if these genomes have high similarity, e.g., represent different strains of the same bacterium. If --mask-level parameter is high, more chains can be extended and further reported by the aligner. If --mask-level is equal to 1, no chains are discarded due to an overlap with another chain, which is critical in a case of multiple similar alignments. —g parameter is set to 2,500 (the default value is 10,000). Minimap2 stops chain elongation if there are no minimizers within -g bp. However, if -g is small, minimap2 may elongate a chain through the large region without any minimizers and report long but erroneous alignment. Using a small -g forces minimap2 to report shorter alignments but with high identity.

1.2.3 Best set of alignments selection

It is unlikely for long contigs to map to the reference genome as a single perfect unambiguous alignment. An alignment software typically reports multiple alignment fragments A = (a\,... ,an) mapped to different positions in the genome, where some alignments may correspond to the same contig fragment mapped to distinct genomic positions. This can be caused by the presence of repetitive regions and transposable elements in the reference genome or by the errors in the assembly and alignment software. To provide accurate and complete mapping information for each contig, QUAST-LG se-

lects a set of non-overlapping alignments A С A that maximizes the total alignment score Score(A') (described below).

A general solution to the problem of how to find an optimal global chain of colinear non-overlapping fragments was proposed by Myers and Miller [46]. In practice, this problem is solved by alignment software using low-level chaining, i.e., joining short matching seeds into larger alignment fragments. For example, MUMmer [47] combines maximal unique matches and minimap2 [44] chains minimizers. In QUAST-LG, a dynamic programming algorithm called BestSetSelection was implemented for high-level chaining, i.e., combining larger alignment fragments into a set.

This algorithm is conceptually similar to the algorithm described in [48] for selecting the best spliced alignment for a mapped transcript. However, it takes into account various misassembly events, implements a set of different penalties, and includes a significant speed-up heuristic which allows to apply the algorithm for large genome assemblies. BestSetSelection is capable of correctly resolving large complex sets of alignments, which are common for eukaryotic genomes assemblies, and produces accurate chaining.

The BestSetSelection algorithm takes as an input the list of contig alignments A sorted according to the position in the contig of the right-most base. The scoring function takes into account the alignment score of each alignment and the locality score between adjacent alignments. Score of a single alignment a is defined as AlignScore(a) = Length(a) • Idy(a), where Idy(a) is reported by the alignment software and is calculated as the percent of matched bases between the alignment and the reference sequence, so Idy is equal to 100% for a perfect match. Given a scored set of alignments S = (ai,... ,ak—i), the score of S' = S U ak is computed as following:

Score(S') = Score(S) + AlignScore(ak) — P(ak-i,ak) — O(ak—1,ak) • Idy(ak) (1.1)

where P(ak—i ,ak) is a penalty and depends on the inconsistency between ak—i and ak in the genome, and O(ak—i,ak) is a number of overlapping bases in a contig between ak—i and ak. Higher penalty values are given for the extensive misassembly events and smaller values are used for local misassemblies (described in Section 1.2.2). The last term in the equation guarantees that the extension of the set with an alignment that fully overlaps in a contig with another alignment from the set is unprofitable.

The BestSetSelection algorithm has a time complexity of O(N2), where N is the number of alignments in A. Usually, this number is relatively small (up to 100), because short alignments are filtered out prior to running the algorithm. However, in the case

A contig (the grey line at the top) has multiple alignments (dark and light green and pink bars in the middle) to the reference genome (the grey line at the bottom), the positions of the alignment are visualized (with blue color in the contig and green/pink in the reference). Alignments A and C (dark green) are solid since they are sufficiently long and do not have large overlaps with other alignments. Alignment B is short, alignments E and F have significant overlaps, so their U niqueLength are not enough to mark them solid. Alignments and D2 are ambiguous and thus cannot be solid (their UniqueLength is equal to 0). Alignments for the best set are colored green (dark green for solid alignments and light green for the rest), the unused alignments are colored light pink. This figure is adapted from [12].

Figure 1.1 — Solid alignments detection.

of large eukaryotic genomes, the number of alignments per contig may reach dozens of thousands. In order to prevent the performance drop, a speed up heuristic procedure was implemented. It is used if the size of A exceeds MaxAlignments threshold (the default value is 100). This procedure guarantees to find the best set of alignments or one of the several sets that maximize the score.

To construct besti (the best alignment set ending at a«) the algorithm iterates through all already computed alignment sets S e BestSets and selects the one that gives the best score for S U ai. However, the distance between the majority of the alignment sets S e BestSets and ai is large, so Score(S U ai) is too small to add S U ai to BestSets. Thus, it is possible to apply a heuristic that allows to iterate only through a small subset of BestSets and significantly reduce the running time.

We define a as a solid alignment if Score(A'Ua) > Score(A') for any A' c A\a, i.e., if its addition to any subset of alignments from A \ a improves its Score. By the definition, all solid alignments from A are included in the resulting best set of alignments (it would be possible to create a set with a higher score otherwise). In particular, besti includes all solid alignments located to the left of ai in the contig. Thus, the algorithm can only iterate through those sets S e BestSets that include the right-most solid alignment before ai. Since BestSets is constructed iteratively, each of these sets will include all solid alignments to the left of ai. This significant reduction of the alignment sets to consider led to a tenfold drop in the running time on the large fragmented

assemblies that allowed us to complete in a reasonable time quality assessment for the human assemblies (see Section 1.3).

The criteria for choosing solid alignments can be deduced from equation (1.1). An alignment a is guaranteed to be solid if

UniqueLength(a) • Idy(a) > m • MaxPenalty (1.2)

where UniqueLength(a) is the length of a part of a that does not overlap with any other alignment, MaxPenalty is the maximal penalty value, m = 1 if a is located on the start/end of the contig and m = 2 otherwise. Figure 1.1 shows an example of the best alignment set and selection of solid alignments.

The described heuristic always identifies the alignment set with the maximal alignment score. Notably, its performance improvement depends on the number of solid alignments in the alignment set. If the input set contains very few solid alignments or, on the contrary, almost all alignments in the set are solid, this heuristic does not improve the running time of the algorithm. However, the benchmarking experiments on various real datasets demonstrated that such cases are almost impossible in practice and the heuristic works sufficiently well.

1.2.4 Handling transposable elements

Eukaryotic genomes usually contain a lot of transposable elements that have accumulated during the evolution [49]. These short variations cause a huge number of differences between the reference genome and the genome actually being assembled that can be reported as assembly errors according to the definition given above. To distinguish between true misassemblies and the ones caused by TEs, QUAST-LG performs an additional check of each relocation and inversion to identify possible TEs.

Figure 1.2 shows several examples where QUAST-LG successfully detects the presence of TEs. Let the threshold for maximal TEs size be X. Figure 1.2A shows the case when TE is present in the reference genome and missing in the assembly. QUAST-LG considered such a case as a relocation. However, if the inconsistency between alignments in the reference and the assembly is lower than X, QUAST-LG reports this case as a local misassembly. Figure 1.2B shows the case when the position of TE in the contig significantly differs from its position in the reference genome.

A

I TE 1 2

5

B

I 2 f 'SB

\ ^

M 2

R

R

C

1 I 2 —I

c

c

R

Ttb]

c

I I TE I

On each subfigure, we plot the reference genome R (top), the contig C (bottom), their matching fragments (blue and green bars for the positions in C and R, respectively) and TEs (violet bars). The inconsistencies in the alignments are shown by arrows and 5 characters. Subfigures demonstrate various cases of TEs presence in the assembly and in the reference genome and are described in details in the text. This figure is adapted from [12].

Original QUAST considers this situation as two misassembly breakpoints since 51 and 52 are usually much larger than X value. QUAST-LG here computes 5 as \S2 — 51\. 5 is equal to the TE length and should be lower than the appropriate X value, so this case is also considered as a possible TE. Figure 1.2C shows the situation where TE is located at the end of the contig, while its location in the reference is significantly away from the neighboring fragment. QUAST-LG cannot confidently consider this case as a TE based on the neighboring fragments as in the previous cases and reports it as a misassembly. An alternative approach would be to classify a TE based on its genomic sequence but is out of the scope of this work, as well as identification of tandem TE insertions and deletions.

The described procedure for TEs identification critically depends on the threshold for TEs size X. The optimal value of X should slightly exceed the length of the largest TE in the studied genome. A user can set this threshold manually if this information is known. Otherwise, we set the X value as 7 kbp by default. To select this value, we performed a following series of experiments.

We measured the number of misassemblies and the number of TEs identified by our procedure in all assemblies of the six benchmark datasets for various X in a range of 1-10 kbp with a step of 1 kbp. Figure 1.3 shows that the number of misassemblies drops significantly with an increase of X until a certain point X = Xc where the curve nearly flattens. Likewise, the number of possible TEs quickly increases from an almost

Figure 1.2 — Detection of false positive assembly errors caused by TEs.

The x-axis shows X value in kbp, the y-axis shows the number of misassemblies (solid lines) and possible TEs (dashed lines). Each line corresponds to one assembly of a corresponding dataset. Yeastpp and Yeast^p correspond to assemblies of S.cerevisiae genome generated using PacBio SMRT and Oxford Nanopore Technologies libraries correspondingly. The rest datasets are described in detail in Section 1.3. This figure is adapted from [12].

Figure 1.3 — Number of misassemblies and possible TEs at different breakpoint thresholds X.