二代测序获得的是bcl格式的原始下机数据，通过bcl2fastq软件 https://support.illumina.com/sequencing/sequencing_software/bcl2fastq-conversion-software.html 可将bcl文件转换成每个样本的fq格式文件，也就是我们常理解的数据拆分。bam文件是由比对软件将质控后的fq格式文件与参考基因组进行比对后的比对信息存储文件。
接下来我们理解下bam文件的内容。参考原文提出的一张经典图片：
上图格式的查看方法为：
samtools的header信息每一行都用‘@’ 符号开头，一般大家不会太关注，但其中的信息对于我们有些生信分析还是很重要的。这里需要重点提一下的是header中的@RG也就是Read group信息，这是在做后续数据分析时专门用于区分不同样本的重要信息。比如测序多条lane获得的bam的合并：如果原来样本的测序深度比较深，一般会按照不同的lane分开比对，最后再合并在一起，那么这个时候你会在这个BAM文件中看到有多个RG，里面记录了不同的lane，甚至测序文库的信息，唯一不变的一定是SM的sample信息，这样合并后才能正确处理。这个合并当然也可以在数据拆分后对rawdata进行cat合并，然后再生成bam文件。
接下来是bam的主体内容record（有时候也叫alignment section，即，比对信息），每 一行代表一条reads，每条reads的信息用tab键进行分隔：
对于每列的解释如下表所示：
森或 比如十进制数据77 = 000001001101 = 1 + 4 + 8 +64，这样就得到了这个FLAG包含的意思：PE read，read比对不上参考序列，它的配对read也同样比不上参考序列，它是read1。
二进制的质量描述见下表：
- MAPQ：比对质量值，这个是大家最为熟悉的比对质量值了。比如说Q30（错配率为0.001），Q20（错配率为0.01），计算公式为：-10logP{错比概率} 。一般结果是这一列的数值是从0到60，且0和60这两个数此颂伍字出现次数最多。
- CIGAR：该标签采用数字和几个字符的组合形象记录樱伍了read比对到参考序列上的细节情况，读起来要比FLAG直观友好许多，只是记录的是不同的信息。比如，一条150bp长的read比对到基因组之后，假如看到它的CIGAR字符串为：33S117M，其意思是说在比对的时候这条read开头的33bp在被跳过了（S），紧接其后的117bp则比对上了参考序列（M）。这里的S代表软跳过（Soft clip），M代表匹配（Match）。N表示可变剪接位置，常见于RNA-seq。H只出现在一条read的前端或末端，但不会出现在中间，S一般会和H成对出现，当有H出现时，一定会有一个与之对应的S出现。CIGAR的标记字符有“MIDNSHP=XB”这10个，分别代表read比对时的不同情况：
详细内容请参考： http://samtools.github.io/hts-specs/SAMv1.pdf
A coordinate system where the rst base of a sequence is one. In this coordinate
system, a region is specied by a closed interval. For example, the region between the 3rd
and the 7th bases inclusive is [37]. The SAM, VCF, GFF and Wiggle formats are using the 1-based coordinate system.
A coordinate system where the rst base of a sequence is zero. In this
coordinate system, a region is specied by a half-closed-half-open interval. For example, the region
between the 3rd and the 7th bases inclusive is [27). The BAM, BCFv2, BED, and PSL formats are using the 0-based coordinate system.
由测序仪测序所得或由原始序列组装所得的DNA/RNA序列
一段连续的序列或者子序列
一段由测序仪测序所得的原始序列。一条Read可能由多个片段组成，在测序数据祥裤中，reads是根据它们被测的顺序来建立索引的。
一个Read单向地比对到参考基因组上，这个比对结果中可以有插入、缺失、跳跃等，但是不能存在“双向”的比对结果，即Read的一段比对到正链参贺滑考基因组、一段匹配到负链，这种方向切换是不允许的，在SAM文件中，线性比对的特性就是：只用一行来记录。
就是当一条Read对比时，比对到了多个区域，但是这些区域并没有重叠的部分，也即由多个“线性比对”结果组成了一个集合，这个集合就组成了一个嵌合比对，嵌合比对中只有一个“线性比对”结果是谨拍简具有代表性的，其余的都以补充的身份出现，嵌合比对的特征就是多个“线性比对”记录中的Read对应的Qname（Read的名字，每个Read只有一个Qname）都是相同的，且这些“线性比对”集合中的每个记录的flag值都是一样的。
无论是上面提到的线性比对还是嵌合比对，只要能够完整的表现出一条Read的对比情况，就是一个Read 比对。
由于序列的重复性，导致一个Read在比对时会被比对到多个区域上，其中只有一个比对质量最好的会被当做比对结果的代表性结果，目前来看，这种决定方式还不是很严谨。
多次比对和嵌合比对是有根本性区别的，多次比对是因为序列本身具有重复性引起的，属于正常比对结果，而嵌合比对更多是由于实验、测序、结构突变、融合、以及其他因素引起的。
一个可能性区间，如果一个碱基被检测正确的正确率是99%，那么错误的概率就是1%=0.01=10（^-2），此时phredscore=-10log10（ 0.01 ）=20，
SAM文件由两部分组成，头部和主体，都以tab分列。
头部内容主要以各种说明为主，比如说明所用软件啦，参考基因组信息啦，排序信息啦等等，下面表格是SAM文件中涉及的一些专有名词的解释。
头部内容说明信息均是以@符合开头的，在头部是以“@说明类型码TAB键TAG：Value” 开始的，每个标签和说明类型码都是由两个字母组成的，下面将列举说明类型码和TAG
HD：代表意思为：SAM文件的开头标志、一般只要出现就会在第一行
HD中的标签（TAG）有：
VN*： 注释版本信息
SO： 比对结果的排序类型：有unknown、unsorted、queryname、coordinate四种排序类型
GO： 比对结果的分组信息：相似的比对结果会被分到一组，这里的分组结果中并不是需要全部进行过排序的，排序的类型有：none (default)、query (alignments are grouped by QNAME)、reference (alignments are grouped by RNAME/POS 三种类型
SS： 子排序类型，比如在某次算法中需要根据coordinate排序，而在每个coordinate排序的结果中又根据QNAME进行了排序，则在SAM文件中就应该表示为：@HD SO:coordinate SS:coordinate:queryname.如果在SO的基础排序中，排序的类型不在之前定义的四种排序类型中时，则SO对应的就应该时unsorted，此时SS就会起主要作用，比如，如果基础排序是根据一个辅助标签MI排序的，之后又根据coordinate排序的，则在SAM头部中就应该表现为@HD SO:unsorted SS:unsorted:MI:coordinate.
SQ： 说明类型码代表的含义为参考序列字典， SQ 的排序决定了比对结果的排序顺序
SQ中又以下标签：
1.SN*：参考序列名称
2.LN*：参考序列长度
3.AH、AN、AS、DS、M5、SP（种族）、TP、UR，这些不常用
例如：@SQ SN:JF-PLAC8_CT_converted LN:88
Read Group，每个Sample都有一个RG ID，一个Sample可以在多个库中进行测序。
RG对应的有以下TAG：
ID* ：Read group的唯一ID
BC：辨别样本或文库的标签序列
SM：样品名
LB：文库名
PU：测序仪
PL：测序平台
CN：产生read的测序中心的名称
PG： 程序 说明类型码,PG中有以下几个标签：
ID*： 程序记录标识，每个程序记录都只有一个ID
PN： 程序名称
CL： 命令行内容（utf-8编码）
PP： 好像是指前一个 @PG-ID：不怎么用这个，有了解的大哥大姐可以帮忙做个备注
DS： 描述
VN： 所用程序版本
以上就是五个说明类型码代表的意思以及常用的标签意义。下面截图作为例子，可以对照上面的内容看一下：
SAM文件主体内容中主要有11列内容，这11列中的内容就如下表：
当数据中这11列内容哪个对应内容有缺失或者无效的话就用‘0' 或 `*'代替。
1.QNAME： Read ID，就简单理解成Read的名称，其中会包含一些测序平台信息
2.FALG： 可以理解为比对结果的标志，可以根据FLAG值筛选比对结果
FLAG值若记不住可以直接在<https://broadinstitute.github.io/picard/explain-flags.html >中进行查询。在SAM文件中出现的flag值是涉及到的value相加得到的值，比如99，177等，可以在上述网站查询
3.RNAME： 理解成比对上的染色体号即可
4.POS： 比对到参考序列上的位置
5.MAPQ： 比对质量值，算法与文中第一部分Phred score算法一致：-10log10（错误率），若值为255表示其比对结果不可用，如果是unmapped read则MAPQ为0
简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report），其以参考序列为基础，使用数字加字母表示比对结果，比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一 个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的；
M”表示 match或 mismatch；
“I”表示 insertion；
“D”表示 deletion （表示的是READ和ref相匹配时，参考基因组中需要deletion的部分，不是READ）；
“N”表示 skipped （跳过这段区域）；
“S”表示 soft clipping （被剪切的序列存在于序列中）；
“H”表示 hard clipping （被剪切的序列不存在于序列中，已经在除去低质量READ的时候被过滤了）；
“P”表示 padding（填充） ；比如参考基因组序列本来为ATTACGAC，read序列为ATAATACGAC，那么，如果在比对的sam结果文件中，有两种比对方式，如下：
1.如果reference是padded reference，也即参考基因组展示为AT**TACGAC，那么**就是代表了ref考虑了read序列的插入情况，但是这种表示情况下CIGAR中就不会再出现p和I（Insertion）了，因为参考基因组中已经考虑了插入的情况了。而如果有其他read跟这个padded ref比对时，对于ref中pad的区域没有序列的话，就会以D（deletion）来表示这个read。
2.另一种情况就是unpadded reference，也即reference是正常的参考基因组（可参考 https://blog.csdn.net/xubo245/article/details/51283022 ），一般情况下，以这一种为主。
“=”表示 match；
“X”表示 mismatch （错配，位置是一一对应的）；
7.RNEXT： mate序列匹配上的染色体号（比如，Read2对应的染色体号），如果第三列和这列都是“*”则说明此Read没有匹配到任何一个染色体，如果第三列有信息，这列是“”或者“=”则代表此Read匹配到第三列对应的染色体上。
8.PNEXT： 该列表示与该reads对应的mate pair reads的比对位置，如果这对pair-end reads比对到同一条reference序列上，在sam文件中reads的id出现2次，Read1比对的第4列等于Read2比对的第8列。同样Read1比对的第8列等于Read2比对的第4列。例如：
第1列（Read id）····第4列（Read1比对位置）····第8列（mate-pair reads比对位置）
22699:1759····124057649····124057667
22699:1759····124057667····124057649
相同的reads id一个来自Read1文件，一个来自Read2文件，第4列和第8列是对应的
9.TLEN： signed observed Template LENgth （可以理解为文库插入片段长度）
如果R1端的read和R2端的read能够mapping到同一条Reference序列上（即第三列RNAME相同），则该列的值表示第8列减去第4列加上第6列的值，R1端和R2端相同id的reads其第九列值相同，但该值为一正一负，R1文件的reads和R2文件的reads，相同id的reads要相对来看。在进行该第列值的计算时，如果取第6列的数值，一定要取出现M的值，S或H的值不能取。
10.序列信息 ，
11.Read质量信息 （ASCII编码）
12.可选区域
举几个例子：NM，MD
NM可以简单的理解为：如果要把READ变为跟reference一致需要几步（一般步骤有三种类型，碱基替换，删除碱基，添加碱基）
MD：简单的理解为匹配结果（跟第10列结合可以看出reference的原序列），下面举个例子：
从上面的图可以看出来这里的NM值为2，MD为3^A42T5，由MD可以知道这个read比对结果中前3个匹配上了，之后需要在read中插入一个A（也就是read相比ref少了一个A），再之后42个碱基是匹配上了，到然后这42个之后的碱基在参考基因组上应该是T，后面5个是匹配上了。这样的话我们根据read的序列就能推断出ref的序列为：CCGATAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCC。
我们检查一下推测的是否正确，来看hg19中这段序列：
可以发现除了开头四个碱基在hg19中是N以外，其他的是匹配一致的（好像是因为有写序列是在端粒中的，所以呈现的序列信息就是N，也有些是因为还没有测出来，所以用N代替。）、
这样再来看的话，就可以知道read相比ref有两处不同，一个是缺了个A，一个是碱基错配（T-A），所以如果read要和ref一致的话需要编辑两次，一次是添加一个A，一次是把A变成T。
如果觉得有用，就帮我点个赞吧^_^!
每次遇到BAM文件flag值都有一些困惑，时间越或运久越迷惑。
在此，针对BAM文件中的flag信息进行梳理和解释：
记录于BAM文件的第2列，以巧芦 bwa 软件比对结果为例
可以使用samtools查询：
samtools view test.bam | cut -f2 | uniq
问题来了，这些数字的意义是什么呢？
首先可以参考SAM/BAM文件的解释文档：
http://samtools.sourceforge.net/SAMv1.pdf
其中，对于FLAG有如下说明：
FLAG: bitwise FLAG. Each bit is explained in the following table:
| Bit | Description |
| 0x1 | template having multiple segments in sequencing |
| 0x2 | each segment properly aligned according to the aligner |
| 0x4 | segment unmapped |
| 0x8 | next segment in the template unmapped |
| 0x10 | SEQ being reverse complemented |
| 0x20 | SEQ of the next segment in the template being reversed |
| 0x40 | the first segment in the template |
| 0x80 | the last segment in the template |
| 0x100 | secondary alignment |
| 0x200 | not passing quality controls |
| 0x400 | PCR or optical duplicate |
| 0x800 | supplementary alignment |
上述0x1, 0x2, …是 十六进制 的数值与十进制的数字截然不同。
对应的十进制数值描述如下：
| 十进制 | 描述 |
| 1 | template having multiple segments in sequencing |
| 2 | each segment properly aligned according to the aligner |
| 4 | segment unmapped |
| 8 | next segment in the template unmapped |
| 16 | SEQ being reverse complemented |
| 32 | SEQ of the next segment in the template being reversed |
| 64 | the first segment in the template |
| 128 | the last segment in the template |
| 256 | secondary alignment |
| 512 | not passing quality controls |
| 1024 | PCR or optical duplicate |
| 2048 | supplementary alignment |
回过头来看，比如16和1024分别是比对到互孝团带补链的片段，对于1024指的是PCR重复片段。
那其他数字的含义呢，他们只是简单数字组合而已，例如：1040是1024 + 16，Read比对到反义链且是一个PCR重复，简单的数字相加而已。
也可以借助 flag解释链接 来解析上述数字的含义，如把1040输入到该网站会返回：
“read reverse strand”和“read is PCR or optical duplicate”。
不过，SAM说明文档中FLAG的代号均使用按位符号显示。bit是信息的基本单元且只有2个数值，1和0。
这谁能搞的懂啊？！
直接用linux bc转换吧：
按下表对10000010000从右到左依次读取：
故BAM flag转换为元字符，轻松获取各种类型BAM flag值背后的信息。
参考资料
https://davetang.org/muse/2014/03/06/understanding-bam-flags/
欢迎分享，转载请注明来源：内存溢出
微信扫一扫
支付宝扫一扫
评论列表（0条）