CeldaCG_Run.Rmd

---
title: "Celda_CG Run"
date: "`r format(Sys.time(), '%B %d, %Y')`"
output: 
 html_document:
     toc: true
     toc_float: true
     number_sections: true
     theme: cosmo
     code_folding: hide
params:
  sce: sce
  L: L
  K: K
  sampleLabel: sampleLabel
  altExpName: altExpName
  useAssay: useAssay
  initialL: initialL
  maxL: maxL
  initialK: initialK
  maxK: maxK
  minCell: minCell
  minCount: minCount
  maxFeatures: maxFeatures
  sceFile: sceFile
  showSession: showSession
  pdf: pdf
---


# Introduction
[Celda](https://github.com/campbio/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.
```{r run_setup}
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.
```{r feature_select, fig.height = 9, fig.width = 18}
# 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.

```{r module_split, fig.height = 9, fig.width = 18}
# 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.

```{r cell_split}
# 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}
```{r dimreduce_tsne, results = "asis"}
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}
```{r dimreduce_umap, results = "asis"}
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)```.
```{r saving}
if (!is.null(sceFile)) {
  saveRDS(sce, sceFile)
}
```


## Methods Summary
Co-clustering of features into modules and cells into subpopulations was performed with the [celda](https://www.biorxiv.org/content/10.1101/2020.11.16.373274v2) 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](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-16-S13-S8). [UMAPs](https://arxiv.org/abs/1602.00370) were generated using the ```celdaUmap``` function. [tSNEs](https://lvdmaaten.github.io/publications/papers/JMLR_2008.pdf) were generated with the ```celdaTsne``` function.


## Session Information
```{r session_info, eval = showSession}
sessionInfo()
```