RNA-seq或者ChIP-seq等等测序的上游分析流程里的比对步骤相信大家都知道,我之前也只是按照各种教程去走流程,并没有仔细的研究过每一步的内容。今天这篇文章学习一下sam文件的格式,以及如何根据read比对的质量来过滤你的sam文件。
一般比对后生成的SAM文件怎么查看里面的内容呢?
$ less -SN *.sam(sam文件的文件名称)
然后会显示如下内容:
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比对到基因组的质量。公式是这样滴:
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