#!/usr/bin/perl -w
use strict;

# BlastN: Basic local alignment sequence tool: nucleotide-to-nucleotide

# Constants
my $true = 1;
my $false = ! $true;
my $gap = "-";
my $bar = "|";

############ Initialization
#
my ($query_name, $query) = get_file_info('DG47.nt');
my ($target_name, $target) = get_file_info('M29081.nt');

print "Query:  $query_name\n",
      "Target: $target_name\n\n";

my $word_length;
set_blast_parameters();

########### Main program
#
for (my $query_word_start = 0;
        $query_word_start <= length($query) - $word_length;
        $query_word_start = $query_word_start+1) {
   my $query_word = substr($query, $query_word_start, $word_length);
   $query_word = qr/$query_word/;
   while ($target =~ /$query_word/g) {
      my $target_word_start = pos($target) - $word_length;
      one_match($target_word_start, $query_word_start);
      pos($target) = $target_word_start+1;
   }
}

######## Initialization subroutines
#
sub get_file_info {
   my ($path) = @_;
   open INPUT, "<$path" or die "Can't open $path: $!\n";
   my $name = <INPUT>;
   chomp $name;
   my $sequence = "";
   while (<INPUT>) {
      chomp;
      $sequence .= uc $_;
   }
   close INPUT;
   return ($name, $sequence);
}

my ($match_reward, $mismatch_penalty, $gap_begin_penalty, $gap_extend_penalty);
sub set_blast_parameters {
   $word_length = 11;
   $match_reward = 1;
   $mismatch_penalty = 2;
   $gap_begin_penalty = 5;
   $gap_extend_penalty = 2;
}

######### Search and print control subroutines
#

sub one_match {
   my ($target_word_start, $query_word_start) = @_;
   my $target_word_end = $target_word_start + $word_length - 1;
   my $query_word_end = $query_word_start + $word_length - 1;
   
   my ($target_max, $query_max, $forward_target, $forward_query)
      = half_match(1, $target_word_end, $query_word_end);

   my ($target_min, $query_min, $backward_target, $backward_query)
      = half_match(-1, $target_word_start, $query_word_start);

   my $target_word = substr($target, $target_word_start, $word_length);
   my $target_hit = reverse($backward_target).$target_word.$forward_target;

   my $query_word = substr($query, $query_word_start, $word_length);
   my $query_hit = reverse($backward_query).$query_word.$forward_query;

   print_hit($target_min, $target_hit, $query_min, $query_hit);

}

# Returns:
#   last index in target of target extension sequence,
#   last index in query of query extension sequence.
#   target extension sequence,
#   query extension sequence,
#
my ($direction, $t, $q);
my ($max_score, $max_i, $max_j, @score, @last_i, @last_j);
sub half_match {
   ($direction, $t, $q) = @_;
   @last_i = ();
   $last_i[0][0] = -1;
   @last_j = ();
   $last_j[0][0] = -1; 
   $max_i = 0;
   $max_j = 0;
   @score = ();
   $score[0][0] = $max_score = $match_reward * $word_length;
   fill_right();
   fill_down();
   my ($xt, $xq) = extensions();
   return (
      $t + $max_i * $direction,
      $q + $max_j * $direction,
      $xt,
      $xq
   );
}

sub extensions {
   my (@xt, @xq);
   my $i = $max_i;
   my $j = $max_j;
   while ($i != 0) {
      my $i0 = $last_i[$i][$j];
      my $j0 = $last_j[$i][$j];
      while ($i > $i0 or $j > $j0) {
         if ($i == $i0) { unshift @xt, $gap }
         else {
            unshift @xt, substr($target, $t + $i*$direction, 1);
            $i = $i-1;
         }
         if ($j == $j0) { unshift @xq, $gap }
         else {
            unshift @xq, substr($query, $q + $j*$direction, 1);
            $j = $j-1;
         }
      }
   }
   return (join("", @xt), join("", @xq));
}

#### Printing subroutines
#

my %seen;
sub print_hit {
   my($target_min, $target_hit, $query_min, $query_hit) = @_;
   $target_min = $target_min+1; # Perl counts from 0
   $query_min = $query_min+1;   # Perl counts from 0

   my $key = "$target_min:$target_hit:$query_min:$query_hit";
   return if $seen{$key};
   $seen{$key} = $true;

   my $start = 0;
   my $max_segment_length = 50;
   while ($start < length($target_hit)) {
      my $segment_length
         = min($max_segment_length, length($target_hit) - $start);
      my $query_segment = substr($query_hit, $start, $segment_length);
      my $target_segment = substr($target_hit, $start, $segment_length);
      my $query_length = count_bases($query_segment);
      my $target_length = count_bases($target_segment);

      print_segment("Query", $query_segment, $query_min, $query_length);
      print_alignment($target_segment, $query_segment);
      print_segment("Target", $target_segment, $target_min, $target_length);
      print "\n";

      $target_min = $target_min + $target_length;
      $query_min = $query_min + $query_length;
      $start = $start + $segment_length;
   }
}

sub count_bases {
   my ($segment) = @_;
   $segment =~ s/$gap//go;
   return length($segment);
}

sub print_segment {
   my ($name, $segment, $start, $length) = @_;
   printf "%-7s %7d %s %7d\n", "$name:", $start, $segment, $start+$length-1;
}

sub print_alignment {
   my ($s_string, $t_string) = @_;
   my @s = split("", $s_string);
   my @t = split("", $t_string);
   my @bars;
   for (my $k = 0; $k < @s; $k = $k+1) {
      if ($s[$k] eq $t[$k] and $s[$k] ne $gap) { push @bars, $bar }
      else { push @bars, " " }
   }
   print " " x 16, @bars, "\n";
}

sub max {
   my ($a, $b) = @_;
   if ($a > $b) { return $a }
   else { return $b }
}

sub min {
   my ($a, $b) = @_;
   if ($a < $b) { return $a }
   else { return $b }
}

################ Searching subroutines
#

sub fill_right {
   my $previous_changed = do_gap_right(0, 0);
   my $i = 1;
   while ($i <= $previous_changed+1) {
      my $this_changed = -1;
      for (my $j = $i; $j <= $previous_changed+1; $j = $j+1) {
         my $changed_score = check_diagonal($i, $j);
         if ($changed_score) {
             $this_changed = max($this_changed, do_gap_right($i, $j));
             do_gap_down($i, $j);
         }
         if ($j <= $previous_changed) {
            # Check vertical
            $changed_score = check_score($i-1, $j, $i, $j, -$gap_begin_penalty);
            if ($changed_score) {
               do_gap_down($i-1, $j);
               $this_changed = max($this_changed, do_gap_right($i, $j));
            }
         }
      }
      $i = $i+1;
      $previous_changed = $this_changed;
   }
}

sub check_diagonal {
   my ($i, $j) = @_;

   my $target_pos = $t + $i*$direction;
   return $false unless 0 <= $target_pos and $target_pos < length($target);

   my $query_pos = $q + $j*$direction;
   return $false unless 0 <= $query_pos and $query_pos < length($query);

   my $change_amount;
   if (substr($target, $target_pos, 1) eq substr($query, $query_pos, 1))
      { $change_amount = $match_reward }
   else { $change_amount = -$mismatch_penalty }

   return check_score($i-1, $j-1, $i, $j, $change_amount);
}

sub fill_down {
   my $previous_changed = do_gap_down(0, 0);
   my $j = 1;
   while ($j <= $previous_changed+1) {
      my $this_changed = -1;
      for (my $i = $j; $i <= $previous_changed+1; $i = $i+1) {
        my $changed_score = check_diagonal($i, $j);
        if ($changed_score) {
            $this_changed = max($this_changed, do_gap_down($i, $j));
            do_gap_right($i, $j);
         }
         if ($i <= $previous_changed) {
            # Check horizontal
            $changed_score = check_score($i, $j-1, $i, $j, -$gap_begin_penalty);
            if ($changed_score) {
               do_gap_right($i, $j-1);
               $this_changed = max($this_changed, do_gap_right($i, $j));
            }
         }
      }
      $j = $j+1;
      $previous_changed = $this_changed;
   }
}

sub check_score {
   my ($from_i, $from_j, $to_i, $to_j, $change_amount) = @_;

   my $old_score = $score[$to_i][$to_j] || 0;

   my $new_score = $score[$from_i][$from_j];
   return $false unless defined($new_score);
   $new_score = $new_score + $change_amount;

   return $false unless $new_score > $old_score;

   record_score_change($from_i, $from_j, $to_i, $to_j, $new_score);
   return $true;
}

sub record_score_change {
   my ($from_i, $from_j, $to_i, $to_j, $new_score) = @_;

   $score[$to_i][$to_j] = $new_score;
   $last_i[$to_i][$to_j] = $from_i;
   $last_j[$to_i][$to_j] = $from_j;

   if ($new_score > $max_score) {
      $max_score = $new_score;
      $max_i = $to_i;
      $max_j = $to_j;
   }

}

sub do_gap_down {
   my ($i, $j) = @_;
   my $new_score = $score[$i][$j] || 0;
   my $penalty = $gap_begin_penalty;
   my $changed = $i;
   while ($true) {
      $new_score = $new_score - $penalty;
      $penalty = $gap_extend_penalty;
      $i = $i+1;
      my $old_score = $score[$i][$j] || 0;
      if ($new_score > $old_score) {
         record_score_change($i-1, $j, $i, $j, $new_score);
         $changed = $i;
      }
      elsif ($new_score-$gap_extend_penalty <= $old_score-$gap_begin_penalty) {
         # We have no further chance of exceeding old scores
         return $changed;
      }
   }
}

sub do_gap_right {
   my ($i, $j) = @_;
   my $new_score = $score[$i][$j] || 00;
   my $penalty = $gap_begin_penalty;
   my $changed = $j;
   while ($true) {
      $new_score = $new_score - $penalty;
      $penalty = $gap_extend_penalty;
      $j = $j+1;
      my $old_score = $score[$i][$j] || 0;
      unless (defined($old_score)) { $old_score = 0 }
      if ($new_score > $old_score) {
         record_score_change($i, $j-1, $i, $j, $new_score);
         $changed = $j;
      }
      elsif ($new_score-$gap_extend_penalty <= $old_score-$gap_begin_penalty) {
         # We have no further chance of exceeding old scores
         return $changed;
      }
   }
}
