BWA是一种能够将低差异度的序列比对到较大的参考基因组上的软件。常常被用来处理质控后的clean reads，通过reads与参考序列比对，生成储存变异信息的sam或bam文件。目前BWA存在3种比对算法，可支持长度在1Mb以内的reads 与参考基因组的比对。
bwa的3种算法：
BWA-backtrack： 专门针对100bp以下的illumina测序数据 (aln/samse/sampe)
BWA-SW：70bp-1Mbp、支持可变剪切(bwasw)
BWA-MEM：70bp-1Mbp、支持可变剪切(mem)
任何一种算法，BWA都需要首先对"参考基因组"建立FM-index (bwa-index)
1. 建立参考基因组的索引
无论选择哪一种比对算法，都需要对待使用的参考基因组建立索引。
bwa index的用法和参数：
bwa index [-p prefix] [-a algoType] <in.db.fasta>
-p #输出数据库的前缀，默认与输入文件名一致
-a #设置参考基因组构建index的算法（有is 和 bwtsw两种）
简要介绍bwa构建参考基因组index的两种算法：
-is用于构建后缀数组的线性时间算法：缺点--比较耗内存（需要参考基因组5.37倍的内存）；优点--速度快。is是构建index的默认算法，但是不适合 “参考基因组” 大于2GB的文件。
-bwtsw使用BWT-SW算法，多用于人类基因组。
2. 三种不同比对算法的使用
2.1 BWA-MEM
在reads中划分seed区域，使用seed序列与建立过索引的参考基因组比对，选择 maximal exact matches(MEMs) 的位置作为reads在参考基因组上的正确比对区域（seed alignment），然后在此基础上通过affine-gap Smith-Waterman算法在比对区域附近延伸（seed extension）。
bwa mem [options] ref.fa reads.fq [mates.fq]
-p pair-end模式参数 ##适应于不同的应用场景：（1）该参数不存在，fastq文件只有reads.fq：reads.fq会被默认为single-end数据；（2）该参数不存在，fastq文件有reads.fq 和mates.fq，软件默认reads.fq 和mates.fq为paired-end mode数据；（3）该参数存在，fastq文件只有reads.fq：read.fq 中的 第 2i-th 和 第 (2i + 1)-th 的 reads 组成一个 read 对 ；（4）该参数存在，astq文件有reads.fq 和mates.fq，软件默认忽略掉mates.fq，认为read.fq 中的 第 2i-th 和 第 (2i + 1)-th 的 reads 组成一个 read 对 。
-M 兼容Picard’s markDuplicates的参数 ##当bam文件中存在一条read比对到参考基因组上不同位置时，Picard的markDuplicates工具不能正常运行，但由于bwa-mem采用的是局部比对的算法，不可避免得会出现， query 序列的多个不同的部位出现各自的最优匹配，导致 reads 有多个最佳匹配位点，因此使用-M参数将shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 插件。
-t INT 线程数（默认为1） ##线程数目越大，运行速度越快
-k INT minSeedLen(最小seed长度,默认为19）## seed正确匹配参考序列长度小于INT时，该匹配删除
-w INT 最大gap长度（默认为100） ## 删除gap大于INT的匹配，不过，需要了解的是：如果gap长度过大，也会由于在记分矩阵中罚分甚高，影响最终匹配分值而被删除
-d INT 停止延伸参数（默认为100）## 最佳匹配的得分与当前匹配得分之差高于|i-j|*A+INT时，终止seed extension（其中i和j分别为query起始位置和reference的当前位置，A为matching得分），该参数的存在可以避免非必要的seed extension。
-r FLOAT seedSplitRatio,启动re-seeding的控制参数 （默认为1.5） ## MEM区域长度大于 minSeedLen*FLOAT时才会启动re-seeding。该数值越大，产生的seeds越少，因而会加快比对速度但是比对准确度却降低了
-c INT maxOcc （默认为10000）## 指定MEM的最大出现频次，倘若在基因组中重复多次出现说明不可信，将被丢掉
-A INT matchScore（默认为1）## 替换记分矩阵中，正确匹配时得分为 INT
-B INT mmPenalty（默认为4）## 替换记分矩阵中，错误匹配时罚分为 INT
-o INT gapOpenPen（默认为6） ## 替换记分矩阵中，gapOpen罚分为 INT
-E INT gapExtPen （默认为1）## 替换记分矩阵中，gapExtension罚分为 INT
-L INT clipPen (默认为5) ## 剪切罚分为INT
-U INT 非配对 read pair 罚分（默认为9） ## BWA-MEM scores an unpaired read pair as scoreRead1+scoreRead2-INT and scores a paired as scoreRead1+scoreRead2-insertPenalty. It compares these two scores to determine whether we should force pairing
-R STR 完整的read group的头部 ## read group 的ID，会被添加到输出文件的每一个read的头部
-T INT 最低mapping分值 （默认为30）## mapping得分越高，read和参考序列越匹配，低于INT的比对被删除
-a 输出single-end 和 unpaired paired-end的 reads ## 这些比对的结果会被标记为次优
-C Append ## 将FASTA/Q第一行中的注释信息输出到SAM文件，如barcode信息（注意格式符合SAM格式，否则程序会报错）
-H sam结果输出中使用hard clipping ’H’ ## 减少冗余（尤其在做长片段比对和细菌人工染色体序列的比对时）
-v INT 控制输出内容（默认为3）## 0: 禁止将stderr输出；1：只输出error信息；2：只输出errors和warnings；3：所有消息均输出；4：在3的基础上，输出更多debugging信息（该模式）
2.2 后记
原本想要把BWA的3种算法都在这篇文章里整理一下，没想到周末想要偷懒，以后再续记吧，毕竟bwa mem是目前用得最多的工具。
PS：
BWT算法的了解链接（https://www.jianshu.com/p/51a7fc499d28）
Smith-Waterman算法的了解链接（https://www.jianshu.com/p/85465f32273a）
SAM&BAM文件了解链接（softclip和hardclip）https://www.jianshu.com/p/3295e179ef14
网友评论