PLSnet <- function(data, genes, ncom=3, qvalue=.1, df.fit=5, ...){
   ## Input:
	 # data - microarray dataset with genes as rows and samples in columns
		 # (no gene names in the first column)
	 # genes - a vector of gene names
	 # ncom  - a number of PLS components (latent variables) in PLS models
	 # qvalue- a cutoff fdr level
	 # df.fit- a number of degrees of freedom used in locfdr for fitting the mixture
	 # 		density
   ## Output:
	 # PLSnet - a matrix of interactions between genes in column 1 and column 2

	library(locfdr) # R package locfdr is required

	data <- scale(t(data)) # scale the data
	n <- ncol(data)
	s <- matrix(0, n, n)

	standr <- function(x) x/as.numeric(sqrt(t(x)%*%x))
	getInter <- function(g, genes, level){
		mm <- outer(genes,genes, paste)
		mm <- mm[upper.tri(mm)]
		mm <- mm[g$fdr<=level]
	
		if(length(mm) == 0)
			return(c(0,0))

		inter <- strsplit(mm, " ")
		int <- matrix(0, nrow=length(inter), 2)
		for(i in 1:length(inter))
			int[i,] <- inter[[i]]
		int
	}

	for(i in 1:n){ # loop through each gene
		X <- data[,-i] 
		y <- data[,i]
		
		c <- matrix(0, n-1, ncom)
		TT <- matrix(0, nrow(data), ncom)
		tX<- X
		for(k in 1:ncom){ # construct ncom PLS components
			c[,k] <- standr(t(tX)%*%y)
			TT[,k] <- tX%*%c[,k]
			tX <- tX-TT[,k]%*%solve(t(TT[,k])%*%TT[,k])%*%t(TT[,k])%*%tX
		}

		b <- solve(t(TT)%*%TT)%*%t(TT)%*%y # least squares solution
		temp <- c%*%b
		if (i == 1) s[i,] <- c(0,temp)
		else ( if (i == n) s[i,] <- c(temp,0)
		else  s[i,] <- c(temp[1:(i-1)],0,temp[i:(n-1)]) )
	}
	s <- 0.5*(s+t(s)) # construct symmetrized coefficient
	ss <- max(s)
	diag(s) <- rep(ss,length(diag(s)))
	s <- s/ss
	s <- s[upper.tri(s)]
	g1 <- locfdr(s, df=df.fit, ...) # perform locfdr testing 
	net <- getInter(g1, genes, qvalue) # extract significant interactions
	net
}