SpatialLeiden - Spatially-aware Leiden clustering

13 Clustering can identify the natural structure that is inherent to measured data. For single-cell 14 omics, clustering finds cells with similar molecular phenotype after which cell types are 15 annotated. Leiden clustering is the algorithm of choice in the single-cell community. 16 However, in the field of spatial omics, Leiden has been considered a non-spatial clustering 17 method. Here, we show that by integrating spatial embeddings Leiden clustering is rendered 18 into a computationally highly performant, spatially aware clustering method that compares 19 well with state-of-the art clustering methods. 20

21 Spatial omics, clustering, Leiden, domains, niches, spatial clustering, spatial biology, 22 bioinformatics 23

24 Single-cell transcriptomics has revolutionised our understanding of cellular heterogeneity by 25 enabling the measurement of gene expression at the individual cell level. However, this high-26 dimensional data poses significant challenges in extracting meaningful biological insights. 27 This can be overcome by grouping cells with similar expression profiles into distinct clusters. 28 By partitioning cells based on transcriptional similarities, clustering facilitates the 29 characterization of cell-type diversity within a heterogeneous cell population. Furthermore, 30 clustering provides a basis for downstream analyses, such as differential expression, 31 trajectory inference, and cell-cell interaction. In single-cell transcriptomics, a variety of 32 clustering algorithms have been used and Leiden clustering has emerged as a performant 33 choice(1). Leiden clustering can be extended to consider multiomics data via the Leiden 34 multiplex functionality(2). 35 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint 36 Progress in spatially resolved omics methods has empowered researchers with the ability to 37 map gene expression in a spatial manner, transcending conventional cell clustering 38 approaches(3). With spatial omics, scientists can discern higher-order tissue structures, 39 termed spatial domains, by integrating spatial information alongside gene expression data. 40 The identification of spatial domains through spatial clustering has emerged as a standard 41 practice in constructing spatial atlases. This is instrumental in visualizing tissue anatomy, 42 delineating tissue spatial continuity, pinpointing domain-specific marker genes, and 43 unravelling domain-dependent molecular regulatory networks. Performance of spatial 44 domain identification improves when leveraging the spatial information compared to non-45 spatial methods(4). 46 47 Leiden clustering has been typically categorised as a “non-spatial” clustering method. 48 However, Leiden multiplex can consider spatial embeddings during clustering, thus 49 rendering prior assumptions of the non-spatial nature of the Leiden algorithm untrue. The 50 Leiden algorithm clusters nodes in a network by optimising a quality function, in a simple 51 case this can be the modularity, which maximises the differences between the actual 52 number of edges in a community and the expected number of such edges under a null 53 model. In single-cell transcriptomics, the cells (nodes) are connected to other cells based on 54 the distance between cells in the gene expression space (edges), usually in a dimensionality 55 reduced latent space. Leiden multiplex enables the user to define an arbitrary number of 56 networks (layers) with the same set of nodes that describe different modalities of edges 57 between the nodes. So, in spatially resolved omics data, the spatial neighbourhood can be 58 encoded by defining a spatial connectivity (based on e.g. Euclidean distance) as the weight 59 of edges between nodes (cells or spots). 60 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint

and Discussion 61 In this study we review how Leiden clustering can utilise spatial information through selection 62 of spatially variable genes (SVGs) instead of highly variable genes (HVGs)(5), spatially-63 aware dimensionality reduction through MUTLISPATI-PCA(6) (msPCA), and explicitly 64 modelling the spatial embedding in the Leiden multiplex clustering (SpatialLeiden) (Figure 65 1a). We demonstrate their application to a 10x Visium spatial transcriptomics dataset of the 66 mouse dorsolateral prefrontal cortex (DLPFC), the most widely used benchmark dataset for 67 spatial clustering methods(7). This dataset consists of spatial gene expression data and 68 histology images of 3 replicate slices from 4 donor mice, together with ground truth 69 annotation of anatomical domains in those tissue samples. We compare the performance of 70 a non-spatially aware and SpatialLeiden clustering to two widely used spatially-aware 71 domain detection tools, SpaGCN and BayesSpace(8,9), and evaluating performance of the 72 tools (Figure 1 b, c). 73 74 While use of SVGs over HVGs yielded only minor improvements, we observed substantial 75 improvement in performance when using spatially aware dimensionality reduction (msPCA) 76 and using SpatialLeiden over non-spatial Leiden, revealing a better representation of the 77 neocortex layering pattern (Figure 1b). We quantitatively evaluated performance of the 78 different clustering strategies using the Adjusted Rand Index (ARI) and Normalised Mutual 79 Information (NMI) score, showing significant improvements of SpatialLeiden over the non-80 spatial Leiden implementation, with performance that was better than SpaGCN and 81 comparable to BayesSpace (Figure 1 c, Supplemental Figure 1-3, Supplemental table 1-82 2) at a fraction of the processing time. SpatialLeiden performed favourably when we 83 compared its performance to other tools in a recent benchmark study, ranking 5th of 15 84 tools(4), (see Supplemental Methods for further details). 85 86 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint 87 Figure 1: SpatialLeiden workflow. (a) Schematic of data processing and modelling steps 88 to enable spatially-aware Leiden clustering. Feature selection is performed by spatially 89 variable genes (SVGs), dimensionality reduction is performed by a spatially aware 90 MULTISPATI-PCA (msPCA), and clustering is performed by the Leiden layer multiplex 91 algorithm with both gene expression and spatial embeddings (SpatialLeiden). (b) Histology 92 and manually annotated neocortex layered domains for the mouse brain DLPFC (slice 93 151673) and spatial domains detected by Leiden, SpatialLeiden, SpaGCN and BayesSpace. 94 (c) Boxplot of Adjusted Rand Index (ARI) and Normalised Mutual Information (NMI) for all 12 95 DLPFC samples. Center line: median; box limits: upper and lower quartiles; whiskers: 1.5× 96 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint interquartile range; dots: outliers; asterisks: significance (FWER, one-sided Wilcoxon signed-97 rank test). 98 99 As with other multi-modal clustering approaches, careful consideration has to be paid to a 100 number of parameters including the definition and weighting of neighbours across 101 modalities, the resolution to be applied to each modality, and the weight of each modality 102 (Supplemental Figure 4). Furthermore, the spatial relationship between cells can be 103 modelled in different ways; a regular grid pattern is suitable for Visium (isometric) and binned 104 Stereo-seq (square), while for imaging-based spatially resolved transcriptomics methods 105 Delaunay triangulation or k-nearest neighbours can be used to define the topological layer. 106 107 To investigate the performance of SpatialLeiden across technologies, tissues and topological 108 modelling we analysed a number of datasets (Stereo-Seq mouse embryo(10), BaristaSeq 109 mouse brain primary cortex, MERFISH mouse brain hypothalamus preoptic area(11), 110 osmFISH mouse brain somatosensory cortex, STARmap mouse brain medial prefrontal 111 cortex(12,13), STARmap1k mouse brain visual cortex) and demonstrated exceptional 112 improvements over non-spatially aware Leiden clustering. SpatialLeiden ranked as the best 113 performing spatial domain clustering methods for Stereo-Seq, MERFISH and osmFISH, 114 while demonstrating top tier performance for the other datasets (Figure 2). For imaging-115 based spatial transcriptomics methods we found that modelling the topological 116 neighbourhood using the 10 k-nearest neighbours generally performed better than using 117 Delaunay triangulation. All non-Stereo-seq dataset were processed within 2 minutes, utilising 118 less than 400 MB of RAM. All Stereo-seq samples were processed within 8 minutes utilising 119 a maximum of 3.5 GB of RAM (Supplemental Figure 5). 120 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint 121 Figure 2: Performance of SpatialLeiden across technologies and tissues: (a) Stereo-122 Seq of the mouse embryo at various development stages; (b) BaristaSeq of mouse primary 123 cortex; (c) MERFISH mouse brain hypothalamus preoptic area; (d) osmFISH of mouse 124 somatosensory cortex; (e) STARmap mouse brain medial prefrontal cortex; and (f) 125 STARmap* of mouse visual cortex. Performance metrics of other tools are taken from Yuan 126 et al 2024(4). SpatialLeiden was run with 5 different random seeds and median results were 127 reported per sample. 128

129 Our results show that the reference implementation of the Leiden algorithm can indeed be 130 used as a spatially-aware clustering algorithm. Subsequent studies that compare spatially-131 aware clustering algorithms should clearly state that they compare to non-spatial 132 implementation of Leiden, rather than misclassifying Leiden as a non-spatial algorithm. 133 We describe the different steps at which spatial awareness can be introduced into the 134 analysis, and our implementation allows easy parameterisation of key considerations for 135 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint modelling gene and spatial modalities. While many spatial domain clustering tools rely on 136 spatially-aware dimensionality reduction approaches this is often followed by non-spatial 137 clustering and we expect these methods to improve with spatially-aware clustering such as 138 SpatialLeiden. 139 140 The same way Leiden became the method of choice for clustering of single-cell data, we 141 believe that SpatialLeiden will become the method of choice for spatial data owing to its 142 efficiency, simplicity, and ease of integration into existing analysis workflows. 143

144 Data processing 145 Data was analysed using python (v3.10.14), Scanpy(14) (v1.10.1), and Squidpy(15) (v1.4.1). 146 Topological neighbourhood graph generation 147 The neighbors of each cell were defined depending on the technology. For datasets with a 148 regular grid the neighbors were defined using squidpy.gr.spatial_neighbors with coord_type 149 ‘grid’ and n_neighs set to 6 or 4 for Visium and Stereo-seq, respectively. For all other 150 datasets the neighbors were defined either using Delaunay triangulation 151 (squidpy.gr.spatial_neighbors with delaunay=True) or using the 10 nearest-neighbors 152 (squidpy.gr.spatial_neighbors with coord_type `generic` and n_neighs set to 10). The 153 untransformed neighbourhood graph ‘spatial_connectivities’ was used as is for regular grids 154 as all connections are equidistant. For Delaunay triangulation and kNN the 155 ‘spatial_distances’ were transformed to ‘spatial_connectivities’ via the following formula: 156 1 − 𝑑 𝑑𝑚𝑎𝑥 157 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint HVG and SVG detection 158 Highly variable genes (HVGs) were detected based on filtered count data 159 (scanpy.pp.highly_variable_genes with flavor ‘seurat_v3’). To detect spatially variable genes 160 (SVGs), first the neighbours were defined as described above and the Moran’s I score was 161 calculated for all genes using squidpy.gr.spatial_autocorr with mode ‘moran’ and selecting 162 the top 3,000 scoring genes. Gene selection was only performed for capture-based spatial 163 transcriptomics technologies (Visium, Stereo-seq). 164 MULTISPATI-PCA 165 We implemented MULTISPATI-PCA(6) in python (https://github.com/HiDiHlabs/multiSPAETI 166 v0.1.0) and used it to perform spatially-aware dimensionality reduction. The topological 167 neighbourhood graph ‘spatial_connectivities’ was used to calculate 30 components 168 (corresponding to the 30 largest eigenvalues) based on the 3,000 HVGs / SVGs (Visium, 169 Stereo-seq) or all genes in the case of imaging-based technologies (STARmap, STARmap*, 170 MERFISH, BaristaSeq, osmFISH) as for normal PCA. 171 Latent neighbourhood graph generation 172 The topological neighbourhood graph in the Visium grid ‘spatial_distances’ from 173 squidpy.gr.spatial_neighbors (as described in SVG detection) is used as the spatial layer for 174 Leiden clustering. To build the neighbourhood graph in latent space of gene expression we 175 first calculated the first 30 principal components based on the top 3,000 variable genes 176 (scanpy.tl.pca). We identified the 15 nearest-neighbours per spot (scanpy.pp.neighbors) 177 based on the PCA or MULTISPATI-PCA results from either the HVGs / SVGs (Visium, 178 Stereo-seq) or for all genes in the case of imaging-based technologies (STARmap, 179 STARmap*, MERFISH, BaristaSeq, osmFISH) and used the resulting ‘connectivities’. 180 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint Non-spatially aware Leiden 181 We used Leiden(1) (https://github.com/vtraag/leidenalg, v0.10.2) as implemented in Scanpy 182 with the default parameters and varied the resolution to achieve the correct number of 183 clusters for each of the DLPFC datasets following the approach of the SpaGCN.search_res 184 function (https://github.com/jianhuupenn/SpaGCN). 185 Spatially-aware Leiden multiplex (SpatialLeiden) 186 We implemented a spatially aware version of 187 Leiden(https://github.com/HiDiHlabs/SpatialLeiden) by using the Layer multiplex(16). An 188 additional graph encoding for the topological neighbourhood of the cells was added as 189 secondary layer in addition to the layer encoding gene expression in latent space. The 190 additional spatial layer was encoded as RBConfigurationVertexPartition as is the default for 191 the scanpy implementation for the latent space graph. The optimal clustering was identified 192 by running the Optimiser.optimise_partition_multiplex from leidenalg until convergence. As 193 only the ratio of the layer_weights is relevant the weight for the gene expression latent space 194 layer was kept at 1 and the weight for the topological neighbourhood was set depending on 195 technology and method of neighbourhood definition. The resolution for the latent space 196 partition was set by running the standard Leiden clustering and identifying the resolution 197 which yields the correct number of clusters. The resolution of the topological partition was 198 then varied to identify the correct number of clusters in the multiplex Leiden using the same 199 approach as described for the standard Leiden method. 200 Implementation and comparison to other spatial clustering 201 algorithms 202 Implementation and comparison to other spatial clustering algorithms is described in the 203 online methods. 204 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint Abbreviations 205 Spatially variable genes (SVGs). 206 Highly variable genes (HVGs). 207 MUTLISPATI-PCA (msPCA). 208 Dorsolateral prefrontal cortex (DLPFC). 209 Adjusted Rand Index (ARI). 210 Normalised Mutual Information (NMI). 211 Declarations 212 Ethics approval and consent to participate 213 Not applicable. 214 Consent for publication 215 Not applicable. 216 Availability of data and materials 217 This study used publicly available spatially resolved transcriptomics data of the mouse brain 218 DLPFC profiled on the 10x Visium platform (http://research.libd.org/spatialLIBD/). Public 219 Stereo-seq, STARmap, STARmap*, MERFISH, osmFISH, and BaristaSeq datasets were 220 downloaded from http://sdmbench.drai.cn/. Performance metrics of other tools presented in 221 Figure 2 were taken from Yuan et al(4). We also provide a repository with Jupyter Notebooks 222 for reproducing all results and figures of this study 223 https://github.com/HiDiHlabs/SpatialLeiden-Study. 224 225 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint Competing interests 226 The authors declare that they have no competing interests. 227 Funding 228 This research has received funding from the Federal Ministry of Education and Research of 229 Germany in the framework of SAGE (project number 031L0265). 230 Authors' contributions 231 NI conceived and designed the study. SS implemented the spatially-aware Leiden. NMB, SS 232 implemented the SpatialLeiden package. NMB implemented MULTISPATI-PCA in python. 233 SS performed code review. NMB performed data analysis. NI, NMB interpreted and 234 analysed results. NI, NMB, SS, RE proofread and corrected the manuscript. All authors 235 contributed to the article and approved the submitted version. 236

237 We thank organisers and participants of the de.NBI BioHackathon SpaceHack 2.0 project in 238 Bielefeld, Germany in December 2023 (Supplemental Table 3). 239

240 1. Traag VA, Waltman L, van Eck NJ. From Louvain to Leiden: guaranteeing well-241 connected communities. Sci Rep. 2019 Mar 26;9(1):5233. 242 2. Bredikhin D, Kats I, Stegle O. MUON: multimodal omics analysis framework. Genome 243 Biol. 2022 Feb 1;23(1):42. 244 3. Vandereyken K, Sifrim A, Thienpont B, Voet T. Methods and applications for single-cell 245 and spatial multi-omics. Nat Rev Genet. 2023 Aug;24(8):494–515. 246 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint 4. Yuan Z, Zhao F, Lin S, Zhao Y, Yao J, Cui Y, et al. Benchmarking spatial clustering 247

with spatially resolved transcriptomics data. Nat Methods. 2024 Mar 15;1–11. 248 5. Chen C, Kim HJ, Yang P. Evaluating spatially variable gene detection methods for 249 spatial transcriptomics data. Genome Biol. 2024 Jan 15;25(1):18. 250 6. Dray S, Saïd S, Débias F. Spatial Ordination of Vegetation Data Using a Generalization 251 of Wartenberg’s Multivariate Spatial Correlation. J Veg Sci. 2008;19(1):45–56. 252 7. Maynard KR, Collado-Torres L, Weber LM, Uytingco C, Barry BK, Williams SR, et al. 253 Transcriptome-scale spatial gene expression in the human dorsolateral prefrontal cortex. 254 Nat Neurosci. 2021 Mar;24(3):425–36. 255 8. Hu J, Li X, Coleman K, Schroeder A, Ma N, Irwin DJ, et al. SpaGCN: Integrating gene 256 expression, spatial location and histology to identify spatial domains and spatially 257 variable genes by graph convolutional network. Nat Methods. 2021 Nov;18(11):1342–51. 258 9. Zhao E, Stone MR, Ren X, Guenthoer J, Smythe KS, Pulliam T, et al. Spatial 259 transcriptomics at subspot resolution with BayesSpace. Nat Biotechnol. 2021 260 Nov;39(11):1375–84. 261 10. Chen A, Liao S, Cheng M, Ma K, Wu L, Lai Y, et al. Spatiotemporal transcriptomic atlas 262 of mouse organogenesis using DNA nanoball-patterned arrays. Cell. 2022 May 263 12;185(10):1777-1792.e21. 264 11. Moffitt JR, Bambah-Mukku D, Eichhorn SW, Vaughn E, Shekhar K, Perez JD, et al. 265 Molecular, spatial, and functional single-cell profiling of the hypothalamic preoptic region. 266 Science. 2018 Nov 16;362(6416):eaau5324. 267 12. Wang X, Allen WE, Wright MA, Sylwestrak EL, Samusik N, Vesuna S, et al. Three-268 dimensional intact-tissue sequencing of single-cell transcriptional states. Science. 2018 269 Jul 27;361(6400):eaat5691. 270 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint 13. Li Z, Zhou X. BASS: multi-scale and multi-sample analysis enables accurate cell type 271 clustering and spatial domain detection in spatial transcriptomic studies. Genome Biol. 272 2022 Aug 4;23(1):168. 273 14. Wolf FA, Angerer P, Theis FJ. SCANPY: large-scale single-cell gene expression data 274 analysis. Genome Biol. 2018 Feb 6;19(1):15. 275 15. Palla G, Spitzer H, Klein M, Fischer D, Schaar AC, Kuemmerle LB, et al. Squidpy: a 276 scalable framework for spatial omics analysis. Nat Methods. 2022 Feb;19(2):171–8. 277 16. Mucha PJ, Richardson T, Macon K, Porter MA, Onnela JP. Community Structure in 278 Time-Dependent, Multiscale, and Multiplex Networks. Science. 2010 May 279 14;328(5980):876–8. 280 281 .CC-BY 4.0 International licensemade available under a (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is The copyright holder for this preprintthis version posted August 23, 2024. ; https://doi.org/10.1101/2024.08.23.609349doi: bioRxiv preprint