基于BWT算法的比对软件原理解析（BWA & Bowtie & Bowtie2）

mapping pipline

remove multiple mapping reads的方法

CHIP-seq: Bowtie2、BWA用的比较多
RNA-seq: Tophat、Bsmap
甲基化：BS-seeker

其中BWA & Bowtie & Bowtie2软件均是基于BWT算法

二代测序的特点：
1. 短 250bp
2. 相对较高的精度 1% = Q30（不懂）

三代测序的特点：
1. 长，有structure variation (这也是为什么要做三代测序算法开发的原因之一)
2. 不稳定

1. pairwise alignment

global---NW
local--SW

1. backtrack算法：
$ bwa aln genome read1.fq > aln_sa1.sai
$ bwa aln genome read2.fq > aln_sa2.sai
$ bwa samse genome aln_sa1.sai read1.fq > aln_se.sam
$ bwa sampe genome aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln_pe.sam

2. 比对算法原理：
BWT算法
Seq1
Seq2
两条序列比对：pairwise alignment方法
全局比对：NM算法 局部比对：SW算法

3. 比对到reference基因组的方法：
1)在seq中取出一个较小的seed(30bp?)
2)通过seed找到ref的index,通过index去和附近的序列做pairwise比对
3)通过seed找ref的index的过程有两个算法（短序列比对到基因组）

- 华大SOAP，MAQ。将基因组打断成小段，将位置和序列存成HASH字典。
- BWA，bowtie算法。解决速度问题。

BWT算法：
raw:ACAACG 添加$
M矩阵：
ACAACG$  平移
$ACAACG
G$ACAAC
CG$ACAA
ACG$ACA
AACG$AC
CAACG$A 
M首字母进行排序
T矩阵：
$ACAACG
AACG$AC
ACAACG$
ACG$ACA
CAACG$A
CG$ACAA
G$ACAAC 
T矩阵第一列为F列,最后一列为L列

F列： $AAACCG
L列： GC$AAAC 
关系：对应L是F的前一个。L排序得到F。单字母相对位置不变（L中的第一个C对应F第一个C，以此类推。）
保留L和L中每一个字母的相对位置，1.L可以推出F。2.根据L、F相对位置可以还原ref

image.png

好处是能够穷举出所有的比对情况，所以可以选择全局最优的结果；最大的缺点是比对的非常慢。

2.将query seq打断成seed，比对ref index经历了两代算法

1. 第一代：华大基因--SOAP MAQ  # 缺点是内存占用大，慢，但是找的准
2. 第二代：BWA，Bowtie，Bowtie2 # 解决了速度慢的问题

image.png

3.BWT 最初用于数据压缩

BWT（Burrows-Wheeler Transform ）

假设原始的序列是
（1）Raw ACAACG # 压缩后能还原，且比对次数尽可能少

第一步，在raw seq中加$符号，并平移，形成一个 raw matrix

image.png

第二步，根据Raw Matrix的首字母进行排序，得到转换矩阵Matrix’，默认$符号排在第一位，

image.png

第一列为First 列，F列
最后一列为Last列，L列

F和L的关系

F是L的前面一列；
对L拍完序以后就是F；
单字母的相对位置不变

所以最后只用保存L列和每个字母的相对位置就可以了，根据L列和每个字母的相对位置可以干两件事情：

推出F列
根据L和F列的相对位置可以完整地还原ref相对位置

例如：第一个是L- $,L-$ 对应F- $；F-$ 的前一个是G，L-G对应F-G；F-G的前一个是L-C，依次类推，得到原来的ref：ACAACG$

image.png

bowtie和bowtie2两个版本之间的区别--gap

1. bowtie/BWA # bowtie只允许mismatch，不允许gap；BWA都允许
2. 用bowtie的话是看不到gap open的

bowtie1 不支持gap open

14bp（high quality）---14bp（low quality of high quality）--8bp（real low quality）
分成三断seed，seed1+seed2比对总共的mismatch <= 2，则继续8bp的比对；如果 > 2 直接放弃后面的比对；

bowtie2 支持gap open

第一步，选择seed区域；
20里面选18---
（18+2）+（18+2）+（18+2）+...+(18+2)
保证一个fragment是20，seed 是18bp
或者，10里面选16--
fragment = 16，overlap = 6，

image.png

那么根据BWT算法，就把拆分的seed mapping到基因组的大概位置；
然后把基因组可能mapping上的那段区域挑出来，和query seq做比对（用NW或者SW算法），因为query seq NW和SW允许gap open

image.png

注意

根据gap或者reads quality罚分得到MAPQ，当MAPQ高于一个阈值是认为可以比对上的，低于阈值就认为比对不上，那如果有20个高于阈值的就认为有20个比对结果。
unique map
1）只有一个seq map上
2）只有一段MAPQ特别高，其他的MAPQ特别小，
这两种情况认为是unique map
只有一个map的结果怎么处理；