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

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

samstat

这个软件对大的bam文件运行经常会报错,就是程序界最出名的segment fault,应该是内存不够。正好,我们把这个大的bam文件根据染色体拆分成了小的bam,可以批量运行一下samstat,对比对文件进行一些统计。我在博客里面写过这个软件的教程,请点击阅读原文查看,或者复制链接(http://www.bio-info-trainee.com/751.html )到浏览器查看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、阿尔的太阳

图文编辑:吃瓜群众

原文发布于微信公众号 - 生信技能树(biotrainee)

原文发表时间:2016-12-31

本文参与腾讯云自媒体分享计划,欢迎正在阅读的你也加入,一起分享。

发表于

我来说两句

0 条评论
登录 后参与评论

相关文章

来自专栏Python中文社区

用Python把告警日志发到微信上

專 欄 ❈RaPoSpectre,Python中文社区专栏作者。 网站:www.rapospectre.com❈ 前言 笔者所在公司项目的告警信息会通过钉钉发...

6079
来自专栏运维前线

阿里&百度&腾讯&facebook&Microsoft&Google开源项目汇总

BAT && YMFT Tencent GitHub地址:https://github.com/Tencent/tinker Tinker是Android的...

3518
来自专栏java工会

在阿里技术一线呆三年,你会遇到哪些故障

在技术面试的时候,我们通常会被问到“你遇到过哪些问题,是怎么解决的”。这个问题就很考验经验了,如果你在一个小作坊呆了很多年,你可能根本就不会遇到这些问题。所以面...

632
来自专栏通信云团队

开发一款即时通讯App,从这几步开始

9114
来自专栏云计算D1net

OpenStack加入Apache顶级项目Cassandra

Apache Cassandra是极高性能、可扩展、分布式NoSQL数据库,使用灵活,简单分区行存储数据模型,可以对商业服务器和跨数据中心进行无单点故障的...

2966
来自专栏知晓程序

好消息!小程序可关联 50 个公众号了!

下面,知晓程序(微信号 zxcx0101)带大家看看,微信官方给大家带来了什么「夜间好消息」。?

802
来自专栏FreeBuf

QQ蠕虫的行为检测方法

作者 Nandisec 选题背景 QQ蠕虫是一种利用QQ等腾讯公司相关产品进行传播的一种特殊蠕虫,该蠕虫的基本原理是利用了QQ帐户的快速登录机制,只要当前...

1938
来自专栏*坤的Blog

提高效率

1712
来自专栏java工会

在阿里技术一线呆三年,你会遇到哪些故障

在技术面试的时候,我们通常会被问到“你遇到过哪些问题,是怎么解决的”。这个问题就很考验经验了,如果你在一个小作坊呆了很多年,你可能根本就不会遇到这些问题。所以面...

842
来自专栏移动开发试验田

【移动开发】市面上主流「移动推送服务」的体验比较

推送服务基本上是每个 App 的刚需,自己也用过许多家推送服务,最近腾讯云上线了一个类似于 firebase 的移动开发平台,上面集成了很多的移动服务,包括推送...

3587

扫码关注云+社区