Введение 3
1. Введение в предметную область 3
2. Процесс секвенирования 4
3. Форматы fasta и fastq 7
4. Выравнивание 10
5. Обзор биоинформатических инструментов 11
Глава 1. Сравнение последовательности ДНК с эталонным геномом. 12
1. Основные принципы работы программы FastqScreen 12
2. Описание работы алайнеров bowtie и bowtie2 и их роль в программе FastqScreen 15
3. Критерий принятия решения 17
3.1. Введение 17
3.2. Описание критерия 19
3.3. Реализация 22
4. Программная реализация FastqScreen 23
4.1. Исходная реализация на языке Perl 23
4.2. Реализация по средствам языка Java 24
4.3. Результаты замеров производительности 27
Глава 2. Применение метода нечеткого поиска к задаче поиска повторяющихся нуклеотидов последовательности ДНК 29
1. Описание форматов SAM/BAM/CRAM 29
2. Описание набора инструментов Picard Tools 32
3. Базовый алгоритм 34
4. Метод нечеткого поиска 36
5. Реализация метода нечеткого поиска 40
Заключение 43
Список литературы 44
1. Введение в предметную область
За последние несколько лет все большую популярность набирает такое направление в науке как биоинформатика. Это обусловлено тем, что в данной области огромную роль играет использование современных компьютерных технологий для анализа и манипуляции большими биологическими данными. В широком смысле биоинформатика — это использование программного обеспечения для решения различных биологических задач, таких как анализ экспериментальных данных, визуализация биологических визуализация, поиск вариаций в геноме. Целью этой науки является не только развитие алгоритмов использования биологических данных, но и подготовка данных для получения наиболее достоверных результатов. Таким образом, биоинформатика включает в себя как естественно-научные дисциплины, такие как генетика, молекулярная биология, так и точные науки — математика, статистика, информатика.
Биоинформатика берет свое начало из биологии, а точнее, молекулярной биологии, одной из задач которой является изучение генетической информации. Генетическая (наследственная) информация хранится в клетках живых существ. Единицами хранения такой информации являются молекулы ДНК, РНК и белки.
Молекула ДНК (дезоксирибонуклеиновая кислота) представляет собой двухцепочечную спиралевидную молекулу, состоящую из нуклеотидных звеньев, которые удерживаются водородными связями. Каждый нуклеотид последовательности ДНК содержит в себе азотиcтое основание одного из четырех видов: аденин (A), гуанин (G), тимин (T) и цитозин (C). В природе пары оснований соединены между собой согласно принципу комплементарности: аденин соединяется только с тимином, гуанин — только с цитозином [1].
Приведем некоторые термины из области биологии:
Ген — это некоторый участок ДНК, представляющий собой последовательность нуклеотидов. Количество генов у организмов различно. Геном — это общая генетическая информация, содержащаяся в клетках организма. Например, человеческий геном состоит из 23 пар хромосом и митохондриальной ДНК. Геномы большинства организмов состоят из ДНК, однако геном некоторых вирусов построен из РНК. Подробную информацию по этому вопросу можно почерпнуть, например, в [1].
2. Процесс секвенирования
Одной из первых задач, которая встает перед исследователями, является получение генетической информации некоторого организма, пригодное для компьютерного представления. Процесс формального описания первичной структуры линейной макромолекулы в виде последовательности мономеров в текстовом виде называется секвенированием. Иными словами, это процесс определения последовательности нуклеотидов, из которых состоят изучаемые данные, и приведение их к цифровому виду. Прибор, с помощью которого экспериментатор из некоторого организма получает текстовый файл для последующей компьютерной обработки, называется секвенатором. Результат работы секвенатора — это текстовый файл, содержащий огромное количество ридов. Рид — это последовательность нуклеотидов, представляющая собой небольшой фрагмент ДНК. Размер этого фрагмента варьируется настройками секвенатора.
В основе процесса секвенирования лежит метод полимеразной цепной реакции, или, сокращенно, ПЦР. Этот метод является мощным инструментом для подготовки данных к формату, в цифровом виде. ПЦР — это метод, при котором происходит амплификация ДНК, то есть увеличение числа копий фрагментов ДНК, с помощью термических и некоторых химических воздействий. Этот метод цикличен и каждая итерация удваивает концентрацию фрагментов, поэтому число копий растет со скоростью , где n — количество циклов. Это происходит для того, чтобы секвенатор смог определить позиции нуклеотидов в фрагменте ДНК.
Перечислим стадии ПЦР:
• Выделение ДНК из биологического материала (пробы).
• Фрагментация, то есть деление ДНК на короткие фрагменты, с помощью ультразвука или ферментов.
• После чего, над полученными фрагментами, производится термическое воздействие, с целью разрыва водородных связей и, как следствие, получения одноцепочечных фрагментов.
Рис. 0.2.1. Фрагмент ДНК с разорванными водородными связями.
• Далее, происходит подбор праймеров. Праймер — это очень короткая последовательность нуклеотидов, которая получается путем добавления комплементарных нуклеотидов на концы фрагмента.
Рис. 0.2.2. Добавление праймеров на концы фрагментов.
• Затем, к фрагменту ДНК добавляется ДНК-полимераза, которая доводит синтез до конца фрагмента.
Рис. 0.2.3. Добавление ДНК-полимеразы.
• Таким образом, мы получили два одинаковых фрагмента и можем повторить процесс.
Существует множество подходов к секвенированию биологических данных: полупроводниковое секвенирование, секвенирование синтезом, секвенирование при помощи нанопор, одномолекульное секвенирование в реальном времени и другие. Самым популярным на данный момент является метод секвенирования по Сэнгеру. Этот метод был изобретен Фредериком Сэнгером в 1977 и пользуется большой популярностью по сей день, поскольку он обладает достаточно высокой точностью. Более детально познакомиться с данным механизмом можно в [2].
Развитие современных технологий положило начало методам секвенирования нового поколения (next-generation sequencing, или сокращенно NGS). К наиболее популярным NGS секвенаторам относятся:
• Illumina. Основан на методе секвенирования синтезом, обрабатывающий короткие риды длиной 100bp1-300bp. Отличительной особенностью Illumina является мостиковая амплификация, которая позволяет производить чтения фрагментов ДНК с обоих концов, что существенно облегчает многие задачи, связанные со сборкой генома. Кроме того этот секвенатор обладает высокой точностью — 99.9% и хорошей производительностью. Однако, имеет достаточно низкую скорость секвенирования: полная сборка генома человека может длиться неделю, что критично для клинических исследований. В настоящее время является наиболее широко используемым.
• IonTorent. Полупроводниковый секвенатор, также работающий с короткими ридами, имеет высокую скорость секвенирования, но обладает более низкой точностью — 99%.
• К более современным технология секвенирования относятся секвенаторы Oxford nanopore и PacBio, основанные на секвенирование при помощи нанопор и одномолекульном секвенирование в реальном времени соответственно. Отличительно особенностью этих подходов является то, что они способны обрабатывать риды длиной до нескольких тысяч оснований. Тем не менее, эти инструменты имеют низкую точность (приблизительно 85%).
3. Форматы fasta и fastq
Результатом работы секвенатора является файл, доступный для чтения компьютеру и даже человеку. Существует два формата представления секвенированных данных: это форматы fasta и fastq. Файлы в данных форматах содержат информацию об изучаемых фрагментах ДНК в виде нуклеотидной последовательности. Отличие этих форматов заключается в том, что в файле fastq, в отличие от fasta, содержится еще и информация о качестве. Рассмотрим эти форматы подробнее.
Формат fasta — это формат представления данных секвенирования в виде последовательности биологической информации совместно с данными о качестве. Ниже представлен пример представления одного рида в формате fasta:
>SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAA
Этот формат достаточно прост и строится следующим образом: первая строка начинается с символа '>' и содержит уникальный идентификатор последовательности, после которого может идти некоторая дополнительная информация; вторая строка содержит саму последовательность нуклеотидов. Для последовательности ДНК буквы A, C, G и T кодируют, соответственно, аденин, цитозин, гуанин и тимин. Кроме того, в файле также может присутствовать буква N, что означает, что секвенатор по какой-либо причине не смог однозначно определить нуклеотид на данной позиции.
Формат fastq очень похож на fasta, но представляет рид с помощью четырех строк:
@SEQ_ID
GATTTGGGGTTCAAAGCAGTATCGATCAAATAGTAAATCCATTTGTTCAA
+
!''*((((***+))%%%++)(%%%%).1***-+*''))**55CCF>>>>>
Первые две строки содержат аналогичную информацию, что и fasta, с тем лишь отличием, что идентификатором начала первой строки является символа '@'. Третья строка начинается с символа '+' и может содержать дополнительную служебную информацию или же оставаться пустой, если такая информация не требуется. Четвертая строка содержит символы, кодирующие уровень качества для каждого из нуклеотидов. Её длина должна совпадать с длиной последовательности нуклеотидов. Рассмотрим подробнее как вычисляется качество каждого нуклеотида.
Уровень качества Q — это некоторое целое число, соответствующее вероятности того, что данный нуклеотид некорректен. Существует два подхода вычисления Q:
Шкала PHRED. В этой шкале каждому нуклеотиду соответствует значение, которое логарифмически зависит от вероятности ошибки:
= ,
где есть вероятность того, что соответствующий нуклеотид ошибочный.
Шкала Solexa. В этом случае уровень качества вычисляется следующим образом:
= ,
Для высокого уровня качества оба этих подхода дают одинаковый результат, но для низкого уровня качества различаются.
Для кодирования символов, характеризующих качество, используется таблица ASCII. ASCII-код символа равен уровню качества Q плюс некоторая константа. В случае шкалы PHRED эта константа равна 33, а для шкалы Salexa — 64. В обоих случаях код ASCII-код символа не должен превышать 127.
4. Выравнивание
Одним из основных инструментов биоинформатики является выравнивание. Выравнивание позволяет сопоставить полученный файл fastq некоторому эталонному референсному геному fasta. Программа, выполняющая выравнивание, называется алайнером (от англ. aligner) или выравнивателем. Существует около сотни алайнеров и многие из них имеют низкий уровень доверия или основаны на алгоритмах, решающих узкий круг задач. Наиболее популярными алгоритмами, решающими задачу выравнивания являются хэш-таблицы, суффиксные деревья и алгоритм слияния-сортировки. В рамках данной работы мы будем рассматривать только алайнеры, основанные на построении суффиксного дерева, а именно, bowtie и bowtie2 [4]. Также, к этому классу программ-выравнивателей относится весьма популярный алайнер BWA.
Bowtie, bowtie2 и BWA являются высокопроизводительными алайнерами, обладающими хорошей точностью и высоким индексом популярности, предназначенные для выравнивания коротких ридов на большой референсный геном. Для проведения процедуры выравнивания рассматриваемым алайнерам обязательно требуется построение FM-индекса. FM-индекс — это сжатый индексный файл, позволяющий производить навигацию по референсу и экономить оперативную память. В основе построения FM-индекса лежит преобразование Барроуза-Уилера [6].
В результате выравнивания мы получаем файлы в формате SAM, BAM или CRAM, с которыми подробнее познакомимся во второй главе настоящей работы. Эти форматы играют ключевую роль в биоинформатических исследованиях и являются наиболее популярными в данной области. Они открывают широкий спектр возможностей по манипуляции, анализу, визуализации биологических данных. Именно эти форматы используются в лабораторных и клинических исследованиях, в диагностике раковых опухолей, наследственных заболеваний и многом другом.
5. Обзор биоинформатических инструментов
Рассмотренные выше технологии являются только подготовительной частью большинства биоинформатических исследований. Наибольший интерес представляет дальнейший анализ данных. В данном разделе мы проведем краткий обзор различных инструментов, которые используются биоинформатиками, и возможностей, которые они предоставляют.
Обычно, анализ секвенированных данных, а именно, файла в формате fastq, начинается с проверки качества полученных данных. Причины плохого качества входных данных могут быть самые разные: ошибки работы прибора, грязные пробирки. Некоторые инструменты позволяют проверить это. Например программа FastQC позволяет провести быстрый обзор файла и сообщить пользователю где могут быть проблемы. Кроме того, этот инструмент может собрать некоторую статистическую информацию о ридах. Еще одним инструментом проверки качества первичных секвенированных данных является программа FastqScreen, которая позволяет получить за достаточно короткий срок информацию о содержащихся в файле примесях.
Дальнейшим этапом исследования обычно является выравнивание, после чего, анализируются полученные файлы в формате SAM/BAM/CRAM. Одними из самых популярных программ для обработки биоинформатических данных являются SamTools и Picard Tools. Эти инструменты позволяют решать огромное количество биоинформатических задач, например, рассчитать покрытие нуклеотидов на референсе, собрать некоторую статистическую информацию о данных, удалить дубликаты и многое другое.
В настоящей работе, мы сфокусируем внимание на следующих задачах:
1. Детальное изучение программы FastqScreen и её оптимизация. Кроме того, мы предложим новое расширение функционала инструмента с целью упрощения его использования.
2. Улучшение одной из утилит Picard Tools путем оптимизации существующего алгоритма обработки данных.
В настоящей работе были изучены основные положения биоинформатики и был проведен ознакомительный анализ некоторых инструментов этой области: FastqScreen и Picard Tools.
Кроме того, был изучен алгоритм работы программы FastqScreen, переписан с языка Perl на Java и оптимизирован. Также, был сформулирован и построен критерий принятия решения, основанный на диаграмме размаха, который согласно какому-то наперед заданному уровню значимости, определяет, может ли полученный файл fastq быть использован в дальнейшем анализе.
Помимо вышесказанного, в данной работе, был улучшен алгоритм работы одного из инструментов Picard Tools. Изначально, утилита EstimateLibraryComplexity работала на посимвольном сравнении строк, что не являлось эффективным решением данной задачи. Поэтому, этот алгоритм был изменен: теперь, алгоритм сравнения строк, основывается на построении хэш-таблицы.
1. David W. Mount, Bioinformatics: Sequence and Genome Analysis, Cold Spring Harbor Laboratory Press, 2004.
2. Д. В. Ребриков, Д. О. Коростин, Е. С. Шубина, В. В. Ильинский, NGS: высокопроизводительное секвенирование, Бином, 2004.
3. Michael Frigge, David C. Hoaglin and Boris Iglewicz, Opening the Box of a Boxplot, The American Statistician 42 (4): 257–262, 1988.
4. Ben Langmead, Cole Trapnell, Mihai Pop and Steven L Salzberg, Ultrafast and memory-efficient alignment of short DNA sequences to the human genome, Genome Biology, 2009.
5. Alexandr Andoni, Nearest Neighbor Search: the Old, the New, and the Impossible, Massachusetts Institute of Technology, 2009.
6. M. Burrows and D. Wheeler, A block sorting lossless data compression algorithm, Technical Report 124, Digital Equipment Corporation, 1994.