#!/usr/bin/perl use strict; # Implementation of Needleman-Wusch # Copyright Dan DeBlasio 2012 # Created 8 Oct 2012 # get input from the user my $string1 = shift @ARGV; my $string2 = shift @ARGV; # Set up DP value and DP direction tables my @DP = ([(0) x length($string1)+1], [(0) x length($string2)+1]); my @DP_dir = ([(0) x length($string1)+1], [(0) x length($string2)+1]); $DP[0][0] = 0; $DP_dir[0][0] = " "; #initialize the top and left column/row for(my $i=1;$i<=length($string2);$i++){ $DP[0][$i] = -$i; $DP_dir[0][$i] = "I"; } for(my $i=1;$i<=length($string1);$i++){ $DP[$i][0] = -$i; $DP_dir[$i][0] = "D"; } #Fill in the DP table in ROW MAJOR order for(my $i=1;$i<=length($string1);$i++){ for(my $j=1;$j<=length($string2);$j++){ #the value for a deletion my $valU = $DP[$i][$j-1] - 1; #the value for an insertion my $valL = $DP[$i-1][$j] - 1; #value for a (mis)match my $valD = $DP[$i-1][$j-1]; #if the value match add 2, otherwise -1 for a mismatch # notice, we need a shortcircut if the indexies are the beginging, cause no char can match the emptystring if((substr($string1,$i-1,1) eq substr($string2,$j-1,1))){ $valD += 2; }else{ $valD -= 1; } #if the (mis)match has the best value if($valD >= $valU && $valD >= $valL){ $DP[$i][$j] = $valD; $DP_dir[$i][$j] = "M"; #otherwise if the insertion is best }elsif($valU >= $valL){ $DP[$i][$j] = $valU; $DP_dir[$i][$j] = "I"; #finally if nothing else then deltion must be best }else{ $DP[$i][$j] = $valL; $DP_dir[$i][$j] = "D"; } } } foreach my $i(0...length($string1)+1){ foreach my $j(0...length($string2)+1){ print "$DP[$i][$j] "; } print "\n"; } foreach my $i(0...length($string1)+1){ foreach my $j(0...length($string2)+1){ print "$DP_dir[$i][$j] "; } print "\n"; } #traceback go get the actual alignment #start at the bottom corner my $j = length($string2); my $i = length($string1); my @S = ("","",""); while($i>0 || $j>0){ #if its a deletion if($DP_dir[$i][$j] eq "I"){ $S[0] = "-" . $S[0]; $S[1] = " " . $S[1]; $S[2] = substr($string2,$j-1,1) . $S[2]; $j--; #if its an insertion }elsif($DP_dir[$i][$j] eq "D"){ $S[0] = substr($string1,$i-1,1) . $S[0]; $S[1] = " " . $S[1]; $S[2] = "-" . $S[2]; $i--; #match }elsif(substr($string1,$i-1,1) eq substr($string2,$j-1,1)){ $S[0] = substr($string1,$i-1,1) . $S[0]; $S[1] = "|" . $S[1]; $S[2] = substr($string2,$j-1,1) . $S[2]; $i--;$j--; #mismatch }else{ $S[0] = substr($string1,$i-1,1) . $S[0]; $S[1] = " " . $S[1]; $S[2] = substr($string2,$j-1,1) . $S[2]; $i--;$j--; } } print "SCORE: ".$DP[length($string1)][length($string2)]."\n"; print "$S[0]\n$S[1]\n$S[2]\n";