#!/usr/bin/perl
# protnear.perl

# tandem gene analysis by proteins (07jul13)
# see below, ~/Desktop/dspp-work/daphwork/daph-tandemgene.info
# FIXME: 1. find max gene-span for alt-tr to filter all overlaps
# FIXME: 2. filter out TE-repeat overlaps (from mRNA gff)
# fixed: 3. had bad $rev value, stats (perl oddity??)

# set blo=dpulexgnomon.blat50
# set bgff=$sc/dpulex4/dpulex_jgi060905_Gnomon.gff.gz
# gzgrep mRNA $bgff | cat - $blo.idchains | perl $td/protnear.perl

my $NEARDIST= 15000; # common with tandy exon choice
my $SAMEDIST= 10;  #  do _isoverlap; but use to screen odd cases 
our $BINSIZE  = 1000 ; #was# 5000;
my(%gr, %gm, %gbe, %galtids, %genepairs, %isscaf,%isfar,%isnear,%issame,
   %nearsteps, %revsteps, %nearpair, %nearrev, %farpair, %scafpair);
my $shownear=0;
my $mrnatype= $ENV{mrnatype} || "mRNA";
my $ngene= $ENV{ngene} || 0;
my $dupgeneok= 1;
my $pctover=0;
my $clustersize= 0; my $clustertop=0;
my @clusize;
my $source="";

use Getopt::Long;

my $optok= GetOptions(
  "mrnatype=s", \$mrnatype, 
  "ngene=s", \$ngene, 
  "source=s", \$source, 
  "dupgeneok!", \$dupgeneok,
  "pctoverlap=i", \$pctover,   
  "NEARDIST=i", \$NEARDIST,  #"NEARSTEPS=i", \$NEARSTEPS,  
  "clustersize=i", \@clusize,  # \$clustersize,  
  "shownear!", \$shownear,  
  );
  
die "usage: cat genes.gff genes.blastp.idchains | perl protnear.perl [ options ] \
  opts: -mrnatype mRNA  -ngene 20000 -nodupgene -neardist $NEARDIST
        -pctoverlap=20  : separate 'same' genes/alttr from near genes
        -shownear : id list for near pairs
        -cluster 10 -cluster 20 : range of cluster sizes
" unless($optok);

$pctover= $pctover/100.0 if($pctover);
if(@clusize) {
  $clustersize= $clusize[0];
  $clustertop= $clusize[1] || 0;
}

=item fixme: subset chains

WBGene00000003|WBGene00000004|WBGene00000006|WBGene00000007|WBGene00000008|WBGene00000009|WBGene00000010
WBGene00000004|WBGene00000006|WBGene00000007|WBGene00000008|WBGene00000010
WBGene00000005|WBGene00000006|WBGene00000007|WBGene00000008|WBGene00000009|WBGene00000010
            ^^^^ the following are essentially subsets of 1st;
            fix to remove these from .idchains
=cut
my ($noloc,$nskip)=(0) x 9;

while(<>){
chomp; 
if(/\t$mrnatype/){ # gff for loc,ids
  my($gr,$gs,$gt,$tb,$te,$gp,$to,@x)=split /\t/; my($d)=m/ID=(\w+)/; 
  $gr{$d}=$gr;  # $gm{$d}= int(($tb+$te)/2); 
  $gbe{$d}= [$tb,$te,$to];
  $source ||= $gs;
  } 
elsif(m/\t/) {} 
elsif(m/\|/) { # blast/blat idchains
  defined %galtids or %galtids= genespans(); 
  my @d=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($gbe{$di}); # count skips

    for my $j ($i+1..$#d) { 
      my $dj= $d[$j];  next unless($gbe{$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 pair?
      $genepairs{$dj}{$di}++;  
      # ^ note we have duplicate IDs across idchain lines; adjust for that how?
      #  no change needed for below counts as 1st genepair near/far is all it uses
      # ? look at idchains using only best match? >> odd result from idchains bestmach
      # look at near/far per contiguous-megabase of genome to account for assembly qual
      # big increase for daphnia w/ scafpair count: many dupl across scaffolds?
      
      #? $nc++; 
      
      if($gr{$di} ne $gr{$dj}) { 
        $nscaf++; $isscaf{$di}++; $isscaf{$dj}++; 
        $scafpair{"$di.$dj"}++; $scafpair{"$dj.$di"}++; 
        } 
      else {
        my ($dist,$rev)= _mindistance( @{$gbe{$di}}, @{$gbe{$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}++;
              }
            } 
          }
        }
                
      } 
    } 
    
  $gmemb+=@d; # this counts some genes > 1 time
#   $gscaf+=scalar(keys %isscaf); $gfar+=scalar(keys %isfar);   ## these are not counts of uniq gene ids; overlap among chains
#   $gnear+=scalar(keys %isnear); $gsame+=scalar(keys %issame); 
  } 
}

OUTPUT: 
  outdist();
  # out1();
  
sub outdist {
  my($gscaf,$gfar,$gnear,$gsame)=(0)x9;
  my($pscaf,$pfar,$pnear)= (0)x9;
  
# 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);
  
  $ngenetotal = $ngene || scalar(keys %gbe); # count all mRNA/prot IDs; fixme for full AA subsets
  print "source=$source; " if $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 " 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) {
    my $cnear =scalar(keys %{$nearsteps{$kb}}); 
    print " kb$kb=$cnear,"
    }
  print "\n";
  foreach my $kb (sort {$a<=>$b} keys %revsteps) {
    my $crev =scalar(keys %{$revsteps{$kb}}); 
    print " rv$kb=$crev," 
    }
  print "\n";
  foreach my $kb (sort {$a<=>$b} keys %revsteps) {
    my $cnear =scalar(keys %{$nearsteps{$kb}}) || 1; 
    my $crev =scalar(keys %{$revsteps{$kb}}); 
    my $rp= sprintf("%.3f",$crev/$cnear);
    print " revp$kb=$rp," # or $crev? both
    }
  
  print "\n";
  
  # print $nskip, ... 
  if($shownear) { print "  near-genes:\n"; map{print $_,"\n"} sort keys %nearpair;  }
 
}


sub out1 { 
  ## for gscaf/gfar/gnear/gsame change to count uniq IDs in each class, not all per chain line (w/ overlap)
  ## provide  result / total-genes (including not in blat/blast result list; from gff)
  #? $nc=scalar(keys %nc); # this is gene IDs returned in blat/blast output < total genes
  my($gscaf,$gfar,$gnear,$gsame)=(0)x9;
  my($pscaf,$pfar,$pnear)= (0)x9;
  
# 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);
    # ^ should instead look at counts/gene or genepair counts ??
    # near, far, scaf pairs counts ?
  $pscaf =scalar(keys %scafpair); 
  $pfar = scalar(keys %farpair);   
  $pnear =scalar(keys %nearpair); 
  $prev =scalar(keys %nearrev); $prev= sprintf("%.3f", $prev/$pnear) if($pnear);
  
  $ngenetotal = $ngene || scalar(keys %gbe); # count all mRNA/prot IDs; fixme for full AA subsets

  print "source=$source; " if $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 " gscaf=$gscaf; gfar=$gfar; gnear=$gnear; gsame=$gsame\n"; # gmemb=$gmemb; 
  print " pscaf=$pscaf; pfar=$pfar; pnear=$pnear; prev-near=$prev \n"; 

  ($gscaf,$gfar,$gnear,$gsame)= 
    map{ sprintf("%.3f", $_/$ngenetotal); } ($gscaf,$gfar,$gnear,$gsame); 
  print "  per-gene: gscaf=$gscaf; gfar=$gfar; gnear=$gnear; gsame=$gsame\n"; 
  ($pscaf,$pfar,$pnear)= 
    map{ sprintf("%.3f", $_/$ngenetotal); } ($pscaf,$pfar,$pnear); 
  print "  per-gene: pscaf=$pscaf; pfar=$pfar; pnear=$pnear; \n"; 

  # print $nskip, ... 
  if($shownear) { print "  near-genes:\n"; map{print $_,"\n"} sort keys %nearpair;  }
}


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

# add pctoverlap option: want to keep or test cases of minor overlap (e.g. UTRs)?


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
  #?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
  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 < 0.50); # 0.33;  #? use pctover here?
    }
  return 0;
}

sub _mindistance { 
  my($tb,$te,$to, $qb,$qe,$qo)= @_;
  ##my $rev= ($to eq "-" or $qo eq "-") and ($to ne $qo); ## BAD
  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 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);
    }
}

#use constant SKIP_ALTTR => 1;

sub genespans {
  my (%locs,%gene,%altr,$gid);
  %altr=(nothing => 1);
  
  foreach my $id (sort keys %gbe) {
    my($tb,$te,$to)= @{$gbe{$id}};
    my $ref= $gr{$id};
    my @bins= ($tb/$BINSIZE .. $te/$BINSIZE);
    foreach my $ib (@bins) { push @{$locs{$ref}{$ib}}, $id; } 
    }
  #return (nothing => 1) if(SKIP_ALTTR);
  
  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 $gbe/tb,te overlap
      my %ids= map{$_,1}@ids; @ids= keys %ids;
      next if(@ids<2);
      my($gb,$ge,$go)=(-1,0,"+");
      my @altids=();
      map{ 
        my($tb,$te,$to)= @{$gbe{$_}};
        if($gb==-1) { ($gb,$ge,$go)= ($tb,$te,$to); push(@altids,$_); }
        elsif ( _isoverlap($gb,$ge,$go, $tb,$te,$to)
          and not _isinside($gb,$ge,$go, $tb,$te,$to)
          ) { 
          #? 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; $gbe{$_}= [$gb,$ge]; } @altids;
      }
      
    }
  }
  return %altr;
}


__END__

# tandem gene analysis by proteins (07jul13)
# find locations of genes/group in .idchains, list near,far,diff-scaff distance counts
## see ~/Desktop/dspp-work/daphwork/daph-tandemgene.info