Re: Sequence alignment

I guess u re talking about the Needleman Wunsch algorithm. I can help you in working out some steps but I think u need to define the gap penalty, match and mismatch score prior to get the alignment.

1. Build an empty matrix where you can calculate the score and also have a pointer for the tracing back so that you can calculate the matches and subsequently get the gaps/mismatches.

2. Filling the matrix and then tracing back the matrix

my @matrix;
$matrix[0][0]{score}   = 0;
$matrix[0][0]{pointer} = "none";
for(my $j = 1; $j <= length($seq1); $j++) {
    $matrix[0][$j]{score}   = $GAP * $j;
    $matrix[0][$j]{pointer} = "left";
}
for (my $i = 1; $i <= length($seq2); $i++) {
    $matrix[$i][0]{score}   = $GAP * $i;
    $matrix[$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 = $matrix[$i-1][$j-1]{score} + $MATCH;
        }
        else {
            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MISMATCH;
        }

        # calculate gap scores
        $up_score   = $matrix[$i-1][$j]{score} + $GAP;
        $left_score = $matrix[$i][$j-1]{score} + $GAP;

        # choose best score
        if ($diagonal_score >= $up_score) {
            if ($diagonal_score >= $left_score) {
                $matrix[$i][$j]{score}   = $diagonal_score;
                $matrix[$i][$j]{pointer} = "diagonal";
            }
        else {
                $matrix[$i][$j]{score}   = $left_score;
                $matrix[$i][$j]{pointer} = "left";
            }
        } else {
            if ($up_score >= $left_score) {
                $matrix[$i][$j]{score}   = $up_score;
                $matrix[$i][$j]{pointer} = "up";
            }
            else {
                $matrix[$i][$j]{score}   = $left_score;
                $matrix[$i][$j]{pointer} = "left";
            }
        }
    }
}
# trace-back

my $align1 = "";
my $align2 = "";

# start at last cell of matrix
my $j = length($seq1);
my $i = length($seq2);

while (1) {
    last if $matrix[$i][$j]{pointer} eq "none"; # ends at first cell o
+f matrix

    if ($matrix[$i][$j]{pointer} eq "diagonal") {
        $align1 .= substr($seq1, $j-1, 1);
        $align2 .= substr($seq2, $i-1, 1);
        $i--;
        $j--;
    }
    elsif ($matrix[$i][$j]{pointer} eq "left") {
        $align1 .= substr($seq1, $j-1, 1);
        $align2 .= "-";
        $j--;
    }
    elsif ($matrix[$i][$j]{pointer} eq "up") {
        $align1 .= "-";
        $align2 .= substr($seq2, $i-1, 1);
        $i--;
    }    
}
$align1 = reverse $align1;
$align2 = reverse $align2;
print "$align1\n";
print "$align2\n";
[download]

For match = 1, mismatch = -1 and gap = -1.

Input

ATGTAGACCTAGATCATGATGACTGATGAT
ATTACCGATGACTGATGACTGATGACTGAT
[download]

Output

ATGTAGACCTA-GA-TCATGA-TGACTGA-TGAT
AT-T--ACCGATGACTGATGACTGA-TGACTGAT
[download]

Comment on Re: Sequence alignment Select or Download Code

Replies are listed 'Best First'.

Re^2: Sequence alignment
by bobosm (Initiate) on May 03, 2016 at 09:40 UTC

How we simply add the code for display max score? I try but do not print correctly.

[reply]

Re^2: Sequence alignment
by Samman_Mahmoud (Initiate) on Oct 28, 2013 at 12:13 UTC

why this code doesn't work with me !!!

[reply]

Re^3: Sequence alignment

by marto (Cardinal) on Oct 28, 2013 at 12:22 UTC

"doesn't work" isn't an error report. You'll have to provide more information when asking for help. See How do I post a question effectively? and PerlMonks for the Absolute Beginner.

[reply]

Re^3: Sequence alignment

by Athanasius (Archbishop) on Oct 28, 2013 at 12:34 UTC

Although marto is quite right, perhaps I can make an educated guess here?

The code as given by snape is incomplete. When I add this before the code:

#! perl
use strict;
use warnings;

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

chomp(my $seq1 = <DATA>);
chomp(my $seq2 = <DATA>);
[download]

and this after the code:

__DATA__
ATGTAGACCTAGATCATGATGACTGATGAT
ATTACCGATGACTGATGACTGATGACTGAT
[download]

I get the same output as reported by snape.

Disclaimer: I have no knowledge of the Needleman Wunsch algorithm (if that’s what this code is implementing), nor of what it’s for. :-)

Hope that helps,

Athanasius <°(((>< contra mundum Iustus alius egestas vitae, eros Piratica,

[reply]
[d/l]
[select]