【直播】我的基因组（19）:根据比对结果来统计测序深度和覆盖度

看来本次直播我的基因组分析流程效果还不错，不少朋友跟着自己动手开始分析全基因组测试数据了，值得表扬。其中有好几个朋友留言向我反映公司给的统计报告里面的覆盖度的问题，如下图：

我在第 8讲中写道，每条染色体的覆盖度都接近于100%，而且测序深度大多在40X以上，不少读者表示看晕了，明明这个条形图显示覆盖度不到60%呀！其实是公司的这个图没有做好，它里面的覆盖度用的是最上面的曲线显示的，条形图是测序深度！！！

测序深度和覆盖度的示意图如下：

那么我们就来自己动手统计一下比对好的sam/bam文件的测序深度和覆盖度吧！

这个统计主要依赖于samtools的depth功能，或者说mpileup功能，输入文件都是sort好bam格式的比对文件。事实上，你如果读samtools的源代码就会发现，其实depth功能调用的就是mpileup的函数。但是mpileup可以设置一系列的过滤参数。而depth命令是纯天然的，所以mpileup的结果一定会小于depth的测序深度。对mpileup，可以不选择-u -f 参数指定参考基因组，因为我们只需要测序深度情况，还有，可以指定-q 1 来过滤掉多比对情况。还可以用-Q来过滤低质量的碱基(base pair),用-A来过滤无法定位的reads，结果如下：

针对这个全基因组位点的统计结果，我们很容易写脚本来计算每条平均的测序深度和各个染色体的覆盖度。

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 ... & ...为命令表明命令后台执行也可以

nohup ... & > *.log 将运行结果写入到一个log文件里面

time 命令可以统计命令运行的时间

这个脚本运行会比较慢，因为是针对整个55G的bam文件。耗时如下：

real 130m34.855s

user 237m14.308s

sys 1m45.692s

结果如下：

其中chromosome的长度在bam文件里面可以看到，用samtools view -H P_jmzeng.final.bam 即可！！！把上面的表格可视化，就是文章最开头的figure。但是很明显公司给我的各个指标均高于我自己算出来的！尤其是其中几个染色体的coverage非常差。当然平均测序深度，我在这里选择的是整个染色体的长度作为分母，对于那些覆盖度很差的染色体，平均测序深度就被被拉低。所以我起初猜想是不是问题出在我用的是samtools的mpileup命令，而不是depth命令！（因为我以前从来没有如此细致的去比较它们的差别，其实这这个命令的确有差别，但是对全基因组层面的统计指标几乎没什么影响）

下面的表格，是我用depth命令的结果，而且平均测序深度我选择被覆盖到的染色体长度作为分母，所以看到测序深度有些许提升，但是有几条染色体的覆盖度堪忧。与公司给我的报告不符合！