作者:麦茬道
审稿:童蒙
编辑:amethyst
随着生物信息数据的爆发式增长,存储生物信息的文件格式也多样化起来,不同的文件格式往往有不同的目的:
为了软件之间的兼容性;
人类可读的用于数据加工、解析和处理的格式,比如.tsv,.csv等;
为了提高计算机计算效率的数据格式,一般可读性不强的二进制文件,比如今天要讲解到的bam文件。
介绍
bam文件是sam文件的二进制格式,sam 文件是Sequence Alignment/Map Format的简写,产生于比对之后的数据输出,记录了比对的具体情况。文件中以tab键分割,包括上下两部分:
- Header section
- Alignments section
01.Header section
- 该部分全部以“@”开头,提供基本的软件版本,参考序列信息,排序信息等
- @HD行:这一行中有各种不同的标识
- 标识“VN”用以说明格式版本
- 标识“SO”用以说明比对排序的情况,有unknown (default)、unsorted、queryname和coordinate,对于coordinate,排序的主键是Alignments section的第三列“RNAME”,其顺序由@SQ行的“SN”标识的顺序定义,次要排序键是Alignments section的第四列“POS”字段。对于RNAME和POS相等的比对,排列顺序则是任意的
- @SQ行的“SN”标签是参考序列说明,它的值主要是用于Alignments section的第三列“RNAME”和第七列“MRNM”比对的记录
- @PG行是使用的程序说明;该行“ID”为程序记录标识符,“PN”为程序名字,“CL”为命令行
- @CO行是任意的说明信息
02.Alignments section
该部分包含了11列必需字段,无效或者没有的字段一般用“0”或者“*”表示。
下面,以数据详细介绍下bam文件的Alignments section中的每一列的含义。
下图是reads比对情况的示意图
- Coor是参考基因组上的位置坐标
- ref为参考序列
- +r001/1和-r001/2组成1对Reads
- +r003和-r003是反向互补关系,其实是一条read
- +r002、+r003和+r004都是单条read比对
下图是将比对情况记录成bam文件格式的样子
图中的Alignments section有6行12列信息详细介绍了6条read的比对情况,其中前11列为必需字段,每列的含义简单汇总成下表。接下来,我们仔细解读下表1中的信息。
第1列:Qname
Read的名字
第2列:FLAG
每一个read的比对情况可以用十进制数字(或者十六进制数字)表示,如果比对情况 有多个,将多个比对情况所代表的十进制数字加和就是这一行的FLAG。
比如,图1中r001的FLAG是99(1+2+32+64),则表示了“该read是pair read中的一个”,“pair read中每个都能够正确比对上”,“该read的mate read的反向互补可以比对上”,“该read是pari read中的read1”;r001的另一个FLAG是147(1+2+16+128),则表示“该read是pair read中的一个”,“pair read中每个都能够正确比对上”,“该read是原read的反向互补”,“该read是pari read中的read2”(也就是说,该read是read2的反向互补序列)。
值得注意的是,r001是pair read,而且都能比对上,所以r001出现了两次,如果r001的read1比对到参考序列的2个地方,r001的名字则会出现三次;如果read1比对上一次,read2没有比对上,r001仍会出现2次,不过,其中一个r001的第三列为“*”;所以pair-end测序,read1文件和read2文件同时mapping,相同reads的id最少出现2次。
另,以下网站可以通过输入FLAG值,直接找出该FLAG是那些FLAG的加和:Decoding SAM flags(https://broadinstitute.github.io/picard/explain-flags.html)
FLAG释义表如下:
第3列:RNAME
比对上的参考序列的名字,该名字出现在Header section的@SQ行的SN标识中,如果该read没有比对上,也就是说该read在参考序列上没有坐标,那么这一列则用“”表示,那么这一行的POS和CIGAR列也会是“”。
第4列:POS
read比对到的参考序列“RNAME”最左侧的位置坐标,也是CIGAR中第一个比对标识“M”对应的最左侧碱基在参考序列的位置,未比对上的read在参考序列中没有坐标,此列标识为“0”。
第5列:MAPQ
比对的质量值,计算方法为比对错误率的-10*log10的值,一般是四舍五入的整数值,如果是255,说明该比对值无效。
第6列:CIGAR
CIGAR标识符表示read中每个碱基的比对情况,主要有以下标识符:
- M: alignment match (can be a sequence match or mismatch)
read上的碱基与参考序列“RNAME”完全匹配,碱基一一对应,包括了正确匹配与错误匹配 - I: insertion to the reference
read上的碱基相对于参考序列“RNAME”有插入现象(如下):
REF: CACGATCA**GACCGATACGTCCGA
READ1: CGATCAGAGACCGATA
CIGAR:6M2I8M
- D: deletion from the reference
read上的碱基相对于参考序列“RNAME”有删除现象(如下):
REF: AGCTAGCATCGTGTCGCCCGTCTAGCATACGC
READ: GTGTAACC-GACTAG
CIGAR:8M1D6M
- N: skipped region from the reference
read上的碱基相对于参考序列“RNAME”存在连续没有比对上的空缺,这些空缺用N来表示,跟“D”相似但远比“D”缺失的更多,这种比对类型也叫“Spliced alignment”,常见cDNA与参考序列比对(如下):"…"表示intron
REF:AGCATCGTGTCGCCCGTCTAGCATACGCATGATCGACTGTCAGCTAGTCAGACTA
READ: GTGTAACCC................................TCAGAATA
CIGAR:9M32N8M
- S: soft clipping (clipped sequences present in SEQ)
- H: hard clipping (clipped sequences NOT present in SEQ)
read的开头或者结尾部分没有比对到参考序列"RNAME”上,但这部分未比对上的连续序列仍保留在sam文件的该read序列中,用“S”来表示;如果未保留,则用“H”表示,也即“hard cliping”(如下所示,也可同图2中r003的比对CIGAR中看出)
REF: AGCTAGCATCGTGTCGCCCGTCTAGCATACGCAT
READ: gggGTGTAACC-GACTAGgggg
CIGAR:3S8M1D6M4S (在sam中存储为GGGGTGTAACCGACTAGGGGG)
CIGAR:3H8M1D6M4H (在sam中存储为GTGTAACCGACTAG)
- P: padding (silent deletion from padded reference)
多条read比对到参考序列的同一位置时,如果不同read单独同该参考序列比对时,参考序列的情况也不同,比如下方READ1同参考序列比对时,“GA”属于插入(6M2I8M),READ2同参考序列比对时,“A”属于插入(4M1I9M ),READ3同参考序列完全匹配(10M),没有插入,但是三条read之前却没有可比性。因此,当参考序列“比对情况包含完整”且序列唯一时,所有read同时进行比对,read3这种原本没有插入却默认插入的比对称之为Padded alignment,这种情况用“P”表示。
REF: CACGATCA**GACCGATACGTCCGA
READ1: CGATCAGAGACCGATA
READ2: ATCA*AGACCGATAC
READ3: GATCA**GACCG
The padded CIGAR are different:
READ1: 6M2I8M
READ2: 4M1P1I9M
READ3: 5M2P5M
- =:sequence match 正确匹配
- X:sequence mismatch 错误匹配
第7列:MRNM
该read的mate read比对上的参考序列的名字,该名字出现在Header section的@SQ行的SN标识中,
- 如果和该read所在行的第三列“RNAME”一样,则用“=”表示,说明这对read比对到了同一条参考序列上;
- 如果mate read没有比对上,第七列则用“*”表示;
- 如果这对read没有比对到同一条参考序列,那么这一列则是mate read所在行第三列的“RNAME”。
第8列:MPOS
该read的mate read比对到的参考序列“RNAME”最左侧的位置坐标,也是mate read CIGAR中第一个比对标识“M”对应的最左侧碱基在参考序列的位置,未比对上的read在参考序列中没有坐标,此列标识为“0”。
第9列:ISIZE
表示pair read完全匹配到同一条参考序列时,两个read之间的长度,可简单理解为测序文库的长度。这个定义有两种情况(虚线表示未比对上的序列,即soft-clipped bases):
- 第一种情况如图4左所示,Read1和Read2比对到同一条参考序列,此时ISIZE即为Read2的最右侧比对坐标减去Read1最左侧比对坐标;
-
第二种情况如图4右所示,由于至今没有明确的定义和共识,因此ISIZE可以是TLEN#1,也可以是TLEN#2,视情况而定。
第10列:SEQ
存储的序列,没有存储,此列则用“*”标识。该序列的长度一定等于CIGAR标识中“M”,“I”,“S”,“=”,“X”标识的碱基长度之和。
第11列:QUAL
序列的每个碱基对应一个碱基质量字符,每个碱基质量字符对应的ASCII码值减去33(Sanger Phred-33 质量值体系),即为该碱基的测序质量得分(Phred Quality Score)。不同Phred Quality Score代表不同的碱基测序错误率,如Phred Quality Score值为20和30分别表示碱基测序错误率为1%和0.1%。
参考资料
- Sam格式详解手册:http://samtools.github.io/hts-specs/SAMv1.pdf
- Sam格式相关文献:https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2723002/
- Sam格式的wiki介绍(着重介绍了比对类型):https://davetang.org/wiki/tiki-index.php?page=SAM
- Sam格式相关的博客:https://blog.csdn.net/genome_denovo/article/details/78712972
- 快速查询FLAG的网站:https://broadinstitute.github.io/picard/explain-flags.html
- ASCII码对应表:https://snipcademy.com/binf/img/lessons/dna-file-formats/ascii.svg
- 关于测序质量体系的博客:https://blog.csdn.net/xxxxy314/article/details/48948891/
网友评论