Fast and accurate short read alignment with Burrows-Wheeler transform
本文主要将生物信息中,BWA的算法设计。
1. BWT变换
假设 T=RUOSHUI
在 T后面加入比T任意字符都小的字符$, X = RUOSHUI$
字符串X循环移位后得到的字符串集合组合成一个N*N的矩阵
0 RUOSHUI$
1 UOSHUI$R
2 OSHUI$RU
3 SHUI$RUO
4 HUI$RUOS
5 UI$RUOSH
6 I$RUOSHU
7 $RUOSHUI
按字符表进行排序后得到矩阵M
0 7 $RUOSHUI
1 4 HUI$RUOS
2 6 I$RUOSHU
3 2 OSHUI$RU
4 0 RUOSHUI$
5 3 SHUI$RUO
6 5 UI$RUOSH
7 1 UOSHUI$R
得到BWT变换字符串L=“ISUU$OHR”;
后缀数组S=[7,4,6,2,0,3,5,1]
2. 后缀数组域
如果,字符串W是X的字串,在X的后缀数组中找到数组域,
R−−(W)
,
R¯¯¯(W)
,分别表述W在后缀数组的区间。
如果W是空字符串,
R−−(W)
=1,
R¯¯¯(W)
=|W|-1
C(a),表示在X的字符串,首次出现在矩阵M的行号,如表
a | $ | H | I | O | R | S | U |
---|---|---|---|---|---|---|---|
C(a) | 0 | 1 | 2 | 3 | 4 | 5 | 6 |
O(a, i ), 表示在L[i]中字符a出现的次数,如表
L | I | S | U | U | $ | O | H | R |
---|---|---|---|---|---|---|---|---|
0 | 1 | 2 | 3 | 4 | 5 | 6 | 7 | |
$ | 0 | 0 | 0 | 0 | 1 | 1 | 1 | 1 |
H | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 1 |
I | 1 | 1 | 1 | 1 | 1 | 1 | 1 | 1 |
O | 0 | 0 | 0 | 0 | 0 | 1 | 1 | 1 |
R | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 1 |
u | 0 | 0 | 1 | 2 | 2 | 2 | 2 | 2 |
3. 精确查找
aW, 公式
R−−(aW)=C(a)+O(a,R−−(W)−1)+1
R¯¯¯(aW)=C(a)+O(a,R¯¯¯(W))
精确求出aW的后缀数域
不妨,求字符串Y的后缀数域,采用上面公式求出。这个方式称为后端查找。
4.非精确查找
非精确查找表示字符串,匹配不超过z个差别。一种是错配,一种出现空穴。采用一种递归的方式计算。具体的方法如下
InexactSearch( W, z )
CalculateD( W )
return inexRecur( W, |W| - 1, z, 1, |X| - 1 )
CalculateD( W )
j = 0
z = 0
for ( i = 0; i < |W| ; ++i ){
if ( W[j,i] is not substring of W ){
z = z + 1
j = i + 1
}
D(i) = z
}
inexRecur( W, i, z, k, l )
if z < D(i)
return NULL
if i < 0
return {[k, l]}
I = NULL
I = I U inexRecur( W, i - 1, z - 1, k , l )
for each b in {A, G, C, T }
k = C(b) + O ( b, k - 1 ) + 1
l = C(b) + O ( b, l )
if k <= l
I = I U inexRecur( W, i , z- 1 , k , l )
if b = W[i]
I = I U inexRecur( W, i - 1, z, k, l )
else
I = I U inexRecur( W, i - 1, z - 1, k, l )
return I