读取带有标签的 BAM 到 tbl / 递归 dplyr::bind_cols 列表列表

Question

我想读入包含一些标签的 BAM 文件，然后将其转换为 tibble 以供进一步处理。

通常，这可以很简单地实现：

library(Rsamtools)
library(tidyverse)
map.info <- c("rname", "strand", "pos")
map.params <- ScanBamParam(what = map.info)
bam <- scanBam(bam.file, param = map.params)

scanBam returns 包含命名向量 rname、strand 和 pos 的列表，可以使用 dplyr::bind_cols(bam) 轻松连接。但是，假设我对 MD 标签感兴趣，我需要执行以下操作：

map.params <- ScanBamParam(what = map.info, tag = c("MD"))
bam <- scanBam(bam.file, param = map.params)

现在，bam 是一个列表的列表，具有命名向量 rname、strand 和 pos，但还有另一个 tag本身就是一个列表，其中有一个名为 MD.

的向量

dplyr::bind_cols 无法处理此嵌套列表列表，并引发错误，但 as.data.frame(bam) 有效。

TL;DR

一个玩具示例归结为核心问题：

> df.list <- list(a = 1:2, b = 3:4, c = 5:6)
> df.nest <- list(a = 1:2, b = 3:4, c = 5:6, d = list( e = 7:8 ))
> dplyr::bind_cols(df.list)
# A tibble: 2 x 3
       a     b     c
   <int> <int> <int>
 1     1     3     5
 2     2     4     6
> dplyr::bind_cols(df.nest)
Error in cbind_all(x) : Argument 4 must be length 2, not 1
> as.data.frame(df.nest)
  a b c e
1 1 3 5 7
2 2 4 6 8

有没有办法在嵌套列表中递归 bind_cols？

暂定答案

受 @mt1022 的回答启发，并且没有进一步检查 Rsamtools 基本代码，尽管格式与玩具示例非常相似，但 scanBam 输出似乎不像玩具示例.

但是，正如我们所知道的那样，以下内容也应该实现完全合并 tibble:

map.info <- c("rname", "strand", "pos")
map.params <- ScanBamParam(what = map.info, tag = c("MD", "NM"))
bam <- scanBam(bam.file, param = map.params)
bam.tbl <- bind_cols(do.call(bind_cols, bam[[1]][c("rname", "strand", "pos")]),
                     do.call(bind_cols, bam[[1]]$tag))

这似乎比我希望（或预期）的更老套，但它确实有效。

基准测试

有三个选项会导致类似的结果：

as.data.frame(bam)
bind_cols(do.call(bind_cols, bam[[1]][map.info]), do.call(bind_cols, bam[[1]]$tag))
bind_cols(lapply(bam, as.data.frame), .id = 'rn')

我认为包含 10762160 行的示例 bam 文件应该可以很好地说明哪种方法最有效。

> length(bam[[1]]$rname)
[1] 10762160
> system.time(for (i in 1:100) as.data.frame(bam))
   user  system elapsed
 70.565  25.821  96.699
> system.time(for (i in 1:100) bind_cols(do.call(bind_cols, bam[[1]][mapI]), do.call(bind_cols, bam[[1]]$tag)))
   user  system elapsed
  0.124   0.020   0.144
> system.time(for (i in 1:100) dplyr::bind_rows(lapply(bam, as.data.frame), .id = 'rn'))
   user  system elapsed
108.091  36.046 144.623

看起来 'clunky'，我猜对 bind_col 的嵌套调用是最快的！

Answer 1

受bind_cols手册中这个例子的启发：

# You can mix vectors and data frames:
bind_rows(
  c(a = 1, b = 2),
  data_frame(a = 3:4, b = 5:6),
  c(a = 7, b = 8)
)

我想我们可以试试：

do.call(bind_cols, df.nest)
# # A tibble: 2 x 4
#       a     b     c     e
#   <int> <int> <int> <int>
# 1     1     3     5     7
# 2     2     4     6     8

编辑：一个真实的例子

library(Rsamtools)
# generate sample data using internal bam file of Rsamtools
which <- RangesList(seq1=IRanges(1000, 2000),
                    seq2=IRanges(c(100, 1000), c(1000, 2000)))
what <- c("rname", "strand", "pos")
param <- ScanBamParam(which=which, what=what, tag = c("NM"))

bamFile <- system.file("extdata", "ex1.bam", package="Rsamtools")
bam <- scanBam(bamFile, param=param)

# convert bam to a data.frame
bam.df <- dplyr::bind_rows(lapply(bam, as.data.frame), .id = 'rn')

#                  rn rname strand  pos NM
# 1    seq1:1000-2000  seq1      +  970  0
# 2    seq1:1000-2000  seq1      +  971  0
# 3    seq1:1000-2000  seq1      +  972  0
# ......

读取带有标签的 BAM 到 tbl / 递归 dplyr::bind_cols 列表列表

Read BAM with tags to tbl / recursive dplyr::bind_cols for list of lists

r

bioinformatics

bioconductor

dplyr

purrr

TL;DR

暂定答案

基准测试

编辑：一个真实的例子