理解SAM文件格式以及过滤sam文件

RNA-seq或者ChIP-seq等等测序的上游分析流程里的比对步骤相信大家都知道，我之前也只是按照各种教程去走流程，并没有仔细的研究过每一步的内容。今天这篇文章学习一下sam文件的格式，以及如何根据read比对的质量来过滤你的sam文件。

一般比对后生成的SAM文件怎么查看里面的内容呢？

$ less -SN *.sam（sam文件的文件名称）

然后会显示如下内容：

头行（header line）以 @ 开始，紧接着一个或两个字母,比如下列代码中的 SQ 表示参考序列信息，SN表示参考序列名称，LN表示参考序列长度，PG表示软件，ID表示项目记录号（唯一），PN表示软件名称，VN表示软件版本，CL表示命令行等等。

SAM比对结果部分有11列，每一列都是不同的信息：

第1列：fastq的read ID
第2列：FLAG，对应的数值如下：
（如果某一个数值不是下面的任意值，那么那个数值就是下面这些数里面几个的和）

1：该read是成对的paired reads中的一个
2：paired reads中每个都正确比对到参考序列上
4：该read没比对到参考序列上
8：与该read成对的matepair read没有比对到参考序列上
16：该read其反向互补序列能够比对到参考序列
32：与该read成对的matepair read其反向互补序列能够比对到参考序列
64：在paired reads中，该read是与参考序列比对的第一条
128：在paired reads中，该read是与参考序列比对的第二条
256：该read是次优的比对结果
512：该read没有通过质量控制
1024：由于PCR或测序错误产生的重复reads
2048：补充匹配的read

比如说，我的比对结果里这一列的值有一个83。那么这个83并不在上述的值里，但是83是1+2+16+64的结果，那么这个read的比对结果的解读就是：
该read是成对read里的一条，该read反向互补序列可以比对到参考基因组上，并且和这read配对的read也能比对到参考基因组上，这条read是这一对read里的第一条。

第3列：染色体名称。如果这列是“ * ”，可以认为这条read没有比对上的序列，则这一行的第四，五，八，九 列是“0”，第六，七列与该列是相同的表示方法。

第4列：比对的位置，从对应上的染色体第1位开始往后计算。没有比对上的，此处为0。

第5列：MAPQ比对质量值。越高说明该read比对到参考基因组上的位置越唯一，例如42。Mapping qulity的计算方法是：Q=-10log10p，Q是一个非负值，p是这个序列不来自这个位点的估计值。假如说一条序列在某个参考序列上找到了两个位点，但是其中一个位点的Q明显大于另一个位点的Q值，这条序列来源于前一个位点的可能性就比较大。Q值的差距越大，独特性越高。如果值为255表示mapping值是不可用的，如果是unmapped read则为0。

第6列：简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report）。其以参考序列为基础，使用数字加字母表示比对结果，match/mismatch、insertion、deletion、skipped region from the reference(表示可变剪接位置)、soft clipping (clipped sequences present in SEQ)、hard clipping (clipped sequences NOT present in SEQ)。对应字母 M、I、D、N、S、H。比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的；例如：36M 表示36个碱基在比对时完全匹配。再比如：如37M1D2M1I，这段字符的意思是37个匹配，1个参考序列上的删除，2个匹配，1个参考序列上的插入。
（NOTE：clipped均表示一条read的序列被分开，之所以被分开，是因为read的一部分序列能匹配到第三列的RNAME序列上，而被分开的那部分不能匹配到RNAME序列上。而H只出现在一条read的前端或末端，但不会出现在中间，S一般会和H成对出现，当有H出现时，一定会有一个与之对应的S出现）

第7列: 这条reads第二次比对的位置。=表示参考序列与reads一模一样，*表示没有完全一模一样的参考序列

第8列: 该列表示与该reads对应的mate pair reads的比对位置（即mate），若无mate,则为0。
（NOTE:mate，在Illuminated中有两种测序技术：paired end sequencing，mate pair sequencing。这两种测序都是测的一个片段的两端，这两端产生的reads被称为mate1，mate2，单末端测序则无mate。）

第9列: 序列模板长度，如果同一个片段都比对上了同一个参考序列，为最左边的碱基位置到最右边的碱基位置（左为正，右为负）。当mate 序列位于本序列上游时该值为负值。不可用时，为0。

第10列: read的序列

第11列: ASCII码格式的序列质量。格式同FASTQ一样。

第12列: 可选的区域。
格式一般差不多是这样的：AS:i:-1 XN:i:0 XM:i:1 XO:i:0 XG:i:0 NM:i:1 MD:Z:35T0 YT:Z:UU
AS:i 匹配的得分
XS:i 第二好的匹配的得分
YS:i mate 序列匹配的得分
XN:i 在参考序列上模糊碱基的个数
XM:i 错配的个数
XO:i gap open的个数
XG:i gap 延伸的个数
NM:i 经过编辑的序列
YF:i 说明为什么这个序列被过滤的字符串
YT:Z 值为UU表示不是pair中一部分(单末端？)、CP(是pair且可以完美匹配)
DP(是pair但不能很好的匹配)、UP(是pair但是无法比对到参考序列上)
MD:Z 代表序列和参考序列错配的字符串

以上就是SAM文件的格式说明。我的这篇文章主要会focus on 在sam文件的第5列：MAPQ。因为后续我想做一个ATAC-seq的练习，那篇文献里方法部分提到他们把sam文件根据MAPQ过滤了一下，所以下面主要是学习MAPQ相关知识点。

参考文章：
1.生信人必会数据格式持续收集-测序原理-数据格式-数据库-生信技能树
2.SAM文件格式介绍 | Public Library of Bioinformatics
3.SAM文件格式说明 | 寂寞先生
4.生信：2：sam格式文件解读https://blog.csdn.net/genome_denovo/article/details/78712972

##############################我是分割线##################################

别人的MAPQ值和你比对出来的MAPQ值能直接拿来比较吗？

这篇文章给了一个很好的解答：
http://www.acgt.me/blog/2014/12/16/understanding-mapq-scores-in-sam-files-does-37-42

以下是这篇文章的一个大概的内容，并没有完全翻译：

序列比对图(SAM)格式文件每一列中都存储了相应的内容。其中，SAM文件的第五列存储比对质量(MAPQ)值。

MAPQ: MAPping Quality. It equals −10 log10 Pr{mapping position is wrong}, rounded to the nearest integer. A value 255 indicates that the mapping quality is not available.

按照上面的公式，如果某一个read的正确比对概率是0.99，那么它的MAPQ值应该是是20，即：-10×log10（1-0.99）。如果正确比对概率是0.999，那么MAPQ的值就是30。所以MAPQ的值取决于你的正确比对的概率。（如果MAPQ值是255，那么这个值不可用）。相反，当正确比对概率趋向于0时，MAPQ的值也趋于0。

在比对read后做的第一件事，就是统计sam文件里MAPQ值的分布。但是也有很多人并没有关心MAPQ值。也许你很相信比对软件输出的sam文件，但是这些分数真的会有很大差异吗？
下面的图是来自两个比对的MAPQ值的分布。下面的图是上面图的放大，可以更清晰地显示0-1之间MAPQ分数的分布:

从这个图里我们能得出什么结论呢？这两个比对有很明显的区别。experiment1最常见的MAPQ得分是42，其次是1。在experiment2中，得分最多的只有37分，其次是0分。实验1基于小鼠数据，实验2使用拟南芥数据。但这可能不是分布不同的原因。小鼠的数据是基于DNase-Seq实验中未配对的Illumina read，a . thaliana的数据来自于全基因组测序中Illumina读取的成对read。然而，这些区别仍然可能不是造成差异的原因。

造成这些MAPQ值分布的不同的真实原因，是实验1利用的bowtie2基因的比对，而实验2利用BWA基因MAPQ值的计算。所以你不应该比较这两个试验的MAPQ值，除非你用的是同一个比对软件。

对于bowtie2比对的sam文件，MAPQ值最大是42；而BWA比对出来的MAPQ值最大是37。

MAPQ的影响因素（参考认识MAPQ）：
（1）基因组重复区域MAPQ会比较低，因为会出现multiple mapping 和 reads聚集的情况；
（2）read 中碱基质量值，低质量值的碱基意味着序列很可能是错误的，错误的序列可能会导致错误的比对，所以MAPQ会低；
（3）比对算法的敏感性，如果比对算法敏感性差，会造成比对错误，MAPQ低；
（4）单双端测序的影响，如果reads两端都可以比对到基因组同一位置，那么比对正确的可能性很大，MAPQ会高。

所以，你需要注意的是：
（1）MAPQ值在不同比对软件的结果是不一样的。
（2）你应该根据你自己的MAPQ值来过滤“真正”不好的比对read。

##############################我是分割线##################################

Bowtie2是如何分配MAPQ值的呢？

上面说了不一样的比对软件，得到的MAPQ值并不一样，你也不能将它们直接拿来比较大小。BWA我就不学习了，因为我主要用的都是bowtie2，如果有同学需要请自行学习。这里我主要搜索一些关于bowtie2比对结果MAPQ值。
这里有一篇文章非常的详细：How does bowtie2 assign MAPQ scores?

下面就来学习一下这篇文章：

比对质量值(MAPQ或MQ)被bowtie2和bwa等软件来评估read比对到基因组的质量。公式是这样滴：

公式里的p代表一个read比对错误的概率

MAPQ的值从0到37,40或者42，取决于你用什么软件进行比对。这里只说一下bowtie2的MAPQ值。bowtie2的MAPQ的值并不用上面的公式来计算。
当我们分析NGS数据时，应该根据MAPQ的值来过滤一下，把低于某一个值的read剔除，但是应该选什么值来作为阈值呢？有人说（一个研究果蝇的同事）应该保留MAPQ值>=30的Read。但是在实际分析中，对于人类细胞系，这个值就不太合适了，因为人类细胞系里SNP，microdeletion，microinsertions，breakpoints等等有很多，导致了比对质量值会偏低。有的人认为MAPQ>=10的read都可以保留，甚至还有人认为只要MAPQ>=1都是可以接受的。还有人说MAPQ值为255的read是unique比对，但这是根据旧的定义来说的，新的定义在SAM官网称MAPQ值为255的read是不可用的。为了搞清楚用bowtie2比对到底应该怎么过滤，这篇文章的作者设计了一个小实验来确定。

这个具体的实验过程我就不详细说了，各种的代码。。。各种公式。。。比较上头，我们来直接看结论吧：

在bowtie2里，第12列的可选字段里，真正的multiread（AS=XS）也可以得到MAPQ=1（如果AS == XS，则认为这个read是真正的multiread，并且MAPQ只能得到0或1），不管这个read比对到基因组多少个位置。当read比对发生0或者1次错配，那么AS=XS将为-6。像这样：

AS:i:-6 XS:i:-6 MAPQ=1

如果有2-5个错配，那么结果是这样的：

AS=XS <= -12 (i.e. -12 to -30.6) MAPQ=0

所以，当你想从你的data里排除“真正的multireads”时，用MAPQ>=2也是可以的。对于高质量的比对结果，MAPQ >= 3代表允许3个错配，MAPQ >= 23代表允许2个错配，MAPQ >= 40允许1个错配，MAPQ >= 42代表允许0个错配。

而在bowtie2里，真正的uniread 可以得到不同的MAPQ值，例如3,8,23,24,40,42。如果你只想保留uniread，那么你就可以只保留MAPQ为上述这些值的reads。比如你可以用下面这个代码：

$ awk '$5 == 3 || $5 == 8 || $5 == 23 || $5 == 24 || $5 == 40 || $5 == 42' file.sam

或者：

$ grep -v XS:i: file.sam

那么如果你想根据某一个MAPQ的值来过滤你的sam文件：

#如果你想把MAPQ小于2的sam文件都丢掉，并转成bam文件
$ samtools view -bSq 2 file.sam > filtered.bam
##-q INT Skip alignments with MAPQ smaller than INT [0]

参考文章：
1.Question: Filtering A Sam File For Quality Scores
2.bowtie2