项目背景
对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)
期待有大神能给我提示（感激~）
网友评论