查看原文
其他

lakeseafly 2018-06-06

当测序后生成的fastq文件比对到到基因序列后,通常会生成一个sam或者bam为扩展名的文件。由于生信分析中,大部分分析都是基于sam格式进行的分析,所以今天我就带大家深入学习一下SAM格式的文件。

认知SAM格式

什么是SAM文件?

SAM的全称是sequence alignment/map format。而BAM就是SAM的二进制文件(B源自binary)。

SAM 格式主要包括两大部分

  1. 标头注释部分(header section)

  2. 比对结果部分(alignment section)

为什么我们需要SAM格式?

SAM格式是用来来支持高通量测序数据分析:

  • 快速查找与坐标重叠的比对。例如,选择与染色体2上的坐标323,567,334重叠的比对。

  • 根据read的属性进行选择和过滤。例如,我们希望能够快速选择能过比对到反向链上的read。

  • 有效地存储数据。例如,从SAM格式转化成BAM格式,单个压缩文件包含所有样本的数据,每个样本都以某种方式标记。

SAM格式的详解

标头注释部分

首先我们看一个例子:

@HD    VN:1.0    SO:unsorted
@SQ    SN:gi|10141003|gb|AF086833.2
|    LN:18959

@PG    ID:bowtie2    PN:bowtie2    
VN:2.2.5    CL:"bowtie2-2.2.5
/bowtie2-align-s --wrapper basic-0
-x 1976.fa -1
SRR1972739_1.fastq -2 SRR19727

39_2.fastq"

标头信息可有可无,都是以@开头,用不同的tag表示不同的信息,主要有:

  • @HD,说明符合标准的版本、对比序列的排列顺序(这里unsorted)

  • @SQ,参考序列说明 (SN:gi|10141003|gb|AF086833.2|)

  • @PG,使用的比对程序说明(这里是bowtie2)

  • LN 是参考序列的长度

比对结果部分

比对结果部分,每一行表示一个read的比对信息,包括11个必须的字段(mandatory fields)和一个可选的字段,字段之间用tag分割。必须的字段有11个,顺序固定,不可用时,根据字段定义,可以为’0‘或者’*‘,这11个字段包括:

第一列: Query Name (QNAME)

这一列代表着比对片段的(template)的编号,例如:

SRR1972739.1

第二列:FLAG

这是一种常用且高效的保存多个布尔特征值的方法。

举个简单的例子: 在 SAM 格式中,当 flag 为 1,也即对应的二进制为 01 时,表示该 read 有多个测序数据 , 一般理解为有双端测序数据 (另一条没被过滤掉), 而 flag 为 2, 也即二进制 10 时, 表示这条 read 的多个片断都有比对结果, 通常理解为双端 reads 都比对上了, 那么就可以推断出 flag 为 3 时, 也即二进制的 11, 表示该 read 有另一端的 read 并且比对成功, 可以看到, 其实就是 01 加 10。

Flag标识对应的情况说明:

简单解析一下常见的几个标识:

要是所得出的标记不是上述列举的数字,比如83=64+16+2+1,就是这几种情况值和。这个网站可以直接通过输入你所得的标记数字,直接告诉其对应的信息:https://broadinstitute.github.io/picard/explain-flags.html

第三列: Reference Name (RNAME)

reference sequence name,实际上就是比对到参考序列上的染色体号。若是无法比对,则是*

第四列: Position (POS)

比对上的位置,注意是从1开始计数,没有比对上,此处为0;

第五列:Mapping Quality (MAPQ)

比对的质量;比对的质量分数,越高说明该read比对到参考基因组上的位置越准确

第六列:Compact Idiosyncratic Gapped Alignment Representation (CIGAR)

CIGAR 代表着简要比对信息表达式,其以参考序列为基础,使用数字加字母表示比对结果。 例如 3S6M1P1I4M

前三个碱基被剪切去除了,然后6个比对上了,然后打开了一 个缺口,有一个碱基插入,最后是4个比对上了。

具体字母意义如下:

“M”表示 match或 mismatch;
“I”表示 insert
“D”表示 deletion
“N”表示 skipped(跳过这段区域)
“S”表示 soft clipping(被剪切的序列存在于序列中)
“H”表示 hard clipping(被剪切的序列不存在于序列中)
“P”表示 padding
“=”表示 match
“X”表示 mismatch(错配,位置是一一对应的)

第七列:MRNM(chr)

下一个片段比对上的参考序列的编号,没有另外的片段,这里是’*‘,同一个片段,用’=‘;

第八列:mate position

下一个片段比对上的位置,如果不可用,此处为0

第九列:ISIZE

Template的长度,mate position大于POS为正,否则为负,中间的不用定义正负,不分区段(single-segment)的比对上,或者不可用时,此处为0;Inferred fragment size.详见Illumina中paired end sequencing 和 mate pair sequencing,是负数,推测应该是两条read之间的间隔(待查证),若无mate则为0;

ISIZE为正,说明amplicon的start 为POS,amplicon 的end为 pos + isize-1

ISIZE为负,说明amplicon的start为mate position,amplicon的end为mate position -isize-1

ISIZE为0,说明无法计算出amlicon的大小

第十列:Sequence

序列片段的序列信息,如果不存储此类信息,此处为’*‘,注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度;就是read的碱基序列,如果是比对到互补链上则是reverse completed

第十一列:ASCII

read质量的ASCII编码。

第十二列:Optional fields

可选的区域格式如:TAG:TYPE:VALUE,其中TAG有两个大写字母组成,每个TAG代表一类信息,每一行一个TAG只能出现一次,TYPE表示TAG对应值的类型,可以是字符串、整数、字节、数组等。

Reference

SAM 格式最权威的解析于介绍:SAM-FORMAT

http://qinqianshan.com/sam-format/

还有更多文章,请移步公众号阅读

如果你生信基本技能已经入门,需要提高自己,请关注上面的生信技能树,看我们是如何完善生信技能,成为一个生信全栈工程师。

如果你是初学者,请关注下面的生信菜鸟团,了解生信基础名词,概念,扎实的打好基础,争取早日入门。

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存