##############################################
############ BLAST PARSING ###################
##############################################




use strict;
use warnings;
use Bio::SearchIO;

#Ferma il processo se gli argomenti allo script sono diversi da 3
die "Usage: $0 <BLAST report file> <number-of-top-hits> <output-file> \n", if (@ARGV!=3 );

my ($infile,$numHits,$outfile) = @ARGV;

print "Parsing BLAST results \n";


my $in = Bio::SearchIO -> new ( -format=>'blast', -file=>$infile   );
								
								
open (OUT, ">$outfile" ) or die "Cannot open $outfile :$!"

#print the header info for tab-delimitated columns
print OUT "query_name\tquery_length\taccession_number\tlength\tdescription\tE value\tbit score\tframe\tquery_start\t";
print OUT "query_end\thit_start\thit_end\tidentical\tHSP_length\tHSP_E value\n";


### Estrazione dei risultati in maniera ricorsiva:

while( my $result = $in -> next_result ) {

			#nome della sequence QUERY:
			print OUT $result-> query_name . "\t";
			
			#lunghezza della QUERY:
			print OUT $result-> query_length;
			
			#"No Hits found" case:
			if ($result->num_hits == 0)  {
			print OUT "\n No hits found \n";
			} else {
			my $count = 0;
			
			
			
			#Processing heach HIT recursively...
			
			while(my $hit = $result -> next_hit) {
			
			print OUTFILE "\t" if ($count > 0);
			 #get the accession number of the hits:
			 print OUTFILE "\t" . $hit->accession . "\t";
							#get the sequence length of the hit:
							print OUT $hit -> length . "\t";
			
			#Get the description of the HIT sequence:
			print OUT $hit -> description . "\t";
			
			#GEt the Evalue of the HIT:
			print OUT $hit -> significance . "\t";
			#Get the score of the HIT:
			print OUT $hit -> bits . "\t";
							
							
			my $hspcount = 0;
				
#Process the top HSP from the top list:

while( my $hsp = $hit -> next_hsp) {

#Metto i separatori perche per raccogliere i dati di ogni singolo HSP (livello più basso) devo dirgli di iniziare scrivere direttamente nella colonna numero 8 lasciando libere le precedenti)
print OUT "\t\t\t\t\t\t\t", if ($hspcount > 0);

		#Get the frame of the query sequence:
		print OUT $hsp->query->frame . "\t";
					
					#Get the start and the end of the QUERY in the alignment:
					print OUT $hsp-> start('query') . "\t" . $hsp->end('query'). "\t";

					#Get the start and the end of the HIT in the alignment:
					print OUT $hsp-> start('hit') . "\t" . $hsp->end('hit'). "\t";
					
					#Get the similarity value:  NOTA BENE: UTILE SOLAMENTE PER BLAST SU SEQUENZE AMINOACIDICHE!!!!
					#printf OUT "%.1f" , ($hsp->frac_conserved * 100);
					#print OUT "%\t";
					
					#Get the similarity value:
					printf OUT "%.1f" , ($hsp->frac_identical * 100);
					print OUT "%\t";
					
					#Get the lenght of the HSP and the E-value:
					print OUT $hsp-> length('total') . "\t";
					print OUT $hsp-> expect . "\t";
					print OUT "\n";   # questo ci dice che deve andare a capo!
					
					
					$hspcount++;
					}
					$count++;
					
					#flow control for the number of hits needed:  @In pratica significa che il secondo argomento dato allo script ($numHits) è il numero di Hit richieste,
					# una volta raggiunto il numero richiesto di Hits per una certa query, lo script passa al 'Result' successivo.
					last if ($count == $numHits);
					
					}
					}
					
}
close OUT;
print "DONE!!!\n";				




#NOTE;
# Un singolo risultato di BLAST è composto da più HIT

# Può accadere che una certa HIT contenga piu di un allineamento: cioè più HSPs

## RESULT-->HIT--->HSPS
### Ogni pacchetto di BIOPERL (Bio::SearchIO) come in questo caso, è caratterizzato da una serie di funzioni proprie che lavorano sugli oggetti creati nel pacchetto stesso

#Ogni oggetto (HSP, Result, HIT) ha le sue proprie funzioni! (vedi www.bioper.org/wiki/HOWTO::SearchIO/)


##USAGE:
# >Perl.exe blast.parsing.pl blast.output.txt 10 blast.output.top10.txt