deleted file mode 100644

@@ -1 +0,0 @@

-load("/Users/maggiewodicka/Desktop/CoGAPS/GWCoGAPS.AP.fixed.Rda")

@@ -1,9 +1,10 @@

 Package: CoGAPS

-Version: 2.7.0

-Date: 2014-08-23

+Version: 2.99.0

+Date: 2018-04-24

 Title: Coordinated Gene Activity in Pattern Sets

 Author: Thomas Sherman, Wai-shing Lee, Conor Kelton, Ondrej Maxian, Jacob Carey,

-    Genevieve Stein-O'Brien, Michael Considine, Maggie Wodicka, John Stansfield, Shawn Sivy, Carlo Colantuoni, Alexander Favorov, Mike Ochs, Elana Fertig

+    Genevieve Stein-O'Brien, Michael Considine, Maggie Wodicka, John Stansfield,

+    Shawn Sivy, Carlo Colantuoni, Alexander Favorov, Mike Ochs, Elana Fertig

 Description: Coordinated Gene Activity in Pattern Sets (CoGAPS)

     implements a Bayesian MCMC matrix factorization algorithm,

     GAPS, and links it to gene set statistic methods to infer biological

@@ -41,4 +42,4 @@ biocViews: GeneExpression, Transcription, GeneSetEnrichment,

     DifferentialExpression, Bayesian, Clustering, TimeCourse, RNASeq, Microarray,

     MultipleComparison, DimensionReduction

 NeedsCompilation: yes

-RoxygenNote: 6.0.1

+RoxygenNote: 5.0.1

@@ -26,6 +26,7 @@ export(reOrderBySet)

 export(reconstructGene)

 export(residuals)

 export(scCoGAPS)

+export(scCoGAPSFromCheckpoint)

 import(doParallel)

 import(foreach)

 import(ggplot2)

@@ -24,8 +24,6 @@

 #' @param fixedPatterns matrix of fixed values in either A or P matrix

 #' @param checkpointInterval time (in seconds) between creating a checkpoint

 #' @param checkpointFile name of the checkpoint file

-#' @param pumpThreshold type of threshold for pump statistic

-#' @param nPumpSamples number of samples used in pump statistic

 #' @param ... keeps backwards compatibility with arguments from older versions

 #' @return list with A and P matrix estimates

 #' @importFrom methods new

@@ -37,8 +35,7 @@ CoGAPS <- function(D, S, nFactor=7, nEquil=1000, nSample=1000, nOutputs=1000,

 nSnapshots=0, alphaA=0.01, alphaP=0.01, maxGibbmassA=100, maxGibbmassP=100,

 seed=-1, messages=TRUE, singleCellRNASeq=FALSE, whichMatrixFixed='N',

 fixedPatterns=matrix(0), checkpointInterval=0, 

-checkpointFile="gaps_checkpoint.out", pumpThreshold="unique",

-nPumpSamples=0, ...)

+checkpointFile="gaps_checkpoint.out", ...)

 {

     # get v2 arguments

     oldArgs <- list(...)

@@ -57,6 +54,14 @@ nPumpSamples=0, ...)

     if (missing(S) & !is.null(oldArgs$unc))

         S <- oldArgs$unc

+    # get pump arguments - hidden for now from user

+    pumpThreshold <- "unique"

+    nPumpSamples <- 0

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

+        pumpThreshold <- list(...)$pumpThreshold

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

+        pumpThreshold <- list(...)$nPumpSamples

+

     # check arguments

     if (class(D) != "matrix" | class(S) != "matrix")

         stop('D and S must be matrices')

@@ -149,7 +154,7 @@ output_atomic=FALSE, fixedBinProbs=FALSE, fixedDomain="N", sampleSnapshots=TRUE,

 numSnapshots=100, alphaA=0.01, nMaxA=100000, max_gibbmass_paraA=100.0,

 alphaP=0.01, nMaxP=100000, max_gibbmass_paraP=100.0, seed=-1, messages=TRUE)

 {

-    #warning('gapsRun is deprecated with v3.0, use CoGAPS')

+    warning('gapsRun is deprecated with v3.0, use CoGAPS')

     CoGAPS(D, S, nFactor=nFactor, nEquil=nEquil, nSample=nSample,

         nOutputs=nOutR, nSnapshots=ifelse(sampleSnapshots,numSnapshots,0),

         alphaA=alphaA, alphaP=alphaP, maxGibbmassA=max_gibbmass_paraA,

@@ -179,7 +184,7 @@ sampleSnapshots=TRUE, numSnapshots=100, alphaA=0.01, nMaxA=100000,

 max_gibbmass_paraA=100.0, alphaP=0.01, nMaxP=100000, max_gibbmass_paraP=100.0,

 seed=-1, messages=TRUE)

 {

-    #warning('gapsMapRun is deprecated with v3.0, use CoGaps')

+    warning('gapsMapRun is deprecated with v3.0, use CoGaps')

     CoGAPS(D, S, nFactor=nFactor, nEquil=nEquil, nSample=nSample,

         nOutputs=nOutR, nSnapshots=ifelse(sampleSnapshots,numSnapshots,0),

         alphaA=alphaA, alphaP=alphaP, maxGibbmassA=max_gibbmass_paraA,

@@ -192,7 +197,7 @@ v2CoGAPS <- function(result, ...)

 {

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

     {

-        #warning('GStoGenes is deprecated with v3.0, see CoGAPS documentation')

+        warning('GStoGenes is deprecated with v3.0, see CoGAPS documentation')

         if (is.null(list(...)$plot) | list(...)$plot)

         {

             plotGAPS(result$Amean, result$Pmean)

@@ -18,16 +18,21 @@

 #' data(SimpSim)

 #' sim_name <- "example"

 #' createGWCoGAPSSets(SimpSim.D, SimpSim.S, nSets=2, sim_name)

-#' result <- GWCoGAPS(sim_name, nFactor=3, nEquil=1000, nSample=1000)

+#' result <- GWCoGAPS(sim_name, nFactor=3, nEquil=200, nSample=200)

 #' @export

 GWCoGAPS <- function(simulationName, nFactor, nCores=NA, cut=NA, minNS=NA, manualMatch=FALSE, consensusPatterns=NULL, ...)

 {

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

+    {

         stop("checkpoint file name automatically set in GWCoGAPS - don't pass this parameter")

-    if (is.null(consensusPatterns)){

+    }

+

+    if (is.null(consensusPatterns))

+    {

         allDataSets <- preInitialPhase(simulationName, nCores)

         initialResult <- runInitialPhase(simulationName, allDataSets, nFactor, ...)

-        if(manualMatch){

+        if (manualMatch)

+        {

             saveRDS(initialResult,file=paste(simulationName,"_initial.rds", sep=""))

             stop("Please provide consensus patterns upon restarting.")

         }

@@ -107,7 +112,9 @@ preInitialPhase <- function(simulationName, nCores)

     # establish the number of cores that we are able to use

     if (is.na(nCores))

+    {

         nCores <- length(allDataSets)

+    }

     registerDoParallel(cores=nCores)

     return(allDataSets)

 }

@@ -139,21 +146,29 @@ postInitialPhase <- function(initialResult, nSets, cut, minNS)

     #run postpattern match function

     if (is.na(cut))

+    {

         cut <- nFactor

+    }

+

     return(patternMatch4Parallel(Ptot=BySet$P, nP=nFactor, nSets=nSets, cnt=cut,

         minNS=minNS, bySet=TRUE))

 }

 runFinalPhase <- function(simulationName, allDataSets, consensusPatterns, nCores, ...)

 {    

-    if(length(dim(consensusPatterns))!=2){stop("consensusPatterns must be a matrix")}

+    if (length(dim(consensusPatterns)) != 2)

+    {

+        stop("consensusPatterns must be a matrix")

+    }

     # find data files if providing consensus patterns

     fileSig <- paste(simulationName, "_partition_[0-9]+.RData", sep="")

     allDataSets <- list.files(full.names=TRUE, pattern=fileSig)

     if (is.na(nCores))

-    nCores <- length(allDataSets)

+    {

+        nCores <- length(allDataSets)

+    }

     registerDoParallel(cores=nCores)

     # generate seeds for parallelization

@@ -1,95 +1,121 @@

-

-#' patternMatch4Parallel

+#' cellMatchR

+#' @export

 #'

 #' @param Atot a matrix containing the total by set estimates of Pmean output from \code{reOrderBySet}

 #' @param nSets number of parallel sets used to generate \code{Atot}

 #' @param cnt  number of branches at which to cut dendrogram

 #' @param minNS minimum of individual set contributions a cluster must contain

+#' @param maxNS maximum of individual set contributions a cluster must contain

 #' @param ignore.NA logical indicating whether or not to ignore NAs from potential over dimensionalization. Default is FALSE.

 #' @param bySet logical indicating whether to return list of matched set solutions from \code{Atot}

-#' @param R2mean 

-#' @param

+#' @param plotDen plot

 #' @param ... additional parameters for \code{agnes}

 #' @return a matrix of concensus patterns by samples. If \code{bySet=TRUE} then a list of the set contributions to each

 #' concensus pattern is also returned.

-#' @export

-#' @seealso \code{\link{fastcluster}}

-#'

-#'

+cellMatchR <- function(Atot, nSets, cnt, minNS=NA, maxNS=NA,

+ignore.NA=FALSE, bySet=FALSE, plotDen=FALSE, ...)

+{

+    if (is.na(minNS))

+    {

+        minNS <- nSets / 2

+    }

+    if (is.na(maxNS))

+    {

+        maxNS <- nSets + minNS

+    }

+    if (!ignore.NA)

+    {

+        if (anyNA(Atot))

+        {

+            warning(paste("Non-sparse matrixes produced. Reducing the number",

+                "of patterns asked for and rerun."))

+        }

+    }

+    else

+    {

+        Atot <- Atot[complete.cases(Atot),]

+    }

-cellMatchR <- function(Atot,nSets, cnt, minNS=NULL, maxNS=NULL, ignore.NA=FALSE, bySet=FALSE, plotDen=FALSE,...){

+    corcut <- function(Atot, minNS, cnt, cluster.method)

+    {

+        corr.dist <- cor(Atot)

+        corr.dist <- 1-corr.dist

-  if(is.null(minNS)){minNS=nSets/2}

-  if(is.null(maxNS)){maxNS=nSets+minNS}

+        clust <- agnes(x=corr.dist, diss=TRUE, cluster.method)

+        #clust=fastcluster::hclust(dist(corr.dist))

+        cut <- cutree(as.hclust(clust),k=cnt)

-  if(ignore.NA==FALSE){if(anyNA(Atot)){

-    warning('Non-sparse matrixes produced. Reducing the number of patterns asked for and rerun.')

-  }}

-  if(ignore.NA==TRUE){Atot<-Atot[complete.cases(Atot),]}

+        cls <- sort(unique(cut))

+        cMNs <- matrix(ncol=cnt,nrow=dim(Atot)[1])

+        colnames(cMNs) <- cls

+        rownames(cMNs) <- rownames(Atot)

-corcut<-function(Atot,minNS,cnt,cluster.method){

-  corr.dist=cor(Atot)

-  corr.dist=1-corr.dist

+        RtoMeanPattern <- list()

+        AByClust <- list()

+        for (i in cls)

+        {

+            if (!is.null(dim(Atot[,cut == i])))

+            {

+                if (dim(Atot[,cut == i])[2] >= minNS)

+                {

+                    cMNs[,i] <- rowMeans(Atot[,cut==i])

+                    AByClust[[i]] <- Atot[,cut==i]

+                    nIN <- sum(cut==i)

+                    RtoMeanPattern[[i]] <- sapply(1:nIN, function(j) 

+                        round(cor(x=Atot[,cut==i][,j], y=cMNs[,i]),3))

+                }

+            }

+        }

-  clust=agnes(x=corr.dist,diss=TRUE,cluster.method)

-  #clust=fastcluster::hclust(dist(corr.dist))

-  cut=cutree(as.hclust(clust),k=cnt)

+        AByClust[sapply(AByClust,is.null)] <- NULL

+        RtoMeanPattern[sapply(RtoMeanPattern,is.null)] <- NULL

+        return(list("RtoMeanPattern"=RtoMeanPattern, "AByClust"=AByClust))

+    }   

-  cls=sort(unique(cut))

-  cMNs=matrix(ncol=cnt,nrow=dim(Atot)[1])

-  colnames(cMNs)=cls

-  rownames(cMNs)=rownames(Atot)

+    cc <- corcut(Atot, minNS, cnt, cluster.method)

-  RtoMeanPattern <- list()

-  AByClust <- list()

-    for(i in cls){

-        if (is.null(dim(Atot[,cut == i]))==TRUE){

-          next

-        } else if(dim(Atot[,cut == i])[2] < minNS){

+    ### split by maxNS

+    indx <- which(unlist(lapply(cc$AByClust, function(x) dim(x)[1] > maxNS)))

+    i <- 1

+    while (length(indx) > 0)

+    { 

+        icc <- corcut(cc$AByClust[[indx[1]]], minNS, 2, cluster.method)

+        if (length(icc[[2]]) == 0)

+        {

+            indx <- indx[-1]

             next

-        } else{

-            cMNs[,i]=rowMeans(Atot[,cut==i])

-            AByClust[[i]] <- Atot[,cut==i]

-            nIN=sum(cut==i)

-            RtoMeanPattern[[i]] <- sapply(1:nIN,function(j) {round(cor(x=Atot[,cut==i][,j],y=cMNs[,i]),3)})

         }

-      }

-    AByClust[sapply(AByClust,is.null)]<-NULL

-    RtoMeanPattern[sapply(RtoMeanPattern,is.null)]<-NULL

-    return(list("RtoMeanPattern"=RtoMeanPattern,"AByClust"=AByClust))

-  }   

-  cc<-corcut(Atot,minNS,cnt,cluster.method)

-

-    ### split by maxNS

-    indx<-which(unlist(lapply(cc$AByClust,function(x) dim(x)[1]>maxNS)))

-    i<-1

-    while(length(indx)>0){ 

-            icc<-corcut(cc$AByClust[[indx[1]]],minNS,2,cluster.method)

-            if(length(icc[[2]])==0){

-              indx<-indx[-1]

-              next

-            } else{

-              cc$AByClust[[indx[1]]]<-icc[[2]][[1]]

-              cc$RtoMeanPattern[[indx[1]]]<-icc[[1]][[1]]

-              if(length(icc[[2]])>1){

+        else

+        {

+            cc$AByClust[[indx[1]]] <- icc[[2]][[1]]

+            cc$RtoMeanPattern[[indx[1]]] <- icc[[1]][[1]]

+            if (length(icc[[2]]) > 1)

+            {

                 cc$AByClust<-append(cc$AByClust,icc[[2]][2])

                 cc$RtoMeanPattern<-append(cc$RtoMeanPattern,icc[[1]][2])

-              } 

-              indx<-which(unlist(lapply(cc$AByClust,function(x) dim(x)[1]>maxNS)))

-            }

+            } 

+            indx <- which(unlist(lapply(cc$AByClust,function(x) dim(x)[1]>maxNS)))

+        }

     }

+    #weighted.mean(AByClustDrop[[1]],RtoMPDrop[[1]])

+    AByCDSWavg <- t(sapply(1:length(cc$AByClust), function(z)

+        apply(cc$AByClust[[z]], 1, function(x)

+            weighted.mean(x, (cc$RtoMeanPattern[[z]])^3))))

+    rownames(AByCDSWavg) <- lapply(1:length(cc$AByClust), function(x)

+        paste("Pattern",x))

+    #scale As

+    Amax <- apply(AByCDSWavg, 1, max)

+    AByCDSWavgScaled <- t(sapply(1:dim(AByCDSWavg)[1], function(x)

+        AByCDSWavg[x,] / Amax[x]))

+    rownames(AByCDSWavgScaled) <- rownames(AByCDSWavg)

-#weighted.mean(AByClustDrop[[1]],RtoMPDrop[[1]])

-AByCDSWavg<- t(sapply(1:length(cc$AByClust),function(z) apply(cc$AByClust[[z]],1,function(x) weighted.mean(x,(cc$RtoMeanPattern[[z]])^3))))

-rownames(AByCDSWavg) <- lapply(1:length(cc$AByClust),function(x) paste("Pattern",x))

-#scale As

-Amax <- apply(AByCDSWavg,1,max)

-AByCDSWavgScaled <- t(sapply(1:dim(AByCDSWavg)[1],function(x) AByCDSWavg[x,]/Amax[x]))

-rownames(AByCDSWavgScaled) <- rownames(AByCDSWavg)

-

-  if(bySet){  

-    return(list("consenusAs"=t(AByCDSWavgScaled),"ABySet"=cc))

-  } else {return(AByCDSWavgScaled)}

-

+    if (bySet)

+    {  

+        return(list("consenusAs"=t(AByCDSWavgScaled),"ABySet"=cc))

+    }

+    else

+    {

+        return(AByCDSWavgScaled)

+    }

 }

@@ -1,42 +1,39 @@

-#' createGWCoGAPSSets

+#' Create Gene Sets for GWCoGAPS

 #'

-#'\code{createGWCoGAPSSets} factors whole genome data into randomly generated sets for indexing;

+#' @details factors whole genome data into randomly generated sets for indexing

 #'

-#'@param data data matrix with unique rownames

-#'@param nSets number of sets for parallelization

-#'@param outRDA name of output file

-#'@param keep logical indicating whether or not to save gene set list. Default is TRUE.

-#'@param path character string indicating were to save resulting data objects

-#'@export

-#'@return list with randomly generated sets of genes from whole genome data

-#'@examples \dontrun{

-#'createGWCoGAPSSet(D,nSets=nSets)

-#'}

-#'

-

-createGWCoGAPSSets <- function(D, S, nSets, simulationName, path="")

+#' @param D data matrix

+#' @param S uncertainty matrix

+#' @param nSets number of sets to partition the data into

+#' @param simulationName name used to identify files created by this simulation

+#' @return simulationName used to identify saved files

+#' @examples

+#' data(SimpSim)

+#' createGWCoGAPSSets(SimpSim.D, SimpSim.S, nSets=2, "example")

+#' @export

+createGWCoGAPSSets <- function(D, S, nSets, simulationName)

 {

     # check gene names

-    if (length(unique(colnames(D))) != length(colnames(D)))

+    if (length(unique(rownames(D))) != length(rownames(D)))

     {

-        warning("Cell identifiers not unique!")

+        warning("Gene identifiers not unique!")

     }

-    # partition data by sampling random sets of cells

+    # partition data by sampling random sets of genes

     genes <- 1:nrow(D)

     setSize <- floor(length(genes) / nSets)

     for (set in 1:nSets)

     {

-        

-        # sample genes

-            sampleSize <- ifelse(set == nSets, length(genes), setSize)

-            geneset <- sample(genes, sampleSize, replace=FALSE)

-            genes <- genes[!(genes %in% geneset)]

+        # sample gene names

+        sampleSize <- ifelse(set == nSets, length(genes), setSize)

+        geneSet <- sample(genes, sampleSize, replace=FALSE)

+        genes <- genes[!(genes %in% geneSet)]

+

         # partition data

-        sampleD <- D[geneset,]

-        sampleS <- S[geneset,]

-        save(sampleD, sampleS, file=paste(path, simulationName, "_partition_", set,

+        sampleD <- D[geneSet,]

+        sampleS <- S[geneSet,]

+        save(sampleD, sampleS, file=paste(simulationName, "_partition_", set,

             ".RData", sep=""));

     }

     return(simulationName)

-}

+}

\ No newline at end of file

@@ -1,9 +1,8 @@

 #' Create Gene Sets for scCoGAPS

+#' @export

 #'

-#' @details factors whole genome data into randomly generated sets for indexing

-#'

+#' @description factors whole genome data into randomly generated sets for indexing

 #' @param D data matrix

-#' @param S uncertainty matrix

 #' @param nSets number of sets to partition the data into

 #' @param simulationName name used to identify files created by this simulation

 #' @param samplingRatio vector of relative quantities to use for sampling celltypes

@@ -12,9 +11,9 @@

 #' @return simulationName used to identify saved files

 #' @examples

 #' data(SimpSim)

-#' createscCoGAPSSets(SimpSim.D, SimpSim.S, nSets=2, "example")

-#' @export

-createscCoGAPSSets <- function(D, nSets, simulationName,samplingRatio=NULL,path="",anotionObj=NULL)

+#' createscCoGAPSSets(SimpSim.D, nSets=2, simulationName="example")

+createscCoGAPSSets <- function(D, nSets, simulationName, samplingRatio=NULL,

+path="", anotionObj=NULL)

 {

     # check gene names

     if (length(unique(colnames(D))) != length(colnames(D)))

@@ -27,16 +26,22 @@ createscCoGAPSSets <- function(D, nSets, simulationName,samplingRatio=NULL,path=

     setSize <- floor(length(cells) / nSets)

     for (set in 1:nSets)

     {

-        

-        if(is.null(samplingRatio)){

-        # sample cell names

+        if (is.null(samplingRatio))

+        {

+            # sample cell names

             sampleSize <- ifelse(set == nSets, length(cells), setSize)

             cellset <- sample(cells, sampleSize, replace=FALSE)

             cells <- cells[!(cells %in% cellset)]

-        } else {

-        if(length(unique(anotionObj))!=length(samplingRatio)){warning("Not all celltypes will be sampled from.")}

-        ct.indx<-lapply(unique(anotionObj),function(x) which(anotionObj == x))

-        cellset<-sample(colnames(D)[ct.indx[[x]]], samplingRatio[x],replace=TRUE)

+        }

+        else

+        {

+            if (length(unique(anotionObj)) != length(samplingRatio))

+            {

+                warning("Not all celltypes will be sampled from.")

+            }

+            ct.indx <- lapply(unique(anotionObj), function(x) which(anotionObj == x))

+            cellset <- lapply(unique(anotionObj), function(x)

+                sample(colnames(D)[ct.indx[[x]]], samplingRatio[x],replace=TRUE))

         }

         # partition data

@@ -45,7 +50,8 @@ createscCoGAPSSets <- function(D, nSets, simulationName,samplingRatio=NULL,path=

         sampleD <- log2(sampleD+1)

         # generate S

         sampleS <- pmax(.1*sampleD, .1)

-        save(sampleD, sampleS, file=paste0(path,simulationName, "_partition_", set,".RData"));

+        save(sampleD, sampleS, file=paste0(path,simulationName, "_partition_",

+            set,".RData"));

     }

     return(simulationName)

 }

@@ -12,56 +12,75 @@

 patternMarkers <- function(Amatrix=NA, scaledPmatrix=FALSE, Pmatrix=NA,

 threshold="all", lp=NA, full=FALSE)

 {

-    if(scaledPmatrix==FALSE)

+    if (scaledPmatrix==FALSE)

     {

-        if(!is.na(Pmatrix)){

-        pscale <- apply(Pmatrix,1,max)   # rescale p's to have max 1

-        Amatrix <- sweep(Amatrix, 2, pscale, FUN="*")   # rescale A in accordance with p's having max 1

+        if(!is.na(Pmatrix))

+        {

+            pscale <- apply(Pmatrix,1,max)   # rescale p's to have max 1

+            Amatrix <- sweep(Amatrix, 2, pscale, FUN="*")   # rescale A in accordance with p's having max 1

+        }

+        else

+        {

+            warning("P values must be provided if not already scaled")

         }

-        else(warning("P values must be provided if not already scaled"))

     }

     # find the A with the highest magnitude

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

-    if(!is.na(lp))

+    if (!is.na(lp))

     {

-        if(length(lp)!=dim(Amatrix)[2]){

+        if (length(lp) != dim(Amatrix)[2])

+        {

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

         }

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

         ssranks <-rank(sstat)

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

-    } else {

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

+    }

+    else

+    {

         # determine which genes are most associated with each pattern

         sstat<-matrix(NA, nrow=nrow(Amatrix), ncol=ncol(Amatrix),dimnames=dimnames(Amatrix))

         ssranks<-matrix(NA, nrow=nrow(Amatrix), ncol=ncol(Amatrix),dimnames=dimnames(Amatrix))#list()

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

-        for(i in 1:nP){

+        for (i in 1:nP)

+        {

             lp <- rep(0,dim(Amatrix)[2])

             lp[i] <- 1

             sstat[,i] <- unlist(apply(Arowmax, 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"){

+        if (threshold=="cut")

+        {

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

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

         }

-        else if (threshold=="all"){

+        else if (threshold=="all")

+        {

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

             gBYp <- list()

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

+            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]

-            })

-        } else{stop("Threshold arguement not viable option")}

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

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

+        }

+        else

+        {

+            stop("Threshold arguement not viable option")

+        }

     }

     if (full)

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

+    {

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

+            "PatternMarkerScores"=sstat))

+    }

     else

+    {

         return("PatternMarkers"=ssgenes.th)

+    }

 }

@@ -1,96 +1,129 @@

 #' patternMatch4Parallel

 #'

-#' @param Ptot a matrix containing the total by set estimates of Pmean output from \code{reOrderBySet}

+#' @param Ptot a matrix containing the total by set estimates of Pmean output

+#' from \code{reOrderBySet}

 #' @param nSets number of parallel sets used to generate \code{Ptot}

 #' @param cnt  number of branches at which to cut dendrogram

 #' @param minNS minimum of individual set contributions a cluster must contain

-#' @param maxNS max of individual set contributions a cluster must contain. default is nSets+minNS

+#' @param maxNS max of individual set contributions a cluster must contain.

+#' default is nSets+minNS

 #' @param cluster.method the agglomeration method to be used for clustering

-#' @param ignore.NA logical indicating whether or not to ignore NAs from potential over dimensionalization. Default is FALSE.

-#' @param bySet logical indicating whether to return list of matched set solutions from \code{Ptot}

+#' @param ignore.NA logical indicating whether or not to ignore NAs from

+#' potential over dimensionalization. Default is FALSE.

+#' @param bySet logical indicating whether to return list of matched set

+#' solutions from \code{Ptot}

 #' @param ... additional parameters for \code{agnes}

-#' @return a matrix of concensus patterns by samples. If \code{bySet=TRUE} then a list of the set contributions to each

+#' @return a matrix of concensus patterns by samples. If \code{bySet=TRUE} then

+#' a list of the set contributions to each

 #' concensus pattern is also returned.

 #' @seealso \code{\link{agnes}}

 #' @export

-patternMatch4Parallel <- function(Ptot, nSets, cnt, minNS=NULL, maxNS=NULL, 

-    cluster.method="complete", ignore.NA=FALSE, bySet=FALSE, ...){

-

-    if (is.null(minNS)){minNS=ceiling(nSets/2)}

-    if (is.null(maxNS)){maxNS=nSets+minNS}

-

+patternMatch4Parallel <- function(Ptot, nSets, cnt, minNS=NA, maxNS=NULL, 

+cluster.method="complete", ignore.NA=FALSE, bySet=FALSE, ...)

+{

+    if (is.na(minNS))

+    {

+        minNS <- ceiling(nSets / 2)

+    }

+    if (is.null(maxNS))

+    {

+        maxNS=nSets+minNS

+    }

     if (ignore.NA==FALSE & anyNA(Ptot))

-        warning('Non-sparse matrixes produced. Reducing the number of patterns asked for and rerun.')

+    {

+        warning(paste("Non-sparse matrixes produced. Reducing the number of",

+            "patterns asked for and rerun."))

+    }

     if (ignore.NA==TRUE)

+    {

         Ptot <- Ptot[complete.cases(Ptot),]

+    }

-    corcut<-function(Ptot,minNS,cnt,cluster.method){

-    # corr dist

-    corr.dist=cor(t(Ptot))

-    corr.dist=1-corr.dist

-    # cluster

-    #library(cluster)

-    clust=agnes(x=corr.dist,diss=TRUE,method=cluster.method)

-    cut=cutree(as.hclust(clust),k=cnt)

-    #save.image(file=paste("CoGAPS.",nP,"P.",nS,"Set.CorrClustCut",cnt,".RData"))

+    corcut <-function(Ptot,minNS,cnt,cluster.method)

+    {

+        # corr dist

+        corr.dist <- cor(t(Ptot))

+        corr.dist <- 1-corr.dist

+        # cluster

+        #library(cluster)

+        clust <- agnes(x=corr.dist, diss=TRUE, method=cluster.method)

+        cut <- cutree(as.hclust(clust), k=cnt)

+        #save.image(file=paste("CoGAPS.", nP, "P.", nS, "Set.CorrClustCut",

+            #cnt,".RData"))

-    cls=sort(unique(cut))

-    cMNs=matrix(nrow=cnt,ncol=dim(Ptot)[2])

-    rownames(cMNs)=cls

-    colnames(cMNs)=colnames(Ptot)

+        cls <- sort(unique(cut))

+        cMNs <- matrix(nrow=cnt, ncol=dim(Ptot)[2])

+        rownames(cMNs) <- cls

+        colnames(cMNs) <- colnames(Ptot)

-    RtoMeanPattern <- list()

-    PByClust <- list()

-    for(i in cls){

-        if (is.null(dim(Ptot[cut == i, ]))==TRUE){

-            next

-        } else if(dim(Ptot[cut == i, ])[1]< minNS){

-            next

-        } else {

-            cMNs[i,]=colMeans(Ptot[cut==i,])

-            PByClust[[i]] <- Ptot[cut==i,]

-            nIN=sum(cut==i)

-            RtoMeanPattern[[i]] <- sapply(1:nIN,function(j) {round(cor(x=Ptot[cut==i,][j,],y=cMNs[i,]),3)})

+        RtoMeanPattern <- list()

+        PByClust <- list()

+        for (i in cls)

+        {

+            if (!is.null(dim(Ptot[cut == i,])))

+            {

+                if (dim(Ptot[cut == i,])[1] >= minNS)

+                {

+                    cMNs[i,] <- colMeans(Ptot[cut==i,])

+                    PByClust[[i]] <- Ptot[cut==i,]

+                    nIN <- sum(cut==i)

+                    RtoMeanPattern[[i]] <- sapply(1:nIN, function(j)

+                        round(cor(x=Ptot[cut==i,][j,], y=cMNs[i,]),3)

+                    )

+                }

+            }

         }

-    }

-    PByClust[sapply(PByClust,is.null)]<-NULL

-    RtoMeanPattern[sapply(RtoMeanPattern,is.null)]<-NULL

-    return(list("RtoMeanPattern"=RtoMeanPattern,"PByClust"=PByClust))

+

+        PByClust[sapply(PByClust,is.null)] <- NULL

+        RtoMeanPattern[sapply(RtoMeanPattern,is.null)] <- NULL

+        return(list("RtoMeanPattern"=RtoMeanPattern, "PByClust"=PByClust))

     }    

-    cc<-corcut(Ptot,minNS,cnt,cluster.method)

+    cc <- corcut(Ptot,minNS,cnt,cluster.method)

     ### split by maxNS

-    indx<-which(unlist(lapply(cc$PByClust,function(x) dim(x)[1]>maxNS)))

-    while(length(indx)>0){

-            icc<-corcut(cc$PByClust[[indx[1]]],minNS,2,cluster.method)

-            if(length(icc[[2]])==0){

-              indx<-indx[-1]

-              next

-            } else{

-                cc$PByClust[[indx[1]]]<-icc[[2]][[1]]

-                cc$RtoMeanPattern[[indx[1]]]<-icc[[1]][[1]]

-                if(length(icc[[2]])>1){

-                    cc$PByClust<-append(cc$PByClust,icc[[2]][2])

-                    cc$RtoMeanPattern<-append(cc$RtoMeanPattern,icc[[1]][2])

-                } 

-                indx<-which(unlist(lapply(cc$PByClust,function(x) dim(x)[1]>maxNS)))

-            }

+    indx <- which(unlist(lapply(cc$PByClust,function(x) dim(x)[1]>maxNS)))

+    while(length(indx) > 0)

+    {

+        icc <- corcut(cc$PByClust[[indx[1]]],minNS,2,cluster.method)

+        if (length(icc[[2]])==0)

+        {

+          indx <- indx[-1]

+          next

+        }

+        else

+        {

+            cc$PByClust[[indx[1]]] <- icc[[2]][[1]]

+            cc$RtoMeanPattern[[indx[1]]] <- icc[[1]][[1]]

+            if (length(icc[[2]]) > 1)

+            {

+                cc$PByClust <- append(cc$PByClust,icc[[2]][2])

+                cc$RtoMeanPattern <- append(cc$RtoMeanPattern,icc[[1]][2])

+            } 

+            indx <- which(unlist(lapply(cc$PByClust, function(x)

+                dim(x)[1] > maxNS)))

+        }

     }

     #weighted.mean

-    PByCDSWavg<- t(sapply(1:length(cc$PByClust),function(z) apply(cc$PByClust[[z]],2,function(x) 

-        weighted.mean(x,(cc$RtoMeanPattern[[z]])^3))))

-    rownames(PByCDSWavg) <- lapply(1:length(cc$PByClust),function(x) paste("Pattern",x))

+    PByCDSWavg <- t(sapply(1:length(cc$PByClust), function(z)

+        apply(cc$PByClust[[z]],2, function(x)

+            weighted.mean(x, (cc$RtoMeanPattern[[z]])^3))))

+    rownames(PByCDSWavg) <- lapply(1:length(cc$PByClust), function(x)

+        paste("Pattern", x))

     #scale ps

     Pmax <- apply(PByCDSWavg,1,max)

-    PByCDSWavgScaled <- t(sapply(1:dim(PByCDSWavg)[1],function(x) PByCDSWavg[x,]/Pmax[x]))

+    PByCDSWavgScaled <- t(sapply(1:dim(PByCDSWavg)[1], function(x)

+        PByCDSWavg[x,] / Pmax[x]))

     rownames(PByCDSWavgScaled) <- rownames(PByCDSWavg)

-    if(bySet){

-        return(list("consenusPatterns"=PByCDSWavgScaled,"PBySet"=cc))

-    } else {

+    if (bySet)

+    {

+        return(list("consenusPatterns"=PByCDSWavgScaled, "PBySet"=cc))

+    }

+    else

+    {

         return("consenusPatterns"=PByCDSWavgScaled)

     }

 }

@@ -12,51 +12,56 @@

 #' concensus pattern is also returned.

 #' @seealso \code{\link{agnes}}

 #' @export

-#' 

-#' 

-

 patternMatch4singleCell <- function(Atot, nSets, cnt, minNS, 

 cluster.method="complete", ignore.NA=FALSE, bySet=FALSE, ...)

 {

-    if (!is.null(minNS))

-        minNS=nSets/2

-

-    if (ignore.NA==FALSE & anyNA(Atot))

-        warning('Non-sparse matrixes produced. Reducing the number of patterns asked for and rerun.')

-    if (ignore.NA==TRUE)

+    if (is.na(minNS))

+    {

+        minNS <- nSets/2

+    }

+    if (!ignore.NA & anyNA(Atot))

+    {

+        warning(paste("Non-sparse matrixes produced. Reducing the number of",

+            "patterns asked for and rerun"))

+    }

+    if (ignore.NA)

+    {

         Atot <- Atot[complete.cases(Atot),]

+    }

     # corr dist

-    corr.dist=cor(Atot)

-    corr.dist=1-corr.dist

+    corr.dist <- cor(Atot)

+    corr.dist <- 1 - corr.dist

     # cluster

     #library(cluster)

-    clust=agnes(x=corr.dist,diss=TRUE,method=cluster.method)

-    cut=cutree(as.hclust(clust),k=cnt)

+    clust <- agnes(x=corr.dist, diss=TRUE, method=cluster.method)

+    cut <- cutree(as.hclust(clust), k=cnt)

     #save.image(file=paste("CoGAPS.",nP,"P.",nS,"Set.CorrClustCut",cnt,".RData"))

     #drop n<4 and get weighted Avg

-    cls=sort(unique(cut))

-    cMNs=matrix(nrow=cnt,ncol=dim(Atot)[1])

-    rownames(cMNs)=cls

-    colnames(cMNs)=colnames(Atot)

+    cls <- sort(unique(cut))

+    cMNs <- matrix(nrow=cnt, ncol=dim(Atot)[1])

+    rownames(cMNs) <- cls

+    colnames(cMNs) <- rownames(Atot)

     RtoMeanPattern <- list()

     PByClust <- list()

     for(i in cls)

     {

-        if (is.null(dim(Atot[,cut == i ]))==TRUE)

+        if (is.null(dim(Atot[,cut == i])))

         {

-            cMNs[i,] <- Atot[,cut == i ]

+            cMNs[i,] <- Atot[,cut == i]

             RtoMeanPattern[[i]] <- rep(1,length(Atot[,cut == i ]))

             PByClust[[i]] <- t(as.matrix(Atot[,cut == i ]))

         }

         else

         {

-            cMNs[i,]=colMeans(Atot[,cut==i])

+            cMNs[i,] <- rowMeans(Atot[,cut==i])

             PByClust[[i]] <- Atot[,cut==i]

-            nIN=sum(cut==i)

-            RtoMeanPattern[[i]] <- sapply(1:nIN,function(j) {round(cor(x=Atot[,cut==i][j,],y=cMNs[i,]),3)})

+            nIN <- sum(cut==i)

+            RtoMeanPattern[[i]] <- sapply(1:nIN, function(j)

+                round(cor(x=Atot[,cut==i][,j], y=cMNs[i,]), 3)

+            )

         }

     }

@@ -65,17 +70,14 @@ cluster.method="complete", ignore.NA=FALSE, bySet=FALSE, ...)

     RtoMPDrop <- list()

     for(i in cls)

     {

-        if (is.null(dim(PByClust[[i]])) == TRUE)

-            next