#!/usr/bin/perl -w
use strict;
use Getopt::Std;

my $TITLE = <<ENDOFTITLE;
#####################################################################
#  SpeedyBee.pl,  v2002-11-22, (c) Sven Rahmann.
#  A (slow) PERL implementation of the FASTA method
#  This script was written for classroom demonstartion only.
#  This is FREE software without any warranty; see the GPL for details.
#####################################################################
ENDOFTITLE

my $USAGE = <<ENDOFUSAGE;
#####################################################################
#  USAGE:          SpeedyBee.pl  [options]  databasefile
#    -h            Help: Print this information; ignore any input
#    -s seqfile    Specify FASTA file with sequences
#                  (default: same as database file)
#    -q qgramlen   Specify q-gram length (default: 2)
#####################################################################
ENDOFUSAGE


#####################################################################

sub fastascores #($query,$q,\%db)
  # Returns hash of scores, one for each sequence in %db
  {
    my $query=shift;    # query sequence (string)
    my $q=shift;        # q-gram length
    my $dbref=shift;    # reference to database sequence hash
    my %scores=();   # hash with return scores

    my ($dbid,$dbseq);
    my ($m,$mq);
    my ($n,$nq);
    my %positions=();  # Hash of position lists indexed by q-grams
    my %dcount=();      # Diagonal counts
    my $sc;
    my ($i,$j,$d);
    my $qgram;

    # Preprocess query
    $m  = length($query);
    $mq = $m-$q+1;
    for ($i=0; $i<$mq; $i++)
      {
        $qgram=substr($query,$i,$q);
        push(@{$positions{$qgram}},$i);
      }

    # Scan each db sequence
    foreach $dbid (keys %$dbref)
      {
        $dbseq = $$dbref{$dbid};
        $n  = length($dbseq);
        $nq = $n-$q+1;
        %dcount=();
        for($j=0; $j<$nq; $j++)
          {
            $qgram=substr($dbseq,$j,$q);
            if (exists($positions{$qgram})) {
              foreach $i (@{$positions{$qgram}}) {
                $d=$j-$i;
                $dcount{$d}++;
              }
            }
          }
        $sc=0;
        foreach $d (keys %dcount)
          { $sc=($dcount{$d}>$sc)?$dcount{$d}:$sc; }
        $scores{$dbid}=$sc;
      }
    return %scores;
  }

#####################################################################

sub readfasta #($filename)
  {
    my $fname=shift;
    my %seqs=();
    my $header = "NIX";
    my $line;
    my $seq = "";

    open FASTA,"<$fname" or die "Cannot open $fname!";
    while (!($header =~ /^>/))  { $header = <FASTA>; }
    chomp($header); $header=~s/^>//;
    while(<FASTA>)
      {
        chomp; $line=$_;
        if ($line =~ /^>/) {
          $seqs{$header}=$seq;
          $header = $line; $header=~s/^>//; $seq="";
        }
        else   {  $seq .= $line;   }
      }
    $seqs{$header}=$seq;
    close FASTA;
    return %seqs;
  }



####################################################################
## MAIN
####################################################################

my %options;
my $sname='';
my $q=2;
my %db;
my %seqs;
my %scores;
my ($key,$k);
my $l;

$TITLE =~ s/\#(  )?//g;
$USAGE =~ s/\#(  )?//g;
$USAGE =~ s/\n//;
print STDERR $TITLE;

getopts("hs:q:", \%options);
if ($options{h} || (!@ARGV)) { print STDERR $USAGE; exit;}
if ($options{s}) { $sname=$options{s}; }
else { $sname=$ARGV[0]; }
if ($options{q}) { $q=$options{q}; }

%db=readfasta($ARGV[0]);
%seqs=readfasta($sname);

foreach $key (sort keys %seqs)
  {
    $l=length($seqs{$key});
    %scores = fastascores($seqs{$key},$q,\%db);
    printf "FASTA-Scores for Sequence '$key' (length $l):\n";
    foreach $k (sort keys %scores)
      { print "   $k => $scores{$k}\n"; }
    printf "\n";
  }
####################################################################