基因组比对专题——sam文件格式详解(2)

在最近的学习工作中，我越来越认识到生物信息学分析重要的基础之一是对文件格式与含义有着深刻的认识。因此，在深入了解比对之前，必须对其中所要用到的文件格式加以研究。本文参考自SAMv1.pdf (samtools.github.io)

本文对SAM文件的比对区域进行了说明，前序知识请参考基因组比对专题--sam文件格式详尽解读 (1) - 简书 (jianshu.com)

比对区域

在SAM格式文件中，比对区域的每一行通常都用来表示与特定参考序列的比对情况。每行由11个或是更多的字段通过TAB ("\t") 彼此分隔组成。前11个字段的内容与顺序被严格要求与定义。若是在这些字段（前11个）字段中，某些信息不可用，那么必须使用"0"或"*"进行站位，选择哪种取决于该字段应该包含的数据类型。下面将分别介绍这11个字段。

1.QNAME

这是在比对区域中每行中第一个记录的字段。全称Query template NAME，翻译为查询模版名称。对于读段或是片段，它们都有独特的QNAME。若是读段或片段具备相同QNAME，则认为它们来自同一模板。如果QNAME字段为"*"时，则该信息不可用或是未指定。

当一个读段可以与参考序列上的多个区域都对应上时（多重比对），或是可能跨越参考序列上多个区域从而形成嵌合比对，即一个读段的不同部分最佳比对结果是参考序列上的一段不连续区域（通俗理解：读段的某部分比对到参考序列的A处，而又有一部分比对到B处），那么这种情况下同一QNAME可以在比对文件中的多行中。

2.FLAG

这是在比对区域中每行中第二个记录的字段。全称Bitwise Flag，翻译为位标志。该字段使用二进制表示法，其中每位（可以形象理解为十位，百位，千位...）的布尔值(0或1)表示改行比对结果是否具备该位置属性。在这里使用了12位来表示比对结果是否分别具备12种属性（这种描述仅针对目前广泛使用的samv1.6版本）。当然，为了节省空间，在存储信息时，通常会将二进制转换为十进制存储。例如，如果在SAM文件中FLAG值为50，那么就意味着二进制的 “000000110010” 换算为的十进制，表示这一行的比对信息具备2，16，32三种值对应的属性。那么现在将对这十二种属性一一加以说明（这里使用十进制）

• 1 模版含有多个读段。即该读段或序列来自于包含了多个读段或是序列的模版。例如在双端测序中，同一DNA片段上生成的读段被认为是来自一个同一模版。

• 2 来自相同模版的读段或是序列都正确比对到了参考序列上（依据软件对正确比对到参考序列的标准而定）。对于双端测序，通常这意味着双端都被比对到相同的序列上，且其相对方向复合预期，且两个读段间的距离符合测序平台与实验设计预期。

• 4 表示该读段或序列未比对成功。即该读段或是序列没有比对到参考序列上的任何区域。

• 8 表示模版中的另一个没有被比对到到参考序列上。这个概念一般用于双端测序中，即一端的读段比对到了参考序列上，但是另一端没有。

• 16 表示SEQ序列被反向互补。我们都知道，DNA序列是由两条互补的链组成，分为正向链和反向链。当序列与基因组比对时，它可以被比对到正向链，也可以是反向链。那么如果SEQ中这个读段记录的序列匹配到基因组的反向链，那么这个序列在方向会被反向，并进行碱基互补替换，从而使其与参考序列的正向链匹配。如果发生这种情况，那么就用该值记录在FLAG中。

• 32 表示SEQ对应模版的下一段读段或序列被反向互补，这可以结合FlAG中16和下面的RNEXT和PNEXT加以理解。即如果模版中的下一段读段经历了形如FLAG=16中的处理，那么就会添加值32

• 64 表示该行比对的读段为模版中的第一个片段

• 128 表示该行比对的读段为模版中的最后的片段

• 256 表示改行比对为次级比对。在比对结果中，一个读段可能会在参考序列上存在多个比对成功区域，这种情况下，比对工具往往会选择一种比对结果作为主要比对（primary alignment），而其余结果作为次级比对（secondary alignment）。这里的次级只是说这些结果不是作为该读段比对的代表性结果，而并非说比对结果很烂。

• 512 表示该比对结果中的读段为通过特定的质量控制而被认为是低质量的。

• 1024 表示该比对结果中的读段可能是由于PCR扩增或是光学错误引起的重复，通常要在分析中去除含有该标志的片段以减少偏差

• 2048 补充比对。表示该比对结果是嵌合比对结果中的一部分，嵌合比对在基因组比对专题--sam文件格式详尽解读 (1) - 简书 (jianshu.com)有所提及。

补充说明：

对于SAM文件，如果存在多个SEQ相同的行，必须只能有一行满足FLAG的值不可分解为256或是2048与其他数字的和。可以看到，256和2048分别表示次级比对和补充比对，如果FLAG值不能分解出它们，那么就代表了该种比对结果既不是次要比对也不是嵌合比对，从而确定该比对结果为主要比对。而每个SEQ只能有一个主要比对。

3.RNAME

这是比对区域当中记录的第3个字段，英文全称为reference name。该字段用于表示QNAME中的读段或是序列比对到参考序列的名称。并且如果该SAM文件中头部区域使用了@SQ（相关知识参考基因组比对专题--sam文件格式详尽解读 (1) - 简书 (jianshu.com)），那么RNAME（除非是“”）必须出现在头部区域中的SN的TAG中，以形成一一对应的关系。当然，可以赋予一个比对不成功的读段或序列一个普通的坐标，例如这可以使其参与到其后的对文件进行排序的过程中。但是，如果RNAME为“”，那么POS字段和CIGAR字段不应当包含任何有效信息。

4.POS

这是比对区域当中记录的第四个字段，英文全称为position，表示位置。该字段表示读段或是序列映射到参考序列的起始位置。这个位置坐标是基于步长为1计数的，并且是以第一个“消耗”参考序列碱基的CIGAR操作的最左侧的比对位置为值。而参考序列的起始坐标设置为1。举例说明：如果GIGAR操作消耗的第一个碱基是参考序列上的坐标为10的碱基，那么POS就记为10。

其中，CIGAR字符串中的一些操作被认为是“消耗”参考序列碱基的，这意味着这些操作涉及到了参考序列上的碱基。例如，匹配（M）、删除（D）、错配（X）和等同（=）操作都是“消耗”参考碱基的操作，因为它们都涉及到了参考序列上的位置。而插入（I）和软裁剪（S）等操作则不“消耗”参考序列碱基，因为它们仅涉及到读段或序列。有关CIGAR字段详细信息参见下文。

对于没有成功比对到参考序列上的读段或序列，那么POS的值将被设置为0。如果POS值为0，那么RNAME和CIGAR区域也无法描述任何有效信息。

5.MAPQ

这是比对区域中记录的第五个字段，英文全称为mapping quality，翻译为比对质量。它的计算方法为：

其中value表示比对质量分数，rate表示比对错误概率。当比对质量不可用时，该字段值会被设置为255

6.CIGAR

这是比对区域当中记录的第六个字段，英文全称为Concise Idiosyncratic Gapped Alignment Report，翻译为简洁的特殊间隙比对报告。说人话就是，通过字符表示查询序列与参考序列的比对结果的具体内容是什么，为了能够清楚知道接下来对字符的解释描述里的内容，我先举一个例子：假设某个CIGAR字符串为10M5I10M，那么这就表示查询的前10个碱基与参考序列相匹配，随后在查询序列中存在一个5碱基的插入片段，紧接着查询序列的10个碱基由于参考序列匹配或不匹配，具体的字符包括哪些呢？如下：

操作符	描述	消耗查询序列	消耗参考序列
M	表示查询序列匹配或不匹配参考序列	是	是
I	向参考序列中插入	是	否
D	从参考序列中删除	否	是
N	从参考序列中跳过	否	是
S	软剪切（在SEQ中呈现的剪切序列）	是	否
H	硬剪切（在SEQ中未呈现剪切的序列）	否	否
P	填充（从填充的参考序列中静默删除）	否	否
=	表示查询序列匹配参考序列	是	是
X	表示查询序列不匹配参考序列	是	是

下面对一些模棱两可的描述给出详细说明

• N：表示从参考序列中跳过的区域，这个操作不消耗查询序列，但消耗参考序列。这种操作通常在mRNA比对中发挥作用，用来表示mRNA比对到参考基因组时跨越的内含子区域

• S：表示在查询序列（读取）的头部或尾部有一部分并未用于与参考序列的对齐。尽管这部分序列出现在读取的SEQ字段中，但在比对评分或分析时被忽略。

• H：硬剪切类似于软剪切，但硬剪切的序列不出现在读取的SEQ字段中。硬剪切通常用于表示已经从读取数据中物理去除的序列。对于S与H来说，H操作只能作为第一个和/或是最后一个出现，而S操作位于CIGAR字符串两端时，只能有硬剪切位于其外侧，例如可以有H10S20M30S10H，但是不能有S10H20M30H10S

• P：这是一种虚拟操作，用于占位而不对应于查询序列或参考序列中的任何实际序列，它通常用于保持多个片段比对时的一致性特别势必对比对到同一位置的不同读段时。例如在处理复杂的结构变异时，需要在对齐的读取之间插入空间以保持它们的相对位置不变。

• M/I/S/=X操作的长度之和必须等于SEQ的长度。

7.RNEXT

这是比对区域当中记录的第七个字段，英文全称为Reference sequence name of the primary alignment of the NEXT read in the template，翻译为模版中下一个读段主要比对的参考序列的名称。该字段用于描述来自同一模版的某一个读段的下一个读段主要比对到的参考序列名称。

如果该读段已经是模版中最后一个读段，那么下一个读段就变成了第一个读段。与RNAME一样，如果使用了@SQ，那么RNEXT字段中的值也要与某个@SQ行中的SN的TAG值相同。当信息不可用时，此字段设为"*"。

如果RNEXT与RNAME相同，则设置为“=”。如果模版中的下一个读段有一个主要比对（即该这种比对结果被认为是最佳比对结果），并且RNEXT不是“=”，那么RNEXT应当与下一个读段或序列的RNAME相同。

此外，如果RNEXT是“*”，则无法给出PNEXT和一些相关FLAG值，如32。

8.PNEXT

这是比对区域当中记录的第八个字段，英文全称为Position of the primary alignment of the NEXT read in the template，翻译为，模版中下一个读段主要比对的参考序列的位置。类似于POS之于RNAME，PNEXT与RNEXT也是这种关系，并且计数方法也想同，不难发现，其值就等同于下一个读段比对结果所在行的POS。如果下一个读段或序列的比对位置信息不可用，那么PNEXT设置为0，此时，不能对RNEXT和一些相关FALG值进行使用，如32。

9.TLEN

这是比对区域当中记录的第九个字段，英文全称为signed observed template length，翻译为带符号的观察到的模版长度。该字段的目的是提供一个快速地方式来了解读段在参考序列上的位置和跨度，特别是对双末端测序结果来说。这里的长度包含了从一端的开始到另一端的结尾，以及他们之间所有序列的长度（|结束-开始 + 1|绝对值）。其正负号代表读段的方向，如果值为正，那么就意味着该读段在双端读取上的位于左侧，为负则是右侧。而程序在计算时，是基于CIGAR字符串的结果，而“S”标注的区域不参与运算。当模版仅为单段或是信息不可用时，这个值将会被设置为0。但是，对于模版比对开始和结束的定义存在争议，故该值的实现以使用的工具为准

10.SEQ

这是比对区域当中记录的第十个字段，英文全称为segment sequence， 翻译为片段序列。它描述了查询序列的本身的序列信息(A,G,C,T...)。当序列信息在这里没有被存储时，通常会使用“*”来代替。否则，序列的长度必须等于CIGAR字段中M/I/S/=/X操作的长度之和。“=”表示基础与参考基础相同。

11.QUAL

这是比对区域中记录的第十一个字段，英文全称为quality，翻译为质量。简单来说，它记录了SEQ字段中每个碱基的质量分数加33后变换的ASCII字符。质量分数的计算方法为:

其中value为质量分数，rate为碱基错误率。将得到的值四舍五入并+33后，对应的ASCII码就是QUAL字段中记载的内容，这种原理与FASTQ文件是一样的。

这里需要注意，如果一个读段或是序列信息不可用，那么QUAL的字段将被设置为"*"。但是如果QUAL字段不为空，那么SEQ序列必须不为空，且QUAL字段长度也必须与SEQ长度相同。

12.可选字段

在结束了上述11个规定字段，随后的字段为可选字段，可以添加更多的比对结果信息，其组织形式一般为TAG:TYPE:VALUE，并与其他格式一样，保持TAB键分隔字段。

一个经典的例子是MD:Z标签，提供了足够的信息来精确地描述哪些位置匹配，哪些位置存在替换错误（即不匹配）。MD标签以一种紧凑的格式描述了参考序列和读取序列之间的匹配情况，包括完美匹配的长度和特定位置的碱基不匹配。

例如，MD:Z:35A64表示前35个碱基完全匹配，第36位上参考序列的碱基是A，但在读取序列中不是（不匹配），后面紧接着是64个碱基的完全匹配。