In steinbaugh/pointillism: Single-cell RNA-seq toolkit

## nolint start
suppressPackageStartupMessages({
    library(goalie)
    library(basejump)
    library(Seurat)
    library(pointillism)
})
prepareTemplate(package = "pointillism")
source("_setup.R")
## knitr arguments for `rmarkdown::render()` looping.
## > opts_chunk[["set"]](
## >     cache.path = paste(
## >         params[["seurat_name"]],
## >         "clustering",
## >         "cache/",
## >         sep = "_"
## >     ),
## >     fig.path = paste(
## >         params[["seurat_name"]],
## >         "clustering",
## >         "files/",
## >         sep = "_"
## >     )
## > )
## nolint end

This workflow is adapted from the following sources:

• Seurat guided clustering tutorial
• Seurat cell-cycle scoring and regression

To identify clusters, the following steps will be performed:

1. Normalization and transformation of the raw gene counts per cell to account for differences in sequencing depth.
2. Identification of high variance genes.
3. Regression of sources of unwanted variation (e.g. number of UMIs per cell, mitochondrial transcript abundance, cell cycle phase).
4. Identification of the primary sources of heterogeneity using principal component (PC) analysis and heatmaps.
5. Clustering cells based on significant PCs (metagenes).

---

Load Seurat object

file <- params[["seurat_file"]]
object <- import(file)
name <- basenameSansExt(file)
rm(file)
assert(
    is(object, "Seurat"),
    is.character(name)
)
invisible(validObject(object))
print(object)

organism <- organism(object)
assert(is.character(organism))
## Quality control features to plot.
features <- c(
    "nCount_RNA",
    "nFeature_RNA",
    "S.Score",
    "G2M.Score"
)
## Vector to use for dimensional reduction plot looping.
reducedDimGroups <- unique(c(
    "ident",
    "sampleName",
    "Phase",
    interestingGroups(object)
))

Initialize Seurat (r seurat_name)

First, let's create a Seurat object using the raw counts from the cells that have passed our quality control filtering parameters. Next, the raw counts are normalized using global-scaling normalization with the NormalizeData() function. This (1) normalizes the gene expression measurements for each cell by the total expression; (2) multiplies this by a scale factor (10,000 by default); and (3) log-transforms the result. Following normalization, the FindVariableFeatures() function is then called, which calculates the average expression and dispersion for each gene, places these genes into bins, and then calculates a z-score for dispersion within each bin. This helps control for the relationship between variability and average expression. Finally, the genes are scaled and centered using the ScaleData() function.

object <- NormalizeData(object)
object <- FindVariableFeatures(object)
object <- ScaleData(object)

Plot variable genes

To better cluster our cells, we need to detect the genes that are most variable within our dataset. We can plot dispersion (a normalized measure of to cell-to-cell variation) as a function of average expression for each gene to identify a set of high-variance genes.

VariableFeaturePlot(object)

Regress out unwanted sources of variation

Your single-cell dataset likely contains "uninteresting" sources of variation. This can include technical noise, batch effects, and/or uncontrolled biological variation (e.g. cell cycle). Regressing these signals out of the analysis can improve downstream dimensionality reduction and clustering [@Buettner2015-ur]. To mitigate the effect of these signals, [Seurat][] constructs linear models to predict gene expression based on user-defined variables.

Cell-cycle scoring

First, we assign each cell a score, based on its expression of G2/M and S phase markers. These marker sets should be anticorrelated in their expression levels, and cells expressing neither are likely not cycling and in G1 phase. We assign scores in the CellCycleScoring() function, which stores S and G2/M scores, along with the predicted classification of each cell in either G2M, S or G1 phase.

data(cellCycleMarkersList, package = "AcidSingleCell")
cellCycleMarkers <- cellCycleMarkersList[[camelCase(organism)]]
assert(is(cellCycleMarkers, "CellCycleMarkers"))
markdownHeader("S phase markers", level = 3L)
sGenes <- as.character(cellCycleMarkers[["s"]][["geneName"]])
assert(hasNoDuplicates(sGenes))
print(sGenes)
markdownHeader("G2/M phase markers", level = 3L)
g2mGenes <- as.character(cellCycleMarkers[["g2m"]][["geneName"]])
assert(hasNoDuplicates(g2mGenes))
print(g2mGenes)

object <- CellCycleScoring(
    object = object,
    g2m.features = g2mGenes,
    s.features = sGenes,
    set.ident = TRUE
)
assignAndSaveData(
    name = paste(params[["seurat_name"]], "preregress", sep = "_"),
    object = object,
    dir = params[["data_dir"]]
)

Here we are checking to see if the cells are grouping by cell cycle. If we don't see clear grouping of the cells into G1, G2M, and S clusters on the PCA plot, then it is recommended that we don't regress out cell-cycle variation. When this is the case, remove S.Score and G2M.Score from the variables to regress (vars_to_regress) in the R Markdown YAML parameters.

cellCycleGenes <- sort(unique(c(sGenes, g2mGenes)))
object <- RunPCA(
    object = object,
    features = cellCycleGenes,
    verbose = FALSE
)
plotPca(
    object = object,
    interestingGroups = "Phase",
    label = FALSE,
    dark = params[["dark"]]
)

Apply regression variables

Here we are regressing out variables of uninteresting variation, using the vars.to.regress argument in the ScaleData() function. When variables are defined in the vars.to.regress argument, [Seurat][] regresses them individually against each gene, then rescales and centers the resulting residuals.

We generally recommend minimizing the effects of variable read count depth (nUMI) and mitochondrial gene expression (mitoRatio) as a standard first-pass approach. If the differences in mitochondrial gene expression represent a biological phenomenon that may help to distinguish cell clusters, then we advise not passing in mitoRatio here.

When regressing out the effects of cell-cycle variation, include S.Score and G2M.Score in the vars.to.regress argument. Cell-cycle regression is generally recommended but should be avoided for samples containing cells undergoing differentiation.

object <- ScaleData(
    object = object,
    vars.to.regress = params[["vars_to_regress"]],
    verbose = FALSE
)

Now that regression has been applied, let's recheck to see if the cells are no longer clustering by cycle. We should now see the phase clusters superimpose.

object <- RunPCA(
    object = object,
    features = cellCycleGenes,
    verbose = FALSE
)
plotPca(
    object = object,
    interestingGroups = "Phase",
    label = FALSE,
    dark = params[["dark"]]
)

Linear dimensionality reduction {.tabset}

Next, we perform principal component analysis (PCA) on the scaled data with RunPCA(). By default, the variable genes are used as input, but can be defined using the features argument. ProjectPCA() scores each gene in the dataset (including genes not included in the PCA) based on their correlation with the calculated components. Though we don't use this further here, it can be used to identify markers that are strongly correlated with cellular heterogeneity, but may not have passed through variable gene selection.

object <- RunPCA(object, verbose = TRUE)

Heatmap

In particular, DimHeatmap() allows for easy exploration of the primary sources of heterogeneity in a dataset, and can be useful when trying to decide which PCs to include for further downstream analyses. Both cells and genes are ordered according to their PCA scores.

invisible(Map(
    f = DimHeatmap,
    dims = seq(from = 1L, to = params[["pc_compute"]], by = 1L),
    MoreArgs = list("object" = object)
))

Determine statistically significant principal components

To overcome the extensive technical noise in any single gene for scRNA-seq data, [Seurat][] clusters cells based on their PCA scores, with each PC essentially representing a "metagene" that combines information across a correlated gene set. Determining how many PCs to include downstream is therefore an important step. To accomplish this, we plot the standard deviation of each PC as an elbow plot with our plotPcElbow() function.

PC selection — identifying the true dimensionality of a dataset — is an important step for [Seurat][], but can be challenging/uncertain. We therefore suggest these three approaches to consider:

1. Supervised, exploring PCs to determine relevant sources of heterogeneity, and could be used in conjunction with GSEA for example.
2. Implement a statistical test based on a random null model. This can be time-consuming for large datasets, and may not return a clear PC cutoff.
3. Heuristic approach, using a metric that can be calculated instantly.

We're using a heuristic approach here, by calculating where the principal components start to elbow. The plots below show where we have defined the principal compoment cutoff used downstream for dimensionality reduction. This is calculated automatically as the larger value of:

1. The point where the principal components only contribute 5% of standard deviation.
2. The point where the principal components cumulatively contribute 90% of the standard deviation.

This methodology is also commonly used for PC covariate analysis on bulk RNA-seq samples.

dims <- params[["dims"]]
if (!is.numeric(dims)) {
    p <- plotPcElbow(object)
    print(p)
    elbow <- attr(p, "elbow")
    dims <- seq_len(elbow)
}

We are using r length(dims) principal components for dimensionality reduction calculations.

Cluster the cells

Seurat applies a graph-based clustering approach, building upon initial strategies in Macosko et al., and further inspired by SNN-Cliq [@Xu2015-je] and PhenoGraph [@Levine2015-hr]. This approach embeds cells in a graph structure - for example a K-nearest neighbor (KNN) graph, with edges drawn between cells with similar feature expression patterns, and then attempt to partition this graph into highly interconnected ‘quasi-cliques’ or ‘communities’.

As in PhenoGraph, Seurat first constructs a KNN graph based on the euclidean distance in PCA space, and refines the edge weights between any two cells based on the shared overlap in their local neighborhoods (Jaccard similarity). This step is performed using the FindNeighbors() function, and takes as input the previously defined dimensionality of the dataset.

To cluster the cells, Seurat next applies modularity optimization techniques such as the Louvain algorithm (default) or SLM [@Blondel2008-rf], to iteratively group cells together, with the goal of optimizing the standard modularity function. The FindClusters() function implements this procedure, and contains a resolution parameter that sets the ‘granularity’ of the downstream clustering, with increased values leading to a greater number of clusters. We find that setting this parameter between 0.4-1.2 typically returns good results for single-cell datasets of around 3K cells. Optimal resolution often increases for larger datasets. The clusters can be found using the Idents() function.

## Note that `resolution` argument supports multiple calculations.
object <- FindNeighbors(
    object = object,
    dims = dims,
    verbose = FALSE
)
object <- FindClusters(
    object = object,
    resolution = params[["resolution_calc"]],
    verbose = FALSE
)

Run non-linear dimensional reduction

Seurat offers several non-linear dimensional reduction techniques, such as UMAP and tSNE, to visualize and explore these datasets. The goal of these algorithms is to learn the underlying manifold of the data in order to place similar cells together in low-dimensional space. Cells within the graph-based clusters determined above should co-localize on these dimension reduction plots. As input to the UMAP and tSNE, we suggest using the same PCs as input to the clustering analysis.

UMAP {.tabset}

[Uniform Manifold Approximation and Projection (UMAP)][UMAP] is a dimension reduction technique that can be used for visualisation similarly to t-SNE, but also for general non-linear dimension reduction. The algorithm is founded on three assumptions about the data:

1. The data is uniformly distributed on Riemannian manifold.
2. The Riemannian metric is locally constant (or can be approximated as such).
3. The manifold is locally connected.

[UMAP][] visualization requires the [Python][] dependency umap-learn. We recommend installing this with [conda][].

object <- RunUMAP(object, dims = dims)

invisible(lapply(
    X = reducedDimGroups,
    FUN = function(group) {
        markdownHeader(group, level = 3L, asis = TRUE, tabset = TRUE)
        lapply(
            X = params[["resolution_calc"]],
            FUN = function(res) {
                resCol <- paste0("RNA_snn_res.", res)
                assert(isSubset(resCol, colnames(object[[]])))
                markdownHeader(resCol, level = 4L, asis = TRUE, tabset = TRUE)
                Idents(object) <- resCol
                show(plotUmap(
                    object = object,
                    interestingGroups = group,
                    dark = params[["dark"]]
                ))
            }
        )
    }
))

tSNE {.tabset}

[Seurat][] continues to use t-distributed stochastic neighbor embedding (t-SNE) as a powerful tool to visualize and explore these datasets. While we no longer advise clustering directly on t-SNE components, cells within the graph-based clusters determined above should co-localize on the t-SNE plot. This is because the t-SNE aims to place cells with similar local neighborhoods in high-dimensional space together in low-dimensional space. As input to the t-SNE, we suggest using the same PCs as input to the clustering analysis.

object <- RunTSNE(object, dims = dims)

invisible(lapply(
    X = reducedDimGroups,
    FUN = function(group) {
        markdownHeader(group, level = 3L, asis = TRUE, tabset = TRUE)
        lapply(
            X = params[["resolution_calc"]],
            FUN = function(res) {
                resCol <- paste0("RNA_snn_res.", res)
                assert(isSubset(resCol, colnames(object[[]])))
                markdownHeader(resCol, level = 4L, asis = TRUE, tabset = TRUE)
                Idents(object) <- resCol
                show(plotTsne(
                    object = object,
                    interestingGroups = group,
                    dark = params[["dark"]]
                ))
            }
        )
    }
))

Pick resolution to use

Pick the resolution to use, if multiple are stashed.

if (!is.null(params[["resolution_use"]])) {
    resCol <- paste0("RNA_snn_res.", params[["resolution_use"]])
    print(resCol)
    Idents(object) <- resCol
}

Cluster quality control

Let's look at the variance in the number of UMI counts (nUMI), gene detection (nGene), and the percentage of mitochondrial gene expression (mitoRatio), to see if there are any obvious cluster artefacts. We can also assess cell cycle batch effects (S.Score, G2M.Score) and any principal component bias toward individual clusters.

## See also `Seurat::FeaturePlot()`.
plotFeature(
    object = object,
    reduction = "UMAP",
    features = features,
    label = FALSE,
    legend = FALSE,
    dark = params[["dark"]]
)

Let's plot the feature specificity of the top principal components.

## Limit to 6 PCs.
featurePcs <- paste0("PC_", head(dims, n = 6L))

plotFeature(
    object = object,
    reduction = "UMAP",
    features = featurePcs,
    label = FALSE,
    legend = FALSE,
    dark = params[["dark"]]
)

Save data

assignAndSaveData(
    name = params[["seurat_name"]],
    object = object,
    dir = params[["data_dir"]]
)

steinbaugh/pointillism documentation built on Oct. 13, 2023, 10:43 p.m.