% Inputs:  1) 2 sequences
%             e.g. Align ATRPLPA STRWLPTK
%          2) Which algorithm to implement: 1 using final cell, 2 using
%             maximum score
% Outputs: 1) Alignment matrix
%          2) Aligned sequences
%
% Notes: 1) Input sequences should have no gaps.
%        2) 1st sequence is located on top & 2nd sequence is located on the
%           side.
%        3) If more than 1 alignment is possible, the algorithm takes only 1
%           alignment, in the following order: Diagonal, Horizontal, Vertical.
%           i.e. if 1 cell has a diagonal arrow & a horizontal arrow, the
%           algorithm will take the diagonal arrow.

function Align(x,y)

% Insert a gap at the beginning of each sequence
x = char(strcat('-', upper(x)));
y = char(strcat('-', upper(y)));

lenx = length(x);
leny = length(y);

% Initialize matrices
S = zeros (lenx, leny);     %S = alignment matrix
A = zeros (lenx, leny);     %A = matrix to store pathway

% Calculations for the gap row & gap column.
for i = 2:lenx
    S(i,1) = S(i-1,1) - 2;
    A(i,1) = 2;             %2 = horizontal arrow
end

for j = 2:leny
    S(1,j) = S(1,j-1) - 2;
    A(1,j) = 3;             %3 = vertical arrow
end
    
% Calculations for the sequences.
for j = 2:leny
    for i = 2:lenx
        H = S(i-1,j) - 2;           %Score calculated horizontally
        V = S(i,j-1) - 2;           %Score calculated vertically
        if x(i) == y(j)
            D = S(i-1,j-1) + 2;     %Score calculated diagonally
        else
            D = S(i-1,j-1) - 1;
        end
        S(i,j) = max(max(H,V),D);   %Take the highest score
        
        % Store the arrow that give the highest score in matrix A
        if S(i,j) == D
            A(i,j) = 1;             %1 = diagonal arrow
        elseif S(i,j) == H
            A(i,j) = 2;             %2 = horizontal arrow
        elseif S(i,j) == V
            A(i,j) = 3;             %3 = vertical arrow
        end
    end
end

disp (S');

disp('There are 2 variations of this algorithm:');
disp('1) Perform traceback from final cell');
disp('2) Perform traceback from cell with maximum score on last row or column');
Algo = input('Please choose 1 or 2: ');

while (Algo ~= 1) & (Algo ~= 2)
    Algo = input('Invalid choice. Please choose 1 or 2: ');
end

pos = lenx+leny;

if Algo == 2
    % Locate cell with maximum score
    for p = 1:lenx-1
        if S(p,j) > S(i,j)
            i = p;
        end
    end

    for q = 1:leny-1
        if S(lenx,q) > S(i,j)
            i = lenx;
            j = q;
        end
    end

    % Align with trailing gaps if necessary
    if i ~= lenx
        for r = lenx:-1:i+1
            newx(pos) = x(r);
            newy(pos) = '-';
            pos = pos-1;
        end
    elseif j ~= leny
        for r = leny:-1:j+1
            newx(pos) = '-';
            newy(pos) = y(r);
            pos = pos-1;
        end
    end
end

% Traceback
while ((i ~= 1) | (j ~= 1)) & (pos ~= 0)
    if A(i,j) == 1          %Diagonal arrow, align x with y
        newx(pos) = x(i);
        newy(pos) = y(j);
        i = i-1;
        j = j-1;
    elseif A(i,j) == 2      %Horizontal arrow, align x with gap
        newx(pos) = x(i);
        newy(pos) = '-';
        i = i-1;
    elseif A(i,j) == 3      %Vertical arrow, align gap with y
        newx(pos) = '-';
        newy(pos) = y(j);
        j = j-1;
    end
    pos = pos-1;
end

disp (' ');
disp (newx);
disp (newy);