#!/usr/local/bin/perl -w
use strict;

# MicroSat: Detects repeated sequences in DNA
#           
# Uses autocorrelation, i.e. matches between bases spaced n units
# apart

################# VARIABLES ###################

  my @matches;                  # @matches(n) contains consecutive matches
                                #   at a distance of n
  my @sequence;                 # Sequence, one base per array element.
  my $trials;                   # number of matches attempted; always
                                #   sequence length - (repeat unit)
  my $unit_min = 2;             # smallest repeat unit
  my $unit_max = 30;            # largest repeat unit


################### MAIN PROGRAM ###############

@sequence = split("", get_genome_sequence('Reptest1.nt'));

Find_matches();
Print_analysis();

################ SUBROUTINES ###################

sub Find_matches {
   for (my $distance = $unit_min; $distance <= $unit_max;
           $distance = $distance+1) {
      $matches[$distance] = Distant_pairs($distance);
   }
}

sub Distant_pairs {
   my ($distance) = @_;
   for (my $j = $distance; $j < @sequence; $j = $j+1) {
      if ($sequence[$j] eq $sequence[$j-$distance]) { 
         $matches[$distance] = $matches[$distance] + 1
      }
   }
}

sub Print_analysis {
  print "Unit  Matches  Trials  Pct\n";
  foreach my $distance ($unit_min..$unit_max) {
    $trials = @sequence - $distance;
    printf "%4d  %7d  %6d  %7.1f\n", $distance, $matches[$distance], $trials,
            $matches[$distance]/$trials*100;
  }
}

sub get_genome_sequence {
   my ($file_path) = @_;
   open GENOME, "<$file_path" or die "Can't open $file_path: $!\n";
   $_ = <GENOME>; # discard first line
   my $sequence = "";
   while (<GENOME>) {
      chomp($_);                  # Remove line break
      s/\r//;                     # And carriage return
      $_ = uc $_;                 # Make the line uppercase;
      $sequence = $sequence . $_; # Append the line to the sequence
   }
   return $sequence;
}
