基因组覆盖率滑动 window

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 的名称（因为它是一个 multi-contig 文件，此 ID 会更改）- 位置（碱基）- 映射的读取数（覆盖率）。

现在我想计算滑动中的覆盖率（第三列）windows； window 大小为 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 和 window，它是计算出来的：

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 的形式从分组操作中 return 任意数量的值。在这种情况下，我们 return 覆盖值的滚动平均值，以及每个 window 中 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 倍。

基因组覆盖率滑动 window

Genome coverage as sliding window

r

bioinformatics

rollapply