#!/usr/bin/env perl
use strict;
use SeqOp;
use File::Basename qw( fileparse );
use Getopt::Long;
# Start with query protein file (one of 500 multifasta chunks)
# BLAT against genome
# pslCDnaFilter -localNearBest=0.2 -maxAligns=3 001.psl 001.filtered.psl
# psl2gff-fast.pl 001.filtered.psl
# sortgff 001.filtered.gff | bedtools merge -d 4000 > 001.filtered.merged.bed 

# Then for each line of 001.filtered.merged.bed, run genewise and format gff properly and uniquely (qname, aln num, $pid);
# Dependencies: SeqOp.pm (from Tyler), genewise, blat, pslCDnaFilter (kent source tree) 

my $doblat = 0;
my $file = "";			#file is the query file
my $genomedir = "";
my $genome = "";
my $TMP = $ENV{TMPDIR};
GetOptions(
	   'q:s'       => \$file,
	   'g|genome:s'      => \$genome
	  );
my ($base,$path,$ext) = fileparse($file,qw(\.fa \.fasta \.mfa));
my $pid = $$;
my $out = "$base.genewise.$pid.gff";
#print STDERR "Writing to $out\n";

if (!-e "$base.raw.psl" || $doblat) {
  print STDERR "blat -maxIntron=30000 -q=prot -t=dnax $genome $file $base.raw.psl\n";
  `blat -maxIntron=30000 -q=prot -t=dnax $genome $file $base.raw.psl`;
}
if (!-e "$base.filtered.psl" || $doblat) {
  print STDERR "pslCDnaFilter -localNearBest=0.2 -maxAligns=3  $base.raw.psl $base.filtered.psl\n";
  `pslCDnaFilter -localNearBest=0.2 -maxAligns=3  $base.raw.psl $base.filtered.psl`;
}

my %region;
print STDERR "psl2gff-fast.pl $base.filtered.psl\n";
`psl2gff-fast.pl $base.filtered.psl`;
print STDERR "sort -k9,9 -k1,1 -k7,7 -k4,5n $base.filtered.gff > $base.filtered.sorted.gff\n";
`sort -k9,9 -k1,1 -k7,7 -k4,5n $base.filtered.gff > $base.filtered.sorted.gff`;
open (IN, "<$base.filtered.sorted.gff") or die "couldn't open $base.filtered.psl\n";
while (<IN>) {
  next if m/^#/;
  my @f = split;
  my $str = substr($f[6],1,1);
  my $seq = $f[0];
  my $p = $f[8];
  my $start = $f[3];
  my $end = $f[4];
  if (!exists($region{"$p-$seq-$str"})) {
    push @{$region{"$p-$seq-$str"}},[$start,$end];
  } else {
    foreach my $aln ( @{$region{"$p-$seq-$str"}}) {
      if ($start < $aln->[1] + 4000) {
	$aln->[1] = $end;
      } else {
	push @{$region{$p-$seq-$str}},[$start,$end];
      }
    }
  }

}
close IN;
open (OUT, ">$out") or die "couldn't open $out\n";
my $count = 1;
my $sbjct = '';
print STDERR "Running genewise on each hit...\n";
my $Q = 0;
my %seen;
my $U = 100000000000;
my $V = 0;
my $SBJCT = 0;;

my $E = 1;
foreach my $prot (keys %region) {
  my ($q,$sbjct,$str)= split("-",$prot);
  foreach my $aln (@{$region{$prot}}) {
    my $strand = '-tfor';
    if ($str eq '-') {
      $strand = '-trev';
    }
    my $alnstart =  $aln->[0]; #start
    my $alnend = $aln->[1];      #end
    my $seq = "$TMP/$sbjct.$pid.fa";
    my $SEQ = "$TMP/$SBJCT.$pid.fa";
    my $qseq = "$TMP/$q.$pid.fa";
    my $QSEQ = "$TMP/$q.$pid.fa";
    print STDERR "$q vs. $sbjct\n";
    if ($q ne $Q) {
      unlink "$QSEQ" if -e "$QSEQ";
      $count = 1;
    }
    $Q=$q;
    if ($sbjct ne $SBJCT) {
      unlink "$SEQ" if -e "$SEQ";
    }
    $SBJCT = $sbjct;
    if (!-e "$TMP/$q.$pid.fa") {
      if (!-e "$file.index") {
	`fastaindex -f $file -i $file.index`;
      }
      print STDERR "fastafetch -f $file -i $file.index -q $q > $qseq\n";
      `fastafetch -f $file -i $file.index -q $q > $qseq`;
    }
  
    if (!-e "$seq") {
      if (!-e "$genome.index") {
	`fastaindex -f $genome -i $genome.index`;
      }
      `fastafetch -f $genome -i $genome.index -q $sbjct > $seq`;
    }

    my $l = SeqOp::fasta_length($seq);
    #print STDERR "length: $l\n";

    $U = $alnstart-2000; if ($U < 1){$U = 1; }
    $V = $alnend+2000;  if ($V > $l){$V = $l;}
    print STDERR "genewise $qseq $seq $strand -u $U -v $V -pseudo -gff | grep cds | \n";
    open GENEWISE, "genewise $qseq $seq $strand -u $U -v $V -pseudo -gff | grep cds | ";
    while (<GENEWISE>) {
      s/\S+genewise-prediction-\S+$/ID=$q.$count.$pid/;
      my @gff = split "\t",$_;
      if ($gff[3]>$gff[4]) {
	my $tmp = $gff[3];
	$gff[3]=$gff[4];
	$gff[4]=$tmp;
      }
      print OUT join("\t",@gff);
    }
    $count++;
  }

}
#unlink "$base.filtered.psl";
#unlink "$base.raw.psl";

close OUT;
print STDERR "done!\n";
