01 背景
Burrow-Wheeler Aligner for short-read alignment (see minimap2 for long-read alignment)
1.1 介绍
BWA是一个用于将DNA序列比对到大型参考基因组(如人类基因组)的软件包。它包含三种算法:BWA-backtrack、BWA-SW和BWA-MEM。第一个算法设计用于 Illumina 测序读取长度最长为100bp,而后两个用于更长序列,范围从70bp到几百万碱基。BWA-MEM和BWA-SW具有类似的特性,如支持长读取和嵌合比对,但通常建议使用BWA-MEM,因为它更快、更准确。对于70-100bp的Illumina读取,BWA-MEM的性能也比BWA-backtrack更好。
对于所有算法,BWA首先需要为参考基因组构建FM索引(使用index命令)。比对算法通过不同的子命令调用:BWA-backtrack使用aln/samse/sampe,BWA-SW使用bwasw,BWA-MEM使用mem。
1.2 可用性
BWA以GPLv3协议发布。最新的源代码可以在 GitHub 免费获取。发布的软件包可以在 SourceForge 下载。获得源代码后,只需使用make进行编译,并将单个可执行文件bwa复制到您想要的目的地。构建BWA所需的唯一依赖是zlib。
自0.7.11版以来,x86_64-linux的预编译二进制文件已经包含在bwakit中。除了BWA之外,这个自洽的软件包还配备了BWA相关的工具和第三方工具,用于正确的BAM到FASTQ转换,比对到ALT染色体,适配器修剪,重复标记,HLA类型鉴定和相关的数据文件。
BWA是一种能够将差异度较小的序列比对到一个较大的参考基因组上的软件包,包括BWA-backtrack, BWA-SW和BWA-MEM三种算法。
1.3 注意
minimap2已经取代了BWA-MEM用于PacBio和Nanopore读取的比对。它保留了所有主要的BWA-MEM特性,但速度大约快了50倍,更加灵活、准确,并且产生更好的碱基级别比对。BWA-MEM2的beta版本已发布,用于短读比对。BWA-MEM2比BWA-MEM快大约两倍,并且输出几乎相同的比对结果。
02 参考
https://github.com/lh3/bwa #官网
03 安装
git clone https://github.com/lh3/bwa.git
cd bwa; make
./bwa index ref.fa
./bwa mem ref.fa read-se.fq.gz | gzip -3 > aln-se.sam.gz
./bwa mem ref.fa read1.fq read2.fq | gzip -3 > aln-pe.sam.gz
04 使用
程序:bwa(通过Burrows-Wheeler变换进行比对)
版本:0.7.17-r1188
联系人:Heng Li <lh3@sanger.ac.uk>
用法: bwa <命令> [选项]
命令: index 对FASTA格式的序列建立索引
mem BWA-MEM算法
fastmap 识别超极大精确匹配
pemerge 合并重叠的配对末端(实验性功能)
aln 带间隙/不带间隙的比对
samse 生成比对结果(单端)
sampe 生成比对结果(配对端)
bwasw 用于长查询的BWA-SW
shm 管理共享内存中的索引
fa2pac 将FASTA格式转换为PAC格式
pac2bwt 从PAC生成BWT
pac2bwtgen 生成BWT的备用算法
bwtupdate 更新.bwt到新格式
bwt2sa 从BWT和Occ生成SA
注意:使用BWA之前,您需要首先使用 `bwa index' 对基因组建立索引。
BWA中有三种比对算法:`mem'、`bwasw' 和 `aln/samse/sampe'。如果您不确定使用哪种,请先尝试 `bwa mem'。
请查阅 `man ./bwa.1' 获取详细的使用说明。
bwa index
用法:bwa index [选项] <in.fasta>
选项: -a STR BWT 构建算法:bwtsw、is 或 rb2 [自动]
-p STR 索引的前缀 [与 fasta 文件名相同]
-b INT bwtsw 算法的块大小(与 -a bwtsw 一起有效)[10000000]
-6 索引文件命名为 <in.fasta>.64.* 而不是 <in.fasta>.*
警告: `-a bwtsw' 不适用于短基因组,而 `-a is' 和 `-a div' 不适用于长基因组。
本质上,短读比对软件
bwa mem
用法:bwa mem [选项] <idxbase> <in1.fq> [in2.fq]
算法选项:
-t INT 线程数 [1]
-k INT 最小种子长度 [19]
-w INT 用于带状比对的带宽 [100]
-d INT 斜对角线 X-dropoff [100]
-r FLOAT 在长于 {-k} * FLOAT 的种子内查找内部种子 [1.5]
-y INT 第3轮种子扩展的种子出现次数 [20]
-c INT 跳过出现次数超过 INT 的种子 [500]
-D FLOAT 删除短于最长重叠链的 FLOAT 分数的链 [0.50]
-W INT 如果种子基数短于 INT,则丢弃链 [0]
-m INT 为每个读取执行最多 INT 轮配对恢复 [50]
-S 跳过配对恢复
-P 跳过配对;除非也使用了 -S
评分选项:
-A INT 序列匹配的分数,该分数调节了选项 -TdBOELU,除非被覆盖 [1]
-B INT 不匹配的惩罚 [4]
-O INT[,INT] 插入和删除的间隙开放惩罚 [6,6]
-E INT[,INT] 间隙扩展惩罚;大小为 k 的间隙成本为 '{-O} + {-E}*k' [1,1]
-L INT[,INT] 5' 和 3' 端修剪的惩罚 [5,5]
-U INT 未配对读取对的惩罚 [17]
-x STR 读取类型。设置 -x 会改变多个参数,除非被覆盖 [null]
pacbio: -k17 -W40 -r10 -A1 -B1 -O1 -E1 -L0 (PacBio 读取到参考序列)
ont2d: -k14 -W20 -r10 -A1 -B1 -O1 -E1 -L0 (Oxford Nanopore 2D 读取到参考序列)
intractg: -B9 -O16 -L5 (种内物种连片到参考序列)
输入/输出选项:
-p 智能配对(忽略 in2.fq)
-R STR 读取分组头行,例如 '@RG\tID:foo\tSM:bar' [null]
-H STR/FILE 如果以 @ 开头,则将 STR 插入到头部;或者插入 FILE 中的行 [null]
-o FILE 输出比对结果的 SAM 文件 [stdout]
-j 将 ALT 染色体视为主要组装的一部分(即忽略 <idxbase>.alt 文件)
-5 对于分割比对,将具有最小坐标的比对视为主要比对
-q 不修改补充比对的 mapQ
-K INT 每批处理 INT 个输入碱基,无论 nThreads 如何(用于可重现性) []
-v INT 详细级别:1=错误,2=警告,3=消息,4+=调试 [3]
-T INT 输出的最低分数 [30]
-h INT[,INT] 如果有 <INT 个比对的分数 >80% 的最大分数,则在 XA 中输出所有 [5,200]
-a 对于单端或未配对的 PE,输出所有比对
-C 将 FASTA/FASTQ 注释附加到 SAM 输出
-V 在 XR 标签中输出参考 FASTA 头部
-Y 对于补充比对使用软剪辑
-M 将较短的分割比对标记为次要比对
-I FLOAT[,FLOAT[,INT[,INT]]]
指定插入大小分布的均值、标准差(如果不存在则为均值的10%),最大值(如果不存在则为均值的4个标准差)和最小值。
仅限 FR 方向。[推断]
注意:请阅读 man 页面以获取命令行和选项的详细描述。
05 常用命令行
bwa index cankao.fasta #默认算法即可
bwa index -a bwtsw draft.asm.fasta #不太合适
bwa mem -SP5M -t 24 cankao.fasta R1.fastq.gz R2.fastq.gz ###快速比对
bwa aln -t 24 cankao.fasta R1.fastq.gz > R1.sai
bwa aln -t 24 cankao.fasta R2.fastq.gz > R2.sai
#一直只能跑一个,内存容易炸
bwa sampe cankao.fasta R1.sai R2.sai R1.fastq.gz R2.fastq.gz > R1_R2.bwa_aln.sam
06 参考文献
-
Li H. and Durbin R. (2009) Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics, 25, 1754-1760. [PMID: 19451168]. (if you use the BWA-backtrack algorithm)
-
Li H. and Durbin R. (2010) Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics, 26, 589-595. [PMID: 20080505]. (if you use the BWA-SW algorithm)
-
Li H. (2013) Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv:1303.3997v2 [q-bio.GN]. (if you use the BWA-MEM algorithm or the fastmap command, or want to cite the whole BWA package)
07 Q&A
7.1 BWA可以处理哪些类型的数据?
BWA可以处理各种类型的DNA序列数据,尽管最佳算法和设置可能会有所不同。以下是推荐的设置:
Illumina/454/IonTorrent的单端读取长度超过~70bp或组装的contigs(成分)长达数百万碱基对,映射到与参考基因组密切相关的:
bwa mem ref.fa reads.fq > aln.sam
Illumina单端读取长度短于~70bp:
bwa aln ref.fa reads.fq > reads.sai; bwa samse ref.fa reads.sai reads.fq > aln-se.sam
Illumina/454/IonTorrent的双端读取长度超过~70bp:
bwa mem ref.fa read1.fq read2.fq > aln-pe.sam
Illumina双端读取长度短于~70bp:
bwa aln ref.fa read1.fq > read1.sai; bwa aln ref.fa read2.fq > read2.sai
bwa sampe ref.fa read1.sai read2.sai read1.fq read2.fq > aln-pe.sam
PacBio子读取或Oxford Nanopore读取到参考基因组:
bwa mem -x pacbio ref.fa reads.fq > aln.sam
bwa mem -x ont2d ref.fa reads.fq > aln.sam
实际上区别不大,对于150bp而言,长度还是用minimap2
对于长度超过~70bp的查询序列,建议使用BWA-MEM,对于各种错误率(或序列差异)都适用。一般来说,BWA-MEM对于长查询序列更容忍错误,因为错过所有seed的几率较小。正如上文所示,使用非默认设置,BWA-MEM可以处理具有超过20%测序错误率的Oxford Nanopore读取。
7.2 为什么一个读取在输出SAM中出现多次?
BWA-SW和BWA-MEM执行局部比对。如果存在染色体易位、基因融合或长缺失,一个跨越断点的读取可能具有两个命中点,在SAM输出中占据两行。使用BWA-MEM的默认设置,只有一行是主要的并且是软剪辑的;其他行带有0x800的SAM标志(辅助比对),并且是硬剪辑的。
7.3 BWA是否适用于总长度超过4GB的参考序列?
是的。自0.6.x版本以来,所有BWA算法都适用于总长度超过4GB的基因组。然而,单个染色体的长度不应超过2GB。
7.4 为什么一个读取对的映射质量很高,但另一个读取的映射质量为零?
这是正确的。映射质量是针对单个读取分配的,而不是针对读取对。有可能一个读取可以被明确地映射,但它的配对读取落在串联重复中,因此无法确定其准确的位置。
7.5 BWA-backtrack比对如何突出染色体的末端?
在内部,BWA将所有参考序列连接成一个长序列。一个读取可能被映射到两个相邻参考序列的连接点。在这种情况下,BWA-backtrack将标记读取为未映射的(0x4),但您将看到位置、CIGAR和所有标签。BWA-SW比对也可能出现类似的问题。BWA-MEM没有这个问题。
7.6 BWA是否适用于GRCh38发布中的ALT contigs?
是的,自0.7.11版本以来,BWA-MEM正式支持对GRCh38+ALT的比对。截至目前,BWA-backtrack和BWA-SW不正确地支持ALT比对。请参阅README-alt.md获取详细信息。简言之,建议使用bwakit生成参考基因组和进行比对。
7.7 我可以直接对GRCh38+ALT运行BWA-MEM而不进行后处理吗?
如果您对ALT contigs的命中不感兴趣,可以在不进行后处理的情况下运行BWA-MEM。这种方式产生的比对结果与不包含ALT contigs的GRCh38的比对结果非常接近。尽管如此,应用后处理有助于减少由于来自ALT contigs分歧部分的读取导致的虚假比对,并且还可以实现HLA分型。建议运行后处理脚本。
比对就是慢,但是准,patient!