#### Analyse a BLAST file

# Expected input:
#
#       (header)
#
#       Query= [query_name] [query_description]
#                 ([query_length] letters)
#
#       (garbage)
#
#       >ContigMMM.revised.geneNNN.protein [subject1_description]
#                Length = [subject1_length]
#
#        Score = [garbage], Expect = [expectation1]
#
#       (garbage)
#
#       >ContigMMM.revised.geneNNN.protein [subject2_description]
#                Length = [subject2_length]
#
#        Score = [discard], Expect = [expectation2]
#
#        (garbage)
#
# NOTE: The Query= line and >Contig line may spill over if the description
#       is too long
#
# NOTE: The >Contig line gives the Subject name: MMM.NNN 
# 
# Output:
#
# query_name query_description query_length subject1_name 
#            subject1_description subject1_length    
#
# NOTE: If there is more than one matching subject, then 
#       subject2_name subject2_description subject2_length is appended 
#            to the end of the line.
#	IF there is no matching subject, then the line ends with query_length
#
# NOTE: There is a line for every query submitted to Blast
#
# NOTE: Each item is separated from the next by a TAB

#### OPEN THE BLAST FILE

my ($blastPath) = '71vsnps.txt';
open BLAST, "<$blastPath" or die "Can't open $blastPath: $!\n";

#### OPEN OUR OUTPUT FILE

my ($resultPath) = '71vsnps-out.tab';
open RESULT, ">$resultPath" or die "Can't open $resultPath: $!\n";

############################## MAIN PROGRAM
#
# We'll accumlate the items of information for each query in the
# array @query_info
#
# Then, for each line of the BLAST file, we check for two things: the start
# of a query, and the start of a match.
#
# When we find the start of a query, we print the accumulated information for
# the previous query. (But when we see the first query there won't be a previous
# query.  What then?  See below at PRINTING THE PREVIOUS QUERY.)
#
# The other possibility is that we might see the start of a match.  If so,
# we record it in @query_info.

my @query_info;

while (<BLAST>) {
   if (/^Query=/) {
      print_previous_query();
      start_new_query();
   }
   if (/^>Contig/) { record_subject 	() }
}

# Now for the last query.  All the others were printed when we were
# about to process the query after them; but there is no query after the
# last query, so we make an extra call to print_previous_query().

print_previous_query();

############# PRINTING THE PREVIOUS QUERY
#
# One quirk here:
#    o We might not have a previous query, so we need to check before printing.

sub print_previous_query {
   if (@query_info) {
      my $result_line = join("\t", @query_info);
      print $result_line, "\n";
      print RESULT $result_line, "\n";
      @query_info = ();
   }
}

sub start_new_query {
   my ($query, $query_description) = /^Query=\s+(\w+)\W+(.+)/;
   unless (defined $query) {
      warn "Didn't find query name on line $.\n";
      return;
   }
   while (<BLAST>) {
      my ($query_length) = /(\d+) letters/;
      if (defined $query_length) {
         @query_info = ($query, $query_description, $query_length);
         return;
      }
   }
}

my ($subject_name, $subject_description, $subject_length, $expectation);

sub record_subject {
   
   find_subject_name();
   unless (defined $subject_name) {
      warn "Didn't find subject_name on line $.\n";
      return;
   }

   find_subject_length();
   unless (defined $subject_length) {
      warn "Didn't find subject_length on line $.\n";
      return;
   }

   find_expectation();
   unless (defined $expectation) {
      warn "Didn't find expectation of chance match on line $.\n";
      return;
   }

   push @query_info, $subject_name, $subject_description, $subject_length,
        $expectation;

}

sub find_subject_name {
   $subject_name = undef;
   if (/^>Contig(\d+)\.revised\.gene(\d+)\.protein\W+(.+)/) { 
     $subject_name = sprintf "%03d.%03d", $1, $2;
     $subject_description = $3;
   }
}

sub find_subject_length {
   $subject_length = undef;
   while (<BLAST>) {
      ($subject_length) = /Length = (\d+)/;
      if (defined $subject_length) {
         return;
      }
   }
}

sub find_expectation {
   $expectation = undef;
   $_ = <BLAST>;
   if (defined $_) { $_ = <BLAST> }
   unless (defined $_) {
      warn "Incomplete query at $.\n";
      return;
   }
   ($expectation) = /Expect = (\S+)/;
}


