#!/usr/bin/perl

use CGI;
use LWP::Simple;
use Bio::Perl;
use strict;

$ENV{'PATH'} = '/var/www/Genomics/Wideloache/WebFiles/blast/bin';

# Standard creation of interface with the input and results web page
my $q = new CGI;
    print $q->header( -type => "text/html", -expires => "now" ), $q->start_html( -title => "EC BLASTer Results" );

open(STDERR, ">&STDOUT");

my $ecNum = $q->param("ec"); 

print "<h2>Blast Results for EC Number: $ecNum</h2><hr>";

my @accNum = getAccNums($ecNum);
my $aaFile = "aaFiles/" . $ecNum . "aa.faa";   
my $outfilename = "outfiles/output" . $ecNum . ".txt";
if (scalar(@accNum) == 0) {
   print "No protein sequences found for EC Number $ecNum.\n";
   exit;
} 
elsif (-e $outfilename) {
print "Blast results from a previous run were found.<br><br>";
#Skip the file writing because it has already been done.
}
else {
   my $output = "";
   my $output2 = "";
   my @seqObjects;
   for (my $i = 0; $i < scalar(@accNum); $i++) { 
      my $seq_object = get_sequence('swiss', $accNum[$i]); 
      push (@seqObjects, $seq_object);     
   }
   write_sequence(">$aaFile", 'fasta', @seqObjects);  
   $output = `blastall -p blastp -i $aaFile -d blast/data/haloAll.faa -e .03 -m 9`; 
   $output2 = `blastall -p blastp -i $aaFile -d blast/data/haloAll.faa -e .03 -m 0`;   
   open (OUT, "> outfiles/output" . $ecNum . ".txt");
   print OUT "$output\n";
   close OUT;
   open (OUT, "> outfiles/output2" . $ecNum . ".txt");
   print OUT "$output2\n";
   close OUT;
}

my $outputRead = "";

open(INFILE, "< outfiles/output2" . $ecNum . ".txt" || die "Error: Output file not found!");
while(my $line = <INFILE>) {
   $outputRead = $outputRead . $line;
}
close INFILE;

$outputRead =~ s/\n//g;
print scalar(@accNum) . " known protein sequences with E.C. number " . $ecNum . " were found and blasted against the <i>Halorhabdus Utahensis</i> protein calls.<br><br>";
if ($outputRead eq "") {
   print "No significant hits were found.";
}
else {
   print "You can find your blast results <a href=\"outfiles/output" . $ecNum . ".txt\"> here </a>. All results with an e-value of <.03 are shown.<br>To see more detailed blast alignments (the full blast output for each protein) click <a href=\"outfiles/output2" . $ecNum . ".txt\"> here </a>.";
}

print end_html;

### Returns an array of swissprot accession numbers for the given EC number
sub getAccNums {
   my $ecNum = shift(@_);
   my $url = "http://expasy.org/enzyme/" . $ecNum;
   my $html = get($url)
     or die "Couldn't fetch page.";
     
   my @accNum;
   
   @accNum = $html =~ m#uniprot/......#g;
   for (my $i = 0; $i < scalar(@accNum); $i++) { 
      $accNum[$i] = substr($accNum[$i],8);
   }
   return @accNum;
}