#! /usr/bin/perl -w

use strict;
use Bio::SeqIO;
use Data::Dumper;

my $over = 1000;

my $file = 'output/alignments.fa';

open(INFO, ">$file");

foreach my $filename (glob("algn/*.fa")) {
	my $gene = $ARGV[0];
	my $result = `bl2seq -p tblastn -i $gene -j $filename`;
	
	my $stream = Bio::SeqIO->new(-file => $filename,
			      -format => 'Fasta');
	my $stream2 = Bio::SeqIO->new(-file => $gene,
			      -format => 'Fasta');

	my @queries = ($result =~ m/Query:.*\n/g);
	my @subjects = ($result =~ m/Sbjct:.*\n/g);
    
	my $seq = $stream->next_seq();
	my $seq2 = $stream2->next_seq();
	

	my $maxQ = -1;
	my $minQ = $seq2->length() + 1;
	my $maxS = -1;
	my $minS = $seq->length() + 1;
   

	#print "\n".$filename;

    for (my $i = 0; $i < @queries; ++$i){
		$queries[$i] =~ m/(\d+)(\D+)(\d+)/;
		my @query = ($1,$2,$3);
		$subjects[$i] =~ m/(\d+)(\D+)(\d+)/;
		my @subject = ($1,$2,$3);
		
		if($query[0] <= $minQ && $subject[0] <= $minS) {
			$minQ = $query[0];
			$minS = $subject[0];
		}
		
		if($query[2] >= $maxQ 
			&& $subject[2] >=  $maxS) {
			$maxQ = $query[2];
			$maxS = $subject[2];
		}
	}

	if($minS - $over > 0) {
		$minS = $minS - $over;
	} else {
		$minS = 1;
	}

	if($maxS + $over <= $seq->length()){
		$maxS = $maxS + $over;
	} else {
		$maxS = $seq->length();
	}

	my $s = $seq->subseq($minS,$maxS);
	#print "\n length".$seq->length();
	#print "\n minQ".$minQ;
	#print "\n minS".$minS;
	#print "\n maxQ".$maxQ;
	#print "\n maxS".$maxS;

	print INFO ">$filename\n$s\n";

	#print  ">$filename\n$result\n";
}

close(INFO);

 
