項目背景
對bam文件按多種自定義的劃窗來統(tǒng)計區(qū)間的測序覆蓋度。
我的需求
對于同一個bam文件，根據(jù)不同分割條件產(chǎn)生的多個bed文件，循環(huán)批量計算測序覆蓋度。
比如，輸入為一個排序后的bam文件（case.sort.bam），以及對基因組進行不同方式的自定義分箱而產(chǎn)生的多個bed文件（100M.bed, 50M_10M.bed, 20M_5M.bed），希望對bed文件做循環(huán)，批量得到不同分箱條件下同一個bam文件的自定義區(qū)間reads數(shù)統(tǒng)計。
準備工作
# 使用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文件已構(gòu)建索引
samtools index xxx.sort.bam
bedtools coverage循環(huán)時出現(xiàn)問題
# 循環(huán)時只有一個樣本可以成功產(chǎn)生結(jié)果
cat window_list.txt|while read window;do (nohup bedtools coverage -a $window.bed -b $sample.sort.bam > $sample.$window.counts.txt.tmp &) ;done
# 產(chǎn)生三個結(jié)果文件，其中兩個文件大小為0，里面沒有內(nèi)容
-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中沒有報錯
我的嘗試
# 不循環(huán)時是正常的
window=50M_10M
sample=case
bedtools coverage -a $window.bed -b $sample.sort.bam > $sample.$window.counts.txt
# 產(chǎn)生結(jié)果：
-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 &
# 產(chǎn)生結(jié)果：
-rw-rw-r-- 1 xpan xpan 18K Aug 10 16:52 case.50M_10M.counts.txt.tmp
# 把每個樣本的各種分箱方式都單獨運行了一遍，全部都能產(chǎn)生正常結(jié)果
為什么單獨運行時正常，循環(huán)時卻只有一個結(jié)果正常呢？
換用samtools bedcov后成功了
# 換用samtools bedcov再嘗試同樣的操作
samtools bedcov $window.bed $sample.sort.bam > $sample.$window.counts.txt.tmp1
# 產(chǎn)生結(jié)果：-rw-rw-r-- 1 xpan xpan 1.3K Aug 10 15:47 case.100M.counts.txt.tmp1
# 加上循環(huán)
cat window_list.txt|while read window;do (nohup samtools bedcov $window.bed $sample.sort.bam > $sample.$window.counts.txt.tmp1 &) ;done
# 產(chǎn)生三個結(jié)果都正常：
-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產(chǎn)生不同結(jié)果
# samtools bedcov產(chǎn)生的結(jié)果只有區(qū)間內(nèi)的reads數(shù)，而且和bedtools coverage計算的結(jié)果出入還不小
# samtools bedcov結(jié)果取前三行展示：
chr1 0 100000000 22300197
chr1 100000000 200000000 22302389
chr1 200000000 248956422 11554278
# bedtools coverage結(jié)果取前三行展示：
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 計算得到的結(jié)果差異非常大！
查了一下之后發(fā)現(xiàn)，samtools bedcov的算法是把區(qū)間內(nèi)每個堿基的測序深度加起來，而bedtools coverage是直接計算區(qū)間內(nèi)的reads數(shù)（見https://www.biostars.org/p/195497/），兩者適用于不同的需求。
那么，我的需求基本被解決了，因為兩種算法都能夠表現(xiàn)測序深度，除以區(qū)間長度就得到了平均覆蓋率。我的最終目的是比較每個區(qū)間的測序深度的相對大小，找到平均覆蓋率顯著增大的區(qū)間。
延伸思考
但是如果我還是想知道開頭的那個問題：對同一個樣本用不同的分箱方式批量計算區(qū)間內(nèi)的reads數(shù)，該怎么實現(xiàn)呢？
另外，bedtools coverage 和 samtools bedcov 這兩個工具都不能選擇線程數(shù)，也不能選擇-O輸出，用起來略有些不舒服。所以，還有其他更合適的工具推薦嗎？
補充說明
有朋友說用bedtools multicov就可以解決問題，那是誤解我的需求了。bedtools multicov可以解決多個bam文件+一個bed文件情況下的批量運行問題，但是我的問題是一個bam文件+多個bed文件情況下如何批量運行？
由于我的shell腳本能力還不太過關(guān)，自己調(diào)試花了大半天也沒有找到問題原因，所以發(fā)出來請大家一起看看（當眾處刑hhh)
期待有大神能給我提示（感激~）
