一些需要知道的概念生物信息学生信上游

理解sam 与bam 文件

2020-06-15  本文已影响0人  Peng_001

参考:碱基矿工
生信学苑

通常我们可能会看到.cram.sam.bam等多种序列比对的文件,其实他们的本质都源于sam。
bam 是sam 则是一种比gz更加高效的压缩算法,可以很好的对sam 文件进行压缩;
而cram则是bam的高压缩格式(相比于BAM有着更高的压缩率,能够节省30%-50%的空间),IO效率比原来的BAM要略差(同样慢30%左右),且CRAM和BAM可以通过samtools或者picard方便地实现互转。

SAM

SAM 文件的全称为SAM(The Sequencing Alignment/Map Format的英文简称),由sanger 定制,是以TAB 为分割符的文本格式。主要应用于测序序列mapping到基因组上的结果表示,也可以表示任意的多重比对的结果。

它是bwa比对输出的标准格式,由于其是纯文本内容,因此缺点便是内容巨大,比如一个人30x全基因组测序的sam大小超过600G,非常容易造成存储爆满。

因此开发者想到了压缩sam 文件这一方法,而bam 格式也在层出不穷的各类压缩文件中脱颖而出。

BAM

正如上面所说,bam的本质其实就是sam 文件的一种压缩形式(因此二者实际上是同一文件的不同格式)。而bam 对于比对的结果也有着非常全面的记录。

除此之外,bwa作者还为bam 专门开发了samtools,更便于处理bam文件,也增加了原有的拓展性。

查看文件

通常来说,sam文件可以像其他文本文件一样直接查看(less等),但bam文件(为二进制文件)则需要借助samtools工具。

samtools view -h E_coli_K12.sorted.markdup.bam|less -SN
# h选项会将bam 的header 信息也一并输出

记录的信息

序列比对,至少应该包含以下的内容与信息。


一般来说,bma 文件分为两个部分,headerrecord。(大部分的组学数据也都是按照这两种内容分割来存储数据的)

header 部分

header 部分信息不会太多,一般来说每一行以@开头,里面主要包含了文件的各类参数信息。

@HD 说明符合标准的版本、对比序列的排列顺序;
@SQ,参考序列说明(如果是人则是和参考序列一致的各个染色体的信息);
@RG,比对序列(reads)的说明;
@PG,使用的程序说明(操作过程与相关参数);
@CO,任意的说明信息(可有可无)。

@RG 主要是用于区分不同样本的重要信息,它是测序数据合并的关键。

record 部分

record ,也叫 alignment section,是序列比对的内容,从header 往下record 部分的每一行都是一条read 的比对信息。并且read的每个信息都会用tab 分开,且一共有12列。



其中有几个重要的信息。

flags

flag 主要记录的是序列的比对信息。实际的flag 应为01 组成的二进制码,一共有12位数字,12位数字分别代表12种0与1 对应的组合情况。(数值取值范围也就是0~2048)


比如下图这个数字

2209 = 100010100001,也就表示为PE测序(1),它的另一条配对read 反向互补后比对到参考序列(6),且为read2片段(8),最后,该read 还可能存在嵌合,也即比对的这部分可能只是来自于其中的一部分序列,supplementary alignment(12)。

一般来说,我们可以自己写脚本,无非就是将原本的十进制数字提取出来,再转化为二进制形式,最后再通过比对每一行的信息,通过 TRUE(1) 还是FALSE(0)的判断,即可得到其表示的结果。比如若转化后的二进制形式的第三位数字为1,则表示该read 没有比对到参考序列上。

另外,python也提供了相应的包可以方便我们的bam 信息处理。

MAPQ 质量值

这个质量值和fastq 格式中的质量值的转化与表示方式一样,都是-10logP(P为错比概率)。它是衡量read 比对质量的一个重要指标。

CIGAR

又称“雪茄值”,全称为,Compact Idiosyncratic Gapped Alignment Report,它和flag 类似,也通过几个字符串,但表示了不同的read比对到参考序列的细节信息。CIGAR的标记字符有MIDNSHP=XB这10个。


比如167H134M表示,比对这条read 的开头167bp被硬跳过了,且紧随其后的134bp比对上了参考序列。另外并非M就是完全匹配,其也包括某些单碱基错配。

其他信息

其他的还包括read 质量值,其实也就是囊括了fastq 文件的信息。有的时候fastq 文件也会存储为uBam(un-mapping BAM)文件,这其实也就是没有做过比对的BAM文件。

另外,除了前11列信息之外,metadata 的内容是不固定的。这点和不同的处理软件输出内容不同也有一定关系。

可视化查看

之前提到,bam 与sam 不同,它是一种压缩后的二进制文件,因此需要使用samtools view进行查看。

此外我们还可以提取指定的某段序列的信息。

samtools view -h E_coli_K12.bam NC_000913.3:1000000-1000200

如果想要可视化,可以使用tview指令。这是samtools 自带的终端工具。

samtools tview --reference in.fa in.bam

在该模式,键盘输入g,就可以调整到对应位置。具体信息可以通过输入?查看。

使用IGV可视化比对结果

http://www.igv.org/

将.bam 文件导入后,就可以可视化查看相关的比对结果。

上一篇下一篇

猜你喜欢

热点阅读