Обнаружил жесткую нехватку информации по биоинформатике в русском сегменте. Не знаю, уж, востребована она или нет, но хочу предоставить на суд читателя вводную часть, которую можно назвать практическая биоинформатика, которой мне очень не хватало для ознакомления с предметом. В этой главе я хочу описать путь, который пришлось пройти мне до настоящего момента, когда я уже не шарахаюсь от фраз: вот вам FASTQ файл и постройте мне bed graph для genome browser. Чтобы в дальнейшем вести разговор об интересном, хочу по диагонали пройтись по определениям и программам первичной обработки данных, без которых трудно говорить на одном языке.
Для начала немного определений. Будем считать, что хромосомы это одномерные координатные оси, от единицы до, примерно, 10е8. Длина оси зависит от длины хромосомы. Каждая точка оси, это целое число. Биологи и химики провели большое количество экспериментов и благодаря им смогли уже с большой точностью (около 90%) описать части хромосом. Эти описания называются аннотацией. Аннотация содержит информацию о длине хромосом, о координатах отдельных участков хромосом, наиболее известны из них ген, intron, exon. Этих участков огромное количество, но основное их свойство заключается в том, что это отрезки или набор отрезков расположены на координатной оси. Одни отрезки могут включать в себя другие или как-либо пересекаться. Вот набор сайтов на котором можно посмотреть аннотации геномов человека, мышки и тому подобных.
Биологи с химиками проводят свои эксперименты и в результате их операций над клетками получают раствор, в котором содержатся относительно небольшие куски ДНК или РНК (не очень хочется вдаваться в подробности разницы или одинаковости, просто последовательность нуклеотидов). Этот раствор пропускают через аппаратуру для секвенирования, выход которой — небольшие строки. Эти строки являются концами кусков ДНК или РНК, находившимися в растворе. Длина полученных из аппаратуры строк составляет всего 36-50 bases (длина строки в нуклеотидах) иногда больше, но на текущий момент вроде не более 200. Эти отрезки, полученные из секвенирующей аппаратуры и определенные, последовательностью нуклеотидов, называются ридами (от английского reads — “считывание”). Стоит отметить, что риды характеризуются только последовательностью нуклеотидов, а не расположением на геноме. Иногда эти последовательности дополнены строкой вероятностей, ставящей в соответствие позиции нуклеотида его вероятность нахождения на этой позиции. FASTA файл — это файл без вероятностей, FASTQ — это файл с вероятностями.
Далее, в зависимости от того, что стало результатом эксперимента — куски ДНК или РНК, проводят одну из двух методик секвенирования ChIP-seq или RNA-seq, соответственно. Более подробно о них рассказано тут http://en.wikipedia.org/wiki/DNA_sequencing.
После того, как дорогие машины по секвенированию отработают и выдадут результат в FAST A/Q файл, необходимо найти получившиеся последовательности в геноме. Для ChIP-seq очень быстрый поиск, возможный в домашних условиях, делает программа bowtie, находила миллионы ридов всего за 5 минут. Т.е. она ищет вхождение строки длиною 36-50 символов, состоящей из, как минимум, четырех буквенного алфавита в строках общей длиною до 10е9. Почему использован оборот как минимум: помимо стандартного использования A/T G/C алфавита часто добавляют символ N, заменяющий любую возможную букву (более детально http://genome.ucsc.edu/FAQ/FAQdownloads.html#download5). У программы множество параметров. Так, например, можно разрешить в строке(риде) ошибку или две (максимум до 3х ошибок). Она может искать не одно вхождение рида в геном, а много. Может хитро сортировать данные, например, если рид с одной ошибкой нашелся однозначно в геноме, а с двумя или тремя ошибками нашелся во многих частях генома, то вывести только первый результат. Этот процесс нахождения ридов в геноме называется мэппингом от английского mapping (еще его называют alignment). Алгоритм такого быстрого поиска очень интересен, но этому можно посвятить отдельную статью или найти ссылку на английскую статью, с сайта разработчика, где рассказывается как алгоритм bzip2 натолкнул их на такое решение. Программ по мэппингу и сайтов, где его делают онлайн, много; по ключевым словам blast, eland +genome в гугле можно найти дополнительную информацию.
Для RNA-seq процедура немного сложнее, для него выполняют ChIP-seq мэппинг, а далее обрабатывают те риды, которые не нашли себя при ChIP-seq. Хорошая программа, которая активно использует bowtie для своей работы, называется tophat. В результате сплайсинга и порожденного им изоформизма части рида могут находиться в разных местах генома. Например, первые 15 символов могут приходиться на один участок генома, а другие 11 — на другой. Это разделение рида на части, приходящееся на конец одного exon и на начало другого, называют splice-junctions. Tophat позволяет их находить, а также определяет новые возможные изоформы существующих генов.
Результатом работы этих программ является sam/bam файл, в котором содержится информация из FAST A/Q файла плюс информация о координатах на соответствующих спиралях хромосомы. Процесс от библиотеки до sam/bam файла, часто называют pipeline процедурой и во многих лабораториях стоит на потоке, поэтому стоит поинтересоваться, что за параметры и версии программ установлены по умолчанию. В общем, это то, на чем заканчивается вступительная часть, далее наступает очередь исследований. Можно отметить, что с этого момента у нас есть данные, находящиеся с некоторой достоверностью в одних и тех же координатах: координатная ось, аннотация, риды с соотвествующими координатами.
Можно приступать к анализу данных от простого выяснения равномерности и непрерывности, сложности, кончая сложными статистическими выкладками, которые помогают разделить риды на определённые группы: шум, нулевой уровень и обогащение. Анализ данных необходим для того, чтобы в дальнейшем можно было обоснованно отбросить те или иные данные.
Если вас заинтересовала вводная часть, то на каждом этапе я могу остановиться подробнее. К сожалению, простая вводная часть уже занимает несколько страниц сухого текста без картинок, поэтому описывать написанные программы и механизмы в этой главе не решился. Меня самого больше всего привлекает последний абзац, где я мельком упомянул про статистику. Хочется освятить существующие библиотеки и механизмы для работы с такими данными. Сюда можно приобщить методы, входящие в datamining (разные типы кластеринга), которые описаны на хабре. Как применить пуассоновское распределение для анализа данных без контроля, как применить сложную цепочку пуассона, ф-теста, для нахождения участков обогащенных ридами на геноме (Diarac delta function)? Полезны ли готовые библиотеки для работы с интервалами boost.intervals, boost.icl?
И, естественно, если эта тема интересна, может кто и подскажет, как и куда копать, в каких вопросах, и дополнит. А может и напишет свое. Ибо решение биологических вопросов без математики и программирования на данном этапе уже точно невозможно. Есть ресурсы англоязычные, где обсуждаются похожие вопросы www.seqanswers.com. Но мне хотелось бы в дальнейшем уйти от описания готовых продуктов и заняться обсуждением правомерности применения математических и статистических методов в этих программах и возможности применения новых методов.
На текущем этапе своей работы мы пытались найти параметры, по которым можно фильтровать риды, как на интересные для исследования, так и на шумовые. Задача довольно нетривиальная при отсутствии контроля. В дальнейшем решено в библиотеку для секвенирования добавлять контроль, который позволит измерять уровень ошибки, но статистика все равно не отпустит.
Для начала немного определений. Будем считать, что хромосомы это одномерные координатные оси, от единицы до, примерно, 10е8. Длина оси зависит от длины хромосомы. Каждая точка оси, это целое число. Биологи и химики провели большое количество экспериментов и благодаря им смогли уже с большой точностью (около 90%) описать части хромосом. Эти описания называются аннотацией. Аннотация содержит информацию о длине хромосом, о координатах отдельных участков хромосом, наиболее известны из них ген, intron, exon. Этих участков огромное количество, но основное их свойство заключается в том, что это отрезки или набор отрезков расположены на координатной оси. Одни отрезки могут включать в себя другие или как-либо пересекаться. Вот набор сайтов на котором можно посмотреть аннотации геномов человека, мышки и тому подобных.
Биологи с химиками проводят свои эксперименты и в результате их операций над клетками получают раствор, в котором содержатся относительно небольшие куски ДНК или РНК (не очень хочется вдаваться в подробности разницы или одинаковости, просто последовательность нуклеотидов). Этот раствор пропускают через аппаратуру для секвенирования, выход которой — небольшие строки. Эти строки являются концами кусков ДНК или РНК, находившимися в растворе. Длина полученных из аппаратуры строк составляет всего 36-50 bases (длина строки в нуклеотидах) иногда больше, но на текущий момент вроде не более 200. Эти отрезки, полученные из секвенирующей аппаратуры и определенные, последовательностью нуклеотидов, называются ридами (от английского reads — “считывание”). Стоит отметить, что риды характеризуются только последовательностью нуклеотидов, а не расположением на геноме. Иногда эти последовательности дополнены строкой вероятностей, ставящей в соответствие позиции нуклеотида его вероятность нахождения на этой позиции. FASTA файл — это файл без вероятностей, FASTQ — это файл с вероятностями.
Далее, в зависимости от того, что стало результатом эксперимента — куски ДНК или РНК, проводят одну из двух методик секвенирования ChIP-seq или RNA-seq, соответственно. Более подробно о них рассказано тут http://en.wikipedia.org/wiki/DNA_sequencing.
После того, как дорогие машины по секвенированию отработают и выдадут результат в FAST A/Q файл, необходимо найти получившиеся последовательности в геноме. Для ChIP-seq очень быстрый поиск, возможный в домашних условиях, делает программа bowtie, находила миллионы ридов всего за 5 минут. Т.е. она ищет вхождение строки длиною 36-50 символов, состоящей из, как минимум, четырех буквенного алфавита в строках общей длиною до 10е9. Почему использован оборот как минимум: помимо стандартного использования A/T G/C алфавита часто добавляют символ N, заменяющий любую возможную букву (более детально http://genome.ucsc.edu/FAQ/FAQdownloads.html#download5). У программы множество параметров. Так, например, можно разрешить в строке(риде) ошибку или две (максимум до 3х ошибок). Она может искать не одно вхождение рида в геном, а много. Может хитро сортировать данные, например, если рид с одной ошибкой нашелся однозначно в геноме, а с двумя или тремя ошибками нашелся во многих частях генома, то вывести только первый результат. Этот процесс нахождения ридов в геноме называется мэппингом от английского mapping (еще его называют alignment). Алгоритм такого быстрого поиска очень интересен, но этому можно посвятить отдельную статью или найти ссылку на английскую статью, с сайта разработчика, где рассказывается как алгоритм bzip2 натолкнул их на такое решение. Программ по мэппингу и сайтов, где его делают онлайн, много; по ключевым словам blast, eland +genome в гугле можно найти дополнительную информацию.
Для RNA-seq процедура немного сложнее, для него выполняют ChIP-seq мэппинг, а далее обрабатывают те риды, которые не нашли себя при ChIP-seq. Хорошая программа, которая активно использует bowtie для своей работы, называется tophat. В результате сплайсинга и порожденного им изоформизма части рида могут находиться в разных местах генома. Например, первые 15 символов могут приходиться на один участок генома, а другие 11 — на другой. Это разделение рида на части, приходящееся на конец одного exon и на начало другого, называют splice-junctions. Tophat позволяет их находить, а также определяет новые возможные изоформы существующих генов.
Результатом работы этих программ является sam/bam файл, в котором содержится информация из FAST A/Q файла плюс информация о координатах на соответствующих спиралях хромосомы. Процесс от библиотеки до sam/bam файла, часто называют pipeline процедурой и во многих лабораториях стоит на потоке, поэтому стоит поинтересоваться, что за параметры и версии программ установлены по умолчанию. В общем, это то, на чем заканчивается вступительная часть, далее наступает очередь исследований. Можно отметить, что с этого момента у нас есть данные, находящиеся с некоторой достоверностью в одних и тех же координатах: координатная ось, аннотация, риды с соотвествующими координатами.
Можно приступать к анализу данных от простого выяснения равномерности и непрерывности, сложности, кончая сложными статистическими выкладками, которые помогают разделить риды на определённые группы: шум, нулевой уровень и обогащение. Анализ данных необходим для того, чтобы в дальнейшем можно было обоснованно отбросить те или иные данные.
Если вас заинтересовала вводная часть, то на каждом этапе я могу остановиться подробнее. К сожалению, простая вводная часть уже занимает несколько страниц сухого текста без картинок, поэтому описывать написанные программы и механизмы в этой главе не решился. Меня самого больше всего привлекает последний абзац, где я мельком упомянул про статистику. Хочется освятить существующие библиотеки и механизмы для работы с такими данными. Сюда можно приобщить методы, входящие в datamining (разные типы кластеринга), которые описаны на хабре. Как применить пуассоновское распределение для анализа данных без контроля, как применить сложную цепочку пуассона, ф-теста, для нахождения участков обогащенных ридами на геноме (Diarac delta function)? Полезны ли готовые библиотеки для работы с интервалами boost.intervals, boost.icl?
И, естественно, если эта тема интересна, может кто и подскажет, как и куда копать, в каких вопросах, и дополнит. А может и напишет свое. Ибо решение биологических вопросов без математики и программирования на данном этапе уже точно невозможно. Есть ресурсы англоязычные, где обсуждаются похожие вопросы www.seqanswers.com. Но мне хотелось бы в дальнейшем уйти от описания готовых продуктов и заняться обсуждением правомерности применения математических и статистических методов в этих программах и возможности применения новых методов.
На текущем этапе своей работы мы пытались найти параметры, по которым можно фильтровать риды, как на интересные для исследования, так и на шумовые. Задача довольно нетривиальная при отсутствии контроля. В дальнейшем решено в библиотеку для секвенирования добавлять контроль, который позволит измерять уровень ошибки, но статистика все равно не отпустит.
