#!/usr/bin/perl
# blastnear.perl

# tandem gene analysis by proteins (07aug)
# combining idchains (blast in) and protnear (idchange, gene gff in); 
# see protnear.perl
# use:  grep -v '^#' dgri_pred6.blastp | sort -k13,13nr | cat dgri_tandy6jmd_all.genes - | \
#  $td/blastnear.perl -blast -near 5000 -bits=150 -norecip -shownear > dgri_tandy6jmd_all.blastnear8 &


use strict;
use constant SHOWREV => 1;

my $NEARDIST= 15000; # common with tandy exon choice
my $SAMEDIST= 10;  #  do _isoverlap; but use to screen odd cases 
my $PCT_ISINSIDE = 0.40; # 0.40; ??
   # for _isinside, gene inside other < pct length of other; problems w/ alttr confusion
   # versus tandem gene models inside over-spread model
our $BINSIZE  = 1000 ; #was# 5000;
my $MIXGROUPS= 0; # segregate or mix predictor groups;
my $ALLPAIRSCORE= 0; # for id pairs output of all,not just nearpair: farpair,scafpair not used

my (%generef,  %genebe, %galtids, %genepairs, %modelerrs,
    %isscaf,%isfar,%isnear,%issame,
    %distclass, # move above %is___ to distclass
    %nearsteps, %revsteps, %nearpair, %nearrev, %farpair, %scafpair);
my $shownear=0;
my $mrnatype= $ENV{mrnatype} || "mRNA";
my $ngene= $ENV{ngene} || 0;
my $dupgeneok= 1;
my $clustersize= 0; my $clustertop=0;
my @clusize;
my $source="";
my $pairfilter="";
my $keepgffattr="";
my $isblat=0; my $isblast=1; 
my $removesubsets=1; my $usebest=0; 
my $eval= 1e-10; my $bitscore=0; 
my $align=0.1; # align=0.5 is too restrictive: Dmoj/Est6 tandems fail
my $skipids="";
my $reciprocal=1;
my $skiplocations=0;
my $groupok=1;
my ($pctalign,$pctbits,$pctover,$debug)=(0,0,0,0);

use Getopt::Long;

my $optok= GetOptions(
  "blast!",\$isblast, "blat!",\$isblat, 
  "eval=s",\$eval, "bitscore=s",\$bitscore, "align=s",\$align,
  "pctbitscore=s",\$pctbits, "pctalign=s",\$pctalign,
  "mrnatype=s", \$mrnatype, 
  "pairfilter=s", \$pairfilter, 
  "keepgffattr=s", \$keepgffattr, 
  "ngene=s", \$ngene, 
  "source=s", \$source, 
  "groupok!", \$groupok,  
  "dupgeneok!", \$dupgeneok,
  "reciprocal!",\$reciprocal,
  "MIXGROUPS!", \$MIXGROUPS, 
#  "skiplocations!",\$skiplocations, # doesnt look useful but to double check total pairs w/o loc filter
  "pctoverlap=i", \$pctover, 
  "NEARDIST=i", \$NEARDIST,  #"NEARSTEPS=i", \$NEARSTEPS,  
  "clustersize=i", \@clusize,  # \$clustersize,  
  "shownear!", \$shownear,  
  "showall!", \$ALLPAIRSCORE,  
  "debug!", \$debug,  
  );

die "usage: cat genes.gff genes.blastp | perl blastnear.perl [ options ] \
  opts: -mrnatype mRNA  -ngene 20000 -neardist $NEARDIST
        -blast|-blat : input format: blast -m 8,9 or blat psl
        -pctoverlap=20  : separate 'same' genes/alttr from near genes
        -pairfilter=xyz , -keepgffattr=xyz : keep match pairs (2nd) or gff attrib ; pattern
        -shownear : id list for near pairs
        -eval=$eval | -bitscore=$bitscore | -align=$align : blast/blat score filter
        -pctbits=50 , -pctalign=50  : remove poor pairs by bitscore/align relative to self-bits/align
        -noreciprocal : dont require recip. blast match
        -skiplocations : dont use genes.gff
        -cluster 10 -cluster 20 : range of cluster sizes
        -mixgroups ; -showall (for shownear > showall)
" unless($optok);

$isblast=0 if($isblat);
$shownear=1 if ($ALLPAIRSCORE); #== showall
$pctover= $pctover/100.0 if($pctover);
$pctbits= $pctbits/100.0 if($pctbits);
$pctalign= $pctalign/100.0 if($pctalign);
if(@clusize) {
  $clustersize= $clusize[0];
  $clustertop= $clusize[1] || 0;
}

my $noreciprocal= ! $reciprocal;
my %skipids=();
# if($skipids && open(F,$skipids)) { while(<F>){chomp; $skipids{$_}=1;} close(F); }

my ($noloc,$nskip, $nchain)=(0) x 9;
my %noloc=(); # by group
my (%nself, %selfsize, %pairs, %palign, %pstarts, %gename);

while(<>){
if(/^\W/) { next; }
chomp; 
my @v= split "\t";

if(/\t$mrnatype\t/ and @v == 9){ # gene/mRNA gff for loc,ids
  my($gr,$gs,$gt,$tb,$te,$gp,$to,$tx, $attr)= @v; # split /\t/; 
  my($id)= m/ID=([^;\s]+)/;  # $attr
  #? collect Parent to help alt-tr filtering? some have (gnomon), some not (genewise)
  next unless( $keepgffattr eq "" or $attr =~ m/$keepgffattr/);
  my($ggroup)= ($groupok) ? $id =~ m/^(\D+)/ : ("all");  
  $generef{$id}= $gr;  # $gm{$id}= int(($tb+$te)/2); 
  $genebe{$id}= [$tb,$te,$to]; # add $gr? need to change sub _xxx() params if so
  my ($gname)= m/Parent=([^;\s]+)/;
  $gename{$id}= $gname || ""; # for JGI diff predictors?
  $source ||= $gs;
  } 

elsif( @v > 11) { # blast nc=13, blat.psl nc=19? no. columns

  if(/_GLEAN_/){ s/[a-z]{4}_GLEAN_/GLEAN_/g; } # fix for _GLEAN_ id change .gff vs .aa

  my ($bbits, $aid, $bid, $abmat, $as, $bs, $beval);
  unless($isblast) {    # .. blat format; use isblat 
    next unless(/^\d/);
    
#     ($abmat,$aid,$as,$bid,$bs)= @v[0,9,10,13,14]; # blat psl format query,subject
#     next if($skipids{$aid} or $skipids{$bid});
#     $bbits= $abmat; # for high score
#     if($aid eq $bid) {
#       $selfsize{$aid}= $abmat;
#       $nself{$aid}= $bbits;
#       next;
#       }
#     next if( $abmat < $as*$align || $abmat < $bs*$align || abs($as-$bs)>20 );

## ... rewrite ...
    my( $mismatches, $rep_matches, $orient, $qstart, $qend, $tstart, $tend);
    my ($blocksizes, $qstarts, $tstarts);
    ( $abmat, $mismatches, $rep_matches, $orient,
      $aid, $as, $qstart, $qend,
      $bid, $bs, $tstart, $tend,
      $blocksizes, $qstarts, $tstarts,
      )= @v[0..2, 8..16, 18..20];

    next if($skipids{$aid} or $skipids{$bid});
    $bbits= $abmat; # for high score
    if($aid eq $bid) {
      $selfsize{$aid}= $abmat;
      $nself{$aid}= $bbits;
      next;
      }
    next if( $abmat < $as*$align || $abmat < $bs*$align || abs($as-$bs)>20 );#?

    my @blocksizes = split( /,/ , $blocksizes );
    my @qstarts    = split( /,/ , $qstarts );
    my @tstarts    = split( /,/ , $tstarts );
    my $npart      = @qstarts;
    for(my $p=0; $p<$npart; $p++) {
      my($qstart,$tstart,$len)= ($qstarts[$p], $tstarts[$p], $blocksizes[$p]);
      my $qend= $qstart+$len;
      my $tend= $tstart+$len;
      $qstart++; $tstart++; # move to 1-origin
      my $ps= join( "\t", $qstart, $qend, $tstart, $tend);
      $pstarts{$aid}{$bid} .= $ps."\n";
    }

  } else {  # .. blast format  
    next unless(/^\w/);
    ($aid,$bid,$beval,$bbits,$abmat)= @v[0,1,10,11,3]; # blast format query,subject
    
    # using prot22_modXX.blastp out got dingbats on aid,bid : clean
    $aid =~ s/\W+$//; $bid =~ s/\W+$//;
    next if($skipids{$aid} or $skipids{$bid});
    if($aid eq $bid) {
      $selfsize{$aid}= $abmat;
      $nself{$aid}= $bbits;
      next;
      }
    # ^ save self size (if prot) here from abmat/align-size
    if($bitscore) { next if($bbits < $bitscore); } else { next if ($beval > $eval); }

    ## add also, for blastp at least, q,s match start,end, to check doubled errors
    ## want to keep all per pair, e.g.  q1 x s1   401 .. 800 x 1 .. 400 ; 1 .. 400 x 1 .. 400
    ## need only qe, se ends now; drop qb,sb for space?? no
    
=item NOTE to Self: Hash of String cuts oodles of memory out of Hash of Hash
    
    Switch to Hash of String INSTEAD OF Hash of Hash:
    this cut mem use by gigabytes for genome blast data set
    
=cut

    #?? save mem on this big hash w/ strings? YES ***
    my $ps= join "\t", @v[6,7,8,9];
    $pstarts{$aid}{$bid} .= $ps."\n";
    
    ## need to watch for b > e reverse
    # my @pstarts= @v[6,7,8,9]; # my($qb,$qe, $sb, $se)
    # $pstarts{$aid}{$bid} or $pstarts{$aid}{$bid}=[];
    # push( @{$pstarts{$aid}{$bid}}, \@pstarts); # do i want to save all this big data set? keep only best scored?

  }

#?  $pairs{$aid}{$bid}++;  # want reciprocal align 
  $pairs{$aid}{$bid}= $bbits unless ($pairs{$aid}{$bid} and $pairs{$aid}{$bid} > $bbits);
  $palign{$aid}{$bid}= $abmat unless ($palign{$aid}{$bid} and $palign{$aid}{$bid} > $abmat);

  }

}


# defined %galtids or 
# ** problem with Alt Tr test here ; mixing predictor groups **
%galtids= genespans() unless($skiplocations); # ok; from %genebe

## FIXME into pairs2near stats
# my $chains= pairs2chains();
# chains2near($chains); # sets hashes: nearsteps, isnear, issame, ...

pairs2near(); # sets hashes: nearsteps, isnear, issame, ...

OUTPUT: 
  outdist();
  # out1();

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

sub filterAlttr {
  my($geneids, $di)= @_;
  
  my($gref,$gb,$ge,$go)=("",-1,0,"+");
  my (@altids,@sepids, %altids, %sepids);
  my $ng= @$geneids;
  
  for (my $i=0; $i<$ng; $i++) {
    my $idi= $geneids->[$i];
    next if($altids{$idi});
    $sepids{$idi}++;
    my($igroup)= ($groupok) ? $idi =~ m/^(\D+)/ : ("all");  
    $gref= $generef{$idi};
    ($gb,$ge,$go)= @{$genebe{$idi}}; 

    # for (my $j=0; $j<$ng; $j++) 
    for (my $j=$i+1; $j<$ng; $j++) 
    {
      next if($j==$i);
      my $idj= $geneids->[$j];
      next if($altids{$idj});
    
      my($jgroup)= ($groupok) ? $idj =~ m/^(\D+)/ : ("all");  
      my $tref= $generef{$idj};
      my($tb,$te,$to)= @{$genebe{$idj}};
      # fixme: ? need to separate predictors: only where group(idi) = group(idj) ??
      
      if ( 
        ( $jgroup eq $igroup) and  # $MIXGROUPS || 
        $tref eq $gref and
        _isoverlap($gb,$ge,$go, $tb,$te,$to) and
        not _isinside($gb,$ge,$go, $tb,$te,$to, 0) #? or keep as alttr ?
        ) { 
        $altids{$idj}++; delete $sepids{$idj};
        }
      else {
        $sepids{$idj}++;
        }
      } # j
    } # i
    
   @sepids= sort{ $pairs{$di}{$b} <=> $pairs{$di}{$a} or $a cmp $b } keys %sepids;
   @altids= sort{ $pairs{$di}{$b} <=> $pairs{$di}{$a} or $a cmp $b } keys %altids;
   return (\@sepids, \@altids);  
    
}



sub pairs2near {

  my($nc,$nscaf,$nfar,$nnear,$nsame)=(0)x9; 

  # my $ng= scalar(@aid); 
  # for my $i (0..$#d)  my $di= $d[$i];   

  # my @aid= sort keys %pairs; #? sort by bitscore instead?
  my @aid= sort{ $nself{$b} <=> $nself{$a} or $a cmp $b } keys %pairs;  
  my %moderrtemp=();
  my %genepairloc=(); 
  
  for my $di (@aid) {
    my($igroup)= ($groupok) ? $di =~ m/^(\D+)/ : ("all");  
    
    #unless($skiplocations) {
    $noloc{$igroup}++ and next unless($genebe{$di}); # count skips
    #}
    
    #? sort other than ID-lexical? by genebe size? pairs score!
    my @bid= sort{ $pairs{$di}{$b} <=> $pairs{$di}{$a} or $a cmp $b } keys %{$pairs{$di}};

# add filter: pct bits ab  / min(nself a, nself b)
# add filter: pct align ab / min(selfsize a, selfsize b) AND check for double-dup size b >> size a

    # filter @bid by below recipro and genebe to count cluster size here
    #($_ ne $di) and  # already removed
    @bid= grep {
           ( exists $genebe{$_} ) #$skiplocations || 
        && ( $pairfilter eq "" || m/$pairfilter/ ) 
        && ( $noreciprocal || $pairs{$_}{$di} )
        && ( !$pctbits  || $pctbits  <  $pairs{$di}{$_} / _min($nself{$_},$nself{$di}) )  
        && ( !$pctalign || $pctalign < $palign{$di}{$_} / _min($selfsize{$_}, $selfsize{$di}) )  
        } @bid; 

#?? need to mark all filtered ids as genepairs{$di}{xxx} ??    
    my($bsep,$balt) = (\@bid, []);
    unless($skiplocations) {
      ($bsep,$balt) = filterAlttr (\@bid, $di);
      @bid= @$bsep;
    }
    
    for my $dj (@$balt) {
      $genepairs{$di}{$dj}++; # store distance instead of count for pairs?    
      $genepairs{$dj}{$di}++;  
      }
      
    my $ncluster= @bid;
    if($clustersize) {
      $nskip++ and next if($clustersize>0 and $ncluster < $clustersize);
      $nskip++ and next if($clustersize<0 and $ncluster > -$clustersize);
      if($clustertop>0) {    
        $nskip++ and next if( $ncluster > $clustertop);
      }
    }
    
    for my $dj (@bid) {
      # next if($dj eq $di); # already removed 
      # next unless( $genebe{$dj} );
      # next unless( $noreciprocal or $pairs{$dj}{$di} );
      # next unless($removesubsets or !$didab{$aid}{$bid});
      
      ## problem with MIXGROUPS here; dont want mixed matches counted in group near,far counts
      ## just want pair ids (dumpgenes) and model error scores
      
      my($jgroup)= ($groupok) ? $dj =~ m/^(\D+)/ : ("all");  
      next unless( $MIXGROUPS || $jgroup eq $igroup);
      my $grp= $igroup; # "$igroup.$jgroup" #??
      
            
      next if($genepairs{$di}{$dj} ); # not $dupgeneok and ()
      # next if($genepairs{$di}{$dj} or $genepairs{$dj}{$di}); # not $dupgeneok and ()
      $genepairs{$di}{$dj}++; # store distance instead of count for pairs?
      # $genepairs{$dj}{$di}++;  
      
      $genepairloc{$di}= [ $generef{$di}, @{$genebe{$di}} ]
        if($genebe{$di}); # saving only di loc here?? drop {$dj} subscript
      
        # one class of tandem model errors: 2+ gene exon sets joined as one
        #  test that palign{di}{dj} > 0.9 * smallsize and bigsize > 1.8 * smallsize
      my $errdoubled = 0;
      my $errskipover= 0;
      
      $errdoubled= ( $selfsize{$di} > 1.8 * $selfsize{$dj}
         && $palign{$di}{$dj} > 0.9 * $selfsize{$dj}) ? 1 : 0;
       # ^ also add test that both 1/2s of 2x prot are matched by same 1x?
       
        # also want other class, dup exons inside gene model
        # can we use genebe to get at gene span diffs?
        # this tests if gene-i spans >> gene-j, but with same size protein, i.e. intron structure changed
        # ** must have evidence of internal exons skipped over to call this err;
        #    other predictors? or better tandy exons
use constant NO_ERRSKIPOVER => 0;
      unless(NO_ERRSKIPOVER or $errdoubled or $skiplocations) {  
        my $ijratio = _sizediff( @{$genebe{$dj}}, @{$genebe{$di}} );  # test returns size-2 / size-1
        # $errskipover= 1 if ($ijratio > 1.9 && $palign{$di}{$dj} > 0.9 * $selfsize{$di});
        $errskipover= 1 if ($ijratio > 1.75 && $palign{$di}{$dj} > 0.75 * $selfsize{$di});
        
        warn "skipov=1 $di.$dj\n" if($debug && $errskipover); # this shows many ; not seen output
        }
      
      # $modelerrs{$igroup}{'doubled'}{$di}++ if($errdoubled);
      $moderrtemp{$igroup}{'doubled'}{$di}{$dj}++ if($errdoubled);
        ## ^ need 2+ other gene matches to assign this err class accurately
        ## need to reduce false doubled calls;
        ## i.e. a 2x prot size matching 1x is often real biology; need more what?
        
      # $modelerrs{$igroup}{'skipover'}{$di}++ if($errskipover);
      $moderrtemp{$igroup}{'skipover'}{$di}{$dj}++ if($errskipover);
      # should mark these with id output to vis. check; see below
      # this misses errors where there isn't a duplicate w/i predictor; compare other predictors?
      # these are counts/pairs not /gene; want latter? e.g. 1 doubled gene can count 4.. times
       ## i.e. 1277 count == 60 genes; 
       
      ## FIXME: need to handle multiple predictors, di >> dj a,b,c
      ## distinguish cases where a-pred has one near long model, b-pred has 2,3 near short models
      ## want to find max number of nearby, distinct genes from all predictors
      ## look at all dj/bid for distance from di/aid and overlap among selves
      
      ## ... if skiplocations all this is useless : count all as Unlinked/scaf ...
      
      if($skiplocations or $generef{$di} ne $generef{$dj}) { 
        $distclass{$grp}{'scaf'}{$di}++;  
        if($ALLPAIRSCORE) { $nearpair{"$di.$dj"}= 399; } # = UNLINKSCORE
        # $distclass{$grp}{'scaf'}{$dj}++;
        } 
      else {
        my ($dist,$rev)= _mindistance( @{$genebe{$di}}, @{$genebe{$dj}} );
        # returns dist:1 == _isinside ; dist:2  == small overlap 

        ### $dist=0 if($dist < 10); # for now; mostly subset gene calls of same larger model
        $dist=0 if($dist < 2); # for now; mostly subset gene calls of same larger model

        # pass dist==2, NO; bad data? as distinct genes w/ small overlap (set w/ pctoverlap)
        # got this when genes MOSTLY overlap ; i.e. same model
        
             # look for cases in @jid of 2+ inside genes matching larger model = tandem mismodel
        if ($dist<=0) { 
          $distclass{$grp}{'same'}{$di}++; 
          # $distclass{$grp}{'same'}{$dj}++;
          if($ALLPAIRSCORE) { $nearpair{"$di.$dj"}= 199; } # = OVERLAPSCORE; skip this one??
          }
        else {
          if($dist>$NEARDIST) { 
             $distclass{$grp}{'far'}{$di}++; 
             # $distclass{$grp}{'far'}{$dj}++;
             $nearsteps{$grp}{999}{$di}++; 
             # $nearsteps{$grp}{999}{$dj}++;
             if($ALLPAIRSCORE) { $nearpair{"$di.$dj"}= 299; } # = FARSCORE
             if($rev) { 
              $revsteps{$grp}{999}{$di}++; 
              # $revsteps{$grp}{999}{$dj}++; 
              }
             } 
          else { 
            $distclass{$grp}{'near'}{$di}++; 
            # $distclass{$grp}{'near'}{$dj}++;
            #? handle special dist == 1 as $distep=0 ?
            my $distep= 1 + int( $dist / 1000);
            
            $distep= 0 if($dist < 30);
            # ^ these are not interesting generally: GeneWise has 10x others for eg Dmoj
            #  they are partial/subset gene calls matching 1/2 of longer gene from other pred
            #  consider as altTr or bad gene call; replace _isinside w/ issame ?
            #  what of some cases where inside calls are better tandem gene models?
            #
            # * what would be good to know is where there are 
            #   tandem 2+ genes from other predictors
            #   inside of di gene model, similar to the larger one (e.g. overextended/bad model)
            
              ## keep this group-less ?            
              ## add score for errdoubled, errskipover ?
            my $pairscore= $distep;
            # .. defer this mark to validate below
            #x $pairscore .=",errD" if($errdoubled); # wait till have 2+ gene matches to mark this
            #x $pairscore .=",errO" if($errskipover);
            $nearpair{"$di.$dj"}= $pairscore; 
                
            #NOT recip here# $nearpair{"$dj.$di"}++; 

            $nearsteps{$grp}{$distep}{$di}++; 
            # $nearsteps{$grp}{$distep}{$dj}++;
            if ($rev) { 
              $revsteps{$grp}{$distep}{$di}++; 
              # $revsteps{$grp}{$distep}{$dj}++; 
              }
            } 
          }
        }
                
      } # dj..
    } # di..


  validmodelerrs(\%genepairloc,\%moderrtemp);
  
  
#   my $nc= scalar(keys %nc);  
#   my $nself= scalar(keys %nself);
#   my $nsing= $nself - $ng;
#   print "# nself=$nself; npaired=$ng; nsingle=$nsing; duplicates: ngroup=$na; nuniqgene=$nc; ngroupedgene=$nb; rmsubsets=$removesubsets \n"; 

}


sub validmodelerrs 
{
  my($genepairloc, $moderrtemp)= @_;
  my (%locs,);
  
  foreach my $di (sort keys %$genepairloc) {
    my ($pgrp)= $di =~ m/^(\D+)/;
    my ($tref,$tb,$te,$to)= @ { $genepairloc->{$di} }; # same for all dj
    my @bins= ($tb/$BINSIZE .. $te/$BINSIZE);
    foreach my $ib (@bins) { push @{$locs{$tref}{$ib}}, $di; }  
    }
  
  my @ggroup= sort keys %$moderrtemp; # distclass;
  foreach my $grp (@ggroup) { 
    foreach my $err ( sort keys %{ $moderrtemp->{$grp} } ) {
    
      foreach my $di ( sort keys %{ $moderrtemp->{$grp}{$err} } ) {
        my @dj = sort keys %{$moderrtemp->{$grp}{$err}{$di}};
        
        # there are mistakes in these calls (less than skipover); 
        # i.e. a 2x prot size matching 1x is often real biology; need more what?
        # adding part count and minlen as below seems to get only real doubles, but
        #  drops count way down to 10s not 100s of errs .. probably right
        #  this also suffers from missing double cases with no shorter dupl. by same predictor
        #  *** use other predictors here also, double of 2x short gene test should be valid across preds
        
        if($err eq 'doubled' and @dj > 0) { # drop to @dj > 0 and use pstarts > 1

          ## pstarts:  all per pair, q1 x s1: 401 .. 800 x 1 .. 400 ; 1 .. 400 x 1 .. 400
          my $isdoubled= 0;
          my %part; my $lastgrp;
          my %grpdoubled=();
          foreach my $dj (@dj) {
          
            ## mixed groups tricky here; shouldnt count part b/n groups
            my ($pgrp)= $dj =~ m/^(\D+)/;
            if($MIXGROUPS) {
              next if($pgrp ne $grp && $pgrp =~ /CGW|CEX|CGM|GMP/); # skip genemappers w/ problem alt-tr
              %part=() if($pgrp ne $lastgrp);
              $lastgrp= $pgrp;
            }
          
            #$pstarts{$di}{$dj} or next;
            #my @pstarts= @{$pstarts{$di}{$dj}}; # want @p > 1 ?
            my $pstarts= $pstarts{$di}{$dj} or next;
            my @pstarts= split "\n",$pstarts;
            
            my ($leni,$lenj)= ($selfsize{$di},$selfsize{$dj});
            next if($leni<$lenj);
            my $minlen=  0.8 * $lenj;
            # next if( $palign{$di}{$dj} < $minlen);
            foreach my $ps (@pstarts) {
               ## my($qe,$se)= @$ps;
               #my($qb,$qe, $sb, $se)= @$ps;
               my($qb,$qe, $sb, $se)= split "\t",$ps;
               
               ($qb,$qe)=($qe,$qb) if($qb>$qe); # rev
               ($sb,$se)=($se,$sb) if($sb>$se); # rev
               next if( 1+$se-$sb < $minlen); ## do i need this instead of palign{i,j} ?
               my $part= int($qe/$se); # expect $qe increase >> $se for doubling; want 2+ parts
               $part{$part}++; # all qe < se == 0; filter any?
               # ^ test for all parts among @dj, doesnt need to be same dj to be doubled?
              }
              
            my $n= scalar(keys %part); # now over all @dj matches
            $isdoubled = $n if ($n>1 && $isdoubled<$n); # some tripled ..
            $grpdoubled{$pgrp}= $n if ($n>1 && $grpdoubled{$pgrp} < $n);
            # no need to test also the align strength again?
            # my $strength= ( $selfsize{$di} > 1.8 * $selfsize{$dj}
            #    && $palign{$di}{$dj} > 0.9 * $selfsize{$dj}) ? 1 : 0;
            
          }
          
          if($isdoubled) {
            $modelerrs{$grp}{'doubled'}++; # drop di here ; below
            my $gname= $gename{$di}||"";
            if($MIXGROUPS) {
              map{ 
                my ($pgrp)= $_ =~ m/^(\D+)/; my $dv= $grpdoubled{$pgrp}||0;  
                $nearpair{"$di.$_"} .= ",errD$dv\t$gname" if($dv);
                } @dj; 
            } else {
              map{ $nearpair{"$di.$_"} .= ",errD$isdoubled\t$gname";} @dj; 
            }
             # ^ for MIXGROUPS, should mark this only per group that qualifies
            }
        }

        #? could/should extend this test to *all* gene locs for $grp using genebe ?
        ## ** this isn't getting real dupl skipovers; none of flagged genes have
        ## inside dupl. of other groups ... why?
        ## part of problem is GWise,Exonr,+ partial/short alt-tr called as dup inside genes
        
        if($err eq 'skipover' and @dj > 0) {
         # need to test if other group predictors have dupl genes internal to this one
         
         ## my $dj= $dj[0];  # same di loc for all dj
         my ($tref,$tb,$te,$to)= @ { $genepairloc->{$di} };
          # drop this big genepairloc hash for genebe ? getting to be a mem pig here
          # see sub genespans() that bins all genebe
          
         # count other group 2+ dups  inside this di loc
         my @bins= ($tb/$BINSIZE .. $te/$BINSIZE);
         my (%dgenes, %dgroup);
         foreach my $ib (@bins) { 
            map{ $dgenes{$_}++;} @{$locs{$tref}{$ib}};
          } 
         foreach my $dg (sort keys %dgenes) {
          my ($pgrp)= $dg=~m/^(\D+)/;
          next if($groupok and $pgrp eq $grp); # skip self ?? NO at least optionally
          next if($pgrp =~ /CGW|CEX|CGM|GMP/); # skip genemappers w/ problem alt-tr
          $dgroup{$pgrp}++; # no. dup genes in di gene span
          }
         my $isdup=0;
         foreach my $pgrp (sort keys %dgroup) { $isdup++ if ($dgroup{$pgrp}>1); }
         
         $isdup=2 if(!$groupok and %dgenes);
         
         if ($isdup>1) {  ## bad here for 1 predictor group; need isdup>0
          $modelerrs{$grp}{'skipover'}++;
          my $gname= $gename{$di}||"";
          if($MIXGROUPS) {
            map{ 
              my ($pgrp)= $_ =~ m/^(\D+)/; my $dv= $dgroup{$pgrp}||0;  
              $nearpair{"$di.$_"} .= ",errO$dv\t$gname" if($dv>1);
              } @dj; 
          } else {
            map{$nearpair{"$di.$_"} .= ",errO$isdup\t$gname";} @dj; # dj may be 0 ?
          }
         }
        }
       
       ## also big 3rd error class/group not seen here; can we add w/ other predictor info?
       #  : genes not called with dupl that others do call; includes skipover, double and misses
       
       
      }
    }
  }


}


sub outopts {
  my @ggroup= sort keys %distclass;
  my $manygroup= @ggroup>1;
  
  my $cc = "# "; ## $dumpgenes ? "#t " : "";
  
  my @flags=();
  push(@flags,pct_align => $pctalign);
  push(@flags,pct_bitscore => $pctbits);
  push(@flags,min_e_value => $eval);
  push(@flags,min_bitscore => $bitscore) if $bitscore;
  push(@flags,pct_overlap => $pctover) if($pctover);
  ## push(@flags,pct_inside => $PCT_ISINSIDE) if($PCT_ISINSIDE);
  push(@flags,noreciprocal => $noreciprocal) if($noreciprocal);

  push(@flags,pairfilter => $pairfilter) if($pairfilter);
  push(@flags,keepgffattr => $keepgffattr) if($keepgffattr);
  push(@flags,mixed_groups => $MIXGROUPS) if($MIXGROUPS);

  my %flags=@flags;
  my $flags= join ",", map{"$_=".$flags{$_}} sort keys %flags;
  
  print $cc,"Tandy-blastnear counts of predictor groups\n"; # $ttype per
  print $cc,"Options: $flags\n";
  # print $cc,"Source : $refgroup\n";
  # print $cc,join("\t","Group        ", "Stat     ",@eqnames),"\n";
}


sub outdist {

  my($gscaf,$gfar,$gnear,$gsame, $gpairs)=(0)x9;
  my($pscaf,$pfar,$pnear)= (0)x9;

  outopts();
  
## add predictor groupings loop .... need all %isxxx{group} ?

# replace this with range of near..far ? see tandynear2.perl
#   $gscaf= scalar(keys %isscaf); 
#   $gfar = scalar(keys %isfar);   ## these count uniq gene ids now
#   $gnear= scalar(keys %isnear); 
#   $gsame= scalar(keys %issame);  ## but gene can be in >1 class here
#   $gpairs=scalar(keys %genepairs);
#?  $prev =scalar(keys %nearrev); $prev= sprintf("%.3f", $prev/$pnear) if($pnear);

  my @grp= sort keys %distclass;
  foreach my $grp (@grp) {
  $gscaf= scalar(keys %{ $distclass{$grp}{'scaf'} } );
  $gfar = scalar(keys %{ $distclass{$grp}{'far'} } );
  $gnear= scalar(keys %{ $distclass{$grp}{'near'} } );
  $gsame= scalar(keys %{ $distclass{$grp}{'same'} } );
 
  my $modelerrs= join ",", 
    map{ my $c= $modelerrs{$grp}{$_}; "$_:$c" } sort keys %{ $modelerrs{$grp} };
    #map{ my $c=scalar(keys %{$modelerrs{$grp}{$_}}); "$_:$c" } 
  
  ## need $grp here :
  my $ngenetotal = $ngene || scalar( grep /^$grp/, keys %genebe ); # count all mRNA/prot IDs; fixme for full AA subsets
  $gpairs=scalar( grep /^$grp/, keys %genepairs);
  my $noloc= $noloc{$grp};
  
  print "source=$grp; " if $grp; # was $source
  print "ngenes=$ngenetotal; npaired=$gpairs; ngroup=$nchain; noloc=$noloc; skip=$nskip; ";
  print "clusters=$clustersize..$clustertop; " if($clustersize);
  print "clusters=all; " unless($clustersize);
  print "neardist=$NEARDIST\n";
  print "modelerrs=$modelerrs\n" if($modelerrs);
  print " gscaf=$gscaf; gfar=$gfar; gnear=$gnear; gsame=$gsame\n";  
  
  # want range of near steps: 1000, .., 5000, 10000, 15000 .. 50000 , far, scaf
  foreach my $kb (sort {$a<=>$b} keys %{$nearsteps{$grp}}) {
    my $cnear =scalar(keys %{$nearsteps{$grp}{$kb}}); 
    print " kb$kb=$cnear,"
    }
  print "\n";
  
  if(SHOWREV) {  
    foreach my $kb (sort {$a<=>$b} keys %{$revsteps{$grp}}) {
      my $crev =scalar(keys %{$revsteps{$grp}{$kb}}); 
      print " rv$kb=$crev," 
      }
    print "\n";
    foreach my $kb (sort {$a<=>$b} keys %{$revsteps{$grp}}) {
      my $cnear =scalar(keys %{$nearsteps{$grp}{$kb}}) || 1; 
      my $crev =scalar(keys %{$revsteps{$grp}{$kb}}); 
      my $rp= sprintf("%.3f",$crev/$cnear);
      print " revp$kb=$rp," # or $crev? both
      }  
    print "\n";
    }
  print "\n";
  } #... end $grp

  # print $nskip, ... 
  if($shownear) { print "  near-genes:\n"; # change to all dupl pairs: ALLPAIRSCORE == showall
    map{
      my $v= $nearpair{$_}; my $pr=1;
#       # NO: filter out cross-pred pairs if we have MIXGROUPS
#       if($MIXGROUPS) {
#         my($gi,$gj)= $_ =~ m/^(\D+)\d+.(\D+)/;
#         $pr=0 unless($gi eq $gj);
#         }
      print $_,"\t$v\n" if($pr);
      } sort keys %nearpair;  
    }
 
}


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

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

sub _isoverlap { 
  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
  #? is this bad?# return 0 if($rev);
  
  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 $over;
}

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

sub _sizediff {
  my($tb,$te,$to, $qb,$qe,$qo, $checkswap)= @_; # test Q / T
  my $tlen= abs($te - $tb) || 1;
  my $qlen= abs($qe - $qb) || 1;  
  if($checkswap and $qlen > $tlen) { #swap
    ($qlen,$tlen)=($tlen,$qlen);
    ($tb,$te,$to, $qb,$qe,$qo)= ($qb,$qe,$qo, $tb,$te,$to);
    }
  return $qlen / $tlen;  
}

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
  # if ($tb <= $qe && $te >= $qb) { return (0,$rev); } # or -?

  #? what of _isinside ?
  my $over= ($tb <= $qe && $te >= $qb) ? 1 : 0;
  
  if($over) {
    my $isin= _isinside($tb,$te,$to, $qb,$qe,$qo, 1);
    return (1,$rev) if($isin); # special case
    }
    
  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 (2,$rev) unless($over); # special case 2; treat different from _isinside?
  }

  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);
    }
}


use constant SKIP_ALTTR => 1; # this alttr filter problematic w/ multiple predictors

sub genespans {
  my (%locs,%gene,%altr,$gid);
  %altr=(nothing => 1);
  return (nothing => 1) if(SKIP_ALTTR);
  
  foreach my $id (sort keys %genebe) {
    my ($pgrp)= $id=~m/^(\D+)/; 
    my $ref= $generef{$id}; $ref.=".".$pgrp;
    my ($tb,$te,$to)= @{$genebe{$id}}; # stick ref here ?
    my @bins= ($tb/$BINSIZE .. $te/$BINSIZE);
    foreach my $ib (@bins) { push @{$locs{$ref}{$ib}}, $id; } 
    }
  
  # *** this is not good for multiple predictors, want only same-predict filter
  # *** tack group onto ref and is ok
  # ??? ignore this alttr filter and use "same/overlap/inside" class for distance test?
  
  foreach my $ref (sort keys %locs) {
    foreach my $ib (sort{$a<=>$b} keys %{$locs{$ref}}) {
      my @ids= @{$locs{$ref}{$ib}}; # all are alt-tr? NO, check $genebe/tb,te overlap
      my %ids= map{$_,1}@ids; 
      @ids= sort keys %ids;
      next if(@ids<2);
      my($gb,$ge,$go,$ggroup)=(-1,0,"+");
      my @altids=();
      ## fixme: dont mark as altr/altids unless same predictor ???
      ## at least overlap of diff-predictors cant be used to extend genespan
      map{ 
        my($tgroup)= $_ =~ m/^(\D+)/;
        my($tb,$te,$to)= @{$genebe{$_}};
        
        if($gb==-1) { ($gb,$ge,$go,$ggroup)= ($tb,$te,$to,$tgroup); push(@altids,$_); }
    #??    elsif($tgroup ne $ggroup) { next; } #*** dont mix predictors in altid set
        elsif ( _isoverlap($gb,$ge,$go, $tb,$te,$to)
          # and not _isinside($gb,$ge,$go, $tb,$te,$to) #? treat as alttr here/same predictor group?
          ) { 
          #? check for _isinside() where big-gn span >> small-gn span
          # have some aberrant predictions with huge span, many genes inside
          $gb= $tb if($tb<$gb);
          $ge= $te if($te>$ge);
          push(@altids,$_);
          }
        } @ids;
      if(@altids>1) {
        map{ $altr{$_}=1; $genebe{$_}= [$gb,$ge]; } @altids;
      }
      
    }
  }
  return %altr;
}


__END__

=item old subs


sub pairs2chains {

  my @aid= sort keys %pairs; 
  my $ng= scalar(@aid); 
  my ($na, $nb, %nc, @chains, %didab,);
  
  foreach my $aid (@aid) { 
    my @bid= sort keys %{$pairs{$aid}}; 
  
    my @cid=(); 
    BLIST: foreach my $bid (@bid) { 
    # next unless($removesubsets or !$didab{$aid}{$bid});
    if($noreciprocal or $pairs{$bid}{$aid} ) { 
      push(@cid, $bid); 
      $didab{$aid}{$bid}= $didab{$bid}{$aid}=1;  # check recip match
      }
     } 
    
    my %bid= map{ $_=>1; } $aid,@cid; 
    @bid= sort keys %bid;
    
    if ($removesubsets and @bid>1) {
      my $nid= @bid;
      foreach my $chain (@chains) {
        my $did=0; map{ $did++ if $bid{$_};  } @$chain;
        @bid=() and last if($did >= $nid);
        }
    }
   
    if (@bid>1) { 
      push(@chains, \@bid);
      # print join("|",@bid),"\n"; 
      $na++; $nb += @bid; @nc{@bid}=1; 
      } 
    } 
  
  my $nc= scalar(keys %nc);  
  my $nself= scalar(keys %nself);
  my $nsing= $nself - $ng;
  print "# nself=$nself; npaired=$ng; nsingle=$nsing; duplicates: ngroup=$na; nuniqgene=$nc; ngroupedgene=$nb; rmsubsets=$removesubsets \n"; 
  return \@chains; # hash by $aid instead ?
}


sub chains2near { # blast/blat idchains
  my ($chains)= @_;
  # my ();
  
  foreach my $chain (@$chains) {
  
  my @d= @$chain; # split /\|/; 
  if(/[a-z]{4}_GLEAN_/){map{s/^[a-z]{4}_//}@d;} # fix for _GLEAN_ id change .gff vs .aa
  
  if($clustersize) {
    $nskip++ and next if($clustersize>0 and @d < $clustersize);
    $nskip++ and next if($clustersize<0 and @d > -$clustersize);
    if($clustertop>0) {    
      $nskip++ and next if( @d > $clustertop);
    }
  }
  
  $nchain++;  # @nc{@d}=1; 
  my($nc,$nscaf,$nfar,$nnear,$nsame)=(0)x9; 

  for my $i (0..$#d) { 
    my $di= $d[$i];  $noloc++ and next unless($genebe{$di}); # count skips

    for my $j ($i+1..$#d) { 
      my $dj= $d[$j];  next unless($genebe{$dj});
      
      next unless($dupgeneok or not($genepairs{$di}{$dj} or $genepairs{$dj}{$di}));
        # no effect here
      $genepairs{$di}{$dj}++; # store distance instead of count for pairs?
      $genepairs{$dj}{$di}++;  
      # ^ note we have duplicate IDs across idchain lines; adjust for that how?
      
      if($generef{$di} ne $generef{$dj}) { 
        $nscaf++; $isscaf{$di}++; $isscaf{$dj}++; 
        $scafpair{"$di.$dj"}++; $scafpair{"$dj.$di"}++; 
        } 
      else {
        my ($dist,$rev)= _mindistance( @{$genebe{$di}}, @{$genebe{$dj}} );
        if ($dist<=0) { $issame{$di}++; $issame{$dj}++; }
        else {
          if($dist>$NEARDIST) { 
             $isfar{$di}++; $isfar{$dj}++; 
             $farpair{"$di.$dj"}++; $farpair{"$dj.$di"}++;
             $nearsteps{999}{$di}++; $nearsteps{999}{$dj}++;
             if($rev) { $revsteps{999}{$di}++; $revsteps{999}{$dj}++; }
             } 
          else { 
            $isnear{$di}++; $isnear{$dj}++; 
            my $distep= 1 + int( $dist / 1000);
            $nearsteps{$distep}{$di}++; $nearsteps{$distep}{$dj}++;
            $nearpair{"$di.$dj"}++; $nearpair{"$dj.$di"}++; 
            if ($rev) { 
              $nearrev{"$di.$dj"}++;  $nearrev{"$dj.$di"}++; 
              $revsteps{$distep}{$di}++; $revsteps{$distep}{$dj}++;
              }
            } 
          }
        }
                
      } # dj..
    } # di..
  } # chains
}


=cut