#!/usr/bin/perl -w

$seq1 = "ATCGAT";
$seq2 = "ATACGT";

$MATCH    =  2;
$MISMATCH = -1;
$GAP      = -2;


# Initialise the score matrix and the trace matrix

$score_matrix[0][0] = 0;
$trace_matrix[0][0] = "STOP";

for ( $row = 1; $row <= length($seq1); $row++ ) {
    $score_matrix[$row][0] = $score_matrix[$row-1][0] + $GAP;
    $trace_matrix[$row][0] = "UP";
}

for ( $column = 1; $column <= length($seq2); $column++ ) {
    $score_matrix[0][$column] = $score_matrix[0][$column-1] + $GAP;
    $trace_matrix[0][$column] = "LEFT";
}


# Fill the score matrix and the trace matrix

for ( $row = 1; $row <= length($seq1); $row++ ) {
    for ( $column = 1; $column <= length($seq2); $column++ ) {
        if ( substr($seq1, $row-1, 1) eq substr($seq2, $column-1, 1) ) {
            $diagonal_score = $score_matrix[$row-1][$column-1] + $MATCH;
        } else {
            $diagonal_score = $score_matrix[$row-1][$column-1] + $MISMATCH;
        }
        $left_score = $score_matrix[$row][$column-1] + $GAP;
        $up_score   = $score_matrix[$row-1][$column] + $GAP;

        if ( ( $diagonal_score >= $left_score ) &&
             ( $diagonal_score >= $up_score ) ) {
            $score_matrix[$row][$column] = $diagonal_score;
            $trace_matrix[$row][$column] = "DIAGONAL";
        } elsif ( $left_score >= $up_score ) {
            $score_matrix[$row][$column] = $left_score;
            $trace_matrix[$row][$column] = "LEFT";
        } else {
            $score_matrix[$row][$column] = $up_score;
            $trace_matrix[$row][$column] = "UP";
        }
    }
}


# Print the score matrix

print "Score matrix:\n      ";
for ( $column=0 ; $column<length($seq2) ; ++$column ) {
    printf  "    " . substr($seq2, $column, 1);
}
print "\n";
for ( $row = 0; $row <= length($seq1); $row++ ) {
    if ( $row==0 ) {
        print " ";
    } else {
        print substr($seq1, $row-1, 1);
    }
    for ( $column = 0; $column <= length($seq2); $column++ ) {
        printf( "%5d",  $score_matrix[$row][$column] );
    }
    print "\n";
}


# Trace back from the bottom-right cell

$aligned1 = "";
$aligned2 = "";

$row = length($seq1);
$column = length($seq2);

while ( $trace_matrix[$row][$column] ne "STOP" ) {
    if ( $trace_matrix[$row][$column] eq "DIAGONAL" ) {
        $aligned1 = substr($seq1, $row-1, 1) . $aligned1;
        $aligned2 = substr($seq2, $column-1, 1) . $aligned2;
        $row--;
        $column--;
    } elsif ( $trace_matrix[$row][$column] eq "LEFT" ) {
        $aligned1 = "-" . $aligned1;
        $aligned2 = substr($seq2, $column-1, 1) . $aligned2;
        $column--;
    } elsif ( $trace_matrix[$row][$column] eq "UP" ) {
        $aligned1 = substr($seq1, $row-1, 1) . $aligned1;
        $aligned2 = "-" . $aligned2;
        $row--;
    }
}

print "\n";
print "$aligned1\n";
print "$aligned2\n";