#!/usr/bin/perl

# genome feature densities / variation (hot spots?)
# count per base window (10kb? 50kb?)
#  cdsbases; genespan; ncds, ngenes, intergene dist; other annot features like prot dupls.

#cat xxx.gff | sort -k1,1 -k4,4n -k5,5nr | grep '^scaffold_..  ' | perl genomedens.pl
#be careful of mixed types: mRNA, CDS sorted by loc will mess up overlap check

BEGIN{ $bins=10000; }

while(<>){
  ($r,$s,$t,$b,$e)=split; 
  $r=~s/scaffold\D*//i; # chr?
  ($lr,$lb,$le)= @{ $last{$t} || [-1,0,0] };
  $ov=($lt eq $t and $r eq $lr and $b < $le and $e > $lb);
  ($lb,$le)=(0,0) if($lr ne $r);
  
  # for CDS, mRNA, ...
  unless($ov)  { 
    $nb=1+$e-$b; 
    $ib=$b - $le - 1; # negative ok?
    $bn= $bins * int(($b+$e)/(2*$bins));
    #? ($id)=m/(?:ID|Parent)=([^;]+)/;

    $nc{$r}{$bn}{$t."ib"} += $ib; # if CDS && id == lid; i.e. introns only
    $nc{$r}{$bn}{$t."nb"} += $nb; 
    $nc{$r}{$bn}{$t}++; 
    $tp{$t}++;
    $tp{$t."nb"}++;
    $tp{$t."ib"}++;
    
    $at="dupp"; ($v)=m/paralog=Omcl\d+,(\d+)/; $v||=0;
    if($v) { $nc{$r}{$bn}{$at}++; $tp{$at}++; }
    $at="dupx"; ($v)=m/tandy=(\w+)/; $v||=0;
    if($v) { $nc{$r}{$bn}{$at}++; $tp{$at}++; }
    }
  
  ($lr,$lb,$le,$lt)=($r,$b,$e,$t);
  $last{$t}=[$r,$b,$e];
}

END{ 

# print "# $eg EST tandem clusters\nscaf\tlocat\tnEST\tnTand\ttandem_clusters\n"; 
@tp= sort keys %tp; # these are columns!
print join("\t","scaf","loc",@tp),"\n";

@r= sort{$a<=>$b}keys %nc; # numeric scaffold# 
foreach $r (@r) {  
  @b= sort{$a<=>$b}keys %{$nc{$r}}; 
  foreach $b (@b) {
    print "$r\t$b";
    foreach $t (@tp) { $nc= $nc{$r}{$b}{$t} || 0;  print "\t",$nc; }
    print "\n";
    } 
  } 
}

# > ! genomedens.ltab