R中的子集SAM / BAM文件答案

【问题标题】：Subset SAM/BAM file in RR中的子集SAM / BAM文件
【发布时间】：2020-09-15 03:29:41
【问题描述】：

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

但是，我只需要读取的一个子集。我有一个 character 向量，其中包含我感兴趣的 qname。

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

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

【问题讨论】：

标签： r subset bioinformatics bioconductor bam

【解决方案1】：

使用 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)

【讨论】：

非常感谢您的详细回答。你所有的建议似乎都是我想要的。当我决定必须有更好的方法时，第三个实际上是我尝试的工作版本。您的第二个解决方案非常复杂，我只需要处理一个 BAM 文件，所以我不必费心等待 30 秒，但很高兴知道未来的工作。您的最后一个建议是我在手册中找到的，但我明确不想写入新文件（在这种情况下，我会直接使用 samtools 而不是编写 R 脚本）。您能否评论一下您的第一个建议相对于我的简单解决方案的好处？
@sg-lecram 你的解决方案很好；主要区别在于 GAlignments 知道范围，因此，例如，如果您对特定的感兴趣区域进一步感兴趣 roi = GRanges(...) 它会直接到 subsetByOverlaps(aln, roi) 或 countOverlaps() 或任何其他非常方便的基于范围的操作。使用 yieldReduce 可能没有真正的成本，大部分时间都花在输入（无论如何都必须完成）和下游处理上； yieldReduce 对于管理内存非常有用，尤其是。结合多个文件的并行处理。

【解决方案2】：

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

【讨论】：