理解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

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 202,607评论 5 476
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 85,047评论 2 379
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 149,496评论 0 335
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 54,405评论 1 273
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 63,400评论 5 364
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 48,479评论 1 281
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 37,883评论 3 395
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 36,535评论 0 256
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 40,743评论 1 295
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 35,544评论 2 319
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 37,612评论 1 329
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 33,309评论 4 318
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 38,881评论 3 306
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 29,891评论 0 19
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,136评论 1 259
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 42,783评论 2 349
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 42,316评论 2 342

推荐阅读更多精彩内容