In WeiSong-bio/roryk-bcbioSinglecell: Single-Cell RNA-Seq Utilities

# Last modified 2018-06-15

# Highly recommended to run this template on an HPC cluster
# stopifnot(detectHPC())

bcbioSingleCell::prepareSingleCellTemplate()
source("_setup.R")

# Load SingleCellExperiment object
bcb_name <- load(params$bcb_file)
bcb <- get(bcb_name, inherits = FALSE)
stopifnot(is(bcb, "SingleCellExperiment"))
invisible(validObject(bcb))

# Vector to use for dimensional reduction plot looping
dim_red_groups <- unique(c(
    "ident", "sampleName", interestingGroups(bcb), "Phase"
))

# Quality control features to plot
features <- c(
    "nUMI",
    "nGene",
    "log10GenesPerUMI",
    "mitoRatio",
    "S.Score",
    "G2M.Score"
)

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

This workflow is adapted from the following sources:

• Satija Lab: Seurat v2 Guided Clustering Tutorial
• Paul Hoffman: 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).

---

Initialize Seurat (r bcb_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 FindVariableGenes() 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.

# S4 object coercion method using `setAs()`, documented in `setAs.R` file. This
# handles gene to symbol conversion and stashes metadata in the `@misc` slot.
# getMethod("coerce", signature(from = "bcbioSingleCell", to = "seurat"))
seurat <- as(bcb, "seurat") %>%
    NormalizeData(
        normalization.method = "LogNormalize",
        scale.factor = 10000
    ) %>%
    FindVariableGenes(
        mean.function = ExpMean,
        dispersion.function = LogVMR,
        do.plot = FALSE
    ) %>%
    ScaleData(model.use = "linear")

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.

VariableGenePlot(seurat)

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. The scaled z-scored residuals of these models are stored in the [e-mail] slot, and are used for dimensionality reduction and clustering.

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 in [e-mail], along with the predicted classification of each cell in either G2M, S or G1 phase.

organism <- metadata(bcb)$organism
cell_cycle_markers <- bcbioSingleCell::cellCycleMarkers[[camel(organism)]]
stopifnot(is.data.frame(cell_cycle_markers))

markdownHeader("S phase markers", level = 3)
s_genes <- cell_cycle_markers %>%
    filter(phase == "S") %>%
    pull("geneName")
print(s_genes)

markdownHeader("G2/M phase markers", level = 3)
g2m_genes <- cell_cycle_markers %>%
    filter(phase == "G2/M") %>%
    pull("geneName")
print(g2m_genes)

saveData(cell_cycle_markers, s_genes, g2m_genes, dir = params$data_dir)

seurat <- CellCycleScoring(
    object = seurat,
    g2m.genes = g2m_genes,
    s.genes = s_genes
)
# Cell-cycle `Phase` column should now be added to `[e-mail]`
assignAndSaveData(
    name = paste(params$seurat_name, "preregress", sep = "_"),
    object = seurat,
    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.

RunPCA(
    object = seurat,
    pc.genes = c(s_genes, g2m_genes),
    do.print = FALSE
) %>%
    plotPCA(
        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.

print(params$vars_to_regress)
seurat <- ScaleData(seurat, vars.to.regress = params$vars_to_regress)

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.

RunPCA(
    object = seurat,
    pc.genes = c(s_genes, g2m_genes),
    do.print = FALSE
) %>%
    plotPCA(
        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 genes in [e-mail] are used as input, but can be defined using the pc.genes 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. The results of the projected PCA can be explored by setting use.full = TRUE for PrintPCA().

seurat <- seurat %>%
    RunPCA(do.print = FALSE) %>%
    ProjectPCA(do.print = FALSE)

PCHeatmap()

In particular, PCHeatmap() 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. Setting cells.use to a number plots the "extreme" cells on both ends of the spectrum, which dramatically speeds plotting for large datasets. Though clearly a supervised analysis, we find this to be a valuable tool for exploring correlated gene sets.

colors <- plasma(n = 3, begin = 0, end = 1)
invisible(mapply(
    FUN = PCHeatmap,
    pc.use = 1:params$pc_compute,
    MoreArgs = list(
        object = seurat,
        col.use = CustomPalette(
            low = colors[[1]],
            mid = colors[[2]],
            high = colors[[3]]
        ),
        do.balanced = TRUE,
        label.columns = FALSE,
        remove.key = TRUE
    ),
    SIMPLIFY = FALSE
))

VizPCA()

invisible(mapply(
    FUN = VizPCA,
    pcs.use = 1:params$pc_compute,
    MoreArgs = list(
        object = seurat,
        do.balanced = TRUE,
        font.size = 1,
        nCol = 1
    ),
    SIMPLIFY = FALSE
))

PrintPCA()

PrintPCA(seurat, pcs.print = 1:params$pc_compute)

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 (bottom left).
2. The point where the principal components cumulatively contribute 90% of the standard deviation (bottom right).

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

# Allow for user-defined PCs to use, otherwise calculate using a PC elbow plot
dims_use <- params$dims_use
if (!is.numeric(params$dims_use)) {
    dims_use <- plotPCElbow(seurat)
}

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

Cluster the cells

Seurat uses a graph-based clustering approach, inspired by SNN-Cliq [@Xu2015-je] and PhenoGraph [@Levine2015-hr]. This approach embeds cells in a graph structure, by default using a K-nearest neighbor (KNN) graph, with edges drawn between cells with similar gene 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 distance). To cluster the cells, it then applies modularity optimization techniques [@Blondel2008-rf], to iteratively group cells together, with the goal of optimizing the standard modularity function.

The FindClusters() function implements the procedure, and contains a resolution argument 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.6-1.2 typically returns good results for single cell datasets of around 3K cells. Optimal resolution often increases for larger datasets. The clusters are saved in the seurat@ident slot.

Regarding the value of the resolution argument, use a value < 1 if you want to obtain fewer clusters.

seurat <- FindClusters(
    object = seurat,
    dims.use = dims_use,
    resolution = params$resolution_calc,
    save.SNN = TRUE,
    force.recalc = TRUE
)

A useful feature in [Seurat][] v2.0 is the ability to recall the parameters that were used in the latest function calls for commonly used functions. For FindClusters(), the authors provide the function PrintFindClustersParams() to print a nicely formatted formatted summary of the parameters that were chosen.

PrintFindClustersParams(seurat)

Run non-linear dimensional reduction

t-SNE {.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, although computing the t-SNE based on scaled gene expression is also supported using the genes.use argument.

seurat <- RunTSNE(
    object = seurat,
    dims.use = dims_use,
    tsne.method = "Rtsne",
    do.fast = TRUE
)

PrintTSNEParams(seurat)

invisible(lapply(
    X = dim_red_groups,
    FUN = function(group) {
        markdownHeader(group, level = 3, asis = TRUE, tabset = TRUE)
        lapply(
            X = params$resolution_calc,
            FUN = function(res) {
                res_col <- paste("res", res, sep = ".")
                markdownHeader(res_col, level = 4, asis = TRUE, tabset = TRUE)
                seurat <- SetAllIdent(seurat, id = res_col)
                show(plotTSNE(
                    object = seurat,
                    interestingGroups = group,
                    dark = params$dark
                ))
            }
        )
    }
))

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][].

```{bash, eval=FALSE} conda install -c conda-forge umap-learn

```r
seurat <- RunUMAP(seurat)

invisible(lapply(
    X = dim_red_groups,
    FUN = function(group) {
        markdownHeader(group, level = 3, asis = TRUE, tabset = TRUE)
        lapply(
            X = params$resolution_calc,
            FUN = function(res) {
                res_col <- paste("res", res, sep = ".")
                markdownHeader(res_col, level = 4, asis = TRUE, tabset = TRUE)
                seurat <- SetAllIdent(seurat, id = res_col)
                show(plotUMAP(
                    object = seurat,
                    interestingGroups = group,
                    dark = params$dark
                ))
            }
        )
    }
))

PCA {.tabset}

Note that t-SNE/UMAP is not PCA! The measurement of distance in a t-SNE plot is difficult to interpret, and is most helpful for the relationships of close neighbors. To better infer separation distance between the putative clusters, let's visualizing using PCA.

invisible(lapply(
    X = dim_red_groups,
    FUN = function(group) {
        markdownHeader(group, level = 3, asis = TRUE, tabset = TRUE)
        lapply(
            X = params$resolution_calc,
            FUN = function(res) {
                res_col <- paste("res", res, sep = ".")
                markdownHeader(res_col, level = 4, asis = TRUE, tabset = TRUE)
                seurat <- SetAllIdent(seurat, id = res_col)
                show(plotPCA(
                    object = seurat,
                    interestingGroups = group,
                    dark = params$dark
                ))
            }
        )
    }
))

Pick resolution to use

Pick the resolution to use, if multiple are stashed.

if (!is.null(params$resolution_use)) {
    seurat <- SetAllIdent(
        object = seurat,
        id = paste("res", params$resolution_use, sep = ".")
    )
}

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()`
plotFeatureTSNE(
    object = seurat,
    features = features,
    dark = params$dark
)

plotFeatureUMAP(
    object = seurat,
    features = features,
    dark = params$dark
)

Principal compoments

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

# Limit to 9 PCs
dim_features <- paste0("PC", head(dims_use, 9))

plotFeatureTSNE(
    object = seurat,
    features = dim_features,
    dark = params$dark
)

plotFeatureUMAP(
    object = seurat,
    features = dim_features,
    dark = params$dark
)

assignAndSaveData(
    name = params$seurat_name,
    object = seurat,
    dir = params$data_dir
)

WeiSong-bio/roryk-bcbioSinglecell documentation built on July 6, 2019, 12:03 a.m.