Fork me on GitHub

BWA 命令详解

发表于   |  阅读量 次   |   分类于   |  

一、BWA 简介

BWA，即Burrows-Wheeler-Alignment Tool。
BWA 的学习主要来自参考网站和使用手册；BWA 是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。它由三个不同的算法：

1
2
3
BWA-backtrack: 是用来比对 Illumina 的序列的，reads 长度最长能到 100bp。
BWA-SW: 用于比对 long-read ，支持的长度为 70bp-1Mbp；同时支持剪接性比对。
BWA-MEM: 都支持较长的read长度，同时都支持剪接性比对（split alignments)，但是BWA-MEM是更新的算法，也更快，更准确，且 BWA-MEM 对于 70bp-100bp 的 Illumina 数据来说，效果也更好些。

二、 BWA 下载和安装

下载地址：sourceforge download

1
2
3
4
5
6
$ wget http://superb-sea2.dl.sourceforge.net/project/bio-bwa/bwa-0.7.12.tar.bz2
$ tar jxf bwa-0.7.12.tar.bz2
$ cd bwa-0.7.12
$ make
$ echo 'PATH=$PATH:/your/bwa/path' >> ~/.bashrc
$ source ~/.bashrc

三、 BWA 使用

在运用这三种算法之前，需要利用 BWA 的 index 命令，构建出参考基因组的 FM-index，而对与上述的三种不同的算法而言，又有不同的命令：

1
2
3
aln/samse/sampe ----> BWA-backtrack   (samse 中的 se 是 single-end 的简写，而 sampe 中的 pe 是 paired-end 的简写）。
bwasw ----> BWA-SW
mem ----> BWA-MEM

3.1 建立索引 index

在进行 reads 的比对前，需要对 fasta 文件构建 FM-index。
用法和参数如下：

1
2
3
4
5
6
7
index   Usage：bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>
        Index database sequence in the FASTA format.
        OPTIONS:
        -P STR  输出数据库的前缀；【默认和输入的文件名一致，输出的数据库在其输入文件所在的文件夹，并以该文件名为前缀。】
        -a [is|bwtsw]   构建index的算法，有两个算法：
                        is  是默认的算法，虽然相对较快，但是需要较大的内存，当构建的数据库大于2GB的时候就不能正常工作了。
                        bwtsw   对于短的参考序列式不工作的，必须要大于等于10MB, 但能用于较大的基因组数据，比如人的全基因组。

根据reference genome data(e.g. ref.fa) 建立 Index File例子：

1
2
$ bwa index ref.fa -p genome
###可以不加-p genome，这样建立索引都是以ref.fa为前缀,以下使用的示例均未使用-p genome###

3.2 mem比对

该算法先使用 MEM(maximal exact matches) 进行 seeding alignments，再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。
BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点。这对 long reads 的比对时比较重要的结果。但是却会和 Picard 的 markDuplicates 程序部兼容。
使用方法如下：

1
mem Usage: bwa mem [options] ref.fa reads.fq [mates.fq]

简单的来讲，mem 使用的 MEMs(maximal exact matches） 进行seedling alignments, 再使用 SW(affine-gap Smith-Waterman 算法）进行seeding extending.
mem 进行局部比对，因此，对于一条序列的不同区域可能会产生多种最优匹配结果， 这对于long reads 来说尤为重要。 有些软件如 Picard’s markDuplicates 跟 mem 的这种剪接性比对不兼容,
在这种情况下，可以使用 –M 选项来将 shorter split hits 标记为次优。
特别说明：
如果 mates.fq 缺省，且参数 –p 未设定，那么 reads.fq 被认为是 single-end;
如果 mates.fq 存在，且参数 –p 未设定，那么 mem 命令会认为 read.fq 和 mates.fq 中的 i-th reads 组成一个read对 (a read pair)，这个模式是常用的 paired-end mode。
如果参数 –p 被设定，那么， mem 命令会认为 read.fq 中的 第 2i-th 和 第 (2i + 1)-th 的 reads 组成一个 read 对 （a read pair），这种方式也被成为交错式的（interleaved paired-end)。 在这种情况下，即使有 mates.fq，也会被忽略。
常用的参数如下：

1
2
3
4
5
6
-t INT  线程数，默认是1。
-M      将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
-p      若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired reads进行比对。若加入此参数：则仅以第1个文件作为输入(输入的文件若有2个，则忽略之)，该文件必须是read1.fq和read2.fa进行reads交叉的数据。
-R STR  完整的read group的头部，可以用 '\t' 作为分隔符， 在输出的SAM文件中被解释为制表符TAB. read group 的ID，会被添加到输出文件的每一个read的头部。
-T INT  当比对的分值比 INT 小时，不输出该比对结果，这个参数只影响输出的结果，不影响比对的过程。
-a      将所有的比对结果都输出，包括 single-end 和 unpaired paired-end的 reads，但是这些比对的结果会被标记为次优。

例子：

1
2
3
$ bwa mem ref.fa reads.fq > mem-se.sam
$ bwa mem ref.fa read1.fq read2.fq > mem-pe.sam
$ bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\"  ref.fa read1.fastq read2.fastq > mem-pe.sam 2> ./mem-pe.log

3.3 align/samse/sampe

用法如下：

1
2
3
4
5
6
7
对于single-read：
bwa aln [options] ref.fa read.fq > aln_sa.sai
bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam
对于pair-reads：
bwa aln [options] ref.fa read1.fq > aln_sa1.sai
bwa aln [options] ref.fa read2.fq > aln_sa2.sai
bwa sampe [options] ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam

3.4 bwasw

对输入的第1个文件的所有序列进行比对。如果输如有 2 个文件，则进行 paired-end 比对，此模式仅对 Illumina 的 short-insert 数据进行比对。在 Paired-end 模式下，BWA-SW依然输出剪接性比对结果，但是这些结果会标记为 not properly paired; 同时如果有多个匹配位点，则不会写入 mate 的匹配位置。
常用参数如下：

1
2
3
bwasw   
Usage: bwa bwasw [ options ] ref.fasta reads.fq [mate.fq] > aln.sam
-t  INT 使用的线程数

例子：

1
2
$ bwa bwasw genome long_read.fq > aln.sam
$ bwa bwasw genome read1.fq read2.fq > aln-pe.sam

3.5 backtrack

经典的 bwa 先使用 aln 命令将单独的 reads 比对到参考序列，再使用 samse 或 sampe 生成 sam 文件。
常用例子：

1
2
3
4
$ bwa aln genome read1.fq > aln_sa1.sai
$ bwa aln genome read2.fq > aln_sa2.sai
$ bwa samse genome aln_sa1.sai read1.fq > aln_se.sam
$ bwa sampe genome aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam

热评文章