This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Need help to add letter into matrix

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";
alignment sequence sequencing matrix perl

I added code markup to your post for increased readability. You can do this by selecting the text and clicking the 101010 button. When you compose or edit a post that button is in your toolbar, see image below:

101010 Button

thank you. I don't know how to make it clear since i am a newbie of this forum.

0 answers

No answers yet.

Log in to answer this question.