设为首页收藏本站
开启辅助访问
切换到窄版

 找回密码
 注册

QQ登录

只需一步,快速开始

搜索
查看: 1896|回复: 0

Smith-Waterman算法的Perl实现

[复制链接]
发表于 2011-10-6 11:05:08 | 显示全部楼层 |阅读模式

  1. #!/usr/bin/perl
  2. # Smith-Waterman  Algorithm
  3. # usage statement
  4. die "usage: {row.content} <sequence 1> <sequence 2>\n" unless @ARGV == 2;
  5. # get sequences from command line
  6. my ($seq1, $seq2) = @ARGV;
  7. # scoring scheme
  8. my $MATCH    =  1; # +1 for letters that match
  9. my $MISMATCH = -1; # -1 for letters that mismatch
  10. my $GAP      = -1; # -1 for any gap
  11. # initialization
  12. my @matrix;
  13. $matrix[0][0]{score}   = 0;
  14. $matrix[0][0]{pointer} = "none";
  15. for(my $j = 1; $j <= length($seq1); $j++) {
  16.    $matrix[0][$j]{score}   = 0;
  17.    $matrix[0][$j]{pointer} = "none";
  18. }
  19. for (my $i = 1; $i <= length($seq2); $i++) {
  20.    $matrix[$i][0]{score}   = 0;
  21.    $matrix[$i][0]{pointer} = "none";
  22. }
  23. # fill
  24. my $max_i     = 0;
  25. my $max_j     = 0;
  26. my $max_score = 0;
  27. for(my $i = 1; $i <= length($seq2); $i++) {
  28.    for(my $j = 1; $j <= length($seq1); $j++) {
  29.        my ($diagonal_score, $left_score, $up_score);
  30.       
  31.        # calculate match score
  32.        my $letter1 = substr($seq1, $j-1, 1);
  33.        my $letter2 = substr($seq2, $i-1, 1);      
  34.        if ($letter1 eq $letter2) {
  35.            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MATCH;
  36.        }
  37.        else {
  38.            $diagonal_score = $matrix[$i-1][$j-1]{score} + $MISMATCH;
  39.        }
  40.       
  41.        # calculate gap scores
  42.        $up_score   = $matrix[$i-1][$j]{score} + $GAP;
  43.        $left_score = $matrix[$i][$j-1]{score} + $GAP;
  44.       
  45.        if ($diagonal_score <= 0 and $up_score <= 0 and $left_score <= 0) {
  46.            $matrix[$i][$j]{score}   = 0;
  47.            $matrix[$i][$j]{pointer} = "none";
  48.            next; # terminate this iteration of the loop
  49.        }
  50.       
  51.        # choose best score
  52.        if ($diagonal_score >= $up_score) {
  53.            if ($diagonal_score >= $left_score) {
  54.                $matrix[$i][$j]{score}   = $diagonal_score;
  55.                $matrix[$i][$j]{pointer} = "diagonal";
  56.            }
  57.            else {
  58.                $matrix[$i][$j]{score}   = $left_score;
  59.                $matrix[$i][$j]{pointer} = "left";
  60.            }
  61.        } else {
  62.            if ($up_score >= $left_score) {
  63.                $matrix[$i][$j]{score}   = $up_score;
  64.                $matrix[$i][$j]{pointer} = "up";
  65.            }
  66.            else {
  67.                $matrix[$i][$j]{score}   = $left_score;
  68.                $matrix[$i][$j]{pointer} = "left";
  69.            }
  70.        }
  71.       
  72.        # set maximum score
  73.        if ($matrix[$i][$j]{score} > $max_score) {
  74.            $max_i     = $i;
  75.            $max_j     = $j;
  76.            $max_score = $matrix[$i][$j]{score};
  77.        }
  78.    }
  79. }
  80. # trace-back
  81. my $align1 = "";
  82. my $align2 = "";
  83. my $j = $max_j;
  84. my $i = $max_i;
  85. while (1) {
  86.    last if $matrix[$i][$j]{pointer} eq "none";
  87.    
  88.    if ($matrix[$i][$j]{pointer} eq "diagonal") {
  89.        $align1 .= substr($seq1, $j-1, 1);
  90.        $align2 .= substr($seq2, $i-1, 1);
  91.        $i--; $j--;
  92.    }
  93.    elsif ($matrix[$i][$j]{pointer} eq "left") {
  94.        $align1 .= substr($seq1, $j-1, 1);
  95.        $align2 .= "-";
  96.        $j--;
  97.    }
  98.    elsif ($matrix[$i][$j]{pointer} eq "up") {
  99.        $align1 .= "-";
  100.        $align2 .= substr($seq2, $i-1, 1);
  101.        $i--;
  102.    }   
  103. }
  104. $align1 = reverse $align1;
  105. $align2 = reverse $align2;
  106. print "$align1\n";
  107. print "$align2\n";
复制代码
您需要登录后才可以回帖 登录 | 注册

本版积分规则

QQ|申请友链|小黑屋|手机版|Archiver|生物信息学论坛 ( 蜀ICP备09031721号  

GMT+8, 2017-2-21 08:49 , Processed in 0.102392 second(s), 21 queries .

Powered by Discuz! X3

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表