提问|如何批量统计区间测序覆盖度？

2019-08-10 本文已影响0人 Steven潘

项目背景

对bam文件按多种自定义的划窗来统计区间的测序覆盖度。

我的需求

对于同一个bam文件，根据不同分割条件产生的多个bed文件，循环批量计算测序覆盖度。

比如，输入为一个排序后的bam文件（case.sort.bam），以及对基因组进行不同方式的自定义分箱而产生的多个bed文件（100M.bed, 50M_10M.bed, 20M_5M.bed），希望对bed文件做循环，批量得到不同分箱条件下同一个bam文件的自定义区间reads数统计。

准备工作

# 使用bedtools对基因组染色体坐标进行窗口划分
# -w定义窗口长度  -s定义划窗步长
bedtools  makewindows -g sizes.genome  -w  100000000 > 100M.bed
bedtools  makewindows -g sizes.genome  -w  50000000 -s 10000000 > 50M_10M.bed
bedtools  makewindows -g sizes.genome  -w  20000000 -s 5000000  > 20M_5M.bed

# bam文件已构建索引
samtools index xxx.sort.bam

bedtools coverage循环时出现问题

# 循环时只有一个样本可以成功产生结果
cat window_list.txt|while read window;do (nohup bedtools coverage -a $window.bed -b $sample.sort.bam > $sample.$window.counts.txt.tmp &) ;done
# 产生三个结果文件，其中两个文件大小为0，里面没有内容
-rw-rw-r-- 1 xpan xpan    0 Aug 10 15:31 case.50M_10M.counts.txt.tmp
-rw-rw-r-- 1 xpan xpan    0 Aug 10 15:31 case.20M_5M.counts.txt.tmp
-rw-rw-r-- 1 xpan xpan 2.3K Aug 10 15:32 case.100M.counts.txt.tmp
# nohup.out中没有报错

我的尝试

# 不循环时是正常的
window=50M_10M
sample=case
bedtools coverage -a $window.bed -b $sample.sort.bam > $sample.$window.counts.txt 
# 产生结果：
-rw-rw-r-- 1 xpan xpan  18K Aug 10 11:50 case.50M_10M.counts.txt

# 单独测试nohup & 语句，也是正常的
nohup bedtools coverage -a $window.bed -b $sample.sort.bam > $sample.$window.counts.txt.tmp &
# 产生结果：
-rw-rw-r-- 1 xpan xpan  18K Aug 10 16:52 case.50M_10M.counts.txt.tmp
# 把每个样本的各种分箱方式都单独运行了一遍，全部都能产生正常结果

为什么单独运行时正常，循环时却只有一个结果正常呢？

换用samtools bedcov后成功了

# 换用samtools bedcov再尝试同样的操作
samtools bedcov $window.bed $sample.sort.bam > $sample.$window.counts.txt.tmp1
# 产生结果：-rw-rw-r-- 1 xpan xpan 1.3K Aug 10 15:47 case.100M.counts.txt.tmp1
# 加上循环
cat window_list.txt|while read window;do (nohup samtools bedcov $window.bed $sample.sort.bam > $sample.$window.counts.txt.tmp1 &) ;done
# 产生三个结果都正常：
-rw-rw-r-- 1 xpan xpan 1.3K Aug 10 15:57  case.100M.counts.txt.tmp1
-rw-rw-r-- 1 xpan xpan  20K Aug 10 15:57  case.20M_5M.counts.txt.tmp1
-rw-rw-r-- 1 xpan xpan  11K Aug 10 15:57  case.50M_10M.counts.txt.tmp1

samtools与bedtools产生不同结果

# samtools bedcov产生的结果只有区间内的reads数，而且和bedtools coverage计算的结果出入还不小
# samtools bedcov结果取前三行展示：
chr1    0   100000000   22300197
chr1    100000000   200000000   22302389
chr1    200000000   248956422   11554278

# bedtools coverage结果取前三行展示：
chr1    0   100000000   182139  5952928 100000000   0.0595293
chr1    100000000   200000000   180032  5266467 100000000   0.0526647
chr1    200000000   248956422   94232   3098065 48956422    0.0632821

samtools bedcov 与 bedtools coverage 计算得到的结果差异非常大！

查了一下之后发现，samtools bedcov的算法是把区间内每个碱基的测序深度加起来，而bedtools coverage是直接计算区间内的reads数（见https://www.biostars.org/p/195497/），两者适用于不同的需求。

那么，我的需求基本被解决了，因为两种算法都能够表现测序深度，除以区间长度就得到了平均覆盖率。我的最终目的是比较每个区间的测序深度的相对大小，找到平均覆盖率显著增大的区间。

延伸思考

但是如果我还是想知道开头的那个问题：对同一个样本用不同的分箱方式批量计算区间内的reads数，该怎么实现呢？

另外，bedtools coverage 和 samtools bedcov 这两个工具都不能选择线程数，也不能选择-O输出，用起来略有些不舒服。所以，还有其他更合适的工具推荐吗？

补充说明

有朋友说用bedtools multicov就可以解决问题，那是误解我的需求了。bedtools multicov可以解决多个bam文件+一个bed文件情况下的批量运行问题，但是我的问题是一个bam文件+多个bed文件情况下如何批量运行？

由于我的shell脚本能力还不太过关，自己调试花了大半天也没有找到问题原因，所以发出来请大家一起看看（当众处刑hhh)

期待有大神能给我提示（感激~）