bwa全名Burrow-Wheeler Aligner，是一款将短测序序列比对到参考基因组的软件。主要包括三种算法，分别是BWA-backtrack, BWA-SW和BWA-MEM，其中BWA-backtrack针对的是illumina 100bp或更短的reads，后两者针对的reads片段较长，跨度范围70bp-1Mb(大片段比对效率很低，可以考虑使用目前三代数据的比对软件)。相比于BWA-SW，BWA-MEM在运行速度及准确度上都有明显的提升，是目前使用最多的比对方法。具体使用方法参照之前的文章，BWA参考基因组比对-序列比对。
今天为大家介绍的是比对结果sam/bam文件CIGAR中的信息。
数据中CIGAR信息
cigar中包含的是比对结果信息，表明了一条reads所有碱基的比对情况，在sam/bam文件的第六列，示例信息如下。
@HD VN:1.6 SO:coordinate
@SQ SN:ref LN:45
r001 99 ref 7 30 8M2I4M1D3M = 37 39 TTAGATAAAGGATACTG *
r002 0 ref 9 30 3S6M1P1I4M * 0 0 AAAAGATAAGGATA *
r003 0 ref 9 30 5S6M * 0 0 GCCTAAGCTAA * SA:Z:ref,29,-,6H5M,17,0;
r004 0 ref 16 30 6M14N5M * 0 0 ATAGCTTCAGC *
r003 2064 ref 29 17 6H5M * 0 0 TAGGC * SA:Z:ref,9,+,5S6M,30,1;
r001 147 ref 37 30 9M = 7 -39 CAGCGGCAT * NM:i:1
CIGAR比对类型
对于cigar信息，sam文件中共定义了9种类型，见下表。
bwa-samtools-bam/sam-cigar
在实际的分析中，一般使用MIDNSHP七种，其中M表示MATCH，I表示INSERTION，D表示DELETION，N表示skipped bases on the reference，S表示SOFT CLIPING，H表示HARD CLIPING，P表示PADDING，以下是相关的示例信息及简要说明。
samtools-cigar-example
(a)最上面是参考基因组，下面给出了reads的比对信息，(b)bam文件中输出的reads比对结果。M、I和D三种类型很简单，都是相对于基因组而言的，对应匹配、插入和缺失；S和H两种类型稍后详细说明；对于P的理解可以看r001和r002两条reads的比对结果，r001插入了两个碱基AG，cigar信息为2I，r002插入了一个G碱基，对应cigar信息为1P1I，这里的P其实是r002相对于r001而言的补位信息，即r001中插入的A碱基在r002中是缺失的，进行补齐；N表示忽略基因组上的序列。
Clipping说明
在Smith-Waterman比对时，短序列的两端可能会出现不能匹配的问题，此时比对结果就会指定为cliping。不同的clipping结果会根据实际的比对情况进行细分，其中S对应soft，h对应hard。
clipping比对结果：
REF: AGCTAGCATCGTGTCGCCCGTCTAGCATACGCATGATCGACTGTCAGCTAGTCAGACTAGTCGATCGATGTG
READ: gggGTGTAACC-GACTAGgggg
上述示例是某一条read的比对结果，其中大写字母表示匹配(非完全匹配，部分碱基是错配的)，-表示缺失，小写字母表示末端为匹配的序列，这部分就是clipping序列。若该read只比对到基因组的这个位置，cigar信息为3S8M1D6M4S，若该序列比对到基因组多个位置，比对的cigar信息为3H8M1D6M4H。S和H除了比对位置的区别以外，在输出数据中的序列也不同，标注为S的序列会显示在bam文件中，标注H的序列则会删除。比如3S8M1D6M4S在bam中输出序列为gggGTGTAACCGACTAGgggg，而3H8M1D6M4H输出的序列为GTGTAACCGACTAG
参考信息
[1] file format：http://samtools.github.io/hts-specs/
[2] The Sequence Alignment/Map format and SAMtools (doi: 10.1093/bioinformatics/btp352)
[3] SAM information：https://davetang.org/wiki/tiki-index.php?page=SAM
网友评论