суббота, 6 февраля 2016 г.

Парадигма картирования seed-and-vote: индексация генома.

В подавляющем большинстве случаев обязательным, неизбежным и неотвратимым шагом в анализе данных 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, что будет обсуждаться в следующих сообщениях. При этом, единожды построенные, они могут использоваться многократно, пока, например, не появится обновление референсной последовательности.

Комментариев нет:

Отправить комментарий