#!/usr/bin/perl -w
####### Program description #######
#
# Title: parse.pl
#
# Author(s): Mattias Ohlsson
#
# Description: 
#  This is a parser for result files produced by the Blastp program
#  for sequence alignment.
#
# List of subroutines:
#  GetSubjects(); 
#  GetHSPs();
#  AnalyzeHSP();
#
# Overall procedure:
#  This program loads a result file from a blastp run. The program does
#  not check that the loaded file is a valid blastp file, only that it
#  is a non-empty file. The file is stored as a single scalar. All the
#  analysis is performed on scalar and not on an array containing lines
#  of the result file. The first step is to split the result into
#  subjects, this is performed in the GetSubjects routine. The program
#  then makes a loop over all found subjects and for each subject all
#  the high scoring pairs are extracted (GetHSP). A second loop is over
#  all HSPs found and where each HSP is analyzed in the AnalyzeHSP
#  routine.  This routine extracts all the information needed in order
#  to select the ones with a percentage of identities within a specified
#  range. The selected ones are printed together with the alignment
#  (optional).
#
# Usage:
#  ./parse.pl {blast result file}
#
###################################
use strict;

#### The main program ####

### Part 1:  Read the input file and check that is non-zero
undef $/;
my $Blast = (<>);
unless ( length $Blast > 0 ) {
    die "Zero length input file\n";
}

### Part 2: Some hardcoded constants
my $Icutlow = 80;
my $Icuthigh = 85;
my $PrintAlign = 1;

### Part 3: The rest
my @subjects = GetSubjects($Blast);

my (@Idents, @Alns);
foreach my $subj ( @subjects ) {

    my ($sID, $rHSPs) = GetHSPs($subj); # All the high scoring pairs

    undef @Idents; 
    undef @Alns;
    for (my $i = 0; $i < @{$rHSPs}; $i++) {
	my $rHPSres = AnalyzeHSP($rHSPs->[$i]);
	my $id = $rHPSres->{'Pid'}; 

	# Check the criteria
	if ( $id >= $Icutlow && $id <= $Icuthigh ) {
	    push @Idents, $id;
	    push @Alns, $rHPSres->{'Aln'};
	}
    } 

    if ( @Idents ) {
	print "======= ID: $sID  =======\n";
	foreach (my $i = 0; $i < @Idents; $i++) {
	    print "=> Identities: $Idents[$i]%\n";
	    if ( $PrintAlign == 1 ) { print $Alns[$i], "\n"; }
	}
    }
}


sub GetSubjects {
    # The blast "file"
    my ($blast) = @_; # Note this is a copy of the supplied blast "file"!

    # Get rid of everything before the ALIGNMENTS line
    $blast =~ s/^.*ALIGNMENTS\n//s;

    # and everthing after the database statement
    $blast =~ s/\s\sDatabase:.*//s;

    # Now split on the >xxx pattern but keep the >xxx itself
    my @subjects = split /(?=>\w{2,3})/, $blast;

    return @subjects;

} # End of GetSubjects


sub GetHSPs {
    # One subject
    my ($subject) = @_; # Note this is a copy of the supplied subject!

    # Extract an ID
    my $id;
    if ( $subject =~ /^>(\w{2,3}\|.*?\|)/ ) {
	$id = $1;
    } else {
	print "Warning: Subject does not start with a known format\n";
	print $subject;
	exit;
    }

    # Now get HSPs by dividing at the Score line
    $subject =~ s/^.*?(?=Score)//s;
    my @hsps =  split /(?=\sScore)/, $subject;

    return($id, \@hsps);

} # End of GetHSPs


sub AnalyzeHSP {
    # The HSP
    my ($hsp) = @_; # Note this is a copy of the supplied hsp!
    
    # This will store the results
    my %HSPres; 
 #Identities = 90/90 (100%), Positives = 90/90 (100%), Gaps = 0/90 (0%)
    # Find the score
    if ( $hsp =~ /Identities\s+=\s+\d+\/\d+\s+\((\d+)%\)/ ) {
	$HSPres{'Pid'} = $1;
    }

    # Get the Alignment
    $hsp =~ s/^.*?(?=Query)//s;
    $hsp =~ s/\n(?=\n)//g;
    $HSPres{'Aln'} = $hsp;
    
    # Now return a reference to the result hash
    return \%HSPres;

} # End of AnalyzeHSP
