# AnalysisYeastFuncat.R
# Example of an application of the R library hcgene to the analysis of the characteristics of the FunCat taxonomy for the yeast.
# October 2007

library(hcgene);

# if VERBOSE is set to TRUE messages about the operation performed are printed on stdout
VERBOSE <- TRUE;

# if this is set to TRUE, postscript files are generated, otherwise graphics are plotted on separated windows
POSTSCRIPT <- FALSE;

if ((POSTSCRIPT) && (!file.exists("ps"))) dir.create("ps");


print(date());
#############################################################################################
######### Universal Funcat tree (not related to any species)   ##############################
#############################################################################################
if (VERBOSE) cat("**************************************** \nAnalysis of the universal FunCat tree. * \n****************************************\n\n");	
# Construction of the  "universal" Funcat tree that collects all the FunCat classes.
if (VERBOSE) cat(" Construction of the  universal Funcat trees ... \n");
gUniversalFuncat <- Do.universal.tree.Funcat();
gUniversalFuncatup <- Do.universal.tree.Funcat.up();

# plotting universal graphs
if (VERBOSE) cat(" Plotting of the  universal Funcat tree ... \n");
if (POSTSCRIPT) 
  postscript("ps/Funcat.Universal.graph.ps",onefile=FALSE) else x11();
Plot.ontology.graph(gUniversalFuncat);
if (POSTSCRIPT) 
	dev.off();
	
if (POSTSCRIPT) 
  postscript("ps/Funcat.Universal.up.graph.ps",onefile=FALSE) else x11();
	Plot.ontology.graph(gUniversalFuncatup);
if (POSTSCRIPT) 
	dev.off();

print("Statistics of the Universal Funcat tree:");
Compute.depth.statistics.GO.graph(gUniversalFuncat, ontology="00");

if (VERBOSE) cat(" Plotting the distribution of the distances from the root ... \n");
# Computing and plotting the distribution of the distances from the root
d<-Compute.distances.from.root(gUniversalFuncat, ontology="00"); 
if (POSTSCRIPT) 
  postscript("ps/Funcat.Universal.hist.depth.ps",onefile=FALSE) else x11();
Plot.hist.Funcat.depth.labels(nodes(gUniversalFuncat));
if (POSTSCRIPT) dev.off();

# statistics, histograms and ecdf of the out degree of the nodes of the universal tree
print("Statistics on the out degree of nodes of the Universal Funcat tree:");
d<-Compute.statistics.Funcat.degree(gUniversalFuncat);	
maxout <- max(d$degree$outDegree);
print(d$outDegree);

if (VERBOSE) cat(" Plotting the distribution of the out degree of the nodes of the universal tree ... \n");
if (POSTSCRIPT) {	
	postscript("ps/Funcat.universal.tree.hist.outdegree.ps",onefile=FALSE);
  hist(d$degree$outDegree, main="",xlab="out degree", breaks=0:maxout, right=FALSE);
  dev.off();
	
	postscript("ps/Funcat.universal.tree.ecdf.outdegree.ps",onefile=FALSE);
	y<-ecdf(d$degree$outDegree);
	plot(y, verticals=TRUE, ylab="", xlab="out degree", main="");
	dev.off();
} else {    
  x11(); hist(d$degree$outDegree, main="",xlab="out degree", breaks=0:maxout, right=FALSE);
	y<-ecdf(d$degree$outDegree);
	x11(); plot(y, verticals=TRUE, ylab="", xlab="out degree", main="");
}	

if (VERBOSE) cat(" Building and plotting subtrees depending on the level ... \n");	
# building and plotting subtrees depending on the level
# tree limited to the nodes with maximum depth equal to 1
nodes<-Select.Funcat.classes.by.depth(gUniversalFuncat,1);
sg <- subGraph(nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.level1.ps",onefile=FALSE)  else x11();
Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=0.9,width=1.2,color="black", fontcolor="black");
if (POSTSCRIPT)
  dev.off();

# tree limited to the nodes with maximum depth equal to 2
nodes<-Select.Funcat.classes.by.depth(gUniversalFuncat,2);
sg <- subGraph(nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.level2.ps",onefile=FALSE) else x11();
#Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=0.9,width=1.2,color="black", fontcolor="black");
Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT) dev.off();

# tree limited to the nodes with maximum depth equal to 2
nodes<-Select.Funcat.classes.by.depth(gUniversalFuncat,3);
sg <- subGraph(nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.level3.ps",onefile=FALSE) else x11();
Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT) dev.off();

if (VERBOSE) cat(" Building and plotting subtrees rooted at specific nodes ... \n");	
# building and plotting specific subtrees of the Funcat hierarchy:
# O1 subtree: METABOLISM
s01.nodes <- Subtree.nodes(gUniversalFuncat, "01");
g01 <- subGraph(s01.nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.01.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g01,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT) dev.off();

# O1.01 subtree: amino acid metabolism
s.nodes <- Subtree.nodes(gUniversalFuncat, "01.01");
g <- subGraph(s.nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.01.01.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT) dev.off();

# 01.02.02 nitrogen metabolism
s.nodes <- Subtree.nodes(gUniversalFuncat, "01.02.02");
g <- subGraph(s.nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.01.02.02.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g,fontsize=12,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
if (POSTSCRIPT) dev.off();

# 30 subtree: CELLULAR COMMUNICATION/SIGNAL TRANSDUCTION MECHANISM
s.nodes <- Subtree.nodes(gUniversalFuncat, "30");
g <- subGraph(s.nodes, gUniversalFuncat);
if (POSTSCRIPT)
  postscript("ps/Funcat.Universal.graph.30.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g,fontsize=8,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="transparent");
if (POSTSCRIPT) dev.off();


#############################################################################################
######### Funcat analysis for the Yeast ####################################################
#############################################################################################
if (VERBOSE) {
  cat("****************************************\n");	
  cat("Analysis of FunCat trees for the yest. *\n");	
  cat("****************************************\n\n");	
}

if (VERBOSE) cat(" Construction of the list of annotated genes in the yeast ... \n"); 
# Construction of the list of annotated genes in the yeast
Yeast.Funcat.specific <- Get.yeast.Funcat.specific.classes();
Yeast.Funcat.general <- Get.yeast.Funcat.all.classes(Yeast.Funcat.specific);
# writing the list of yeast annotated genes to a file
Write.gene.classes.associations(Yeast.Funcat.general,filename="yeast.gene.Funcat.general.txt");

if (VERBOSE) cat(" Construction of the table of multilabels associated with each gene in the yeast ... \n"); 
# Construction of the table of multilabels associated with each gene in the yeast
Yeast.Funcat.Table <- Build.Funcat.Table.labels(Yeast.Funcat.general, store=TRUE, file="Yeast.Funcat.Table.data");
# writing the table to a file
write.table(Yeast.Funcat.Table, file="Yeast.gene.Funcat.general.table.txt");


# basic statistics
print("Statistics of Funcat gene labels for all annotated genes:");
l<-Compute.statistics.gene.labels(Yeast.Funcat.general);
print(l);

print("Statistics of Funcat classes for all the annotated genes:");
l<-Compute.statistics.functional.classes(Yeast.Funcat.Table[,2:ncol(Yeast.Funcat.Table)]); # the first variable refers to the dummy "00" class
print(l);

if (VERBOSE) cat(" Plotting the distribution of the cardinality of gene labels (number of classes for each gene) ... \n"); 
# plotting the histogram of the cardinality of gene labels (number of classes for each gene)
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.general.labels.histogram.ps",onefile=FALSE) else x11();
Do.table.hist.cardinality.labels (Yeast.Funcat.general, title="");
if (POSTSCRIPT)
  dev.off();

# plotting the ecdf of the cardinality of gene labels (number of classes for each gene)	
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.general.labels.ecdf.ps",onefile=FALSE) else x11();
Do.ecdf.cardinality.labels (Yeast.Funcat.general, title="");
if (POSTSCRIPT)
  dev.off();
	
	
l <- Count.examples.per.class(Yeast.Funcat.Table);

if (VERBOSE) cat(" Plotting the distribution of the cardinality of Funcat classes ... \n"); 
# Histogram of the cardinality of Funcat classes
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.general.per.class.hist.ps") else x11();
Plot.histogram.gene.per.class(Yeast.Funcat.Table);	
if (POSTSCRIPT)
  dev.off();

# Ecdf of the cardinality of Funcat classes
if (POSTSCRIPT) 
   postscript("ps/Yeast.Funcat.general.per.class.distr.ps") else x11();
Plot.distribution.gene.per.class(l[[1]], card=1:5813)
if (POSTSCRIPT)
  dev.off();

if (VERBOSE) cat(" Construction of the trees of the classes with annotated genes for the yeast ... \n"); 
# Construction of the graph of the classes with annotated genes for the yeast
Yeast.Funcat.classes <- names(Yeast.Funcat.Table);
cat("Num. Funcat classes for the yeast: ", length(Yeast.Funcat.classes), "\n");

gYeastFuncat <- subGraph(Yeast.Funcat.classes, gUniversalFuncat);
gYeastFuncatup <- subGraph(Yeast.Funcat.classes, gUniversalFuncatup);

if (VERBOSE) cat(" Plotting Yeast Funcat graphs ... \n"); 
# plotting Yeast Funcat graphs
if (POSTSCRIPT) 
  postscript("ps/Yeast.Funcat.graph.ps",onefile=FALSE) else x11();
	Plot.ontology.graph(gYeastFuncat);
if (POSTSCRIPT) 
	dev.off();
	
if (POSTSCRIPT) 
  postscript("ps/Yeast.Funcat.up.graph.ps",onefile=FALSE) else x11();
	Plot.ontology.graph(gYeastFuncatup);
if (POSTSCRIPT) 
	dev.off();

if (VERBOSE) cat(" Computation of some statistics on the depth of the Yeast Funcat graph ... \n"); 
# Computation of some statistics on the "depth" of the Yeast Funcat graph
print("Statistics of the Yeast Funcat tree:");
Compute.depth.statistics.GO.graph(gYeastFuncat, ontology="00");

if (VERBOSE) cat(" Computing and plotting the distribution of the distances from the root ... \n"); 
# Computing and plotting distribution of the distances from the root
d<-Compute.distances.from.root(gYeastFuncat, ontology="00"); 
if (POSTSCRIPT) 
  postscript("ps/Yeast.Funcat.hist.depth.ps",onefile=FALSE) else x11();
hist(d, main="",xlab="distances from the root");
if (POSTSCRIPT) 
  dev.off();

# statistics, histograms and ecdf of the out degree of the nodes of the Yeast Funcat graph
print("Statistics on the out degree of nodes of the Universal Funcat tree:");
d<-Compute.statistics.Funcat.degree(gYeastFuncat);	
maxout <- max(d$degree$outDegree);
print(d$outDegree);

if (VERBOSE) cat(" Plotting the distribution of the out degree of the nodes of the Yeast Funcat graph ... \n"); 
if (POSTSCRIPT) {	
	postscript("ps/Yeast.Funcat.hist.outdegree.ps",onefile=FALSE);
  hist(d$degree$outDegree, main="",xlab="out degree", breaks=0:maxout, right=FALSE);
  dev.off();
	
	postscript("ps/Yeast.Funcat.ecdf.outdegree.ps",onefile=FALSE);
	y<-ecdf(d$degree$outDegree);
	plot(y, verticals=TRUE, ylab="", xlab="out degree", main="");
	dev.off();
} else {    
  x11(); hist(d$degree$outDegree, main="",xlab="out degree", breaks=0:maxout, right=FALSE);
	y<-ecdf(d$degree$outDegree);
	x11(); plot(y, verticals=TRUE, ylab="", xlab="out degree", main="");
}	

if (VERBOSE) cat(" Construction and plotting of subtrees of Yeast Funcat depending on the depth of nodes ... \n"); 
# building and plotting subtrees of Yeast Funcat depending on the level
# tree limited to the nodes with depth 1
nodes<-Select.Funcat.classes.by.depth(gYeastFuncat,1);
sg <- subGraph(nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.level1.ps",onefile=FALSE) else x11();
Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=0.9,width=1.2,color="black", fontcolor="black");
if (POSTSCRIPT)
  dev.off();

# tree limited to the nodes with depth 2
nodes<-Select.Funcat.classes.by.depth(gYeastFuncat,2);
sg <- subGraph(nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.level2.ps",onefile=FALSE) else x11();
Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT)
  dev.off();

# tree limited to the nodes with depth 2
nodes<-Select.Funcat.classes.by.depth(gYeastFuncat,3);
sg <- subGraph(nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.level3.ps",onefile=FALSE) else x11();
Pretty.plot.graph(sg,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT)
  dev.off();

if (VERBOSE) cat(" Construction and plotting of subtrees of Yeast Funcat rooted at specific nodes ... \n"); 
# building and plotting specific subtrees of the Yeast Funcat hierarchy:
# O1 subtree: METABOLISM
s01.nodes <- Subtree.nodes(gYeastFuncat, "01");
g01 <- subGraph(s01.nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.01.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g01,fontsize=12,fillcolor="lightgreen",height=1.2,width=3.2,color="black", fontcolor="transparent");
if (POSTSCRIPT)
  dev.off();

# O1.01 subtree: amino acid metabolism
s.nodes <- Subtree.nodes(gYeastFuncat, "01.01");
g <- subGraph(s.nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.01.01.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g,fontsize=12,fillcolor="lightgreen",height=1.2,width=4.2,color="black", fontcolor="transparent");
if (POSTSCRIPT)
  dev.off();
	
# 01.02.02 nitrogen metabolism
s.nodes <- Subtree.nodes(gYeastFuncat, "01.02.02");
g <- subGraph(s.nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.01.02.02.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g,fontsize=12,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
if (POSTSCRIPT)
  dev.off();

# 30 subtree: CELLULAR COMMUNICATION/SIGNAL TRANSDUCTION MECHANISM
s.nodes <- Subtree.nodes(gYeastFuncat, "30");
g <- subGraph(s.nodes, gYeastFuncat);
if (POSTSCRIPT)
  postscript("ps/Yeast.Funcat.graph.30.ps",onefile=FALSE) else x11();
Pretty.plot.graph(g,fontsize=8,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
if (POSTSCRIPT)
  dev.off();

if (VERBOSE) cat(" Construction and plotting subtrees with classes selected according to the  number of examples ... \n"); 
### Building and plotting subtrees with classes selected according to the  number of examples
# 1. Selection of classes with more than 20 genes
nodes <- Select.functional.classes.by.cardinality(Yeast.Funcat.Table, 20); 
# 2. Building  the corresponding graph
gYeast.card.20 <- subGraph(nodes, gYeastFuncat);

# 1. Selection of classes with more than 50 genes
nodes <- Select.functional.classes.by.cardinality(Yeast.Funcat.Table, 50); 
# 2. Building  the corresponding graph
gYeast.card.50 <- subGraph(nodes, gYeastFuncat);

# 1. Selection of classes with more than 100 genes
nodes <- Select.functional.classes.by.cardinality(Yeast.Funcat.Table, 100); 
# 2. Building  the corresponding graph
gYeast.card.100 <- subGraph(nodes, gYeastFuncat);

# 1. Selection of classes with more than 200 genes
nodes <- Select.functional.classes.by.cardinality(Yeast.Funcat.Table, 200); 
# 2. Building  the corresponding graph
gYeast.card.200 <- subGraph(nodes, gYeastFuncat);

# 1. Selection of classes with more than 400 genes
nodes <- Select.functional.classes.by.cardinality(Yeast.Funcat.Table, 400); 
# 2. Building  the corresponding graph
gYeast.card.400 <- subGraph(nodes, gYeastFuncat);

if (POSTSCRIPT) {
	postscript("ps/Yeast.Funcat.graph.card20.ps",onefile=FALSE);
	Plot.ontology.graph(gYeast.card.20, fcolor="white");
	dev.off();
	postscript("ps/Yeast.Funcat.graph.card50.ps",onefile=FALSE);
	Plot.ontology.graph(gYeast.card.50, fcolor="white");
	dev.off();
	postscript("ps/Yeast.Funcat.graph.card100.ps",onefile=FALSE);
	Plot.ontology.graph(gYeast.card.100, fcolor="white");
	dev.off();
	postscript("ps/Yeast.Funcat.graph.card200.ps",onefile=FALSE);
	Pretty.plot.graph(gYeast.card.200,fontsize=12,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
	dev.off();
	postscript("ps/Yeast.Funcat.graph.card400.ps",onefile=FALSE);
	Pretty.plot.graph(gYeast.card.400,fontsize=14,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
	dev.off();
} else {  
  x11(); Plot.ontology.graph(gYeast.card.20, fcolor="white");
  x11(); Plot.ontology.graph(gYeast.card.50, fcolor="white");
  x11(); Plot.ontology.graph(gYeast.card.100, fcolor="white");
	x11(); Pretty.plot.graph(gYeast.card.200,fontsize=12,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
	#Plot.ontology.graph(gYeast.card.200, fcolor="white");
	x11(); Pretty.plot.graph(gYeast.card.400,fontsize=14,fillcolor="lightgreen",height=1.2,width=5.2,color="black", fontcolor="black");
}

print(date());