高级数据结构与算法（一）：Smith-Waterman算法及matlab实现

1.概念

史密斯-沃特曼算法（Smith-Waterman algorithm）是一种进行局部序列比对（相对于全局比对）的算法，该算法的目的不是进行全序列的比对，而是找出两个序列中具有高相似度的片段。可简称为SW算法。

2.过程

2.1算法目的

寻找具有最高相似性的局部序列，进行局部序列匹配。

2.2基本假设

比对的两序列为：

则两序列的长度分别为len(A) = n，Len(B)=m；
s(a,b)：字符a和字符b的相似分数；
H：匹配分数矩阵
= 2 ：一个空位罚分，这里设置为2（可根据需要设置）

2.3算法步骤

1.初始化算法分数矩阵H，使行i表示字符，列j表示字符；

2.计算矩阵中每一项的：

3.回溯，从矩阵H中分数最大的一项开始：
若，则回溯到左上角单元格
若，回溯到左上角、上边、左边中值最大的单元格，若有相同最大值的单元格，优先级按照左上角、上边、左边的顺序

4.根据回溯路径，写出匹配字符串：
若回溯到左上角单元格，将添加到匹配字串A‘，将添加到匹配字串B’；
若回溯到上边单元格，将添加到匹配字串A’，将_添加到匹配字串B’；
若回溯到左边单元格，将_添加到匹配字串A’，将添加到匹配字串B’。

5.得到局部最优匹配序列，结束

3.代码实现

clc;clear;
x1=input("请输入序列一:",'s');
x2=input("请输入序列二:",'s');
W1=2;
H=zeros(length(x2)+1,length(x1)+1);
for i=2:length(x2)+1
    for j=2:length(x1)+1
        if(x1(j-1)==x2(i-1))
            s=3;
        else
            s=-3;
        end
        a=H(i-1,j-1)+s;
        H(i,j)=max([a,H(i-1,j)-W1,H(i,j-1)-W1,0]);
    end
end
[n1,m1]=find(H==max(max(H)'));
n(1)=n1;
m(1)=m1;
i=1;
while(~(n(i)==1||m(i)==1))
    if(x1(m(i)-1)==x2(n(i)-1))
        n(i+1)=n(i)-1;
        m(i+1)=m(i)-1;
        i=i+1;
        continue;
    end
    c=H(n(i)-1,m(i)-1);
    if(c<H(n(i),m(i)-1)&&H(n(i),m(i)-1)>=H(n(i)-1,m(i)))
        n(i+1)=n(i);
        m(i+1)=m(i)-1;
    elseif(H(n(i)-1,m(i))>c&&H(n(i)-1,m(i))>H(n(i),m(i)-1))
        n(i+1)=n(i)-1;
        m(i+1)=m(i);
    else
        n(i+1)=n(i)-1;
        m(i+1)=m(i)-1;
    end
    
    i=i+1;
end
n=n(1:end-1);
m=m(1:end-1);
for i=1:length(n)
    if(H(n(1,i),m(1,i))==0)
        n(2,i)=1;
        m(2,i)=1;
    else
        n(2,i)=0;
        m(2,i)=0;
    end
end
n2=n(1,find(n(2,:)==0))-1;
m2=m(1,find(m(2,:)==0))-1;
for i=length(n2):-1:1
    A(-i+length(n2)+1)=x1(m2(i));
    B(-i+length(n2)+1)=x2(n2(i));
    if(i~=length(n2))
        if(n2(i+1)==n2(i))
            B(-i+length(n2)+1)='_';
        elseif(m2(i+1)==m2(i))
            A(-i+length(n2)+1)='_';
        end
    end
end
X1(1)=string(num2str(0));
x11='abcdefghlmn';
for i=1:length(x1)
    X1(i+1)=strcat(x1(i),x11(i));
end
X2(1)=string(num2str(0));

for i=1:length(x2)
    X2(i+1)=strcat(x2(i),x11(i));
end

h=heatmap(X1,X2,H);
set(h,'Position',[0.1,0.1,0.7,0.8]);
colormap(cool);
for i=1:length(n2)-1
    annotation('arrow',[0.7*(m2(i)+0.5)/(length(x1)+1)+0.1,0.7*(m2(i+1)+0.5)/(length(x1)+1)+0.1],[0.8*(length(x2)+0.5-n2(i))/(length(x2)+1)+0.1,0.8*(length(x2)+0.5-n2(i+1))/(length(x2)+1)+0.1],'Color','w');
end
disp('比对结果为:');
disp(A);
disp(B);











        

4.结果