Seurat Curated Workflow

Irzam Sarfraz

Source: vignettes/articles/seurat_curated_workflow.Rmd

Introduction

singleCellTK integrates functions from Seurat [1][2][3][4] package in an easy to use streamlined workflow using both the shiny user interface as well as the R console. The shiny application contains a separate tab that lets the users run the steps of the workflow in a sequential manner with ability to visualize through interactive plots from within the application. On the R console, the toolkit offers wrapper functions that use the SingleCellExperiment [5] object as the input and the output. All computations from the wrapper functions are stored within this object for further manipulation.

To view detailed instructions on how to use the workflow, please select ‘Interactive Analysis’ for using the workflow in shiny application or ‘Console Analysis’ for using these methods on R console from the tabs below:

Workflow Guide

In this tutorial example, we illustrate all the steps of the curated workflow and focus on the options available to manipulate and customize the steps of the workflow as per user requirements. This tutorial takes the real-world scRNAseq dataset as an example, which consists of 2,700 Peripheral Blood Mononuclear Cells (PBMCs) collected from a healthy donor, namingly PBMC3K. This dataset is available from 10X Genomics and can be found on the 10X website. To initiate the Seurat workflow, click on the ‘Curated Workflows’ from the top menu and select Seurat:


NOTE: Before heading into the next steps, we assume that users have already loaded SCTK, imported and QC’ed the PBMC3K data, following the Import and QC Tutorial.

1. Normalize Data
Assuming that the data has been uploaded via the Upload tab of the toolkit, the first step for the analysis of the data is the Normalization of data. For this purpose, any assay available in the uploaded data can be used against one of the three methods of normalization available through Seurat i.e. LogNormalize, CLR (Centered Log Ratio) or RC (Relative Counts).


  1. Open up the ‘Normalize’ tab by clicking on it.
  2. Select the assay to normalize from the dropdown menu.
  3. Select the normalization method from the dropdown menu. Available methods are LogNormalize, CLR or RC.
  4. Set the scaling factor which represents the numeric value by which the data values are multiplied before log transformation. Default is set to 10000.
  5. Press the ‘Normalize’ button to start the normalization process.

2. Highly Variable Genes
Identification of the highly variable genes is core to the Seurat workflow and these highly variable genes are used throughout the remaining workflow. Seurat provides three methods for variable genes identification i.e. vst (uses local polynomial regression to fit a relationship between log of variance and log of mean), mean.var.plot (uses mean and dispersion to divide features into bins) and dispersion (uses highest dispersion values only).


  1. Open up the ‘Highly Variable Genes’ tab.
  2. Select method for computation of highly variable genes from vst, mean.var.plot and dispersion.
  3. Input the number of genes that should be identified as variable. Default is 2000.
  4. Press ‘Find HVG’ button to compute the variable genes.
  5. Once genes are computed, select number of the top most variable genes to display in (6).
  6. Displays the top most highly variable genes as selected in (5).
  7. Graph that plots each gene and its relationship based upon the selected model in (2), and identifies if a gene is highly variable or not.

3. Dimensionality Reduction
Seurat workflow offers PCA or ICA for dimensionality reduction and the components from these methods can be used in the downstream analysis. Moreover, several plots are available for the user to inspect the output of the dimensionality reduction such as the standard ‘PCA Plot,’ ‘Elbow Plot,’ ‘Jackstraw Plot’ and ‘Heatmap Plot.’



  1. Open up the ‘Dimensionality Reduction’ tab.
  2. Select a sub-tab for either ‘PCA’ or ‘ICA’ computation. Separate tabs are available for both methods if the user wishes to compute and inspect both separately.
  3. Input the number of components to compute. Default value is 50.
  4. Select the plots that should be computed with the overall processing. The standard ‘PCA Plot’ will be computed at all times, while the remaining can be turned off if not required.
  5. Input the number of features against which a ‘Heatmap’ should be plotted. Only available when ‘Compute Heatmap’ is set to TRUE.
  6. Specify seed value.
  7. Press the ‘Run PCA’ button to start processing.
  8. Select the number of computed components that should be used in the downstream analysis e.g. in ‘tSNE/UMAP’ computation or with ‘Clustering.’ If ‘Elbow Plot’ is computed, a suggested value will be indicated that should be preferred for downstream analysis.
  9. The plot area from where all computed plots can be viewed by the user.
  10. Heatmap plot has various options available for the users to customize the plot. Since a plot is computed against each component, a user-defined number of components can be selected against. Moreover, for viewing quality, a number of columns can be selected in which the plots should be shown.

4. tSNE/UMAP
‘tSNE’ and ‘UMAP’ can be computed and plotted once components are available from ‘Dimensionality Reduction’ tab.


  1. Open up the ‘tSNE/UMAP’ tab.
  2. Select ‘tSNE’ or ‘UMAP’ sub-tab.
  3. Select a reduction method. Only methods that are computed previously in the ‘Dimensionality Reduction’ tab are available.
  4. Set perplexity tuning parameter.
  5. Set seed value.
  6. Information displayed to the user that how many components from the selected reduction method will be used. This value can only be changed from the ‘Dimensionality Reduction’ tab.
  7. Press ‘Run tSNE’ or ‘Run UMAP’ button to start processing.
  8. ‘tSNE’ or ‘UMAP’ plot depending upon the selected computation.

5. Clustering
Cluster labels can be generated for all cells/samples using one of the computed reduction method. Plots are automatically re-computed with cluster labels. The available algorithms for clustering as provided by Seurat include original Louvain algorithm, Louvain algorithm with multilevel refinement and SLM algorithm.


  1. Open up the ‘Clustering’ tab.
  2. Select a previously computed reduction method.
  3. Select clustering algorithm from original Louvain algorithm, Louvain algorithm with multilevel refinement and SLM algorithm
  4. Set resolution parameter value for the algorithm. Default is 0.8.
  5. Set if singletons should be grouped to nearest clusters or not. Default is TRUE.
  6. Information displayed to the user that how many components from the selected reduction method will be used. This value can only be changed from the ‘Dimensionality Reduction’ tab.
  7. Press ‘Find Clusters’ button to start processing.
  8. Re-computed plots with cluster labels. Only those plots are available that have previously been computed.

6. Find Markers
‘Find Markers’ tab can be used to identify and visualize the marker genes using on of the provided visualization methods. The tab offers identification of markers between two selected phenotype groups or between all groups and can be decided at the time of the computation. Furthermore, markers that are conserved between two phenotype groups can also be identified. Visualizations such as Ridge Plot, Violin Plot, Feature Plot and Heatmap Plot can be used to visualize the individual marker genes.


1. Select if you want to identify marker genes against all groups in a biological variable or between two pre-defined groups. Additionally, users can select the last option to identify the marker genes that are conserved between two groups. 2. Select phenotype variable that contains the grouping information. 3. Select test used for marker genes identification. 4. Select if only positive markers should be returned. 5. Press “Find Markers” button to run marker identification. 6. Identified marker genes are populated in the table. 7. Filters can be applied on the table.


8. Filters allow different comparisons based on the type of the column of the table.


9. Table re-populated after applying filters.


10. Heatmap plot can be visualized for all genes populated in the table (9) against all biological groups in the selected phenotype variable.


11. To visualize each individual marker gene through gene plots, they can be selected by clicking on the relevant rows of the table.


12. Selected marker genes from the table are plotted with gene plots.

7. Downstream Analysis
Once all steps of the Seurat workflow are completed, users can further analyze the data by directly going to the various downstream analysis options (Differential Expression, Marker Selection & Pathway Analysis) from within the Seurat workflow.


All methods provided by SCTK for Seurat workflow use a SingleCellExperiment object both as an input and output.

Using a sample dataset:

library(singleCellTK)

# Load filtered pbmc3k data
sce <- readRDS("tutorial_pbmc3k_qc.rds")

print(sce)
## class: SingleCellExperiment 
## dim: 32738 2627 
## metadata(2): sctk assayType
## assays(2): counts decontXcounts
## rownames(32738): MIR1302-10 FAM138A ... AC002321.2 AC002321.1
## rowData names(3): ENSEMBL_ID Symbol_TENx Symbol
## colnames(2627): pbmc3k_AAACATACAACCAC-1 pbmc3k_AAACATTGAGCTAC-1 ...
##   pbmc3k_TTTGCATGAGAGGC-1 pbmc3k_TTTGCATGCCTCAC-1
## colData names(29): Sample Barcode ... decontX_contamination
##   decontX_clusters
## reducedDimNames(2): decontX_UMAP QC_UMAP
## mainExpName: NULL
## altExpNames(0):

1. Normalize Data
Once raw data is uploaded and stored in a SingleCellExperiment object, runSeuratNormalizeData() function can be used to normalize the data. The method returns a SingleCellExperiment object with normalized data stored as a new assay in the input object.

Parameters to this function include useAssay (specify the assay that should be normalized), normAssayName (specify the new name of the normalized assay, defaults to "seuratNormData"), normalizationMethod (specify the normalization method to use, defaults to "LogNormalize") and scaleFactor (defaults to 10000).

sce <- runSeuratNormalizeData(inSCE = sce, useAssay = "decontXcounts", normAssayName = "seuratNormData", normalizationMethod = "LogNormalize", scaleFactor = 10000)

2. Highly Variable Genes
Highly variable genes can be identified by first using the runSeuratFindHVG() function that computes that statistics against a selected HVG method in the rowData of input object. The genes can be identified by using the getTopHVG() function. The variable genes can be visualized using the plotSeuratHVG() method.
Parameters for runSeuratFindHVG() include useAssay (specify the name of the assay to use), method (specify the method to use for variable genes computation, defaults to "vst"), hvgNumber (number of variable features to identify, defaults) and createFeatureSubset (define a name for the subset of selected top variable features).

sce <- runSeuratFindHVG(inSCE = sce, useAssay = "decontXcounts", method = "vst", hvgNumber = 2000, createFeatureSubset = "hvf")

# Print names of top 10 variable features
print(getTopHVG(inSCE = sce, method = "vst", hvgNumber = 10))
##  [1] "GNLY"    "LYZ"     "IGLL5"   "S100A9"  "PPBP"    "FTL"     "PF4"    
##  [8] "FTH1"    "FCER1A"  "HLA-DRA"
# Plot variable features with top 10 labeled
plotSeuratHVG(sce, labelPoints = 10)

3. Dimensionality Reduction
PCA or ICA can be computed using the runSeuratPCA() and runSeuratICA() functions respectively. Plots can be visualized using plotSeuratReduction(), plotSeuratElbow(), plotSeuratJackStraw() (must previously be computed by runSeuratJackStraw()) and runSeuratHeatmap().

sce <- runSeuratPCA(inSCE = sce, useAssay = "seuratNormData", reducedDimName = "pca", nPCs = 50, seed = 42, scale = TRUE, useFeatureSubset = "hvf")

# Plot PC1 vs PC2 plot
plotSeuratReduction(inSCE = sce, useReduction = "pca")

# Plot Elbowplot
plotSeuratElbow(inSCE = sce)
# Compute JackStraw
sce <- runSeuratJackStraw(inSCE = sce, useAssay = "seuratNormData", dims = 50)

# Plot JackStraw
plotSeuratJackStraw(inSCE = sce, dims = 50)

# Compute and plot first 4 dimensions and 30 top features in a heatmap
runSeuratHeatmap(inSCE = sce, useAssay = "seuratNormData", useReduction = "pca", nfeatures = 30, dims = 4)

## NULL

4. tSNE/UMAP
runSeuratTSNE() and runSeuratUMAP() can be used to compute tSNE/UMAP statistics and store into the input object. Parameters to both functions include inSCE (input SCE object), useReduction (specify the reduction to use i.e. "pca" or "ica"), reducedDimName (name of this new reduction) and dims (number of dims to use). Additionally, method specific parameters can be used to fine tune the algorithm. plotSeuratReduction() can be used to visualize the results.

sce <- runSeuratTSNE(inSCE = sce, useReduction = "pca", reducedDimName = "seuratTSNE", dims = 10, perplexity = 30, seed = 1)

# Plot TSNE
plotSeuratReduction(sce, "tsne")

sce <- runSeuratUMAP(inSCE = sce, useReduction = "pca", reducedDimName = "seuratUMAP", dims = 10, minDist = 0.3, nNeighbors = 30, spread = 1, seed = 42)
# Plot UMAP
plotSeuratReduction(sce, "umap")

5. Clustering
runSeuratFindClusters() function can be used to compute the clusters, which can later be plotted through the plotSeuratReduction() method with cluster labels. The parameters to the function include inSCE (input SCE object), useAssay (name of the assay if no reduction to be used), useReduction (specify which reduction to use i.e. "pca" or "ica"), dims (number of dims to use), the algorithm (either "louvain", "multilevel" or "SLM") and resolution (defaults to 0.8).

sce <- runSeuratFindClusters(inSCE = sce, useReduction = "pca", resolution = 0.8, algorithm = "louvain", dims = 10) 
## Modularity Optimizer version 1.3.0 by Ludo Waltman and Nees Jan van Eck
## 
## Number of nodes: 2627
## Number of edges: 94434
## 
## Running Louvain algorithm...
## Maximum modularity in 10 random starts: 0.8366
## Number of communities: 10
## Elapsed time: 0 seconds

plotSeuratReduction() can then be used to plot all reductions previously computed with cluster labels:

plotSeuratReduction(sce, "pca", showLegend = TRUE)

plotSeuratReduction(sce, "tsne", showLegend = TRUE)

plotSeuratReduction(sce, "umap", showLegend = TRUE)

6. Find Markers
Marker genes can be identified using the runSeuratFindMarkers() function. This function can either use one specified column from colData of the input object as a group variable if all groups from that variable are to be used (allGroup parameter) or users can manually specify the cells included in one group vs cells included in the second group (cells1 and cells2 parameter).

sce <- runSeuratFindMarkers(inSCE = sce, allGroup = "Seurat_louvain_Resolution0.8")

# Fetch marker genes table
markerGenes <- metadata(sce)[["seuratMarkers"]] 

# Order by log fold change and p value
markerGenes <- markerGenes[order(-markerGenes$avg_log2FC, markerGenes$p_val),]

head(markerGenes)
##        gene.id         p_val avg_log2FC pct.1 pct.2     p_val_adj cluster1
## 4005      FTH1 4.997420e-109   278.7536 1.000 0.989 1.636055e-104        4
## 10628   TMSB10  2.717889e-03   276.9933 0.981 0.996  1.000000e+00        8
## 2218       IGJ  1.929933e-43   246.5411 0.139 0.010  6.318216e-39        2
## 10821  HLA-DRA  1.218767e-22   216.9737 1.000 0.557  3.989998e-18        9
## 10200     ACTB  3.659473e-07   184.3126 1.000 0.993  1.198038e-02        8
## 10813 HLA-DPB1  4.213390e-24   172.2514 1.000 0.516  1.379380e-19        9
##       cluster2
## 4005       all
## 10628      all
## 2218       all
## 10821      all
## 10200      all
## 10813      all

The marker genes identified can be visualized through one of the available plots from ridge plot, violin plot, feature plot, dot plot and heatmap plot. All marker genes visualizations can be plotted through the wrapper function plotSeuratGenes(), which must be supplied the SCE object (markers previously computed), name of the scaled assay, type of the plot (available options are "ridge", "feature", "violin", "dot" and "heatmap"), features that should be plotted (character vector) and the grouping variable that is available in the colData slot of the input object. An additional parameter ncol decides in how many columns should the visualizations be plotted.

plotSeuratGenes(inSCE = sce, useAssay = "seuratNormData", plotType = "ridge", features = metadata(sce)[["seuratMarkers"]]$gene.id[1:4], groupVariable = "Seurat_louvain_Resolution0.8", ncol = 2, combine = TRUE)

References

[1]
A. Butler, P. Hoffman, P. Smibert, E. Papalexi, and R. Satija, “Integrating single-cell transcriptomic data across different conditions, technologies, and species,” Nature Biotechnology, vol. 36, no. 5, pp. 411–420, Jun. 2018, doi: 10.1038/nbt.4096.
[2]
T. Stuart et al., “Comprehensive Integration of Single-Cell Data,” Cell, vol. 177, no. 7, pp. 1888–1902.e21, Jun. 2019, doi: 10.1016/j.cell.2019.05.031.
[3]
R. Satija, J. A. Farrell, D. Gennert, A. F. Schier, and A. Regev, “Spatial reconstruction of single-cell gene expression data,” Nature Biotechnology, vol. 33, pp. 495–502, 2015, doi: 10.1038/nbt.3192.
[4]
Y. Hao et al., “Integrated analysis of multimodal single-cell data,” Cell, 2021, doi: 10.1016/j.cell.2021.04.048.
[5]
R. A. Amezquita et al., “Orchestrating single-cell analysis with Bioconductor,” Nature Methods, vol. 17, no. 2, pp. 137–145, Feb. 2020, doi: 10.1038/s41592-019-0654-x.