理解sam 与bam 文件
通常我们可能会看到.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 文件分为两个部分,header
与record
。(大部分的组学数据也都是按照这两种内容分割来存储数据的)
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可视化比对结果
将.bam 文件导入后,就可以可视化查看相关的比对结果。