Entering edit mode
6.5 years ago
danghoanghamy
•
0
Hi, I have a script use to sequence alignment, but when running it just only show the matrix table without letter of the sequence. I want to add letter in but i don't know how.
Hi, I have a script
use strict;
use warnings;
use List::Util 'max';
STDOUT->autoflush;
my ($seq1, $seq2) = qw/ ACGGAT AGATAA /;
my ($len1, $len2) = map length, $seq1, $seq2;
my @seq1 = $seq1 =~ /./g;
my @seq2 = $seq2 =~ /./g;
my @matrix1;
$matrix1[$_][0] = -$_ for 0 .. $len2;
$matrix1[0][$_] = -$_ for 0 .. $len1;
for my $x ( 1 .. $len2 ) { # Rows (bases of sequence 2)
for my $y ( 1 .. $len1 ) { # Columns (bases of sequence 1)
my $match = $seq1[$y-1] eq $seq2[$x-1];
my @scores = (
$matrix1[$x-1][$y] - 1, # Up
$matrix1[$x][$y-1] - 1, # Left
$match ? $matrix1[$x-1][$y-1] + 1 : $matrix1[$x-1][$y-1] - 1, # Diagonal
);
$matrix1[$x][$y] = max @scores;
}
}
print "******\n";
print "MATRIX\n";
print "******\n";
for my $i ( 0 .. $len2 ) {
my $row = $matrix1[$i];
print join(' ', map { sprintf '%3d', $_ } @$row), "\n";
}
# my $seq1='ACGGAT';
# my $seq2='AGATAA';
my $MATCH = 1; # +1 for letters that match
my $MISMATCH = -1; # -1 for letters that mismatch
my $GAP = -1; # -1 for any gap
my @matrix2;
$matrix2[0][0]{score} = 0;
$matrix2[0][0]{pointer} = "none";
for(my $j = 1; $j <= length($seq1); $j++) {
$matrix2[0][$j]{score} = $GAP * $j;
$matrix2[0][$j]{pointer} = "left";
}
for (my $i = 1; $i <= length($seq2); $i++) {
$matrix2[$i][0]{score} = $GAP * $i;
$matrix2[$i][0]{pointer} = "up";
}
# fill
for(my $i = 1; $i <= length($seq2); $i++) {
for(my $j = 1; $j <= length($seq1); $j++) {
my ($diagonal_score, $left_score, $up_score);
# calculate match score
my $letter1 = substr($seq1, $j-1, 1);
my $letter2 = substr($seq2, $i-1, 1);
if ($letter1 eq $letter2) {
$diagonal_score = $matrix2[$i-1][$j-1]{score} + $MATCH;
}
else {
$diagonal_score = $matrix2[$i-1][$j-1]{score} + $MISMATCH;
}
# calculate gap scores
$up_score = $matrix2[$i-1][$j]{score} + $GAP;
$left_score = $matrix2[$i][$j-1]{score} + $GAP;
# choose best score
if ($diagonal_score >= $up_score) {
if ($diagonal_score >= $left_score) {
$matrix2[$i][$j]{score} = $diagonal_score;
$matrix2[$i][$j]{pointer} = "diagonal";
}
else {
$matrix2[$i][$j]{score} = $left_score;
$matrix2[$i][$j]{pointer} = "left";
}
} else {
if ($up_score >= $left_score) {
$matrix2[$i][$j]{score} = $up_score;
$matrix2[$i][$j]{pointer} = "up";
}
else {
$matrix2[$i][$j]{score} = $left_score;
$matrix2[$i][$j]{pointer} = "left";
}
}
}
}
# trace-back
my $align1 = "";
my $align2 = "";
while (1) {
last if $matrix2[$len2][$len1]{pointer} eq "none"; # ends at first cell of matrix
if ($matrix2[$len2][$len1]{pointer} eq "diagonal") {
$align1 .= substr($seq1, $len1-1, 1);
$align2 .= substr($seq2, $len2-1, 1);
$len2--;
$len1--;
}
elsif ($matrix2[$len2][$len1]{pointer} eq "left") {
$align1 .= substr($seq1, $len1-1, 1);
$align2 .= "-";
$len1--;
}
elsif ($matrix2[$len2][$len1]{pointer} eq "up") {
$align1 .= "-";
$align2 .= substr($seq2, $len2-1, 1);
$len2--;
}
}
$align1 = reverse $align1;
$align2 = reverse $align2;
print "********\n";
print "SEQUENCE\n";
print "********\n";
print "$align1\n";
print "$align2\n";