#!/usr/local/bin/perl -w

use strict;
use diagnostics;
use warnings;
use Getopt::Long;

###############################################################################
#Tiling
###############################################################################
use Bio::SearchIO;
use Bio::Search::Tiling::MapTiling;
use Bio::Tools::GFF;

my $input = $ARGV[0];
my $type  = $ARGV[1];

my $out   = new Bio::Tools::GFF( -gff_version => 3 );

my $blio  = Bio::SearchIO->new(
	-format => 'blast',
	-file   => $input
);

my $other = $type eq 'query' ? 'hit' : 'query';

my %seen;

while ( my $result = $blio->next_result ) {
	my $countHit = 0;
	while ( my $hit = $result->next_hit ) {
		$countHit++;

		#Tiling
		my $tiling = Bio::Search::Tiling::MapTiling->new($hit);
		my @contexts = $tiling->contexts($type);

		#TEST
		my $bestContext = pop @contexts;
		my $bestIdentities = $tiling->identities($type, 'exact', $bestContext);
		my @bestContexts = ($bestContext);
		for my $context ( @contexts ) {
			my $identities = $tiling->identities($type, 'exact', $context);
			#print "current identities=$identities in context=$context for type=$type\n";
			if ($identities > $bestIdentities) {
				$bestIdentities = $identities;
				@bestContexts = ($context);
			}
			elsif($identities == $bestIdentities){
				push(@bestContexts, $context);
			}
		}
		#print "best identities=$bestIdentities in context=$bestContext for type=$type\n";
		@contexts = @bestContexts;

		my $countTypeContext = 0;
		for my $typeContext ( @contexts ) {
			$countTypeContext++;
    		my @hspTiled = $tiling->next_tiling($type, $typeContext);

    		my @otherContexts = $tiling->contexts($other);
			my $countOtherContext = 0;
    		for my $otherContext ( @otherContexts ) {
    			$countOtherContext++;
    			#print "otherContext=$otherContext for type=$other\n";
    			my $fracAligned     = $tiling->frac_aligned    (-type => $other, -action => 'exact', -context => $otherContext, -denom => 'total') * 100;
   				my $percentIdentity = $tiling->percent_identity(-type => $other, -action => 'exact', -context => $otherContext);

   				#next if ($fracAligned < 50);
				#next if ($percentIdentity < 50);

    			my ($hspFeatTab, $min, $max, $st, $seqid, $name) = getHSPFeatTab(\@hspTiled, $hit, $result, "$countHit.$countTypeContext.$countOtherContext");
    			print "#####################TILED#####################\n";
   				makeGFF($hit, $hspFeatTab, $st, $min, $max, $type, $name, "$countHit.$countTypeContext.$countOtherContext", $seqid, $percentIdentity, $fracAligned);
    		}
		}
		#end tiling

		#test old behaviour
		my @hsps = $hit->hsps();
		my ($hspFeatTab, $min, $max, $st, $seqid, $name) = getHSPFeatTab(\@hsps, $hit, $result, $countHit);
    	print "***************UNTILED*****************\n";
    	makeGFF($hit, $hspFeatTab, $st, $min, $max, $type, $name, $countHit, $seqid, "percentIdentity", "fracAligned");

	}#boucle Hits
}#boucle Results

sub getHSPFeatTab {

	my ($hspTab, $hit, $result, $countHit) = @_;
	my ($name, %min, %max, @hspFeatTab, $st, $seqid);
	my @hspTabFeatures;
	for my $hsp ( @$hspTab ) {
		#make HSP feature
		my $feature = new Bio::SeqFeature::Generic;
		my ( $proxyfor, $otherf );
		if ( $type eq 'query' ) {
			( $proxyfor, $otherf ) = ( $hsp->query, $hsp->hit );
			$name ||= $hit->name;
		}
		else {
			( $otherf, $proxyfor ) = ( $hsp->query, $hsp->hit );
			$name ||= $result->query_name;
		}
		$proxyfor->score( $hit->bits ) unless ( $proxyfor->score );
		$min{$type}  = $proxyfor->start unless defined $min{$type}  && $min{$type}  < $proxyfor->start;
	   	$max{$type}  = $proxyfor->end   unless defined $max{$type}  && $max{$type}  > $proxyfor->end;
	   	$min{$other} = $otherf->start   unless defined $min{$other} && $min{$other} < $otherf->start;
	   	$max{$other} = $otherf->end     unless defined $max{$other} && $max{$other} > $otherf->end;
	   	#get strand
		my $strand;
		if ( $otherf->strand && $proxyfor->strand ) {
			$strand = $proxyfor->strand * $otherf->strand;
		}
		elsif ( $otherf->strand && !$proxyfor->strand ) {
			$strand = $otherf->strand;
		}
		elsif ( !$otherf->strand && $proxyfor->strand ) {
			$strand = $proxyfor->strand;
		}
		else {
			$strand = "+";
		}
		#set HSP feat values
	   	$feature->location($proxyfor->location,$otherf->location);
		$feature->add_tag_value( 'Parent', "$name:$countHit" );
	   	$feature->add_tag_value( 'Target', $name );
		$feature->add_tag_value( 'Target', $otherf->start );
		$feature->add_tag_value( 'Target', $otherf->end );
		$feature->strand($strand);
		$feature->source_tag( $proxyfor->source_tag );
		$feature->score( $proxyfor->score );
		$feature->frame( $proxyfor->frame );
		$feature->seq_id( $proxyfor->seq_id );
		$feature->primary_tag("match_part");

		push( @hspFeatTab, $feature );
		$seqid = $proxyfor->seq_id;
		$st = $proxyfor->source_tag;
	}
	return(\@hspFeatTab, \%min, \%max, $st, $seqid, $name);
}

sub makeGFF {

	my ($hit, $hspFeatTab, $st, $min, $max, $type, $name, $countHit, $seqid, $percentIdentity, $fracAligned) = @_;
	#make Hit feature
	my $matchf = Bio::SeqFeature::Generic->new(
												-start       => $min->{$type},
												-end         => $max->{$type},
												-strand      => $hit->strand($type),
												-primary_tag => 'match_set',
												-source_tag  => $st,
												-score       => $hit->bits,
												-seq_id      => $seqid
	);
	$matchf->add_tag_value( 'ID',                 "$name:$countHit" );
	$matchf->add_tag_value( 'Name',               $name );
	$matchf->add_tag_value( 'EValue',             $hit->significance );
	$matchf->add_tag_value( 'AlignLength',        $hit->length_aln($type) );
	$matchf->add_tag_value( 'GapNumber',          $hit->gaps($type) );
	$matchf->add_tag_value( 'PercentageIdentity', $percentIdentity );
	$matchf->add_tag_value( 'FractionAligned',    $fracAligned );
	#print match_set
	$out->write_feature($matchf);
	#print match_part
	for my $hspFeat (@$hspFeatTab) {
		$out->write_feature($hspFeat);
	}
}
exit;