function align()
% Demonstration of DNA alignment with the Needleman-Wunsch dynamic
% programming algorithm.

seq2 = 'GAATTCAGTTA';
seq1 = 'GGATCGA';

n = length(seq1);
m = length(seq2);

score = zeros(n+1, m+1);
w = -2;
back = cell(n+1, m+1);

% fill the score matrix
for i=2:n+1
    for j=2:m+1
        sij = -1;
        if seq1(i-1) == seq2(j-1)
            sij = 2;
        end
        a = [score(i-1,j-1)+sij, score(i,j-1)+w, score(i-1,j)+w];
        [smax, imax] = max(a);
        score(i,j) = smax;
        switch imax
            case 1
                back{i,j} = [i-1, j-1];
            case 2
                back{i,j} = [i, j-1];
            case 3
                back{i,j} = [i-1, j];
        end
    end
end

score
bar3(score);
title('Score Matrix');
set(gca,'XTickLabel', str2cell(seq2));
set(gca,'YTickLabel', str2cell(seq1));
zlabel('Score');

% back tracing
i = n+1;
j = m+1;
align = cell(1,max([n, m])-1);
align{1} = [i, j];
k = 2;
while i > 2 && j > 2
   ij = back{i,j};
   align{k} = ij;
   k = k + 1;
   i = ij(1);
   j = ij(2);
end

iold = 0;
jold = 0;
l = 1;
for k=length(align):-1:1
    ij = align{k};
    i = ij(1);
    j = ij(2);
    if i ~= iold
        c1 = seq1(i-1);
    else
        c1 = '_';
    end
    if j ~= jold
        c2 = seq2(j-1);
    else
        c2 = '_';
    end
    iold = i;
    jold = j;
    aseq1(l) = c1;
    aseq2(l) = c2;
    l = l + 1;
end

aseq1
aseq2

% draw the backtrace matrix
n = length(align)-1;
x = zeros(1,n);
y = zeros(1,n);
u = zeros(1,n);
v = zeros(1,n);
for k=1:n
    from = align{k};
    y(k) = from(1);
    x(k) = from(2);
    to = align{k+1};
    v(k) = to(1) - y(k);
    u(k) = to(2) - x(k);
end
figure;
pcolor(score);
colormap(gray);
hold on;
quiver(x,y,u,v,0.25,'LineWidth',3,'Color','r');
title('Backtrace');
set(gca,'XTickLabel', str2cell(seq2));
set(gca,'YTickLabel', str2cell(seq1));
colorbar;

end

function c = str2cell(s)
% Convert a string to a cell array of size 1 by N, where N is the
% string length. Each cell contains one character of the string.
c = cell(1,length(s)+1);
c(1) = cellstr(' ');
for i=1:length(s)
    c(i+1) = cellstr([s(i)]);
end
end