转自公众号：universebiologygirl
http://mp.weixin.qq.com/s?__biz=MzU1NDkzOTk2MQ==&mid=2247485204&idx=1&sn=f114e247aa39f0002d7bb4b1ea627593

BAM/SAM文件两三事

写在前面的

上期推文主要解读了BAM文件的Flags值，通过一个perl程序快速批量翻译Flags值。flags详细记录了reads的比对信息。是由一堆2的n次方（n=0~11）加和组成。例如某条reads经过比对后的flags是83, 83= 1+2+16+64转化为2进制为1010011， 则代表这是一条双端测序的Read1，该reads及其对应reads都是正常比对且当前reads是反向互补后比对上的参考基因组。而perl程序就是基于这个基本思路，通过构建10进制转2进制进行Flags值计算的。虽然已经有flags值计算器了。但是通过编写这些程序会加深对BAM文件的理解。从而加深对生物信息学数据的理解。

CIGAR值的解读

CIGAR是Compact Idiosyncratic Gapped Alignment Report的首字母缩写，是也是记录了reads比对到参考基因组上的字符串。只不过这个字符串是由数字+字母组合成的字符串。与Flags值不同的是CIGAR记录的是每条reads内部每个部分的比对情况。而Flags值则是记录每条reads总体的比对情况。与Flags值相比，CIGAR值对于reads比对的细节记录的更详细。那么，在对CIGAR值有了初步的了解以后，我们继续深入了解这个值的意义。由于CIGAR值是由数字+字母构成的字符串，那么下面这张来自于碱基旷工的图则清晰的阐明了CIGAR值中出现的每个字母的比对含义。

举例而言，如果一条250bp的reads在BAM文件中的CIGAR值记录为50S100M10D100M，那么则表示这条reads的前50bp被跳过，中间100M比对上了参考基因组，后100bp也比对上了参考基因组。但是中间100M和后面100M则在比对参考基因组时出现了一个10bp的gap。也就是参考基因组上有10bp是这条reads没有的，参考基因组这10bp两端各100bp的序列恰好与reads有匹配，也就是这10bp对于reads来说是一个潜在的gap。而这个被跳过是什么意思呢？我自己用自己的数据测试了一下，有两种情况，一种是这部分序列比对不上参考基因组，另外一种情况是它可以比对到参考基因组但是和reads剩余部分比对到基因组上的距离相比相差太远。而经常出现于RNA-seq中的N一般也就几十到几百bp那么仔细想想，会不会出现50S10D200M的情况呢？其实是不存在的，因为这条reads的前50bp被跳过参考基因组。后200M匹配到参考基因组，如果前面还有10D缺失，那么BAM文件另外一列的比对position值将改变，因为比对position值指的是reads比对Match到参考基因组的第一个位置。所以CIGAR值的标记字符不会像随机排列组合那样所有组合情况都有。下面我们继续追踪这样的问题，对于reads sequence来说，我们更希望直观明了地找到这些reads哪些部分是M，哪些部分是S，哪些部分又是D，D可以用-表示。接下来我们将通过写程序的方式找到每个reads标记字符对应的子序列。我们思考这样的情况，对于D的出现，它只能出现在子序列之间，且D所表示的区域在reads上是无子序列的，只能填-充当gap。对于N，在reads中也是不存在的。对于S,M,I这些序列reads中是实际存在的，由于标识字符前对应的数字就是实际长度，所以只要正确找到这些子序列的在 reads上的起始位点。就可以正确找到这些子序列。显然，reads中D或者N的存在是影响寻找SMI等的起始位点的。如果D或者N出现，那么后面的SMI的起始位点计算则不能把D或N的长度算入。那么这样的模型我们可以用动态规划算法，令stRi代表第Ri个reads子序列的起始位点，由于D或N不可能直接出现在reads起始部位，我们仅Ri-1的标识字符情况，用公式表示则为

上面模型主要阐述为如果某条子序列前面出现了一个D或者N，那么它实际的起始位点则是它在D或N前面的对应子序列结束的位置后，如果它前一个子序列是SIMH，则它的实际起始位点则在上一个子序列结束位置后，起始对于H来说，它已经将reads截断了，后面不可能有序列了，所以下面的程序起始多余了一个判断。
同样编写perl语言完成这个程序。

定义一个CIRAG_count子程序，该子程序主要用于计算每个子序列的起始位点

sub CIRAG_count {

这个就是利用perl程序提取的子序列

总结
核心目的主要是理解CIGAR值的实际生物信息学含义。了解动态规划算法的一种简单运用。
参考
碱基旷工
SAMtools手册
作者信息
熊东彦，中国科学院武汉病毒研究所在读研究生，擅长方向R语言编程，Perl语言编程，转录组分析，宏基因组分析。近期推文：生物信息学实用小技巧。
往期精彩推文
【BioLeetCode】生信编程挑战与github项目维护邀请
医学遗传学数据库小结分析

                                                        
《生物信息学分析实用小技巧（六）：BAM/SAM中的雪茄值》来自互联网公开内容，收录仅供学习使用，如侵权请联系删除。本文URL：https://www.ezixuan.com/1020478.html