Introduction to robustSingleCell

We used two replicates of CD44+ T cell data sets from Ciucci et al. 2019 1 as an example to demonstrate the use of robustSingleCell. The analysis requires at least 8G of memory on slurm 2 high performance computing workload manager (for example, you can start by requesting srun --pty -p <partition> --mem=8G -t 1:00:00 bash to start an interactive session).

We first download the raw 10X data from GEO using GEOquery, which can be obtained using the following command if not already installed:


The two datasets LCMV1, LCMV2 will be downloaded into TMPDIR. Each folder will contain the matrix.mtx, gene.tsv, barcode.tsv files as in 10X genomics format.


We cluster each dataset separately to account for dataset-specific technical and biological differences. Then, we measure the transcriptional similarity and divergence between the clusters identified in the two datasets using correlation analysis.

Individual analysis of LCMV1 and LCMV2

First, we set up the directory where the results of the analysis will be stored.

LCMV1 <- initialize.project(datasets = "LCMV1", 
                          origins = "CD44+ cells",
                          experiments = "Rep1",
                          data.path = file.path(tempdir(), "LCMV"),
                          work.path = file.path(tempdir(), "LCMV/LCMV_analysis")) function reads the data in 10X genomics format and performs quality filtering as described in Magen et al 2019 3. We randomly downsampled the datasets to 1000 cells to shorten the simplify this example.

LCMV1 <-, subsample = 500)

Next, we identify highly variable genes for the following PCA and clustering analyses. We also compute the activation of gene sets of interest, such as cell cycle genes, for confounder correction.

LCMV1 <- get.variable.genes(LCMV1) 
exhaustion_markers <- c('Pdcd1', 'Cd244', 'Havcr2', 'Ctla4', 'Cd160', 'Lag3', 'Tigit', 'Cd96')
LCMV1 <- add.confounder.variables(LCMV1,
    ribosomal.score = ribosomal.score(LCMV1),
    mitochondrial.score = mitochondrial.score(LCMV1),
    cell.cycle.score = cell.cycle.score(LCMV1),
    Exhaustion = controlled.mean.score(LCMV1, exhaustion_markers))

Figure 1 shows the mitochondrial score versus number of UMIs, pre and post filtering.

Fig 1. Mitochondrial genes score vs. number of UMIs for pre (top) and post (bottom) quality control filtering.

The PCA function performs multiple simulation analyses of shuffled data to determine the appropriate number of PCs. You can also run each simulation in parallel using the option local = F.

LCMV1 <- PCA(LCMV1, local = T)

We then perform clustering analysis for a range of clustering resolutions. The analysis is repeated multiple times over shuffled data to estimate the appropriate clustering resolution and control for false discovery of clusters. At the end of the clustering, the function will prompt you to choose an optimal clustering resolution. We choose 0.05 for our KNN ratio, which is the smallest value tested with mdlrty/mean.shfl > 2.

LCMV1 <- cluster.analysis(LCMV1, local = T)
Fig 2. Bar plot shows the clustering modularity of the original data versus shuffled data across multiple clustering resolutions. Numbers on top represent the fold change of original versus shuffled analysis for each resolution.

We select the appropriate resolution, typically the one where there is more than two (2) fold change modularity difference relative to the shuffled analysis.

The summarize function which performs differential expression analysis, computes tSNE and visualizes the results in the analysis folder. After differential expression analysis, get.cluster.names assigns clusters with names using a customized set of marker genes which users should adapt to their own data.

types = rbind(
                data.frame(type='CD4', gene = c("Cd4")),
summarize(LCMV1, local = T)
LCMV1_cluster_names <- get.cluster.names(LCMV1, types, min.fold = 1.0, max.Qval = 0.01)
LCMV1 <- set.cluster.names(LCMV1, names = LCMV1_cluster_names)
summarize(LCMV1, local = T)

Figure 3 shows violin plots indicating the activation of the cell cycle genes.

Fig 3. Violin plot pf cell cycle score.

Figure 4 places individual cells on a two dimensional grid corresponding to the scores of the first two PCs (note that the PCA figures are created in the next step via summarize function below).

Fig 4. Single cells placement on a 2D grid corresponding to the first two PCs.

The genes driving the PCs are visualized in figure 5 according to the PCA loadings after removing the lowly ranked genes.

Fig 5. Top ranked genes contribution to PC1 and PC2 scores.

The average expression of genes driving the PCs can be visualized as a heatmap visualized in figure 6 according to the PCA loadings after removing the lowly ranked genes.

Fig 6. Heatmap shows loadings of the first PC.

Figure 7 shows the tSNE visualization of the cells, color coded by cluster assignment.

Fig 7. t-SNE plot colored by cluster assignment.