BAM 文件解读

发表于 分类于 ， 阅读次数： 本文字数： 4.2k 阅读时长 ≈ 4 分钟

Introduction

了解生信的都知道，BAM 是一种测序数据比对信息存储的一种文件格式，类似的还有 SAM 格式；可以这么理解：BAM 是 SAM 经过压缩之后的一种文件格式，其压缩率可以减少很大的同级别的存储空间。

BAM 文件中存放了大量序列比对相关的信息，读懂 BAM 才能更好的理解解析生物大数据。

Getting Started

sam 文件是纯文本文件；bam 文件是 sam 文件的压缩格式；bam 文件的大小只有原来 sam 的 1/6（大约）；

bam文件共有两部分：header+record

BAM-header

1
2
3
4
5
6
7
8
9
10
## 查看 header
$ samtools view -h test.bam | head
## output
@HD	VN:1.6	SO:coordinate
@SQ	SN:chrM	LN:16571
@SQ	SN:chr1	LN:249250621
@SQ	SN:chr2	LN:24319937 3
.......
@RG	ID:sampleid 	PL:Illumina	LB:lib	SM:sampleid   ## bwa -R 给定的参数；如：bwa mem -t 4 -M -k 30 -R '@RG\tID:sampleid\tPL:Illumina\tLB:lib\tSM:sampleid' ....
@PG	ID:bwa	PN:bwa	......
  • 解释：
    • HD：VN-版本号；SO-排列方式
    • SQ：SN-参考序列序号；LN-长度
    • RG(read group)：ID-样本信息；
    • PG：ID-比对工具；以及比对时运行的命令

BAM-record(比对信息)

每一行表示每一条read的比对情况

1
2
## 查看比对信息
$ samtools view test.bam | le

  • 解读：

    每一行的结构如下：

    QNAME - flags - position - MAPQ - CIGAR - mate information - read sequence - quality scores - metadata

    • read name:
      • QNAME(Query template NAME): fastq 文件中的 read id；即对应测序fastq数据中read的名称
    • flags:
      • flag(bitwise FLAG): 二进制的比对信息位；一共有12位，每一位都代表了一个特定信息
    • position:
      • RNAME(Reference sequence NAME): 参考基因组对应的染色体号；若为 * 表示没有比对到参考基因组; RNAME中的R表示Reference
      • POS(1-based leftmost mapping POSition): 比对位置，从对应染色体的第一位开始往后计算
    • MAPQ(MAPping Quality): 比对质量值(Mapping Quality)；值越大，表示比对错误的概率越低；如 MAPQ>30 -- 表示错比概率低于0.001（千分之一）；计算公式：-10logP{错比概率}
    • CIGAR(Concise Idiosyncratic Gapped Alignment Report): 比对信息，用数字和字符的组合形象的记录read比对到参考序列上的细节
    • Mate information:
      • RNEXT(Reference name of the mate/next read): 配对read所比对到的染色体（仅PE测序的数据才有）; ‘=’ -- 表示前后是位于同一个片段；‘*’ -- 表示没有对应的片段信息
      • PNEXT(Position of the mate/next read): 下一个片段比对上的位置；配对read所比对到的位置（仅PE测序的数据才有）；若没有相关信息则为0
      • TLEN(observed Template LENgth): 插入片段长度（仅PE测序的数据才有）；TLEN的绝对值等于模板序列的映射末端与模板序列映射起点之间的距离（end-start+1; 包括两端）；【比对上的碱基不包括soft-clip碱基】，值的正负代表着比对上的正负链；单链，值为0
    • read sequence:
      • SEQ(segment SEQuence): read序列；若为 ‘*’ 表示该序列没有存储，否则会有对应的序列片段，长度 和CIGAR中的 M / I / S / = / X 是对应的。
    • quality scores:
      • QUAL(ASCII of Phred-scaled base QUALity): read 质量值
    • metadata: 元信息，从第12列开始往后都是metadata，一般包括RG信息，missmatch信息，二次比对信息等；在不同的处理软件中输出会有所不同，不固定

BAM-record-flags

FLAG(十进制) Bitwise(二进制序列) 十六进制 内容描述
1 000000000001 0x1 Pair-End(PE)测序，0表示Single-End(SE)测序
2 000000000010 0x2 正常比对（PE测序的两个read都能正确比对到参考序列上）；若是PE测序，表示PE的两条read之间的比对距离没有明显偏离插入片段长度
4 000000000100 0x4 表示该read没有比对到参考基因组上
8 000000001000 0x8 PE测序的另一个配对read没有比对到参考序列（不是当前read）
16 000000010000 0x10 反向互补后比对到参考序列，即比对到负链
32 000000100000 0x20 PE测序的另一个配对read反向互补后比对到参考序列
64 000001000000 0x40 PE测序read1；说明该read是与参考序列比对的第一条
128 000010000000 0x80 PE测序read2
256 000100000000 0x100 二次比对，说明该read在参考基因组上比对了多个位置，当前比对位置是次佳比对位置，通常需要过滤掉；但可能存在生物信息
512 001000000000 0x200 低于过滤阈值
1024 010000000000 0x400 PCR重复序列(来自测序文库构建过程)或者光学重复(来自测序过程)；即 dup
2048 100000000000 0x800 说明该read可能存在嵌合，这个比对的部分只是来自其中的一部分序列

77 = 000001001101 = 1 + 4 + 8 +64，这样就得到了这个FLAG包含的意思：PE read，read比对不上参考序列，它的配对read也同样比不上参考序列，它是read1

当然，如果你希望自己在程序中写一段处理FLAG的代码，那么显然是不会像我们这个例子那样去分解这个整数的，多麻烦啊！那么该如何做呢？其实也很简单，比如我们要获得其中某个位（假设第N位）的值——只需要将这个FLAG值和2的N次方做与的运算即可。在与运算时，FLAG值首先会被转换成一串二进制序列（如77=000001001101），而2的N次方除了第N位是1之外，其它的都是0，“与”了之后其它信息就会被屏蔽掉。比如，我们想知道该read是否比对上了参考序列，那么只需要计算FLAG & 4 的值就行了，如果结果是1那么就是比对上了，如果是0则代表没有比上。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
## samtools 根据 flags 值对reads进行筛选
$ samtools flags 10
## output
0xa	10	PROPER_PAIR,MUNMAP

## 
Flags:
	0x1		PAIRED        .. paired-end (or multiple-segment) sequencing technology
	0x2		PROPER_PAIR   .. each segment properly aligned according to the aligner
	0x4		UNMAP         .. segment unmapped
	0x8		MUNMAP        .. next segment in the template unmapped
	0x10	REVERSE       .. SEQ is reverse complemented
	0x20	MREVERSE      .. SEQ of the next segment in the template is reversed
	0x40	READ1         .. the first segment in the template
	0x80	READ2         .. the last segment in the template
	0x100	SECONDARY     .. secondary alignment
	0x200	QCFAIL        .. not passing quality controls
	0x400	DUP           .. PCR or optical duplicate
	0x800	SUPPLEMENTARY .. supplementary alignment

BAM-record-CIGAR

更加直观的体现出read比对到参考序列的比对情况，如插入，缺失，跳跃等

标记字符 内容描述
M 匹配（包含完全匹配和单碱基错配）
I 序列插入
D 序列缺失
N 跳过参考序列，常见于RNA-Seq数据的比对；N 表示可变剪切位置
S 软跳过(soft-clip)，跳过read中的部分序列，不会改变read的长度
H 硬跳过(hard-clip)，直接剪切掉read中分布序列，会改变read的长度；H 只会出现在reads的两端，不会出现在中间，且一般与 S 成对出现
P Padding；和 N 类似，是read比对时中间跳过参考序列的部分区域
= 完全匹配
X 序列错配
B BAM_CBACK，一个曾经适用于CG测序技术特点的标记（华大BGI SEQ测序系列），作用和 N 类似，跳过一定长度的参考序列
* 显示该条 read 完全错配；即与参考基因组不存在匹配的位置

重点说明的是：M -- 并不能表示reads完全比对上参考基因组序列，不能百分百认为其中没有miss-match；多态性碱基或单碱基错配也是用M标记

  • Hard Clip：存在的本意，是减少BAM文件序列的冗余度，比如有一条read，它能比对到A，B两个地方，在A地方，是60M90S，在B地方是60H90M，此时一条read其实已经在A位置有了完整的序列信息，在B位置的信息其实是冗余的，所以在B位置可以引入Hard Clip这样一个标记形式，就能把B位置的序列标记为secondary。
  • 软连接最终会在后面的序列中保留对应的序列，而硬连接会在后面的序列中直接删掉该片段

References

[1] https://www.jianshu.com/p/11275987b24b

[2] https://www.jianshu.com/p/705330b383c3