参考：
踏踏实实做技术：BWA，Bowtie，Bowtie2的比对算法推导
mapping pipline
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-对应F-
的前一个是G，L-G对应F-G；F-G的前一个是L-C，依次类推，得到原来的ref：ACAACG$
image.png
image.png
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的结果怎么处理；
网友评论