【4.1.3】序列比对到参考基因组--bwa

Mapping 阅读量：次

一、简介

BWA，即Burrows-Wheeler-Alignment Tool。BWA 是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。它由三个不同的算法：

  • BWA-backtrack: 是用来比对 Illumina 的序列的，reads 长度最长能到 100bp。-
  • BWA-SW: 用于比对 long-read ，支持的长度为 70bp-1Mbp；同时支持剪接性比对。
  • BWA-MEM: 推荐使用的算法，支持较长的read长度，同时支持剪接性比对（split alignments)，但是BWA-MEM是更新的算法，也更快，更准确，且 BWA-MEM 对于 70bp-100bp 的 Illumina 数据来说，效果也更好些。 对于上述三种算法，首先需要使用索引命令构建参考基因组的索引，用于后面的比对。所以，使用BWA整个比对过程主要分为两步，
  1. 第一步建索引，
  2. 第二步使用BWA MEM进行比对。

bwa的使用需要两中输入文件：

Reference genome data（fasta格式 .fa, .fasta, .fna）
Short reads data (fastaq格式 .fastaq, .fq)

二、BWA的安装

a.下载BWA http://bio-bwa.sourceforge.net/bwa.shtml

#解压
$ tar -jxvf bwa-*.tar.bz2
$ cd bwa-*；

# 编译BWA
$ make

$ echo 'PATH=$PATH:/path/bwa--*' >> ~/.bashrc
$ source ~/.bashrc

或

git clone https://github.com/lh3/bwa.git
cd bwa
make

安装好之后，会出现一个名为bwa的可执行文件，输入下面命令可以查看帮助信息

./bwa

Program: bwa (alignment via Burrows-Wheeler transformation)
Version: 0.7.17-r1188
Contact: Heng Li <<a href="mailto:lh3@sanger.ac.uk" _href="mailto:lh3@sanger.ac.uk">lh3@sanger.ac.uk</a>>
Usage:   bwa <command> [options]
Command: index         index sequences in the FASTA format
         mem           BWA-MEM algorithm
         fastmap       identify super-maximal exact matches
         pemerge       merge overlapping paired ends (EXPERIMENTAL)
         aln           gapped/ungapped alignment
         samse         generate alignment (single ended)
         sampe         generate alignment (paired ended)
         bwasw         BWA-SW for long queries
         shm           manage indices in shared memory
         fa2pac        convert FASTA to PAC format
         pac2bwt       generate BWT from PAC
         pac2bwtgen    alternative algorithm for generating BWT
         bwtupdate     update .bwt to the new format
         bwt2sa        generate SA from BWT and Occ
Note: To use BWA, you need to first index the genome with `bwa index'.
      There are three alignment algorithms in BWA: `mem', `bwasw', and
      `aln/samse/sampe'. If you are not sure which to use, try `bwa mem'
      first. Please `man ./bwa.1' for the manual.

三、 BWA 的使用

bwa软件的作用是将序列比对到参考基因组上，在比对之前，首先需要对参考基因组建立索引。

建立索引 index

在进行 reads 的比对前，需要对 fasta 文件构建 FM-index。

用法和参数如下：

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

bwa index 指令更多的用法及 options，通过bwa index 命令来查看

# 根据reference genome data(e.g. ref.fa) 建立 Index File：
$ bwa index ref.fa -p genome        # 可以不加-p genome，这样建立索引都是以ref.fa为前缀

构建出参考基因组的 FM-index，建立好参考基因组之后，就可以进行比对了。不同的比对算法，命令不同。

aln/samse/sampe ----> BWA-backtrack (samse 中的 se 是 single-end 的简写，而 sampe 中的 pe 是 paired-end 的简写）。

bwasw ----> BWA-SW
mem ----> BWA-MEM

3.1 BWA-MEM 算法

该算法先使用 MEM(maximal exact matches) 进行 seeding alignments，再使用 SW(affine-gap Smith-Waterman) 算法进行 seeds 的延伸。BWA–MEM 算法执行局部比对和剪接性。可能会出现 query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点。有些软件如 Picard’s markDuplicates 跟 mem 的这种剪接性比对不兼容,在这种情况下，可以使用 –M 选项来将 shorter split hits 标记为次优。

对应的子命令为mem, 基本用法如下

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

Options:

　　-t INT 线程数，默认是1，增加线程数，会减少运行时间。
　　-M 将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。
　　-p 若无此参数：输入文件只有1个，则进行单端比对；若输入文件有2个，则作为paired 
reads进行比对。若加入此参数：则仅以第1个文件作为输入(会忽略第二个输入序列文件，把第一个文件当做单端测序的数据进行比对)，该文件必须是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，但是这些比对的结果会被标记为次优。
　　-Y 对数据进行soft clipping, 当错配或者gap数过多比对不上时，会对序列进行切除，这里的切除并只是在比对时去掉这部分序列，最终输出结果中序列还是存在的，所以称为soft clipping。

特别说明：

  • 如果 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，也会被忽略。

示例：

# SE （single end）
$ bwa mem ref.fa reads.fq > mem-se.sam

# PE（paired end）
$ 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

对于超长读长的reads,比如PacBio和Nanopore测序仪产生的reads, 用法如下

$ bwa mem -x pacbio ref.fa reads.fq > aln.sam
$ bwa mem -x ont2d ref.fa reads.fq > aln.sam

上述的代码中, ref.fa指的是参考基因组索引的名字，对于上面提供的小鼠的例子而言，参考基因组索引的名字为mm10.fasta, 注意是不包含后缀的。

3.2 BWA-backtrack 算法

对应的子命令为aln/samse/sample

Usage: bwa aln [options] ref.fa read.fq > aln_sa.sai 
Options:
　　-o int：允许出现的最大gap数。
　　-e int：每个gap允许的最大长度。
　　-d int：不允许在3’端出现大于多少bp的deletion。
　　-i int：不允许在reads两端出现大于多少bp的indel。
　　-l int：Read前多少个碱基作为seed，如果设置的seed大于read长度，将无法继续，最好设置在25-35，与-k 2 配合使用。
　　-k int：在seed中的最大编辑距离，使用默认2，与-l配合使用。
　　-t int：要使用的线程数。
　　-R int：此参数只应用于pair end中，当没有出现大于此值的最佳比对结果时，将会降低标准再次进行比对。增加这个值可以提高配对比对的准确率，但是同时会消耗更长的时间，默认是32。
　　-I int：表示输入的文件格式为Illumina 1.3+数据格式。
　　-B int：设置标记序列。从5’端开始多少个碱基作为标记序列，当-B为正值时，在比对之前会将每个read的标记序列剪切，并将此标记序列表示在BC SAM 标签里，对于pair end数据，两端的标记序列会被连接。
　　-b ：指定输入格式为bam格式。

用法如下：

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

# 对于single-read
$ bwa aln [options] ref.fa read.fq > aln_sa.sai                　#寻找 SA coordinates
$ bwa samse [options] ref.fa aln_sa.sai read.fq > aln-se.sam    # 转换SA coordinates输出为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

如果希望多线程运行，在其中加入 -t这个参数，另外-f这个参数可以指定结果输出文件，如:

$ bwa aln -c -t 3 -f aln_sa1.sai ref.fa  read1.fq

３.3 BWA-SW 算法

对应的子命令为bwasw, 基本用法如下

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

例子：

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

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

利用bwa进行比对的例子：

1）对参考基因组（reference）构建索引

echo "bwa index starts@"`date` && \
cd ref && \
bwa index -a bwtsw genome.fa && \
echo "bwa index ends@"`date`

生成的文件有以下几种类型：

genome.fasta.ann
genome.fasta.pac
genome.fasta.amb
genome.fasta.bwt
genome.fasta.sa

2）用bwa mem 进行比对

echo "bwa mem starts@"`date` && \
cd ref && \
bwa mem -t 4 -M -R "\@RG\tID:{library}\tLB:{library}\tPL:Illumina\tPU:{sample}\tSM:{sample}\" \
                /home/data/genome.fasta \
                R_1.fastq R_2.fastq \
                > bwa.mem.sam \
                2> ./bwa.log&& \
echo "bwa mem ends@"`date`
fai：是对ref基因组文件建的索引，方便软件快速随机读取基因组序列
sai：是将fastq比对后出来的文件，用于最后输出比对结果sam文件的

四、我的案例

参考资料

药企，独角兽，苏州。团队长期招人，感兴趣的都可以发邮件聊聊：tiehan@sina.cn
个人公众号，比较懒，很少更新，可以在上面提问题，如果回复不及时，可发邮件给我： tiehan@sina.cn