2) Clustering and cell type annotation
September 25, 2026
This script integrates the samples into a joint UMAP graph and estimates clusters. Marker genes (= upregulated differentially expressed genes) are calculated for each cluster and provided to the user for annotation of the cell types. The clustering resolution can be adjusted based on the desired level of granularity, allowing either a rough or a more detailed cell-type annotation. Initially 2 different cluster granularities are provided to the user.
1 Setup
#Load helper functions
source("/people/nrq364/NeuSiC/scRNA_helper_new.R")
library(magrittr)
library(dplyr)
library(conos)
library(pagoda2)
library(qs2)
library(ggplot2)
library(ggrastr)
library(openxlsx)
library(tidyverse)2 Read in data
Read in filtered list of count matrices generated in notebook 1) Pre-processing and QC filtering.
cms <- qs_read("cms_filtered.qs2", nthreads = 10)Table with dimensions of count matrices: genes x cells. As well as total number of cells.
sapply(cms, dim)## control_1 control_2 control_4 control_5 control_6 control_7 treated_1 treated_2 treated_3 treated_4 treated_6 treated_7
## [1,] 32285 32285 32285 32285 32285 32285 32285 32285 32285 32285 32285 32285
## [2,] 898 2341 1579 790 3468 2505 1080 754 1045 2198 2695 3032# sum of all cells
sum(sapply(cms, function(x) dim(x)[2]))## [1] 223853 Embed cells in a joint UMAP graph using conos
Create a vector with sample names.
names <- paste0(names(cms), "_")Use a helper function (defined in the script scRNA_helper_new.R) to embed cells from all samples in a joint UMAP graph and perform first clustering. This helper function performs pagoda2 pre-processing and normalization of each count matrix, then creates a conos object and integrates the samples into a joint UMAP graph. Alignment strength can be adjusted to force stronger integration of samples.
con <- quickConos(cms,
names,
n.cores.p2=10,
n.cores.con=20, get.tsne = F, alignment.strength=0)
con <- con$conInvestigate how cells cluster on the UMAP. Are there any batches visible (clustering by sample, by condition)? If yes, alignment.strength parameter can be increased. The embedding used here was created with alignment.strength = 0.2.
Plot UMAP coloured by sample.
data <- con$plotGraph(groups=con$getDatasetPerCell())
data <- data$data
ggplot(data, aes(x = x, y =y, color= Group))+ geom_point(size=0.2, alpha=0.5)+ theme_light()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ xlab("")+ylab("") + theme(axis.ticks = element_blank(), axis.text =element_blank())+labs(color = "Sample") +guides(color = guide_legend(override.aes = list(size=3, alpha=1)))+scale_color_manual(values=viridisLite::turbo(12))
#ggsave("img/2_clustering/fig0.png")
Plot UMAP coloured by cluster.
con$plotGraph() + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ xlab("")+ylab("") + theme(axis.ticks = element_blank(), axis.text =element_blank())
#ggsave("img/2_clustering/fig1.png")
4 Metadata table
Generate one table/tibble with sample, condition, depth, mito.fraction, doublet.scores ordered by barcodes/cell names of conos object.
4.1 Metadata from conos
From conos object we save conos_id, sample, condition, and cell barcode in a table/tibble.
clusters <- con$clusters$leiden$groups
# initiate the data frame
df <- tibble(con_id = names(clusters), cluster = as.integer(as.character(clusters)))
df %<>% separate(con_id, into = c("sample", "barcode"), sep = "!!", remove=F)
df %<>% mutate(condition = sub("_[0-9]$", "", sample)) # generate condition column by replacing the "_x" part with "" (nothing)
#df %<>% mutate(barcode = sub("^[^_]+_[0-9]+_", "", barcode))
# add a "order" column to see in case rows get mixed up. Necessary to keep same order as in conos object.
df %<>% mutate(order = row_number())
head(df)4.2 Metadata from crm object
We add doublet scores, depth (UMI counts per cell), and mito gene fraction from the CRMetrics object to the metadata table.
Get doublet scores and add to metadata table
scrublet <- crm$doublets$scrublet$result$scores
names(scrublet) <- row.names(crm$doublets$scrublet$result)
new_column <- tibble(con_id = row.names(crm$doublets$scrublet$result),scrublet = crm$doublets$scrublet$result$scores)
df %<>% left_join(new_column, by = "con_id")Get depth and add to metadata table.
depth <- crm$getDepth()
new_column <- tibble(con_id = names(depth), UMI_counts = depth)
df %<>% left_join(new_column, by = "con_id")Get mito gene fraction and add to metadata table.
mito <- crm$getMitoFraction()
new_column <- tibble(con_id = names(mito), mito_gene_fraction = mito)
df %<>% left_join(new_column, by = "con_id")head(df)5 UMAPs coloured by metadata
We colour the UMAP by different variables from the metadata to evaluate the embedding, clustering and QC filtering.
# How to use the metadata to colour UMAP: generate a named vector feed into "groups" parameter
setNames(df$desired_metadata_column, df$con_id)UMAP coloured by condition
condition <- setNames(df$condition, df$con_id)
data <- con$plotGraph(groups=condition)
data <- data$data
ggplot(data, aes(x = x, y =y, color= Group))+ geom_point(size=0.2, alpha=0.5)+ theme_light()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ xlab("") + ylab("") + theme(axis.ticks = element_blank(), axis.text =element_blank())+labs(color = "Condition") +guides(color = guide_legend(override.aes = list(size=3, alpha=1))) +scale_color_manual(values=rainbow(2))
#ggsave("img/2_clustering/fig2.png")
UMAP coloured by doublet scores
doubletscores <- setNames(df$scrublet, df$con_id)
data <- con$plotGraph()
data <- data$data
ggplot(data, aes(x = x, y =y, color=doubletscores))+ geom_point(size=0.2, alpha=0.5)+ theme_light()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ xlab("")+ylab("") + theme(axis.ticks = element_blank(), axis.text =element_blank())+ scale_color_viridis_b(option = "viridis", direction = -1) # change _b to _c for continuous coloring scale
#ggsave("img/2_clustering/fig3.png")
UMAP coloured by UMI counts.
depth <- setNames(df$UMI_counts, df$con_id)
data <- con$plotGraph()
data <- data$data
ggplot(data, aes(x = x, y =y, color=depth))+ geom_point(size=0.2, alpha=0.5)+ theme_light()+ theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ xlab("")+ylab("") + theme(axis.ticks = element_blank(), axis.text =element_blank())+ scale_color_viridis_c(option = "viridis", direction = -1) # change _b to _c for continuous coloring scale
#ggsave("img/2_clustering/fig4.png")
UMAP can be coloured by various metadata to confirm QC filtering was successful and we do not see clusters of low quality cells. Save conos object if satisfied with sample integration. Otherwise adjust alignment strength parameter.
6 Stacked barplots
Plot composition of samples per cluster. We expect that each sample is represented in each cluster and not that one cluster is composed of one/a few samples.
plotClusterBarplots(con, groups = con$clusters$leiden$groups , show.entropy = F, show.size = F) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_discrete(name = "condition") + theme(plot.margin = margin(0.5,0,0,2, "cm"))
#ggsave("img/2_clustering/fig5.png")
Plot composition of
conditions per cluster.
condition <- setNames(df$condition, df$con_id)
plotClusterBarplots(con, groups = con$clusters$leiden$groups , show.entropy = F, show.size = F, sample.factor = condition) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + scale_fill_discrete(name = "condition") + theme(plot.margin = margin(0.5,0,0,2, "cm"))
#ggsave("img/2_clustering/fig6.png")
7 Refine clustering
If needed, the leiden clustering algorithm can be rerun to refine the clusters (change resolution parameter). It is not performed here.
con$findCommunities(method = leiden.community, resolution = 4, min.group.size = 50)Plot new clusters in embedding and also in a table.
con$plotGraph(groups = con$clusters$leiden$groups %>% factor, title="")
table(con$clusters$leiden$groups %>% factor)Get the clustering factor from conos object.
leiden_number_of_clusters <- con$clusters$leiden$groups %>% factor8 Cluster marker genes
Calculate upregulated genes per cluster (= marker genes), use only control cells, because treatment can alter marker gene expression. Use the leiden factor saved above or within the conos object and extract only the ctrl cells.
leiden_xx <- con$clusters$leiden$groups %>% factor
leiden_xx_ctrl <- leiden_xx[grep("control", names(leiden_xx))]
length(leiden_xx_ctrl)Then, get marker genes per cluster using only control cells. Determine differential expressed genes, comparing each group against all others using Wilcoxon rank sum test.
# uses pagoda2's "getDifferentialGenes()"
de_conos <- con$getDifferentialGenes(groups=leiden_xx_ctrl, append.auc = T, upregulated.only = T)Ordering the genes by AUC prioritizes genes that distinguish the respective cluster well from the rest.
de_conos %<>% lapply(function(x) {x %>% arrange(desc(AUC)) %>% head()})head(de_conos$`1`)The following parameters are reported:
- M = log2 fold change
- Z = adjusted Z score, with positive values indicating higher
expression in a given group compare to the rest
- PValue = significance corresponding to the test
- PAdj = adjusted P-value
- AUC = discrimination between cluster and other cells
- Specificity, Precision
- ExpressionFraction = sum expression of the gene within the cluster
divided by the total expression of this gene
Export top 30 genes for users to explore. Transform the list of marker genes into one big table by adding the column “Clusters”. Order by AUC and save csv file with top 30 genes.
library(data.table)
de_conos <- lapply(de_conos, function(x) {x %>% arrange(desc(AUC)) %>% head(30)})
write.xlsx(de_table, file = "top30_marker_byAUC.csv")8.1 Plot marker genes on UMAP
If user gave a list of expected cell types and marker genes those genes can be plotted here.
#plot several genes
genes <- c("Gad1", "Gad2", "Lhx2", "Ascl1")
plots <- lapply(genes, function(gene) {
con$plotGraph(gene = gene, title = gene, size = 0.2) + theme(legend.position = "right", legend.title = element_text("Expression"))+ scale_color_viridis_c(direction = -1)
})
cowplot::plot_grid(plotlist = plots)
#ggsave("img/2_clustering/fig7.png")
9 Annotation
Cell type annotation is manually performed by the user who is the
expert in the respective tissue. CyteTypeR, a multi-agent AI system
where specialized agents collaborate on marker analysis, literature
evidence, and Cell Ontology mapping, can be used as guidance. MapMyCells,
which uses the Allen Institutes taxonomies, is currently tested as an
alternative.
The final annotation is stored as a factor “anno” and is also added to
the metadata table
new_column <- tibble(con_id = names(anno), annotation = anno)
df %<>% left_join(new_column, by = "con_id")10 Useful functions for refining and annotating clusters and subsetting cms
11 Final plots
Final UMAP coloured by cluster It is possible to define custom colours of clusters.
colours <- c( "#FF0000" ,"#FF6000", "#FFBF00", "#DFFF00" ,"#00FFFF" , "#20FF00", "#00FF40", "#DF00FF","#80FF00","#009FFF", "#0040FF","#FF0060" , "#8000FF", "#00FF9F", "#FF00BF","#2000FF" ) %>% setNames(c("Cajal_Retzius_cells", "CGE_progenitors","Cortical_hem" , "Dorsal_IPC" , "Immature_GABAergic_N","Highly_prolif_RG_S" , "Immature_PN", "Immature_PN_Tbr1","MGE_progenitors","Microglia",
"Dorsal_RG" , "Anteriomedial_RG", "Lateral_RG", "PIC_immature_PN", "Roof_plate","Highly_prolif_RG_G2M"))11.1 UMAPs
UMAP without cell type annotation, so user can add cell type names.
con$plotGraph(groups = anno, color.by="cluster",mark.groups=F, size=0.2) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ scale_colour_manual(values=colours)
#ggsave("img/2_clustering/fig8.png")
UMAP with cell type names
con$plotGraph(groups=anno, color.by="cluster",mark.groups=T, size=0.2) + theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ scale_colour_manual(values=colours)
#ggsave("img/2_clustering/fig9.png")