搜索
查看: 3064|回复: 0

【直播】我的基因组(十三):了解sam格式比对结果

[复制链接]

103

主题

133

帖子

860

积分

版主

Rank: 7Rank: 7Rank: 7

积分
860
发表于 2017-2-2 09:49:31 | 显示全部楼层 |阅读模式
本帖最后由 zckoo007 于 2017-2-2 09:52 编辑

【直播】我的基因组(十三):了解sam格式比对结果

很抱歉这么久都没有推直播给大家了!每到年终的这个时候,小编真的是忙的找不到北大/(ㄒoㄒ)/~~。这一周的直播应该会每天一更啦,希望大家可以跟着一起学习,不要脱粉哦~~

另外还要谢谢大家在生信菜鸟团困难时候的帮助!再次表达感激之情!


言归正传,十一讲中将我们主要讲了如何将下机数据比对到参考基因组中。但是很多人对比对结果却是一头雾水。那我们现在来了解一下Sam格式的比对结果吧!

比对工具到现在已经多如牛毛了,见列表:[url=]https://en.wikipedia.org/wiki/List_of_sequence_alignment_software[/url] 。但是能被大多数人熟知的,就是bowtie和bwa(我们在十一讲中用的才是bwa),它们把测序数据比对到参考基因组之后,都会生成一个sam格式的文件。随后的大部分分析都是基于sam格式进行的分析,虽然Jimmy多次强调这些基础知识的重要性需要大家私下自学。但是由于这个sam文件实在是太重要了!!!所以,不得不亲自抽出一讲来说说它,后面也会基于此写十多篇文章:



目录
14-把bam文件给按照染色体给分割成小文件
15-提取未比对的测序数据
16-提取多比对的测序数据
17-提取左右端测序数据比对到不同染色体的PE reads
18-去除PCR的duplication情况
19-根据比对结果来统计测序深度和覆盖度
20-覆盖度累积曲线

因为这个是基础,如果你后面的十几篇有不理解的,请回头来再仔细看看sam文件的定义!

当然,不仅是这些分析是基于对sam文件的理解,我只是举几个例子,大家千万要熟练使用sam格式的比对结果,最权威的定义见:[url=]https://samtools.github.io/hts-specs/SAMv1.pdf[/url]

记住,我们的双端测序的数据,一个paired reads,有左右两端两条reads,所以在sam文件里面会有且只有两条记录,除非你设置特殊参数,允许输出多比对情况。


a7507f723939262a9fafa6802a4a0e58.jpg

上面是一个典型的PEreads输出的sam比对结果,反正必须要有的就是下面11列,其中第3和第7列,可以用来判断某条reads是否比对成功到了基因组的染色体,左右两条reads是否比对到同一条染色体。而第1,10,11列可以提取出来还原成我们的测序数据fastq格式的。第9列是我们建库的时候打断的片段长度,本次是PE150的数据,打断成350bp,所以这里应该是350个字符左右,但如果是RNA-seq数据,就不一样了。


1ca4dbdff020d64a6ca0709c88298ab0.jpg

其中第二列flag是比较反人类的,一般人用不了二进制,有网页可以帮助你:[url=]http://picard.sourceforge.net/explain-flags.html[/url]。我们的sam里面第二列是下面这些二级制转为十进制后的和!


c939b99c5f345f6c67598c0aaa2f5f49.jpg

然后第6列CIGAR是比较重要的,解释如下,其中M并不是说match,所以我们的PE 150的reads,大部分都会是150M,但是并不代表着跟参考序列一模一样。其中S/H是比较特殊的,很难讲清楚,但是大部分情况下用不到。(soft-clipping碱基是指一条reads未匹配上当前基因组位置的部分,如果有多个reads在这种情况并且这些reads的soft-clipping碱基都能够比对在基因组另一位置,那么就可能存在SV)


2dc62989da20bdbe33736b1268f24d26.jpg

第5列,比对结果的质量值,也是因工具而异。
a. Match score: Score awarded for a base in a sequence matching a base in another sequence
b. Alignment score: Cumulative score of the bases of a sequence matching the bases of another sequence (more this score, better the alignment, if all else equal)
c. Mapping Quality score: Probability that the shorter sequence is mapped to the right spot on the longer sequence.

如果定义某条reads比对的质量值是一个非常复杂的问题,我也没办法说清楚,感兴趣的朋友可以去查看 [url=]http://biofinysics.blogspot.com/2014/05/how-does-bowtie2-assign-mapq-scores.html[/url]

但是需要记住,质量值越高这个比对越可信,如果质量值为0,可能是该序列在参考基因组有多种定位的可能性。

最后,一般来说,sam文件肯定是大于11列的,后面多余的列是各种各样的 tag。而且只要是你开发了一个比对工具,你就可以定义一堆tag,这个并没有公认的标准,因为sam文件的定义就是前面的11列,后面的tag是随心所欲的!

但是一般RG代表着你的sam文件比对来自于哪个样本的fastq程序结果。NM这个tag是编辑距离,大概就是你的reads如果想转变成参考基因组,需要改变多少个碱基,如果编辑距离是0才说明你的这个150bp长度的序列跟参考基因组一模一样。

MD这个tag里面写明了,你的序列跟参考基因组不同在哪里,比如下面的截图里面的,我的某个位点相比参考基因组来说,就变成了G,而其余的碱基都是一样的。

AS和XS在两个标签貌似没什么用,以后再说吧。

如果你用的bowtie或者hisat等其它比对工具,还会有更多的稀奇古怪的tag,学无止境呀!


0eb0499ba21579ca6c881d708ffdcea9.jpg

文:Jimmy、吃瓜群众

图文编辑:吃瓜群众






上一篇:【直播】我的基因组(十二):先粗略看看几个基因吧
下一篇:【直播】我的基因组(十四):bam文件给按照染色体给分割...
基因组,转绿组,肿瘤信息,生物统计,Python, Linux.
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树 ( 粤ICP备15016384号  

GMT+8, 2020-6-1 13:14 , Processed in 0.027687 second(s), 30 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.