部分转自生信菜鸟团公众号
1.什么是SAM文件?
SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件(B源自binary)。
SAM 格式主要包括两大部分:
1.标头注释部分(header section)
2.比对结果部分(alignment section)
2.为什么我们需要SAM格式?
SAM格式是用来来支持高通量测序数据分析:
(1):快速查找与坐标重叠的比对。例如,选择与染色体2上的坐标323,567,334重叠的比对。
(2):根据read的属性进行选择和过滤。例如,我们希望能够快速选择能过比对到反向链上的read。
(3):有效地存储数据。例如,从SAM格式转化成BAM格式,单个压缩文件包含所有样本的数据,每个样本都以某种方式标记。
3.SAM格式的详解
标头注释部分
标头信息可有可无,都是以@开头,用不同的tag表示不同的信息
比对结果部分
每一列表示一个read的比对信息,包括11个必须的字段和一个可选字段,字段之间用tag分割。必须的字段有11个,顺序固定。这11个字段包括:
第一列: Query Name (QNAME):
这一列代表着比对片段的(template)的编号
第二列:FLAG:
这是一种常用且高效的保存多个布尔特征值的方法。
举个简单的例子: 在 SAM 格式中,当 flag 为 1,也即对应的二进制为 01 时,表示该 read 有多个测序数据 , 一般理解为有双端测序数据 (另一条没被过滤掉), 而 flag 为 2, 也即二进制 10 时, 表示这条 read 的多个片断都有比对结果, 通常理解为双端 reads 都比对上了, 那么就可以推断出 flag 为 3 时, 也即二进制的 11, 表示该 read 有另一端的 read 并且比对成功, 可以看到, 其实就是 01 加 10。
一般flag值不需要自己去算,直接将flag值导入网站即可
http://broadinstitute.github.io/picard/explain-flags.html
所有flag对应值的含义
1 : 代表这个序列采用的是PE双端测序
2: 代表这个序列和参考序列完全匹配,没有插入缺失
4: 代表这个序列没有mapping到参考序列上
8: 代表这个序列的另一端序列没有比对到参考序列上,比如这条序列是R1,它对应的R2端序列没有比对到参考序列上
16:代表这个序列比对到参考序列的负链上
32 :代表这个序列对应的另一端序列比对到参考序列的负链上
64 : 代表这个序列是R1端序列, read1;
128 : 代表这个序列是R2端序列,read2;
256: 代表这个序列不是主要的比对,一条序列可能比对到参考序列的多个位置,只有一个是首要的比对位置,其他都是次要的
512: 代表这个序列在QC时失败了,被过滤不掉了(# 这个标签不常用)
1024: 代表这个序列是PCR重复序列(#这个标签不常用)
2048: 代表这个序列是补充的比对(#这个标签具体什么意思,没搞清楚,但是不常用)
第三列: Reference Name (RNAME):
reference sequence name,实际上就是比对到参考序列上的染色体号。若是无法比对,则是*
第四列: Position (POS):
比对上的位置,注意是从1开始计数,没有比对上,此处为0
第五列:Mapping Quality (MAPQ):
比对的质量;比对的质量分数,越高说明该read比对到参考基因组上的位置越准确
第六列:Compact Idiosyncratic Gapped Alignment Representation (CIGAR)
CIGAR 代表着简要比对信息表达式,其以参考序列为基础,使用数字加字母表示比对结果
例如 3S6M1P1I4M
前三个碱基被剪切去除了,然后6个比对上了,然后打开了一 个缺口,有一个碱基插入,最后是4个比对上了。
这里的总长度对应的就是测出来的一条序列的长度,如果是PE100,这里加起来就应该是100,如果是PE150,这里加起来就是150,这里的长度和第10列的长度是一致的
第七列:RNEXT:
双端测序中下一个reads比对的参考系列的名称。“*”是完全没有比对上,“=”代表完全比对
第3和第7列,可以用来判断某条reads是否比对成功到了基因组的染色体,左右两条reads是否比对到同一条染色体
第八列:PNEXT:
如果是双端测序,是指另一端匹配到参考基因组的位置,如果设置为0,那么该列不可用
第九列:TLEN Template的长度
最左边得为正,最右边的为负,中间的不用定义正负,不分区段(single-segment)的比对上,或者不可用时,此处为0
区别于第6列和第10列是对应测出来的序列的长度。这里第9列的长度是对应插入片段的长度,insert size,也就是建库时,将DNA片段打断成的长度。
第十列:Sequence:
序列片段的序列信息,如果不存储此类信息,此处为’*‘,注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度;就是read的碱基序列,如果是比对到互补链上则是reverse completed。
就是测序的reads序列信息
第十一列:ASCII:
read质量值
其实很容易发现,如果将第1,10,11列提取出来的话,就能还原成我们常见的fastq格式信息。
第十二列:Optional fields:
可选的区域
格式如:TAG:TYPE:VALUE,其中TAG有两个大写字母组成,每个TAG代表一类信息,每一行一个TAG只能出现一次,TYPE表示TAG对应值的类型,可以是字符串、整数、字节、数组等。
备注:
看一下KPGP-00001这个韩国人的测序reads比对到hg38的其中一个lane的sam格式部分信息:
B80KJTABXX:4:1:1404:2065#CTAGTTAT 163 chr7 50962731 60 90M = 50963137 496 AGAAAATTATTTAAATGACCCGAGCCTCACCTTCAACATGAGGAACATCAT
ATTCCCTTTGATAAAATGTTGCTGGTGTAAGTGCTCCAT gggfgfggeggdgggadegggdegegeeggeegedggegggeggegedgggedgggfggeceeggaedgcgggggecgaQ_`X``BBBBB NM:i:0 MD:Z:90 MC:Z:90M
AS:i:90 XS:i:20 RG:Z:KPGP-00001_L6
可以看出这个是用的PE90测序,测序read长度为90bp,建库打断成约490bp,这个read名称是B80KJTABXX:4:1:1404:2065#CTAGTTAT,flag值是163,代表着
image.pngreads是比对到7号染色体,比对的位置是50962731,比对的质量值是60,"90M"意味着90个碱基都match(当然可能是mismatch),“=”意味着双端测序的另一条read也比对上,并且是比对到同一个片段,另一条read比对的位置是
50963137 ,这条read的序列信息是“AGAAAATTATTTAAATGACCCGAGCCTCACCTTCAACATGAGGAACATCAT
ATTCCCTTTGATAAAATGTTGCTGGTGTAAGTGCTCCAT ”
对应质量值分ASCII值为“gggfgfggeggdgggadegggdegegeeggeegedggegggeggegedgggedgggfggeceeggaedgcgggggecgaQ_`X``BBBBB ”
以上。
网友评论