July 28th, 2013

GenomicRanges

Дисклеймер: кратко основные вещи по Bioconductor GenomicRanges

library("GenomicRanges")
 
### Сборка объекта GenomicRanges с нуля
GenRange <- GRanges(
  seqnames = Rle(c("chrX", "chrY", "chrZ"), c(1, 1, 2)), # задание "хромосом" и числа блоков на каждой как Rle
  ranges = IRanges(start = c(1, 3, 19, 30), end = c(7, 5, 35, 40), names = c('One', 'Two', 'XX', 'ZZ')), # задание координат и имён блоков как IRanges
  strand = Rle(strand(c("+", "-", "*")), c(1, 1, 2)), # задание ориентации цепей блоков с указанием "числа" как Rle
 
  score = c(10, 100, 500, 400), # Опциональная колонка численных значений
  lolID = c('d', 'g', 't', 'y') # Произвольная колонка
)
GenRange
 
GenRange2 <- GRanges(
  seqnames = Rle(c("chrX", "chrY", "chrZ"), c(1, 1, 2)), # задание "хромосом" и числа блоков на каждой как Rle
  ranges = IRanges(start = c(1, 3, 19, 25), end = c(7, 5, 35, 30), names = c('One', 'Two', 'XX', 'DD')), # задание координат и имён блоков как IRanges
  strand = Rle(strand(c("+", "-", "*")), c(1, 1, 2)), # задание ориентации цепей блоков с указанием "числа" как Rle
 
  score = c(10, 100, 500, 400), # Опциональная колонка численных значений
  lolID = c('d', 'g', 't', 's') # Произвольная колонка
)
GenRange2
 
### Базовый вывод (и ввод) из объекта
seqnames(GenRange) # Rle-факторный вывод числа блоков и хромосом
seqlengths(GenRange) <- c(1e6, 1.2e6, 9e5) # Вывод и ВВОД длин хромосом
seqlengths(GenRange2) <- c(1e6, 1.2e6, 9e5)
strand(GenRange) # Rle-факторный вывод цепей
length(GenRange) # вывод числа блоков
names(GenRange) # вывод имён блоков
ranges(GenRange) # вывод диапазонов в виде IRanges
mcols(GenRange) # вывод колонок метаданных в IRangesовский DataFrame
start(GenRange) # получение вектора всех стартов
end(GenRange) # ...и всех концов
width(GenRange) # всех длин
GenRange[4] # вызов , скажем, 4 строки (блока)
 
  # манипуляции с компонентами
rep(GenRange[2], 3) # повторение блоков
rev(GenRange) # полное инвертирование порядка
window(GenRange, start = 3, end = 4) # поднабор по номерам строк (блоков)
seqselect(GenRange, start = c(1,2), end = c(2,3)) # почти то-же самое с поддержкой конкатенации, только нихрена не удобно.
 
  # преобразования
    # Intra range
EndFlanksGenRange <- flank(GenRange, 2, start = F) # получение набора обрезков длиной 2 После краёв блоков (исключены статы)
  ranges(EndFlanksGenRange); ranges(GenRange) # сравнение
ShiftedGenRange <- shift(GenRange, 10) # получение смещённого на 10 GRanges
  ranges(ShiftedGenRange); ranges(GenRange) # смещение на + 10
ResizedGenRange <- resize(GenRange, 5, fix = 'start') # переприсвоение длины, стартовая точка задана стартом (это дефолт)
  ranges(ResizedGenRange); ranges(GenRange) # все сжались с конца
    # Inter range, имена стираются
reduce(GenRange) # слияние пересекающихся диапазонов, много внутренних параметров, чувствителен к цепи
Gaps <- gaps(GenRange) # получение набора всех непокрытых участков, требует длин хромосом, возвращает всю цепь-связанную кашу
  ranges(Gaps) # эээ... гэпы с цепью?...
DisjoinedGenRange <- disjoin(GenRange) # выщепление пересекающихся участков блоков
  ranges(DisjoinedGenRange); ranges(GenRange) # непересекающиеся оставлены, 2 пересекающихся расщеплены на непересекающиеся куски и 1 участок пересечения
Coverage <- coverage(GenRange) # получение зачений покрытия в данном случае по хромосомам в виде RleList
 
### манипуляции со многими GRanges
union(GenRange, GenRange2) # общие блоки
intersect(GenRange, GenRange2) # пересекающиеся блоки
setdiff(GenRange, GenRange2) # отличающеся блоки
 
### получение и работа с объектами GenomicRangeList
  # сборка с нуля
gr1 <- GRanges(seqnames = "chr1", ranges = IRanges(2, 4), strand = "+", score = 50, lolID = 'o'); seqlengths(gr1) <- 1e6
gr2 <- GRanges(seqnames = "chr1", ranges = IRanges(2, 6), strand = "*", score = 20, lolID = 'm'); seqlengths(gr2) <- 1e6
GRangeLN <- GRangesList("ex1" = gr1, "ex2" = gr2)
GRangeLN
 
  # Расщепление GRanges
GRangeL <- split(GenRange, rep(1:2, each = 2)) # получение GenomicRagneList объекта из 2 элементов по 2 блока в каждом
GRangeL[2] # второй элемент, 
GRangeL[[2]] # он-же, просто корректнее
GRangeL[1:2, "score"] # поднабор элементов и колонок метаданных
GRange <- c(GRangeL[[1]], GRangeL[[2]]) # склеиваем обратно в объект GRanges
GRange == GenRange #[1]  TRUE TRUE TRUE TRUE # ОНО!
 
  # базовый ввод-вывод идентичен GRanges кроме
elementLengths(GRangeLN)
GRangeLN$ex1
GRangeLN["ex2"]
  # манипуляции с диапазонами аналогичны
 
  # смешанная работа
ChromRangeL <- split(GenRange, names(GenRange)) # расщепление GRange в GRangeList по хромосомам
addGR <- GenRange # копирование
addGR[3] <- ChromRangeL[[1]] # переприсвоение третьего блока, элемент GRL нормально присваивается для GR
addGR # имена не меняются, переписан 3 блок, содержимое как у 1го
 
### эквивалентные циклам операции
sapply(GRangeLN, length) # возвращается вектор длин
lapply(GRangeLN, length) # возвращается список длин
mapply(c, GRangeLN, GRangeL) # не работает
Map(c, GRangeLN, GRangeL) # не работает
endoapply(c, GRangeLN, GRangeL) # не работает
mendoapply(c, GRangeLN, GRangeL) # не работает
 
### пересечения интервалов
overlaps <- findOverlaps(gr2, GRangeLN) # нахождения пересечений блоков в диапазонах и списках диапазонов, класс Hits
countOverlaps(gr1, GRangeLN) # подсчёт пересечений, возврат вектора
subsetByOverlaps(gr1, GRangeLN) # возврат пересечения в виде GRanges
match(GRangeLN, gr1) # устарело
gr1 %in% GRangeLN # устарело

Использование Rsamtools

Дисклеймер: это описание использования пакета Rsamtools на фоне реального анализа данных Ribosome Profiling

Адаптация samtools для R является основным средством импорта данных картирования NGS в R. Сопряжение со статистическим языком обеспечивает функционал, недоступный для самого samtools.

Сделано несколько методов чтения BAM файлов, различаются и возвращаемые ими объекты, но общий ввод основан на указании референсных файлов или содержащих множество BAM файлов  папок - объекты типа BamFile.
Метод scanBam возвращает (загруженный в память) объект класса list со своей сложной структурой,
метод readBamGappedAlignments из группы методов ReadBamGapped возвращает класс GenomicRanges (точнее GappedAlignments), что из пакета GenomicRanges, а не Rsamtools,
близкий ему метод readBamGappedAlignmentPairs возвращает класс GappedAlignmentPairs из пакета GenomicRanges,
метод readBamGappedReads принадлежит пакету ShortRead и возвращает класс GappedReads.

Главный метод чтения Rsamtools это scanBam с использованием sacnBamParam - указание диапазонов выборочного чтения BAM файла (вот они - Big Data): указание which - где читать, задаются Rle диапазоны и what - что читать, какие колонки из BAM файла оставлять. Серии файлов составляются с помощью list.files с использованием паттернов, не все методы Rsamtools поддерживают чтение серий файлов. Примечательно что методы GenomicRanges типа summarize работают и для list BAM и для серий BAM файлов.

library(Rsamtools)
### задание референса BamFile
mRNAalignd <- "/media/data/profiling_data_processing/mRNA-exp.sort.bam" # выровненные риды общего транскриптома
mRNAindx <- "/media/data/profiling_data_processing/mRNA-exp.sort.bam.bai" # индекс выравнивания
filebam <- BamFile(file = mRNAalignd, index = mRNAindx) # объект-ссылка типа BamFile
### перебор разных методов чтения BAM файлов
scannedBam <- scanBam(filebam) # файл в виде списка, тормозит, жрёт много памяти, зато загружает ВСЁ в виде иерархического списка
readedGAlign <- readBamGappedAlignments(filebam) # файл в виде GRanges с 8 колонками без метаданных (цепь-длина-старт-конец, гэпы и т.д.)
readedGAlignPair <- readBamGappedAlignmentPairs(filebam) # этот метод для paired-end! он не прочитал нормально файл, но должен быть аналог для paired
readedGappedReads <- readBamGappedReads(filebam) # не видно разницы между GappedAlignments - такой-же GRanges
### рассмотрение и преобразование основного формата хранения BAM Rsamtools
sapply(scannedBam[[1]], class) # получение типов элементов элемента листа
lst <- lapply(names(scannedBam[[1]]), function(elt) { # применение функции к элементу списка с выводом в результате списка
do.call(c, unname(lapply(scannedBam, "[[", elt))) # развоваричвание одного элемента списка в список со многимим элементами
}
)
names(lst) <- names(scannedBam[[1]]) # задание имён элементов списка
dataBam <- do.call("DataFrame", lst) # трансформация в data frame, очень ресурсоёмко
### фрагментарное чтение BAM файла
testRangeList <- RangesList(NC_000913.2 = IRanges(50000, 100000)) # пробный диапазон в виде RangeList с одним элементом, имя должно совпадать с хромосомой
partBamParam <- ScanBamParam(what = c("rname", "strand","pos"), which = testRangeList) # сборка параметров для чтения BAM
partBam <- scanBam(filebam, param = partBamParam) # читается
### чтение серий файлов и summarize из GenomicRanges
fls <- list.files(path = "/media/data/profiling_data_processing", pattern = "*mRNA.*sort.bam$", full = T) # создание списка файлов
BamFilesList <- BamFileList(fls, index = character()) # сборка листа файлов
testRange <- GRanges(seqnames = 'NC_000913.2', ranges = IRanges(start = 500000, end = 1000000), strand = '+') # тестовый GRange, GRL почему-то не заработал
summary <- summarizeOverlaps(testRange, BamFilesList) # саммари в объект
result <- assays(summary) # неудачный тестовый GRanges, ничего просяняющего
### остальные фукции Rsamtools - стандартные манипуляции с BAM и VCF файлами и построение референсов - BamViews

Colored with dumpz.org