項目背景
對bam文件按多種自定義的劃窗來統(tǒng)計區(qū)間的測序覆蓋度。
我的需求
對于同一個bam文件剃氧，根據(jù)不同分割條件產生的多個bed文件敏储，循環(huán)批量計算測序覆蓋度。
比如朋鞍，輸入為一個排序后的bam文件（case.sort.bam）已添，以及對基因組進行不同方式的自定義分箱而產生的多個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文件已構建索引
samtools index xxx.sort.bam
bedtools coverage循環(huán)時出現(xiàn)問題
# 循環(huán)時只有一個樣本可以成功產生結果
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中沒有報錯
我的嘗試
# 不循環(huán)時是正常的
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
# 把每個樣本的各種分箱方式都單獨運行了一遍，全部都能產生正常結果
為什么單獨運行時正常禾怠，循環(huán)時卻只有一個結果正常呢返奉？
換用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
# 加上循環(huán)
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產生的結果只有區(qū)間內的reads數(shù)，而且和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 計算得到的結果差異非常大吗氏！
查了一下之后發(fā)現(xiàn)芽偏，samtools bedcov的算法是把區(qū)間內每個堿基的測序深度加起來，而bedtools coverage是直接計算區(qū)間內的reads數(shù)（見https://www.biostars.org/p/195497/）弦讽，兩者適用于不同的需求污尉。
那么，我的需求基本被解決了往产，因為兩種算法都能夠表現(xiàn)測序深度被碗，除以區(qū)間長度就得到了平均覆蓋率。我的最終目的是比較每個區(qū)間的測序深度的相對大小仿村，找到平均覆蓋率顯著增大的區(qū)間锐朴。
延伸思考
但是如果我還是想知道開頭的那個問題：對同一個樣本用不同的分箱方式批量計算區(qū)間內的reads數(shù)，該怎么實現(xiàn)呢蔼囊？
另外焚志，bedtools coverage 和 samtools bedcov 這兩個工具都不能選擇線程數(shù)衣迷，也不能選擇-O輸出，用起來略有些不舒服酱酬。所以壶谒，還有其他更合適的工具推薦嗎？
補充說明
有朋友說用bedtools multicov就可以解決問題膳沽，那是誤解我的需求了佃迄。bedtools multicov可以解決多個bam文件+一個bed文件情況下的批量運行問題，但是我的問題是一個bam文件+多個bed文件情況下如何批量運行贵少？
由于我的shell腳本能力還不太過關呵俏，自己調試花了大半天也沒有找到問題原因，所以發(fā)出來請大家一起看看（當眾處刑hhh)
期待有大神能給我提示（感激~）
