#!/usr/bin/perl
# overbestgene2.perl

=item overbestgene2
  
  special case of gene gff overlap detection:
   - input one gene/cds gff file with quality scores
      and many overlapping predictions (e.g. exonerate predicts)
   - output best by score non-overlapping subset genes + cds
     (some overlap possible, want all mostly distinct predictions)
     
 gzcat exonerate-gldmoj.gff3.gz | $td/overbestgene2.perl -in stdin > exonerate-gldmoj-best.gff &
 #collect_gff=1660943, ngene=763282   + 6735 geneids with terepeat > about same...
 #done kept=23000, skipped=744001
    vs kept=15038, skipped=748244 for simple filter (includes terepeat filter) 
    
    
=cut

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

use constant SAMEBASE => 0; # for _sameloc, slop allowed in loca == locb
# use constant { ACT_DROP=>1, ACT_KEEP=>2, ACT_MARK=>3, ACT_MARK_WITH_ID=>4 };
# use constant { kOVERLAP=>1, kSAMELOC=>2, kNEARLOC=>3, kBESTGENE=>4, kBESTEXONS=>5, };
# 
# use constant BESTEXON_SCORE => 2;

our $BINSIZE  = 50000 ; #was# 1000; want large steps for this case?
our $NEARDIST =  500; # was 15k; needs to be < BINSIZE

our $debug=1;
my ($showskips,$overlaps,$overlaplist,%geneexons,$scoretype,$input,$itype,$action,$actid,$ok,$mark);
# my ($overtype,$typeover,$sametypes,$pctover,)= (kBESTGENE,"",1,0);
my $exontypes="CDS,exon,match_part,HSP";
my $mrnatypes="mRNA,match"; # fixme ?? need this
my $badattr="";
# my $correct_partial_hsp=1;
$scoretype= "score";
my $pctover=0;

my $optok= GetOptions(
  "input=s", \$input,  
#  "overlaps=s", \$overlaps, 
#  "typeover=s", \$typeover, 
  "exontypes=s", \$exontypes, 
  "badattr=s", \$badattr, 
  "mrnatypes=s", \$mrnatypes, 
  "pctover=i", \$pctover, 
#  "NEARDIST=i", \$NEARDIST, 
  "BINSIZE=i", \$BINSIZE, 
#  "mark=s", \$mark, 
#  "itype=s", \$itype,  
  "scoretype=s", \$scoretype, # return not ID= but other attribute or score
#  "action=s", \$action, 
  "skipshow!", \$showskips, 
  "debug!", \$debug, 
  );

die "
usage:  perl overbestgene2 
  -exontypes $exontypes  -mrnatypes $mrnatypes
  -input genes.gff|stdin > bestgenes.gff
" unless($optok);
# -act keep|mark|markid -mark=bestgene

$pctover= $pctover/100.0 if($pctover);
$exontypes =~ s/[,; ]+/|/g;
$mrnatypes =~ s/[,; ]+/|/g;

my $inh= *STDIN;
$ok = ($input =~ /.gz$/) ? open($inh,"gunzip -c $input |") 
      : ($input =~ /^(stdin|-)/) ? $inh= *STDIN
      : open($inh,$input);
die "bad -input=$input" unless($ok);


my ($genebins, $genes, $exons,) = collect_gff($inh); 

my %keepgene=(); # print
my %skipgene=(); # optionally list ..
my($nskip, $nkeep)= (0,0);

foreach my $ref (sort keys %$genebins) {
  my @bins= sort{$a <=> $b} keys %{$genebins->{$ref}};
  foreach my $ib (@bins) { # each 1000 bases?? too many steps?
    my @didexons=();
    my @sgenes = sort _sort_refscoreloc @ { $genebins->{$ref}{$ib} };
    foreach my $gref (@sgenes) {
      my($ref,$src,$typ,$tb,$te,$tscore,$to,$tph,$tattr,$gid)= @$gref;
      my $rexons= $exons->{$gid};
      # next unless($rexons); # are we missing exon-gene id links?
      
      if($keepgene{$gid}) {
        push(@didexons, @$rexons);
      } else {
        my $nover=0; my $xbad=0;
        foreach my $ex (@$rexons) {
          $nover++ if (overlaps($ex->[3], $ex->[4], \@didexons));
          $xbad++  if ($badattr && $ex->[8] =~ m/$badattr/);
          }
        if($nover > 0 || $xbad > 0) { 
          $skipgene{$gid}++; $nskip++;  
          }  
        else { 
          $keepgene{$gid}++; $nkeep++;
          push(@didexons, @$rexons);
          putgene( $gref, $rexons); 
          }
        }
      
      }
      
    }
}

if($showskips) {
  print "# skipped genes ",(".") x 50," \n";
  foreach my $gid (sort keys %skipgene) {  
    putgene( $genes->{$gid}, undef, "skip=1");  # $exons->{$gid},
    }
}

warn"#done kept=$nkeep, skipped=$nskip\n" if $debug;

#.......................................

sub _sort_refscoreloc # $overlaps == mrnatypes
{
# ($ref,$src,$typ,$tb,$te,$tscore,$to,$tph,$tattr,$gid)
# fixed: ?? with ref 1st this isn't working; big score not at front
  return ($a->[0] cmp $b->[0]) # ref
      || ($b->[5] <=> $a->[5]) # score
      || ($a->[3] <=> $b->[3]) # begin
      || ($b->[4] <=> $a->[4]); # end

#   return # ok for score
#          ($b->[5] <=> $a->[5]) # score
#       || ($a->[3] <=> $b->[3]) # begin
#       || ($b->[4] <=> $a->[4]); # end
}


sub putgene {
  my($gref,$exons,$flag)= @_;
  foreach my $ex ($gref, @$exons) {
    my($ref,$src,$typ,$tb,$te,$tscore,$to,$tph,$tattr,$gid)= @$ex;
    $tattr .=";$flag" if($flag);
    print join("\t", $ref, $src, $typ, $tb, $te, $tscore, $to, $tph, $tattr),"\n";
  }
}


sub _min { return ($_[1] < $_[0]) ? $_[1] : $_[0]; }
sub _max { return ($_[1] > $_[0]) ? $_[1] : $_[0]; }

sub overlaps {
  my($tb, $te, $locs) = @_;
  foreach my $rloc (@$locs) {
    ### my ($lb,$le,$lid)= @{$rloc}[0,1,3];
    my ($lb,$le,$lid)= @{$rloc}[3,4,9];
    my $over= ($tb <= $le && $te >= $lb) ? 1 : 0;
    # add option to screen out trival 5% UTR overlaps .. pctoverlap=i
    if($over and $pctover) {
      my ($bb,$be)= ( _max($tb,$lb), _min($te,$le) );
      my $maxo= abs($be - $bb);
      my $leno= _min( abs($le - $lb), abs($te - $tb)) || 1;
      $over = 0 if $maxo/$leno < $pctover;
      }
    return 1 if($over);
    }
  return 0;
}

sub collect_gff
{
  my($gff)= @_;
  my (%genebins, %genes, %exons);
  my ($ng,$nr)= (0,0);
  while(<$gff>){
    next unless(/^\w/); chomp;
    my $line= $_;
    my($ref,$src,$typ,$tb,$te,$tscore,$to,$tph,$tattr)= split"\t";
    $nr++;
    
    $tscore=0 if($tscore eq ".");
    my $score= $tscore; # want this to select; w/ field choice?
    if($scoretype =~ /^score/i) { $score=$tscore; }
    elsif($tattr =~ m/\b$scoretype=([^;]+)/) { $score=$1; }
    
    my($gid); 
    if($tattr =~ m/\bID=([^;]+)/) {  $gid=$1; }
    elsif($tattr =~ m/\bParent=([^;]+)/) {  $gid=$1; }
    unless(defined $gid) { $gid = "N".$ng; }
    
    my $rloc= [$ref,$src,$typ,$tb,$te,$score,$to,$tph,$tattr,$gid]; # change to string; save mem
    ##? only need access these:  $gid, $tb, $te, $tscore
    # $rloc= [ $ref, $tb, $te, $tscore, $gid, $line]; 
    
    if($typ =~ /^($mrnatypes)$/) { # $overtype == kBESTEXONS and 
      $genes{$gid}= $rloc; $ng++; # should be uniq/id ?
      my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE));
      foreach my $ib (@bins) {
        push( @{$genebins{$ref}{$ib}}, $rloc); #? or store $gid
        }
    } elsif($typ =~ /^($exontypes)$/) { #$overtype == kBESTEXONS and 
      push( @{$exons{$gid}}, $rloc);
      } 
    }
  warn"#collect_gff=$nr, ngene=$ng\n" if $debug;
  return (\%genebins, \%genes, \%exons);
}