In this tutorial, we demonstrate how to perform vertical integrative clustering (iCluster) analysis of single cell RNA-seq and ATAC-seq data using the iCluster2b function, which implements an updated iCluster method originally reported by Shen, Wang and Mo. 2023. We use a paired scRNA-seq and scATAC-seq data measured for 10,412 PBMCs for iCluster analysis. The paired data are publicly available on the 10x website, which include the following files.

The paired data are used by Seurat to illustrate the Weighted Nearest Neighbor (WNN) method for integrative clustering of single cell multi-omics data. We will follow the tutorial WNN Analysis to process the paired scRNA-seq and scATAC-seq data.

library(iClusterPlus)
library(lattice)
library(gplots)
library(Seurat)
library(SeuratData)
library(cowplot)
library(dplyr)
library(Rphenograph)
library(umap)
library(mclust)
library(MatrixGenerics)
library(Matrix)
library(EnsDb.Hsapiens.v86)
library(ggplot2)
library(GenomeInfoDb) #standardChromosomes
library(Signac) #v1.16.0
library(Azimuth)
library(Silhouette)
library(cluster)
library(knitr)

The iCluster analysis of single cell multi-omics data can be carried out in the following 4 steps.

Step 1: Processing data using the Seurat package

We follow the Seurat’s WNN Analysis tutorial to process the scRNA-seq and scATAC-seq data. The processing workflow includes filtering, normalization, feature selection, scaling and dimension reduction for each dataset.

# The following are Seurat's codes to process the scRNA-seq and scATACT-seq data
# read scRNA-seq and scATAC-seq data from the data folder of working directory
inputdata.10x <- Read10X_h5("../pbmc_data/pbmc_granulocyte_sorted_10k_filtered_feature_bc_matrix.h5")

# extract RNA and ATAC data
rna_counts <- inputdata.10x$`Gene Expression`
atac_counts <- inputdata.10x$Peaks

# Create Seurat object
pbmc <- CreateSeuratObject(counts = rna_counts)
pbmc[["percent.mt"]] <- PercentageFeatureSet(pbmc, pattern = "^MT-")

# Add in the ATAC-seq data and only use peaks in standard chromosomes
grange.counts <- StringToGRanges(rownames(atac_counts), sep = c(":", "-"))
grange.use <- seqnames(grange.counts) %in% standardChromosomes(grange.counts)
atac_counts <- atac_counts[as.vector(grange.use), ]
annotations <- GetGRangesFromEnsDb(ensdb = EnsDb.Hsapiens.v86)
seqlevelsStyle(annotations) <- 'UCSC'
genome(annotations) <- "hg38"

# read ATAC Per fragment information file in the data folder of working directory
frag.file <- "../pbmc_data/pbmc_granulocyte_sorted_10k_atac_fragments.tsv.gz"

# create chromosome assay
chrom_assay <- CreateChromatinAssay(
   counts = atac_counts,
   sep = c(":", "-"),
   genome = 'hg38',
   fragments = frag.file,
   min.cells = 10,
   annotation = annotations
 )
pbmc[["ATAC"]] <- chrom_assay

#basic QC based on the number of detected molecules for RNA, ATAC and mitochondrial percentage.
#VlnPlot(pbmc,features=c("nCount_RNA", "nCount_ATAC", "percent.mt"), ncol=3,log=TRUE,pt.size = 0) + NoLegend()

pbmc <- subset(
  x = pbmc,
  subset = nCount_ATAC < 7e4 &
    nCount_ATAC > 5e3 &
    nCount_RNA < 25000 &
    nCount_RNA > 1000 &
    percent.mt < 20
)

# RNA analysis 
DefaultAssay(pbmc) <- "RNA"
options(future.globals.maxSize= 850*1024^2) #calculate the size for an 850mb limit:
pbmc <- SCTransform(pbmc, verbose = FALSE) %>% RunPCA() %>% RunUMAP(dims = 1:50, reduction.name = 'umap.rna', reduction.key = 'rnaUMAP_')

# ATAC analysis
# Seurat excludes the first dimension as this is typically correlated with sequencing depth
DefaultAssay(pbmc) <- "ATAC"
pbmc <- RunTFIDF(pbmc)
pbmc <- FindTopFeatures(pbmc, min.cutoff = 'q0')
pbmc <- RunSVD(pbmc) #reduction.name = "lsi
pbmc <- RunUMAP(pbmc, reduction = 'lsi', dims = 2:50, reduction.name = "umap.atac", reduction.key = "atacUMAP_")

# cell annotation using Azimuth package
DefaultAssay(pbmc) <- "RNA"
pbmc <- RunAzimuth(pbmc, reference = "pbmcref")

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 scATAC-seq). As a result, the latent factor matrix can be used for sample clustering and visualization.

# Select the top 20% most variable regions and merge the neighboring regions with gap width <=10000 bps
atac <- mergeRegion(peak.ct=pbmc[["ATAC"]]$data,rm.prob=0.8,rm.isolate=FALSE,min.gapwidth=10000)
dim(atac)
## [1] 15441 10412
# we use the normalized and scaled RNA-seq data (top 3000 genes) and top 10% most variable ATAC-seq features for iCluster analysis
dim(pbmc[['SCT']]$scale.data)
## [1]  3000 10412
pbmcList <- list(RNA=t(pbmc[['SCT']]$scale.data),ATAC=t(as.matrix(atac)))
rm(atac)

# PCA variance plot 
pcaVarPlot(pbmcList,K=100)

From the 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. Here, we select K=50 as the number of latent factors for iCluster modeling.

# since features have been selected, we set lambda to a small value to minimize feature selection
date()
## [1] "Thu Nov  5 23:15:05 2026"
fit2 = iCluster2b(xList=pbmcList,K=50,lambda=list(0.1,0.1),method=c("lasso","lasso"))
date()
## [1] "Thu Nov  5 23:16:32 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(fit2$meanZ, k=30)
##   Finding nearest neighbors...DONE ~ 3.395 s
##   Compute jaccard coefficient between nearest-neighbor sets...DONE ~ 1.392 s
##   Build undirected graph from the weighted links...DONE ~ 0.278 s
##   Run louvain clustering on the graph ...DONE ~ 0.414 s
##   Return a community class
##   -Modularity value: 0.8725663 
##   -Number of clusters: 20
tabclst <- table(membership(ic.graph[[2]]))

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

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

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

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
ic.p1 <- DimPlot(pbmc,reduction="ic.umap",group.by="iClusters",label=TRUE,label.size=4,repel=TRUE) + 
  NoLegend() + ggtitle("iClusters") + labs(x="UMAP_1",y="UMAP_2")

ic.p2 <- DimPlot(pbmc, reduction="ic.umap",group.by="predicted.celltype.l2",label=TRUE,label.size=4, repel=TRUE) + 
  NoLegend() + ggtitle("Cell types") + labs(x="UMAP_1",y="UMAP_2")

ic.p1 + ic.p2

  1. Heatmap of cell clusters generated by iCluster analysis. We show the heatmap of cell clusters generated by iCluster analysis using the top 10% of the RNA and ATAC features. The ranking of RNA and ATAC 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)

fit2$clusters <- pbmc$iClusters
ic.hm <- iClusterPlus::plotHeatmap(fit=fit2,xList=pbmcList,type=c("gaussian","gaussian"), 
                     threshold=c(0.10, 0.10),feature.order=c(T,T),feature.scale=c(T,T),
                     dist.method = "euclidean", hclust.method="ward.D",col.scheme = col.scheme,width=10,
                     sparse=c(T,T),cap=c(0.99,0.99))

For comparison, we perform integrative 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).

pbmc <- FindMultiModalNeighbors(pbmc, reduction.list = list("pca", "lsi"), dims.list = list(1:50, 2:50))
pbmc <- RunUMAP(pbmc, nn.name = "weighted.nn", reduction.name = "wnn.umap", reduction.key = "wnnUMAP_")
pbmc <- FindClusters(pbmc, graph.name = "wsnn", algorithm = 3, verbose = FALSE)

We visualize the cell clusters generated by Seurat WNN analysis, alongside cell annotations

p3 <- DimPlot(pbmc, reduction="wnn.umap",group.by="seurat_clusters", label=TRUE,label.size=4,repel=TRUE) + 
  NoLegend() + ggtitle("Seurat_clusters")

p4 <- DimPlot(pbmc, reduction="wnn.umap",group.by="predicted.celltype.l2",label=TRUE,label.size=4,repel=TRUE) + 
  NoLegend() + ggtitle("Cell types")

p3 + p4

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.

# comparison using normalized mutual information (NMI)
compare(pbmc$predicted.celltype.l2, pbmc$iClusters, method = "nmi") 
## [1] 0.7623989
compare(pbmc$predicted.celltype.l2, pbmc$seurat_clusters, method = "nmi") 
## [1] 0.7315697
# Comparison using adjusted rand index 
compare(pbmc$predicted.celltype.l2, pbmc$iClusters, method = "adjusted.rand") 
## [1] 0.6155965
compare(pbmc$predicted.celltype.l2, pbmc$seurat_clusters, method = "adjusted.rand") 
## [1] 0.596371
# Comparison using average silhouette width
sil.ic <- silhouette(as.numeric(pbmc$iClusters), dist(Embeddings(pbmc, reduction = "ic.umap")))
summary(sil.ic)$avg.width 
## [1] 0.428286
sil.wnn <- silhouette(as.numeric(pbmc$seurat_clusters), dist(Embeddings(pbmc, reduction = "wnn.umap")))
summary(sil.wnn)$avg.width
## [1] 0.4177851

Table of cell types vs. iClusters

# table of annotated cell types and iClusters 
icTable <- table(pbmc$predicted.celltype.l2,pbmc$iClusters)
knitr::kable(icTable)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20
B intermediate 1 1 3 1 0 0 346 0 0 0 25 10 2 0 0 0 0 0 0 0
B memory 0 1 0 0 0 0 105 0 0 0 0 1 0 0 0 0 0 0 0 0
B naive 2 0 1 0 0 0 14 0 0 0 318 16 0 0 0 0 0 0 0 0
CD14 Mono 2123 3 0 2 1 559 1 1 31 0 0 280 46 1 1 0 0 0 1 0
CD16 Mono 11 0 0 0 0 2 0 0 370 0 0 3 0 0 0 0 0 0 0 0
CD4 Naive 0 731 0 18 0 0 0 0 0 9 0 1 0 1 0 0 0 0 0 0
CD4 Proliferating 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0
CD4 TCM 0 660 1343 47 8 0 0 1 0 8 1 7 1 40 2 0 4 0 1 0
CD4 TEM 0 0 10 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
CD8 Naive 0 20 1 1000 17 0 0 2 0 358 1 0 0 0 0 0 3 0 0 0
CD8 TCM 0 2 23 11 87 0 0 0 0 4 0 1 0 0 1 0 0 0 0 0
CD8 TEM 0 1 5 1 530 0 0 19 0 1 0 7 0 2 1 0 32 0 0 0
cDC1 0 0 0 0 0 0 0 0 0 0 0 0 6 0 0 0 0 0 0 0
cDC2 0 0 0 0 0 0 0 0 0 0 0 0 124 0 0 0 0 0 0 0
dnT 0 0 6 2 2 0 0 0 0 2 0 0 0 0 0 0 0 0 1 0
gdT 0 4 0 1 69 0 0 9 0 1 0 1 0 0 5 0 12 0 0 0
HSPC 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 25 0 0
ILC 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0
MAIT 0 0 5 1 1 0 0 1 0 0 1 1 0 0 105 0 0 0 0 0
NK 1 0 0 0 36 0 2 380 0 0 0 7 0 0 2 0 0 0 0 0
NK Proliferating 0 0 0 0 1 0 0 3 0 0 0 0 0 0 0 0 0 0 10 0
NK_CD56bright 0 0 0 0 0 0 0 22 0 0 0 0 0 0 0 0 0 0 0 0
pDC 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 106 0 0 0 0
Plasmablast 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 18
Treg 0 3 1 1 0 0 0 0 0 0 0 1 0 116 0 0 1 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] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] future_1.75.0             knitr_1.51               
##  [3] cluster_2.1.8.3           Silhouette_0.9.6         
##  [5] Azimuth_0.5.0             shinyBS_0.65.0           
##  [7] Signac_1.16.0             GenomeInfoDb_1.46.2      
##  [9] EnsDb.Hsapiens.v86_2.99.0 ensembldb_2.34.0         
## [11] AnnotationFilter_1.34.0   GenomicFeatures_1.62.0   
## [13] AnnotationDbi_1.72.0      Biobase_2.70.0           
## [15] GenomicRanges_1.62.1      Seqinfo_1.0.0            
## [17] IRanges_2.44.0            S4Vectors_0.48.1         
## [19] BiocGenerics_0.56.0       generics_0.1.4           
## [21] MatrixGenerics_1.22.0     matrixStats_1.5.0        
## [23] mclust_6.1.3              umap_0.2.10.0            
## [25] Rphenograph_0.99.1        igraph_2.3.3             
## [27] ggplot2_4.0.3             dplyr_1.2.1              
## [29] cowplot_1.2.0             pbmcsca.SeuratData_3.0.0 
## [31] pbmcref.SeuratData_1.0.0  bmcite.SeuratData_0.3.0  
## [33] SeuratData_0.2.2.9002     Seurat_5.5.1             
## [35] SeuratObject_5.4.0        sp_2.2-3                 
## [37] gplots_3.3.0              lattice_0.22-9           
## [39] iClusterPlus_1.49.1       irlba_2.3.7              
## [41] Matrix_1.7-6             
## 
## loaded via a namespace (and not attached):
##   [1] fs_2.1.0                          ProtGenerics_1.42.0              
##   [3] spatstat.sparse_3.2-0             bitops_1.1-0                     
##   [5] DirichletMultinomial_1.52.0       TFBSTools_1.48.0                 
##   [7] httr_1.4.8                        RColorBrewer_1.1-3               
##   [9] backports_1.5.1                   tools_4.5.2                      
##  [11] sctransform_0.4.3                 R6_2.6.1                         
##  [13] DT_0.34.0                         lazyeval_0.2.3                   
##  [15] uwot_0.2.4                        rhdf5filters_1.22.0              
##  [17] withr_3.0.3                       gridExtra_2.3.1                  
##  [19] progressr_1.0.0                   cli_3.6.6                        
##  [21] spatstat.explore_3.8-2            fastDummies_1.7.6                
##  [23] shinyjs_2.1.1                     labeling_0.4.3                   
##  [25] sass_0.4.10                       S7_0.2.2                         
##  [27] spatstat.data_3.1-9               ggridges_0.5.7                   
##  [29] pbapply_1.7-4                     askpass_1.2.1                    
##  [31] Rsamtools_2.26.0                  foreign_0.8-91                   
##  [33] dichromat_2.0-1                   parallelly_1.48.0                
##  [35] BSgenome_1.78.0                   rstudioapi_0.19.0                
##  [37] RSQLite_3.53.3                    BiocIO_1.20.0                    
##  [39] gtools_3.9.5                      ica_1.0-3                        
##  [41] spatstat.random_3.5-1             car_3.1-5                        
##  [43] googlesheets4_1.1.2               abind_1.4-8                      
##  [45] lifecycle_1.0.5                   yaml_2.3.12                      
##  [47] carData_3.0-6                     SummarizedExperiment_1.40.0      
##  [49] glmGamPoi_1.22.0                  rhdf5_2.54.1                     
##  [51] SparseArray_1.10.10               Rtsne_0.17                       
##  [53] grid_4.5.2                        blob_1.3.0                       
##  [55] promises_1.5.0                    shinydashboard_0.7.3             
##  [57] crayon_1.5.3                      pwalign_1.6.0                    
##  [59] miniUI_0.1.2                      beachmat_2.26.0                  
##  [61] cigarillo_1.0.0                   KEGGREST_1.50.0                  
##  [63] pillar_1.11.1                     rjson_0.2.23                     
##  [65] future.apply_1.20.2               codetools_0.2-20                 
##  [67] fastmatch_1.1-8                   glue_1.8.1                       
##  [69] spatstat.univar_3.2-0             data.table_1.18.4                
##  [71] vctrs_0.7.3                       png_0.1-9                        
##  [73] spam_2.11-4                       cellranger_1.1.0                 
##  [75] gtable_0.3.6                      cachem_1.1.0                     
##  [77] xfun_0.60                         S4Arrays_1.10.1                  
##  [79] mime_0.13                         survival_3.8-9                   
##  [81] gargle_1.6.1                      RcppRoll_0.3.2                   
##  [83] fitdistrplus_1.2-6                ROCR_1.0-12                      
##  [85] nlme_3.1-170                      bit64_4.8.2                      
##  [87] RcppAnnoy_0.0.23                  bslib_0.12.0                     
##  [89] rpart_4.1.27                      KernSmooth_2.23-26               
##  [91] otel_0.2.0                        colorspace_2.1-3                 
##  [93] Hmisc_5.2-6                       SeuratDisk_0.0.0.9021            
##  [95] seqLogo_1.76.0                    DBI_1.3.0                        
##  [97] nnet_7.3-21                       tidyselect_1.2.1                 
##  [99] bit_4.6.0                         compiler_4.5.2                   
## [101] curl_7.1.0                        htmlTable_2.5.0                  
## [103] hdf5r_1.3.12                      DelayedArray_0.36.1              
## [105] plotly_4.12.1                     rtracklayer_1.70.1               
## [107] checkmate_2.3.4                   scales_1.4.0                     
## [109] caTools_1.18.4                    lmtest_0.9-40                    
## [111] rappdirs_0.3.4                    stringr_1.6.0                    
## [113] digest_0.6.39                     goftest_1.2-3                    
## [115] presto_1.0.0                      spatstat.utils_3.2-4             
## [117] rmarkdown_2.31                    XVector_0.50.0                   
## [119] base64enc_0.1-6                   htmltools_0.5.9                  
## [121] pkgconfig_2.0.3                   sparseMatrixStats_1.22.0         
## [123] fastmap_1.2.0                     rlang_1.3.0                      
## [125] htmlwidgets_1.6.4                 UCSC.utils_1.6.1                 
## [127] DelayedMatrixStats_1.32.0         shiny_1.14.0                     
## [129] farver_2.1.2                      jquerylib_0.1.4                  
## [131] zoo_1.9-0                         jsonlite_2.0.0                   
## [133] BiocParallel_1.44.0               VariantAnnotation_1.56.0         
## [135] RCurl_1.98-1.19                   magrittr_2.0.5                   
## [137] Formula_1.2-6                     dotCall64_1.2                    
## [139] patchwork_1.3.2                   Rhdf5lib_1.32.0                  
## [141] Rcpp_1.1.2                        reticulate_1.46.0                
## [143] stringi_1.8.9                     MASS_7.3-66                      
## [145] plyr_1.8.9                        listenv_1.0.0                    
## [147] ggrepel_0.9.8                     deldir_2.0-4                     
## [149] Biostrings_2.78.0                 splines_4.5.2                    
## [151] tensor_1.5.1                      BSgenome.Hsapiens.UCSC.hg38_1.4.5
## [153] ggpubr_1.0.0                      spatstat.geom_3.8-2              
## [155] ggsignif_0.6.4                    RcppHNSW_0.7.0                   
## [157] reshape2_1.4.5                    TFMPvalue_1.0.0                  
## [159] XML_3.99-0.23                     evaluate_1.0.5                   
## [161] biovizBase_1.58.0                 JASPAR2020_0.99.10               
## [163] httpuv_1.6.17                     RANN_2.6.2                       
## [165] tidyr_1.3.2                       openssl_2.4.2                    
## [167] purrr_1.2.2                       polyclip_1.10-7                  
## [169] scattermore_1.2                   broom_1.0.13                     
## [171] xtable_1.8-8                      restfulr_0.0.17                  
## [173] RSpectra_0.16-2                   rstatix_1.1.0                    
## [175] later_1.4.8                       googledrive_2.1.2                
## [177] viridisLite_0.4.3                 tibble_3.3.1                     
## [179] memoise_2.0.1                     GenomicAlignments_1.46.0         
## [181] globals_0.19.1