序列比对BWA（Docker篇）

BWA，即Burrows-Wheeler-Alignment Tool。BWA 是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包。首先通过BWT(Burrows-Wheeler Transformation，BWT压缩算法)为大型参考基因组建立索引，然后将reads比对到基因组。特点是快速、准确、省内存。它有三个不同的算法：

BWA-backtrack: 是用来比对 Illumina 的序列的，reads 长度最长能到 100bp。

BWA-SW: 用于比对 long-read ，支持的长度为 70bp-1Mbp；同时支持剪接性比对。

BWA-MEM: 推荐使用的算法，支持较长的read长度，同时支持剪接性比对（split alignments)，但是BWA-MEM是更新的算法，也更快，更准确，且 BWA-MEM 对于 70bp-100bp 的 Illumina 数据来说，效果也更好些。

对于不熟悉Linux的人来说，安装BWA并不太容易，那有没有方法可以省去繁琐的安装步骤而直接使用BWA呢？答案是有的！Docker工具可以完美的决这个问题！

下载镜像并运行

安装好Docker后，搜索我们omicsclass提供的docker镜像，其中便有安装好BWA软件的镜像（下图均可左右拖动）。

$ docker search omicsclass

NAME DESCRIPTION STARS OFFICIAL AUTOMATED

omicsclass/gene-family gene-family analysis docker image 2

omicsclass/rnaseq RNA-seq analysis docker image build by omics… 1

omicsclass/bwa BWA v0.7.17 build by omicsclass 0

omicsclass/samtools samtools v1.10 build by omicsclass 0

omicsclass/isoseq3 isoseq3 v3.3.0 build by omicsclass 0

omicsclass/biocontainer-base Biocontainers base Image centos7 0

omicsclass/blast-plus blast+ v2.9.0 0

omicsclass/blastall legacy blastall v2.2.26 0

omicsclass/sratoolkit SRAtoolkit v2.10.3 and aspera v3.9.9.177872 0

下载镜像。

$ docker pull omicsclass/bwa

Using default tag: latest

latest: Pulling from omicsclass/bwa

ab5ef0e58194: Already exists

417469905807: Already exists

ed09842cc19f: Already exists

f860268ff83f: Already exists

f87dd41136a6: Already exists

90091b4f5d91: Already exists

4d81fa76e64d: Pull complete

Digest: sha256:b2afe73fa4abbd4527d349ecd94ccafb9e06f719d5c12b3335b6bc1c8c397a3d

Status: Downloaded newer imageforomicsclass/bwa:latest

下载完成后可以检查一下。

$ docker images

REPOSITORY TAG IMAGE ID CREATED SIZE

omicsclass/samtools latest 9373e18781bf 3 weeks ago 2.04GB

omicsclass/bwa latest b39d68bfa8e5 3 weeks ago 2.1GB

omicsclass/blast-plus latest 0220cac51a6e 3 weeks ago 2.55GB

进入虚拟机。

$ docker run --rm -v /d/bwa:/work -it omicsclass/bwa:latest##这里我使用的Docker是windows Toolbox版本

######################################################

# 欢迎使用组学大讲堂提供的docker镜像 #

# 问题交流请访问：www.omicsclass.com #

######################################################

Linux新手建议学习课程：

--> https://www.omicsclass.com/article/702

搭建实验室生信分析平台与docker使用详情见课程：

--> https://www.omicsclass.com/article/1181

BWA使用（BWA-MEM算法）

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

用法和参数如下：

Usage：bwa index [ –p prefix ] [ –a algoType ]

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，建立好参考基因组之后，就可以进行比对了。对应的子命令为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。

$ 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

对于超长读长的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指的是参考基因组索引的名字。

更多生物信息课程：

1. 文章越来越难发？是你没发现新思路，基因家族分析发2-4分文章简单快速，学习链接：基因家族分析实操课程、基因家族文献思路解读

2. 转录组数据理解不深入？图表看不懂？点击链接学习深入解读数据结果文件，学习链接：转录组（有参）结果解读；转录组（无参）结果解读

3. 转录组数据深入挖掘技能-WGCNA，提升你的文章档次，学习链接：WGCNA-加权基因共表达网络分析

4. 转录组数据怎么挖掘？学习链接：转录组标准分析后的数据挖掘、转录组文献解读

5.微生物16S/ITS/18S分析原理及结果解读、OTU网络图绘制、cytoscape与网络图绘制课程

6. 生物信息入门到精通必修基础课，学习链接：linux系统使用、perl入门到精通、perl语言高级、R语言画图

7. 医学相关数据挖掘课程，不用做实验也能发文章，学习链接：TCGA-差异基因分析、GEO芯片数据挖掘、GSEA富集分析课程、TCGA临床数据生存分析、TCGA-转录因子分析、TCGA-ceRNA调控网络分析

8.其他课程链接：二代测序转录组数据自主分析、NCBI数据上传、二代测序数据解读。