#!/usr/bin/perl
# tandytab2.perl
# collate results of tandynear2/d*_*-pexons.tandynear41mb.txt, other for R stats table

use strict;
use warnings;

my @in= @ARGV; # need fnames for species, group

my @idstats= qw( found nxfound nexon   align   cdslen);
my($spp,$ggroup,@flds,@sameflds,@idflds,%stats,%info);

foreach my $in (@in) {

  if($in =~ m/pexons.tandynear/) {
    ($spp,$ggroup)= $in =~ m,([a-zA-Z]+)_(\w+)-pexons,;
    die "#bad spp=$spp; group=$ggroup from $in \n" unless($spp and $ggroup);
    
    $info{$spp}{$ggroup}{file}= $in;
    open(IN,$in);
    while(<IN>){ 
      chomp; 
      if(s/^#t\s+//) { 
#t Group                Stat    Same    Inside  Near15k Near30k Near45k Far
#t all                  count   31664   7678    3146    174     96      1014
#t all                  freq    1.000   0.242   0.099   0.005   0.003   0.032
#t all                  found   31664   84      189     15      15      34
        my($gall,$stat,@v)=split; 
        if(s/^Options:\s*//) {  $info{$spp}{$ggroup}{options}= $_;  }
        elsif($stat =~ /^Stat/) { @flds=@v; @sameflds= map{ $_."_found"; } @flds; }
        elsif($stat =~ /^count/){
          map{ my($k,$v)= ($_,shift @v); $stats{$spp}{$ggroup}{$k}=$v; } @flds;
          }
        elsif($stat =~ /^found/){
          map{ my($k,$v)= ($_,shift @v); $stats{$spp}{$ggroup}{$k}=$v; } @sameflds;
          }

        }
        
      elsif(m/^\w/) { # geneid stats .. collect and add sums? to %stats
# # id quality table for near-duplicate genes, class=matched
# group   gid     loc     near    found   anyfnd  nxfound nexon   align   cdslen  marks
## group: all; Near15k=1263
# GM_NCBI_GNO_    32000328        scaffold_16:53619-55695 Near15k 0       0       1       5       167          1749    dmelhsp=1

        my @v=split; 
        if(/^group/){ @idflds= @v; }
        else {
          my %v= map{ $_, shift(@v) } @idflds;
          foreach my $idk (@idstats) {
            $stats{$spp}{$ggroup}{$idk} || 0;
            $stats{$spp}{$ggroup}{$idk} += $v{$idk};
            }
          # collect stats for found,nxfound,nexon,align,cdslen
          # By pred-group, near state
          }
        }
      }
  }
  
}

# collectstats();
dumpstats();

sub dumpstats {
  my $didhead;
  
  my @info;
  my @spp= sort keys %stats;
  foreach my $spp (@spp) {
    my @ggroup= sort keys %{$stats{$spp}};
    foreach my $gr (@ggroup) {
    
      push @info, join ", ", $spp,$gr, 
        map { $_ ."=". $info{$spp}{$gr}{$_} } sort keys %{ $info{$spp}{$gr} };
      
      my @sfld= sort keys %{ $stats{$spp}{$gr} }; #%sr; ## %{$stats{$spp}{$gr}};
      my @vals= map{ $stats{$spp}{$gr}{$_} } @sfld;  

      print join("\t","species","ggroup",@sfld),"\n" unless($didhead++);
      print join("\t",$spp,sprintf("%10s",$gr),@vals),"\n";
      }
  }
 
  print "\n#info: \n"; foreach (@info) { print "# $_\n"; }
}

__END__

R stats:

setwd("~/Desktop/dspp-work/genomesoft/tandy/")

pe <- read.table("pexons.tandynear4.txt", header=T, strip.white=T)
# pe <- pe1

ggroup <- "Gnomon"
pelev <- grep("NCBI",as.character(pe$ggroup),value=T) 

ggroup <- "GleanR"
pelev <- grep("GLEANR|JGI|FBtr",as.character(pe$ggroup),value=T) 

  sppset <- c("dmel","dsec","dyak","dere","dana","dpse","dwil","dgri","dappulx")
  spat <- paste(sppset,collapse="|")
  
  pefac <- factor(pe$ggroup, levels=pelev, ordered=T,exclude=NULL)
  
  pea <-  pe[!is.na(pefac),c(-1,-2)]
  rown <- pe[!is.na(pefac),1] # species
  rownames(pea) <- rown
  pea <- pea[sppset,]

  predtitle <- paste("Pred. exons/species, predictor",ggroup)
  ylabel <- "exon count"

  pea$align <- pea$align/pea$cdslen
  pea$nxfound <- pea$nxfound/pea$nexon
  
  tparm <- c( 6,8,10) # near_found ,15 = found ~= nearcount; 13,17 == align,nxfound   ; 2,4, = far, inside
  ymax <- 0.04 # w/ above as relative Same, max is 0.03+ for Gnomon, 0.015 for GleanR

  tparm <- c( 5,7,9) # near ; lacks dyak odd spike > dyak tandems were found at higher rate than others?
  ymax <- 0.15 # w/ above as relative Same, max is 0.03+ for Gnomon, 0.015 for GleanR
  
  pea[,tparm] <- pea[,tparm] / pea$Same_found
  
  tcols <- c("lightblue", "mistyrose", "lightcyan", "lavender", 2:20) [1:length(tparm)]
  
  bpx <- barplot( as.matrix(t( pea[,tparm])),ylab=ylabel, main=predtitle,
    ylim=c(0,ymax), beside=T, legend=F, col = tcols ) #  

  legend( x="topleft", colnames(pea[,tparm]), lwd=3, col=tcols)


#........
  # tparm <- grep("_found",as.character(colnames(pea))) #   2  4  6  8 10 12

# ymax <- max(pea[,tparm])

# layout(as.matrix(1))
# par(mar=c(4,4,3,2)+0.1)

# sort by sppset

> colnames(pea)
 [1] "Far"           "Far_found"     "Inside"        "Inside_found"  "Near15k"      
 [6] "Near15k_found" "Near30k"       "Near30k_found" "Near45k"       "Near45k_found"
[11] "Same"          "Same_found"    "align"         "cdslen"        "found"        
[16] "nexon"         "nxfound"