Single-cell RNAseq Clustering
Introduction
This tool is designed to perform the clustering of cells coming from Single-Cell RNA Sequencing (scRNA-seq) data. That is, it aims to find groups of cells that share similar expression patterns, which should correspond to the same cell type or state. Prior to the clustering, this tool performs the preprocessing of the data in order to make it suitable for the clustering algorithm (Figure 1). This application is based on the widely-used Seurat package.
Please cite Seurat as:
Figure 1. General workflow for the single-cell RNA seq clustering analysis
Run single-cell RNA sequencing clustering
In order to perform the Single Cell Clustering, a Count Matrix object has to be opened. It can be loaded from different formats by going to transcriptomics → Load → Single Cell RNA-Seq Count Matrix or generated from FASTQ sequencing files with the Single Cell RNA-Seq Quantification tool available in transcriptomics → Single Cell RNA-Seq → Single Cell RNA-Seq Quantification.
Once loaded, go to the Side Panel → Actions → Clustering (Figure 2).
Figure 2. Launch clustering analysis from a scRNA-Seq Count Matrix Side Panel.
Configuration: Preprocessing.
These steps are meant to prepare data for the clustering analysis (Figure 3).
Normalization
Normalization aims to reduce differences in gene expression due to technical variation. This step tries to ensure that the observed heterogeneity within cells is due to biological reasons, rather than technical biases.
Normalize Data. Check this option to perform the normalization step. Extremely recommendable unless your data is already normalized or you are sure that your data does not need it.
Normalization Method. Which normalization method to use:
Regularized Negative Binomial Regression. This option uses the SCTransform normalization method designed by the Seurat team. This method applies a generalized linear model for each gene, with counts as the response and sequencing depth as the explanatory variable. This is used to learn the parameters of the model that depend on a gene’s average expression. Then, a negative binomial regression is applied and the Person residuals (the difference between the expected and the observed expression) are treated as the normalized expression levels. For more details about the statistical method please see Hafemeister and Satija 2019.
Log Normalization. Feature (gene) counts for each cell are divided by the total counts for that cell and multiplied by a scale factor. This is then natural-log transformed using log1p.
Relative Counts. Feature counts for each cell are divided by the total counts for that cell and multiplied by a scale factor. No log transformation is applied.
Centered Log Ratio Transformation (CLR). For each feature, the CLR-transformation is defined as the logarithm of the feature counts for one cell, divided by the geometric mean of all counts for that cell.
Data Adjustment
These steps are necessary whether the normalization is done or not, in order to prepare data for the dimensional reduction step.
High Variable Genes. Number of high variable genes to keep for further analysis. Keeping only those genes reduces both computational cost and non-relevant signals.
Scale Data. If checked, it scales features to have unit variance.
Center Data. If checked, it centers features to have zero mean.
Data Correction
Regress Out Mitochondrial Genes. Check this option to reduce the heterogeneity between cells associated with the expression of mitochondrial genes. This could prevent grouping cells during the clustering step that present a higher expression of those genes. In some cases, this information could not be informative since the higher expression of mitochondrial genes could be caused by technical reasons (e.g. because the mRNA has been leaked during cell manipulation) rather than by biological reasons. In order to perform this analysis, a file with a list of mitochondrial genes must be provided with one gene per line (Figure 4).
Regress Out Cell Cycle Genes. Check this option to reduce the heterogeneity between cells associated with the expression of cell cycle genes. This could prevent grouping cells during the clustering step that are in the same developmental stage, independently of their type. In some cases, these differences in cell cycle may be uninformative but, in other cases, they could be treated as indicative of proliferating cell populations which can be different across treatment conditions, for example. So whether to check this option or not will depend on the dataset under study and the target of the experiment. In order to perform this analysis, a file with cell cycle genes in one column and the cell cycle phase they belong to (S or G2/M) must be provided. One gene per line and columns separated by a tab (Figure 5).
Dimensional Reduction
Common single-cell RNAseq analysis implies from hundreds to thousands of cells and tens of thousands of genes. That is, it is really high dimensional data, so it is advisable to reduce the dimensionality of the dataset in order to reduce the computational cost of further analysis. The most widely used method for dimensional reduction is Principal Component Analysis (PCA). Keeping only the firsts Principal Components for further analysis reduces the dimensionality of the data while keeping the heterogeneity of the dataset, since they explain the largest amount of variance present in the sample. For more details please see "Orchestrating Single-Cell Analysis", 2020.
Principal Components. Number of Principal Components to compute from the Principal Component Analysis. It must be smaller than the number of cells present in the count table.
Figure 3. Preprocessing wizard page.
Figure 4. Mitochondrial genes file example.
Figure 5. Cell cycle genes file example.
Configuration: Multi-sample Data Integration.
This is an additional step needed when the input count table contains data coming from multiple samples e.g. from wild-type and mutant organisms, from control and stimulated samples, etc. All these situations could cause changes in gene expression that could make a joint analysis of all the data difficult, with cells clustering both by condition and by cell type (Figure 6-A). In order to avoid that, the Multi-sample Data Integration step aims to integrate scRNA-seq datasets by identifying common cell types based on common sources of variation. As a result, this step enables the identification of shared populations across datasets (Figure 6-B and C) and thus further downstream comparative analyses. For more details about the integration algorithm please read the papers from the Seurat Team: Butler et al. 2018 and Stuart et al. 2019.
The integration algorithm consists of two steps: it firsts computes a dimensional reduction using Canonical Correlation Analysis (CCA) and then performs a mutual nearest neighbors (MNNs) algorithm to identify shared subpopulations across datasets. In this way, it finds “anchors” to integrate the datasets, that is, equivalent cells across datasets.
The next parameters modify the behavior of the integration algorithm. In case the data integration has not been successful, it is advisable to modify these parameters to try to improve the results. See Figure 7-A for a “bad” and Figure 7-B for a “good” integrations.
Integration Factor. Choose the condition to integrate datasets by. That is, the condition that you don’t want to interfere during the clustering of cells. If no integration is to be done, select “None”.
N. Dimensions for Integration. The number of dimensions to use from the Canonical Correlation Analysis for the integration step.
K Anchor. How many neighbors (k) to use during the MNNs when picking anchors.
K Filter. How many neighbors (k) to use during the MNNs when filtering anchors.
K Score. How many neighbors (k) to use during the MNNs when scoring anchors.
K Weight. How many neighbors (k) to use during the MNNs when weighing anchors.
The Experimental Design specified with the tool “Experimental Desing” available on the Count Matrix' Side Panel will be shown. If the Count Matrix only contains one sample, the Multi-sample Data Integration parameters will be desabled. In order to perform Data Integration, perform a “Merge Single Cell Counts” first.
Figure 6. UMAP plots of 10,799 cells coming from wild type (WT) and Rb mutant Drosophila eye disc, prior to (A) and post (B) alignment. After alignment, cells across conditions are grouped together based on shared cell type, allowing for a single joint clustering (C) to detect 11 populations.
Figure 7. Multi-sample Data Integration configuration wizard page.
Configuration: Clustering.
This step performs the actual cell clustering. It groups cells with similar expression patterns, which should correspond to the same cell type. For this analysis, Seurat uses a graph-based clustering algorithm. For more information, please see "Orchestrating Single-Cell Analysis", 2020.
Clustering
These parameters affect the clustering algorithm (Figure 8).
Define Dimensions by. How to decide the number of dimensions to use from the dimensional reduction (PCA) for the clustering step.
Elbow Point. This option automatically decides the number of dimensions to use. The assumption is that each of the top PCs capturing biological signals should explain much more variance than the remaining PCs. Thus, there should be a sharp drop in the percentage of variance explained from the last “biological” PC to the others. This is the so-called “Elow Point”.
Manual. With this option, it is necessary that the user specifies the number of PCs to use in the “Number of Dimensions” parameter.
It should be taken into account that strong biological variation in the early PCs will shift the elbow point, potentially excluding weaker (but still interesting) variation in the next PCs immediately following it. So it could be interesting to repeat the analysis increasing the number of dimensions established by the Elbow Point to see if it improves the results.
Number of Dimensions. The number of dimensions to use from the dimensional reduction for the clustering step. This option is only available if the “Manual” option from the “Dimensions Choice” is selected.
k-value. The number of neighboring cells computed during the clustering algorithm. Normally, a greater k-value would produce a smaller number of clusters, and vice versa.
Resolution. This parameter determines how "fine" the clustering is: values above 1.0 would produce a larger number of clusters, and vice versa.
UMAP Configuration
These parameters affect only the UMAP visualization. The Uniform Manifold Approximation and Projection (UMAP) method is a non-linear dimensionality reduction technique (similar to PCA), but this time it is used to plot the high-dimensional single-cell data into two dimensions. The UMAP is used for visualization purposes because it represents the variability of single-cell data more accurately than the PCA. For more information, please refer to "Orchestrating Single-Cell Analysis", 2020.
Point’s Minimum Distance. This controls how tightly are the points (cells) compressed together.
Point’s Spread. This controls how expanded the points are. In combination with the minimum distance, this determines how clustered the points are.
Figure 8. Clustering configuration wizard page.
Results
Once the input counts have been processed and analyzed via the "scRNA-seq Clustering '' tool, a new tab is opened containing the Clustering Results (Figure 9). This new tab contains one row for each of the obtained clusters and the number of cells composing them.
It also generates a Summary Report (Figure 10):
The “Data Overview” table shows some general statistics about the data.
The “Experimental Design” section informs about the conditions present in the Experimental Design, if previously specified.
The “Clustering Results” table shows the total number and a link to a list of cells belonging to each cluster. Before the table, there’s one line specifying the number of dimensions from the PCA used for the clustering step.
The “Analysis Parameters” table shows the parameters used for the clustering analysis.
In addition, two charts are generated:
Elbow Plot (Figure 11). It shows the amount of variance (given by the standard deviation) explained by successive PCs. This helps to decide how many principal components to use for the clustering algorithm, in case you want to re-run the clustering with different parameters.
UMAP (Figure 12). Each dot is a cell, colored by the cluster assigned by the graph-based clustering algorithm. The X and Y axes are the values in the UMAP coordinate 1 and 2 of each cell, respectively.
Figure 9. Clustering Results main viewer.
Figure 10. Part of the Summary Report.
Figure 11. Elbow Plot example.
Figure 12. UMAP plot example.
Side Panel Features
Summary Report
It shows the Summary Report previously explained in the above “Results” section (Figure 10).
Extract Cluster(s) Counts
It may be the case that one or more big clusters have been obtained during the analysis, so it could be desirable to extract them and re-run the clustering to obtain sub-clusters of cells within it. That could make it possible to find more specific cell types.
With the Extract Cluster(s) feature, it is possible to specify one or more clusters to extract by checking them in the “Cluster(s) to Extract” option of the wizard (Figure 12). The result is a new tab with a normalized count table containing only the cells that belong to the selected cluster(s). A new clustering analysis can be performed with it.
Differential Expression
Perform Differential Expression Analysis with the Clustering results. For more detailed information, please go to the corresponding section (link).
Figure 13. Extract Cluster(s) wizard.
Charts
UMAP
It shows the UMAP plot. It has different coloring options available in the wizard (Figure 14):
Color by Metadata. With this option, the same UMAP as in Figure 11 will be obtained. It is possible to choose in the chart’s Side Panel coloring cells by the obtained clusters or by the factors in the experimental design (if any).
Color by Gene Expression. Cells will be colored according to the expression level of the specified gene (Figure 15).
Expression Profile
With this feature, it is possible to see the expression levels of known gene markers across the different clusters in order to identify putative cell types. To this end, a Bubble Plot is generated with clusters in rows and the specified genes in columns (Figure 16). The size of the dot represents the percentage of cells expressing the gene, that is, the percentage of cells that have a gene expression level greater than 0. The color represents the average gene expression in that cluster. You can configure the following options in the wizard (Figure 17):
Input genes.
Input Genes. You can specify here which genes to plot. The gene name or ID should correspond to the one in the input count table used during the clustering analysis. They can be specified in two ways:
Text: write the genes to plot in the “Genes List” text box, one per line.
File: specify a file containing the genes to plot, one gene per line.
Plot Options.
This affects the visualization.
Scale Gene Expression. When checked, it applies the Z-Score transformation to scale average gene expression across clusters. It allows the visualization of both highly and lowly expressed genes on the same color scale. It should be noted that this may exaggerate the results, but it is still advisable if you are going to plot genes with different expression level ranges. If unchecked, it colors the dots by the raw average gene expression of each cluster.
Figure 14. UMAP wizard with different coloring options.
Figure 15. UMAP with cells colored by the expression of the gene TFR2.
Figure 16. Bubble Plot example.
Figure 17. Expression Profile Wizard.
Export
Export Cell Metadata
It generates a file containing one row per each cell, containing the Cell ID, the cluster label assigned, and the sample of origin (Figure 17). Columns are separated by a tab.
Figure 18. Cell Metadata file.
Context menu options
Rename Cluster
Once you know the cell type of one cluster, it could be interesting to rename it so it has a more descriptive name. In order to do that, go to the cluster column and click the mouse right button. It will open a list of options. On that list, click on “Rename Cluster”. A new wizard will appear to type the new name. The result is the same project but with the renamed cluster (Figure 19). It has to be saved in order to conserve the modification.
Figure 19. Rename a cluster (column) from a Seurat clustering object.