#!/usr/bin/perl
# R stats convert

@gvals= qw(npaired gscaf gfar gnear);
print "dl <- c() \n";
while(<DATA>){
  s/^\s+//;
  unless(/^\w/) {
    if(%vals) { 
      print "species <- \"".$vals{source}."\"\n";
      print "famsize <- \"".$vals{clusters}."\"\n";
      print "gvals <- c( ", join(", ", map{ "$_=$vals{$_}" } @gvals) ," );\n";
      @kb= map{ $k="kb".$_; $v=$vals{$k}||0; "$k=$v"; } (1..20);
      print "neard <- c(",join(",",@kb,),")\n";
      print "dl <- rbind(dl, c(species=species,famsize=famsize,gvals,neard))\n";
      print "#\n";
      #last if($n++>3);
    }
    %vals=(); 
    next;
    }
  chomp;
  @v= split /\s*[;,]\s*/;
  map{ ($k,$v)=split"="; $vals{$k}=$v; }@v;
}
print "df <- as.data.frame(dl)\n";

=item protnear table

cat dmoj1_gleanr*.genes dmoj1_gleanr*.idchains | \
perl $td/protnear.perl -near 20000 -clust 10 -clust 80

##!/bin/tcsh
## protnear.sh

set td=/bio/bio-grid/dpulex/prots/

foreach da (cele1_modCE_aa dpulex1_gnomon_aa dgri1_gnomon_aa dmel4_modDM_aa mmus1_modMM_aa)
  echo "# ==============[ $da ]=========================="
  echo "# protein near steps cluster size 2-6"
  gzcat $da.genes.gz $da.blastp150.idchains.gz | perl $td/protnear.perl -clus 1 -clust 5 -near 20000
  echo "#"
  echo "# protein near steps cluster size 10-80"
  cat $da.genes.gz $da.blastp150.idchains.gz | perl $td/protnear.perl -clus 10 -clust 80 -near 20000
  echo
end

foreach da (dgri1_genewise_aa dgri1_gleanr_aa dgri1_gnomon_aa)
  echo "# ==============[ $da ]=========================="
  echo "# protein near steps cluster size 2-6"
  gzcat $da.genes.gz $da.blastp150.idchains.gz | perl $td/protnear.perl -clus 1 -clust 5 -near 20000
  echo "#"
  echo "# protein near steps cluster size 10-80"
  gzcat $da.genes.gz $da.blastp150.idchains.gz | perl $td/protnear.perl -clus 10 -clust 80 -near 20000
  echo
end

=item plots

famsize <- "10..80"
famsize <- "1..5"
famsize <- "all"

source <- "GleanR"  
source <- "Genewise"  # 
source <- "Gnomon"
species <- "Dros-7" #"Dsec"


dfset <- df$famsize==famsize & df$source==source
nd <- df[dfset,8:27]
rownames(nd) <- paste( df[dfset,1],df[dfset,2],sep="-")
nx <- 1000 * as.integer( gsub("kb","",colnames(nd[1,])) )

nd <- df[df$famsize==famsize,7:26]
rownames(nd) <- df[df$famsize==famsize,1]
nx <- 1000 * as.integer( gsub("kb","",colnames(nd[1,])) )

# damn r
# no way!! for (j in 1:5) { for (i in 1:20) nd[j,i] <- as.numeric(as.character(nd[j,i])) }
nd2 <- matrix(as.numeric(as.matrix(nd)),nrow=nrow(nd))
rownames(nd2) <- rownames(nd)
colnames(nd2) <- colnames(nd)
nd <- nd2
# rownames(nd) <- c("C.elegans","Daphnia","Dros.gri","Dros.mel","Mouse")

ymax <- 500 # 1600

plot( nx, nd[2,],  type="n", 
 main=paste("Tandem distances, family size",famsize),
 xlab="Tandem Distance (bp)",  log="x", ylim=c(0,ymax),
 ylab="N. Duplicate genes")

nr <- nrow(nd) 
for (j in 1:nr) 
 lines(nx, nd[j,], type="b",pch=j+1,col=j+1)

nr1 <- nr+1 
legend("topright", rownames(nd), col=2:nr1, pch=2:nr1, lwd = c(1,1), bty="n")

#..... proportion of total dupl. found

td <- df[df$famsize==famsize,4:6]
td2 <- matrix(as.numeric(as.matrix(td)), nrow=nrow(td))
tn <- apply(td2,1,sum)

ymax <- 0.10
plot( nx, nd[2,],  type="n", 
 main=paste("Tandem distances, family size",famsize),
 xlab="Tandem Distance (bp)",  log="x", ylim=c(0,ymax),
 ylab="% of Duplicate genes")

nr <- nrow(nd) 
for (j in 1:nr) 
 lines(nx, nd[j,]/tn[j], type="b",pch=j+1,col=j+1)

nr1 <- nr+1 
legend("topright", rownames(nd), col=2:nr1, pch=2:nr1, lwd = c(1,1), bty="n")

#..... predictor bar plot, 1000 dist vs all gnear/gfar/gdup ??

# x steps : species
# x-bars  : predictors
#   ..      -- at 1k, at 5k..20k, at > 20k

lsp <- levels(df$species)
lpd <- levels(df$source)

famsize <- "1..5"
famsize <- "10..80"
famsize <- "all"

dfset <- df$famsize==famsize 

nd <- c()
for (ipd in lpd) { 
dfseti <- dfset & df$source == ipd
ndi <-  matrix(as.numeric(as.matrix(df[dfseti, 8])),ncol=1) 
ndk <-  matrix(as.numeric(as.matrix(df[dfseti, 13:27])),ncol=15) 
ndk <-  matrix( apply(ndk, 1, sum),ncol=1) 

ndf <-  matrix(as.numeric(as.matrix(df[dfseti, 5:6])),ncol=2) # gscaf,gfar
ndf <-  matrix( apply(ndf, 1, sum),ncol=1) 
ndf <-  ndf/10 # scale down to match others for plot (???)

rownames(ndi) <- df[dfseti,1]
colnames(ndi) <- paste( (df[dfseti,2])[1], "1k", sep="-")
colnames(ndk) <- paste( (df[dfseti,2])[1], "10k", sep="-")
colnames(ndf) <- paste( (df[dfseti,2])[1], "Far", sep="-")
nd <- cbind( nd, ndi, ndk, ndf )
}

# by predictor main X group: YES, see spp trends for like-dmel effects
cols <- heat.colors(nrow(nd)) # diff color range set / predictor

barplot(nd,beside=T,col=cols,
  main = paste("Dmel-effect Predictor x Tandem genes, clusters",famsize),
  ylab = "N. Duplicate genes",
  xlab = "Predictor-Distance"
  )
legend("topleft", rownames(nd), horiz=T, col=cols, lwd = c(2,2), bty="n")


#...........


#.......  
plot( nx, nd[2,],  type="n", 
  main="Tandem duplicates",
  xlab="Distance",  
  ylab="N. Genes")
for (j in 1:5) 
  lines(nx, nd[j,], type="b",pch=j+1,col=j+1)
  
legend("topright", rownames(nd), col=2:6, pch=2:6, lwd = c(1,1), bty="n")

=cut


__DATA__

# ==============[ dgri1_genewise_aa ]==========================
# protein near steps cluster size 2-6
source=Dgri_GeneWise; ngenes=19258; npaired=2212; ngroup=996; noloc=0; skip=322; clusters=1..5; neardist=20000
 gscaf=1779; gfar=542; gnear=461; gsame=28
 kb1=139, kb2=115, kb3=56, kb4=55, kb5=44, kb6=28, kb7=34, kb8=24, kb9=16, kb10=13, kb11=6, kb12=2, kb13=10, kb14=6, kb15=12, kb16=6, kb17=12, kb18=12, kb19=2, kb999=542,
 rv1=70, rv2=61, rv3=16, rv4=18, rv5=18, rv6=8, rv7=16, rv8=10, rv9=8, rv11=4, rv13=4, rv14=2, rv16=2, rv17=6, rv18=6, rv19=2, rv999=223,
 revp1=0.504, revp2=0.530, revp3=0.286, revp4=0.327, revp5=0.409, revp6=0.286, revp7=0.471, revp8=0.417, revp9=0.500, revp11=0.667, revp13=0.400, revp14=0.333, revp16=0.333, revp17=0.500, revp18=0.500, revp19=1.000, revp999=0.411,
#
# protein near steps cluster size 10-80
source=Dgri_GeneWise; ngenes=19258; npaired=662; ngroup=147; noloc=0; skip=1171; clusters=10..80; neardist=20000
 gscaf=662; gfar=487; gnear=160; gsame=12
 kb1=55, kb2=46, kb3=30, kb4=19, kb5=21, kb6=24, kb7=16, kb8=16, kb9=12, kb10=14, kb11=23, kb12=6, kb13=8, kb14=2, kb15=9, kb16=12, kb18=8, kb19=4, kb20=2, kb999=487,
 rv1=33, rv2=20, rv3=16, rv4=14, rv5=9, rv6=6, rv7=8, rv8=9, rv9=2, rv10=2, rv11=8, rv12=2, rv13=2, rv14=2, rv15=2, rv16=2, rv19=2, rv20=2, rv999=238,
 revp1=0.600, revp2=0.435, revp3=0.533, revp4=0.737, revp5=0.429, revp6=0.250, revp7=0.500, revp8=0.562, revp9=0.167, revp10=0.143, revp11=0.348, revp12=0.333, revp13=0.250, revp14=1.000, revp15=0.222, revp16=0.167, revp19=0.500, revp20=1.000, revp999=0.489,

# ==============[ dgri1_gleanr_aa ]==========================
# protein near steps cluster size 2-6
source=Dgri_GleanR; ngenes=16901; npaired=4208; ngroup=1833; noloc=0; skip=881; clusters=1..5; neardist=20000
 gscaf=3556; gfar=689; gnear=912; gsame=48
 kb1=388, kb2=228, kb3=133, kb4=90, kb5=72, kb6=50, kb7=40, kb8=34, kb9=14, kb10=16, kb11=18, kb12=18, kb13=10, kb14=18, kb15=4, kb16=10, kb17=6, kb18=10, kb19=13, kb20=8, kb999=689,
 rv1=219, rv2=121, rv3=83, rv4=48, rv5=48, rv6=30, rv7=30, rv8=16, rv9=10, rv10=10, rv11=10, rv12=14, rv13=8, rv14=14, rv15=2, rv16=4, rv17=2, rv18=2, rv19=2, rv20=2, rv999=486,
 revp1=0.564, revp2=0.531, revp3=0.624, revp4=0.533, revp5=0.667, revp6=0.600, revp7=0.750, revp8=0.471, revp9=0.714, revp10=0.625, revp11=0.556, revp12=0.778, revp13=0.800, revp14=0.778, revp15=0.500, revp16=0.400, revp17=0.333, revp18=0.200, revp19=0.154, revp20=0.250, revp999=0.705,
#
# protein near steps cluster size 10-80
source=Dgri_GleanR; ngenes=16901; npaired=2768; ngroup=502; noloc=0; skip=2212; clusters=10..80; neardist=20000
 gscaf=2757; gfar=1354; gnear=947; gsame=40
 kb1=374, kb2=250, kb3=224, kb4=174, kb5=202, kb6=117, kb7=147, kb8=97, kb9=104, kb10=75, kb11=110, kb12=91, kb13=56, kb14=69, kb15=69, kb16=77, kb17=74, kb18=62, kb19=44, kb20=44, kb999=1354,
 rv1=260, rv2=170, rv3=134, rv4=122, rv5=114, rv6=85, rv7=103, rv8=64, rv9=74, rv10=51, rv11=68, rv12=64, rv13=42, rv14=46, rv15=47, rv16=47, rv17=42, rv18=46, rv19=26, rv20=33, rv999=1069,
 revp1=0.695, revp2=0.680, revp3=0.598, revp4=0.701, revp5=0.564, revp6=0.726, revp7=0.701, revp8=0.660, revp9=0.712, revp10=0.680, revp11=0.618, revp12=0.703, revp13=0.750, revp14=0.667, revp15=0.681, revp16=0.610, revp17=0.568, revp18=0.742, revp19=0.591, revp20=0.750, revp999=0.790,

# ==============[ dgri1_gnomon_aa ]==========================
# protein near steps cluster size 2-6
source=Dgri_Gnomon; ngenes=17922; npaired=3868; ngroup=1674; noloc=0; skip=738; clusters=1..5; neardist=20000
 gscaf=3051; gfar=645; gnear=917; gsame=210
 kb1=362, kb2=210, kb3=167, kb4=104, kb5=93, kb6=54, kb7=48, kb8=36, kb9=24, kb10=10, kb11=21, kb12=37, kb13=14, kb14=14, kb15=8, kb16=10, kb17=15, kb18=16, kb19=16, kb20=8, kb999=645,
 rv1=210, rv2=124, rv3=97, rv4=54, rv5=52, rv6=32, rv7=32, rv8=21, rv9=8, rv10=6, rv11=6, rv12=16, rv13=10, rv14=8, rv15=2, rv17=2, rv18=10, rv19=4, rv20=4, rv999=449,
 revp1=0.580, revp2=0.590, revp3=0.581, revp4=0.519, revp5=0.559, revp6=0.593, revp7=0.667, revp8=0.583, revp9=0.333, revp10=0.600, revp11=0.286, revp12=0.432, revp13=0.714, revp14=0.571, revp15=0.250, revp17=0.133, revp18=0.625, revp19=0.250, revp20=0.500, revp999=0.696,
#
# protein near steps cluster size 10-80
source=Dgri_Gnomon; ngenes=17922; npaired=2240; ngroup=494; noloc=0; skip=1918; clusters=10..80; neardist=20000
 gscaf=2227; gfar=1335; gnear=924; gsame=64
 kb1=414, kb2=238, kb3=212, kb4=155, kb5=184, kb6=111, kb7=147, kb8=93, kb9=111, kb10=92, kb11=112, kb12=98, kb13=48, kb14=68, kb15=71, kb16=72, kb17=65, kb18=53, kb19=50, kb20=51, kb999=1335,
 rv1=265, rv2=147, rv3=126, rv4=100, rv5=102, rv6=68, rv7=96, rv8=55, rv9=73, rv10=56, rv11=63, rv12=63, rv13=34, rv14=40, rv15=47, rv16=44, rv17=35, rv18=32, rv19=28, rv20=32, rv999=1048,
 revp1=0.640, revp2=0.618, revp3=0.594, revp4=0.645, revp5=0.554, revp6=0.613, revp7=0.653, revp8=0.591, revp9=0.658, revp10=0.609, revp11=0.562, revp12=0.643, revp13=0.708, revp14=0.588, revp15=0.662, revp16=0.611, revp17=0.538, revp18=0.604, revp19=0.560, revp20=0.627, revp999=0.785,

#