#!/usr/local/bin/perl

use vars qw($opt_e $opt_h);
use strict;
use Getopt::Std;

getopts('e:h');

help() if (!$opt_e || $opt_h);

#-------------------------------------------------------------#
# Main
my %sub_seq;
my $query;
my $flag;
my $seq;
my $read_on;

while(<>) {

	
	if (!$read_on) {
	
		$read_on = 1 if (/^CONVERGED/);
		next;
	}
	
	if (/^>[^\|]+\|[^\|]+\|(\S+)/) {
	
		if ($flag) {

			print_seq("$query\_$sub_seq{$query}",$seq);
		}
		
		$flag = undef;
		$seq = undef;
		
		$sub_seq{$1}++;
		$query = $1;
		next;
	}
	if (/Expect\s*=\s*(\S+)/) {

		if ($seq && $flag) {

			print_seq("$query\_$sub_seq{$query}",$seq);
			$sub_seq{$query}++;
			$flag = undef;
			$seq = undef;
		}
		
		$flag = (scalar($1) <= scalar($opt_e))? 1:undef;
		next;
	}
	if (/^\s*Sbjct:\s*\d+\s+(\S+)/) {

		$seq .= $1;
	}
}

#-------------------------------------------------------------#
# Print sequence in fasta format
sub print_seq() {

	my ($q,$s) = @_;
	
	$s =~ s/[-]//g;
		
	print ">$q\n";
	
	while (length($s)>0) {
		
		my $out = substr($s,0,60);
		substr($s,0,60) = '';
		print "$out\n";
	}
	print "\n";
}

#-------------------------------------------------------------#
# Help
sub help {

print <<EOF;

extract_psi_seq.pl [options] psi_blast_file

where options are:
	
	-e <float> Threshold E-val (recover sequences with e-val <= this value)
	-h         This help

EOF
exit(0);
}
