четверг, 16 июня 2016 г.

Загрузка файла формата BAM в рабочее пространство R среды.

BAM файлы являются основными итоговыми файлами, генерируемыми элайнерами. Они же являются основным источником информации при проведении разнообразных вариантов постанализа данных высокопроизводительного секвенирования (идентификации сайтов простого нуклеотидного полиморфизма, анализа экспрессии генов, сборки транскриптов и т. д.). Поэтому нет ничего удивительного, что для работы с BAM файлами разработано разнообразное программное обеспечение, написанное как на разных языках, так и обладающее разными возможностями. Провести сравнительный анализ такого программного обеспечения дело трудоемкое и не совсем целесообразное, а вот рассмотреть базовые возможности R/Bioconductor по чтению BAM файлов и манипуляций над ними вполне посильно.
Ключевой R/Bioconductor библиотекой функций, предназначенных для работы с BAM файлами, является пакет Rsamtools. Пакет Rsamtools может быть установлен с официального сайта проекта Bioconductor с помощью стандартной процедуры:

source("https://bioconductor.org/biocLite.R")
biocLite("Rsamtools")

На этом же сайте можно найти стандартное руководство к пакету и две виньетки, а на видеохостинге YouTube можно посмотреть ряд видеоуроков, например, введение в Rsamtools от Kasper Hansen или использование функции pileup от Nathaniel Hayden.

Загрузка BAM файла в рабочее пространство R среды.
С практической точки зрения доступ к данным, хранящимся в BAM файле, лучше всего создавать через объект класса BamFile. Объект класса BamFile является не самим BAM файлом, загруженным в рабочее пространство R среды, а лишь ссылкой на этот файл и файл с индексами. Использование такого объекта вместо прямой загрузки самого BAM файла позволяет обрабатывать данные из BAM файла с помощью разных методов без использования больших объемов оперативной памяти и без повторных загрузок в R среду файла с индексами. Ниже представлен код для создания такого объекта:

rm(list = ls())
library(Rsamtools)
dataBAM = BamFile(file = "/media/gvv/NGS_Illumina_MiSeq/Files_BAM/fusGene.bam",
        index = "/media/gvv/NGS_Illumina_MiSeq/Files_BAM/fusGene.bam.bai")
print(dataBAM)
class: BamFile
path: /media/gvv/NGS_Illumina_MiSeq/Files_BAM/fusGene.bam
index: /media/gvv/NGS_Illumina_MiSeq/Files_BAM/fusGene.bam.bai
isOpen: FALSE
yieldSize: NA
obeyQname: FALSE
asMates: FALSE
qnamePrefixEnd: NA
qnameSuffixStart: NA

Чтение раздела “Заголовок”.
Содержание раздела “Заголовок” можно просмотреть с помощью функции scanBamHeader(). Эта функция имеет два аргумента:
q  files – объект класса BamFile (или же путь к BAM файлу и его название, если объект BamFile не создавался).
q what – элемент заголовка, который должен быть отображен в консоли. Может принимать два значения: "targets" (будет выведена информация о референс-последовательностях) и "text" (будет выведена вся остальная информация, содержащаяся в заголовке).

scanBamHeader(files = dataBAM, what = c("targets", "text"))
$targets
     chr8
103888433
$text
$text$`@HD`
[1] "VN:1.0"        "SO:coordinate"
$text$`@SQ`
[1] "SN:chr8"        "LN:103888433\r"
$text$`@PG`
[1] "ID:subread"           "PN:subread"           "VN:Rsubread 1.20.1\r"

Чтение раздела “Картирование”.
Как уже обсуждалось в сообщениях Парадигма картирования seed-and-vote: итоговые файлы в формате SAMи “Парадигма картирования seed-and-vote: итоговые файлы в формате BAM” каждая запись (строка) раздела “Картирование” SAM/BAM файла содержит 11 обязательных полей, расположение и название которых постоянны и, кроме того, может содержать ряд дополнительных необязательных полей, в том числе предложенных самим пользователем.
Основной функцией, позволяющей работать с данными, содержащимися в разделе “Картирование”, является функция scanBam(). У этой функции только три аргумента:
q file – объект класса BamFile (или же путь к BAM файлу и его название, если объект BamFile не создавался).
q index – путь к BAI файлу и его название.
q param – параметры сканирования BAM файла.
Аргумент param является ключевым у функции scanBam(), так как именно через этот аргумент можно задать, что же должно извлекаться из BAM файла. Значение аргумента param задается с помощью функции низкого уровня ScanBamParam(), которая, в свою очередь, имеет восемь аргументов.
q which – объект класса GRanges или класса RangesList (или же какой-либо другой объект, который может быть преобразован в объект указанных классов), содержащий координаты интересующего исследователя интервала референс-последовательности. Функция ScanBamParam() преобразует этот объект в объект класса IRangesList, длина которого равна количеству заданных интервалов. Если на функцию scanBam() через аргумент which подан такой объект, то из BAM файла будут извлекаться только те записи, которые перекрываются с указанным интервалом. Если же этот аргумент не задействован, то будут извлекаться все записи, вне зависимости от их координат.
q what – название (или символьный вектор с несколькими названиями) обязательного поля (или полей), которое необходимо извлечь из BAM файла. Если этот аргумент не задан, то ни одно из обязательных полей не будет извлечено. Если нужно извлечь все поля, то можно подать на аргумент специальный “помощник конструктора” scanBamWhat().
q flag – чтение поля FLAG. Значение этого аргумента может быть задано с помощью функции низкого уровня scanBamFlag(). Эта функция имеет одиннадцать аргументов: isPaired, isProperPair, isUnmappedQuery, hasUnmappedMate, isMinusStrand, isMateMinusStrand, isFirstMateRead, isSecondMateRead, isSecondaryAlignment, isNotPassingQualityControls и isDuplicate. Каждый из этих аргументов может принимать одно из трех логических значений: TRUE, FALSE или NA (при чтении BAM файла соответствующий FLAG игнорируется и загружаются все строки вне зависимости от его значения).
q simpleCigar – чтение поля CIGAR. Может принимать одно из двух логических значений: если TRUE, то будут загружены только те строки, для которых поле CIGAR имеет значение либо * (значение поля не определено), либо M (все нуклеотиды рида полностью совпадают с референс-последовательностью), а если FALSE, то будут загружены все строки.
q reverseComplement – логический аргумент, позволяющий экстрагировать из BAM файла обратные комплементарные последовательности тех ридов, которые картированы по минус цепи. По умолчанию имеет значение FALSE.
q tag – название ярлыка (или символьный вектор с названиями нескольких ярлыков) необязательного поля (или полей). При наличии у этого аргумента актуального значения из BAM файла будут извлечены только те записи, которые имеют соответствующее поле, причем вместе со значениями этого поля.
q tagFilter – именованный объект класса list, содержащий список ярлыков необязательных полей и их значений, которые необходимо извлечь из BAM файла. При наличии у этого аргумента актуального значения из BAM файла будут извлечены только те записи, которые имеют не просто соответствующее поле (или поля), а еще и с указанными значениями.
q mapqFilter – целочисленное положительное минимальное значение качества картирования ридов. Если поле MAPQ (а точнее поле mapq, так как речь идет о BAM файле) имеет значение, меньшее чем подано на аргумент mapqFilter, то такая запись не будет извлекаться из BAM файла.
Теперь мы можем выразить все сказанное выше в виде одного консолидированного R кода:

rm(list = ls())
library(Rsamtools)
dataBAM = BamFile(file = "/media/gvv/NGS_Illumina_MiSeq/Files_BAM/fusGene.bam",
        index = "/media/gvv/NGS_Illumina_MiSeq/Files_BAM/fusGene.bam.bai")
which = GRanges(seqnames = Rle("der8"), ranges = IRanges(start = 91960000, end = 91960100))
what = c("pos", "isize")
flag = scanBamFlag(isPaired = TRUE, isUnmappedQuery = FALSE, isMinusStrand = TRUE)
param = ScanBamParam(which = which,
                                          what = what,
                                          flag = flag,
                                          simpleCigar = FALSE,
                                          reverseComplement  = FALSE,
                                          tag = c("NH", "NM"),
                                          tagFilter = list(NM = c(0, 1, 2)),
                                          mapqFilter = 50)
res = scanBam(file = dataBAM, param = param)

Первоначально мы создали объект dataBAM класса BamFile. В качестве примера был взят заранее подготовленный небольшой BAM файл, содержащий RNA-Seq риды гибридного онкогена RUNX1-RUNX1T1, картированные по деривату 8-й хромосомы с помощью элайнера Subjunc.
Далее с помощью стандартной функции GRanges() был задан интересующий нас интервал which референс-последовательности, прописано, какие именно поля what и значения flag нам нужны.
После этого были прописаны все остальные параметры сканирования. В частности, через параметры сканирования был создан запрос на два необязательных поля NH (количество записей, где фигурирует текущий рид; number of hits) и NM (количество нуклеотидов рида, не соответствующих референс-последовательности (за исключением “свисающих” нуклеотидов); number of mismatches), причем через аргумент tagFilter было введено ограничение извлекать из BAM файла только те записи, которые по полю NM имеют значения 0, 1 или 2 (т. е. содержат не более 2 мисматчей). Кроме того, через аргумент mapqFilter было введено ограничение не извлекать записи, где качество картирования ниже 50.
Только после этого мы подали на функцию scanBam() объект dataBAM и параметры сканирования. Следует обратить внимание на то, что при подаче на эту функцию не пути к BAM файлу и его названия, а объекта класса BamFile, прописывать аргумент index не нужно, так как информация о BAI файле уже есть в объекте класса BamFile. В противном случае будет выдана ошибка:
Warning message:
In scanBam(file = dataBAM, index = dataBAM, param = param):
  'index' ignored for scanBam,BamFile-method
Полученный в результате работы функции scanBam() объект res является объектом класса list. Длина этого объекта равна количеству заданных интервалов. Содержание этого объекта показано ниже (для компактности показаны только восемь записей).

res
$`chr8:91960000-91960100`
$`chr8:91960000-91960100`$pos
  [1] 91959902 91959903 91959914 91959917 91959917 91959918 91959918 91959919
$`chr8:91960000-91960100`$isize
  [1] -809 -286 -179 -339 -188 -280 -142 -184
$`chr8:91960000-91960100`$tag
$`chr8:91960000-91960100`$tag$NH
  [1] 1 1 1 1 1 1 1 1
$`chr8:91960000-91960100`$tag$NM
  [1] 0 0 0 0 1 0 0 0
class(res)
[1] "list"

Этот объект содержит всю запрошенную выше информацию. Очевидно, что вариантов запроса может быть множество – все определяется той задачей (или задачами), которая стоит перед исследователем. Следовательно, и содержание объекта res будет различаться. Для дальнейшей работы с ним можно использовать те же методы, что и при работе с любым объектом класса list.

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

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