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

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

測序深度和覆蓋度的示意圖如下：

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

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

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

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 ... & ...為命令 表明命令后臺執(zhí)行  也可以

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

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

這個腳本運(yùn)行會比較慢,，因為是針對整個55G的bam文件。耗時如下：

real130m34.855s

user237m14.308s

sys1m45.692s

結(jié)果如下：

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

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

先不要急著怪公司,，請聽下回分解