@@ -16,6 +16,7 @@ Maintainer: Elana J. Fertig <ejfertig@jhmi.edu>,

 Depends:

     R (>= 3.5.0)

 Imports:

+    BiocParallel,

     cluster,

     data.table,

     methods,

@@ -37,7 +38,8 @@ biocViews: GeneExpression, Transcription, GeneSetEnrichment,

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

     MultipleComparison, DimensionReduction

 NeedsCompilation: yes

-RoxygenNote: 6.0.1

+RoxygenNote: 6.1.0

+Encoding: UTF-8

 Collate:

     'class-CogapsParams.R'

     'CoGAPS.R'

@@ -12,7 +12,6 @@ export(computeGeneGSProb)

 export(getMeanChiSq)

 export(getParam)

 export(patternMarkers)

-export(patternMatch)

 export(plotPatternMarkers)

 export(plotResiduals)

 export(reconstructGene)

@@ -24,6 +23,8 @@ exportClasses(CogapsResult)

 importClassesFrom(S4Vectors,Annotated)

 importClassesFrom(S4Vectors,character_OR_NULL)

 importClassesFrom(SingleCellExperiment,LinearEmbeddingMatrix)

+importFrom(BiocParallel,SnowParam)

+importFrom(BiocParallel,bplapply)

 importFrom(Rcpp,evalCpp)

 importFrom(SummarizedExperiment,assay)

 importFrom(cluster,agnes)

@@ -53,7 +53,7 @@ supported <- function(file)

 #' resultC <- CoGAPS(GIST.D, params)

 #' @importFrom methods new is

 #' @importFrom SummarizedExperiment assay

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

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

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

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

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

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

         "checkpointInterval"=checkpointInterval,

         "checkpointInFile"=checkpointInFile,

         "transposeData"=transposeData,

-        "whichMatrixFixed"="" # internal parameter

+        "whichMatrixFixed"=NULL # internal parameter

     )

     allParams <- parseExtraParams(allParams, list(...))

@@ -76,13 +76,14 @@ checkpointInFile=NULL, transposeData=FALSE, ...)

         stop("unsupported file extension for data")

     # check uncertainty matrix

-    if (is(data, "character") & !is(uncertainty, "character"))

+    if (is(data, "character") & !is.null(uncertainty) & !is(uncertainty, "character"))

         stop("uncertainty must be same data type as data (file name)")

     if (is(uncertainty, "character") & !supported(uncertainty))

         stop("unsupported file extension for uncertainty")

     if (!is(data, "character") & !is.null(uncertainty) & !is(uncertainty, "matrix"))

         stop("uncertainty must be a matrix unless data is a file path")

-    checkDataMatrix(data, uncertainty, params)

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

+        checkDataMatrix(data, uncertainty, params)

     # convert data to matrix

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

@@ -137,6 +138,7 @@ checkpointOutFile="gaps_checkpoint.out", checkpointInterval=1000,

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

 {

     params@distributed <- "single-cell"

+    params@singleCell <- TRUE

     CoGAPS(data, params, nThreads, messages, outputFrequency, uncertainty,

         checkpointOutFile, checkpointInterval, checkpointInFile, transposeData,

         ...)

@@ -7,69 +7,87 @@

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

 {

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

+    {

+        internal <- ifelse(is(data, "character"), cogaps_cpp_from_file, cogaps_cpp)

+        raw <- internal(data, allParams, uncertainty, set, fixedMatrix)

+        new("CogapsResult", Amean=raw$Amean, Asd=raw$Asd, Pmean=raw$Pmean,

+            Psd=raw$Psd, seed=raw$seed, meanChiSq=raw$meanChiSq)    

+    }

+

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

     # run Cogaps normally on each subset of the data

-    initialResult <- foreach(i=1:allParams$gaps@nSets) %dopar%

-    {

-        cogaps_cpp(data, allParams, uncertainty, sets[[i]])

-    }

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

+        allParams=allParams, uncertainty=uncertainty)

     # match patterns in either A or P matrix

     consensusMatrix <- findConsensusMatrix(initialResult, allParams)

+    allParams$gaps@nPatterns <- ncol(consensusMatrix)

     # set fixed matrix

     allParams$whichMatrixFixed <- ifelse(allParams$gaps@distributed

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

-    # run all subsets with the same fixed matrix

-    finalResult <- foreach(i=1:allParams$gaps@nSets) %dopar%

-    {

-        cogaps_cpp(data, allParams, uncertainty, sets[[i]], consensusMatrix)

-    }

+    # ru final phase with fixed matrix

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

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

     # get result 

-    resultList <- stitchTogether(finalResult, allParams)

-    resultList$seed <- allParams$modelParams@seed

-    resultList$meanChiSq <- sum(sapply(finalResult, function(r) r$meanChiSq))

-    resultList$diagnostics <- list()

-    return(resultList)

+    return(stitchTogether(finalResult, allParams))

 }

+#' get number of rows from supported file name or matrix

+#' @param data either a file name or a matrix

+#' @return number of rows

 #' @importFrom data.table fread

+#' @importFrom tools file_ext

 nrow_helper <- function(data)

 {

-    if (class(data) == "character")

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

     {

-        switch(file_ext(data),

+        return(switch(tools::file_ext(data),

             "csv" = nrow(data.table::fread(data, select=1)),

             "tsv" = nrow(data.table::fread(data, select=1)),

             "mtx" = as.numeric(data.table::fread(data, nrows=1, fill=TRUE)[1,1])

-        )

+        ))

     }

     return(nrow(data))

 }

+#' get number of columns from supported file name or matrix

+#' @param data either a file name or a matrix

+#' @return number of columns

 #' @importFrom data.table fread

+#' @importFrom tools file_ext

 ncol_helper <- function(data)

 {

-    if (class(data) == "character")

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

     {

-        switch(file_ext(data),

+        return(switch(tools::file_ext(data),

             "csv" = ncol(data.table::fread(data, nrows=1)) - 1,

             "tsv" = ncol(data.table::fread(data, nrows=1)) - 1,

             "mtx" = as.numeric(data.table::fread(data, nrows=1, fill=TRUE)[1,2])

-        )

+        ))

     }

     return(ncol(data))

-}        

+}

+#' partition genes/samples into subsets

+#' @description either genes or samples or partitioned depending on the type

+#' of distributed CoGAPS (i.e. genome-wide or single-cell)

+#' @param data either file name or matrix

+#' @param allParams list of all CoGAPS parameters

+#' @return list of sorted subsets of either genes or samples

 createSets <- function(data, allParams)

 {

-    total <- ifelse(allParams$gaps@distributed == "genome-wide",

+    total <- ifelse(xor(allParams$transposeData, allParams$gaps@distributed == "genome-wide"),

         nrow_helper(data), ncol_helper(data))

     setSize <- floor(total / allParams$gaps@nSets)

@@ -78,13 +96,17 @@ createSets <- function(data, allParams)

     for (n in 1:(allParams$gaps@nSets - 1))

     {

         selected <- sample(remaining, setSize, replace=FALSE)

-        sets[[n]] <- selected

+        sets[[n]] <- sort(selected)

         remaining <- setdiff(remaining, selected)

     }

-    sets[[allParams$gaps@nSets]] <- remaining

+    sets[[allParams$gaps@nSets]] <- sort(remaining)

     return(sets)

 }

+#' find the consensus pattern matrix across all subsets

+#' @param result list of CogapsResult object from all runs across subsets

+#' @param allParams list of all CoGAPS parameters

+#' @return matrix of consensus patterns

 findConsensusMatrix <- function(result, allParams)

 {

     if (allParams$gaps@distributed == "genome-wide")

@@ -93,13 +115,11 @@ findConsensusMatrix <- function(result, allParams)

         patterns <- do.call(cbind, lapply(result, function(x) x@featureLoadings))

     comb <- expand.grid(1:allParams$gaps@nSets, 1:allParams$gaps@nPatterns)

-    rownames(patterns) <- paste(comb[,1], comb[,2], sep=".")

+    colnames(patterns) <- paste(comb[,1], comb[,2], sep=".")

     return(patternMatch(patterns, allParams))

 }

 #' Match Patterns Across Multiple Runs

-#' @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

@@ -115,7 +135,7 @@ patternMatch <- function(allPatterns, allParams)

 {

     cc <- corcut(allPatterns, allParams)

-    ### split by maxNS

+    # split by maxNS

     indx <- which(sapply(cc$PatsByClust, function(x) ncol(x) > allParams$gaps@maxNS))

     while (length(indx) > 0)

     { 

@@ -124,22 +144,27 @@ patternMatch <- function(allPatterns, allParams)

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

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

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

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

-

+        if (length(icc$PatsByClust) > 1)

+        {

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

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

+        }

         indx <- which(sapply(cc$PatsByClust, function(x) ncol(x) > allParams$gaps@maxNS))

     }

     # create matrix of mean patterns - weighted by coefficient of determination

-    PatsByCDSWavg <- t(sapply(1:length(cc$PatsByClust), function(z)

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

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

-    return(PatsByCDSWavg)

+    PatsByCDSWavg <- sapply(1:length(cc$PatsByClust), function(z)

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

+    colnames(PatsByCDSWavg) <- lapply(1:length(cc$PatsByClust), function(x) paste("Pattern", x))

     # scale

-    return(apply(PatsByCDSWavg, 1, function(row) row / max(row)))

+    return(apply(PatsByCDSWavg, 2, function(col) col / max(col)))

 }

+#' cluster patterns together

+#' @param allPatterns matrix of all patterns across subsets

+#' @param allParams list of all CoGAPS parameters

+#' @return patterns listed by which cluster they belong to

 #' @importFrom cluster agnes

 #' @importFrom stats cutree as.hclust cor

 corcut <- function(allPatterns, allParams)

@@ -166,27 +191,31 @@ corcut <- function(allPatterns, allParams)

     return(list("RtoMeanPattern"=RtoMeanPattern, "PatsByClust"=PatsByClust))

 }

+#' concatenate final results across subsets

+#' @param result list of CogapsResult object from all runs across subsets

+#' @param allParams list of all CoGAPS parameters

+#' @return list with all CoGAPS output

 stitchTogether <- function(result, allParams)

 {

-    if (allParams$modelParams@distributed == "genome-wide")

+    if (allParams$gaps@distributed == "genome-wide")

     {

-        consensus <- result[[1]]@Pmean

-        return(list(

-            "Amean" = do.call(rbind, lapply(result, function(x) x@featureLoadings)),

-            "Asd"   = do.call(rbind, lapply(result, function(x) x@featureStdDev)),

-            "Pmean" = consensus,

-            "Psd"   = matrix(0, nrow=nrow(consensus), ncol=ncol(consensus))

-        ))

+        consensus <- result[[1]]@sampleFactors

+        Amean <- do.call(rbind, lapply(result, function(x) x@featureLoadings))

+        Asd   <- do.call(rbind, lapply(result, function(x) x@featureStdDev))

+        Pmean <- consensus

+        Psd   <- matrix(0, nrow=nrow(consensus), ncol=ncol(consensus))

     }

     else

     {

-        consensus <- result[[1]]@Amean

-        return(list(

-            "Amean" = consensus,

-            "Asd"   = matrix(0, nrow=nrow(consensus), ncol=ncol(consensus)),

-            "Pmean" = do.call(rbind, lapply(result, function(x) x@sampleFactors)),

-            "Psd"   = do.call(rbind, lapply(result, function(x) x@sampleStdDev))

-        ))

+        consensus <- result[[1]]@featureLoadings

+        Amean <- consensus

+        Asd   <- matrix(0, nrow=nrow(consensus), ncol=ncol(consensus))

+        Pmean <- do.call(rbind, lapply(result, function(x) x@sampleFactors))

+        Psd   <- do.call(rbind, lapply(result, function(x) x@sampleStdDev))

     }

+

+    return(list("Amean"=Amean, "Asd"=Asd, "Pmean"=Pmean, "Psd"=Psd,

+        "seed"=allParams$gaps@seed,

+        "meanChiSq"=sum(sapply(result, function(r) r@metadata$meanChiSq))))

 }

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

 # Generated by using Rcpp::compileAttributes() -> do not edit by hand

 # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

-cogaps_cpp_from_file <- function(data, allParams, uncertainty, indices = NULL, fixedMatrix = NULL) {

+cogaps_cpp_from_file <- function(data, allParams, uncertainty = NULL, indices = NULL, fixedMatrix = NULL) {

     .Call('_CoGAPS_cogaps_cpp_from_file', PACKAGE = 'CoGAPS', data, allParams, uncertainty, indices, fixedMatrix)

 }

-cogaps_cpp <- function(data, allParams, uncertainty, indices = NULL, fixedMatrix = NULL) {

+cogaps_cpp <- function(data, allParams, uncertainty = NULL, indices = NULL, fixedMatrix = NULL) {

     .Call('_CoGAPS_cogaps_cpp', PACKAGE = 'CoGAPS', data, allParams, uncertainty, indices, fixedMatrix)

 }

@@ -5,7 +5,7 @@

 #' @importClassesFrom S4Vectors Annotated

 #' @importClassesFrom SingleCellExperiment LinearEmbeddingMatrix

 setClass("CogapsResult", contains="LinearEmbeddingMatrix", slots=c(

-    sampleStdDev = "ANY",   # Psd transpose

+    sampleStdDev = "ANY",   # Psd

     featureStdDev = "ANY"   # Asd

 ))

@@ -13,15 +13,16 @@ setClass("CogapsResult", contains="LinearEmbeddingMatrix", slots=c(

 #' @return initialized CogapsResult object

 #' @importFrom methods callNextMethod

 setMethod("initialize", "CogapsResult",

-function(.Object, Amean, Pmean, Asd, Psd, seed, meanChiSq, diagnostics, ...)

+function(.Object, Amean, Pmean, Asd, Psd, seed, meanChiSq, diagnostics=list(),

+...)

 {

     .Object@featureLoadings <- Amean

-    .Object@sampleFactors <- t(Pmean)

+    .Object@sampleFactors <- Pmean

     if (!is.null(Asd))

         .Object@featureStdDev <- Asd

     if (!is.null(Psd))

-        .Object@sampleStdDev <- t(Psd)

+        .Object@sampleStdDev <- Psd

     .Object@metadata[["seed"]] <- seed

     .Object@metadata[["meanChiSq"]] <- meanChiSq

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

 \alias{CoGAPS}

 \title{CoGAPS Matrix Factorization Algorithm}

 \usage{

-CoGAPS(data, params = new("CogapsParams"), nThreads = NULL,

+CoGAPS(data, params = new("CogapsParams"), nThreads = 1,

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

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

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

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

 \name{binaryA}

 \alias{binaryA}

 \alias{binaryA,CogapsResult-method}

-\alias{binaryA}

 \title{binary heatmap for standardized feature matrix}

 \usage{

 binaryA(object, threshold = 3)

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

 \name{calcCoGAPSStat}

 \alias{calcCoGAPSStat}

 \alias{calcCoGAPSStat,CogapsResult-method}

-\alias{calcCoGAPSStat}

 \title{calculate gene set statistics}

 \usage{

 calcCoGAPSStat(object, GStoGenes, numPerm = 500)

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

 \name{calcGeneGSStat}

 \alias{calcGeneGSStat}

 \alias{calcGeneGSStat,CogapsResult-method}

-\alias{calcGeneGSStat}

 \title{probability gene belongs in gene set}

 \usage{

 calcGeneGSStat(object, GStoGenes, numPerm, Pw = rep(1,

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

 \name{calcZ}

 \alias{calcZ}

 \alias{calcZ,CogapsResult-method}

-\alias{calcZ}

 \title{compute z-score matrix}

 \usage{

 calcZ(object, which = "feature")

@@ -4,14 +4,14 @@

 \name{computeGeneGSProb}

 \alias{computeGeneGSProb}

 \alias{computeGeneGSProb,CogapsResult-method}

-\alias{computeGeneGSProb}

 \title{compute gene probability}

 \usage{

 computeGeneGSProb(object, GStoGenes, numPerm = 500, Pw = rep(1,

   ncol(object@featureLoadings)), PwNull = FALSE)

-\S4method{computeGeneGSProb}{CogapsResult}(object, GStoGenes, numPerm = 500,

-  Pw = rep(1, ncol(object@featureLoadings)), PwNull = FALSE)

+\S4method{computeGeneGSProb}{CogapsResult}(object, GStoGenes,

+  numPerm = 500, Pw = rep(1, ncol(object@featureLoadings)),

+  PwNull = FALSE)

 }

 \arguments{

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

new file mode 100644

@@ -0,0 +1,19 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/DistributedCogaps.R

+\name{corcut}

+\alias{corcut}

+\title{cluster patterns together}

+\usage{

+corcut(allPatterns, allParams)

+}

+\arguments{

+\item{allPatterns}{matrix of all patterns across subsets}

+

+\item{allParams}{list of all CoGAPS parameters}

+}

+\value{

+patterns listed by which cluster they belong to

+}

+\description{

+cluster patterns together

+}

new file mode 100644

@@ -0,0 +1,20 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/DistributedCogaps.R

+\name{createSets}

+\alias{createSets}

+\title{partition genes/samples into subsets}

+\usage{

+createSets(data, allParams)

+}

+\arguments{

+\item{data}{either file name or matrix}

+

+\item{allParams}{list of all CoGAPS parameters}

+}

+\value{

+list of sorted subsets of either genes or samples

+}

+\description{

+either genes or samples or partitioned depending on the type

+of distributed CoGAPS (i.e. genome-wide or single-cell)

+}

new file mode 100644

@@ -0,0 +1,19 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/DistributedCogaps.R

+\name{findConsensusMatrix}

+\alias{findConsensusMatrix}

+\title{find the consensus pattern matrix across all subsets}

+\usage{

+findConsensusMatrix(result, allParams)

+}

+\arguments{

+\item{result}{list of CogapsResult object from all runs across subsets}

+

+\item{allParams}{list of all CoGAPS parameters}

+}

+\value{

+matrix of consensus patterns

+}

+\description{

+find the consensus pattern matrix across all subsets

+}

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

 \name{getMeanChiSq}

 \alias{getMeanChiSq}

 \alias{getMeanChiSq,CogapsResult-method}

-\alias{getMeanChiSq}

 \title{return chi-sq of final matrices}

 \usage{

 getMeanChiSq(object)

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

 \name{getParam}

 \alias{getParam}

 \alias{getParam,CogapsParams-method}

-\alias{getParam}

 \title{get the value of a parameter}

 \usage{

 getParam(object, whichParam)

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

 \title{Constructor for CogapsResult}

 \usage{

 \S4method{initialize}{CogapsResult}(.Object, Amean, Pmean, Asd, Psd, seed,

-  meanChiSq, diagnostics, ...)

+  meanChiSq, diagnostics = list(), ...)

 }

 \value{

 initialized CogapsResult object

new file mode 100644

@@ -0,0 +1,17 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/DistributedCogaps.R

+\name{ncol_helper}

+\alias{ncol_helper}

+\title{get number of columns from supported file name or matrix}

+\usage{

+ncol_helper(data)

+}

+\arguments{

+\item{data}{either a file name or a matrix}

+}

+\value{

+number of columns

+}

+\description{

+get number of columns from supported file name or matrix

+}

new file mode 100644

@@ -0,0 +1,17 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/DistributedCogaps.R

+\name{nrow_helper}

+\alias{nrow_helper}

+\title{get number of rows from supported file name or matrix}

+\usage{

+nrow_helper(data)

+}

+\arguments{

+\item{data}{either a file name or a matrix}

+}

+\value{

+number of rows

+}

+\description{

+get number of rows from supported file name or matrix

+}

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

 \name{patternMarkers}

 \alias{patternMarkers}

 \alias{patternMarkers,CogapsResult-method}

-\alias{patternMarkers}

 \title{compute pattern markers statistic}

 \usage{

 patternMarkers(object, threshold = "all", lp = NA, full = FALSE)

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

 \name{plotResiduals}

 \alias{plotResiduals}

 \alias{plotResiduals,CogapsResult-method}

-\alias{plotResiduals}

 \title{plot of residuals}

 \usage{

 plotResiduals(object, data, uncertainty = NULL)

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

 \name{reconstructGene}

 \alias{reconstructGene}

 \alias{reconstructGene,CogapsResult-method}

-\alias{reconstructGene}

 \title{reconstruct gene}

 \usage{

 reconstructGene(object, genes = NULL)

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

 \name{setDistributedParams}

 \alias{setDistributedParams}

 \alias{setDistributedParams,CogapsParams-method}

-\alias{setDistributedParams}

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

 \usage{

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

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

 \name{setParam}

 \alias{setParam}

 \alias{setParam,CogapsParams-method}

-\alias{setParam}

 \title{set the value of a parameter}

 \usage{

 setParam(object, whichParam, value)

new file mode 100644

@@ -0,0 +1,19 @@

+% Generated by roxygen2: do not edit by hand

+% Please edit documentation in R/DistributedCogaps.R

+\name{stitchTogether}

+\alias{stitchTogether}

+\title{concatenate final results across subsets}

+\usage{

+stitchTogether(result, allParams)

+}

+\arguments{

+\item{result}{list of CogapsResult object from all runs across subsets}

+

+\item{allParams}{list of all CoGAPS parameters}

+}

+\value{

+list with all CoGAPS output

+}

+\description{

+concatenate final results across subsets

+}

@@ -23,21 +23,7 @@ public:

     Archive(const std::string &path, std::ios_base::openmode flags)

         : mStream(path.c_str(), std::ios::binary | flags)

-    {

-        /*if (flags == ARCHIVE_WRITE)

-        {

-            *this << ARCHIVE_MAGIC_NUM;

-        }

-        else if (flags == ARCHIVE_READ)

-        {

-            uint32_t magicNum = 0;

-            *this >> magicNum;

-            if (magicNum != ARCHIVE_MAGIC_NUM)

-            {

-                Rcpp::Rcout << "warning: invalid checkpoint file" << std::endl;

-            }

-        }*/

-    }

+    {}

     void close() {mStream.close();}

@@ -80,4 +66,4 @@ public:

     friend Archive& operator>>(Archive &ar, boost::random::mt11213b &val) { return readFromArchive(ar, val); }

 };

-#endif

+#endif // __COGAPS_ARCHIVE_H__

@@ -7,28 +7,32 @@

 // these are helper functions for converting matrix/vector types

 // to and from R objects

-static Matrix convertRMatrix(const Rcpp::NumericMatrix &rmat)

+static Matrix convertRMatrix(const Rcpp::NumericMatrix &rmat, bool transpose=false)

 {

-    Matrix mat(rmat.nrow(), rmat.ncol());

-    for (unsigned i = 0; i < mat.nRow(); ++i)

+    unsigned nr = transpose ? rmat.ncol() : rmat.nrow();

+    unsigned nc = transpose ? rmat.nrow() : rmat.ncol();

+    Matrix mat(nr, nc);

+    for (unsigned i = 0; i < nr; ++i)

     {

-        for (unsigned j = 0; j < mat.nCol(); ++j)

+        for (unsigned j = 0; j < nc; ++j)

         {

-            mat(i,j) = rmat(i,j);

+            mat(i,j) = transpose ? rmat(j,i) : rmat(i,j);

         }

     }

     return mat;

 }

 template <class Matrix>

-static Rcpp::NumericMatrix createRMatrix(const Matrix &mat)

+static Rcpp::NumericMatrix createRMatrix(const Matrix &mat, bool transpose=false)

 {

-    Rcpp::NumericMatrix rmat(mat.nRow(), mat.nCol());

-    for (unsigned i = 0; i < mat.nRow(); ++i)

+    unsigned nr = transpose ? mat.nCol() : mat.nRow();

+    unsigned nc = transpose ? mat.nRow() : mat.nCol();

+    Rcpp::NumericMatrix rmat(nr, nc);

+    for (unsigned i = 0; i < nr; ++i)

     {

-        for (unsigned j = 0; j < mat.nCol(); ++j)

+        for (unsigned j = 0; j < nc; ++j)

         {

-            rmat(i,j) = mat(i,j);

+            rmat(i,j) = transpose ? mat(j,i) : mat(i,j);

         }

     }

     return rmat;

@@ -60,6 +64,7 @@ unsigned getNumPatterns(const Rcpp::List &allParams)

         ar >> nPatterns;

         ar.close();

     }

+    return nPatterns;

 }

 std::vector<unsigned> getSubsetIndices(const Rcpp::Nullable<Rcpp::IntegerVector> &indices)

@@ -80,7 +85,7 @@ bool processDistributedParameters(const Rcpp::List &allParams)

         GAPS_ASSERT(d == "genome-wide" || d == "single-cell");

         return d == "genome-wide";

     }

-    return true;

+    return false; // will be ignored anyways

 }

 // this is the main function that creates a GapsRunner and runs CoGAPS

@@ -118,11 +123,12 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix)

     else // no checkpoint, populate from given parameters

     {

         // set fixed matrix

-        if (fixedMatrix.isNotNull());

+        if (fixedMatrix.isNotNull())

         {

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

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

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

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

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

         }

         // set parameters that would be saved in the checkpoint

@@ -147,9 +153,9 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix)

     GapsResult result(runner.run());

     return Rcpp::List::create(

         Rcpp::Named("Amean") = createRMatrix(result.Amean),

-        Rcpp::Named("Pmean") = createRMatrix(result.Pmean),

+        Rcpp::Named("Pmean") = createRMatrix(result.Pmean, true),

         Rcpp::Named("Asd") = createRMatrix(result.Asd),

-        Rcpp::Named("Psd") = createRMatrix(result.Psd),

+        Rcpp::Named("Psd") = createRMatrix(result.Psd, true),

         Rcpp::Named("seed") = runner.getSeed(),

         Rcpp::Named("meanChiSq") = result.meanChiSq,

         Rcpp::Named("diagnostics") = Rcpp::List::create()

@@ -161,7 +167,7 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix)

 // [[Rcpp::export]]

 Rcpp::List cogaps_cpp_from_file(const Rcpp::CharacterVector &data,

 const Rcpp::List &allParams,

-const Rcpp::Nullable<Rcpp::CharacterVector> &uncertainty,

+const Rcpp::Nullable<Rcpp::CharacterVector> &uncertainty=R_NilValue,

 const Rcpp::Nullable<Rcpp::IntegerVector> &indices=R_NilValue,

 const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix=R_NilValue)

 {

@@ -170,14 +176,13 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix=R_NilValue)

     {

         unc = Rcpp::as<std::string>(Rcpp::CharacterVector(uncertainty));

     }

-

     return cogapsRun(Rcpp::as<std::string>(data), allParams, unc, indices, fixedMatrix);

 }

 // [[Rcpp::export]]

 Rcpp::List cogaps_cpp(const Rcpp::NumericMatrix &data,

 const Rcpp::List &allParams,

-const Rcpp::Nullable<Rcpp::NumericMatrix> &uncertainty,

+const Rcpp::Nullable<Rcpp::NumericMatrix> &uncertainty=R_NilValue,

 const Rcpp::Nullable<Rcpp::IntegerVector> &indices=R_NilValue,

 const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix=R_NilValue)

 {

@@ -186,7 +191,6 @@ const Rcpp::Nullable<Rcpp::NumericMatrix> &fixedMatrix=R_NilValue)

     {

         unc = convertRMatrix(Rcpp::NumericMatrix(uncertainty));

     }

-

     return cogapsRun(convertRMatrix(data), allParams, unc, indices, fixedMatrix);

 }

@@ -40,9 +40,12 @@

                 gaps_stop();                                            \

             }                                                           \

         } while(0)

+

+    #define DEBUG_PING gaps_printf("here %s %d\n", __FILE__, __LINE__);

 #else

     #define GAPS_ASSERT(cond) do {} while(0)

     #define GAPS_ASSERT_MSG(cond, msg) do {} while(0)

+    #define DEBUG_PING   

 #endif

 #endif

\ No newline at end of file

@@ -79,7 +79,6 @@ void GapsRunner::setCheckpointInterval(unsigned interval)

 GapsResult GapsRunner::run()

 {

-#if 0

     mStartTime = bpt_now();

     // calculate appropiate number of threads if compiled with openmp

@@ -114,7 +113,6 @@ GapsResult GapsRunner::run()

             runOnePhase();

             break;

     }

-#endif

     GapsResult result(mStatistics);

     result.meanChiSq = mStatistics.meanChiSq(mASampler);

     return result;    

@@ -161,6 +159,7 @@ void GapsRunner::updateSampler(unsigned nA, unsigned nP)

         {

             mPSampler.sync(mASampler);

         }

+        GAPS_ASSERT(mASampler.internallyConsistent());

     }

     if (mFixedMatrix != 'P')

@@ -171,10 +170,8 @@ void GapsRunner::updateSampler(unsigned nA, unsigned nP)

         {

             mASampler.sync(mPSampler);

         }

+        GAPS_ASSERT(mPSampler.internallyConsistent());

     }

-

-    GAPS_ASSERT(mASampler.internallyConsistent());

-    GAPS_ASSERT(mPSampler.internallyConsistent());

 }

 void GapsRunner::displayStatus()

@@ -193,6 +190,7 @@ void GapsRunner::displayStatus()

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

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

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

+        gaps_flush();

     }

 }

@@ -326,34 +326,41 @@ float GibbsSampler<T, MatA, MatB>::getAvgQueue() const

 template <class T, class MatA, class MatB>

 bool GibbsSampler<T, MatA, MatB>::internallyConsistent()

 {

-    Atom a = mDomain.front();

-    float current = a.mass;

-    uint64_t row = impl()->getRow(a.pos);

-    uint64_t col = impl()->getCol(a.pos);

-

-    while (mDomain.hasRight(a))

+    if (mDomain.size() > 0)

     {

-        a = mDomain.right(a);

-        if (row != impl()->getRow(a.pos) || col != impl()->getCol(a.pos))

+        Atom a = mDomain.front();

+        float current = a.mass;

+        uint64_t row = impl()->getRow(a.pos);

+        uint64_t col = impl()->getCol(a.pos);

+

+        while (mDomain.hasRight(a))

         {

-            float matVal = mMatrix(row, col);

-            if (std::abs(current - matVal) > 0.1f)

+            a = mDomain.right(a);

+            if (row != impl()->getRow(a.pos) || col != impl()->getCol(a.pos))

             {

-                gaps_printf("mass difference detected at row %lu, column %lu: %f %f\n",

-                    row, col, current, matVal); 

-                return false;

+                float matVal = mMatrix(row, col);

+                if (std::abs(current - matVal) > 0.1f)

+                {