#!/usr/bin/perl
# R stats convert
# tandynear2rplot.pl

use strict;
# use warnings;
use Getopt::Long;


my $title="Predictor x Tandem genes, tBLASTn";
my $options="";

my @gvals= qw(npaired gscaf gfar gnear); #  edoubled eskipover : drop?
my @skeys =  qw( ngenes npaired gsame gscaf gfar gnear); # add below:  kb1 kb2 kb3
my @sppkeys= qw(Dsec Dsim Dyak Dere Dana Dpse Dwil Dvir Dmoj Dgri);
my %spcode=( CG => "Dmel",
  GM => "Dsec", GD => "Dsim", 'GE' => "Dyak", GG => "Dere", 
  GF => "Dana", GA => "Dpse",  Gfixme => "Dper",
  GK => "Dwil", GI => "Dmoj", GJ => "Dvir",  GH => "Dgri",
  );

use constant OLDROUT => 0;
       
my $kbrange= "1..15";
my $speciesin= "none";
my (%vals, $species, %savetab, $tabsave, $input, $asource, $athead);
my $tabsave=1;

my $optok= GetOptions(
  "input=s", \$input, 
  "title=s", \$title, 
  "options=s", \$options, 
  "kbrange=s", \$kbrange, 
  "species=s", \$speciesin, 
  "tabsave!", \$tabsave,
);

die "usage: tandynear2rplot -input=file|stdin -species=$species -kbrange=$kbrange < in > out
e.g. gzcat blastnear/*blastnear22c.gz | grep -v '^(G.|d...)_' | tandynear2rplot 
or,  head d???/*snapdmel.blastnear | tandynear2rplot.pl > snapdmel.blastnear.rdata 
" unless($optok);

my @kbrange= eval $kbrange; # (1..15); # (1..20);
push @skeys, map{"kb$_"} @kbrange;

my $inids=0;

print "# R-data: source('this.rdata')\n";
# if(OLDROUT){
# print "dl <- c() \n";
# }

my $inh=*STDIN;
if($input && -f $input) { open(IN,$input); $inh=*IN; }

while(<$inh>){  # want (<>) ??
  next if($inids and /^(G\w|d\w\w\w)_/);
  s/^\s+//;
  $inids=0 if(/=|#/);
  $inids=1 and next if(/near-genes:/);
  
  if(/^#/){ 
    $species= $speciesin unless($athead); 
    if(/Options: (.+)$/) { my $op=$1; $options .= ", $op" unless($options =~ /$op/); }
    # Options: min_bitscore=150,min_e_value=1e-10,pct_align=0.3,pct_bitscore=0,pct_overlap=0.2
  } 
#   elsif(m,^==\> ([a-zA-Z]+),) { # ==> dpse/dpse_dgilsno.blastnear <==  head output
#      my $sptag=$1;
#      # if(!$species or $speciesin =~ /^none$/)
#      $species= $sptag;
#      $athead++;
#   }
  unless(/^\w/) {
    if(%vals and $vals{source}) { 

      #? collect species from $sptag
      my $sptag=""; if( $vals{source} =~ s/^(G\w|d\w\w\w)_//) { $sptag=$1; }
      if($sptag and (!$species or $speciesin =~ /^none$/)) { $species= $sptag; }
      
      $species = $spcode{$species} || $species;
      $vals{species} ||= $species || $speciesin;
      $vals{species} =~ s/\d$//; 
      $vals{species}=~s/^(.)/\u$1/;
      $species = $vals{species};
      push(@sppkeys,$species) unless(grep(/$species/,@sppkeys));

      $vals{clusters} ||= "all";
      $vals{source} =~ s/_$//;
      $vals{source} =~ s/GLEANR/GleanR/;
      $vals{source} =~ s/GLEAN/SKIP/;
      $vals{source} =~ s/EISE_CGW/Wise/; # Genewise too long
      $vals{source} =~ s/NCBI_GNO/Gnomon/;
      $vals{source} =~ s/BREN_NSC/Nscan/;
      $vals{source} =~ s/DGIL_SNO/Snap/;
      
      unless ($vals{source} =~ m/SKIP/) {  
        $asource= $vals{source} unless($asource);
        my($src,$spp, $clus)= @vals{ qw(source species clusters) }; # group factors
        my @tab= @vals{ @skeys };
        $savetab{$src}{$clus}{$spp}= \@tab;

# if(OLDROUT) {      
#         # defer print til get species?
#         print "species <- \"".$vals{species}."\" # $sptag\n";
#         my $src= "\"".$vals{source}."\"";
#         my $fam= "\"".$vals{clusters}."\"";
#         print "gvals <- c( ", join(", ", map{ my $v=$vals{$_}||0; "$_=$v" } @gvals) ," );\n";
#         my @kb= map{ my $k="kb".$_; my $v=$vals{$k}||0; "$k=$v"; } @kbrange;
#         print "neard <- c(",join(",",@kb,),")\n";
#         print "dl <- rbind(dl, c(species=species,source=$src,famsize=$fam,gvals,neard))\n";
#         print "#\n";
# }        
        }
    }
    
    %vals=(); 
    @vals{ @skeys }= (0) x scalar(@skeys); #? init
    if(/\==\[ (\w+) \]/) { $vals{species}= $species= $1; }
    elsif(m,\==\>\s+(\w+),) { # ==> dpse/dpse_dgilsno.blastnear <==  head output
     $vals{species}= $species= $1;
     $athead++;
    }
    next;
    }
  chomp;
  if(s/^modelerrs=/modelerrs; /){ s/(\w+):/e$1=/g;}
  if(/\=\d/) {
    my @v= split /\s*[;,]\s*/;
    map{ my($k,$v)=split"="; $vals{$k}=$v||0; }@v;
    }
}

# if(OLDROUT) {      
# print "df <- as.data.frame(dl)\n";
# }

my $specieslist= '"'.join('", "',@sppkeys).'"';
$options .= "\nin=$input" if($input);

print_Rdata(); # print_table();
print_Rrun();
print_Rfuncs();


sub print_Rdata {
  my $cc=""; # "# "
  print "dtable <-paste(\"\n";
  print $cc,join("\t", "source", "species", "famsize", @skeys), "\n";
  foreach my $src (sort keys %savetab) {
    #print "\n";
    foreach my $clus (sort keys %{$savetab{$src}} ) { #? drop clus
    foreach my $spp (@sppkeys) {
      next unless($savetab{$src}{$clus}{$spp});
      my @v= @{$savetab{$src}{$clus}{$spp}};
      print $cc,join( "\t", $src, $spp, $clus, @v), "\n";
      }
    }
    }
  print "\")\n\n";
  print "df <- read.table(textConnection(dtable),header=T)\n";
}

sub print_Rrun {
  print <<"EOR";
#..... tandem gene predictor bar plots ..........
title <- "$title"
sublab <- "$options"
species <- c($specieslist)
asource<- "$asource"
EOR

  print <<'EOR';
cat("
#... data table ....
nd <- tandynear_xtab(df)
#... plots ....
nd.barpanels(nd, title=title)
#or# nd.barplot(nd, title=title)

#......................................
")
EOR
}

sub print_Rfuncs 
{
  print <<'EOR';
tandynear_xtab <- function(df,lpd=NULL,dfset=NULL) {
  if(is.null(dfset)) dfset <- 1 # df$famsize==famsize 
  if(is.null(lpd)) lpd <- levels(df$source)
  maxd <- ncol(df)
  i4 <- grep("kb4",colnames(df)) # need more info on $kbrange
  cmat <- c()
  for (ipd in lpd) { 
    iset <- dfset & df$source == ipd
    predname <- as.character(df$source[iset])[1]
    k5 <- matrix( as.numeric(as.matrix( df[iset, i4:(i4+4) ])),ncol=5)
    far <- matrix( as.numeric(as.matrix( df[iset, c("gfar","gscaf") ])),ncol=2)
    all <- matrix( as.numeric(as.matrix( df[iset, c("gfar","gscaf","gnear") ])),ncol=3)
    mat <- data.frame(
      k1=as.numeric(as.character(df$kb1[iset])),
      k2=as.numeric(as.character(df$kb2[iset])),
      k5=apply( k5, 1, sum),
      Far=apply( far, 1, sum)/5, # /5 or /2?
      # All=apply( all, 1, sum)/4, #?
      row.names=df$species[iset]
      )
    colnames(mat) <- paste(predname,c("1k","2k","5k","Far/5"),sep="-") #colnames(mat)
    cmat <- cbind(cmat,as.matrix(mat))
    }
  return(cmat)
  }


nd.barplot <- function(nd, col=NULL, log="", legloc="topleft",
  title="Predictor x Tandem genes", sub=NULL,
  ylab="N. Duplicate genes", xlab="Predictor-Distance"
 ) {
  if(is.null(col)) col <- heat.colors(nrow(nd)) # diff color range set / predictor
  if(is.null(sub) && !is.null(sublab)) sub <- paste("Options:",sublab)
  sf <- factor(rownames(nd), levels=species, ordered=T)
  spo <- as.character( sf[order(sf)])
  barplot(nd[spo,], beside=T, col=col, log=log,
    main = title, sub  = sub, ylab = ylab, xlab = xlab )  
  legend( legloc, cex=0.9, spo, horiz=T, col=col, lwd = c(2,2), bty="n" )
  }

nd.barpanels <- function(nd, ylim=NULL, col=NULL, title="Predictor x Tandem genes") {
  nc <- length( grep(asource,colnames(nd)) )  
  if(nc<1) nc <-1  
  np <- trunc( ncol(nd) / nc)
  if(is.null(ylim)) ylim <- c(0,max(nd))
  sf <- factor(rownames(nd), levels=species, ordered=T)
  spo <- as.character( sf[order(sf)])
  if(is.null(col)) col <- heat.colors(nrow(nd)) # diff color range set / predictor
  nl <- 2 * round((np+1.9)/2)
  nlf <- layout( matrix( 1:nl, ncol=2, byrow = T) )
  layout.show(nlf)
  
  par(mar=c(1.1,1.1,1.1,1.1))
  par(cex=1)
  plot( 1:10, 1:10, axes=F, xlab="", ylab="", type="n")
  text( 5,6, title, cex=1.3)
  subtxt <- sub(",pct",",\n pct",sublab)
  text( 5,2, paste("Options:",subtxt), cex=0.9)
  
  legend( "right", cex=0.9, spo, horiz=F, col=col, lwd = c(2,2), bty="o" )
  
  for (i in 1:np) {
   par(mar=c(4.1,4.1,2.1,2))
   par(cex=0.9)
   iset <- (i-1)*nc + c(1:nc)  # 1:4 = num sub-bar vals/predictor
   barplot( nd[spo, iset], 
     beside=T, col=col, ylim=ylim,
     ylab="N. Duplicate genes",
     xlab="Predictor-Distance"
     )  
   }
   
  }

EOR
}


#........ end

__END__

splow <- c("dmel","dsec", "dsim", "dyak", "dere", "dana", "dper", "dpse",  "dwil", "dvir", "dmoj", "dgri")


set dspp=("dmel" "dsec"  "dsim"  "dyak"  "dere"  "dana"  "dper"  "dpse"   "dwil"  "dvir"  "dmoj"  "dgri")