This tutorial demonstrates how to perform vertical integrative clustering (iCluster) analysis of single cell multi-omics data. The iCluster2b function in the iClusterPlus package will be used as an example. The iCluster2b function implements an updated iCluster method originally reported by Shen, Wang and Mo. 2023. The CITE-seq dataset (Stuart, Butler et al., Cell 2019) consisting of scRNA-seq and ADT (antibody-derived tag) data of 30,672 cells from human bone marrow is used for the demonstration. The dataset can be obtained by installing the R package SeuratData.

We use the R scripts of the Seurat’s tutorial Weighted Nearest Neighbor Analysis to process the scRNA-seq and ADT data. The iCluster analysis of single cell multi-omics data can be carried out in 4 steps as described in the following.

library(iClusterPlus)
library(lattice)
library(gplots)
library(Seurat)
library(SeuratData)
library(cowplot)
library(dplyr)
library(Rphenograph)
library(umap)
library(mclust)
library(Silhouette)
library(cluster)
library(knitr)

Step 1: Processing data using the Seurat package

The scRNA-seq and ADT data are processed separately, which includes normalization, feature selection, scaling and dimension reduction.

InstallData("bmcite")
bm <- LoadData(ds = "bmcite")
DefaultAssay(bm) <- 'RNA'
bm <- NormalizeData(bm) %>% FindVariableFeatures() %>% ScaleData() %>% RunPCA()

# process ADT data; all the 25 features are used for dimensional reduction
DefaultAssay(bm) <- 'ADT'
VariableFeatures(bm) <- rownames(bm[["ADT"]])
bm <- NormalizeData(bm, normalization.method = 'CLR', margin = 2) %>% ScaleData() %>% RunPCA(reduction.name = 'apca')

Step 2: Perform iCluster modeling

The iCluster modeling generates a latent factor matrix (meanZ), which captures the major sources of variation of the multi-omics data (scRNA-seq and ADT). As a result, the latent factor matrix can be used for sample clustering and visualization.

# input data for iCluster2b 
BMlist <- list(ADT=t(bm[['ADT']]$scale.data),RNA=t(bm[['RNA']]$scale.data))

# PCA variance plot 
pcaVarPlot(BMlist,K=50)

From the PCA variance plot, we can see that the explained variance starts to level off around the 10th PC. For iCluster modeling of single cell multi-omics data, we suggest selecting the number PCs based on the point at which the explained variance levels off. In general, using the first 10 to 20 PCs is sufficient to distinguish major cell types. Using higher-numbered PCs (e.g., 40-50) may help identify rare cell populations and capture more subtle heterogeneity. As an example, we select K=20 as the number of latent factors for iCluster modeling. It can be seen that the clustering results are quite consistent for a range of numbers (e.g., K = 20 - 30).

# since features have been selected, we set lambda to a small value to minimize feature selection
date()
## [1] "Thu Nov  5 15:27:43 2026"
bm.ic = iCluster2b(xList=BMlist,K=20,lambda=list(0.1,0.1),method=c("lasso","lasso"))
date()
## [1] "Thu Nov  5 15:28:00 2026"

Step 3: Clustering of cells using the latent factor matrix

We use the Louvain method implemented in the Rphenograph package to perform cell clustering analysis. Then, we project the cells to a 2-dimensional space using UMAP for visualization.

## clustering of cells using Louvain method and the latent factor matrix generated by iCluster
ic.graph <- Rphenograph(bm.ic$meanZ, k=30)
##   Finding nearest neighbors...DONE ~ 4.312 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 5.182 s
##   Build undirected graph from the weighted links...DONE ~ 0.88 s
##   Run louvain clustering on the graph ...DONE ~ 1.315 s
##   Return a community class
##   -Modularity value: 0.9018461 
##   -Number of clusters: 25
tabclst <- table(membership(ic.graph[[2]]))

## order the clusters based on cluster size
clstlevels <- as.numeric(names(tabclst)[order(tabclst,decreasing = TRUE)])
bm.iClusters <- factor(membership(ic.graph[[2]]),levels=clstlevels,labels=1:length(clstlevels))
names(bm.iClusters) <- rownames(BMlist$ADT)

# Projecting cells to 2-dimensional space using UMAP
ic.umap <- umap(bm.ic$meanZ)
rownames(ic.umap$layout) <- rownames(BMlist$ADT)
colnames(ic.umap$layout) <- c("umap1","umap2")

## put umap layout and cluster IDs into the Seurat object bm 
bm[["ic.umap"]] <- CreateDimReducObject(embeddings = ic.umap$layout, key = "ic_umap")
bm$iCluster <- bm.iClusters

For comparison, we perform unsupervised clustering analysis using Seurat’s WNN method. For more details of WNN analysis, we refer users to the Seurat’s tutorial (Weighted Nearest Neighbor Analysis).

# Generate WNN graph 
bm <- FindMultiModalNeighbors(bm, reduction.list = list("pca", "apca"),dims.list <- list(1:30, 1:18), k.nn = 20,
   knn.graph.name = "wknn",snn.graph.name = "wsnn",modality.weight.name = "RNA.weight")
# multimodal neighbor matrix used for UMAP can be accessed using bm[['weighted.nn']]
# WNN graph can be accessed using bm[["wknn"]]
# SNN graph used for clustering can be accessed at bm[["wsnn"]]

bm <- RunUMAP(bm, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
bm <- FindClusters(bm, graph.name = "wsnn", algorithm = 3, resolution = 2, verbose = FALSE)

Step 4: Visualization of cell clusters

We use Seurat’s function DimPlot to visualize the cells on a 2-dimensional UMAP space.

  1. Visualization of cell clusters generated by iCluster analysis, alongside cell annotations.
p1 <- DimPlot(bm, reduction = "ic.umap", group.by = "iCluster", label = TRUE, label.size = 2.5, repel = TRUE) + NoLegend() + ggtitle("iClusters") + labs(x="UMAP_1",y="UMAP_2")
p2 <- DimPlot(bm, reduction = "ic.umap", group.by = "celltype.l2", label = TRUE, label.size = 2.5, repel = TRUE) + NoLegend() + ggtitle("Cell types") + labs(x="UMAP_1",y="UMAP_2")
p1 + p2

  1. Heatmap of cell clusters generated by iCluster analysis. We show the heatmap of cell clusters generated by iCluster analysis using the top 25% features. The ranking of ADT and RNA features are based on the sum of the absolute value of the parameter (beta) in a decreasing order.
col.scheme <- alist()
col.scheme[[1]] <- bluered(256)
col.scheme[[2]] <- bluered(256)

bm.ic$clusters <- bm.iClusters
ic.hm <- plotHeatmap(fit=bm.ic,xList=BMlist,type=c("gaussian","gaussian"), 
            threshold=c(0.25, 0.25),feature.order=c(T,T),feature.scale=c(F,F),
            dist.method = "euclidean", hclust.method="ward.D",col.scheme = col.scheme,width=10,
            sparse=c(F,T),cap=c(0.99,0.99))

  1. Visualization of cell clusters generated by Seurat WNN analysis, alongside cell annotations
p3 <- DimPlot(bm, reduction = "wnn.umap", group.by = "seurat_clusters", label = TRUE, label.size = 2.5, repel = TRUE) + NoLegend() + ggtitle("Seurat")
p4 <- DimPlot(bm, reduction = "wnn.umap", group.by = "celltype.l2", label = TRUE, label.size = 2.5, repel = TRUE) + NoLegend() + ggtitle("celltype")
p3 + p4

  1. Heatmap of cell clusters generated by Seurat WNN analysis. The ranking of ADT and RNA top 25% features are based on the sum of the absolute value of the iCluster parameter (beta) in a decreasing order.
bm.ic$clusters <- bm$seurat_clusters
wnn.hm <- plotHeatmap(fit=bm.ic,xList=BMlist,type=c("gaussian","gaussian"), 
            threshold=rep(0.25, 2),feature.order=c(T,T),feature.scale=c(F,F),
            dist.method = "euclidean", hclust.method="ward.D",col.scheme = col.scheme,width=10,
            sparse=c(F,T),cap=c(0.99,0.99))

Finally, we compare the integrative clustering results of iCluster and Seurat WNN using normalized mutual information (NMI), adjusted rand index (ARI) and average silhouette width (ASW). In general, higher index values indicate better clustering performance.

compare(bm$celltype.l2, bm$iCluster, method = "nmi")
## [1] 0.8464881
compare(bm$celltype.l2, bm$seurat_clusters, method = "nmi")
## [1] 0.8841955
# Comparison using adjusted rand index 
compare(bm$celltype.l2, bm$iCluster, method = "adjusted.rand")
## [1] 0.7788463
compare(bm$celltype.l2, bm$seurat_clusters, method = "adjusted.rand")
## [1] 0.718272
# Comparison using average silhouette width
sil.ic <- silhouette(as.numeric(bm$iCluster), dist(Embeddings(bm, reduction = "ic.umap")))
summary(sil.ic)$avg.width  
## [1] 0.4790479
sil.wnn <- silhouette(as.numeric(bm$seurat_clusters), dist(Embeddings(bm, reduction = "wnn.umap")))
summary(sil.wnn)$avg.width 
## [1] 0.2809133

#Table of predicted cell types vs. iClusters

icTable <- table(bm$celltype.l2,bm.iClusters)
knitr::kable(icTable)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
CD14 Mono 0 4127 0 0 2232 0 1 0 0 1 12 0 0 4 6 2 5 0 0 3 0 28 0 40 25
CD16 Mono 0 0 0 0 12 0 0 0 0 0 0 0 0 0 420 0 0 0 0 0 0 1 0 0 0
CD4 Memory 221 0 0 3085 0 0 2 1 0 5 0 0 0 0 0 0 5 0 0 0 0 22 19 0 0
CD4 Naive 4345 1 1 56 0 0 0 1 0 1 0 0 0 0 0 0 1 0 1 0 0 73 20 0 0
CD56 bright NK 3 1 0 3 0 0 9 0 103 1 0 0 0 0 0 0 1 0 0 0 0 2 20 0 0
CD8 Effector_1 0 0 2 0 0 0 552 0 2 14 0 0 0 0 0 0 0 0 0 0 0 3 4 0 0
CD8 Effector_2 0 1 3 1 1 0 294 0 0 5 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0
CD8 Memory_1 0 0 10 0 0 0 406 0 0 21 0 0 0 0 0 0 0 0 0 0 0 0 7 0 0
CD8 Memory_2 0 1 1 1 0 0 503 0 1 37 1 0 0 1 0 0 3 0 0 0 0 6 0 0 0
CD8 Naive 1 1 3889 0 1 0 62 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 20 0 0
cDC2 0 13 0 0 120 0 0 0 0 0 0 0 1 347 0 0 0 0 0 0 0 1 0 0 0
gdT 0 0 0 2 0 0 13 0 0 349 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0
GMP 0 23 0 0 4 0 0 0 0 0 573 8 0 2 0 0 136 0 0 2 0 0 0 0 0
HSC 1 0 0 0 0 0 0 0 0 1 0 247 0 0 0 12 0 0 0 0 0 2 0 0 0
LMPP 1 0 0 0 0 0 0 0 0 0 15 264 0 1 0 8 2 0 1 0 0 0 0 0 0
MAIT 0 1 0 1 0 0 26 0 0 488 0 0 0 0 0 0 0 0 0 0 0 1 3 0 0
Memory B 0 0 1 0 0 229 0 1389 1 1 0 0 0 0 0 0 0 0 0 0 0 8 1 0 0
Naive B 0 0 0 0 0 1756 0 117 0 0 0 0 0 0 0 0 0 0 8 0 0 15 0 0 4
NK 0 1 1 0 0 0 6 0 1237 8 1 0 0 0 0 0 8 0 0 0 0 4 1 0 0
pDC 0 0 0 0 0 0 0 0 0 0 0 3 0 6 0 1 0 318 0 0 0 0 0 0 0
Plasmablast 0 1 0 0 0 0 0 7 0 0 0 0 0 0 0 0 1 0 0 1 211 1 1 0 0
Prog_B 1 0 0 0 0 0 0 0 1 0 0 0 31 0 0 0 0 0 0 111 1 0 0 1 0 0
Prog_B 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 4 0 125 0 0 0 0 0 0
Prog_DC 0 2 0 0 6 0 0 0 0 0 4 2 0 91 0 0 155 2 0 1 0 0 0 0 0
Prog_Mk 0 0 0 0 0 0 0 0 0 0 0 2 0 0 0 119 0 0 0 0 0 3 0 0 0
Prog_RBC 1 0 0 0 2 0 2 1 0 0 0 1 468 0 0 223 1 0 0 213 0 1 2 0 0
Treg 163 0 1 126 0 0 1 0 0 1 0 0 0 0 0 0 2 0 0 0 0 3 0 0 0
Session Info 
sessionInfo()
## R version 4.5.2 (2025-10-31)
## Platform: aarch64-apple-darwin20
## Running under: macOS Sequoia 15.8.1
## 
## Matrix products: default
## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] parallel  stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] future_1.75.0            knitr_1.51               cluster_2.1.8.3         
##  [4] Silhouette_0.9.6         mclust_6.1.3             umap_0.2.10.0           
##  [7] Rphenograph_0.99.1       igraph_2.3.3             ggplot2_4.0.3           
## [10] dplyr_1.2.1              cowplot_1.2.0            pbmcsca.SeuratData_3.0.0
## [13] pbmcref.SeuratData_1.0.0 bmcite.SeuratData_0.3.0  SeuratData_0.2.2.9002   
## [16] Seurat_5.5.1             SeuratObject_5.4.0       sp_2.2-3                
## [19] gplots_3.3.0             lattice_0.22-9           iClusterPlus_1.49.1     
## [22] irlba_2.3.7              Matrix_1.7-6            
## 
## loaded via a namespace (and not attached):
##   [1] RColorBrewer_1.1-3     rstudioapi_0.19.0      jsonlite_2.0.0        
##   [4] magrittr_2.0.5         spatstat.utils_3.2-4   farver_2.1.2          
##   [7] rmarkdown_2.31         vctrs_0.7.3            ROCR_1.0-12           
##  [10] spatstat.explore_3.8-2 askpass_1.2.1          rstatix_1.1.0         
##  [13] htmltools_0.5.9        broom_1.0.13           Formula_1.2-6         
##  [16] sass_0.4.10            sctransform_0.4.3      parallelly_1.48.0     
##  [19] KernSmooth_2.23-26     bslib_0.12.0           htmlwidgets_1.6.4     
##  [22] ica_1.0-3              plyr_1.8.9             plotly_4.12.1         
##  [25] zoo_1.9-0              cachem_1.1.0           mime_0.13             
##  [28] lifecycle_1.0.5        pkgconfig_2.0.3        R6_2.6.1              
##  [31] fastmap_1.2.0          fitdistrplus_1.2-6     shiny_1.14.0          
##  [34] digest_0.6.39          patchwork_1.3.2        tensor_1.5.1          
##  [37] RSpectra_0.16-2        ggpubr_1.0.0           labeling_0.4.3        
##  [40] progressr_1.0.0        spatstat.sparse_3.2-0  httr_1.4.8            
##  [43] polyclip_1.10-7        abind_1.4-8            compiler_4.5.2        
##  [46] withr_3.0.3            backports_1.5.1        S7_0.2.2              
##  [49] carData_3.0-6          fastDummies_1.7.6      ggsignif_0.6.4        
##  [52] MASS_7.3-66            openssl_2.4.2          rappdirs_0.3.4        
##  [55] gtools_3.9.5           caTools_1.18.4         tools_4.5.2           
##  [58] lmtest_0.9-40          otel_0.2.0             httpuv_1.6.17         
##  [61] future.apply_1.20.2    goftest_1.2-3          glue_1.8.1            
##  [64] nlme_3.1-170           promises_1.5.0         grid_4.5.2            
##  [67] Rtsne_0.17             reshape2_1.4.5         generics_0.1.4        
##  [70] gtable_0.3.6           spatstat.data_3.1-9    tidyr_1.3.2           
##  [73] data.table_1.18.4      car_3.1-5              spatstat.geom_3.8-2   
##  [76] RcppAnnoy_0.0.23       ggrepel_0.9.8          RANN_2.6.2            
##  [79] pillar_1.11.1          stringr_1.6.0          spam_2.11-4           
##  [82] RcppHNSW_0.7.0         later_1.4.8            splines_4.5.2         
##  [85] survival_3.8-9         deldir_2.0-4           tidyselect_1.2.1      
##  [88] miniUI_0.1.2           pbapply_1.7-4          gridExtra_2.3.1       
##  [91] scattermore_1.2        xfun_0.60              matrixStats_1.5.0     
##  [94] stringi_1.8.9          yaml_2.3.12            evaluate_1.0.5        
##  [97] codetools_0.2-20       tibble_3.3.1           cli_3.6.6             
## [100] uwot_0.2.4             xtable_1.8-8           reticulate_1.46.0     
## [103] jquerylib_0.1.4        dichromat_2.0-1        Rcpp_1.1.2            
## [106] globals_0.19.1         spatstat.random_3.5-1  png_0.1-9             
## [109] spatstat.univar_3.2-0  dotCall64_1.2          bitops_1.1-0          
## [112] listenv_1.0.0          viridisLite_0.4.3      scales_1.4.0          
## [115] ggridges_0.5.7         purrr_1.2.2            crayon_1.5.3          
## [118] rlang_1.3.0