【直播】我的基因组 33：用samstat软件对sam文件做统计

在此之前，我不止一次强调过QC的重要性，对全基因组测序等以找variation为主的分析流程来说，不仅仅是对测序数据的QC，还有比对之后的sam/bam文件也需要QC，最后找出的variation文件也需要QC。在我们生信技能树的论坛上面有详细的说明，见： () 。之前曾耗费了10讲来解析sam格式，就是基于我个人的知识背景来对比对结果进行QC，而事实上，有很多成熟的软件就可以完成全套流程，当然，所谓的全套也只是基于软件作者的知识背景。我这里讲挑出两个读者来信咨询的最多的软件来简单讲解一下吧！

samstat

这个软件对大的bam文件运行经常会报错，就是程序界最出名的segment fault，应该是内存不够。正好，我们把这个大的bam文件根据染色体拆分成了小的bam，可以批量运行一下samstat，对比对文件进行一些统计。我在博客里面写过这个软件的教程，请点击阅读原文查看，或者复制链接（ ）到浏览器查看SAMStat软件使用说明书。

这是一个C源码程序，所以安装的时候需要./configure,make,make install 三步曲即可，如果没有root权限就指定一个环境变量来编译这个源码，具体参照博客，如果报错自己学会解决！

用法非常简单，因为这个软件做得事情不多，也就是跟我们之前耗费了10讲来解析sam格式的比对文件差不多！

命令如下：

  samstat P_jmzeng.final.REF_*bam

默认对每一个输入的bam文件，都是会输出一个网页版的统计QC报告的，上面的命令会把所有染色体的bam文件都输入，但事实上这个软件对某几条染色体还是有限制，可以换成是

ls P_jmzeng.final.REF_*bam  |xargs samstat  

来运行，避免某些染色体报错，后面的都不能运行了。

但事实上，这个代码仍然会因为某个染色体的bam文件太大而中断，不知道shell里面有没有try这个类似功能的函数。

我简单Google了一下，shell的确没有这个功能，但是有程序猿给出了模拟的解决方案，我测试了一下：

ls P_jmzeng.final.REF_*bam  |while read id

do

{ # try

samstat  $id

} || { # catch

# save log for exception

echo "This bam file is too big : $id  "

}

done

还不错，可以把所有文件该跑出来的结果都跑出来了，结果没什么好解读的，这个软件其实很烂的，我一般都不用了。

samtools stats 

而且samtools更新的版本本身就包括了一系列的统计和绘图功能！

现在比较好的就是：

samtools stats  test.bam  >test.stats

plot-bamstats -p test  test.stats

因为我安装samtools的时候把它添加到了环境变量，所以我不需要用全路径来使用该程序，直接调用即可，plot-bamstats就是安装samtools软件时候附带的一个功能。

有时候，觉得生信工程师真心不容易，什么都得学一点什么都得会一点，但是最重要的，是要有解决问题的勇气和感觉。(现在划重点，本讲最重要的不是samstat软件怎么用，因为实在是太烂了这个软件，也不是samtools的stats功能怎么用，虽然出图很漂亮，统计的也很清晰。)

下一讲我们来重点谈谈Qualimap这个软件。

文：Jimmy、阿尔的太阳

图文编辑：吃瓜群众