一、参数说明
T=ACGTCTCGAGACGT
|T|=14
T[i]=第i个碱基
T[i,j]=第i到第j个碱基的字符串
Ti 整个的字符串
S:S(i)是第i小的数组的位置
B[i]=尾缀数组
C(a)共四个值,分别为C(A)C(C)C(G)C(T) C(A)={0<=i<=n-1:T[i]<A} 的个数
O(a,i)=从B[]表从0~i中,a 的个数,大小为4*N
P·W表示:字符串P和W的直接连接
Pa 表示:字符串P和碱基a 的连接
P(上横线)表示P的补
二、算法举例
问题:reference:ACGTCTC; read:GTC;
由read的T开始去匹配,那么在bi-interval模型中他是怎么实现的呢?
先对reference取反:TGCAGAG,再前后调换位置:GAGACGT
reference 右移:
- ACGTCTCGAGACGT
- TACGTCTCGAGACG
- GTACGTCTCGAGAC
- CGTACGTCTCGAGA
- ACGTACGTCTCGAG
- GACGTACGTCTCGA
- AGACGTACGTCTCG
- GAGACGTACGTCTC
- CGAGACGTACGTCT
- TCGAGACGTACGTC
- CTCGAGACGTACGT
- TCTCGAGACGTACG
- GTCTCGAGACGTAC
- CGTCTCGAGACGTA
排序:
- ACGTACGTCTCGAG
- ACGTCTCGAGACGT
- AGACGTACGTCTCG
- CGAGACGTACGTCT
- CGTACGTCTCGAGA
- CGTCTCGAGACGTA
- CTCGAGACGTACGT
- GACGTACGTCTCGA
- GAGACGTACGTCTC
- GTACGTCTCGAGAC
- GTCTCGAGACGTAC
- TACGTCTCGAGACG
- TCGAGACGTACGTC
- TCTCGAGACGTACG
BWT表建立完毕,从read(GTC)的中间碱基T开始匹配,其interval为:[12,1,3]。
- Backward匹配:
算法分析:
b从0~5表示$,0,1,2,3,4,N。先算出k上面有多少b,加上O(b)记为Kb,同事算出在[k,k+s]这个范围内,有多少b,记为Sb。
L递推求,L0与W的L相同,L4=L0+S0,L3=L4+S4,L2=L3+S3......
针对具体问题:T:bi-interval:[12,1,3] a:G
K0=C($~0)+O(0,11)=0 ; S0=0
K1=C(A~1)+O(1,11)=1+3=4 ; S1=O(1,14)- O(1,11)=0(不存在)
K2=C(C~2)+O(2,,11)=4+3=7 ; S2=O(2,13)- O(2,11)=1 (有1个)
K3=C(G~3)+O(3,11)=8+2=10 ; S3=O(3,13)-O(3,11)=2(有2个)
K4=C(T~4)+O(4,11)=12+3=15 ; S4=0 (不存在)
K5~没有奇异值,不算。
L0=L=1;
L4=L0+S0=1+0=1
L3=L4+S4=1(即A右边有C)
L2=L3+S3=1+2=3~即T的反变换:A右边有G的是第3行
L1=L2+S2=3+1=4~即T的反变换:A右边有T的是第4行。但s4=0,也是不存在
根据输入的a=G,则得到的bi-interval[k3,l2,s3]=[10,1,2]
- Forward匹配:
Forward匹配要做的就是利用backward 的公式,改变输入参数,得到forward(向右)匹配的bi-interval结果。
1.首先:将a=C变成a(补)=G,
2.再将原bi-interval的k和l换位置,得到[1,10,2]。
这样就得到了一个新的backward匹配工作。
针对具体问题:GT(补)=AC,它的bi-interval是:[1,10,2] ; a=C->a(补)=G
L0=C($~0)+O(0,0)=0 ; S0=0
L1=C(A~1)+O(1,0)=1+0=1 ; S1=O(1,2)- O(1,0)=0(不存在)
L2=C(C~2)+O(2,,0)=4+0=4 ; S2=O(2,2)- O(2,0)=0 (不存在)
L3=C(G~3)+O(3,0)=8+0=8 ; S3=O(3,2)-O(3,0)=1(有1个)
L4=C(T~4)+O(4,0)=12+0=12 ; S4=O(4,2)-O(4,0)=1 (有1个)
L5~没有奇异值,不算。
K0=K=10;
k4=K0+S0=10+0=10(即GT右边有A)
k3=k4+S4=10+1=11(即GT右边有C)
K2=k3+S3=11+1=12~S2=0所以不存在
K1=K2+S2=12+0=12~S1=0所以不存在
根据输入的a=G,输出的时候,要调换K和L 的位置,则得到的bi-interval[k3,l3,s3]=[11,8,1]
匹配流程走完。