#!/usr/bin/perl
# geneclust.perl
# usage: cat scaffold.tandy.gff | perl geneclust.perl -ngene=3
# old usage: grep match scaffold.tandy.gff | perl -n geneclust.perl 3
# count of rich regions: with 3+ genes, each with large size (measure against gene size?)
# ?? check, remove overlapping regions, or matches, of lower quality? -- how many are there?

## see this result
## cat $em/dmel4d/*/*_exons_tandy6j.gff | grep HSP | perl -ne\
# '@v=split;  print join("\t",$v[1]),"\n" if(1);' | sort | uniq -c |\
# perl -ne'($n,$t)=split; next unless($t=~/tandy/); $tn{$t}=$n; \
# END{ $ne=$tn{"tandy.equal"}; @t=sort keys %tn; foreach $t (@t) { printf "%5d\t%5.3f\t%s\n",$tn{$t},$tn{$t}/$ne,$t;}}'


use Getopt::Long;

my $SKIPOVERLAP=1;
my $nc= $ENV{nc} || 2;
my $optok= GetOptions(
  "ngene=i", \$nc,
  );

my %b;  
while(<>){
  next unless(/^\w/);
  chomp; 
  my($ref,$src,$typ,$tb,$te,@gff)= split"\t";  
  next unless($typ =~ /^(match|HSP)/);
  
  # src should be 'tandy.near,far,...', tandyweak.near,far,.. also attr tclass=near,far;
  # my ($tand,$eqclass)= split /\./, $src;
  my($eqclass)= ( m/tclass=(\w+)/ ) ? $1 : ""; #unless($eqclass)
  my($terepeat)= ( m/terepeat=(\w+)/ ) ? $1 : 0;
  
  # per gene, gid
  # m/info=([^;]+)/ or next; 
  if( m/info=([^;]+)/) { # == $typ eq "match"
    my $b=$1; my @b=split",",$b; %b= map { split ":" }@b;
  }
  
  my $ishsp= ($typ eq "HSP");
  my($sid)=m/^(\w+)/;  $sid=~s/scaffold_/s/;
  my($cid)= m/cid=(\d+)/;  my($gx)= m/gid=(\d+)/;
  my $scid= "$sid.c$cid";  my $gid= "$scid.g$gx";
  
  $criteria="fullgenes>0 or pctalign>33 and not(overlap or TErepeat or lowquality)"; # nexons > 1 ??
  next if ($b{lowquality}>0 or $terepeat>0 or $src =~ /tandyweak/);
  next unless($b{fullgenes} > 0 or $b{pctalign}>33 );
  
  my $bin= 10000 * int(($tb+$te)/20000);
  if (!$ishsp && $SKIPOVERLAP && $overlapped{$ref.$bin}++) {
    #? want to swap earlier for better gene? how?  
    my @oldloc= @ { $saveloc{$ref.$bin} || [] };
    my $ov= 0;
    foreach my $rloc (@oldloc) {
      $ov=1 and last if _isoverlap( $rloc->[1], $rloc->[2], $tb, $te);
      }
    $skipov++ and next if($ov);
    }
  push( @{$saveloc{$ref.$bin}}, [$ref,$tb,$te]);
  
  
  # check, filter by nexons, pctalign, fullgenes ?
  $gs{ngene}{$scid}++ unless($ishsp);
  $gs{geneid}{$scid}{$gid}++ unless($ishsp);
  
  $gs{nexon}{$scid}++ if($ishsp);
  # this now counts/HSP
  if($eqclass && $ishsp) { $gs{eqnear}{$scid}++ if($eqclass =~ /^near$/); }
  #else { $gs{eqnear}{$scid}++ if($b{eqnear} > 0); }
  
  if($eqclass && $ishsp) { $gs{eqfar}{$scid}++ if($eqclass =~ /^far$/); }
  #else { $gs{eqfar }{$scid}++ if($b{eqfar } > 0); }

  if($eqclass && $ishsp) { $gs{eqsame}{$scid}++ if($eqclass =~ /^equal$/); }
  #else { $gs{eqsame}{$scid}++ if($b{eqsame} > 0); }

  $gs{duponly }{$scid}++  #?? fixme / drop?
    if(!$ishsp && $b{eqsame} == 0 && ($b{eqnear} > 0 || $b{eqfar } > 0));
}


END{ 
  my (%ns, %ts, %ngn);
  
  @gn= sort keys %{$gs{ngene}};
  @g3= grep { $gs{ngene}{$_} >= $nc } @gn;
  my $ngnt=0;
  
  foreach $gn (@g3) {
    $nonly{eqnear}++ if ($gs{eqnear}{$gn} > 0 && ! $gs{eqfar }{$gn});
    $nonly{eqfar }++ if ($gs{eqfar }{$gn} > 0 && ! $gs{eqnear}{$gn});
    #my $ngn= scalar(keys %{$gs{geneid}{$gn}}); #??
    #$gs{ngeneid}{$gn}= $ngn;
    $ngnt += $gs{ngene}{$gn};
    }
  
  print "# criteria: $criteria\n";
  print "# gene cluster (n>=$nc) counts\n";
  print "total regions: ",scalar(@g3),"\t";
  foreach $eq (qw(eqnear eqfar eqsame duponly)) {
    foreach $gn (@g3) {
      $ns{$eq}++ if ($gs{$eq}{$gn} > 0);
      $ts{$eq}+= $gs{$eq}{$gn};
      $ngn= $gs{ngene}{$gn}; # = $gs{ngeneid}{$gn};
      $ngn{$eq}+= $ngn if ($gs{$eq}{$gn} > 0); #???
      }
    print "$eq: ",$ns{$eq},"\t";  
    }
  print"\n"; 
  print "genes/region: $ngnt\t";
  foreach $eq (qw(eqnear eqfar eqsame duponly)) {
    $ns= $ngn{$eq} || 0; 
    print "$eq: ",$ns,"\t";  
    }
  print"\n"; 
  print "only   class regions: \t";
  foreach $eq (qw(eqnear eqfar)) {
    $ns= $nonly{$eq}; 
    print "$eq: ",$ns,"\t";  
    }
  print"\n"; 
  print "average HSPs/region: \t";
  foreach $eq (qw(eqnear eqfar eqsame duponly)) {
    $ns= $ns{$eq} || 1; 
    $ts= sprintf("%.1f",$ts{$eq} / $ns);
    print "$eq: ",$ts,"\t";  
    }
  print"\n"; 
  print "# Overlapping regions skipped: $skipov\n";
  }

sub _isoverlap { 
  my($gb,$ge, $qb,$qe, $qid)= @_;
  return ($gb <= $qe && $ge >= $qb) ? 1 : 0;
}
 
# ID=td_c1g2.GI_NCBI_GNO_32126697;
# tclass=equal;cid=1;gid=1586;
# nexons=1;altexons=1;info=eqfar:0,eqnear:0,eqsame:0,
#   fullgenes:2,geneids:3,maxgenes:1,methods:3,nexons:1,pctalign:5,repeats:0,tandems:0;
# bestids=GI_NCBI_GNO_32126697