deleted file mode 100644

@@ -1 +0,0 @@

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

@@ -1,16 +1,18 @@

 Package: CoGAPS

-Version: 2.99.0

-Date: 2014-08-23

+Version: 2.99.2

+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

     process activity.  It can be used to perform sparse matrix factorization on

     any data, and when this data represents biomolecules, to do gene set

     analysis.

-Maintainer: Elana J. Fertig <ejfertig@jhmi.edu>

+Maintainer: Elana J. Fertig <ejfertig@jhmi.edu>,

+    Thomas D. Sherman <tomsherman159@gmail.com>

 Depends:

     R (>= 3.4.0),

     Rcpp (>= 0.11.0)

@@ -1,7 +1,6 @@

 # Generated by roxygen2: do not edit by hand

 export(CoGAPS)

-export(CoGapsFromCheckpoint)

 export(GWCoGAPS)

 export(GWCoGapsFromCheckpoint)

 export(binaryA)

@@ -9,16 +8,21 @@ export(buildReport)

 export(calcCoGAPSStat)

 export(calcGeneGSStat)

 export(calcZ)

+export(cellMatchR)

 export(computeGeneGSProb)

 export(createGWCoGAPSSets)

+export(createscCoGAPSSets)

 export(gapsMapRun)

 export(gapsRun)

+export(patternMarkers)

+export(patternMatch4Parallel)

 export(plotAtoms)

 export(plotDiag)

 export(plotGAPS)

 export(plotP)

 export(reconstructGene)

 export(residuals)

+export(scCoGAPS)

 import(doParallel)

 import(foreach)

 import(ggplot2)

@@ -34,8 +34,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

@@ -47,8 +45,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(...)

@@ -67,6 +64,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')

@@ -115,7 +120,6 @@ nPumpSamples=0, ...)

 #' @param path path to checkpoint file

 #' @param checkpointFile name for future checkpooints made

 #' @return list with A and P matrix estimates

-#' @export

 CoGapsFromCheckpoint <- function(D, S, path, checkpointFile=NA)

 {

     if (is.na(checkpointFile))

@@ -129,7 +133,7 @@ CoGapsFromCheckpoint <- function(D, S, path, checkpointFile=NA)

 #'  compiler/version was used, which compile time options were enabled, etc...

 #' @return display builds information

 #' @examples

-#'  CoGAPS::displayBuildReport()

+#'  CoGAPS::buildReport()

 #' @export

 buildReport <- function()

 {

@@ -166,7 +170,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,

@@ -196,7 +200,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,

@@ -209,7 +213,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)

@@ -9,6 +9,8 @@

 #' @param nCores number of cores for parallelization. If left to the default NA, nCores = nSets.

 #' @param cut number of branches at which to cut dendrogram used in patternMatch4Parallel

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

+#' @param manualMatch logical indicating whether or not to stop after initial phase for manual pattern matching

+#' @param consensusPatterns fixed pattern matrix to be used to ensure reciprocity of A weights accross sets 

 #' @param ... additional parameters to be fed into \code{gapsRun} and \code{gapsMapRun}

 #' @return list of A and P estimates

 #' @seealso \code{\link{gapsRun}}, \code{\link{patternMatch4Parallel}}, and \code{\link{gapsMapRun}}

@@ -16,18 +18,30 @@

 #' 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, ...)

+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")

-    allDataSets <- preInitialPhase(simulationName, nCores)

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

-    matchedPatternSets <- postInitialPhase(initialResult, length(allDataSets), cut, minNS)

-    save(matchedPatternSets, file=paste(simulationName, "_matched_ps.RData", sep=""))

-    finalResult <- runFinalPhase(simulationName, allDataSets, matchedPatternSets, ...)

-    return(postFinalPhase(finalResult, matchedPatternSets))

+    }

+

+    if (is.null(consensusPatterns))

+    {

+        allDataSets <- preInitialPhase(simulationName, nCores)

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

+        if (manualMatch)

+        {

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

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

+        }

+        matchedPatternSets <- postInitialPhase(initialResult, length(allDataSets), cut, minNS)

+        save(matchedPatternSets, file=paste(simulationName, "_matched_ps.RData", sep=""))

+        consensusPatterns <- matchedPatternSets[[1]]

+    } 

+    finalResult <- runFinalPhase(simulationName, allDataSets, consensusPatterns, nCores, ...)

+    return(postFinalPhase(finalResult, consensusPatterns))

 }

 #' Restart a GWCoGaps Run from Checkpoint

@@ -35,15 +49,23 @@ GWCoGAPS <- function(simulationName, nFactor, nCores=NA, cut=NA, minNS=NA, ...)

 #' @inheritParams GWCoGAPS

 #' @return list of A and P estimates

 #' @importFrom utils file_test

+#' @examples

+#' data(SimpSim)

+#' sim_name <- "example"

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

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

+#' result <- GWCoGapsFromCheckpoint(sim_name, 2)

 #' @export

-GWCoGapsFromCheckpoint <- function(simulationName, nCores=NA, cut=NA, minNS=NA, ...)

+GWCoGapsFromCheckpoint <- function(simulationName, nCores, cut=NA, minNS=NA, ...)

 {

     # find data files

     allDataSets <- preInitialPhase(simulationName, nCores)

     # figure out phase from file signature

-    initialCpts <- list.files(full.names=TRUE, pattern=paste(simulationName, "_initial_cpt_[0-9]+.out", sep=""))

-    finalCpts <- list.files(full.names=TRUE, pattern=paste(simulationName, "_final_cpt_[0-9]+.out", sep=""))

+    initialCpts <- list.files(full.names=TRUE, pattern=paste(simulationName,

+        "_initial_cpt_[0-9]+.out", sep=""))

+    finalCpts <- list.files(full.names=TRUE, pattern=paste(simulationName,

+        "_final_cpt_[0-9]+.out", sep=""))

     if (length(finalCpts))

     {

@@ -51,10 +73,12 @@ GWCoGapsFromCheckpoint <- function(simulationName, nCores=NA, cut=NA, minNS=NA,

         {

             # load data set and shift values so gene minimum is zero

             load(allDataSets[[i]])

-            sampleD <- sweep(sampleD, 1, apply(sampleD, 1, function(x) pmin(0,min(x))))

+            sampleD <- sweep(sampleD, 1, apply(sampleD, 1, function(x)

+                pmin(0,min(x))))

             # run CoGAPS with fixed patterns

-            cptFileName <- paste(simulationName, "_final_cpt_", i, ".out", sep="")

+            cptFileName <- paste(simulationName, "_final_cpt_", i, ".out",

+                sep="")

             CoGapsFromCheckpoint(sampleD, sampleS, cptFileName)

         }

         load(paste(simulationName, "_matched_ps.RData", sep=""))

@@ -62,7 +86,9 @@ GWCoGapsFromCheckpoint <- function(simulationName, nCores=NA, cut=NA, minNS=NA,

     else if (file_test("-f", paste(simulationName, "_matched_ps.RData", sep="")))

     {

         load(paste(simulationName, "_matched_ps.RData", sep=""))

-        finalResult <- runFinalPhase(simulationName, allDataSets, matchedPatternSets, ...)

+        consensusPatterns<-matchedPatternSets[[1]]

+        finalResult <- runFinalPhase(simulationName, allDataSets,

+            consensusPatterns, nCores, ...)

     }

     else if (length(initialCpts))

     {

@@ -71,21 +97,27 @@ GWCoGapsFromCheckpoint <- function(simulationName, nCores=NA, cut=NA, minNS=NA,

         {

             # load data set and shift values so gene minimum is zero

             load(allDataSets[[i]])

-            sampleD <- sweep(sampleD, 1, apply(sampleD, 1, function(x) pmin(0,min(x))))

+            sampleD <- sweep(sampleD, 1, apply(sampleD, 1, function(x)

+                pmin(0,min(x))))

             # run CoGAPS from checkpoint

-            cptFileName <- paste(simulationName, "_initial_cpt_", i, ".out", sep="")

+            cptFileName <- paste(simulationName, "_initial_cpt_", i, ".out",

+                sep="")

             CoGapsFromCheckpoint(sampleD, sampleS, cptFileName)

         }

-        matchedPatternSets <- postInitialPhase(initialResult, length(allDataSets), cut, minNS)

-        save(matchedPatternSets, file=paste(simulationName, "_matched_ps.RData", sep=""))

-        finalResult <- runFinalPhase(simulationName, allDataSets, matchedPatternSets, ...)

+        matchedPatternSets <- postInitialPhase(initialResult,

+            length(allDataSets), cut, minNS)

+        save(matchedPatternSets, file=paste(simulationName, "_matched_ps.RData",

+            sep=""))

+        consensusPatterns<-matchedPatternSets[[1]]

+        finalResult <- runFinalPhase(simulationName, allDataSets,

+            consensusPatterns, nCores, ...)

     }

     else

     {

         stop("no checkpoint files found")

     }

-    return(postFinalPhase(finalResult, matchedPatternSets))

+    return(postFinalPhase(finalResult, consensusPatterns))

 }

 preInitialPhase <- function(simulationName, nCores)

@@ -96,7 +128,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)

 }

@@ -128,18 +162,36 @@ 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, matchedPatternSets, ...)

-{

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

+{    

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

+    }

+    registerDoParallel(cores=nCores)

+

     # generate seeds for parallelization

     nut <- generateSeeds(chains=length(allDataSets), seed=-1)

     # final number of factors

-    nFactorFinal <- nrow(matchedPatternSets[[1]])

+    nFactorFinal <- nrow(consensusPatterns)

     # run fixed CoGAPS

     finalResult <- foreach(i=1:length(allDataSets)) %dopar%

@@ -150,7 +202,7 @@ runFinalPhase <- function(simulationName, allDataSets, matchedPatternSets, ...)

         # run CoGAPS with fixed patterns

         cptFileName <- paste(simulationName, "_final_cpt_", i, ".out", sep="")

-        CoGAPS(sampleD, sampleS, fixedPatterns=matchedPatternSets[[1]],

+        CoGAPS(sampleD, sampleS, fixedPatterns=consensusPatterns,

             nFactor=nFactorFinal, seed=nut[i], checkpointFile=cptFileName,

             whichMatrixFixed='P', ...)

@@ -158,11 +210,10 @@ runFinalPhase <- function(simulationName, allDataSets, matchedPatternSets, ...)

     return(finalResult)

 }

-postFinalPhase <- function(finalResult, matchedPatternSets)

+postFinalPhase <- function(finalResult, consensusPatterns)

 {

-    Aresult <- postFixed4Parallel(finalResult, matchedPatternSets[[1]])

-    finalResult <- list("Amean"=Aresult$A, "Asd"=Aresult$Asd,

-        "Pmean"=matchedPatternSets, "PbySet"=matchedPatternSets[["PBySet"]])

+    Aresult <- postFixed4Parallel(finalResult, consensusPatterns)

+    finalResult <- list("Amean"=Aresult$A, "Asd"=Aresult$Asd,"Pmean"=consensusPatterns)

     class(finalResult) <- append(class(finalResult), "CoGAPS")

     return(finalResult)

 }

@@ -11,8 +11,8 @@

 #' @param nullGenes logical indicating gene adjustment

 #' @return gene similiarity statistic

 #' @examples

-#' data('SimpSim')

-#' calcGeneGSStat(SimpSim.result$Amean, SimpSim.result$Asd, GStoGenes=GSets[[1]],

+#' data("SimpSim")

+#' calcGeneGSStat(SimpSim.result$Amean, SimpSim.result$Asd, GSGenes=GSets[[1]],

 #' numPerm=500)

 #' @export

 calcGeneGSStat  <- function(Amean, Asd, GSGenes, numPerm, Pw=rep(1,ncol(Amean)),

@@ -67,8 +67,8 @@ nullGenes=FALSE)

 #' @return A vector of length GSGenes containing the p-values of set membership

 #' for each gene containined in the set specified in GSGenes.

 #' @examples

-#' data('SimpSim')

-#' computeGeneGSProb(SimpSim.result$Amean, SimpSim.result$Asd, GStoGenes=GSets[[1]],

+#' data("SimpSim")

+#' computeGeneGSProb(SimpSim.result$Amean, SimpSim.result$Asd, GSGenes=GSets[[1]],

 #' numPerm=500)

 #' @export

 computeGeneGSProb <- function(Amean, Asd, GSGenes, Pw=rep(1,ncol(Amean)),

new file mode 100644

@@ -0,0 +1,121 @@

+#' 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 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.

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

+    }

+

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

+    {

+        corr.dist <- cor(Atot)

+        corr.dist <- 1-corr.dist

+

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

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

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

+

+        cls <- sort(unique(cut))

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

+        colnames(cMNs) <- cls

+        rownames(cMNs) <- rownames(Atot)

+

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

+                }

+            }

+        }

+

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

+            {

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

+        }

+    }

+

+    #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)

+    }

+}

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

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

     }

     return(simulationName)

-}

+}

\ No newline at end of file

new file mode 100644

@@ -0,0 +1,57 @@

+#' Create Gene Sets for scCoGAPS

+#' @export

+#'

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

+#' @param D data 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

+#' @param anotionObj vector of same length as number of columns of D 

+#' @param path character string indicating were to save resulting data objects. default is current working dir

+#' @return simulationName used to identify saved files

+#' @examples

+#' data(SimpSim)

+#' 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)))

+    {

+        warning("Cell identifiers not unique!")

+    }

+

+    # partition data by sampling random sets of cells

+    cells <- 1:ncol(D)

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

+    for (set in 1:nSets)

+    {

+        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 <- lapply(unique(anotionObj), function(x)

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

+        }

+

+        # partition data

+        sampleD <- D[,cellset]

+        #log transform 

+        sampleD <- log2(sampleD+1)

+        # generate S

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

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

+            set,".RData"));

+    }

+    return(simulationName)

+}

@@ -8,70 +8,79 @@

 #' @param full logical indicating whether to return the ranks of each gene for each pattern

 #' @return By default a non-overlapping list of genes associated with each \code{lp}. If \code{full=TRUE} a data.frame of

 #' genes rankings with a column for each \code{lp} will also be returned.

+#' @export

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

-    pmax<-apply(Amatrix, 1, max)

-    # determine which genes are most associated with each pattern

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

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

-    nP=dim(Amatrix)[2]

-    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")

         }

-        for (i in 1:nP)

-        {

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

-            ssranks[order(sstat),i] <- 1:length(sstat)

-            ssgenes[,i]<-names(sort(sstat,decreasing=FALSE))

-        }

+        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

     {

-        for(i in 1:nP)

+        # 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)

         {

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

             lp[i] <- 1

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

-            ssranks[order(sstat),i] <- 1:length(sstat)

-            ssgenes[,i]<-names(sort(sstat,decreasing=FALSE))

+            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")

+        {

+            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")

+        {

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

+        }

+        else

+        {

+            stop("Threshold arguement not viable option")

         }

     }

-

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

-        #geneThresh <- apply(sweep(ssranks,1,t(apply(ssranks, 1, min)),"-"),2,function(x) which(x==0))

-        #ssgenes.th <- lapply(geneThresh,names)

-    }

-    else if (threshold=="all")

+    if (full)

     {

-        pIndx<-apply(ssranks,1,which.min)

-        gBYp <- lapply(sort(unique(pIndx)),function(x) names(pIndx[pIndx==x]))

-        ssgenes.th <- lapply(1:nP, function(x) ssgenes[which(ssgenes[,x] %in% gBYp[[x]]),x])

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

+            "PatternMarkerScores"=sstat))

     }

     else

     {

-        stop("Threshold arguement not viable option")

-    }

-

-    if (full)

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

-    else

         return("PatternMarkers"=ssgenes.th)

+    }

 }

@@ -1,129 +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 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}}

-patternMatch4Parallel <- function(Ptot, nSets, cnt, minNS, 

+#' @export

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

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

 {

-    if (!is.null(minNS))

-        minNS=nSets/2

-

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

+    }

-    # 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"))

-    #drop n<4 and get weighted Avg

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

+        RtoMeanPattern <- list()

+        PByClust <- list()

+        for (i in cls)

         {

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

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

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

-        }

-        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)})

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

+                    )

+                }

+            }

         }

-    }

-    #drop n<minNS 

-    PByClustDrop <- list()

-    RtoMPDrop <- list()

-    for(i in cls)

-    {

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

-            next

-        if (dim(PByClust[[i]])[1] < minNS)

-        {

-            next

-        }

-        else

-        {

-            #indx <- which(RtoMeanPattern[[i]]>.7,arr.ind = TRUE)

-            PByClustDrop <- append(PByClustDrop,list(PByClust[[i]]))

-            RtoMPDrop <- append(RtoMPDrop,list(RtoMeanPattern[[i]]))

-        }

-    }

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

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

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

+    }    

-    ### split by corr  (build in drop if below minNS)

-    PByCDS <- list()

-    RtoMPDS <- list()

-    for (j in 1:length(PByClustDrop))

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

     {

-        if (is.null(dim(PByClustDrop[[j]]))==TRUE)

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

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

         {

-            next

+          indx <- indx[-1]

+          next

         }

-        if (dim(PByClustDrop[[j]])[1]<minNS+nSets)

-        {

-            PByCDS <- append(PByCDS,PByClustDrop[j])

-            RtoMPDS <- append(RtoMPDS,RtoMPDrop[j])

-        }

-        if (dim(PByClustDrop[[j]])[1]>=minNS+nSets)

+        else

         {

-            corr.distPBCD=cor(t(PByClustDrop[[j]]))

-            corr.distPBCD=1-corr.distPBCD

-            clustPBCD=agnes(x=corr.distPBCD,diss=TRUE,method="complete")

-            cutPBCD=cutree(as.hclust(clustPBCD),k=2)

-            g1 <- PByClustDrop[[j]][cutPBCD==1,]

-            PByCDS <- append(PByCDS,list(g1))

-            RtoMPDS <- append(RtoMPDS,list(sapply(1:dim(g1)[1],function(z) round(cor(x=g1[z,],y=colMeans(PByClustDrop[[j]][cutPBCD==1,])),3))))

-            g2 <- PByClustDrop[[j]][cutPBCD==2,]

-            if (is.null(dim(g2)[1])==FALSE)

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

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

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

             {

-                PByCDS <- append(PByCDS,list(g2))

-                RtoMPDS <- append(RtoMPDS,list(sapply(1:dim(g2)[1],function(z) round(cor(x=g2[z,],y=colMeans(PByClustDrop[[j]][cutPBCD==2,])),3))))

-            }

+                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(PByClustDrop[[1]],RtoMPDrop[[1]])

-    PByCDSWavg<- t(sapply(1:length(PByCDS),function(z) apply(PByCDS[[z]],2,function(x) weighted.mean(x,(RtoMPDS[[z]])^3))))

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

+    #weighted.mean