看來本次直播我的基因組分析流程效果還不錯(cuò)，不少朋友跟著自己動(dòng)手開始分析全基因組測(cè)試數(shù)據(jù)了，值得表揚(yáng)。其中有好幾個(gè)朋友留言向我反映公司給的統(tǒng)計(jì)報(bào)告里面的覆蓋度的問題，如下圖：

我在第 8講中寫道，每條染色體的覆蓋度都接近于100%，而且測(cè)序深度大多在40X以上，不少讀者表示看暈了，明明這個(gè)條形圖顯示覆蓋度不到60%呀！其實(shí)是公司的這個(gè)圖沒有做好，它里面的覆蓋度用的是最上面的曲線顯示的，條形圖是測(cè)序深度?。?！

測(cè)序深度和覆蓋度的示意圖如下：

那么我們就來自己動(dòng)手統(tǒng)計(jì)一下比對(duì)好的sam/bam文件的測(cè)序深度和覆蓋度吧！

這個(gè)統(tǒng)計(jì)主要依賴于samtools的depth功能，或者說mpileup功能，輸入文件都是sort好bam格式的比對(duì)文件。事實(shí)上，你如果讀samtools的源代碼就會(huì)發(fā)現(xiàn)，其實(shí)depth功能調(diào)用的就是mpileup的函數(shù)。但是mpileup可以設(shè)置一系列的過濾參數(shù)。而depth命令是純天然的，所以mpileup的結(jié)果一定會(huì)小于depth的測(cè)序深度。對(duì)mpileup，可以不選擇-u -f 參數(shù)指定參考基因組，因?yàn)槲覀冎恍枰獪y(cè)序深度情況，還有，可以指定-q 1 來過濾掉多比對(duì)情況。還可以用-Q來過濾低質(zhì)量的堿基(base pair),用-A來過濾無法定位的reads，結(jié)果如下：

針對(duì)這個(gè)全基因組位點(diǎn)的統(tǒng)計(jì)結(jié)果，我們很容易寫腳本來計(jì)算每條平均的測(cè)序深度和各個(gè)染色體的覆蓋度。

nohup time samtools mpileup P_jmzeng.final.bam |perl -alne '{$pos{$F[0]}++;$depth{$F[0]}+=$F[3]} END{print "$_\t$pos{$_}\t$depth{$_}" foreach sort keys %pos}' 1>coverage_depth.txt 2>coverage_depth.log &

nohup ... & ...為命令 表明命令后臺(tái)執(zhí)行  也可以

nohup ... & > *.log 將運(yùn)行結(jié)果 寫入到一個(gè)log文件里面

time 命令 可以統(tǒng)計(jì) 命令 運(yùn)行的時(shí)間

這個(gè)腳本運(yùn)行會(huì)比較慢，因?yàn)槭轻槍?duì)整個(gè)55G的bam文件。耗時(shí)如下：

real130m34.855s

user237m14.308s

sys1m45.692s

結(jié)果如下：

其中chromosome的長度在bam文件里面可以看到，用samtools view -H P_jmzeng.final.bam 即可！??！把上面的表格可視化，就是文章最開頭的figure。但是很明顯公司給我的各個(gè)指標(biāo)均高于我自己算出來的！尤其是其中幾個(gè)染色體的coverage非常差。當(dāng)然平均測(cè)序深度，我在這里選擇的是整個(gè)染色體的長度作為分母，對(duì)于那些覆蓋度很差的染色體，平均測(cè)序深度就被被拉低。所以我起初猜想是不是問題出在我用的是samtools的mpileup命令，而不是depth命令?。?span style="color: rgb(0, 0, 0);">因?yàn)槲乙郧皬膩頉]有如此細(xì)致的去比較它們的差別，其實(shí)這這個(gè)命令的確有差別，但是對(duì)全基因組層面的統(tǒng)計(jì)指標(biāo)幾乎沒什么影響）

下面的表格，是我用depth命令的結(jié)果，而且平均測(cè)序深度我選擇被覆蓋到的染色體長度作為分母，所以看到測(cè)序深度有些許提升，但是有幾條染色體的覆蓋度堪憂。與公司給我的報(bào)告不符合！

先不要急著怪公司，請(qǐng)聽下回分解