简介
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整个比对过程主要分为两步,第一步建索引,第二步使用BWA MEM进行比对。
bwa安装
直接使用mamba安装即可
mamba install bwa
建立参考基因组索引
bwa index [ –p prefix ] [ –a algoType ] <in.db.fasta>
用法说明
-
-p
STR 输出数据库的前缀;【默认和输入的文件名一致,输出的数据库在其输入文件所在的文件夹,并以该文件名为前缀。】 -
-a [algoType]
构建index的算法,有以下两个选项:-
-a is
是默认的算法,虽然相对较快,但是需要较大的内存,当构建的数据库大于2GB的时候就不能正常工作了; -
-a bwtsw
对于短的参考序列式不工作的,必须要大于等于10MB, 但能用于较大的基因组数据,比如人的全基因组。
-
示例代码
bwa index GRCh38.primary_assembly.genome.fa
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, 基本用法如下
bwa mem [options] ref.fa reads.fq [mates.fq]
参数说明
-
-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,也会被忽略。
示例代码
# single end
bwa mem ref.fa reads.fq > mem-se.sam
# 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
参考
https://www.jianshu.com/p/19f58a07e6f4
https://blog.csdn.net/weixin_42192188/article/details/132286609