#!/usr/bin/perl
# tdprotnear.perl

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

=item problems

 * alt-tr versus nearby duplicates: is that still a problem?
    - gnomon has some alt-tr, more for near-dmel
    - genewise has lots (40k mrna/13k genes for dvir?, 20k for others)
 * gnomon has utrs that expand gene size: use CDS overlap instead?
    - problem for Inside/Outside/Overlap classifying based on gene span
    
 
=cut

use constant SAMEBASE => 10; # for _sameloc, slop allowed in loca == locb
our $BINSIZE  = 1000 ; #was# 5000;
our $MINEVAL  = 1e-5;  our $MINBITSCORE = 0;
our $MINALIGN = 0.1; # 0.5; 
our $NEARDIST = 15000; # was 15k
our $NEARSTEPS= 3;
our $EQSAME=1;  our $EQINSIDE=3; our $EQOUTSIDE=4; our $EQOVERLAP=2; our $EQPARTOVER=5;  
our $EQNEAR=-2; our $EQFAR=-999; our $EQNONE=-1;  

my $groupok=1;
my $pctover=0;
my $pctinside=50;
my $mrnatype= $ENV{mrnatype} || "mRNA";
my $exontype= $ENV{exontype} || "CDS";
my $BY_GENE= 1;  
my ($clustersize, $clustertop, $debug, $dumpgenes, $nskip, $npoor, $useexons)
  = (0) x 20;
my ($ok, $tandemidpairs, $blastptable, $sourcegff, $othergff,
    @clusize, @neargroup, %ggroup, %genespan, %exons,  $skipids, $refgroup, $keeponly);

my $optok= GetOptions(
  "tandemidpairs=s", \$tandemidpairs, 
  "sourcegff=s", \$sourcegff, 
  "othergff=s", \$othergff, 
  "refgroup=s", \$refgroup, 
  "blastptable=s", \$blastptable, 
    "eval=s",\$MINEVAL, "bitscore=s",\$MINBITSCORE, "align=s",\$MINALIGN,
  "skipids=s", \$skipids, # == +keepids
  "mrnatype=s", \$mrnatype,"exontype=s", \$exontype,  
#  "keeponly=s", \$keeponly, # regex for id-group 
#  "NEARDIST=i", \$NEARDIST,  "NEARSTEPS=i", \$NEARSTEPS,  
  "neargroup=i", \@neargroup,   
  "pctoverlap=i", \$pctover,   
  "pctinside=i", \$pctinside,   
  "clustersize=i", \@clusize,   
  "dumpgenes=s", \$dumpgenes, 
  "debug=i", \$debug,  
  );

die "usage: perl tdprotnear.perl [ options ] 
     -tandemidpairs dsec-protnear15k.ids : as ID1.ID2 list of tandem genes 
  OR -blastptable dsec-pred6.blastp.gz : blastp table of matches
        -eval 1e-5 or -bitscore 100
     -source maingene.gff.gz : contains tandemid gene group 
     -other  allpredictors.gff.gz : other gene preds to compare to 
     opts: -skipids=idlist ; -pctoverlap=75 ; -pctinside=50 ; -dumpgenes=any|found
           -neargroup=9000 -neargroup=15000 : limit range of tandem distance
" unless($optok && ($tandemidpairs or $blastptable) && $sourcegff && $othergff);

# needs input fix for cluster size among all pairs, not just protnear15k.ids ..
#      -cluster=10 -cluster=80 : min,max dupl source genes
# ?? input blastp table and find pairs from that?


$pctover= $pctover/100.0 if($pctover);
$pctinside= ($pctinside) ? $pctinside/100.0 : 0.5; ## must have
if(@clusize) {
  $clustersize= $clusize[0];
  $clustertop= $clusize[1] || 0;
}

my @eqclass= ($EQSAME, $EQINSIDE, $EQOUTSIDE, $EQOVERLAP, $EQPARTOVER, $EQNONE); # no Near, Far yet
  ## (map{$EQNEAR * $_}(1..$NEARSTEPS))
  ##, $EQFAR); # same, near, far
my @eqnames= ("Same", "Inside", "Outside", "Overlap", "Partov", "None"); 
  ## (map{ my $kb= int((1 + $NEARDIST * $_)/1000); "Near".$kb."k"; } (1..$NEARSTEPS))
  ## ,"Far"); # same, near, far

my %eqnames;
for (my $i=0; $i<@eqclass; $i++) { $eqnames{ $eqclass[$i] }= $eqnames[$i]; }

# want tandem id pairs here instead ? see below
my %skipids=();
my $keepids = 0;
if($skipids) { 
  $keepids = ($skipids =~ s/^[\+\^]//) ? 1 : 0;
  open(F,$skipids); while(<F>){chomp; s/\s.*$//; $skipids{$_}++ if(/^\w/);} close(F); 
  }

#........ source duplicate pairs IN ............
my $nidpairs=0;
my %pairs=();

inputTandemidpairs($tandemidpairs, \%pairs, ) if($tandemidpairs); # return \$pairs ?

inputBlastpTable($blastptable, \%pairs ) if ($blastptable);

#........ end source duplicate pairs IN ............



#........ source gff IN ............
# add tandem cluster size filter here ?
# add neardist filter here?

my $nsourcin=0;
%genespan=(); 
$ok= ($sourcegff=~/\.gz$/) 
  ? open(IN,"gunzip -c $sourcegff |")
  : open(IN,$sourcegff);
die "open $sourcegff" unless($ok);

while(<IN>){
  next unless(/^\w/);
  chomp; my @v=split "\t";
  my($tref,$ftype,$tb,$te,$eval,$tor,$attr) = @v[0,2,3,4,5,6,8]; 
  my $gid;
  if($attr =~ m/ID=([^;]+)/) { $gid=$1; }
  elsif($attr =~ m/Parent=([^;]+)/) { $gid=$1; }
  
  my($ggroup)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); 
  next unless($ggroup =~ m/$refgroup/);
  $nskip++ and next unless( exists $pairs{$gid} ); # ignore other source genes

## clustersize is not good with current input id pairs: limited to 15k near tandems
## need all dupl. prots to count cluster sizes, then restrict by distance.
  if($clustersize) {
    my $ndups= scalar( keys %{ $pairs{$gid} } );
    my $skipit=0;
    $skipit=1 if($clustersize>0 and $ndups < $clustersize);
    $skipit=1 if($clustersize<0 and $ndups > -$clustersize);
    $skipit=1 if($clustertop>0 and $ndups > $clustertop);
    if($skipit) {
      delete $pairs{$gid};
      $nskip++ and next;
      }    
  }
  
  if($ftype =~ m/$mrnatype/) {
    $nsourcin++;  
    $genespan{$gid} = [$tb,$te,$tor,$tref,$gid]; # add CDS here?
    }
  elsif($useexons and $ftype =~ m/$exontype/) {
    push( @{$exons{$gid}}, [$tb,$te,$tor]);
    }
    
} close(IN);


# remove source genes outside of neardist range: need min,max distance
if(@neargroup>1) {
  my($mind,$maxd)= @neargroup; ($mind,$maxd)=($maxd,$mind) if($mind>$maxd); 
  foreach my $ga (sort keys %pairs) {
    foreach my $gb (sort keys %{$pairs{$ga}}) {
      my ($dist,$rev)= _mindistance( @{$genespan{$ga}}[0,1,2], @{$genespan{$gb}}[0,1,2] );
      $nskip++ and delete $pairs{$ga}{$gb} if($dist < $mind or $dist > $maxd);
      }
    my @tleft= sort keys %{$pairs{$ga}};
    delete $pairs{$ga} unless(@tleft);
    }
}

my %refgrouplocs=();
# if($refgroup) 
foreach my $gid (sort keys %genespan) {
  my($ggroup)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); 
  next unless($ggroup =~ m/$refgroup/);
  my ($tb,$te,$tor,$tref)= @{$genespan{$gid}}; 
  my @bins= (int($tb/$BINSIZE) .. int($te/$BINSIZE));
  foreach my $ib (@bins) { push @{$refgrouplocs{$ib}}, $gid; } 
  }

#........ end source gff IN ............



#........ other genes gff IN ............

%ggroup=(); # this is main result: other-genes x ref genes overlap types
my %predcount=(); #??
my %sourcecount=(); #??

$ok= ($othergff=~/\.gz$/) 
  ? open(IN,"gunzip -c $othergff |")
  : open(IN,$othergff);
die "open $othergff" unless($ok);

my $notherin=0;
while(<IN>) {
  next unless(/^\w/);
  chomp; my @v=split"\t"; 

  my($gid,$tref,$pctid,$alen,$tb,$te,$tor,$eval,$bits,$attr);

  ($tref,$tb,$te,$eval,$tor,$attr) = @v[0,3,4,5,6,8]; 
  if($attr =~ m/ID=([^;]+)/) {$gid=$1; }
  elsif($attr =~ m/Parent=([^;]+)/) { $gid=$1; }
  $tor="" if($tor eq ".");
  $notherin++;
  
  # can we filter out any other AltTr here by overlaps? esp. EISE_CGW
  # is input other GFF sorted by big-gene location?  sort -k4,4n -k5,5nr
  
  my($ggroup)= ($groupok) ? $gid =~ m/^(\D+)/ : ("all"); # BAD for non-grouping-prefixed IDs 
  next if($refgroup and $ggroup =~ /$refgroup/);
  
  ## add overlapfilter marks
  # my $overmarks= $attr; ##"" . join ",",@v[12..$#v];
  # $npoor++ and next if($overmarks and $overmarkdrop and $overmarks =~ m/$overmarkdrop/);
  # $npoor++ and next if($overmarkkeep and $overmarks !~ m/$overmarkkeep/);

  if(%skipids){
    if($keepids) {
      $nskip++ and next unless($skipids{$gid});
    } else {
      $nskip++ and next if($skipids{$gid}); # gene/transcript id here; or both? (exonid)
    }
  }

  my $notinrefgroup= 0;
  if($refgroup) {  
    my @sabins= (int($tb/$BINSIZE) .. int($te/$BINSIZE));
    my @samelocs=();
    foreach my $ib (@sabins) {  
      push @samelocs, @ { $refgrouplocs{$ib} } if($refgrouplocs{$ib});
      }

    my %didloc=();
    my $found=0; my $olap=0;  my $refid="";
    my $eq= $EQNONE; my $isrev= 0;
    foreach my $sgid (@samelocs) {  
      $genespan{$sgid} or next;
      next if ($didloc{$sgid.$gid}++);
      my ($gb,$ge,$gor,$gref)= @{$genespan{$sgid}};
      next unless($gref eq $tref);
      next unless(_isoverlap($gb,$ge,$tb,$te));
      
      $olap= 1; $refid= $sgid; 
      $eq= $EQOVERLAP; 
      if( _sameloc( $gref,$gb,$ge,$gor, $tref,$tb,$te,$tor) ) {
        $eq= $EQSAME;  
        
      } elsif (_isinside($gb,$ge,$gor, $tb,$te,$tor)) {
        $eq= $EQINSIDE; # means Other is inside Source
        
      } elsif (_isinside( $tb,$te,$tor, $gb,$ge,$gor)) {
        $eq= $EQOUTSIDE; # means Source is inside Other
        
      } elsif (! _isoverlapfull($gb,$ge,$gor, $tb,$te,$tor)) {
        $eq= $EQPARTOVER; 
      }

      $ggroup{$ggroup}{source}{$sgid}{$gid}= $eq;
      $ggroup{$ggroup}{other}{$gid}{$sgid}= $eq;

      ## these not useful
      ## $predfound{$ggroup}{$eq}{$gid}++ if $found;
      $predcount{$ggroup}{$eq}{$gid}++; #?? see tandynear2.perl
      $sourcecount{$ggroup}{$eq}{$sgid}++;  
      
      # group/refgene/othergene = overlap type,
      
      # check all matches from samelocs; 
      # NO: last; 
      
      ## dont care about near here
#       } else {
#         my $dist= _mindistance2($qb,$qe,$tb,$te);  # _isnear2($qb,$qe,$tb,$te);
#         $isrev= 0; if($qor && $tor) { $isrev = ($qor ne $tor) ? 1 : 0; }
#         $eq= $EQFAR; 
#         for (my $iq=1; $iq<=$NEARSTEPS; $iq++) {
#           if($dist <= $iq * $NEARDIST) { $eq= $iq * $EQNEAR; last; }
#           }
#       }

      }
      
    $genespan{$gid} = [$tb,$te,$tor,$tref,$gid] if($olap); # save other loc for dumpgenes
    $notinrefgroup=1 unless($olap);
  }

}
close(IN);

warn "# $refgroup: idpairs=$nidpairs; ngenes=$nsourcin; nother=$notherin; poor_ids=$npoor; skipped_ids=$nskip; \n" if($debug); ## && $npoor + $nskip>0);

#........ end other genes gff IN ............


print_bygene();  # output




#..........  main subs ......................

# this is much slower read than prepared dup ids table .. premake?
sub inputBlastpTable {
  my($tandemidpairs,$pairs,)=@_;
  $pairs={} unless(ref $pairs);
  
  die "inputBlastpTable: need -refgroup=source_predictor" unless($refgroup);
  my $ok= ($tandemidpairs=~/\.gz$/) 
    ? open(IN,"gunzip -c $tandemidpairs |")
    : open(IN,$tandemidpairs);
  die "open $tandemidpairs" unless($ok);
  my $isblast=1;
  my $isblatpsl=0;
  my %nself=(); #? save
  
  while(<IN>){  
    next unless(/^\w/);
    chomp; my @v= split "\t";
    my ($agroup,$bgroup,$aid,$bid,$bbits,$beval,$abmat,$as,$bs);
    
    if($isblatpsl) {    # .. blat format
      next unless(/^\d/);
      ($abmat,$aid,$as,$bid,$bs)= @v[0,9,10,13,14]; # blat psl format query,subject
      $bbits= $abmat; # for high score
      next if( $abmat < $as*$MINALIGN || $abmat < $bs*$MINALIGN || abs($as-$bs)>20 );
      
    } else {  # .. blast format  
      ($aid,$bid,$beval,$bbits)= @v[0,1,10,11]; # blast format query,subject
      $aid =~ s/\W+$//; $bid =~ s/\W+$//;
      if($MINBITSCORE) { next if($bbits < $MINBITSCORE); } 
      else { next if ($beval > $MINEVAL); }
    }

    ## only want refgroup same-predictor matches here?
    $nself{$aid}= $bbits and next if($aid eq $bid);
    next if($skipids{$aid} or $skipids{$bid});
    ($agroup)= $aid =~ m/^(\D+)/;
    ($bgroup)= $bid =~ m/^(\D+)/;
    next unless($agroup eq $refgroup and $agroup eq $bgroup);
    #? instead collect all HSP in pairs and separate out elsewhere?

    $pairs->{$aid}{$bid}= $bbits 
      unless (defined $pairs->{$aid}{$bid} and $pairs->{$aid}{$bid} > $bbits);
    ## $pairs->{$aid}{$bid}= $pairs->{$bid}{$aid}= $bbits; 
      # is this recip? NO; need to check both a/b scores
    $nidpairs++;
     
  } close(IN);
  return $pairs; # $nidpairs ?
}


sub inputTandemidpairs {
  my($tandemidpairs,$pairs,)=@_;
  $pairs={} unless(ref $pairs);
  
  my $ok= ($tandemidpairs=~/\.gz$/) 
    ? open(IN,"gunzip -c $tandemidpairs |")
    : open(IN,$tandemidpairs);
  die "open $tandemidpairs" unless($ok);
  
  # fixme for refgroup defined and idlist has many others
  while(<IN>){  
    if(m/^(\w+)[\.\t](\w+)/) { 
    my($ga,$gb)=($1,$2);
    #my($agroup)= $ga =~ m/^(\D+)/;
    #my($bgroup)= $gb =~ m/^(\D+)/;
    
    unless($refgroup) {
      ($refgroup)= $ga =~ m/^(\D+)/;
    }
    $pairs->{$ga}{$gb}= 1 if($ga =~ m/$refgroup/);
    $pairs->{$gb}{$ga}= 1 if($gb =~ m/$refgroup/); # is this recip?
    $nidpairs++;
    } 
  } close(IN);
  return $pairs; # $nidpairs ?
}


sub print_bygene {

  my @ggroup= sort keys %ggroup;
  my $manygroup= @ggroup>1;
  
  my $ttype = ($BY_GENE) ? "gene" : "exon";
  ## $ttype .=  ($COUNTKEYS) ? " ids" : " matches";
  my $cc = $dumpgenes ? "#t " : "";
  
  my @flags=();
  # push(@flags,reference_group => $refgroup) if($refgroup);
  # push(@flags,keep_only => $keeponly) if($keeponly);
  # push(@flags,overlap_filter => $overmarkfilter) if($overmarkfilter);    
  # push(@flags,by_gene => $BY_GENE);
  # push(@flags,count_type => (($BY_GENE)? "gene" : "exon"));
  # push(@flags,count_by => (($COUNTKEYS)? "ids" : "matches"));
  # push(@flags,min_align => $MINALIGN);
  # push(@flags,min_e_value => $MINEVAL);
  push(@flags,pct_overlap => $pctover) if($pctover);
  push(@flags,pct_inside => $pctinside) if($pctinside);

  push(@flags,tandem_pairs => $tandemidpairs);
  push(@flags,source_genes => $sourcegff);
  push(@flags,target_genes => $othergff);

  push(@flags,clusters => "$clustersize..$clustertop") if($clustersize);
  push(@flags,clusters => "all") unless($clustersize);
  push(@flags,distance_limit => "$neargroup[0]..$neargroup[1]") if(@neargroup);

  my %flags=@flags;
  my $flags= join ",", map{"$_=".$flags{$_}} sort keys %flags;
  
  print $cc,"Tandy-prot count of $ttype per predictor group\n";
  print $cc,"Options: $flags\n";
  print $cc,"Source : $refgroup\n";
  print $cc,join("\t","Group        ", "Stat     ",@eqnames),"\n";
  
  my ($pc,$pv,$tsame);
  foreach my $ggroup (@ggroup) {
    my $sgroup= sprintf("%-15s", substr($ggroup,0,14));
    my %ofound=(); # other found count classes?
    my %odone=(); my %tandone=();
    
#       $ggroup{$ggroup}{source}{$sgid}{$gid}= $eq;
#       $ggroup{$ggroup}{other}{$gid}{$sgid}= $eq;
    ## bring in source pairs here
    ## should use %pairs keys for @sgenes instead of %ggroup, to get missing cases

    # my @sgenes= sort keys %{ $ggroup{$ggroup}{source} }; 
    ## pairs can now contain cross-predictor matches ... want only both of refgroup?
    ##     next unless($agroup eq $refgroup and $agroup eq $bgroup);

    my @sgenes= sort keys %pairs; 
    my $nsourcegenes= 0; # scalar(@sgenes);
    my $ntandemgenes= 0;
    
    foreach my $gid (@sgenes) {
      next if($refgroup && $gid !~ m/$refgroup/);
      $nsourcegenes++;
      my @tandems= sort keys %{ $pairs{$gid} }; # recall we recip mapped these: dup count
      my @others = sort keys %{ $ggroup{$ggroup}{source}{$gid} };
      ## can we count sgenes,tandems w/ no other overlaps ? == EQNONE
      ## always have others keys here, but below tandems may be missing others
      
      foreach my $oid (@others) {
        next if($odone{same}{$gid.$oid}++);
        my $eq= $ggroup{$ggroup}{source}{$gid}{$oid} || $EQNONE;
        $ofound{same}{$eq}++; 
        }

      foreach my $tid (@tandems) {
        next if($refgroup && $tid !~ m/$refgroup/);
        $ntandemgenes++ unless($tandone{$tid}++);
        @others= sort keys %{ $ggroup{$ggroup}{source}{$tid} };
        push(@others, "Missing") unless @others;
        foreach my $oid (@others) {
          next if($odone{tandem}{$tid.$oid}++);
          my $eq= $ggroup{$ggroup}{source}{$tid}{$oid} || $EQNONE;
          $ofound{tandem}{$eq}++; 
          }
       }
      }

    my @ogenes= sort keys %{ $ggroup{$ggroup}{other} }; 
    my $nothergenes= scalar(@ogenes);
    foreach my $oid (@ogenes) {
      my @sogenes = sort keys %{ $ggroup{$ggroup}{other}{$oid} };
      ## here want to know how many tandems an other covers
      ## total count wont do ?
      
      foreach my $gid (@sogenes) {
        next if($odone{other}{$oid.$gid}++);
        my $eq= $ggroup{$ggroup}{other}{$oid}{$gid} || $EQNONE;
        $ofound{other}{$eq}++; 
        }

      }

    my @same  = map { $ofound{same}{$_}   || 0; } @eqclass;
    my @tandem= map { $ofound{tandem}{$_} || 0; } @eqclass;
    my @other = map { $ofound{other}{$_}  || 0; } @eqclass;
      
#     my @scounts= map { 
#       $pc= scalar( keys %{$sourcecount{$ggroup}{$_}} );
#       #$pc=0; map{ $pc += $_ } values %{$sourcecount{$ggroup}{$_}};
#       $pc;
#       } @eqclass;
#     my @ocounts= map { 
#       $pc= scalar( keys %{$predcount{$ggroup}{$_}} );
#       # $pc=0; map{ $pc += $_ } values %{$predcount{$ggroup}{$_}};
#       $pc;
#       } @eqclass;

    ## ** with %odone correction above for already counted pairs, 
    ##  source>other == tandem>other but for extra Missing/None case of tandem>other
    
    print $cc,join("\t",$sgroup,"source.genes",$nsourcegenes),"\n";  
    print $cc,join("\t",$sgroup,"tandem.genes",$ntandemgenes),"\n";  
    print $cc,join("\t",$sgroup,"other.genes", $nothergenes),"\n";  
    print $cc,join("\t",$sgroup,"source>other",@same),"\n";  
    print $cc,join("\t",$sgroup,"tandem>other",@tandem),"\n";  
    #fixme.now same as @same# print $cc,join("\t",$sgroup,"other>source",@other),"\n"; 
    
  #  # these counts are about same; ocounts =~ ofound{other}
  #   print $cc,join("\t",$sgroup,"source.count",@scounts),"\n";  
  #   print $cc,join("\t",$sgroup,"other.count",@ocounts),"\n";  
  
    print "\n"; 
    
    dumpgenes($ggroup,\@sgenes,) if($dumpgenes);
    
  }
  
}

sub dumpgenes {
  my($ggroup,$sgenes,)= @_;
  my %ofound;
  my $cc="";
  
  print "#i other genes overlapping source+tandem genes \n";
  print "#i", join("\t","Class",
    "Gene", "Location", "Other", "OverlapType", ),"\n";
  my $ig=0;
  
  foreach my $sgid (@$sgenes) {
    next if($refgroup && $sgid !~ m/$refgroup/);
    $ig++;
    my @tandems= sort keys %{ $pairs{$sgid} }; # recall we recip mapped these: dup count
    my @others = sort keys %{ $ggroup{$ggroup}{source}{$sgid} };
    my ($gb,$ge,$gor,$gref)= @{$genespan{$sgid}};
    push(@others,"Missing") unless(@others); # show no others...
    
    foreach my $oid (@others) {
      my $eq= $ggroup{$ggroup}{source}{$sgid}{$oid} || $EQNONE;
      # $ofound{same}{$eq}++; 
      # my ($ob,$oe,$oor,$oref)= @{$genespan{$oid}};
      print $cc, join("\t","source.$ig",
        $sgid, "$gref:$gb-$ge:$gor",
        $oid,  $eqnames{$eq},
        ),"\n";
      }

    my $it=0;
    foreach my $tid (@tandems) {
      next if($refgroup && $tid !~ m/$refgroup/);
      $it++;
      my @tothers= sort keys %{ $ggroup{$ggroup}{source}{$tid} };
      push(@tothers,"Missing") unless(@tothers); # show no others...
      my ($tb,$te,$tor,$tref)= 
        (exists $genespan{$tid}) ? @{$genespan{$tid}} : ("") x 10;
      foreach my $oid (@tothers) {
        my $eq= $ggroup{$ggroup}{source}{$tid}{$oid} || $EQNONE;
        # $ofound{tandem}{$eq}++; 

        print $cc, join("\t","tandem.$ig.$it",
            $tid,  "$tref:$tb-$te:$tor",
            $oid,  $eqnames{$eq},
            ),"\n";
            
        }
     }
     
    }
  print "\n";
}

#.......

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

sub _isoverlap {
  my($gb,$ge, $qb,$qe)= @_; 
  return ($gb <= $qe && $ge >= $qb) ? 1 : 0;
}

sub _isoverlapfull {
  my($tb,$te,$to, $qb,$qe,$qo)= @_;
  if($pctover) {
    my ($bb,$be)= ( _max($tb,$qb), _min($te,$qe) );
    my $maxo= abs($be - $bb);
    my $leno= _min( abs($qe - $qb), abs($te - $tb)) || 1;
    return ($maxo/$leno < $pctover) ? 0 : 1;
  }
  return 1; 
}

sub _sameloc {  
  my($ar,$ab,$ae,$aor, $br,$bb,$be,$bor)= @_;
  return ( $aor eq $bor &&  # strand/orient
          (abs($ab-$bb) < SAMEBASE) &&
          (abs($ae-$be) < SAMEBASE) ) ? 1 : 0;
}


sub _isinside {  #assume overlap? check for g-span >> q-span
  my($tb,$te,$to, $qb,$qe,$qo)= @_;
  if($qb > $tb and $qe < $te) {
    my $tlen= abs($te - $tb) || 1;
    my $qlen= abs($qe - $qb);  
    return 1 if( $qlen/$tlen < $pctinside); # 0.50 # 0.33;  #? use pctover here?
    }
  return 0;
}

sub _mindistance { 
  my($tb,$te,$to, $qb,$qe,$qo)= @_;
  my $rev= ($to and $qo and ($to eq "-" or $qo eq "-") and ($to ne $qo)) ? 1 : 0; # reversed not counted as overlap

  #? what of _isinside ?
  my $over= ($tb <= $qe && $te >= $qb) ? 1 : 0;
  if($over and $pctover) {
    my ($bb,$be)= ( _max($tb,$qb), _min($te,$qe) );
    my $maxo= abs($be - $bb);
    my $leno= _min( abs($qe - $qb), abs($te - $tb)) || 1;
    $over = 0 if $maxo/$leno < $pctover;
    return (1,$rev) unless($over);
  }

  if($over) { return (0,$rev); }
  else { # assume gb,ge, qb,qe are ordered
    my $bd= abs($tb - $qe); # g above q
    my $ed= abs($qb - $te); # g below q ; dont need gb - qb, ge - qe test
    return ($ed < $bd) ? ($ed,$rev) : ($bd,$rev);
    }
}

__END__

=item inputs

1. source gene id list
  $dp-protnear15k.id1 or $dp-protnear15k.ids (paired) : this one mixes predictors

2. source gene gff (add CDS gff?):
    $dp-$predictor-genes.note.gff.gz
    
3. target predictor gene gff (add CDS gff??)

=item outputs

  -- counts
  -- any marked gene gff?

=item notes

*** Measuring tandem gene model errors ***
Types: missed all, one model overlaps all, match source gene/miss tandem, 
      overlaps include: a. all dupl. exons matched, b. dupl. exons skipped
Protein near gene sets: for each predictor w/ duplicates, near-group, find overlap
	of genes, CDS-exons (?) among all/subset other predictors.  Count these:
	 1. substantial agreement (same and dupl. genes, mostly same exons)
	 2. (== 1) same-gene overlap/no-overlap (missed source gene?)
   3. extended-overlap (other has one gene/exon overlap for source and dupl genes)
   		3a. all dupl. exons matched; 3b. dupl. exons skipped
   4. tandem-missed (other found same source gene, but no model for tandem)
   5. ?
   
foreach source predictor gene id 
  - get source gene GFF (and CDS GFF?) == SG1
    - get all source tandem genes (foreach nearstep or for nearregion?) == @SDn
  - overlapfilter target predictor (gene and/or CDS) GFF == OG1, @ODn
  - count: SG1 x OG1 [1,2], SDn x OG1 [3 a,b], SDn x ODn [4], ? SG1 x ODn [5?]
  
Start with limited protnear geneid set as above:   $dp-protnear15k.id1
  
$td/tdprotnear.perl -debug=1 -neardist=1000 -nearstep=15 \
 -skipids=+$dp-protnear15k.id1  \
 xxxx? \
 > & ${dp}-tdprotnear1.txt &



=cut