@@ -23,7 +23,7 @@ exportClasses(CogapsResult)

 importClassesFrom(S4Vectors,Annotated)

 importClassesFrom(S4Vectors,character_OR_NULL)

 importClassesFrom(SingleCellExperiment,LinearEmbeddingMatrix)

-importFrom(BiocParallel,SnowParam)

+importFrom(BiocParallel,MulticoreParam)

 importFrom(BiocParallel,bplapply)

 importFrom(Rcpp,evalCpp)

 importFrom(SummarizedExperiment,assay)

@@ -85,6 +85,11 @@ checkpointInFile=NULL, transposeData=FALSE, ...)

     if (!is(data, "character"))

         checkDataMatrix(data, uncertainty, params)

+    # check single cell parameter

+    if (!is.null(allParams$gaps@distributed))

+        if (allParams$gaps@distributed == "single-cell" & !allParams$gaps@singleCell)

+            warning("running single-cell CoGAPS with singleCell=FALSE")

+

     # convert data to matrix

     if (is(data, "data.frame"))

         data <- data.matrix(data)

@@ -132,7 +137,7 @@ checkpointInFile=NULL, transposeData=FALSE, ...)

 #' @inheritParams CoGAPS

 #' @return CogapsResult object

 #' @importFrom methods new

-scCoGAPS <- function(data, params=new("CogapsParams"), nThreads=NULL,

+scCoGAPS <- function(data, params=new("CogapsParams"), nThreads=1,

 messages=TRUE, outputFrequency=500, uncertainty=NULL,

 checkpointOutFile="gaps_checkpoint.out", checkpointInterval=1000,

 checkpointInFile=NULL, transposeData=FALSE, ...)

@@ -151,7 +156,7 @@ checkpointInFile=NULL, transposeData=FALSE, ...)

 #' @inheritParams CoGAPS

 #' @return CogapsResult object

 #' @importFrom methods new

-GWCoGAPS <- function(data, params=new("CogapsParams"), nThreads=NULL,

+GWCoGAPS <- function(data, params=new("CogapsParams"), nThreads=1,

 messages=TRUE, outputFrequency=500, uncertainty=NULL,

 checkpointOutFile="gaps_checkpoint.out", checkpointInterval=1000,

 checkpointInFile=NULL, transposeData=FALSE, ...)

@@ -7,8 +7,8 @@

 #' @param allParams list of all parameters used in computation

 #' @param uncertainty uncertainty matrix (same supported types as data)

 #' @return list

-#' @importFrom BiocParallel bplapply SnowParam

-distributedCogaps <- function(data, allParams, uncertainty)

+#' @importFrom BiocParallel bplapply MulticoreParam

+distributedCogaps <- function(data, allParams, uncertainty, BPPARAM=NULL)

 {

     FUN <- function(set, data, allParams, uncertainty, fixedMatrix=NULL)

     {

@@ -21,10 +21,11 @@ distributedCogaps <- function(data, allParams, uncertainty)

     # randomly sample either rows or columns into subsets to break the data up

     set.seed(allParams$gaps@seed)

     sets <- createSets(data, allParams)

-    snow <- SnowParam(workers=length(sets), type="SOCK")

-

+    if (is.null(BPPARAM))

+        BPPARAM <- BiocParallel::MulticoreParam(workers=length(sets))

+    

     # run Cogaps normally on each subset of the data

-    initialResult <- bplapply(sets, FUN, BPPARAM=snow, data=data,

+    initialResult <- bplapply(sets, FUN, BPPARAM=BPPARAM, data=data,

         allParams=allParams, uncertainty=uncertainty)

     # match patterns in either A or P matrix

@@ -36,7 +37,7 @@ distributedCogaps <- function(data, allParams, uncertainty)

         == "genome-wide", "P", "A")

     # ru final phase with fixed matrix

-    finalResult <- bplapply(sets, FUN, data=data, BPPARAM=snow,

+    finalResult <- bplapply(sets, FUN, BPPARAM=BPPARAM, data=data, 

         allParams=allParams, uncertainty=uncertainty, fixedMatrix=consensusMatrix)

     # get result 

@@ -27,26 +27,8 @@ NULL

 #' Simulated data

 #' @docType data

-#' @name SimpSim.D

-#' @usage SimpSim.D

-NULL

-

-#' Simulated data

-#' @docType data

-#' @name SimpSim.S

-#' @usage SimpSim.S

-NULL

-

-#' Simulated data

-#' @docType data

-#' @name SimpSim.A

-#' @usage SimpSim.A

-NULL

-

-#' Simulated data

-#' @docType data

-#' @name SimpSim.P

-#' @usage SimpSim.P

+#' @name SimpSim.data

+#' @usage SimpSim.data

 NULL

 #' CoGAPS result from running on simulated data

@@ -57,20 +39,14 @@ NULL

 #' Sample GIST gene expression data from Ochs et al. (2009)

 #' @docType data

-#' @name GIST.D

-#' @usage GIST.D

+#' @name GIST.data_frame

+#' @usage GIST.data_frame

 NULL

 #' Sample GIST gene expression data from Ochs et al. (2009)

 #' @docType data

-#' @name GIST.S

-#' @usage GIST.S

-NULL

-

-#' Simulated dataset to quantify gene set membership in the GIST dataset 

-#' @docType data

-#' @name GIST.GeneSets

-#' @usage GIST.GeneSets

+#' @name GIST.matrix

+#' @usage GIST.matrix

 NULL

 #' CoGAPS result from running on GIST dataset

@@ -79,12 +55,6 @@ NULL

 #' @usage GIST.result

 NULL

-#' Simulated dataset to quantify gene set membership.

-#' @docType data

-#' @name GSets

-#' @usage GSets

-NULL

-

 #' Gene sets defined by transcription factors defined from TRANSFAC.

 #' @docType data

 #' @name tf2ugFC

@@ -79,6 +79,8 @@ setValidity("CogapsParams",

             "random seed must be an integer greater than zero"

         if (object@minNS <= 1 | object@minNS %% 1 != 0)

             "minNS must be an integer greater than one"

+        if (object@nSets <= 1 | object@nSets %% 1 != 0)

+            "minNS must be an integer greater than one"

     }

 )

@@ -104,14 +106,16 @@ setGeneric("setParam", function(object, whichParam, value)

 #'

 #' @description these parameters  are interrelated so they must be set together

 #' @param object an object of type CogapsParams

-#' @param cut a distributed CoGAPS parameter 

-#' @param minNS a distributed CoGAPS parameter 

-#' @param maxNS a distributed CoGAPS parameter 

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

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

+#' pattern matching

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

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

 #' @return the modified params object

 #' @examples

 #'  params <- new("CogapsParams")

 #'  params <- setDistributedParams(3, 2, 4)

-setGeneric("setDistributedParams", function(object, cut=NULL,

+setGeneric("setDistributedParams", function(object, nSets, cut=NULL,

 minNS=NULL, maxNS=NULL)

     {standardGeneric("setDistributedParams")})

@@ -57,7 +57,7 @@ function(object, whichParam, value)

         object@maxGibbsMassA <- value

         object@maxGibbsMassP <- value

     }

-    else if (whichParam %in% c("cut", "minNS", "maxNS"))

+    else if (whichParam %in% c("nSets", "cut", "minNS", "maxNS"))

     {

         stop("please set this parameter with setDistributedParams")

     }

@@ -73,19 +73,15 @@ function(object, whichParam, value)

 #' @aliases setDistributedParams

 #' @importFrom methods slot

 setMethod("setDistributedParams", signature(object="CogapsParams"),

-function(object, cut, minNS, maxNS)

+function(object, nSets, cut, minNS, maxNS)

 {

-    object@cut <- cut

-    object@minNS <- minNS

-    object@maxNS <- maxNS

+    object@nSets <- nSets

-    if (is.null(object@cut))

-        object@cut <- object@nPatterns

-    if (is.null(object@minNS))

-        object@minNS <- ceiling(object@nSets / 2)

-    if (is.null(object@maxNS))

-        object@maxNS <- object@minNS + object@nSets

+    object@cut <- ifelse(is.null(cut), object@nPatterns, cut)

+    object@minNS <- ifelse(is.null(minNS), ceiling(object@nSets / 2), minNS)

+    object@maxNS <- ifelse(is.null(maxNS), object@minNS + object@nSets, maxNS)

+    validObject(object)

     return(object)

 })

@@ -15,18 +15,18 @@ function(object)

 #' @importFrom grDevices rainbow

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

 {

-    nSamples <- nrow(object@sampleFactors)

-    nFactors <- ncol(object@sampleFactors)

+    nSamples <- nrow(x@sampleFactors)

+    nFactors <- ncol(x@sampleFactors)

     colors <- rainbow(nFactors)

     xlimits <- c(0, nSamples + 1)

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

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

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

     for (i in 1:nFactors)

     {

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

-        points(x=1:nSamples, y=object@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,

Binary files a/data/GIST.RData and b/data/GIST.RData differ

Binary files a/data/SimpSim.RData and b/data/SimpSim.RData differ

deleted file mode 100644

@@ -1,12 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/Package.R

-\docType{data}

-\name{GIST.GeneSets}

-\alias{GIST.GeneSets}

-\title{Simulated dataset to quantify gene set membership in the GIST dataset}

-\usage{

-GIST.GeneSets

-}

-\description{

-Simulated dataset to quantify gene set membership in the GIST dataset

-}

similarity index 79%

rename from man/GIST.S.Rd

rename to man/GIST.data_frame.Rd

@@ -1,11 +1,11 @@

 % Generated by roxygen2: do not edit by hand

 % Please edit documentation in R/Package.R

 \docType{data}

-\name{GIST.S}

-\alias{GIST.S}

+\name{GIST.data_frame}

+\alias{GIST.data_frame}

 \title{Sample GIST gene expression data from Ochs et al. (2009)}

 \usage{

-GIST.S

+GIST.data_frame

 }

 \description{

 Sample GIST gene expression data from Ochs et al. (2009)

similarity index 83%

rename from man/GIST.D.Rd

rename to man/GIST.matrix.Rd

@@ -1,11 +1,11 @@

 % Generated by roxygen2: do not edit by hand

 % Please edit documentation in R/Package.R

 \docType{data}

-\name{GIST.D}

-\alias{GIST.D}

+\name{GIST.matrix}

+\alias{GIST.matrix}

 \title{Sample GIST gene expression data from Ochs et al. (2009)}

 \usage{

-GIST.D

+GIST.matrix

 }

 \description{

 Sample GIST gene expression data from Ochs et al. (2009)

deleted file mode 100644

@@ -1,12 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/Package.R

-\docType{data}

-\name{GSets}

-\alias{GSets}

-\title{Simulated dataset to quantify gene set membership.}

-\usage{

-GSets

-}

-\description{

-Simulated dataset to quantify gene set membership.

-}

deleted file mode 100644

@@ -1,12 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/Package.R

-\docType{data}

-\name{SimpSim.A}

-\alias{SimpSim.A}

-\title{Simulated data}

-\usage{

-SimpSim.A

-}

-\description{

-Simulated data

-}

deleted file mode 100644

@@ -1,12 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/Package.R

-\docType{data}

-\name{SimpSim.P}

-\alias{SimpSim.P}

-\title{Simulated data}

-\usage{

-SimpSim.P

-}

-\description{

-Simulated data

-}

deleted file mode 100644

@@ -1,12 +0,0 @@

-% Generated by roxygen2: do not edit by hand

-% Please edit documentation in R/Package.R

-\docType{data}

-\name{SimpSim.S}

-\alias{SimpSim.S}

-\title{Simulated data}

-\usage{

-SimpSim.S

-}

-\description{

-Simulated data

-}

similarity index 75%

rename from man/SimpSim.D.Rd

rename to man/SimpSim.data.Rd

@@ -1,11 +1,11 @@

 % Generated by roxygen2: do not edit by hand

 % Please edit documentation in R/Package.R

 \docType{data}

-\name{SimpSim.D}

-\alias{SimpSim.D}

+\name{SimpSim.data}

+\alias{SimpSim.data}

 \title{Simulated data}

 \usage{

-SimpSim.D

+SimpSim.data

 }

 \description{

 Simulated data

@@ -4,7 +4,7 @@

 \alias{distributedCogaps}

 \title{CoGAPS Distributed Matrix Factorization Algorithm}

 \usage{

-distributedCogaps(data, allParams, uncertainty)

+distributedCogaps(data, allParams, uncertainty, BPPARAM = NULL)

 }

 \arguments{

 \item{data}{File name or R object (see details for supported types)}

@@ -6,19 +6,23 @@

 \alias{setDistributedParams,CogapsParams-method}

 \title{set the value of parameters for distributed CoGAPS}

 \usage{

-setDistributedParams(object, cut = NULL, minNS = NULL, maxNS = NULL)

+setDistributedParams(object, nSets, cut = NULL, minNS = NULL,

+  maxNS = NULL)

-\S4method{setDistributedParams}{CogapsParams}(object, cut = NULL,

+\S4method{setDistributedParams}{CogapsParams}(object, nSets, cut = NULL,

   minNS = NULL, maxNS = NULL)

 }

 \arguments{

 \item{object}{an object of type CogapsParams}

-\item{cut}{a distributed CoGAPS parameter}

+\item{nSets}{number of sets to break data into}

-\item{minNS}{a distributed CoGAPS parameter}

+\item{cut}{number of branches at which to cut dendrogram used in

+pattern matching}

-\item{maxNS}{a distributed CoGAPS parameter}

+\item{minNS}{minimum of individual set contributions a cluster must contain}

+

+\item{maxNS}{maximum of individual set contributions a cluster can contain}

 }

 \value{

 the modified params object

@@ -127,7 +127,6 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix)

         {

             GAPS_ASSERT(!Rf_isNull(allParams["whichMatrixFixed"]));

             std::string which = Rcpp::as<std::string>(allParams["whichMatrixFixed"]);

-            gaps_printf("%s %c\n", which.c_str(), which[0]);

             runner.setFixedMatrix(which[0], convertRMatrix(Rcpp::NumericMatrix(fixedMatrix), which[0]=='P'));

         }

@@ -174,22 +174,63 @@ void GapsRunner::updateSampler(unsigned nA, unsigned nP)

     }

 }

+// sum coef * log(i) for i = 1 to total, fit coef from number of atoms

+// approximates sum of number of atoms (stirling approx to factorial)

+// this should be proportional to total number of updates

+static double estimatedNumUpdates(double current, double total, float nAtoms)

+{

+    double coef = nAtoms / std::log(current);

+    return coef * std::log(std::sqrt(2 * total * gaps::algo::pi)) +

+        total * coef * std::log(total) - total * coef;

+}

+

+double GapsRunner::estimatedPercentComplete() const

+{

+    double nIter = static_cast<double>(mCurrentIteration);

+    double nAtomsA = static_cast<double>(mASampler.nAtoms());

+    double nAtomsP = static_cast<double>(mPSampler.nAtoms());

+    

+    if (mPhase == 'S')

+    {

+        nIter += mMaxIterations;

+    }

+

+    double totalIter = 2.0 * static_cast<double>(mMaxIterations);

+

+    double estimatedCompleted = estimatedNumUpdates(nIter, nIter, nAtomsA) + 

+        estimatedNumUpdates(nIter, nIter, nAtomsP);

+

+    double estimatedTotal = estimatedNumUpdates(nIter, totalIter, nAtomsA) + 

+        estimatedNumUpdates(nIter, totalIter, nAtomsP);

+

+    return estimatedCompleted / estimatedTotal;

+}

+

 void GapsRunner::displayStatus()

 {

     if (mPrintMessages && mOutputFrequency > 0 && ((mCurrentIteration + 1) % mOutputFrequency) == 0)

     {

         bpt::time_duration diff = bpt_now() - mStartTime;

-        unsigned elapsedSeconds = static_cast<unsigned>(diff.total_seconds());

+        double nSecondsCurrent = diff.total_seconds();

+        double nSecondsTotal = nSecondsCurrent / estimatedPercentComplete();

+

+        unsigned elapsedSeconds = static_cast<unsigned>(nSecondsCurrent);

+        unsigned totalSeconds = static_cast<unsigned>(nSecondsTotal);

+

+        unsigned elapsedHours = elapsedSeconds / 3600;

+        elapsedSeconds -= elapsedHours * 3600;

+        unsigned elapsedMinutes = elapsedSeconds / 60;

+        elapsedSeconds -= elapsedMinutes * 60;

-        unsigned hours = elapsedSeconds / 3600;

-        elapsedSeconds -= hours * 3600;

-        unsigned minutes = elapsedSeconds / 60;

-        elapsedSeconds -= minutes * 60;

-        unsigned seconds = elapsedSeconds;

+        unsigned totalHours = totalSeconds / 3600;

+        totalSeconds -= totalHours * 3600;

+        unsigned totalMinutes = totalSeconds / 60;

+        totalSeconds -= totalMinutes * 60;

-        gaps_printf("%d of %d, Atoms: %lu(%lu), ChiSq: %.0f, elapsed time: %02d:%02d:%02d\n",

+        gaps_printf("%d of %d, Atoms: %lu(%lu), ChiSq: %.0f, Time: %02d:%02d:%02d / %02d:%02d:%02d\n",

             mCurrentIteration + 1, mMaxIterations, mASampler.nAtoms(),

-            mPSampler.nAtoms(), mASampler.chi2(), hours, minutes, seconds);

+            mPSampler.nAtoms(), mASampler.chi2(), elapsedHours, elapsedMinutes,

+            elapsedSeconds, totalHours, totalMinutes, totalSeconds);

         gaps_flush();

     }

 }

@@ -61,6 +61,7 @@ private:

     void runOnePhase();

     void updateSampler(unsigned nA, unsigned nP);

+    double estimatedPercentComplete() const;

     void displayStatus();

     void createCheckpoint();

@@ -156,9 +156,9 @@ inline std::string buildReport()

 #endif

 #if defined( __GAPS_AVX__ )

-    std::string simd = "AVX enabled\n";

+    std::string simd = "SIMD: AVX instructions enabled\n";

 #elif defined( __GAPS_SSE__ )

-    std::string simd = "SSE enabled\n";

+    std::string simd = "SIMD: SSE instructions enabled\n";

 #else

     std::string simd = "SIMD not enabled\n";

 #endif

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

 ---

 title: "CoGAPS - Coordinated Gene Association in Pattern Sets"

-author: "Genevieve Stein-O'Brien, Thomas Sherman, Hyejune Limb, Elana Fertig"

+author: "Thomas Sherman, Genevieve Stein-O'Brien, Hyejune Limb, Elana Fertig"

 date: "`r doc_date()`"

 bibliography: References.bib

 vignette: >

@@ -17,15 +17,23 @@ library(CoGAPS)

 # Introduction

-Coordinated Gene Association in Pattern Sets (CoGAPS) infers underlying patterns

-in a matrix of measurements that can be interpreted as arising from the

-multiplication of two lower dimensional matrices. 

-

-The Markov chain Monte Carlo (MCMC) matrix factorization that infers patterns

-also infers the extent to which individual genes belong to these patterns. The

-CoGAPS algorithm extends GAPS to infer the coordinated activity in sets of genes

-for each of the inferred patterns based upon (5) and to refine gene set

-membership based upon (2).

+Coordinated Gene Association in Pattern Sets (CoGAPS) is a technique for latent

+space learning in gene expression data. CoGAPS is a member of the

+Nonnegative Matrix Factorization (NMF) class of algorithms. NMFs factorize a

+data matrix into two related matrices containing gene weights, the

+Amplitude (A) matrix, and sample weights, the Pattern (P) Matrix. Each column

+of A or row of P defines a feature and together this set of features defines

+the latent space among genes and samples, respectively. In NMF, the values of

+the elements in the A and P matrices are constrained to be greater than or

+equal to zero. This constraint simultaneously reflects the non-negative nature

+of gene expression data and enforces the additive nature of the resulting

+feature dimensions, generating solutions that are biologically intuitive to

+interpret @SEUNG_1999.

+

+CoGAPS has two extensions that allow it to scale up to large data sets,

+Genome-Wide CoGAPS (GWCoGAPS) and Single-Cell CoGAPS (scCOGAPS). This package

+presents a unified R interface for all three methods, with a parallel,

+efficient underlying implementation in C++.

 # Installing CoGAPS

@@ -37,7 +45,7 @@ source("https://bioconductor.org/biocLite.R")

 biocLite("CoGAPS")

 ```

-The most up to date version of *CoGAPS* can be installed directly from the 

+The most up-to-date version of *CoGAPS* can be installed directly from the 

 *FertigLab* Github Repository:

 ```{r eval=FALSE}

@@ -49,8 +57,8 @@ devtools::install_github("FertigLab/CoGAPS")

 ```

 There is also an option to install the development version of *CoGAPS*, 

-although this is not guaranteed to be stable and should not be used for

-for publications.

+while this version has the latest experimental features, it is not guaranteed

+to be stable.

 ```{r eval=FALSE}

 ## Method 1 using biocLite

@@ -60,36 +68,48 @@ biocLite("FertigLab/CoGAPS", ref="develop", dependencies = TRUE, build_vignettes

 devtools::install_github("FertigLab/CoGAPS", ref="develop")

 ```

-# Running CoGAPS

+# Package Overview

-Here we walk through running *CoGAPS* on the *GIST* data set from Ochs et

-al. (2009), which contains gene expression data from gastrointestinal

-stromal tumor cell lines treated with Gleevec. The matrix contains bulk gene

-expression data with 1363 genes and 9 samples.

+We first give a walkthrough of the package features using a simple, simulated

+data set. In later sections we provide two example workflows on real data

+sets.

-## Running with Default Parameters

+## Running CoGAPS with Default Parameters

-The only required argument to *CoGAPS* is the data set. The simplest call to

-*CoGAPS* is the following:

+The only required argument to `CoGAPS` is the data set. This can be a `matrix`,

+`data.frame`, `SummarizedExperiment`, `SingleCellExperiment` or the path of a

+file (`tsv`, `csv`, `mtx`) containing the data.

-```{r eval=FALSE}

+```{r}

 # load data

-data(GIST)

+data(SimpSim)

 # run CoGAPS

-result <- CoGAPS(GIST.D)

+CoGAPS(SimpSim.data)

 ```

+While CoGAPS is running it periodically prints status messages. For example,

+`20000 of 25000, Atoms: 2932(80), ChiSq: 9728, time: 00:00:29 / 00:01:19`. This

+message tells us that CoGAPS is at iteration 20000 out of 25000 for this phase,

+and that 29 seconds out of an estimated 1 minute 19 seconds have passed. It 

+also tells us the size of the atomic domain which is a core component of the 

+algorithm but can be ignored for now. Finally, the ChiSq value tells us how

+closely the A and P matrices reconstruct the original data. In general, we want

+this value to go down - but it is not a perfect measurment of how well CoGAPS

+is finding the biological processes contained in the data. CoGAPS also prints

+a message indicating which phase is currently happening. There are two phases

+to the algorithm - *Calibration* and *Sampling*.

+

 ## Setting Parameters

 ### Model Parameters

-Most of the time we'll want to set some parameters before running *CoGAPS*.

-Parameters are managed with a *CogapsParams* object. This object will

+Most of the time we'll want to set some parameters before running CoGAPS.

+Parameters are managed with a `CogapsParams` object. This object will

 store all parameters needed to run CoGAPS and provides a simple interface for

 viewing and setting the parameter values.

-```{r eval=FALSE}

+```{r}

 # create new parameters object

 params <- new("CogapsParams")

@@ -105,219 +125,178 @@ getParam(params, "nPatterns")

 ``` 

 Once we've created the parameters object we can pass it along with our data to

-*CoGAPS*.

+`CoGAPS`.

-```{r eval=FALSE}

+```{r}

 # run CoGAPS with specified model parameters

-result <- CoGAPS(GIST.D, params)

+CoGAPS(SimpSim.data, params)

 ```

 ### Run Configuration Options

-The *CogapsParams* class manages the model parameters - i.e. the parameters

-that affect the result. There are also a few parameters that are passed directly

-to *CoGAPS* that control things like displaying the status of the run.

+The `CogapsParams` class manages the model parameters - i.e. the parameters

+that affect the result. There are also a few parameters that are passed

+directly to `CoGAPS` that control things like displaying the status of the run.

-```{r eval=FALSE}

+```{r}

 # run CoGAPS with specified output frequency

-result <- CoGAPS(GIST.D, params, outputFrequency=200)

+CoGAPS(SimpSim.data, params, outputFrequency=250)

 ```

-There are several other arguments that are passed directly to *CoGAPS* which

+There are several other arguments that are passed directly to `CoGAPS` which

 are covered in later sections.

-# Visualizing Output

+## Visualizing Output

-*CoGAPS* returns a *CogapsResult* object that can be passed on to the analysis

-and plotting functions provided in the package. By default, the *plot* function

-displays the patterns across samples.

+`CoGAPS` returns a `CogapsResult` object that can be passed on to the analysis

+and plotting functions provided in the package. By default, the `plot` function

+displays how the patterns vary across the samples. (Note that we pass the 

+`nIterations` parameter here directly, this is allowed for any parameters in 

+the `CogapsParams` class and will always take precedent over the values given

+in `params`).

-```{r eval=FALSE}

-# show result summary

-result

+```{r}

+# store result

+result <- CoGAPS(SimpSim.data, params, nIterations=5000, outputFrequency=2500)

 # plot CogapsResult object returned from CoGAPS

 plot(result)

 ```

-## Pattern Markers

-

-# Selecting Appropiate Number of Patterns

+In the example workflows we'll explore some more analysis functions provided in

+the package.

-Selecting the best value for *nPatterns* is the most difficult part of the

-analysis. For starters, there is not one "best" value for the number of

-patterns - various numbers of patterns can capture various levels of

-granularity in the data. To further complicate the problem, there's not a

-clear way to compare runs for different numbers of patterns.

+## Running CoGAPS in Parallel

-Here we show the simplest approach of selecting dimensionality by plotting the

-error and selecting the least number of patterns that sufficiently reduce the

-error. We also introduce another way to pass parameters - any parameter in the 

-*CogapsParams* class can be passed by name directly to the *CoGAPS* function,

-overwriting the value contained in `params`.

+Non-Negative Matrix Factorization algorithms typically require long computation

+times and CoGAPS is no exception. In order to scale CoGAPS up to the size of

+data sets seen in practice we need to take advantage of modern hardware

+and parallelize the algorithm.

-```{r eval=FALSE}

-# define the range of patterns we are searching over

-pattern_range <- c(3,5,8)

+### Multi-Threaded Parallelization

-# run CoGAPS with each value in range

-resultList <- lapply(pattern_range, function(p) CoGAPS(GIST.D, params, nPatterns=p, nIterations=3000, outputFrequency=2500))

+The simplest way to run CoGAPS in parallel is to provide the `nThreads`

+argument to `CoGAPS`. This allows the underlying algorithm to run on multiple

+threads and has no effect on the mathematics of the algorithm i.e. this is

+still standard CoGAPS. The precise number of threads to use depends on many

+things like hardware and data size. The best approach is to play around with

+different values and see how it effects the estimated time.

-# plot chi-sq values for each run

-chisq <- sapply(resultList, function(result) getMeanChiSq(result))

-plot(pattern_range, chisq)

+```{r}

+CoGAPS(SimpSim.data, nIterations=10000, outputFrequency=5000, nThreads=1)

+CoGAPS(SimpSim.data, nIterations=10000, outputFrequency=5000, nThreads=4)

 ```

-# Running CoGAPS in Parallel

-

-## Single Process Parallelization

+Note this method relies on CoGAPS being compiled with OpenMP support, use

+`buildReport` to check.

-*CoGAPS* can be run across multiple cores by setting the *nCores* parameter.

-This is different from (Running Distributed CoGAPS)[running-distributed-cogaps].

-The algorithm will run in the normal fashion, without splitting up the data.

-

-```{r eval=FALSE}

-params <- setParam(params, "nCores", 4)

+```{r}

+cat(CoGAPS::buildReport())

 ```

-## Running Distributed CoGAPS

-

-For extremely large datasets (greater than a few thousand rows or columns), it

-is much more efficient to break up the data into random subsets and perform

-*CoGAPS* on each subset in parallel - stitching the results back together

-at the end. The full explanation of this method can be seen here

- (CITE GWCOGAPS).

-

-To run *CoGAPS* in this distributed way there are two distinct functions.

-Genome-Wide CoGAPS - *GWCoGAPS* and Single Cell CoGAPS - *scCoGAPS*. 

-*GWCoGAPS* is used to break the genes (rows) into subsets and *scCoGAPS* is

-used to break the samples (columns - cells in single cell data) into subsets.

-Both methods follow a similiar workflow and the underlying method is 

-essentially the same.

-

-The first step is breaking the data into subsets. This is done as a separate

-step where all subsets are then saved back to a file. In the case where the

-data set is so large it requires a machine with large memory, this allows the

-user to only pay for the machine for enough time to split the data. The actual 

-algorithm should use a reasonable amount of memory on each machine during

-the distributed computation.

-

-```{r eval=FALSE}

-# splitting data into subsets for GWCoGAPS

-gw_sim_name <- createGWCoGAPSSets(as.matrix(GIST.D), as.matrix(GIST.S), nSets=4, "gw_example")

+### Distributed CoGAPS (GWCoGAPS/scCoGAPS)

+

+For large datasets (greater than a few thousand genes or samples) the

+multi-threaded parallelization isn't enough. It is more efficient to break up

+the data into subsets and perform CoGAPS on each subset in parallel, stitching 

+the results back together at the end. The CoGAPS extensions, GWCOGAPS and 

+scCoGAPS, each implement a version of this method (CITE).

+

+In order to use these extensions, some additional parameters are required.

+`nSets` specifies the number of subsets to break the data set into. `cut`,

+`minNS`, and `maxNS` control the process of matching patterns across subsets

+and in general should not be changed from defaults. More information about

+these parameters are available in the original papers (CITE). These parameters

+need to be set with a different function than `setParam` since they depend

+on each other. Here we only set `nSets` (always required), but we have the

+option to pass the other parameters as well.

+

+```{r}

+params <- setDistributedParams(params, nSets=3)

+```

-# splitting data into subsets for scCoGAPS

-sc_sim_name <- createscCoGAPSSets(t(as.matrix(GIST.D)), nSets=4, "sc_example")

+Setting `nSets` requires balancing available hardware and run time against the

+size of your data. In general, `nSets` should be less than or equal to the

+number of nodes/cores that are available. If that is true, then the more subsets

+you create, the faster CoGAPS will run - however, some robustness can be lost

+when the subsets get too small. The general rule of thumb is to set `nSets`

+so that each subset has between 1000 and 5000 genes or cells. We will see an

+example of this on real data in the next two sections.

+

+Once the distributed parameters have been set we can call CoGAPS either by

+setting the `distributed` parameter or by using the provided wrapper functions.

+The following calls are equivalent:

+

+```{r}

+# genome-wide CoGAPS

+CoGAPS(SimpSim.data, params, distributed="genome-wide")

+GWCoGAPS(SimpSim.data, params)

+

+# single-cell CoGAPS

+CoGAPS(SimpSim.data, params, distributed="single-cell", singleCell=TRUE)

+scCoGAPS(SimpSim.data, params)

 ```

-Notice that we pass the transpose of the *GIST* data to *scCoGAPS* since

-the normal dataset is 1363 x 9 and we want a data set with a large number of

-samples.

+Notice that we also set the parameter `singleCell=TRUE`. This makes some 

+adjustments to the algorithm to account for the sparsity in single-cell data.

+The `scCoGAPS` wrapper automatically sets this parameter for us.

-It is neccesary to name the simulation so that the computation portion knows

-what the subset files are called. For convenience the name provided to the

-subset function is returned back so that it can be saved.

+The parallel backend for this computation is managed by the package `BiocParallel`

+and there is an option for the user to specifiy which backend they want. See the

+[Additional Features](#setting-parallel-backend-for-gwcogapssccogaps)

+section for more information.

-Running the computation is just as easy as running normal *CoGAPS*. The results

-will be in the exact same format as *CoGAPS*, only the intermediate computation

-is different. Notice here that no *params* object is accepted - all parameters

-must be passed by name.

+In general it is preferred to pass a file name to `GWCoGAPS`/`scCoGAPS` since

+otherwise the entire data set must be copied across multiple processes which

+will slow things down and potentially cause an out-of-memory error. We will

+see examples of this in the next two sections.

-```{r eval=FALSE}

-GWCoGAPS(gw_sim_name, nPatterns=3, nIterations=1000)

-scCoGAPS(sc_sim_name, nPatterns=3, nIterations=1000)

-```

+# Workflow Example - Bulk Data

-# Supported Data Formats

+THIS SECTION UNDER CONSTRUCTION

-CoGAPS supports both R objects and file names as input. Any of these data types

-can be passed directly to *CoGAPS* as the *data* parameter.

+Here we walk through running CoGAPS on the *GIST* data set from @OCHS_2009,

+which contains gene expression data from gastrointestinal

+stromal tumor cell lines treated with Gleevec. The matrix contains bulk gene

+expression data with 1363 genes and 9 samples.

-R object

+# Workflow Example - Single Cell Data

-* matrix

-* data.frame

-* SummarizedExperiment

-* SingleCellExperiment

+THIS SECTION UNDER CONSTRUCTION

-File