Библиотека Rsubread содержит
два готовых к использованию элайнера: Subread и Subjunc. Первый из них проводит только локальное
картирование и подходит как для работы с DNA-seq, так и RNA-seq ридами.
Однако в случае с RNA-seq ридами этот элайнер не может картировать те риды,
которые попадают на области стыков двух экзонов. Для таких случаев нужно
глобальное картирование, при котором допускается, что разные части некоторых
ридов на уровне генома могут иметь разное происхождение. Такое картирование способен
провести элайнер Subjunc, стратегия работы которого схематически показана
на рисунке.
 |
Схема глобального картирования RNA-seq ридов, основанного на
парадигме seed-and-vote
(источник: Liao Y., Smyth G. K., Shi W. The Subread aligner: fast, accurate andscalable read mapping
|
Элайнер Subjunc сканирует каждый рид дважды. Во время первого сканирования
из рида извлекаются суб-риды и картируются относительно референсного генома по хеш-таблице. Для
каждого рида выбирается два геномных региона, получивших максимальное
количество “голосов” (детали картирования ридов с помощью подхода seed-and-vote описаны в разделе “Парадигма картирования seed-and-vote: индексация генома”). Далее принимается решение, включать ли эти регионы
в предварительную таблицу экзон-экзонных переходов. Положительное решение принимается при выполнении двух
условий: 1) между этими двумя регионами присутствуют донорный (“GT”) и акцепторный (“AG”) сайты сплайсинга; 2) суммарная длина геномных регионов L1 и L2, перекрываемых суб-ридами, участвовавшими в “голосовании”, равна длине L фрагмента рида, из которого произошли такие суб-риды. Если решение
положительное, то в предварительную таблицу экзон-экзонных переходов записываются
координаты предполагаемых сайтов сплайсинга. После этого создается еще одна таблица – таблица “якорных” мест для
каждого рида. Для этого из двух геномных регионов, по которым картирован рид,
выбирается самый лучший – тот, который получил наибольшее количество “голосов”,
- и записывается его крайняя левая координата. В эту же таблицу записывается
крайняя левая координата той части рида, которая картирована по самому “успешному”
геномному региону, а также размер самого региона.
Во время второго
сканирования проводится перебор всех возможных вариантов картирования рида в привязке
к “якорным” местах, в том числе рассматривается возможность картирования рида
целиком внутри экзона. При этом если рассматриваются варианты экзон-экзонных
переходов, то к несоответствиям между нуклеотидной последовательностью рида и
геномного региона предъявляются более строгие требования, чем при картировании
рида целиком внутри экзона. В конечном итоге, после перебора всех вариантов
картирования и выбора самых оптимальных для каждого из ридов, из первоначальной
таблицы предполагаемых сайтов сплайсинга исключаются те варианты, которые не
получают подтверждения.
Аргументы
функции subjunc.
Все описанные процедуры по глобальному картированию RNA-seq ридов могут быть проделаны с
помощью функции subjunc, которая является составной частью библиотеки Rsubread. Эта функция имеет несколько аргументов, которые
описаны ниже.
q index – базовое название файлов с индексами (с указанием пути и
папки, где записаны все эти файлы).
q readfile1 – название файла (и путь к нему, если он не находится в рабочей папке),
содержащего риды, предназначенные для картирования. Этот файл может быть в
формате FASTA или FASTQ. В том случае, если картировать нужно спаренные риды (paired-end reads),
аргументу присваивается название файла, содержащего первые риды.
q readfile2 – название файла, содержащего вторые риды в случае спаренных
ридов. По умолчанию этому аргументу присвоено значение NULL.
q input_format – формат файла, содержащего
риды. По умолчанию принимается формат gzFASTQ. Кроме того, функция может работать с форматами FASTQ, FASTA (или
его архивами), SAM и BAM.
q output_format – формат файла, в
котором будут записаны картированные риды. Аргумент может принимать два
значения – BAM (по
умолчанию) и SAM.
q output_file – название файла, в
котором будут записаны картированные риды.
q nsubreads –
количество суб-ридов, извлекаемых из каждого рида. По умолчанию этот аргумент
имеет значение 14, однако для сокращения времени, необходимого для картирования
большого количества ридов, авторы библиотеки Rsubread предлагают уменьшить этот аргумент до 10.
q TH1 – минимальное
количество картированных суб-ридов, необходимое для признания успешного
картирования рида (первого рида, если картируются спаренные риды). Хороший
результат (в плане точности и чувствительности) достигается тогда, когда
количество таких суб-ридов не меньше 30% от всего количества суб-ридов,
извлеченных из рида.
q TH2 –
минимальное количество картированных суб-ридов, необходимое для успешного
картирования второго рида (если картируются спаренные риды).
q maxMismatches – максимально
допустимое количество мис-матчей (несовпадений нуклеотидов между ридом и
референсной последовательностью). По умолчанию этот аргумент имеет значение 3.
q nthreads – количество
потоков, задействованных при проведении параллельных вычислений. По умолчанию
этому аргументу присваивается значение 1.
q indels – максимальноe количество инсерций и/или
делеций, допустимое во время картирования данного рида. По умолчанию этот
аргумент имеет значение 5.
q complexIndels – логический аргумент,
определяющий, будут ли детектироваться сложные инсерции/делеции. Под сложными
инсерциями/делециями подразумевается множество коротких близко расположенных
друг относительно (на расстоянии до 1 нуклеотида) инсерций/делеций,
локализованных в небольшом геномном регионе. По умолчанию аргументу присвоено значение FALSE.
q phredOffset – число, добавляемое к балам Фреда (Phred scores) при переводе буквенной кодировки качества прочтения нуклеотидов (ASCII letters-based base-calling quality scores) в цифровую. Этот аргумент может принимать
только два значения – 33 и 64. По умолчанию аргументу присвоено значение 33,
что применимо для файлов FASTQ, полученных на большинстве геномных секвенаторов
(например, Illumina HiSeq 2500).
q unique – логический
оператор, определяющий, будут ли сообщаться только уникально картированные
риды. По умолчанию этот аргумент имеет значение TRUE.
q nBestLocations – максимальное
количество равнозначных мест картирования, сообщаемое для рида. Допускаются
значения от 1 (по умолчанию) до 16 (включительно). Следует отметить, что
аргумент unique
имеет приоритет над аргументом nBestLocations.
q minFragLength – минимальный размер
сиквенированного фрагмента, на котором синтезирован рид. По умолчанию этому
аргументу присвоено значение 50.
q maxFragLength – максимальный
размер сиквенированного фрагмента, на котором синтезирован рид. По умолчанию этот
аргумент имеет значение 600.
q PE_orientation – ориентация ридов в паре. Аргумент может иметь одно
из трех значений: fr, ff и rf (где f обозначает смысловую, а r антисмысловую
цепь). По умолчанию ему присвоено значение fr – один рид из пары принадлежит
смысловой, а второй антисмысловой цепи сиквенированного фрагмента.
q nTrim5 – количество
нуклеотидов, удаляемых с 5’-конца каждого из ридов. По умолчанию аргумент имеет
значение 0.
q nTrim3 –
количество нуклеотидов, удаляемых с 3’-конца каждого из ридов. По умолчанию
аргумент имеет значение 0.
q readGroupID – идентификатор,
добавляемый к названию группы ридов и каждому риду в этой группе в конечном
файле. По умолчанию этому аргументу присвоено значение NULL.
q readGroup – название
в формате tag:value, добавляемое к группе ридов в конечном файле. По умолчанию
этот аргумент имеет значение NULL.
q color2base – логический
оператор, позволяющий провести конвертацию формата записи нуклеотидной
последовательности ридов из цветового в цифровой и сохранить риды в таком
формате в конечном файле. По умолчанию этому аргументу присвоено значение FALSE.
q reportAllJunctions –
логический оператор, позволяющий идентифицировать по данным RNA-seq любые сплайсинговые события без привязки к каноническим
сайтам сплайсинга. По умолчанию этот аргумент имеет значение FALSE.
Консолидированный
R код.
Ниже
представлен обобщенный R код по глобальному картированию RNA-seq ридов с
помощью функции subjunc. В качестве примера мы взяли пару FASTQ файлов SRR1145838,
содержащих риды после спаренного RNA-seq. Эти файлы выложены в свободный доступ в
репозитории European Nucleotide Archive.
rm(list = ls())
setwd("/media/gvv/NGS_Illumina_MiSeq/")
library(Rsubread)
GlobalAlign = subjunc(index = "Reference_Genomes/Index_HG19/",
readfile1
= "Files_FASTQ/SRR1145838_1.fastq",
readfile2
= "Files_FASTQ/SRR1145838_2.fastq",
input_format = "FASTQ",
output_format = "BAM",
output_file = "Files_BAM/SRR1145838.bam",
nsubreads = 14,
TH1 = 3,
TH2 = 3,
maxMismatches = 3,
nthreads = 8,
indels = 5,
phredOffset = 33,
unique = TRUE,
nBestLocations = 1,
minFragLength = 50,
maxFragLength = 600,
PE_orientation = "fr",
nTrim5 = 0,
nTrim3 = 0,
readGroupID = NULL,
readGroup = NULL,
color2base = FALSE,
reportAllJunctions = FALSE)
По мере выполнения кода в консоли будет идти
отображение хода расчетов, что в сокращенном виде показано на нижеследующем
рисунке. Время, затрачиваемое на выполнение кода, во многом зависит от
параметров, которые были использованы при индексации референсного генома,
вычислительных возможностей компьютера, количества ридов в исходном FASTQ файле
(или файлах, если речь идет о спаренных ридах), количества суб-ридов,
извлекаемых из каждого рида, а также значения аргумента reportAllJunctions (если
этому аргументу присвоено значение TRUE, то это
существенно увеличивает объем вычислений). После выполнения кода в рабочую
директорию будет записано три файла: 1) файл в формате BED, содержащий
информацию о сплайсинговых событиях (экзон-экзонных переходах); 2) файл в
формате INDEL, содержащий список идентифицированных инсерций и
делеций; 3) файл в формате BAM, содержащий риды и результаты их картирования
относительно референсного генома. Кроме того, если аргумент reportAllJunctions = TRUE, то
генерируется и сохраняется один дополнительный файл в формате TXT,
содержащий информацию об идентифицированных на уровне РНК межгенных сшивках. Более
подробно структура и содержание всех этих файлов мы рассмотрим в следующих сообщениях.
 |
| Окно вывода хода расчетов при работе функции subjunc |
Комментариев нет:
Отправить комментарий