The easiest way to use scigenex is to perform the following steps using the Seurat R package:
The resulting object can be used as input to SciGeneX. You can also provide a normalized matrix as input.
For this tutorial, we’ll be using a peripheral blood mononuclear cell (PBMC) scRNA-seq dataset available from the 10x Genomics website. This dataset contains 2700 individual cells sequenced on the Illumina NextSeq 500 and can be downloaded from 10X Genomics web site or via the SeuratData
library.
In this step, we’ll carry out the classic pre-processing steps of the tutorial. Please refer to this tutorial for more information. If you have already pre-processed your data with Seurat, or if you have a normalized count matrix as input, you can skip this step.
library(ggplot2)
library(patchwork)
library(Seurat, quietly = TRUE)
library(SeuratData)
InstallData("pbmc3k")
data("pbmc3k")
We next run the classical steps of the seurat pipeline. For more information, you can check Seurat website here.
# Quality control
pbmc3k[["percent.mt"]] <- PercentageFeatureSet(pbmc3k, pattern = "^MT-")
pbmc3k <- subset(pbmc3k, subset = percent.mt < 5 & nFeature_RNA > 200)
# Normalizing
pbmc3k <- NormalizeData(pbmc3k)
# Identification of highly variable genes
pbmc3k <- FindVariableFeatures(pbmc3k, selection.method = "vst", nfeatures = 2000)
# Scaling data
pbmc3k <- ScaleData(pbmc3k, features = rownames(pbmc3k), verbose = FALSE)
# Perform principal component analysis
pbmc3k <- RunPCA(pbmc3k, features = VariableFeatures(object = pbmc3k), verbose = FALSE)
# Cell clustering
pbmc3k <- FindNeighbors(pbmc3k, dims = 1:10, verbose = FALSE)
pbmc3k <- FindClusters(pbmc3k, resolution = 0.5, verbose = FALSE)
# Dimension reduction
pbmc3k <-suppressWarnings(RunUMAP(pbmc3k, dims = 1:10, verbose = FALSE))
DimPlot(pbmc3k, reduction = "umap")
In this section, we use the previously generated Seurat object as input to run the main SciGeneX commands. This command executes the main algorithm which will:
First, we’ll load the SciGeneX library. To limit the verbosity of the SciGeneX functions, we’ll set the verbosity level to zero (which will disable information and debug messages).
library(scigenex)
scigenex::set_verbosity(0)
Then we call successively the select_genes()
function which will select informative genes (i.e co-regulated), then the gene_clustering()
function which will call MCL and partition the dataset into gene modules.
# Select informative genes
pbmc_scigenex <- select_genes(pbmc3k,
k = 50,
distance_method = "pearson",
which_slot = "data",
row_sum=2)
# Run MCL
pbmc_scigenex <- gene_clustering(pbmc_scigenex,
s = 5,
threads = 2,
inflation = 2)
The object produced is a ClusterSet
objet that is a S4 object that is intented to store gene modules.
isS4(pbmc_scigenex)
## [1] TRUE
pbmc_scigenex
## An object of class ClusterSet
## Name: s8gZ5uYQ33
## Memory used: 80388912
## Number of cells: 2643
## Number of informative genes: 1854
## Number of gene clusters: 163
## This object contains the following informations:
## - data
## - gene_clusters
## - top_genes
## - gene_clusters_metadata
## - gene_cluster_annotations
## - cells_metadata
## - dbf_output
## - parameters
## * distance_method = pearson
## * k = 50
## * noise_level = 5e-05
## * fdr = 0.005
## * row_sum = 2
## * no_dknn_filter = FALSE
## * seed = 123
## * keep_nn = FALSE
## * k_mcl_graph = 5
## * output_path = /var/folders/zy/wl3dj2_n76zfc8sdvny1q06c0000gn/T//RtmpUE4KyV
## * name = s8gZ5uYQ33
## * inflation = 2
There are various methods associated with the ClusterSet
objects.
## [1] "[" "%in%"
## [3] "clust_names" "clust_size"
## [5] "cluster_set_to_xls" "cluster_stats"
## [7] "col_names" "dim"
## [9] "enrich_go" "gene_cluster"
## [11] "grep_clust" "module_quality_scores"
## [13] "nclust" "plot_clust_enrichments"
## [15] "plot_ggheatmap" "plot_markers_to_clusters"
## [17] "rename_clust" "reorder_clust"
## [19] "reorder_genes" "row_names"
## [21] "show" "top_by_go"
## [23] "top_genes" "viz_enrich"
## [25] "which_clust"
The current object contains 1854 informative genes, 2643 samples and 163 gene modules.
nrow(pbmc_scigenex)
## [1] 1854
ncol(pbmc_scigenex)
## [1] 2643
nclust(pbmc_scigenex)
## [1] 163
At this stage, several modules need to be filtered, as many of them may be singletons. Interestingly, the ClusterSet
class implements the indexing operator/function (“[”). The first argument/dimension of the indexing function corresponds to the cluster stored in the object. The second dimension corresponds to the cell/spot names. As an example, we can simply store gene modules whose size (*i.e. number of genes) is greater than 7 using the following code. The result is an object containing gene modules 40.
pbmc_scigenex <- pbmc_scigenex[clust_size(pbmc_scigenex) > 7, ]
nclust(pbmc_scigenex)
## [1] 40
It may also be important to filter out gene based on dispersion. Several parameters can be computed for each cluster using the cluster_stats()
function.
plot_cluster_stats(cluster_stats(pbmc_scigenex)) +
ggplot2::theme(axis.text.y = ggplot2::element_blank(),
axis.text.x = element_text(angle=45, vjust = 0.5),
panel.grid = element_blank())
Here we will select gene modules based on standard deviation (> 0.1) and rename the cluster (from 1 to the number of clusters):
pbmc_scigenex <- pbmc_scigenex[cluster_stats(pbmc_scigenex)$sd_total > 0.1, ]
pbmc_scigenex <- rename_clust(pbmc_scigenex)
nclust(pbmc_scigenex)
## [1] 18
Then we check the statistics again.
plot_cluster_stats(cluster_stats(pbmc_scigenex)) +
ggplot2::theme(axis.text.y = ggplot2::element_blank(),
axis.text.x = element_text(angle=45, vjust = 0.5),
panel.grid = element_blank())
Clusters of genes are stored in the gene_clusters
slot. One can access the gene names from a cluster using the get_genes()
command. By default, all genes are returned.
# Extract gene names from the 5th gene cluster
genes_module_5 <- get_genes(pbmc_scigenex, cluster = 5)
head(genes_module_5)
## [1] "PTCRA" "CLDN5" "C2orf88" "ACRBP" "SPOCD1" "TAL1"
One can also access the gene to cluster mapping using get_genes()
.
head(gene_cluster(pbmc_scigenex))
## CST3 TYROBP AIF1 LST1 LYZ FTL
## 1 1 1 1 1 1
tail(gene_cluster(pbmc_scigenex))
## PRR7-AS1 GPR153 CTD-2210P24.4 C6orf229 RP11-279F6.3
## 18 18 18 18 18
## MRVI1-AS1
## 18
Visualizing a heatmap containing all cells and modules can be time-consuming and require significant memory resources. A first alternative is to examine gene modules individually or a subset of modules. Gene modules can be visualized using the plot_heatmap()
the plot_ggheatmap()
functions. The former is primarily intended to provide an interactive visualization based on the iheatmapr
library, and allows easy browsing of results and zooming in on particular regions of the heatmap. This second solution leads to a diagram that is more easily customized, as it is based on the ggplot framework. Here, we use plot_ggheatmap()
to view the first 4 clusters. Note that here, we choose to order the columns/cells on the results of Seurat::FindClusters
.
plot_ggheatmap(pbmc_scigenex[1:4, ],
use_top_genes = FALSE,
ident=Idents(pbmc3k)) + ggtitle("Cluster 1 to 4")
However, an alternative is to extract the most representative genes from each group. This can be achieved using the top_genes()
function. This function stores the identifiers of these representative genes in the top_genes
slot of the ClusterSet
object. The get_genes()
function is used to access the top_genes
slot.
pbmc_scigenex <- top_genes(pbmc_scigenex)
genes_cl5_top <- get_genes(pbmc_scigenex, cluster = 5, top = TRUE)
genes_cl5_top
## [1] "SPOCD1" "CLDN5" "C2orf88" "PTCRA"
## [5] "ITGB3" "TAL1" "ACRBP" "ESAM"
## [9] "TNS1" "HIST1H2BJ" "ATP9A" "RP11-359I18.5"
## [13] "PDZK1IP1" "C19orf33" "ALOX12" "PLA2G4C"
## [17] "LNP1" "TSC22D1" "SYTL4" "TGFB1I1"
Both plot_heatmap()
and plot_ggheatmap()
support the use_top_genes
argument:
plot_ggheatmap(pbmc_scigenex,
use_top_genes = TRUE,
ident=Seurat::Idents(pbmc3k)) + ggtitle("All clusters (top genes)") +
theme(strip.text.y = element_text(size=4))
A very interesting feature of SciGeneX is its ability to display gene expression levels in cells/spots using interactive heatmaps. With this function, the user can interactively evaluate expression levels in selected cells or groups, or over the whole dataset. However, it is generally advisable, when using all clusters, to restrict the analysis using top_genes=TRUE
. Here, we’ll display expression levels in cells from clusters 3 to 5.
plot_heatmap(pbmc_scigenex[3:5, ],
use_top_genes = TRUE,
cell_clusters =Seurat::Idents(pbmc3k))