@@ -331,20 +331,60 @@ function(object, threshold, lp, axis)

     ## keep only a subset of markers for each pattern depending on the type of threshold

     if (threshold == "cut") # all markers which achieve minimum rank

     {

-        rankCutoff <- sapply(1:ncol(markerRanks), function(patternIndex)

-        {

-            patternRank <- markerRanks[,patternIndex]

-            return(max(patternRank[patternRank == apply(markerRanks, 1, min)]))

-        })

-        markersByPattern <- lapply(1:ncol(markerRanks), function(patternIndex)

-            rownames(markerRanks)[markerRanks[,patternIndex] <= rankCutoff[patternIndex]])

-        names(markersByPattern) <- colnames(markerRanks)

+      simplicityGENES <- function(As, Ps) {

+        # rescale p's to have max 1

+        pscale <- apply(Ps,1,max)

+        

+        # rescale A in accordance with p's having max 1

+        As <- sweep(As, 2, pscale, FUN="*")

+        

+        # find the A with the highest magnitude

+        Arowmax <- t(apply(As, 1, function(x) x/max(x)))

+        pmax <- apply(As, 1, max)

+        

+        # determine which genes are most associated with each pattern

+        ssl <- matrix(NA, nrow=nrow(As), ncol=ncol(As),

+                      dimnames=dimnames(As))

+        for (i in 1:ncol(As)) {

+          lp <- rep(0, ncol(As))

+          lp[i] <- 1

+          ssl.stat <- apply(Arowmax, 1, function(x) sqrt(t(x-lp)%*%(x-lp)))

+          ssl[order(ssl.stat),i] <- 1:length(ssl.stat)

+        }

+        

+        return(ssl)

+        

+      }

+      simGenes <- simplicityGENES(As=object@featureLoadings,

+                                  Ps=object@sampleFactors)

+      

+      patternMarkers <- list()

+      

+      nP <- ncol(simGenes)

+      

+      for (i in 1:nP) {

+        

+        sortSim <- names(sort(simGenes[,i],decreasing=F))

+        

+        geneThresh <- min(which(simGenes[sortSim,i] > 

+                                  apply(simGenes[sortSim,],1,min)))

+        

+        markerGenes <- sortSim[1:geneThresh]

+        

+        markerGenes <- unique(markerGenes)

+        

+        patternMarkers[[i]] <- markerGenes

+        

+        markersByPattern <- patternMarkers

+        

+      }

     }

     else if (threshold == "all") # only the markers with the lowest scores

     {

         patternsByMarker <- colnames(markerScores)[apply(markerScores, 1, which.min)]

         markersByPattern <- sapply(colnames(markerScores), USE.NAMES=TRUE, simplify=FALSE,

             function(pattern) rownames(markerScores)[which(patternsByMarker==pattern)])

+        

     }

     return(list(

         "PatternMarkers"=markersByPattern,

@@ -77,7 +77,7 @@ plot.CogapsResult <- function(x, groups=NULL, ...)

     ylimits <- c(0, max(x@sampleFactors) * 1.1)

     plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude", xlab="", axes=FALSE) 

     axis(1, labels=FALSE, at=1:length(grpset))

-    text(x = seq(1, 12, by=1), par("usr")[3]-0.15, labels = grpset, srt = 60, pos = 1, xpd = TRUE)

+    text(x = seq(1, length(grpset), by=1), par("usr")[3]-0.15, labels = grpset, srt = 60, pos = 1, xpd = TRUE)

     axis(2)

     for (i in 1:nFactors)

@@ -64,25 +64,50 @@ function(object)

 #' @export

 #' @importFrom graphics plot legend lines points

 #' @importFrom grDevices rainbow

-plot.CogapsResult <- function(x, ...)

+plot.CogapsResult <- function(x, groups=NULL, ...)

 {

+  

+  if(!is.null(groups)) {

+    grpset <- unique(groups)

+    groups <- table(groups)

+    ngroups <- length(grpset)

+    nFactors <- ncol(x@sampleFactors)

+    colors <- rainbow(nFactors)

+    xlimits <- c(0, ngroups + 1)

+    ylimits <- c(0, max(x@sampleFactors) * 1.1)

+    plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude", xlab="", axes=FALSE) 

+    axis(1, labels=FALSE, at=1:length(grpset))

+    text(x = seq(1, 12, by=1), par("usr")[3]-0.15, labels = grpset, srt = 60, pos = 1, xpd = TRUE)

+    axis(2)

+    

+    for (i in 1:nFactors)

+    {

+      lines(x=1:length(groups), y=x@sampleFactors[groups,i], col=colors[i])

+      points(x=1:length(groups), y=x@sampleFactors[groups,i], col=colors[i], pch=i)

+    }

+    

+    legend("top", paste("Pattern", 1:nFactors, sep = ""), pch = 1:nFactors,

+           lty=1, cex=0.8, col=colors, bty="y", ncol=5)

+  }

+  else{

     nSamples <- nrow(x@sampleFactors)

     nFactors <- ncol(x@sampleFactors)

     colors <- rainbow(nFactors)

     xlimits <- c(0, nSamples + 1)

     ylimits <- c(0, max(x@sampleFactors) * 1.1)

-

+    

     plot(NULL, xlim=xlimits, ylim=ylimits, ylab="Relative Amplitude",

-        xlab="Samples")

-

+         xlab="Samples")

+    

     for (i in 1:nFactors)

     {

-        lines(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i])

-        points(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i], pch=i)

+      lines(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i])

+      points(x=1:nSamples, y=x@sampleFactors[,i], col=colors[i], pch=i)

     }

-

+    

     legend("top", paste("Pattern", 1:nFactors, sep = ""), pch = 1:nFactors,

-        lty=1, cex=0.8, col=colors, bty="y", ncol=5)

+           lty=1, cex=0.8, col=colors, bty="y", ncol=5)

+  }

 }

 #' @rdname getFeatureLoadings-methods

@@ -437,10 +437,12 @@ function(object, GStoGenes, numPerm, Pw, PwNull)

 #' @param data the original data as a matrix

 #' @param patternPalette a vector indicating what color should be used

 #' for each pattern

+#' @param patternMarkers pattern markers to be plotted, as generated by the 

+#' patternMarkers function

 #' @param sampleNames names of the samples to use for labeling 

 #' @param samplePalette  a vector indicating what color should be used

 #' for each sample

-#' @param colDenogram logical indicating whether to display sample denogram

+#' @param colDendrogram logical indicating whether to display sample dendrogram

 #' @param heatmapCol pallelet giving color scheme for heatmap

 #' @param scale character indicating if the values should be centered and scaled

 #' in either the row direction or the column direction, or none. The

@@ -451,8 +453,8 @@ function(object, GStoGenes, numPerm, Pw, PwNull)

 #' @importFrom gplots bluered

 #' @importFrom gplots heatmap.2

 #' @importFrom stats hclust as.dist cor

-plotPatternMarkers <- function(object, data, patternPalette, sampleNames,

-samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)

+plotPatternMarkers <- function(object, data, patternMarkers, patternPalette, sampleNames,

+samplePalette=NULL, heatmapCol=bluered, colDendrogram=TRUE, scale="row", ...)

 {

     if (is.null(samplePalette))

         samplePalette <- rep("black", ncol(data))

@@ -462,10 +464,10 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)

     names(patternCols) <- unlist(patternMarkers)

     for (i in 1:length(patternMarkers))

     {

-        patternCols[patternMarkers[[i]]] <- patternPalette[i]

+        patternCols[unlist(patternMarkers$PatternMarkers[[i]])] <- patternPalette[i]

     }

-    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers),],...),

+    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers$PatternMarkers),],),

         symkey=FALSE,

         symm=FALSE, 

         scale=scale,

@@ -473,9 +475,9 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)

         distfun=function(x) as.dist((1-cor(t(x)))/2),

         hclustfun=function(x) stats::hclust(x,method="average"),

         Rowv=FALSE,

-        Colv=colDenogram,

+        Colv=colDendrogram,

         trace='none',

-        RowSideColors = as.character(patternCols[unlist(patternMarkers)]),

+        RowSideColors = as.character(patternCols[unlist(patternMarkers$PatternMarkers)]),

         labCol= sampleNames,

         cexCol=0.8,

         ColSideColors=as.character(samplePalette),

@@ -97,8 +97,6 @@ function(object)

 #' @rdname getAmplitudeMatrix-methods

 #' @aliases getAmplitudeMatrix

-#' The method is same as getFeatureLoadings

-#' Alternate name for getFeatureLoadings

 setMethod("getAmplitudeMatrix", signature(object="CogapsResult"),

 function(object)

 {

@@ -117,8 +115,6 @@ function(object)

 #' @rdname getPatternMatrix-methods

 #' @aliases getPatternMatrix

-#' The method is same as getSampleFactors

-#' Alternate name for getSampleFactors

 setMethod("getPatternMatrix", signature(object="CogapsResult"),

 function(object)

 {

@@ -93,6 +93,20 @@ function(object)

     object@featureLoadings

 })

+

+

+#' @rdname getAmplitudeMatrix-methods

+#' @aliases getAmplitudeMatrix

+#' The method is same as getFeatureLoadings

+#' Alternate name for getFeatureLoadings

+setMethod("getAmplitudeMatrix", signature(object="CogapsResult"),

+function(object)

+{

+    object@featureLoadings

+})

+

+

+

 #' @rdname getSampleFactors-methods

 #' @aliases getSampleFactors

 setMethod("getSampleFactors", signature(object="CogapsResult"),

@@ -101,6 +115,17 @@ function(object)

     object@sampleFactors

 })

+#' @rdname getPatternMatrix-methods

+#' @aliases getPatternMatrix

+#' The method is same as getSampleFactors

+#' Alternate name for getSampleFactors

+setMethod("getPatternMatrix", signature(object="CogapsResult"),

+function(object)

+{

+    object@sampleFactors

+})

+

+

 #' @rdname getMeanChiSq-methods

 #' @aliases getMeanChiSq

 setMethod("getMeanChiSq", signature(object="CogapsResult"),

@@ -315,7 +315,7 @@ function(object, sets, whichMatrix, numPerm, ...)

 {

     ## do this for backwards compatibility

     if (!is.null(list(...)$GStoGenes))

-        sets <- GStoGenes

+        sets <- list(...)$GStoGenes

     ## check input arguments

     if (!is(sets, "list"))

         stop("sets must be a list of either measurements or samples")

@@ -56,7 +56,6 @@ function(object)

     nFeatures <- nrow(object@featureLoadings)

     nSamples <- nrow(object@sampleFactors)

     nPatterns <- ncol(object@featureLoadings)

-

     print(paste("CogapsResult object with", nFeatures, "features and", nSamples,

         "samples"))

     print(paste(nPatterns, "patterns were learned"))

@@ -133,7 +132,6 @@ function(object)

 {

     if (!is.null(object@metadata$unmatchedPatterns))

         return(object@metadata$unmatchedPatterns)

-

     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")

     return(NULL)

 })

@@ -145,7 +143,6 @@ function(object)

 {

     if (!is.null(object@metadata$clusteredPatterns))

         return(object@metadata$clusteredPatterns)

-

     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")

     return(NULL)

 })

@@ -157,7 +154,6 @@ function(object)

 {

     if (!is.null(object@metadata$CorrToMeanPattern))

         return(object@metadata$CorrToMeanPattern)

-

     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")

     return(NULL)

 })

@@ -169,7 +165,6 @@ function(object)

 {

     if (!is.null(object@metadata$subsets))

         return(object@metadata$subsets)

-

     message("this result was not generated with a call to GWCoGAPS or scCoGAPS")

     return(NULL)

 })

@@ -179,28 +174,14 @@ function(object)

 setMethod("calcZ", signature(object="CogapsResult"),

 function(object, whichMatrix)

 {

-    if (whichMatrix=="featureLoadings")

-    {

-        if (sum(object@featureStdDev==0) > 0)

-        {

-            warning("zeros detected in the standard deviation matrix")

-            object@featureStdDev[object@featureStdDev==0] <- 1e-6

-        }

-        return(object@featureLoadings / object@featureStdDev)

-    }

-    else if (whichMatrix=="sampleFactors")

-    {

-        if (sum(object@sampleStdDev==0) > 0)

-        {

-            warning("zeros detected in the standard deviation matrix")

-            object@sampleStdDev[object@sampleStdDev==0] <- 1e-6

-        }

-        return(object@sampleFactors / object@sampleStdDev)

-    }

-    else

-    {

+    if (!(whichMatrix %in% c("featureLoadings", "sampleFactors")))

         stop("whichMatrix must be either 'featureLoadings' or 'sampleFactors'")

-    }

+    mean <- if (whichMatrix == "sampleFactors") object@sampleFactors else object@featureLoadings

+    stddev <- if (whichMatrix == "sampleFactors") object@factorStdDev else object@loadingStdDev

+    if (sum(stddev == 0) > 0)

+        warning("zeros detected in the standard deviation matrix")

+    stddev[stddev == 0] <- 1e-6

+    return(mean / stddev)    

 })

 #' @rdname reconstructGene-methods

@@ -384,11 +365,11 @@ function(object, GStoGenes, numPerm, Pw, nullGenes)

     if (nullGenes)

     {

         ZD <- object@featureLoadings[setdiff(row.names(object@featureLoadings), GStoGenes),] /

-            object@featureStdDev[setdiff(row.names(object@featureLoadings), GStoGenes),]

+            object@loadingStdDev[setdiff(row.names(object@featureLoadings), GStoGenes),]

     }

     else

     {

-        ZD <- object@featureLoadings[GStoGenes,]/object@featureStdDev[GStoGenes,]

+        ZD <- object@featureLoadings[GStoGenes,]/object@loadingStdDev[GStoGenes,]

     }

     outStats <- apply(sweep(ZD,2,gsStat,FUN="*"),1,sum) / (sum(gsStat))

     outStats <- outStats / apply(ZD,1,sum)

@@ -258,127 +258,73 @@ function(object, data, uncertainty)

         lwid=c(1,10), lhei=c(1, 4, 1.2 ), main="Heatmap of Residuals")

 })

-patternMarkersCalculation <- function(mat, threshold, lp)

+unitVector <- function(n, length)

 {

+    vec <- rep(0, length)

+    vec[n] <- 1

+    return(vec)

+}

+

+#' @rdname patternMarkers-methods

+#' @aliases patternMarkers

+setMethod("patternMarkers", signature(object="CogapsResult"),

+function(object, threshold, lp, axis)

+{

+    ## check inputs to the function

     if (!(threshold %in% c("cut", "all")))

-        stop("invalid thresold argument")

-    if (!is.na(lp) & length(lp) != ncol(mat))

+        stop("threshold must be either 'cut' or 'all'")

+    if (!is.na(lp) & length(lp) != ncol(object@featureLoadings))

         stop("lp length must equal the number of patterns")

-    

-    nMeasurements <- nrow(mat)

-    nPatterns <- ncol(mat)

-    rowMax <- t(apply(mat, 1, function(x) x / max(x)))

+    if (!(axis %in% 1:2))

+        stop("axis must be either 1 or 2")

+    ## need to scale each row of the matrix of interest so that the maximum is 1

+    resultMatrix <- if (axis == 1) object@featureLoadings else object@sampleFactors

+    normedMatrix <- t(apply(resultMatrix, 1, function(row) row / max(row)))

+    ## handle the case where a custom linear combination of patterns was passed in

     if (!is.na(lp))

     {

-        sstat <- apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp)))

-        ssranks <- rank(sstat)

-        ssgenes.th <- names(sort(sstat,decreasing=FALSE,na.last=TRUE))

+        markerScores <- apply(normedMatrix, 1, function(row) sqrt(sum((row-lp)^2)))

+        markersByPattern <- names(sort(markerScores, decreasing=FALSE, na.last=TRUE))

+        return(list(

+            "PatternMarkers"=markersByPattern,

+            "PatternMarkerRanks"=rank(markerScores),

+            "PatternMarkerScores"=markerScores

+        ))

     }

-    else

-    {

-        # determine which measurements are most associated with each pattern

-        sstat <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=dimnames(mat))

-        ssranks <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=dimnames(mat))

-        ssgenes <- matrix(NA, nrow=nMeasurements, ncol=nPatterns, dimnames=NULL)

-        for (i in 1:nPatterns)

+    ## default pattern marker calculation, each pattern has unit weight

+    markerScores <- sapply(1:ncol(normedMatrix), function(patternIndex)

+        apply(normedMatrix, 1, function(row)

         {

-            lp <- rep(0,nPatterns)

-            lp[i] <- 1

-            sstat[,i] <- unlist(apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp))))

-            ssranks[,i] <- rank(sstat[,i])

-            ssgenes[,i] <- names(sort(sstat[,i], decreasing=FALSE, na.last=TRUE))

-        }

-        if (threshold=="cut")

-        {

-            geneThresh <- sapply(1:nPatterns, function(x) min(which(ssranks[ssgenes[,x],x] > apply(ssranks[ssgenes[,x],],1,min))))

-            ssgenes.th <- sapply(1:nPatterns, function(x) ssgenes[1:geneThresh[x],x])

-        }

-        else if (threshold=="all")

+            lp <- unitVector(patternIndex, ncol(normedMatrix))

+            return(sqrt(sum((row-lp)^2)))

+        })

+    )

+    markerRanks <- apply(markerScores, 2, rank)

+    colnames(markerScores) <- colnames(markerRanks) <- colnames(normedMatrix)

+    ## keep only a subset of markers for each pattern depending on the type of threshold

+    if (threshold == "cut") # all markers which achieve minimum rank

+    {

+        rankCutoff <- sapply(1:ncol(markerRanks), function(patternIndex)

         {

-            pIndx <- unlist(apply(sstat, 1, which.min))

-            gBYp <- list()

-            for (i in sort(unique(pIndx)))

-            {

-                gBYp[[i]] <- sapply(strsplit(names(pIndx[pIndx==i]),"[.]"),function(x) x[[1]][1])

-            }

-            ssgenes.th <- lapply(1:max(sort(unique(pIndx))), function(x)

-                ssgenes[which(ssgenes[,x] %in% gBYp[[x]]),x])

-        }

+            patternRank <- markerRanks[,patternIndex]

+            return(max(patternRank[patternRank == apply(markerRanks, 1, min)]))

+        })

+        markersByPattern <- lapply(1:ncol(markerRanks), function(patternIndex)

+            rownames(markerRanks)[markerRanks[,patternIndex] <= rankCutoff[patternIndex]])

+        names(markersByPattern) <- colnames(markerRanks)

     }

-    return(list("PatternMarkers"=ssgenes.th, "PatternRanks"=ssranks,

-        "PatternMarkerScores"=sstat))       

-}

-

-#' @rdname patternMarkers-methods

-#' @aliases patternMarkers

-setMethod("patternMarkers", signature(object="CogapsResult"),

-function(object, threshold, lp)

-{

-    patternMarkersCalculation(object@featureLoadings, threshold, lp)

-})

-

-#' @rdname amplitudeMarkers-methods

-#' @aliases amplitudeMarkers

-setMethod("amplitudeMarkers", signature(object="CogapsResult"),

-function(object, threshold, lp)

-{

-    patternMarkersCalculation(object@sampleFactors, threshold, lp)

-})

-

-#' heatmap of original data clustered by pattern markers statistic

-#' @export

-#' @docType methods

-#' @rdname plotPatternMarkers-methods

-#'

-#' @param object an object of type CogapsResult

-#' @param data the original data as a matrix

-#' @param patternPalette a vector indicating what color should be used

-#' for each pattern

-#' @param sampleNames names of the samples to use for labeling 

-#' @param samplePalette  a vector indicating what color should be used

-#' for each sample

-#' @param colDenogram logical indicating whether to display sample denogram

-#' @param heatmapCol pallelet giving color scheme for heatmap

-#' @param scale character indicating if the values should be centered and scaled

-#' in either the row direction or the column direction, or none. The

-#' default is "row". 

-#' @param ... additional graphical parameters to be passed to \code{heatmap.2}

-#' @return heatmap of the \code{data} values for the \code{patternMarkers}

-#' @seealso  \code{\link{heatmap.2}}

-#' @importFrom gplots bluered

-#' @importFrom gplots heatmap.2

-#' @importFrom stats hclust as.dist cor

-plotPatternMarkers <- function(object, data, patternPalette, sampleNames,

-samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)

-{

-    if (is.null(samplePalette))

-        samplePalette <- rep("black", ncol(data))

-

-    ### coloring of genes

-    patternCols <- rep(NA, length(unlist(patternMarkers)))

-    names(patternCols) <- unlist(patternMarkers)

-    for (i in 1:length(patternMarkers))

+    else if (threshold == "all") # only the markers with the lowest scores

     {

-        patternCols[patternMarkers[[i]]] <- patternPalette[i]

+        patternsByMarker <- colnames(markerScores)[apply(markerScores, 1, which.min)]

+        markersByPattern <- sapply(colnames(markerScores), USE.NAMES=TRUE, simplify=FALSE,

+            function(pattern) rownames(markerScores)[which(patternsByMarker==pattern)])

     }

-

-    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers),],...),

-        symkey=FALSE,

-        symm=FALSE, 

-        scale=scale,

-        col=heatmapCol,

-        distfun=function(x) as.dist((1-cor(t(x)))/2),

-        hclustfun=function(x) stats::hclust(x,method="average"),

-        Rowv=FALSE,

-        Colv=colDenogram,

-        trace='none',

-        RowSideColors = as.character(patternCols[unlist(patternMarkers)]),

-        labCol= sampleNames,

-        cexCol=0.8,

-        ColSideColors=as.character(samplePalette),

-        rowsep=cumsum(sapply(patternMarkers,length))

-    )

-}

+    return(list(

+        "PatternMarkers"=markersByPattern,

+        "PatternMarkerRanks"=markerRanks,

+        "PatternMarkerScores"=markerScores

+    ))

+})

 #' @rdname calcCoGAPSStat-methods

 #' @aliases calcCoGAPSStat

@@ -394,20 +340,18 @@ function(object, sets, whichMatrix, numPerm, ...)

         stop("sets must be a list of either measurements or samples")

     ## do a permutation test

     zMatrix <- calcZ(object, whichMatrix)

-    pvalUpReg <- sapply(sets,

-        function(thisSet)

+    pvalUpReg <- sapply(sets, function(thisSet)

+    {

+        lessThanCount <- rep(0, ncol(zMatrix))

+        actualZScore <- colMeans(zMatrix[rownames(zMatrix) %in% thisSet,])

+        for (n in 1:numPerm)

         {

-            numValues <- length(thisSet)

-            lessThanCount <- rep(0, ncol(zMatrix))

-            actualZScore <- colMeans(zMatrix[rownames(zMatrix) %in% thisSet,])

-            for (n in 1:numPerm)

-            {

-                permutedIndices <- sample(1:nrow(zMatrix), size=numValues, replace=FALSE)

-                permutedZScore <- colMeans(zMatrix[permutedIndices,])

-                lessThanCount <- lessThanCount + (actualZScore < permutedZScore)

-            }

-            return(lessThanCount / numPerm)

-        })

+            permutedIndices <- sample(1:nrow(zMatrix), size=length(thisSet), replace=FALSE)

+            permutedZScore <- colMeans(zMatrix[permutedIndices,])

+            lessThanCount <- lessThanCount + (actualZScore < permutedZScore)

+        }

+        return(lessThanCount / numPerm)

+    })

     ## calculate other quantities of interest, return a list

     pvalDownReg <- 1 - pvalUpReg # this is techinically >=, but == should rarely happen

     activityEstimate <- 1 - 2 * pvalUpReg

@@ -482,4 +426,57 @@ function(object, GStoGenes, numPerm, Pw, PwNull)

     return(finalStats)

 })

+#' heatmap of original data clustered by pattern markers statistic

+#' @export

+#' @docType methods

+#' @rdname plotPatternMarkers-methods

+#'

+#' @param object an object of type CogapsResult

+#' @param data the original data as a matrix

+#' @param patternPalette a vector indicating what color should be used

+#' for each pattern

+#' @param sampleNames names of the samples to use for labeling 

+#' @param samplePalette  a vector indicating what color should be used

+#' for each sample

+#' @param colDenogram logical indicating whether to display sample denogram

+#' @param heatmapCol pallelet giving color scheme for heatmap

+#' @param scale character indicating if the values should be centered and scaled

+#' in either the row direction or the column direction, or none. The

+#' default is "row". 

+#' @param ... additional graphical parameters to be passed to \code{heatmap.2}

+#' @return heatmap of the \code{data} values for the \code{patternMarkers}

+#' @seealso  \code{\link{heatmap.2}}

+#' @importFrom gplots bluered

+#' @importFrom gplots heatmap.2

+#' @importFrom stats hclust as.dist cor

+plotPatternMarkers <- function(object, data, patternPalette, sampleNames,

+samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)

+{

+    if (is.null(samplePalette))

+        samplePalette <- rep("black", ncol(data))

+

+    ### coloring of genes

+    patternCols <- rep(NA, length(unlist(patternMarkers)))

+    names(patternCols) <- unlist(patternMarkers)

+    for (i in 1:length(patternMarkers))

+    {

+        patternCols[patternMarkers[[i]]] <- patternPalette[i]

+    }

+    gplots::heatmap.2(as.matrix(data[unlist(patternMarkers),],...),

+        symkey=FALSE,

+        symm=FALSE, 

+        scale=scale,

+        col=heatmapCol,

+        distfun=function(x) as.dist((1-cor(t(x)))/2),

+        hclustfun=function(x) stats::hclust(x,method="average"),

+        Rowv=FALSE,

+        Colv=colDenogram,

+        trace='none',

+        RowSideColors = as.character(patternCols[unlist(patternMarkers)]),

+        labCol= sampleNames,

+        cexCol=0.8,

+        ColSideColors=as.character(samplePalette),

+        rowsep=cumsum(sapply(patternMarkers,length))

+    )

+}

@@ -18,6 +18,34 @@ createCogapsResult <- function(returnList, allParams)

         diagnostics = append(returnList$diagnostics,

             list("params"=allParams$gaps, "version"=utils::packageVersion("CoGAPS")))

     )

+

+    # label snapshots

+    if (allParams$nSnapshots > 0)

+    {

+        freq <- floor(allParams$gaps@nIterations / allParams$nSnapshots)

+        labels <- seq(freq, allParams$gaps@nIterations, freq)

+        if (allParams$snapshotPhase == 'equilibration' | allParams$snapshotPhase == 'all')

+        {

+            for (n in 1:allParams$nSnapshots)

+            {

+                dimnames(res@metadata$equilibrationSnapshotsA[[n]]) <- dimnames(res@featureLoadings)

+                dimnames(res@metadata$equilibrationSnapshotsP[[n]]) <- dimnames(res@sampleFactors)

+            }

+            names(res@metadata$equilibrationSnapshotsA) <- labels

+            names(res@metadata$equilibrationSnapshotsP) <- labels

+        }

+        if (allParams$snapshotPhase == 'sampling' | allParams$snapshotPhase == 'all')

+        {

+            for (n in 1:allParams$nSnapshots)

+            {

+                dimnames(res@metadata$samplingSnapshotsA[[n]]) <- dimnames(res@featureLoadings)

+                dimnames(res@metadata$samplingSnapshotsP[[n]]) <- dimnames(res@sampleFactors)

+            }

+            names(res@metadata$samplingSnapshotsA) <- labels

+            names(res@metadata$samplingSnapshotsP) <- labels

+        }

+    }        

+    

     validObject(res)

     return(res)

 }

@@ -232,16 +260,16 @@ function(object, data, uncertainty)

 patternMarkersCalculation <- function(mat, threshold, lp)

 {

+    if (!(threshold %in% c("cut", "all")))

+        stop("invalid thresold argument")

+    if (!is.na(lp) & length(lp) != ncol(mat))

+        stop("lp length must equal the number of patterns")

+    

     nMeasurements <- nrow(mat)

     nPatterns <- ncol(mat)

     rowMax <- t(apply(mat, 1, function(x) x / max(x)))

-    

     if (!is.na(lp))

     {

-        if (length(lp) != nPatterns)

-        {

-            warning("lp length must equal the number of columns of the Amatrix")

-        }

         sstat <- apply(rowMax, 1, function(x) sqrt(t(x-lp) %*% (x-lp)))

         ssranks <- rank(sstat)

         ssgenes.th <- names(sort(sstat,decreasing=FALSE,na.last=TRUE))

@@ -276,10 +304,6 @@ patternMarkersCalculation <- function(mat, threshold, lp)

             ssgenes.th <- lapply(1:max(sort(unique(pIndx))), function(x)

                 ssgenes[which(ssgenes[,x] %in% gBYp[[x]]),x])

         }

-        else

-        {

-            stop("invalid thresold argument")

-        }

     }

     return(list("PatternMarkers"=ssgenes.th, "PatternRanks"=ssranks,

         "PatternMarkerScores"=sstat))       

@@ -151,17 +151,22 @@ function(object)

 setMethod("calcZ", signature(object="CogapsResult"),

 function(object, whichMatrix)

 {

-

     if (whichMatrix=="featureLoadings")

     {

         if (sum(object@featureStdDev==0) > 0)

-            stop("zeros detected in the standard deviation matrix")

+        {

+            warning("zeros detected in the standard deviation matrix")

+            object@featureStdDev[object@featureStdDev==0] <- 1e-6

+        }

         return(object@featureLoadings / object@featureStdDev)

     }

     else if (whichMatrix=="sampleFactors")

     {

         if (sum(object@sampleStdDev==0) > 0)

-            stop("zeros detected in the standard deviation matrix")

+        {

+            warning("zeros detected in the standard deviation matrix")

+            object@sampleStdDev[object@sampleStdDev==0] <- 1e-6

+        }

         return(object@sampleFactors / object@sampleStdDev)

     }

     else

@@ -383,7 +383,7 @@ function(object, sets, whichMatrix, numPerm, ...)

     pvalDownReg <- 1 - pvalUpReg # this is techinically >=, but == should rarely happen

     activityEstimate <- 1 - 2 * pvalUpReg

     return(list(

-        'twoSidedPValue'=pmin(pvalDownReg, pvalUpReg),

+        'twoSidedPValue'=pmax(pmin(pvalDownReg, pvalUpReg), 1 / numPerm),

         'GSUpreg'=pvalUpReg,

         'GSDownreg'=pvalDownReg,

         'GSActEst'=activityEstimate

@@ -151,12 +151,23 @@ function(object)

 setMethod("calcZ", signature(object="CogapsResult"),

 function(object, whichMatrix)

 {

+

     if (whichMatrix=="featureLoadings")

+    {

+        if (sum(object@featureStdDev==0) > 0)

+            stop("zeros detected in the standard deviation matrix")

         return(object@featureLoadings / object@featureStdDev)

+    }

     else if (whichMatrix=="sampleFactors")

+    {

+        if (sum(object@sampleStdDev==0) > 0)

+            stop("zeros detected in the standard deviation matrix")

         return(object@sampleFactors / object@sampleStdDev)

+    }

     else

+    {

         stop("whichMatrix must be either 'featureLoadings' or 'sampleFactors'")

+    }

 })

 #' @rdname reconstructGene-methods

@@ -352,8 +363,6 @@ function(object, sets, whichMatrix, numPerm, ...)

     ## check input arguments

     if (!is(sets, "list"))

         stop("sets must be a list of either measurements or samples")

-    if (sum(object@featureStdDev==0) > 0)

-        stop("zeros detected in the standard deviation matrix")

     ## do a permutation test

     zMatrix <- calcZ(object, whichMatrix)

     pvalUpReg <- sapply(sets,

@@ -149,20 +149,14 @@ function(object)

 #' @rdname calcZ-methods

 #' @aliases calcZ

 setMethod("calcZ", signature(object="CogapsResult"),

-function(object, which)

+function(object, whichMatrix)

 {

-    if (which == "feature")

-    {

-        object@featureLoadings / object@featureStdDev

-    }

-    else if (which == "sample")

-    {

-        object@sampleFactors / object@sampleStdDev

-    }

+    if (whichMatrix=="featureLoadings")

+        return(object@featureLoadings / object@featureStdDev)

+    else if (whichMatrix=="sampleFactors")

+        return(object@sampleFactors / object@sampleStdDev)

     else

-    {

-        stop("\"which\" must be either \"feature\" or \"sample\"")

-    }

+        stop("whichMatrix must be either 'featureLoadings' or 'sampleFactors'")

 })

 #' @rdname reconstructGene-methods

@@ -348,118 +342,43 @@ samplePalette=NULL, heatmapCol=bluered, colDenogram=TRUE, scale="row", ...)

 #' @rdname calcCoGAPSStat-methods

 #' @aliases calcCoGAPSStat

-#' @importFrom stats runif

 #' @importFrom methods is

 setMethod("calcCoGAPSStat", signature(object="CogapsResult"),

-function(object, GStoGenes, numPerm)

+function(object, sets, whichMatrix, numPerm, ...)

 {

-    # test for std dev of zero, possible mostly in simple simulations

+    ## do this for backwards compatibility

+    if (!is.null(list(...)$GStoGenes))

+        sets <- GStoGenes

+    ## check input arguments

+    if (!is(sets, "list"))

+        stop("sets must be a list of either measurements or samples")

     if (sum(object@featureStdDev==0) > 0)

-    {

-        object@featureStdDev[object@featureStdDev==0] <- 1e-6

-    }

-

-    # calculate Z scores

-    zMatrix <- calcZ(object)

-

-    # check input arguments

-    if (!is(GStoGenes, "data.frame") && !is(GStoGenes, "list") && !is(GStoGenes,"GSA.genesets"))

-    {

-        stop("GStoGenes must be a data.frame,GSA.genesets, or list with format specified in the users manual.")

-    }

-

-    if (is(GStoGenes, "GSA.genesets"))

-    {

-        names(GStoGenes$genesets) <- GStoGenes$geneset.names

-        GStoGenes <- GStoGenes$genesets

-    }

-

-    if (is(GStoGenes, "list"))

-    {

-        GStoGenesList <- GStoGenes

-    }

-    else

-    {

-        GStoGenesList <- list()

-        for (i in 1:dim(GStoGenes)[2])

+        stop("zeros detected in the standard deviation matrix")

+    ## do a permutation test

+    zMatrix <- calcZ(object, whichMatrix)

+    pvalUpReg <- sapply(sets,

+        function(thisSet)

         {

-            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   <- list()

-    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 (!any(names(results)==label))

-        {

-            for (p in 1:numPatt)