#! /usr/bin/perl -w

##### Program description #######
#
# Title:     bp_ex5.pl
# Author(s): Mattias Ohlsson
# Description:
#   A Perl program that can run a remote blast at NCBI. The 
#   module Tools::Run::RemoteBlast is used.
#
##################################

use strict;
use Bio::SeqIO;
use Bio::Tools::Run::RemoteBlast;

# Get the query sequence
my $Seq_in = Bio::SeqIO->new (-file => '1TEN.fasta', 
			      -format => 'fasta');
my $query = $Seq_in->next_seq(); 


# Remote-blast "factory object" creation
my $factory = Bio::Tools::Run::RemoteBlast->new(
						-prog       => 'blastp',
						-data       => 'nr',
						-expect     => '1e-6',
						-readmethod => 'SearchIO');

# Here we submit the query sequence to the Blast server
$factory->submit_blast($query);

print STDERR "waiting...";

# Loop over all possible remote id's
while ( my @rids = $factory->each_rid ) {
    # printf "Debug: No rids %d\n", scalar @rids;

   foreach my $rid ( @rids ) {
      my $rc = $factory->retrieve_blast($rid);
      
      unless ( ref $rc ) {
	 
	 # retrieve_blast returns -1 on error
	 if ( $rc < 0 ) {
	    $factory->remove_rid($rid);
	 }
	 
	 # retrieve_blast returns 0 on 'job not finished'
	 sleep 5;
	 print STDERR ".";
	 
      } else {
	 
	 # Blast done
	 $factory->remove_rid($rid);
	 my $result = $rc->next_result();
	 
	 # Save the output
	 my $filename = 'blast.out';
	 $factory->save_output($filename);
	 
	 # Analyse the blast result
	 print "Database: ", $result->database_name(), "\n";
	 print "\nQuery Name: ", $result->query_name(), "\n\n";
	 while ( my $hit = $result->next_hit ) {
	    while ( my $hsp = $hit->next_hsp ) {
	       if ( $hsp->percent_identity >= 95 && $hsp->percent_identity <= 99 ) {
		  print "======= ID: ", $hit->name, "=======\n";
		  printf "=> Identities: %5.2f%%\n\n", $hsp->percent_identity;
		  
		  my $q = $hsp->query_string; 
		  my @q = ($q =~ /(.{1,60})/g);
		  my $sQ = $hsp->start('query');
		  
		  my $s = $hsp->hit_string;
		  my @s = ($s =~ /(.{1,60})/g);
		  my $sS = $hsp->start('subject');
		  
		  my $hom = $hsp->homology_string;
		  my @hom = ($hom =~ /(.{1,60})/g);
		  
		  for ( my $i = 0; $i < scalar @q; $i++ ) {
		     printf "Query: %4d $q[$i] %4d\n", $sQ, $sQ+length($q[$i])-1;
		     printf "       %4s $hom[$i]\n", '';
		     printf "Subj:  %4d $s[$i] %4d\n\n", $sS, $sS+length($s[$i])-1;
		  }
	       }
	    }
	 }
      }
   }
}
