什么是BAM

BAM是目前基因数据分析中最通用的比对数据存储格式，它既适合于短read也适合于长read，最长可以支持128Mbp的超大read！除了后缀是.bam之外，有些同学可能还会看到.cram，甚至.sam后缀的文件，其实它们一个是BAM的高压缩格式(.cram)——IO效率比原来的BAM要略差；另一个是BAM的纯文本格式（.sam）。当然格式都是一样的，因此为了描述上的清晰，我下面都统一用BAM。

BAM文件格式

其实一开始它的名字是SAM（The Sequencing Alignment/Map Format的英文简称），第一次出现的时候，它是bwa比对软件的标准输出文件，但原生的SAM是纯文本文件，十分巨大（比如：一个人30x全基因组测序的sam大小超过600G），非常容易导致存储空间爆满！为了解决这个问题，bwa的开发者李恒（lh3）设计了一种比gz更加高效的压缩算法，对其进行压缩，这就是我们所说的BAM，它的文件大小差不多只有原来的1/6。

在2007年，NGS技术刚刚兴起之时，各类短序列比对软件层出不穷，输出格式也是各有特点，各家各有一套，并没有什么真正的标准可言，可以说那是一个谁都说我最好的时期。
但逐渐的，研究者们发现BAM格式对Mapping信息的记录是最全面的，用起来也是最灵活的。bwa的作者还为BAM文件开发了一个非常好用的工具包——Samtools，使得人们对BAM文件的处理变得十分便利，拓展性也变得非常强，后来还有类似于IGV等专门支持BAM的工具也越来越多，因此它就逐渐成为了主流。

现在基本上所有的比对数据都是用BAM格式存储的，俨然已经成为了业内的默认标准。

在2013年，研究者们还专门将Samtools的处理核心剥离出来，并将其打包成为一个专门用于处理高通量数据的API——htslib，除了C语言版本之外还有Java和Python版本，这些在github上都能直接找到。后续许多与NGS数据处理有关的工具基本都会使用这个API进行相关功能的开发，可见其影响力。

ok，背景的介绍就先到此为止了，我们回归主题。下面这个图是我从一份刚刚完成比对的bam文件中截取出来的内容：

由于屏幕所限，无法把全部的内容都包含进来，特别是header信息，贴在这里仅是为了让还没见过BAM文件的同学们能够对它有一个总体的感觉。

如果是SAM文件，同时你也熟悉linux操作的话，直接在linux终端用less打开即可（注意：不要试图在本地使用文本编辑器，如vim等直接打开文件，会撑死机子的），但如果我们要查看的是BAM，那么必须通过Samtools（可以到samtools的网站下载并安装）。

$ less -SN in.sam          # 打开sam文件
$ samtools view -h in.bam  # 打开bam文件

BAM文件分为两个部分：header和record。这里额外说一句，许多NGS组学数据的存储格式都是由header和record两部分组成的。

以上例子，在samtools view中加上-h参数目的是为了同时把它的header输出出来，如果没有这个参数，那么header默认是不显示的。

image

header内容不多，也不会太复杂，每一行都用‘@’ 符号开头，里面主要包含了版本信息，序列比对的参考序列信息，如果是标准工具（bwa，bowtie，picard）生成的BAM，一般还会包含生成该份文件的参数信息（如上图），@PG标签开头的那些这里需要重点提一下的是header中的@RG也就是Read group信息，这是在做后续数据分析时专门用于区分不同样本的重要信息。它的重要性还体现在，如果原来样本的测序深度比较深，一般会按照不同的lane分开比对，最后再合并在一起，那么这个时候你会在这个BAM文件中看到有多个RG，里面记录了不同的lane，甚至测序文库的信息，唯一不变的一定是SM的sample信息，这样合并后才能正确处理。

接下来重点要说的是BAM的核心：record（有时候也叫alignment section，即，比对信息）。这是我们通常所说的序列比对内容，每一行都是一条read比对信息，它的记录看起来是这样的：

我这里借用了网上的一张图片来辅助说明，recoed中的每一个信息都是用制表符tab分开的。

下面我们就来仔细瞧瞧这里的每一个信息分别都是什么。

image

以上，前11列是所有BAM文件中都必须要有的信息，而且从描述中我们也能够比较清楚地知道其所代表的含义。但其中，有几个信息实在太重要了，以至于我认为有必要对其进行详细说明。

flags

对于双端比对的数据，生成的BAM文件中，R1端序列和R2端序列的标识符是一样的，之前一直不知道如何根据bam文件区分哪条序列是R1端，哪条序列是R2端，昨天仔细研究了一下，原来代表R1端和R2端的信息都存储在flag中，即bam文件的第二列；

在bam文件格式中定义了各种flag代表的意思

/*! @abstract the read is paired in sequencing, no matter whether it is mapped in a pair */

define BAM_FPAIRED 1

/*! @abstract the read is mapped in a proper pair */

define BAM_FPROPER_PAIR 2

/*! @abstract the read itself is unmapped; conflictive with >BAM_FPROPER_PAIR */

define BAM_FUNMAP 4

/*! @abstract the mate is unmapped */

define BAM_FMUNMAP 8

/*! @abstract the read is mapped to the reverse strand */

define BAM_FREVERSE 16

/*! @abstract the mate is mapped to the reverse strand */

define BAM_FMREVERSE 32

/*! @abstract this is read1 */

define BAM_FREAD1 64

/*! @abstract this is read2 */

define BAM_FREAD2 128

/*! @abstract not primary alignment */

define BAM_FSECONDARY 256

/*! @abstract QC failure */

define BAM_FQCFAIL 512

/*! @abstract optical or PCR duplicate */

define BAM_FDUP 1024

/*! @abstract supplementary alignment */

define BAM_FSUPPLEMENTARY 2048</pre>

1 ： 代表这个序列采用的是PE双端测序

2： 代表这个序列和参考序列完全匹配，没有插入缺失

4： 代表这个序列没有mapping到参考序列上

8： 代表这个序列的另一端序列没有比对到参考序列上，比如这条序列是R1,它对应的R2端序列没有比对到参考序列上

16：代表这个序列比对到参考序列的负链上

32 ：代表这个序列对应的另一端序列比对到参考序列的负链上

64 ： 代表这个序列是R1端序列， read1;

128 : 代表这个序列是R2端序列，read2；

256： 代表这个序列不是主要的比对，一条序列可能比对到参考序列的多个位置，只有一个是首要的比对位置，其他都是次要的

512： 代表这个序列在QC时失败了，被过滤不掉了（# 这个标签不常用）

1024: 代表这个序列是PCR重复序列（#这个标签不常用）

2048: 代表这个序列是补充的比对（#这个标签具体什么意思，没搞清楚，但是不常用）

上面的这几个标签都是2的n次方，这样的数列有一个特点，就是随机挑选其中的几个，它们的和是唯一的，比如

65 只能是1 和 64 组成，代表这个序列是双端测序，而且是read1

所以在bam文件中的第二列，即flag列的值代表这条序列符合上述所有条件的值的和，所以根据这个flag我们可以确定这条序列究竟是read1 还是read2

在新版的samtools中提供了一个flags 功能，可以查看flag 代表的含义，比如下面的sam文件

NB500986:26:HGJ2VBGXX:3:13403:17250:17986 419 chr1 12013 1 131M = 12016 131 TCTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAG
CACTGGAGTGGAGTTTTCCTGTGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCAT AAAAAEEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEEEEEEAEEEEAEEEEEEEEEEEEEEEEEEEE
EEEEEEEEEEEEEEEEEEEEEEEEE/EE/EE<EEEEEEEEEEEEEEEEEEEEEEEEEEEEEE AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:131 YT:Z:UU XS:A:+ NH:i:4 CC:Z:
chr15 CP:i:102519027 HI:i:0 NB500986:26:HGJ2VBGXX:3:13403:17250:17986 339 chr1 12016 1 128M = 12013 -131 GACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAGCAC
TGGAGTGGAGTTTTCCTGTGGAGAGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCAT EEEEEEEEEEEEEEEEEEEEEEEEEEAEEEEEEE<AEEEEEEEEEEEAEEEAEEEEEEEEEEEAEEEEE
EEEEEEEEEEEEEEEEEEEEEAAEEAEEEEEAEEEEEEEEEE/EEEEEEEEEEEAAAAA AS:i:0 XN:i:0 XM:i:0 XO:i:0 XG:i:0 NM:i:0 MD:Z:128 YT:Z:UU XS:A:+ NH:i:4 CC:Z:
chr15 CP:i:102519027 HI:i:0

一共有两个序列，他们的标识符是一样的，那怎么区分那个是R1端，那个是R2端呢？

直接根据flag 查看

samtools flags 419
0x1a3 419 PAIRED,PROPER_PAIR,MREVERSE,READ2,SECONDARY

samtools flags 339
0x153 339 PAIRED,PROPER_PAIR,REVERSE,READ1,SECONDARY</pre>

根据上面给结果可以看出，flag为419代表PAIREED(1), PROPER_PAIR(2),MREVRSE(32),READ2(128),SECONDARY(256)，其中

PAIRED 代表这条序列采用双端测序， 其值为1；

PROPER_PAIR 代表这条序列完全匹配， 其值为2；

MREVRSE 代表这条序列对应的另一端序列比对到参考序列的负链上，其值为32；

READ2 代表这条序列是R2端序列，其值为128

SECONDARY 代表这条序列不是primary alignment, 其值为256

1 + 2 +32 +128 + 256 = 419

所以419对应的序列是R2端序列， 同理，339代表的序列是R1段序列

搞清楚了序列究竟是R1端还是R2端之后，还有一点值得注意的地方，samtools 还提供了一个bam2fq 的功能，根据bam文件抽取比对时采用的fastq序列

在bam文件中，第10列代表的是序列，第11列代表的是序列的质量，根据第二列的flag还可以确定序列是R1还是R2,

但是当序列本身比对到参考序列的负链时，即flag 包含16时，bam文件中记录的序列是原始序列的反向互补序列，而且质量值也是反向的，所以根据这样的序列还原时要小心一点，以上面flag为339的序列为例，

因为339包含了REVERSE,所以对应的序列应该是第10列序列的反向互补序列，碱基质量值为第11列的反向序列，

CIGAR

CIGAR是CompactIdiosyncraticGappedAlignmentReport的首字母缩写，称为“雪茄”字符串。
作为一个字符串，它用数字和几个字符的组合形象记录了read比对到参考序列上的细节情况，读起来要比FLAG直观友好许多，只是记录的是不同的信息。比如，一条150bp长的read比对到基因组之后，假如看到它的CIGAR字符串为：33S117M，其意思是说在比对的时候这条read开头的33bp在被跳过了（S），紧接其后的117bp则比对上了参考序列（M）。这里的S代表软跳过（Soft clip），M代表匹配（Match）。CIGAR的标记字符有“MIDNSHP=XB”这10个，分别代表read比对时的不同情况：

除了最后‘=XB’非常少见之外，其它的标记符通常都会在实际的BAM文件中碰到。另外，对于M还是再强调一次，CIGAR中的M，不能觉得它代表的是匹配就以为是百分百没有任何miss-match，这是不对的，多态性碱基或者单碱基错配也是用M标记！

，MAPQ，比对质量值

这个值同样非常重要，它告诉我们的是这个read比对到参考序列上这个位置的可靠程度，用错误比对到该位置的概率值（转化为Phred scale）来描述：-10logP{错比概率}。因此MAPQ（mapping quality）值大于30就意味着错比概率低于0.001（千分之一），这个值也是我们衡量read比对质量的一个重要因子。

剩下的几列在上面的格式表中描述的也比较清楚，基本没有过于隐藏的信息，因此我就不打算再一一细说了，如果大家依然有困惑可以到后台留言。

此外，细心的同学可能也已经发现了：fastq的所有信息都被涵盖到了BAM文件中了，包括比对不上的read也在，因此获得了BAM其实也等于获得了所有的read。而且，fastq有时也会被转换成一种uBam文件，指的就是un-mapping BAM——没有做过比对的BAM文件。它相比于Fastq可以用metadata存储更多有用的信息，不过这不是我们这篇文章想说的内容。

最后，还是再说明一次：BAM文件中除了必须的前11列信息之外，不同的BAM文件中后面记录metadata的列是不固定的，在不同的处理软件中输出时也会有所不同，我们也可以依据实际的情况增删不同的metadata信息。

使用samtools view查看BAM文件

BAM文件由于是特殊的二进制格式，因此没办法通过文本的形式直接打开，要用samtools的view功能在终端上进行查看（上文也已经说到这里在进行系统补充），如：

$ samtools view in.bam

如果不想从头开始看，希望快速地跳转到基因组的其它位置上，比如chr22染色体，那么可以先用samtools index生成BAM文件的索引(如果已经有索引文件则不需该步骤)，然后这样操作：

$ samtools index in.bam  # 生成in.bam的索引文件in.bam.bai
$ samtools view in.bam chr22            # 跳转到chr22染色体
$ samtools view in.bam chr22:16050103   # 跳转到chr22:16050103位置
$ samtools view in.bam chr22:16050103-16050103  # 只查看该位置

IGV或者samtools tview查看比对情况

以上，我基本上列举了我们会在终端上如何查看BAM文件的几个最常用操作。但如果你想更直观查看的BAM文件，IGV是目前最好的一个选择，但仅适合于文件还比较小的情况，效果如下：

如果你的BAM文件很大，都超过了你的本地电脑磁盘了，你还是想看该怎么办？你有两个选择：

第一，把你想查看的那部分区域用samtools view提取出来，生成一份小一些的BAM，然后下载下来，在导入到IGV中。

$ samtools view -h in.bam chr22:16050103-16050203 | samtools view -Sb - > small.bam

第二，不下载，直接在终端用samtools tview进行查看。samtools tview有类似于IGV的功能，虽然体验会稍差一些。

$ samtools tview --reference hg38.fa in.bam

在该模式下，按下键盘‘g’后，会跳出一个Goto框，在里面输入想要调整过去的位置，就行了，比如：

image

按下esc键则可以取消。另外，为了节省空间，加快查询效率，read中与参考序列相同的部分被用一串串不同颜色的点表示，只留下miss-match的碱基和发生indel变异的区域。其中圆点表示正链比对，逗号表示负链比对。不同的颜色代表不同的比对质量值：白色>=30，黄色20-29，绿色10-19，蓝色0-9。如果你还想知道的其他的功能，可以在tview模式里按下“?”问号，就会弹出类似下面这样的帮助窗口，然后按照指引做就行了。

虽然看起来不如IGV体验那样好，功能也比较单一（仅可以查看比对情况），但可贵之处在于可以在终端里面直接操作，当需要快速查看某个位置的比对情况时，操作效率非常高。而如果要退出该模式，也非常简单，按下q键就可以了。

---------------------------------------------------------------------------------------------------------------------------------------------------I`m a line ! Thanks !-------------------------------------------------------------------------------------------------------------------------------------------------------------------------

参考：
https://www.cnblogs.com/xudongliang/p/5437850.html
https://luohao-brian.gitbooks.io/gene_sequencing_book/content/di-5-8282-li-jie-bing-cao-zuo-bam-wen-jian.html