# goblstats9g2.R
# 30 dec 06 ... for gene-detail GO associations (gene6-g2)
# ---- 29 Dec 06: *** Gene-Detail GO Associations (g2) give best phylo separation ***

# setwd("/Users/gilbertd/active/gmod04work/eugenomes/eugenomestats/egstats9")
# setwd("/Users/gilbertd/Desktop/dspp-work/eugenomestats/egstats9")

version <- "9g3"  
gvar <- "genes6"; gvlab <- "Gene matches (HSP groups) "
GOcat <- "GOdetail" ## FIXME ; g2
#datafn <- "go2-spp-"
datafn <- "go3-spp-"

source("../goassocplot9.R")
# uses library("gmodels") # CrossTable

chilim <- 1.5 # default ** this is best for needs ?
# chilim <- 2 # should use this for cele-dpulx-dros compares to reduce choice overload

gocats <- read.gocat(file="../gocatall.tab")
golabs <- list(P="GO_Process",F="GO_Function",C="GO_CellLocation")

protdbs <- c( "modMM", "modDM", "modCE", "modSC")
protmat <- as.matrix( c( "Mouse", "Fruitfly", "Worm", "Yeast")) ; row.names(protmat) <- protdbs
proteomes <- protmat[protdbs,]

ceorder   <- c("cele", "dpul", "dmel","dpse","dvir");
flyorder  <- c("dmel","dsim","dsec","dyak","dere","dana","dper","dpse","dwil","dmoj","dvir","dgri")
fly1order <- c("dmel","dsim","dyak","dere","dana","dpse","dwil","dmoj","dvir","dgri")
allorder  <- c("cele", "dpul", flyorder)

modorder <- c( "modDM", "modCE", "modMM",)

# 1letter abbrevs?
sppabbrev <- as.matrix( c( "C", "W", "m", "s", "c", "y", "e", "a", "r", "p", "w", "j", "v", "g"))
row.names(sppabbrev) <- allorder
sppabbrev2 <- as.matrix( c( "Ce", "Wf", "Dm", "Ds", "Dc", "Dy", "De", "Da", "Dr", "Dp", "Dw", "Dj", "Dv", "Dg"))
row.names(sppabbrev2) <- allorder
sppnames <- as.matrix( c( "C.elegans", "Daphnia_pulex", "Dro.mel","Dro.sim","Dro.sec","Dro.yak","Dro.ere","Dro.ana","Dro.per","Dro.pse","Dro.wil","Dro.moj","Dro.vir","Dro.gri"))
row.names(sppnames) <- allorder

mods2 <-  modorder # c("modDM") # modorder[-1]
# for (mod in modorder) {
for (mod in  mods2) {

cat("#------ Working on ",mod,"---------------------\n")

gospp <- read.delim(fn<-paste(datafn,mod,".tab",sep=""))
cat("# read ",fn, " dim=",dim(gospp),"\n")

iorder <- 1
for (iorder in 1:2) { ## 1:2 

if(iorder==1) { spporder1 <- flyorder; chilim <- 1.5 }
else { spporder1 <- ceorder; chilim <- 2.0 }
spp1 <- spporder1[1]
nspp <- length(spporder1)

dxtabs <- gospp[,spporder1]
row.names(dxtabs) <- gospp[,1]  # GO ids now

proteomes <- protmat[mod,]

gdir <- paste(gvar,"-",spp1,nspp,"-",tolower(proteomes),"-v",version,"/",sep="")
linkpage <- paste("index-detail.html",sep="")
crossfile <- paste(gdir,"index-summary",sep="")

cat("# output:",gdir,"\n")
dir.create(gdir)

plotheader <- paste(
"<h1>Gene Function Associations in Species Genomes</h1>\n",
##"<h1>Gene matches x Species associations for GO Categories</h1>\n",
" <h2>Deviations in gene matches for GO categories by species genomes. </h2>\n",
" These deviations <i>suggest</i> where species may differ in biological function groups.<br>\n",
" Target species genomes: ",paste(sppnames[spporder1,],collapse=", "),"<br>\n",
" Source proteome: ",paste(proteomes,collapse=", "),
" [version ",version,"]","<br>\n",
" Click on term for detailed information, or see <a href='",linkpage,"'>here</a>. <br>\n",
sep="")

png()
goassocplot(dxtabs, chilimit=chilim, crossfile=crossfile, htmlheader=plotheader )

# chilimit=1.5, << good or stronger?

}  # iorder

}  # for mod