#!/usr/bin/env perl

use warnings;
use strict;

#DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD

=pod

=head1 NAME

seqdiff.pl - Compare two versions of UniProtKB FASTA sequence files and adjust
             the corresponding set of PolyPhen-2 precomputed files to accomodate
             the differences found

=head1 SYNOPSIS

seqdiff.pl [options] seqfile1 seqfile2

where options are:

  -p PATH     specify path to the precomputed directory tree; default is to
              use PRECOMPATH configuration variable retrieved from the current
              PolyPhen-2 installation configuration

  -r          apply any changes found to precomputed files; default is to show
              the differences but do not touch the files

  -v          verbose mode, prints pairs of entries with mismatched sequences

  -h          display extended help

=head1 DESCRIPTION

Input files should be FASTA sequence collections in the PolyPhen-2 .seq variant of
FASTA format, i.e., with single-line unwrapped sequences (files in this format are
created by uniprot.pl script). Any differences between the two sets will be reported:
entries in the first file missing from the second one, entries added in the second
file and entries present in both files but with non-identical sequences. If -r option
was specified, all files corresponding to the entries missing from the second FASTA
file or with non-identical sequences will be deleted from the current set of PolyPhen-2
precomputed files found under the directory tree pointed to by the PRECOMPATH
configuration variable (see config/paths.cnf file).

=head1 AUTHORS

Ivan Adzhubey

=head1 SUBVERSION

 $LastChangedDate: 2012-12-21 14:45:27 -0500 (Fri, 21 Dec 2012) $
 $LastChangedRevision: 416 $
 $LastChangedBy: ivan $

=cut

#DDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDDD

use Getopt::Std;
use Pod::Usage;

BEGIN {
  $Getopt::Std::STANDARD_HELP_VERSION = 1;
  use FindBin qw($RealBin);
  BEGIN {
    unless ($ENV{'PPH'}) {
      my $binpath = $RealBin;
      $binpath =~ s|/[^/]+$||;
      $ENV{'PPH'} = $binpath;
    }
  }
  use lib "$ENV{PPH}/perl";
}

use PPH::Config;
use PPH::Align;
use PPH::Seq;

my %opts;
pod2usage() unless @ARGV;
getopts('hrvp:', \%opts) or pod2usage();
pod2usage(-verbose=>2) if $opts{'h'};
pod2usage('ERROR: Missing input file(s)') unless @ARGV > 1;
my $seqfile1 = shift;
my $seqfile2 = shift;

-f $seqfile1 or die "File does not exist: $seqfile1\n";
-f $seqfile2 or die "File does not exist: $seqfile2\n";

my $path;
if      ($opts{'p'}) {
  $path = $opts{'p'};
} elsif (defined $CONFIG{'PRECOMPATH'} && length $CONFIG{'PRECOMPATH'}) {
  $path = $CONFIG{'PRECOMPATH'};
} else {
  die "Fatal: PRECOMPATH directory undefined\n";
}

-d $path or die "PRECOMPATH directory not found: $path\n";

if ($opts{'r'}) {
  print "## PRECOMPATH=$path\n";
} else {
  print "## PRECOMPATH=$path (DRY-RUN)\n";
}

my @subdirs_missing;
foreach my $subdir (qw{alignments blastfiles lock profiles structures}) {
  -d "$path/$subdir" || push @subdirs_missing, $subdir;
}
die('PRECOMPATH subdirectories missing: ' . join(' ', @subdirs_missing) . "\n") if @subdirs_missing;

# Load a list of accessions from first file
my %accessions;
open(my $SEQ1, '<', $seqfile1) or die "Can't open file: $seqfile1\n";
my $out = do { local $/; <$SEQ1> };
close($SEQ1);
defined $out && length $out or die "Empty or unreadable file: $seqfile1\n";
my $i = 0;
foreach my $entry (split /^>/m, $out) {
  next unless length $entry;
  my ($head, $seq) = split /\n/m, $entry, 2;
  my ($def, $dbid, $acc, $name, $desc) = PPH::Align::_parse_defline($head);
  $accessions{$acc} = $i;
  $i++;
}

my %seen;
my %missing;
my %altered;
$i = 0;
foreach my $acc (sort { $accessions{$a} <=> $accessions{$b} } keys %accessions) {
  my ($acc1, $name1, $desc1, $seq1) = PPH::Seq::fetch_FASTA($acc, $seqfile1);
  warn("Failed to fetch entry $acc from file: $seqfile1\n"), next unless defined $acc1;
  my ($acc2, $name2, $desc2, $seq2) = PPH::Seq::fetch_FASTA($acc, $seqfile2);
  if (defined $acc2) {
    if ($seq1 eq $seq2) {
      ;
    } else {
      $altered{$acc} = 1;
      print "Altered: $acc\n";
      print "(1) $desc1\n$seq1\n(2) $desc2\n$seq2\n" if $opts{'v'};
      delete_files($acc) if $opts{'r'};
    }
    $seen{$acc} = 1;
  } else {
    $missing{$acc} = 1;
    print "Missing: $acc\n";
    delete_files($acc) if $opts{'r'};
  }
  $i++;
}

# Load and scan second file for the list of accessions added
%accessions = ();
open(my $SEQ2, '<', $seqfile2) or die "Can't open file: $seqfile2\n";
$out = do { local $/; <$SEQ2> };
close($SEQ2);
defined $out && length $out or die "Empty or unreadable file: $seqfile2\n";
$i = 0;
foreach my $entry (split /^>/m, $out) {
  next unless length $entry;
  my ($head, $seq) = split /\n/m, $entry, 2;
  my ($def, $dbid, $acc, $name, $desc) = PPH::Align::_parse_defline($head);
  $accessions{$acc} = $i;
  $i++;
}

my %added;
if (%accessions) {
  foreach my $acc (sort { $accessions{$a} <=> $accessions{$b} } keys %accessions) {
    unless ($seen{$acc}) {
      $added{$acc} = 1;
      print "Added:  $acc\n";
    }
  }
}

print "## Totals:\n";
print "##    Missing: ", (%missing ? scalar(keys %missing) : 0), "\n";
print "##    Altered: ", (%altered ? scalar(keys %altered) : 0), "\n";
print "##    Added:   ", (%added   ? scalar(keys %added)   : 0), "\n";

sub delete_files {
  my $acc = shift;
  my $path = $CONFIG{'PRECOMPATH'};
  my @rm = qw(/bin/rm -fv);
  foreach my $subdir (qw{alignments blastfiles lock profiles structures}) {
    my $basename = "$path/$subdir/$acc";
    if      ($subdir eq 'alignments') {
      system @rm, $basename.'.aln';
      system @rm, $basename.'.msa.fa';
    } elsif ($subdir eq 'blastfiles') {
      system @rm, $basename.'.blast';
    } elsif ($subdir eq 'lock') {
      system @rm, $basename.'.lock';
    } elsif ($subdir eq 'profiles') {
      system @rm, $basename.'.prf';
    } elsif ($subdir eq 'structures') {
      system @rm, $basename.'.pdb.blast';
    }
  }
}
