Map with BWA-MEM(BWA比对,mem方法)

分析模块,输入构建好BWA索引的参考序列(由分析模块“BWA reference index”生成),以及测序原始数据(fastq文件)。这里,推荐测序原始数据,先通过分析模块“Trimmomatic PE/SE”进行去接头污染和质量控制。

分析模块,利用bwa mem算法(适合中、长序列),将测序原始数据(fastq文件)比对回参考序列,生成sam格式的比对结果文件。

注:sam文件为文本格式文件,可用Windows写字板打开查看。这里,建议先将sam文件转成bam文件,bam文件排序之后,利用IGV软件打开,查看比对结果。

IGV安装和使用,包含Windows桌面版和iPad版,官方网站提供了详细的文档。参考网站:(http://www.broadinstitute.org/igv/)。


BWA输入测序数据分两种模式(SE/PE):

选择Paired时,分析模块会对双末端数据进行比对,用户需提供2个fastq原始数据文件,分别对应DNA片段的左端和右端测序结果。

选择Single时,分析模块会对单末端数据进行比对,用户需提供1个fastq原始数据文件。


输入:


1、构建好BWA索引的参考序列文件,由分析模块 "BWA reference index" 生成


2、fastq格式的测序原始数据文件。

示例:

@ecoli_section_9642_10205_3:0:0_1:0:0_0/1

CCCGCCATCTCTTGCAGAAGCGCCTGTTGCTGTACATGGTGCATTCGCATCCCCATCCCTACGCGGCTTC

+

2222222222222222222222222222222222222222222222222222222222222222222222


输出:


sam格式的比对结果文件。

示例:

@SQ          SN:ecoli_section       LN:30000

@RG         ID:ecoli_A1      SM:ecoli_A1

@PG          ID:bwa     PN:bwa    VN:0.7.12-r1039

ecoli_section_9642_10205_3:0:0_1:0:0_0       83     ecoli_section     10136       60     70M =       9642         -564         GAAGCCGCGTAGGGATGGGGATGCGAATGCACCATGTACAGCAACAGGCGCTTCTGCAAGAGATGGCGGG        2222222222222222222222222222222222222222222222222222222222222222222222     RG:Z:ecoli_A1 XT:A:U     NM:i:1         SM:i:37    AM:i:37   X0:i:1       X1:i:0       XM:i:1      XO:i:0       XG:i:0       MD:Z:56G13

ecoli_section_9642_10205_3:0:0_1:0:0_0       163  ecoli_section     9642         60     70M =       10136       564         CAGATGTCGATCCCATCGCCATTGCTGAAATCGGCATCCTTACCAGGGCAATTTTCGTGTTCCGTTTCAT        2222222222222222222222222222222222222222222222222222222222222222222222     RG:Z:ecoli_A1 XT:A:U     NM:i:3         SM:i:37    AM:i:37   X0:i:1       X1:i:0       XM:i:3      XO:i:0       XG:i:0       MD:Z:30G2C4A31

ecoli_section_3782_4324_1:0:0_3:0:0_1         83     ecoli_section     4255         60     70M =       3782         -543         GAGCCCGAATAAACGGTCTCAGCCAGTTAAGGCACTCCGACAAGAAATTAGCCCAGCTAATGGGGCATTG        2222222222222222222222222222222222222222222222222222222222222222222222     RG:Z:ecoli_A1 XT:A:U     NM:i:3         SM:i:37    AM:i:37   X0:i:1       X1:i:0       XM:i:3      XO:i:0       XG:i:0       MD:Z:2T51C2A12


SAM格式说明:

SAM是一种序列比对格式标准,由Sanger制定,是以TAB为分割符的文本格式,主要应用于测序片段(Query)比对到参考序列上的结果保存。

关于SAM格式的详细介绍,参考文档:(https://samtools.github.io/hts-specs/SAMv1.pdf)。或者Google搜索SAM format

SAM文件分为两部分:注释信息部分(header section)和比对结果部分(alignment section)。注释信息部分属于可选部分,以@开头,用不同的字母表示不同的信息。主要有:@HD,表明SAM文件的格式版本、比对结果的排序方式;@SQ,表明参考序列的信息,如染色体名称,序列长度等;@RG,表明测序序列(Query)的信息,如文库名、样品名、测序平台等;@PG,表明使用的程序信息,如程序名称、版本等;@CO,任意的注释信息。

比对结果部分(alignment section),每一行表示对应Query的一个比对结果信息。包括11个必须的字段(mandatory fields)和一个可选的字段,字段之间用TAB进行分隔。必须的字段有11个,顺序固定,如果字段的信息不可用,视具体字段用’0‘或’*‘进行表示, 11个字段从左到右,依次为:

1、QNAME,测序序列(Query)的名称;

2、FLAG,位标识,用数字表示Query的比对情况,每一个数字代表一种比对情况,这里的值是所有符合情况的数字相加。数字的具体含义,参考如下所示的FLAG说明;

l   |0x0001 (1)|         the read is paired in sequencing

l   |0x0002 (2)|         the read is mapped in a proper pair

l   |0x0004 (4)|         the query sequence itself is unmapped

l   |0x0008 (8)|         the mate is unmapped

l   |0x0010 (16)|       query being reverse complemented

l   |0x0020 (32)|       mate being reverse complemented

l   |0x0040 (64)|       the read is the first read in a pair

l   |0x0080 (128)|     the read is the second read in a pair

l   |0x0100 (256)|     the alignment is not primary

l   |0x0200 (512)|     the read fails quality controls

l   |0x0400 (1024)|   the read is either a PCR or an optical duplicate

3、RNAME,比对上的参考序列名称。对于没有比对上的结果,用’*‘表示;

4、POS,比对上的参考序列位置坐标,从1开始计数(1-based leftmost position)。对于没有比对上的结果,用’0‘表示;

5、MAPQ,比对的质量值;

6、CIGAR,简要比对信息表达式(CIGAR,Compact Idiosyncratic Gapped Alignment Report),其以参考序列为基础,使用数字加字母表示比对结果。比如,3S6M1P1I4M,前三个碱基被剪切去除了,然后6个比对上了,然后打开了一个缺口,有一个碱基插入,最后4个比对上了。对于没有比对上的结果,用’*‘表示。每个字母的具体含义,参考如下所示的说明;

l   M      alignment match (can be a sequence match or mismatch)

l   I        insertion to the reference

l   D       deletion from the reference

l   N       skipped region from the reference

l   S        soft clipping (clipped sequences present in SEQ)

l   H       hard clipping (clipped sequences NOT present in SEQ)

l   P        padding (silent deletion from padded reference)

l   =        sequence match

l   X       sequence mismatch

7、RNEXT,PE测序另一端read比对上的参考序列名称。如果信息不可用,用’*‘表示。如果比对上的序列与RNAME相同,用’=‘表示;

8、PNEXT,PE测序另一端read比对上的参考序列位置坐标,如果信息不可用,用’0‘表示;

9、TLEN,PE测序数据,根据比对结果推断的插入片段大小,如果是SE测序或者信息不可用,用’0‘表示;

10、SEQ,Query对应的序列信息,Query序列与参考序列具有同样的方向;

11、QUAL,Query序列的测序质量值信息。

可选字段(optional fields),格式如:TAG:TYPE:VALUE,每个TAG代表一类信息, TYPE表示TAG对应值的类型,可以是字符串、整数、字节、数组等,一些TAG及其含义如下所示:

l   NM    Edit distance

l   MD    Mismatching positions/bases

l   AS     Alignment score

l   BC     Barcode sequence

l   X0     Number of best hits

l   X1     Number of suboptimal hits found by BWA

l   XN    Number of ambiguous bases in the referenece

l   XM    Number of mismatches in the alignment

l   XO    Number of gap opens

l   XG    Number of gap extentions

l   XT     Type: Unique/Repeat/N/Mate-sw

l   XA    Alternative hits; format: (chr,pos,CIGAR,NM;)*

l   XS     Suboptimal alignment score

l   XF     Support from forward/reverse alignment

l   XE     Number of supporting seeds

 注:X开头的TAGBWA特有的TAG


绝大多数情况下,使用BWA的默认参数进行比对即可。如需进行个性的调整,更改对应的参数即可。


分析模块引用了BWA v0.7.12软件(http://bio-bwa.sourceforge.net/)。


相关文献如下所示:

Li H. and Durbin R. (2009) Fast and accurate short read alignment with Burrows-Wheeler Transform. Bioinformatics, 25:1754-60. [PMID: 19451168]

Li H. and Durbin R. (2010) Fast and accurate long-read alignment with Burrows-Wheeler Transform. Bioinformatics, Epub. [PMID: 20080505]