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.
Комментариев нет:
Отправить комментарий