#!/usr/bin/perl -w

=head1 NAME

migrate_annotations.pl - Migrate annotations (e.g. gene locations) from one version of genome to another 

=head1 SYNOPSIS

migrate_annotations.pl -r reference.fasta -g reference.gff -n new.fasta [-v input_gff_version] [-help] [-debug]

 Required arguments:
    -r        Reference fasta file
    -g        Reference gff file
    -n        New fasta file

 Options:
    -v        GFF file version (default is 3)
    -help -h  Help

 Examples:
 Migrate genes locations based on sequence.fasta to ones based on new_sequence.fasta
	migrate_annotations.pl -r sequence.fasta -g genes.gff -n new_sequence.fasta > new_genes.gff
     
 Migrate genes locations using a GFF version 2 file
	migrate_annotations.pl -r sequence.fasta -g genes.gff -v 2 -n new_sequence.fasta > new_genes.gff

=head1 DESCRIPTION

Take two sequences and a GFF file that corresponds to one sequence. BLAST align the two sequences.
If there is a close match, convert positions in the GFF file to positions in the second sequence.

NOTE: Fasta files must have the same sequences in the same order!!!

=head1 AUTHOR

Lance Parsons <lparsons@princeton.edu>

=head1 LICENSE

This script is licensed by the Simplified BSD License
See LICENSE.TXT and <http://www.opensource.org/licenses/bsd-license>

Copyright (c) 2011, Lance Parsons
All rights reserved.

=cut

use strict;

# Blast binaries should be in PATH
$ENV{'PATH'} = $ENV{'PATH'} . ':/Users/lparsons/blast-2.2.21/bin';

use Bio::SeqIO;
use Bio::Tools::GFF;
use Bio::Tools::Run::StandAloneBlast;
use Getopt::Long;
use English;
use Carp;
use Pod::Usage;

my $needsHelp = '';
my $refseq_file;
my $refgff_file;
my $newseq_file;
my $gff_version = 3;
my $debug       = 0;

my $options_okay = &Getopt::Long::GetOptions(
	'refseq|r=s'     => \$refseq_file,
	'refgff|g=s'     => \$refgff_file,
	'newseq|n=s'     => \$newseq_file,
	'gffversion|v=i' => \$gff_version,
	'debug'          => \$debug,
	'help|h'         => \$needsHelp
);

# Check to make sure options are specified correctly and files exist
&check_opts();

my @params = (
	program => 'blastn',
	'-F'    => 'F',
);

my $refseq = Bio::SeqIO->new(
	-file   => "<$refseq_file",
	-format => 'fasta'
);

my $newseq = Bio::SeqIO->new(
	-file   => "<$newseq_file",
	-format => 'fasta'
);

# specify input via -fh or -file
my $gffio =
  Bio::Tools::GFF->new( -file => "$refgff_file", -gff_version => $gff_version );
my @features;

# Scan through GFF file and store features in array
while ( my $feature = $gffio->next_feature() ) {
	push( @features, $feature );
}
$gffio->close();

if ($debug) {
	print STDERR (
"Feature name\tSubject_Id\tPercent_id\tAlignment_Length\tIdentical\tGap_Openings\tQuery_Start\tQuery_End\tOriginal_Start\tOriginal_End\tNew_Start\tNew_End\tStrand\n"
	);
}

print "##gff-version 3\n";
while ( my $seq = $refseq->next_seq ) {
	my $seq2 = $newseq->next_seq();

	foreach my $feature (@features) {
		if ( $feature->seq_id() eq $seq->id() ) {

			my $annotationCollection = $feature->annotation();
			my @gene_annotations     =
			  $annotationCollection->get_Annotations('gene');
			my $feature_name = $gene_annotations[0];
			my $refsubseq    = $seq->trunc( $feature->start, $feature->end );

			# Run blast
			my $factory = Bio::Tools::Run::StandAloneBlast->new(@params);

			# SearchIO object
			my $hit_reported = 0;
			my $hit_count    = 0;

			my $blast_report = $factory->bl2seq( $refsubseq, $seq2 );
			while ( my $result = $blast_report->next_result ) {
				while ( my $hit = $result->next_hit ) {
					while ( my $hsp = $hit->next_hsp ) {
						if (
							(
								$hsp->length('total') >=
								$refsubseq->length() * .95
							)
							&& $hsp->strand('hit') == 1
						  )
						{
							$hit_count++;
							my ( $sbeg, $send ) = $hsp->range('hit');
							my ( $qbeg, $qend ) = $hsp->range('query');

							if ($debug) {
								printf STDERR (
"%s\t%s\t%f\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\n",
									$feature_name,
									$hit->name,
									$hsp->percent_identity,
									$hsp->length('total'),
									$hsp->num_identical(),
									$hsp->gaps,
									$qbeg,
									$qend,
									$feature->start,
									$feature->end,
									$sbeg,
									$send,
									$hsp->strand('hit')
								);
							}
							my $newFeature = $feature;
							$newFeature->seq_id( $seq2->id() );
							$newFeature->start($sbeg);
							$newFeature->end($send);
							my $gffio =
							  Bio::Tools::GFF->new( -gff_version => 3 );
							printf( "%s\n", $newFeature->gff_string($gffio) );
						}
					}
				}
			}
			if ( $hit_count != 1 ) {
				printf STDERR (
					"%d hits found for:\n %s!\n\n",
					$hit_count, $feature->gff_string()
				);
			}
		}
	}
}

# Check for problem with the options or if user requests help
sub check_opts {
	if ($needsHelp) {
		pod2usage( -verbose => 2 );
	}

	if ( !$options_okay ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message => "Error specifying options."
		);
	}
	if ( !-e $refseq_file ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message => "Cannot read fasta file: '$refseq_file'!\n"
		);
	}
	if ( !-e $refgff_file ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message => "Cannot read gff file: '$refgff_file'!\n"
		);
	}
	if ( !-e $newseq_file ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message => "Cannot read fasta file: '$newseq_file'!\n"
		);
	}
}
