calcCoGAPSStat <- function (Amean, Asd, GStoGenes, numPerm=500) {
 
   # calculate Z scores
   zMatrix <- calcZ(Amean,Asd)
   
   # compute the p-value for each gene set belonging to each pattern
 
   # check input arguments
   if (!is(GStoGenes, "data.frame") && !is(GStoGenes, "list")) {
     stop("GStoGenes must be a data.frame or list with format specified in the users manual.")
   }
 	
   if (is(GStoGenes, "list")) {
     GStoGenesList <- GStoGenes
   } else {
     GStoGenesList <- list()
     for (i in 1:dim(GStoGenes)[2]) {
       GStoGenesList[[as.character(colnames(GStoGenes)[i])]] <- as.character(unique(GStoGenes[,i]))
     }
   }
   
   # get dimensions
   numGS   <- length(names(GStoGenesList))
   numPatt <- dim(zMatrix)[2]
   numG    <- dim(zMatrix)[1]+0.9999
 
   # initialize matrices
   statsUp   <- matrix(ncol = numGS, nrow = numPatt)
   statsDown <- matrix(ncol = numGS, nrow = numPatt)
   actEst    <- matrix(ncol = numGS, nrow = numPatt)
   results   <- new.env()
   zPerm     <- matrix(ncol=numPerm,nrow=numPatt)
 
   # do permutation test
   for (gs in 1:numGS) {
     genes <- GStoGenesList[[names(GStoGenesList)[gs]]]
     index <- rownames(zMatrix) %in% genes
     zValues <- zMatrix[index,1]
     numGenes <- length(zValues)
     label <- as.character(numGenes)
     
     if (!exists(label,env=results)) {
       for (p in 1:numPatt) {
         for (j in 1:numPerm) {
           temp <- floor(runif(numGenes,1,numG))
           temp2 <- zMatrix[temp,p]
           zPerm[p,j] <- mean(temp2)
         }
       }
       assign(label,zPerm,env=results)
     }
   }
 
   # get p-value
   for (p in 1:numPatt) {
 
     for (gs in 1:numGS) {
 
       genes <- GStoGenesList[[names(GStoGenesList)[gs]]]
       index <- rownames(zMatrix) %in% genes
       zValues <- zMatrix[index,p]
       zScore <- mean(zValues)
 
       numGenes <- length(zValues)
       label <- as.character(numGenes)
 
       permzValues <- get(label, env=results)[p,]
       ordering <- order(permzValues)
       permzValues <- permzValues[ordering]
 
       statistic <- sum(zScore > permzValues)
       statUpReg <- 1 - statistic/length(permzValues) 
       statsUp[p,gs] <- max(statUpReg, 1/numPerm)
 
       statistic <- sum(zScore < permzValues)
       statDownReg <- 1 - statistic/length(permzValues) 
       statsDown[p,gs] <- max(statDownReg,1/numPerm)
 
       activity <- 1 - 2*max(statUpReg, 1/numPerm)
       actEst[p,gs] <- activity
     }
 
   }
 
   # format output
   colnames(statsUp) <- names(GStoGenesList)
   colnames(statsDown) <- names(GStoGenesList)
   colnames(actEst) <- names(GStoGenesList)
 
   rownames(statsUp) <- colnames(zMatrix)
   rownames(statsDown) <- colnames(zMatrix)
   rownames(actEst) <- colnames(zMatrix)
 
   assign('GSUpreg', statsUp, env=results)
   assign('GSDownreg', statsDown, env=results)
   assign('GSActEst', actEst, env=results)
   
   return(results)
 }