基因组覆盖作为滑动窗口

Question

我使用bwa mem 算法将我的读取映射到我的程序集，并使用samtools depth 提取每个碱基的读取数（= 覆盖率）。生成的文件如下：

1091900001  1   236
1091900001  2   245
1091900001  3   265
1091900001  4   283
1091900001  5   288
1091900002  1   297
1091900002  2   312
1091900002  3   327
1091900002  4   338
1091900002  5   348

包含三列：contig 的名称（因为它是一个多 contig 文件，所以此 ID 会发生变化）- 位置（碱基）- 映射的读取数（覆盖率）。

现在我想计算滑动窗口的覆盖率（第三列）；窗口大小为 3，滑动量为 2 作为平均值 - 每个重叠群（第一列）。

我想使用zoo包的rollapply函数。

require(zoo)
cov <- read.table("file",header=FALSE, sep="\t", na.strings="NA", dec=".", strip.white=TRUE)
library(reshape) #loads the library to rename the column names
cov<-rename(cov,c(V1="Chr", V2="locus", V3="depth")) #renames the header
rollapply(cov$depth, width = 3, by = 2, FUN = mean, align = "left")

但这当然没有考虑重叠群。另外，我的预期输出应该包括 contig-info 和窗口，它是计算出来的：

1091900001  1   3   248.6667
1091900001  3   5   278.6667
1091900002  1   3   312.0000
1091900002  3   5   337.6667

在R 中是否有一种简单的方法可以做到这一点？

Answer 1

以下是使用dplyr 函数group_by 和do 执行此操作的方法：

library(dplyr)

cov %>% 
  group_by(Chr) %>% 
  do(
    data.frame(
      window.start = rollapply(.$locus, width=3, by=2, FUN=min, align="left"),
      window.end = rollapply(.$locus, width=3, by=2, FUN=max, align="left"),
      coverage = rollapply(.$depth, width=3, by=2, FUN=mean, align="left")
      )
    )

# # A tibble: 4 x 4
# # Groups:   Chr [2]
#          Chr window.start window.end coverage
#        <int>        <int>      <int>    <dbl>
# 1 1091900001            1          3 248.6667
# 2 1091900001            3          5 278.6667
# 3 1091900002            1          3 312.0000
# 4 1091900002            3          5 337.6667

do 允许您以 data.frame 的形式从分组操作中返回任意数量的值。在这种情况下，我们返回覆盖率值的滚动平均值，以及来自每个窗口的 locus 中的 min 和 max 值。

编辑：

如果您的数据集很大，您最好使用data.table 执行计算。如果您以前没有见过它的语法，它的语法有点难以理解，但它可以在更大数据的分组操作中提供显着的速度改进。以下是您使用data.table 进行操作的方式：

library(data.table)    

setDT(cov)
cov[, .(
      window.start = rollapply(locus, width=3, by=2, FUN=min, align="left"),
      window.end = rollapply(locus, width=3, by=2, FUN=max, align="left"),
      coverage = rollapply(depth, width=3, by=2, FUN=mean, align="left")
      ),
    .(Chr)]

根据您提供的示例行，以下是dplyr 和data.table 方法的基准测试结果（以毫秒为单位）：

# dplyr:
      min       lq     mean   median       uq      max neval
 7.811753 8.685976 10.10268 9.243551 10.42691 144.5274  1000

# data.table:
      min       lq     mean  median       uq      max neval
 1.924472 2.105459 2.510832 2.30479 2.685706 8.848451  1000

因此，在示例数据中，data.table 选项平均快了大约 4 倍。