Closed:Need help to add letter into matrix
0
0
Entering edit mode
6.5 years ago

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 • 102 views
ADD COMMENT
This thread is not open. No new answers may be added
Traffic: 1968 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6