#########################
# Under R 2.5.1 version #
#########################







####################################
# loading of the required packages #
####################################

library(pls)
library(GeneNet)
library(graph)
library(Rgraphviz)

####################################
# Run this commands for example    #
####################################
example = function() 
{
   #####################################
   # Cut and paste the following lines #
   #####################################
   data(ecoli)
   # Apply PLS-PC network using 3 components on ecoli data
   # Plot the resulting graph
   result = PLSPC.network(ecoli)
   #####################
   # End cut and paste #
   #####################
}

PLSPC.network = function(X, nbre = 3, method = "kernelpls", viz = 1){

#########
# INPUT #
#########

# X : Microarray experiments (n * p)
#         n = number of experiments
#         p = number of genes
# 
# nbre = Number of PLS components to include in the PLSPC network (default = 3 components)
#
# method = implementation of the PLS algorithm (default = Kernel PLS)
#
# viz = visualization (default = 1 -> visualization desired)  

##############################
# Output : list of variables #
##############################
#
#      pcor  : PLSPC matrix 
#      graph : data.frame with the significant edges
#      gr    : a graph object
#      coef  : PLS matrix coefficients
#
             
if (viz ==1){
X = X[ , as.logical(1-duplicated(colnames(X)))]
X = X[ , which(colnames(X) != "NA")]
}

if ((nbre > qr(scale(X, center = TRUE, scale = TRUE))$rank) && (is.null(nbre) == FALSE)){ 
print(paste("Number of PLS components must be < ", qr(scale(X, center = TRUE, scale = TRUE))$rank))
}

else{
##################################################                                     
# Construction of the Partial Correlation matrix #
#            Package PLS required                #
##################################################
pcor = coef = diag(ncol(X))
for (i in 1:ncol(X)){
	temp = mvr(X[, i]~X[, -i], ncomp = nbre, method = method, scale = TRUE, stripped = TRUE)
  coef[i, -i] = as.data.frame(temp$coefficients)[ , nbre]
}

for (i in 1:ncol(X)){
	 for (j in i:ncol(X)){
	   	if (i==j) pcor[i, j] = 1
		  else {
        if(sign(coef[j, i]) != sign(coef[i, j])) pcor[i, j] = 0
        else{
			   pcor[i, j] = sign(coef[j, i])*sqrt(coef[i, j]*coef[j, i])	
        }
      }
    pcor[j, i] = pcor[i, j]
	 }
  }
}

if (viz == 1){
#########################################
# Local fdr multiple testing on the PLS #
#########################################
test.results = ggm.test.edges(pcor, plot = FALSE)
signif = test.results$prob > 0.8                                                                 

if (sum(signif)>1){

####################################
# Visualization of the PLSPC network #
####################################

node.labels = colnames(X)
gr = ggm.make.graph(test.results[signif, ], node.labels, drop.singles = FALSE)
x11(); plot(gr,"twopi")
return(list(pcor = pcor, graph = test.results[signif, ], gr = gr, coef = coef))}
else{print("No significant interaction !");return(list(pcor = pcor, coef = coef))}
}
else{return(list(pcor = pcor, coef = coef))}
}