In campbio/celda: CEllular Latent Dirichlet Allocation

Introduction

Celda package (Cellular Latent Dirichlet Allocation) performs co-clustering of features into modules and cells into subpopulations using count data generated by single-cell genomic platforms. The results can then used to explore the different combination of modules that define each cell population.

Run celda_CG

Setup

The following code loads required libraries, reads in parameters, checks input variables, and sets up display options for the Celda_CG Run report.

require(singleCellTK)
require(celda)
require(knitr)
require(gridExtra)
require(ggplot2)

sce <- params$sce
L <- params$L
K <- params$K
sampleName <- params$sampleName
altExpName <- params$altExpName
useAssay <- params$useAssay
initialL <- params$initialL
maxL <- params$maxL
initialK <- params$initialK
maxK <- params$maxK
minCell <- params$minCell
minCount <- params$minCount
sceFile <- params$sceFile
maxFeatures <- params$maxFeatures
pdf <- params$pdf
showSession <- params$showSession

# Set up Rmarkdown variables
tab3 <- "### %s {-}
"
tab4 <- "#### %s {-}
"
space <- "

"

dev <- ifelse(isTRUE(pdf), c("png"), c("png", "pdf"))
knitr::opts_chunk$set(
  echo = TRUE,
  cache = FALSE,
  cache.lazy = FALSE,
  # don't do lazy-loading for big objects
  cache.comments = FALSE,
  fig.align = "center",
  fig.keep = "all",
  dev = dev
)

Selection of features

Before clustering, features that do not have at least r minCount counts in at least r minCell cells are excluded. If the number of features is more than r maxFeatures, then Seurat's VST function is used to select the r maxFeatures most variable features. While Celda can handle features with many zero counts, lowering the number of features can reduce computational time.

# Select features with minimum counts across a minimum number of cells
sce <- selectFeatures(
  sce,
  minCount = minCount,
  minCell = minCell,
  altExpName = altExpName,
  useAssay = useAssay
)

# Use Seurat to find highly variable features if the number of rows is still
# more than maxFeatures
varFilter <-
  ifelse(nrow(altExp(sce, altExpName)) > maxFeatures, TRUE, FALSE)
if (varFilter) {
  temp.sce <- sce
  temp.sce <-
    seuratFindHVG(inSCE = temp.sce, useAssay = useAssay)
  o <-
    head(
      order(
        rowData(temp.sce)$seurat_variableFeatures_vst_varianceStandardized,
        decreasing = TRUE
      ),
      n = maxFeatures
    )
  altExp(sce, altExpName) <-
    subsetSCERows(temp.sce,
                  index = o,
                  returnAsAltExp = FALSE)
}

Determining the number of modules (L)

Two stepwise splitting procedures are implemented in celda to help determine the number of modules (L) and the number of cell populations (K), respectively. The recursiveSplitModule function fits different celda models for a range of L values between r initialL to r maxL. The first model is fit with r paste0("L = ", initialL). Then the celda_G model is used to split each module into two new modules and the likelihood is re-calculated. The split that produced the best overall likelihood out of all splits is used for the next model with L+1 modules. Perplexity can be used as a measure of "goodness of fit" for discrete Bayesian models. The perplexity is calculated for each choice of L. The Rate of Perplexity Change (RPC) can be used to determine a better "elbow" in the perplexity plot. The elbow represents a dramatic shift in the amount that the perplexity decreases with each new module and is a good place to start looking at possible choices of L. Note that sometimes going beyond the elbow can be helpful if certain modules should be further split up after visual inspection.

# Run recursiveModuleSplit to identify modules for different L
moduleSplit <-
  recursiveSplitModule(
    sce,
    initialL = initialL,
    maxL = maxL,
    sampleLabel = sampleLabel,
    altExpName = altExpName,
    useAssay = useAssay
  )

# Plot perplexity to help choose L
sep <- ifelse(maxL > 100, 10, 5)
p1 <-
  plotGridSearchPerplexity(moduleSplit, altExpName = altExpName, sep = sep)
p2 <- plotRPC(moduleSplit, altExpName = altExpName,  sep = sep)
grid.arrange(p1, p2,  ncol = 2)

Determining the number of cell populations (K)

The recursiveSplitCell function fits different celda models for a range of K values from r initialK to r maxK. The number of modules is set to L and the module labels from the recursiveSplitModule output are used for initialization. Similarly, the first model is fit with r paste0("K = ", initialK). Then the celda_C model is used to split each cell population into two new cell populations and the likelihood is re-calculated. The split that produced the best overall likelihood out of all splits is used for the next model with K+1 cell populations. Perplexity and RPC are calculated as described in the previous section. The elbow can be used as a good starting point for possible choices of K. Lastly, the final model is selected and the modules and cells are reordered using hierarchical clustering so that more similar modules and cell populations will have more similar values of L and K, respectively. Different choices for K are also visualized in reduced dimensional plots in the next section.

# Select number of gene modules
moduleSplitSelect <-
  subsetCeldaList(moduleSplit, params = list(L = L))
initial.modules <- celdaModules(moduleSplitSelect)

# Split cell modules
cellSplit <-
  recursiveSplitCell(
    sce,
    initialK = initialK,
    maxK = maxK,
    yInit = initial.modules,
    sampleLabel = sampleLabel,
    altExpName = altExpName,
    useAssay = useAssay,
    reorder = FALSE)

# Show results of clustering genes (perplexity plot)
p3 <-
  plotGridSearchPerplexity(cellSplit, altExpName = altExpName, sep = 5) +
  theme(legend.position = "bottom")
p4 <- plotRPC(cellSplit, altExpName = altExpName, sep = 5) +
  theme(legend.position = "bottom")
grid.arrange(p3, p4,  ncol = 2)

# Select the final model and reorder L/K values
sce <- subsetCeldaList(cellSplit, params = list(K = K))
sce <- reorderCelda(sce)

Dimensionality Reduction {.tabset .tabset-fade}

Reduced dimensional 2-D plots created by algorithms such as tSNE and UMAP are useful for visualizing the relationship between cells. Each point on the plot represents a single cell. Cells closer together on the plot have more similar expression profiles across all genes. The following tSNEs and UMAPs are colored by different celda models with different choices of K. This plots can be used to help select a final solution.

tSNE {.tabset .tabset-fade}

sce <- celdaTsne(sce, useAssay = useAssay, altExpName = altExpName)
tsne <- reducedDim(altExp(sce, altExpName), "celda_tSNE")
for (i in seq.int(initialK, maxK)) {
  cat(sprintf(tab4, paste0("K = ", i)))
  sce.temp <- subsetCeldaList(cellSplit, params = list(K = i))
  print(
    plotDimReduceCluster(
      celdaClusters(sce.temp),
      dim1 = tsne[, 1],
      dim2 = tsne[, 2],
      labelClusters = TRUE
    )
  )
  cat(space)
}

UMAP {.tabset .tabset-fade}

sce <- celdaUmap(sce, useAssay = useAssay, altExpName = altExpName)
umap <- reducedDim(altExp(sce, altExpName), "celda_UMAP")
for (i in seq.int(initialK, maxK)) {
  cat(sprintf(tab4, paste0("K = ", i)))
  sce.temp <- subsetCeldaList(cellSplit, params = list(K = i))
  print(
    plotDimReduceCluster(
      celdaClusters(sce.temp),
      dim1 = umap[, 1],
      dim2 = umap[, 2],
      labelClusters = TRUE
    )
  )
  cat(space)
}

Save the final model

The SCE object with the final celda model will be saved to the file: r basename(sceFile).

if (!is.null(sceFile)) {
  saveRDS(sce, sceFile)
}

Methods Summary

Co-clustering of features into modules and cells into subpopulations was performed with the celda package using the celda_CG model. Features that did not have at least r minCount counts in at least r minCell cells were excluded in the analysis. r if(isTRUE(varFilter)) paste0("The top ", maxFeatures, " variable features were selected with the Seurat VST method using the singleCellTK package.") The number of modules was identified using the recursiveSplitModule function with r paste0("initialL = ", initialL) and r paste0("maxL = ", maxL). The number of cell populations was identified using the recursiveSplitCell function with r paste0("initialK = ", initialK) and r paste0("maxK = ", maxK). Final values of r paste0("L = ", L) and r paste0("K = ", K) were chosen for final analysis based on the curve for the Rate of Perplexity Change. UMAPs were generated using the celdaUmap function. tSNEs were generated with the celdaTsne function.

Session Information

sessionInfo()

campbio/celda documentation built on April 5, 2024, 11:47 a.m.