1

我有一个包含大量读取的 BAM 文件。scanBam我可以使用from将它加载到 R 中Rsamtools。

但是，我只需要读取的一个子集。我有一个character带有我感兴趣的 qnames 的向量。

scanBam返回一个包含 1 个元素的列表，该列表包含 13 个元素，其中包含所有数千次读取的数据。

如何通过qname保留结构来子集这个对象？我无法在手册或网上找到任何内容。

4

2 回答 2

2

使用 GenomicAlignments::readGAlignments 输入数据可能更方便，包括通过指定 param=ScanBamParam(what="qname") 作为参数的 qname。然后，您可以使用 %in% 作为子集。这是一个更完整的示例，使用ExperimentData包之一

library(GenomicAlignments)
library(RNAseqData.HNRNPC.bam.chr14)

fname <- RNAseqData.HNRNPC.bam.chr14_BAMFILES[1]    
want <- c("ERR127306.11930623", "ERR127306.24720935",
    "ERR127306.23706011", "ERR127306.22418829", "ERR127306.13372247",
    "ERR127306.20686334", "ERR127306.11412145", "ERR127306.4711647",
    "ERR127306.7479582", "ERR127306.12737243")
aln <- readGAlignments(fname, param=ScanBamParam(what="qname"))
aln[mcols(aln)$qname %in% want]

BAM 文件当然很大，而 qnames 是其中很大一部分；以块的形式遍历文件通常是有意义的。这通过 yieldReduce 启用（在当前的 Rsamtools 中），其中向 BamFile 提供设置为合理（例如，1M）读取次数的 yieldSize，一个 MAP 函数来输入一大块数据并处理它（例如，过滤不需要的读取)，一个（可选）REDUCE 函数来连接结果，以及一个（可选）DONE 函数来指示迭代何时完成。解决方案看起来像（yieldSize 是人为的小，以便用示例数据进行说明）：

bfl <- BamFile(fname, yieldSize=100000)  ## larger, e.g., 1M-5M
MAP <- function(bfl, want) {
    ## message("tick")
    aln <- readGAlignments(bfl, param=ScanBamParam(what="qname"))
    if (length(aln) == 0)
        NA                          # semaphore -- DONE
    else
        aln[mcols(aln)$qname %in% want]
}
REDUCE <- c
DONE <- function(x) identical(x, NA)
result <- yieldReduce(bfl, MAP, REDUCE, DONE, want=want)

可以使用 scanBam 采用类似的方法，但数据结构（列表列表）处理起来更复杂：

x <- scanBam(fname, param=ScanBamParam(what=c("qname", "pos")))
keep <- lapply(lapply(x, "[[", "qname"), "%in%", want)
result <- Map(function(elts, keep) {
    lapply(elts, "[", keep)
}, x, keep)

这也可以与 yieldReduce 一起使用。

如果您有兴趣使用过滤后的读取创建一个新的bam 文件，那么

filter_factory <- function(want) {
    list(KeepQname = function(x) x$qname %in% want)
}
filter <- FilterRules(filter_factory(want))
dest <- filterBam(fname, tempfile(), filter=filter,
                  param=ScanBamParam(what="qname"))
readGAlignments(dest)
于 2014-06-13T13:11:43.823 回答
1

我最终使用了subset(DataFrame(scanBam(bam_file)[[1]]),qname %in% qname_vector). 这不会保持完全相同的结构（列表列表），但所有信息都保留并易于访问。

于 2014-06-13T15:53:20.233 回答