生物信息学分析实用小技巧(五):BAM/SAM两三事

转自公众号:universebiologygirl
http://mp.weixin.qq.com/s?__biz=MzU1NDkzOTk2MQ==&mid=2247485110&idx=1&sn=a087717670d8ce219709068f55a2b7d1

BAM/SAM文件两三事

写在前面的

上期推文主要讲了Linux系统下三种语言环境的管理,R、Perl和Python作为生物信息学中常用的编程语言,但是这三门编程语言的环境在系统中的管理有不相同的地方。Linux系统自带了Perl和Python,但是Linux系统只能支持一个Perl的存在,对于其他如qiime类的软件需要用perl最好把这些软件安装在一个独立的虚拟环境中,对于python由于系统支持多个版本的python所以只需要知道如何快速切换python版本就行。对于R语言,由于大部分系统并不自带这个语言,所以只要选择一个合适的版本安装就行。

BAM/SAM文件的细致解读

        这期的推文,主要讲述BAM/SAM文件的解读,咱们公众号已经解读过这些文件了。这次我主要从文件的应用角度来讲,一次肯定是讲不完的,为此我专门以这个文件的内容为主体,分多期推文讲述。主要目的在于让大家再次加深对这些文件的映像,同时具有独立的编程能力完成一些分析。其实,BAM文件的爸爸是SAM文件,SAM文件是Sequencing Alignment/Map Format的英文首字母缩写。由于早期NGS分析中序列比对软件多种多样,它们输出的文件也是五花八门。后来李恒大神开发了bwa,这个软件输出的结果是SAM文件。正因为SAM文件的出现使得后来在NGS分析中序列比对的比对信息储存格式标准化。所以不管你现在是做RNAseq还是WGS亦或是其他组学分析如宏基因组等,只要涉及到对高通量数据进行序列比对,就会出现SAM文件。但是,SAM文件因其占用磁盘空间很大,非常不利于数据的储存。李恒大神由在SAM的基础上开发了BAM储存的方式,储存信息和SAM一样,但是大大节约了储存空间。例如当你使用如下命令查看BAM文件时

samtools view xxx.bam | head -n 100

生物信息学分析实用小技巧(五):BAM/SAM两三事

        这个结果看上去十分复杂,其实我们公众号已经网上都有资料讲过这里是按照如下信息分隔的

生物信息学分析实用小技巧(五):BAM/SAM两三事

                                                     图片来源于网络

        光看上图的信息也是相对较多的。今天我们从flags信息讲起。这个flags是SAM/BAM文件record部分的第二列。flags详细记录了reads的比对信息。仔细观察发现,flags是一系列数字如99,83,147等。这里数字实际上是由一堆2的n次方(n=0~11)加和组成。而这部分数字都经过二进制转化得到一个二进制数字来表示比对的信息。具体信息如下图所示

        例如某条reads经过比对后的flags是83, 83= 1+2+16+64转化为2进制为1010011, 则代表这是一条双端测序的Read1,该reads及其对应reads都是正常比对且当前reads是反向互补后比对上的参考基因组。
        联系到实际的例子:例如我们在宏基因组测序分析中有时候会去掉宿主基因组。那么使用比对软件对宏基因组数据比对后会使用命令samtools -bF4 提取未比对上的reads,而比对上参考基因组的reads则被过滤掉了。samtools则是自动进行了flags值的检测,只要检测到某条reads的比对信息中有4那么就判断出这条reads是未比对上参考基因组的。除了这个情况外,有时候我们在检查比对结果时更希望了解reads比对的更多信息。当然还是以flags值为例,如果我们一个值一个值去计算它等于多少,就比较麻烦。如果我们设计一个程序输入BAM文件的比对信息,输出flags值对应的人类语言信息。在这里我们采用perl语言为基础编写这个程序。要解析reads比对情况,把10进制数值转化为2进制数值更方便。10进制整数转化为2进制,就是除2,得到设为结果逆时向上即为2进制整数,列如27,2进制为11011,用图表示则是


程序写好后进行实际操作

samtools view BAM.bam | head -n 10000 > BAM.info.txt
perl Describe_BAMinfo.pl BAM.info.txt > reads_alignmentdescribe.txt

这样每条reads比对的信息就一目了然了。

写在后面的

这次主要采用理论联系实际,以BAM文件的内容的flags值为中心,通过生物信息学定义加以基本数学知识进制转化来挖掘reads比对信息。利用这些信息,我们都可以自己提取比对上或者未必对上的reads了。这次主要利用perl语言编程。关于perl语言编程我后期专门通过几篇推文来讲述。

参考

碱基旷工-知乎
samtools帮助文档

往期精彩推文

SAM文件相关的进阶知识点(中)

SAM文件相关的进阶知识点(上)

作者信息

熊东彦, 中国科学院武汉病毒研究所在读研究生,擅长方向:转录组分析、宏基因组分析、R语言编程、Perl语言编程。近期推文生物信息学分析小技巧

生物信息学分析实用小技巧(五):BAM/SAM两三事》来自互联网公开内容,收录仅供学习使用,如侵权请联系删除。本文URL:https://www.ezixuan.com/1020480.html

(0)
上一篇 2023年 1月 27日 上午9:48
下一篇 2023年 1月 27日 上午10:03