inst/doc/CoGAPS.R

In CoGAPS: Coordinated Gene Activity in Pattern Sets

## ----include=FALSE, cache=FALSE-----------------------------------------------
library(CoGAPS)
library(BiocParallel)

## -----------------------------------------------------------------------------
packageVersion("CoGAPS")

## ----eval=FALSE---------------------------------------------------------------
#  source("https://bioconductor.org/biocLite.R")
#  biocLite("CoGAPS")

## ----eval=FALSE---------------------------------------------------------------
#  ## Method 1 using biocLite
#  biocLite("FertigLab/CoGAPS", dependencies = TRUE, build_vignettes = TRUE)
#  
#  ## Method 2 using devtools package
#  devtools::install_github("FertigLab/CoGAPS")

## ----eval=FALSE---------------------------------------------------------------
#  ## Method 1 using biocLite
#  biocLite("FertigLab/CoGAPS", ref="develop", dependencies = TRUE, build_vignettes = TRUE)
#  
#  ## Method 2 using devtools package
#  devtools::install_github("FertigLab/CoGAPS", ref="develop")

## -----------------------------------------------------------------------------
# load data
data(GIST)

# run CoGAPS (low number of iterations since this is just an example)
CoGAPS(GIST.matrix, nIterations=1000)

## -----------------------------------------------------------------------------
# create new parameters object
params <- new("CogapsParams")

# view all parameters
params

# get the value for a specific parameter
getParam(params, "nPatterns")

# set the value for a specific parameter
params <- setParam(params, "nPatterns", 3)
getParam(params, "nPatterns")

## -----------------------------------------------------------------------------
# run CoGAPS with specified model parameters
CoGAPS(GIST.matrix, params, nIterations=1000)

## -----------------------------------------------------------------------------
# run CoGAPS with specified output frequency
CoGAPS(GIST.matrix, params, nIterations=1000, outputFrequency=250)

## -----------------------------------------------------------------------------
# run CoGAPS
result <- CoGAPS(GIST.matrix, params, messages=FALSE, nIterations=1000)

# get the mean ChiSq statistic over all samples
getMeanChiSq(result)

# get the version number used to create this result
getVersion(result)

# get the original parameters used to create this result
getOriginalParameters(result)

## -----------------------------------------------------------------------------
as(result, "LinearEmbeddingMatrix")

## -----------------------------------------------------------------------------
# store result
result <- CoGAPS(GIST.matrix, params, nIterations=5000, messages=FALSE)

# plot CogapsResult object returned from CoGAPS
plot(result)

## -----------------------------------------------------------------------------
CoGAPS(GIST.matrix, nIterations=10000, outputFrequency=5000, nThreads=1, seed=5)
CoGAPS(GIST.matrix, nIterations=10000, outputFrequency=5000, nThreads=4, seed=5)

## -----------------------------------------------------------------------------
cat(CoGAPS::buildReport())

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

## -----------------------------------------------------------------------------
# need to use a file with distributed cogaps
GISTCsvPath <- system.file("extdata/GIST.csv", package="CoGAPS")

# genome-wide CoGAPS
GWCoGAPS(GISTCsvPath, params, messages=FALSE, nIterations=1000)

# genome-wide CoGAPS
CoGAPS(GISTCsvPath, params, distributed="genome-wide", messages=FALSE, nIterations=1000)

# single-cell CoGAPS
scCoGAPS(GISTCsvPath, params, messages=FALSE, transposeData=TRUE, nIterations=1000)

# single-cell CoGAPS
CoGAPS(GISTCsvPath, params, distributed="single-cell", messages=FALSE, transposeData=TRUE, nIterations=1000)

## -----------------------------------------------------------------------------
if (CoGAPS::checkpointsEnabled())
{
    # our initial run
    res1 <- CoGAPS(GIST.matrix, params, checkpointInterval=100, checkpointOutFile="vignette_example.out", messages=FALSE)

    # assume the previous run crashes
    res2 <- CoGAPS(GIST.matrix, checkpointInFile="vignette_example.out", messages=FALSE)

    # check that they're equal
    all(res1@featureLoadings == res2@featureLoadings)
    all(res1@sampleFactors == res2@sampleFactors)
}

## -----------------------------------------------------------------------------
# run CoGAPS with custom uncertainty
data(GIST)
result <- CoGAPS(GIST.matrix, params, uncertainty=GIST.uncertainty, messages=FALSE, nIterations=1000)

## -----------------------------------------------------------------------------
# run CoGAPS with serial backend
scCoGAPS(GISTCsvPath, params, BPPARAM=BiocParallel::SerialParam(), messages=FALSE, transposeData=TRUE, nIterations=1000)

## -----------------------------------------------------------------------------
# sampling with weights
anno <- sample(letters[1:5], size=nrow(GIST.matrix), replace=TRUE)
w <- c(1,1,2,2,1)
names(w) <- letters[1:5]
params <- new("CogapsParams")
params <- setAnnotationWeights(params, annotation=anno, weights=w)
result <- GWCoGAPS(GISTCsvPath, params, messages=FALSE, nIterations=1000)

## -----------------------------------------------------------------------------
# running cogaps with given subsets
sets <- list(1:225, 226:450, 451:675, 676:900)
params <- new("CogapsParams")
params <- setDistributedParams(params, nSets=length(sets))
result <- GWCoGAPS(GISTCsvPath, params, explicitSets=sets, messages=FALSE, nIterations=1000)

## -----------------------------------------------------------------------------
# run GWCoGAPS (subset data so the displayed output is small)
params <- new("CogapsParams")
params <- setParam(params, "nPatterns", 3)
params <- setDistributedParams(params, nSets=2)
result <- GWCoGAPS(GISTCsvPath, params, messages=FALSE, nIterations=1000)

# get the unmatched patterns from each subset
getUnmatchedPatterns(result)

# get the clustered patterns from the set of all patterns
getClusteredPatterns(result)

# get the correlation of each pattern to the cluster mean
getCorrelationToMeanPattern(result)

# get the size of the subsets used
sapply(getSubsets(result), length)

## -----------------------------------------------------------------------------
# initial run
result <- GWCoGAPS(GISTCsvPath, messages=FALSE, nIterations=1000)

# custom matching process (just take matrix from first subset as a dummy)
consensusMatrix <- getUnmatchedPatterns(result)[[1]]

# run with our custom matched patterns matrix
params <- CogapsParams()
params <- setFixedPatterns(params, consensusMatrix, 'P')
GWCoGAPS(GISTCsvPath, params, explicitSets=getSubsets(result), nIterations=1000)

## -----------------------------------------------------------------------------
sessionInfo()