#!/usr/bin/perl -w

=head1 NAME

blast_compare.pl - BLAST sequences with matching feature_ids, output a list of differences

=head1 SYNOPSIS

blast_compare.pl -ref refseq.fasta -new newseq.fasta [options]

 Required arguments:
    -ref      Reference sequnces (in fasta format or tab delimited format)
    -new      New sequnces (in fasta format or tab delimited format)

 Options:
    -type     AA (amino acids) or DNA (nucleotide bases), default is AA
    -synonyms Synonyms file that specifies sequence_ids that are synonymous
    -diff     Output format suitable for diff
    -output   Output file basename
    -help -h  Help

 Examples:

 Output a BLAST like report of difference for each sequence in reference.fasta and new.fasta
    blast_compare.pl -ref refernece.fasta -new new.fasta
	     
 Output a list of differences, one per line between the nucleotide sequences in reference.txt
 and new.txt (tab delimited files)
    blast_compare.pl -ref refernece.txt -new new.txt -diff -type DNA -output output_prefix

=head1 DESCRIPTION

Read two files,	either fasta files or tab delimited files in format:
    seq_id	feature_id	sequence

Uses BLAST to compare sequences with matching feature_id and reports all differences 
between the two sequences using a standard format (reference residue-position-new residue, e.g. A100T or C200R).
Insertions and deletions were reported using a dash, with the position determined on the reference sequence.

An option (-type) allows the user to specify whether DNA or AA sequences are being compared.  
By default, output is a BLAST alignment of the two sequences, which provids a visual display of differences.  
An option (-diff) generats a list of the differences between two strains, one per line (e.g. UL42: E331A), 
for each coding sequence in the genome.

=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'} . ':/usr/local/blast-2.2.24/bin';

use Bio::Tools::Run::StandAloneBlast;
use File::Basename;
use Getopt::Long;
use English;
use Carp;
use Pod::Usage;

my $needsHelp = '';
my $ref_filename;
my $new_filename;
my $diff_output = 0;
my $newseq_file;
my $seq_type = 'AA';
my $synonyms_filename;
my $output_base;

my $options_okay = &Getopt::Long::GetOptions(
	'ref|r=s'      => \$ref_filename,
	'new|n=s'      => \$new_filename,
	'diff|d'       => \$diff_output,
	'type|t=s'     => \$seq_type,
	'output|o=s'   => \$output_base,
	'synonyms|s=s' => \$synonyms_filename,
	'help|h'       => \$needsHelp
);

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

my @suffixlist = ( '.fasta', '.fa', '.txt' );
unless ($output_base) {
	$output_base =
	    basename( $ref_filename, @suffixlist ) . "_vs_"
	  . basename( $new_filename, @suffixlist );
}

# Read in files to hash
my $ref_seqs = &read_seqs($ref_filename);
my $new_seqs = &read_seqs($new_filename);

my ( $cmp_file, $homology_file );
unless ($diff_output) {
	open( $cmp_file, ">", $output_base . ".compare_table.txt" )
	  or croak(
		"Unable to write to file $output_base.compare_table.txt: $OS_ERROR");
	printf $cmp_file (
		"%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\n",
		"Gene(strand)",
		"Size Ref Seq",
		"Size New Seq",
		"% Ref Aligned",
		"% Ref Different",
		"% Identity",
		"Hit Length",
		"# Identical",
		"# Gaps",
		"start of match (new)",
		"end of match (new)",
		"start of match (ref)",
		"end of match (ref)",
		"# Changes",
		"# Change Events",
		"Change Summary",
	);
	open( $homology_file, ">", $output_base . ".homology.txt" )
	  or croak("Unable to write to file $output_base.homology.txt: $OS_ERROR");
}

my $synonyms;
if ($synonyms_filename) {
	$synonyms = read_synonyms($synonyms_filename);
}

# For each sequence in file1, find in file2
while ( my ( $k, $aa ) = each %$ref_seqs ) {

	my $s;
	if ( defined( $synonyms->{$k} ) ) {
		$s = $synonyms->{$k};
		if ( !$new_seqs->{$k} ) {
			$k = $s;
		}
	}
	if ( $new_seqs->{$k} ) {
		my $aa2 = $new_seqs->{$k};
		my ( $refseq, $newseq, $blast_program );
		if ( $seq_type eq 'AA' ) {
			$refseq = Bio::Seq->new(
				-display_id => "ref",
				-seq        => $aa,
				-alphabet   => 'protein'
			);
			$newseq = Bio::Seq->new(
				-display_id => "new",
				-seq        => $aa2,
				-alphabet   => 'protein'
			);
			$blast_program = 'blastp';
		} elsif ( $seq_type eq 'DNA' ) {
			$refseq = Bio::Seq->new(
				-display_id => "ref",
				-seq        => $aa,
				-alphabet   => 'dna'
			);
			$newseq = Bio::Seq->new(
				-display_id => "new",
				-seq        => $aa2,
				-alphabet   => 'dna'
			);
			$blast_program = 'blastn';
		} else {
			croak "Unknown sequence type '$seq_type', must be AA or DNA.\n";
		}
		my $blast_report = &blast( $refseq, $newseq, $blast_program );
		my $reported = 0;
		while ( my $result = $blast_report->next_result ) {
			while ( my $hit = $result->next_hit ) {
				while ( my $hsp = $hit->next_hsp ) {
					if ($reported) {
						next;
					}
					$reported++;
					my @diff =
					  &diff_summary( $hsp->query_string, $hsp->homology_string,
						$hsp->hit_string );
					if ($diff_output) {
						foreach my $d (@diff) {
							printf( "%s\t%s\n", $k, $d );
						}
					} else {
						my $diffs = join( ", ", @diff );
						my $event_cnt = &noncontig_count( \@diff );
						my ( $sbeg, $send ) = $hsp->range('hit');
						my ( $qbeg, $qend ) = $hsp->range('query');
						printf $cmp_file (
"%s\t%d\t%d\t%.2f\t%.2f\t%.2f\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%s\n",
							$k,
							$refseq->length(),
							$newseq->length(),
							( $hsp->length('total') / $refseq->length() ) * 100,
							(
								(
									$hsp->length('total') -
									  $hsp->num_identical()
								) / $refseq->length()
							  ) * 100,
							$hsp->percent_identity,
							$hsp->length('total'),
							$hsp->num_identical(),
							$hsp->gaps,
							$sbeg, $send, $qbeg, $qend,
							scalar(@diff),
							$event_cnt,
							$diffs
						);
						print $homology_file ("$k\n");
						print $homology_file ( $hsp->hit_string . "\n" );
						print $homology_file ( $hsp->homology_string . "\n" );
						print $homology_file ( $hsp->query_string . "\n\n" );
					}
				}
			}
		}
		if ( !$reported ) {
			printf $cmp_file (
"%s\t%d\t%d\t%.2f\t%.2f\t%.2f\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%d\t%s\n",
				$k,
				$refseq->length(),
				$newseq->length(),
				0,
				100,
				0,
				0,
				0,
				0,
				0,
				0,
				0,
				0,
				0,
				""
			);
		}
	} else {
		printf( "%s\tNo matching sequence found in %s!\n\n", $k,
			$new_filename );
	}
}

sub noncontig_count {
	my $diffr = shift;
	my @diff  = @$diffr;
	my $count = 0;
	my $pp    = -2;
	for ( my $i = 0 ; $i < scalar(@diff) ; $i++ ) {
		$diff[$i] =~ /\d+/;
		my $p = $&;
		if ( $p > $pp + 1 ) {
			$count++;
		}
		$pp = $p;
	}
	return $count;
}

sub diff_summary {
	my @diff     = ();
	my $ref      = shift;
	my $homology = shift;
	my $new      = shift;

	my $ref_position;
	for ( my $i = 0 ; $i < length($homology) ; $i++ ) {
		my $h = substr( $homology, $i, 1 );
		my $r = substr( $ref,      $i, 1 );
		my $n = substr( $new,      $i, 1 );
		if ( $r ne "-" ) {
			$ref_position++;
		}
		if ( $h ne " " && $h ne "+" ) {
			next;
		} else {
			push( @diff, $r . ($ref_position) . $n );
		}
	}
	return @diff;
}

sub blast {
	my $seq1          = shift;
	my $seq2          = shift;
	my $blast_program = shift;
	my @params        = (
		'program' => $blast_program,
		'-F'      => 'F',
		'--task'  => 'dc-megablast',

		#'-G'      => '6',
		#'-E'      => '2',
		'-r' => '2',
		'-q' => '-3',
	);
	my $factory = Bio::Tools::Run::StandAloneBlast->new(@params);
	my $blast_report = $factory->bl2seq( $seq1, $seq2 );
	return $blast_report;
}

sub read_seqs {
	my $filename = shift;
	open( my $file, "<", $filename )
	  or croak "Can't open '$filename': $OS_ERROR";
	my %seqs;
	if ( $filename =~ /\.fa(sta)?$/ ) {
		my $inseq =
		  Bio::SeqIO->new( -file => "<$filename", -format => 'fasta' );
		while ( my $seq = $inseq->next_seq ) {
			$seqs{ $seq->id } = $seq->seq;
		}
	} else {
		while ( my $line = <$file> ) {
			chomp($line);
			my ( $seq, $feature, $aa ) = split( /\t/, $line );
			$seqs{$feature} = $aa;
		}
	}
	return \%seqs;
}

sub read_synonyms {
	my $filename = shift;
	open( my $file, "<", $filename )
	  or croak "Can't open '$filename': $OS_ERROR";
	my $synonyms;
	while ( my $line = <$file> ) {
		chomp($line);
		my ( $orig, $syn ) = split( /\t/, $line, 2 );
		$synonyms->{$orig} = $syn;
	}
	return $synonyms;
}

# 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 $ref_filename ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message =>
			  "Cannot read reference sequence file: '$ref_filename!'\n"
		);
	}
	if ( !-e $new_filename ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message => "Cannot read new sequence file: '$new_filename!'\n"
		);
	}
	if ( $synonyms_filename && ( !-r $synonyms_filename ) ) {
		pod2usage(
			-exitval => 2,
			-verbose => 1,
			-message => "Cannot read synonyms file: $synonyms_filename!\n"
		);
	}
}
