Global & Local Alignment

[asylu3]

Global

[src]sub computeMatrix;
sub printAllAlignments;
sub printMatrixResult;

# string sequence A and B
$ssa = 'GTATCTT';
$ssb = 'TGAGT';

# convert string sequence to array
@sa = split('',$ssa);
@sb = split('',$ssb);

# half of the substitution metrix
$sm = {
A => {A=> 6},
T => {A=>-1, T=> 4},
C => {A=>-2, T=>-1, C=> 5},
G => {A=>-3, T=> 1, C=>-2, G=>3}};
# generate the other half
while (my ($x,$X) = each %$sm) {
while (my ($y,$z) = each %$X) {
$sm->{$y}{$x} = $z if ($x ne $y);
}
}
$gap = -4; # gap penalty

# initialize the alignment matrix and the direction matrix
$m->[0][0] = 0;
$d->[0][0] = [0,0,0]; # diag,horz,vert

computeMatrix;
printMatrixResult;

# find all alignments
@ra = ();
@rb = ();
$i = 0; # the alignment number
printAllAlignments scalar(@sa), scalar(@sb), 0, @ra, 0, @rb;
print "n";

# -------------------------------------------
sub computeMatrix {
# fill in the first row (horizontal direction [1])
for (my $x=1; $x<=@sa; $x++) {
$m->[$x][0] = $gap * $x;
$d->[$x][0] = [0,1,0];
}

# fill in the first column (vertical direction [2])
for (my $y=1; $y<=@sb; $y++) {
$m->[0][$y] = $gap * $y;
$d->[0][$y] = [0,0,1];
}

# fill in the rest of the matrices
for (my $x=1; $x<=@sa; $x++) {
for (my $y=1; $y<=@sb; $y++) {

# compute all three directions
my $horz = $m->[$x-1][$y] + $gap;
my $vert = $m->[$x][$y-1] + $gap;
my $diag = $m->[$x-1][$y-1]
+ $sm->{$sa[$x-1]}{$sb[$y-1]};

# choose the max among the three directions
my $max = $diag;
$max = $horz if $max < $horz;
$max = $vert if $max < $vert;
$m->[$x][$y] = $max;
$d->[$x][$y] = [0,0,0];
$d->[$x][$y][0] = 1 if $max == $diag;
$d->[$x][$y][1] = 1 if $max == $horz;
$d->[$x][$y][2] = 1 if $max == $vert;
}
}
}

# -------------------------------------------
sub printMatrixResult {
# print the alignment matrix
print "Matrix:n Sa ";
# print the header row
foreach my $x (@sa) {
print "[ $x ]";
}
print "n";
for (my $y=0; $y<=@sb; $y++) {
# print the header column
if ($y>0) {
print "[ $sb[$y-1] ]";
} else {
print " Sb ";
}
# print the values
for (my $x=0; $x<=@sa; $x++) {
if (!$d->[$x][$y][0] && !$d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[ '; }
if (!$d->[$x][$y][0] && !$d->[$x][$y][1] && $d->[$x][$y][2]) { print '[| '; }
if (!$d->[$x][$y][0] && $d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[- '; }
if (!$d->[$x][$y][0] && $d->[$x][$y][1] && $d->[$x][$y][2]) { print '[+ '; }
if ( $d->[$x][$y][0] && !$d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[ '; }
if ( $d->[$x][$y][0] && !$d->[$x][$y][1] && $d->[$x][$y][2]) { print '[V '; }
if ( $d->[$x][$y][0] && $d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[> '; }
if ( $d->[$x][$y][0] && $d->[$x][$y][1] && $d->[$x][$y][2]) { print '[* '; }
printf '%3d]', $m->[$x][$y];
}
print "n";
}
}

# -------------------------------------------
sub printAllAlignments {
my $x = shift @_;
my $y = shift @_;
my @ra = splice(@_,0,shift);
my @rb = splice(@_,0,shift);
if ($d->[$x][$y][0]) { #while see arrow
my @a = @ra;
my @b = @rb;
unshift(@a,$sa[$x-1]);
unshift(@b,$sb[$y-1]);
printAllAlignments $x-1, $y-1, scalar(@a), @a, scalar(@, @b;
}#no arrow
if ($d->[$x][$y][1]) {
my @a = @ra;
my @b = @rb;
unshift(@a,$sa[$x-1]);
unshift(@b,'-');
printAllAlignments $x-1, $y, scalar(@a), @a, scalar(@, @b;
}
if ($d->[$x][$y][2]) {
my @a = @ra;
my @b = @rb;
unshift(@a,'-');
unshift(@b,$sb[$y-1]);
printAllAlignments $x, $y-1, scalar(@a), @a, scalar(@, @b;
}
if (!$d->[$x][$y][0] && !$d->[$x][$y][1] && !$d->[$x][$y][2]) {
# print out the answer
$i++; # global variable
$scr = 0; # score for each alignment
print "nAlignment: $inSa:";
foreach my $x (@ra) { print " $x"; }
print "nSb:";
foreach my $y (@rb) { print " $y"; }
print "nScr";
for (my $z=0; $z<@ra; $z++) {
if (($ra[$z] eq '-') || ($rb[$z] eq '-')) {
$scr += $gap;
printf '%4d', $gap;
} else {
$scr += $sm->{$ra[$z]}{$rb[$z]};
printf '%4d',$sm->{$ra[$z]}{$rb[$z]};
}
}
print " = $scr";
}
}
[/src]

Local
[src]sub computeMatrix;
sub printAllAlignments;
sub printMatrixResult;

# string sequence A and B
$ssa = 'GTACGTG';
$ssb = 'AGCGA';

# convert string sequence to array
@sa = split('',$ssa);
@sb = split('',$ssb);

# half of the substitution metrix
$sm = {
A => {A=> 6},
T => {A=>-1, T=> 4},
C => {A=>-2, T=>-1, C=> 5},
G => {A=>-3, T=> 1, C=>-2, G=>3}};
# generate the other half
while (my ($x,$X) = each %$sm) {
while (my ($y,$z) = each %$X) {
$sm->{$y}{$x} = $z if ($x ne $y);
}
}
$gap = -4; # gap penalty

# initialize the alignment matrix and the direction matrix
$m->[0][0] = 0;
$d->[0][0] = [0,0,0]; # diag,horz,vert

$top = computeMatrix;
printMatrixResult;

$i = 0; # the alignment number

# find top value
for (my $p=0; $p<=@sa; $p++) {
for (my $q=0; $q<=@sb; $q++) {
if ($m->[$p][$q] == $top) {
# find all alignments
@ra = @sa;
@rb = @sb;
@rc = ();
splice @ra, 0, $p;
splice @rb, 0, $q;
while (scalar(@ra)<scalar(@rb)) { push @ra, '.'; }
while (scalar(@rb)<scalar(@ra)) { push @rb, '.'; }
while (scalar(@rc)<scalar(@ra)) { push @rc, '.'; }
printAllAlignments $p, $q, scalar(@rc), @ra, @rb, @rc;
}
}
}

# -------------------------------------------
sub computeMatrix {
# fill in the first row (horizontal direction [1])
for (my $x=1; $x<=@sa; $x++) {
$m->[$x][0] = 0;
$d->[$x][0] = [0,0,0];
}

# fill in the first column (vertical direction [2])
for (my $y=1; $y<=@sb; $y++) {
$m->[0][$y] = 0;
$d->[0][$y] = [0,0,0];
}

my $top = 0; # top score overall

# fill in the rest of the matrices
for (my $x=1; $x<=@sa; $x++) {
for (my $y=1; $y<=@sb; $y++) {

# compute all three directions
my $horz = $m->[$x-1][$y] + $gap;
my $vert = $m->[$x][$y-1] + $gap;
my $diag = $m->[$x-1][$y-1]
+ $sm->{$sa[$x-1]}{$sb[$y-1]};

# choose the max among the three directions
my $max = 0;
$max = $diag if $max < $diag;
$max = $horz if $max < $horz;
$max = $vert if $max < $vert;
$m->[$x][$y] = $max;
$d->[$x][$y] = [0,0,0];
$d->[$x][$y][0] = 1 if $max == $diag;
$d->[$x][$y][1] = 1 if $max == $horz;
$d->[$x][$y][2] = 1 if $max == $vert;
$top = $ max if $top < $max; # remember the max
}
}
return $top;
}

# -------------------------------------------
sub printMatrixResult {
# print the alignment matrix
print "Matrix:n Sa ";
# print the header row
foreach my $x (@sa) {
print "[ $x ]";
}
print "n";
for (my $y=0; $y<=@sb; $y++) {
# print the header column
if ($y>0) {
print "[ $sb[$y-1] ]";
} else {
print " Sb ";
}
# print the values
for (my $x=0; $x<=@sa; $x++) {
if (!$d->[$x][$y][0] && !$d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[ '; }
if (!$d->[$x][$y][0] && !$d->[$x][$y][1] && $d->[$x][$y][2]) { print '[| '; }
if (!$d->[$x][$y][0] && $d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[- '; }
if (!$d->[$x][$y][0] && $d->[$x][$y][1] && $d->[$x][$y][2]) { print '[+ '; }
if ( $d->[$x][$y][0] && !$d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[ '; }
if ( $d->[$x][$y][0] && !$d->[$x][$y][1] && $d->[$x][$y][2]) { print '[V '; }
if ( $d->[$x][$y][0] && $d->[$x][$y][1] && !$d->[$x][$y][2]) { print '[> '; }
if ( $d->[$x][$y][0] && $d->[$x][$y][1] && $d->[$x][$y][2]) { print '[* '; }
printf '%3d]', $m->[$x][$y];
}
print "n";
}
}

# -------------------------------------------
sub printAllAlignments {
my $x = shift @_;
my $y = shift @_;
my $r = shift @_;
my @ra = splice(@_,0,$r);
my @rb = splice(@_,0,$r);
my @rc = splice(@_,0,$r);
if ($d->[$x][$y][0]) {
my @a = @ra;
my @b = @rb;
my @c = @rc;
unshift(@a,$sa[$x-1]);
unshift(@b,$sb[$y-1]);
unshift(@c,$sm->{$sa[$x-1]}{$sb[$y-1]});
printAllAlignments $x-1, $y-1, scalar(@c), @a, @b, @c;
}
if ($d->[$x][$y][1]) {
my @a = @ra;
my @b = @rb;
my @c = @rc;
unshift(@a,$sa[$x-1]);
unshift(@b,'-');
unshift(@c,$gap);
printAllAlignments $x-1, $y, scalar(@c), @a, @b, @c;
}
if ($d->[$x][$y][2]) {
my @a = @ra;
my @b = @rb;
my @c = @rc;
unshift(@a,'-');
unshift(@b,$sb[$y-1]);
unshift(@c,$gap);
printAllAlignments $x, $y-1, scalar(@c), @a, @b, @c;
}
if (!$d->[$x][$y][0] && !$d->[$x][$y][1] && !$d->[$x][$y][2]) {
for (my $u=$x; $u>0; $u--) { unshift @ra, $sa[$u-1]; }
for (my $v=$y; $v>0; $v--) { unshift @rb, $sb[$v-1]; }
while (scalar(@ra)<scalar(@rb)) { unshift @ra, '.'; }
while (scalar(@rb)<scalar(@ra)) { unshift @rb, '.'; }
while (scalar(@rc)<scalar(@ra)) { unshift @rc, '.'; }
# print out the answer
$i++; # global variable
$scr = 0; # score for each alignment
print "nAlignment: $inSa:";
foreach my $x (@ra) { print " $x"; }
print "nSb:";
foreach my $y (@rb) { print " $y"; }
print "nScr";
foreach my $z (@rc) {
if ($z eq '.') {
print " .";
} else {
printf '%4d', $z;
$scr += $z;
}
}
print " = $scrn";
}
}
[/src]