В подавляющем
большинстве случаев обязательным, неизбежным и неотвратимым шагом в анализе
данных NGS является картирование (или выравнивание) полученных
ридов относительно референсного генома (или иной, интересующей исследователя,
нуклеотидной последовательности). В свою очередь этому шагу предшествует
индексация референсного генома: преобразование большого генетического текста в
компактную хеш-таблицу (hash table), позволяющую быстро определять геномное
происхождение ридов из .fastq файла.
Хеш-таблица имеет
два столбца: ключ (key) и значение (value). Ключом в такой таблице является олигонуклеотид, извлеченный из
референсной последовательности, а значением будут уникальные геномные
координаты или список геномных координат (в случае повторяющихся олигонуклеотидов),
по которым можно найти олигонуклеотид в геноме. При картировании первоначально из
рида извлекаются фрагменты генетического текста, по хеш-таблице определяется их
вероятное геномное происхождение, и только после этого по рабочему алгоритму
программы выравнивания (элайнера) принимается решение о геномном происхождении
всего рида.
В настоящее время
преобладающими по разнообразию и популярности являются алгоритмы картирования,
основанные на парадигме seed-and-extend. Это не обозначает, что такой подход является
единственно возможным и единственно правильным. Он имеет как свои плюсы, так и
свои минусы и именно последние стимулируют появление новых подходов. Одним из
таких новых подходов является парадигма картирования seed-and-vote.
Эта
парадигма была разработана австралийскими исследователями Yang Liao, Gordon K. Smyth и Wei Shi, и впервые
описана в 2013 году в статье “Liao Y., Smyth G. K., Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. // Nucleic Acids Res. – 2013 May
1;41(10):e108. doi: 10.1093/nar/gkt214”. Суть
идеи австралийцев довольно проста и схематически проиллюстрирована на нижеследующем
рисунке. На первом этапе из рида извлекаются суб-риды – непрерывающиеся
фрагменты целевой последовательности длиной в 16 нуклеотидов.
Методология извлечения суб-ридов (со сдвигом в три нуклеотида или же со сдвигом в
один нуклеотид) зависит от того, какое значение было задано аргументу gappedIndex функции buildindex на этапе
индексации генома (описание аргументов этой функции дано ниже). |
Картирование ридов, основанное на парадигме seed-and-vote
(источник: Liao Y., Smyth G. K., Shi W. The Subread aligner: fast, accurate and
scalable read mapping
by
seed-and-vote. // Nucleic Acids Res. – 2013 May 1;41(10):e108. doi:
10.1093/nar/gkt214) |
На втором этапе
каждый суб-рид проверяется по хеш-таблице и получает координаты
его геномного происхождения. При этом возможны три неоднозначные ситуации: 1)
суб-рид отсутствует в хеш-таблице, что наиболее вероятно обусловлено его неинформативностью
(объяснение дано ниже в описании аргументов функции buildindex); 2) суб-рид может быть локализован в хеш-таблице только при условии
редактирования его последовательности (из-за наличия одного или нескольких
несовпадающих нуклеотидов, причем допускаемое количество несовпадений может
контролироваться пользователем); 3) суб-риду присваиваются не уникальные геномные координаты, а список из
геномных координат. Последняя ситуация обусловлена несколькими причинами.
Во-первых, суб-риды имеют длину всего в 16 нуклеотидов, что допускает множественность локализации в столь большом
геноме, как геном человека, в силу случайного совпадения. Во-вторых, в геноме
имеется множество коротких или длинных повторяющихся последовательностей, что
неизбежно порождает неоднозначность геномного происхождения для многих
суб-ридов. В-третьих, как уже отмечалось, на этапе картирования допускается
незначительное редактирование суб-ридов, что также может привести к
неоднозначной геномной локализации таких коротких последовательностей.
В конечном итоге к
началу третьего этапа картирования для каждого рида мы получаем сложную
картинку распределения его суб-ридов по геному: в геноме будет найдено
множество мест, куда попадает разное количество суб-риды. Что бы окончательно
определиться с происхождением рида на третьем этапе геномные интервалы, куда
попал хотя бы один суб-рид, сравниваются между собой по количеству суб-ридов (проходит
«голосование») и в качестве истинного места происхождения рида выбирается тот
генетический локус, который получает максимум «голосов». Не сложно догадаться,
что не такими уж и редкими будут ситуации, когда два или больше мест в геноме
получат одинаковый максимум «голосов». Это будут случаи множественного
картирования, управление которыми обеспечивается с помощью соответствующих
аргументов функции выравнивания.
Но прежде
чем провести картирование нам нужно получить хеш-таблицу. Для ее создания мы
можем воспользоваться функцией buildindex из
библиотеки функций Rsubread, разработанной сотрудниками Отдела биоинформатики Университета
Мельбурна (Мельбурн, Австралия) Wei Shi и Yang Liao. Об этой библиотеке уже
говорилось в разделе “Идентификация сайтов SNP по данным NGS: Использование R функции exactSNP”,
поэтому мы опустим описание путей поиска и установки этого программного
продукта на локальный компьютер, а сразу же перейдем к обсуждению особенностей работы
с данной функцией и ее и возможностей.
Аргументы
функции buildindex.
Функция buildindex
имеет семь аргументов.
q basename – базовое название (например, “Index_HG19”), которое будет использовано как префикс для
обозначения всех создаваемых файлов с индексами. Здесь же
целесообразно прописать путь и папку, куда будут записаны все эти файлы (если
она не задана как рабочая).
q reference – название файла (и путь к нему, если он не находится в рабочей папке),
содержащего сиквенс целевого генома в формате FASTA.
q gappedIndex – логический
оператор, позволяющий выбирать один из двух способов извлечения суб-ридов из
референсной последовательности. Если данному аргументу присвоено значение TRUE (по
умолчанию), то суб-риды будут извлечены из референсной последовательности со сдвигом
в три нуклеотида, если же FALSE –
то со сдвигом в один нуклеотид. Последний вариант требует больше времени на
расчеты и больше аппаратных ресурсов, но он же дает более точные результаты по
картированию.
q indexSplit – логический
оператор, задающий способ сегментации хеш-таблицы. Если
данному аргументу присвоено значение FALSE, то создается одна цельная хеш-таблица. В противном случае (по
умолчанию indexSplit = TRUE) хеш-таблица будет разбита на множество сегментов (или блоков индексов).
Количество сегментов зависит от значения аргумента memory, размера референсного генома, а также значения
аргумента gappedIndex. При картировании ридов каждый из таких сегментов
будет загружаться в оперативную память компьютера и использоваться как единое
целое.
q memory – размер оперативной
памяти (в мегабайтах), отводимой под работу по созданию хеш-таблицы. По
умолчанию отводится 8000 Мбайт. Заданный аргументом размер оперативной памяти
будет задействован и на этапе картирования ридов, при этом чем больше
оперативной памяти может быть использовано, тем быстрее будет проходить
картирование.
q TH_subread – порог
частоты встречаемости суб-рида в референсном геноме, выше которого он относится
к неинформативным (часто встречаемым) и исключается из хеш-таблицы, а ниже или
равно которому – к информативным и сохраняется в хеш-таблице. По умолчанию этот
аргумент имеет значение 24, что, по мнению авторов библиотеки Rsubread, является
оптимальным для большинства случаев индексации референсных геномов и
картирования ридов.
q сolorspace – логический оператор, позволяющий задавать формат кодировки
нуклеотидной последовательности. По умолчанию этому аргументу присвоено
значение FALSE. В этом случае нуклеотидная
последовательность имеет прямую цифровую кодировку, где каждому из четырех
вариантов нуклеотидов присваивается соответствующая цифра (так называемая «base space» кодировка). Такой способ
кодировки используется в подавляющем большинстве технологий NGS. Если же этому аргументу
присвоено значение TRUE,
то это позволяет работать с последовательностями, записанными в цветовом
формате (так называемая «color space» кодировка). Такой формат используется в технологии
секвенирования SOLiD (Sequencing by Oligonucleotide Ligation and Detection), разработанной
компанией Life Technologies и предполагает кодировку последовательности с помощью цифр, соответствующих
четырем цветам – 0 (синий), 1 (зеленый), 2 (желтый) и 3 (красный).
Каждый из этих цветов определяет динуклеотид, а не отдельный нуклеотид (например,
динуклеотид TA соответствует красному цвету и цифре 3). Более подробно о самой
технологии секвенироания SOLiD можно узнать на сайте компании Thermo Fisher Scientific Inc., которая занимается
производством и реализацией оборудования, реагентов и расходных материалов для
данной технологии. Кроме того, цветовой формат кодировки нуклеотидных
последовательностей обсуждается в ряде аналитических работ (например, «Overview of sequencing technology platforms” издательства Springer), посвященных NGS, а так же на форумах, связанных
с биоинформатикой (например, Biostars).
Консолидированный
R код.
Ниже представлен
обобщенный R
код для индексации референсного генома с помощью функции buildindex. В
качестве примера мы взяли референсную сборку GRCh37/GCA_000001405.1
генома человека, сиквенс которой можно получить в формате .2bit по адресу http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips
rm(list = ls())
setwd("/media/gvv/NGS_Illumina_MiSeq/")
library(Rsubread)
INDEX = buildindex(basename = "Reference_Genomes/Index_HG19/",
reference = "Reference_FASTA_files/hg19.fa",
gappedIndex = TRUE,
indexSplit = FALSE,
memory = 7000,
TH_subread = 24,
colorspace = FALSE)
На выполнение кода
уходит 30-40 минут – в зависимости от вычислительных возможностей компьютера.
После завершения расчетов в рабочую директорию будет записано 4 файла данных (при
условии, что аргумент indexSplit = FALSE) с префиксом, заданным аргументом basename: .array, .tab, .reads и .files, - а так же один лог-файл. Общий размер этих файлов составит около 5,7
Гбайт. Полученные индексы могут быть использованы для картирования ридов с
помощью элайнеров Subread и Subjunc, что будет обсуждаться в следующих сообщениях. При этом, единожды
построенные, они могут использоваться многократно, пока, например, не появится обновление
референсной последовательности.
Комментариев нет:
Отправить комментарий