Characterizing spatial gene expression heterogeneity in spatially resolved single-cell transcriptomic data with nonuniform cellular densities

Overview of MERINGUE

Given a set of spatial positions such as those corresponding to single cells, MERINGUE first represents these cells as neighborhoods using Voronoi tessellation (Fig. 1A). In Voronoi tessellation, planes are partitioned into neighborhoods where a neighborhood for a cell consists of all points closer to that cell than any other (Okabe et al. 1992). Cells are then considered adjacent if their neighborhoods share an edge. For biological interpretability, we further require adjacent cells to be within a certain spatial distance in space to accommodate realistic length scales of cellular interactions. This neighborhood representation of cells accommodates varying neighborhood sizes and distances between cells and thus can characterize cell types and tissues with nonuniform densities. We also find that such neighborhood adjacency relationships to be more stable than k-nearest-neighbor or k-mutual-nearest-neighbor relationships, because such relationships require k to be specified beforehand and a single k value may not be appropriate for all densities and regions within a spatially resolved data set (Supplemental Fig. S2A). MERINGUE encodes these adjacency relationships using a binary adjacency weight matrix W with a weight of 1 if two data sets are adjacent and 0 otherwise (Fig. 1A). Such adjacency relationships are not restricted to 2D and thus can accommodate 3D information, such as from imaging of multiple slices of the same tissue or 3D volumetric imaging of a tissue block, if available (Lee et al. 2015; Wang et al. 2018).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 1.

Overview of MERINGUE. (A) MERINGUE encodes spatial relationships among spatially resolved data sets, such as cells, using a binary adjacency weight matrix W. (Left) Two cells are considered adjacent if their neighborhoods inferred from Voronoi tessellation share an edge. (Right) The binary adjacency weight matrix W is visualized by plotting all cells in space with a red line connecting cells if cells are adjacent and no line otherwise. (B) MERINGUE identifies genes with spatially heterogeneous expression using W to compute and evaluate the significance of a spatial autocorrelation index I for each gene. When a gene's expression magnitude between spatially adjacent cells are highly correlated, I will be positive, indicative of spatial heterogeneity. Three simulated gene expression profiles are provided to illustrate examples of high and low spatial heterogeneity with red indicating high expression and blue indicating low expression. (C) MERINGUE groups identified spatially heterogeneous genes into primary spatial patterns by computing a spatial cross-correlation index for every gene pair. (Top) The resulting spatial cross-correlation matrix is used to construct a hierarchical dendrogram. (Bottom) Dynamic tree cutting is applied to partition genes into patterns. Groups of genes are Z-scored and averaged, with empty regions filled in using Akima interpolation to visualize final patterns. (D) MERINGUE identifies gene expression patterns that may be indicative of putative cell–cell communication using an inter-cell-type weight matrix W_ict between two cell types, which can then be used to compute an inter-cell-type spatial cross-correlation index iSCI between two genes (Top left). Two cell types, A and B, are shown as green triangles and orange squares, respectively. W_ict is visualized with a red line for cells of cell type A spatially adjacent to cells of cell type B. (Top right) Cell type–A cells express gene A at variable levels, whereas cell type–B cells express gene B at variable levels, with red indicating high expression and blue indicating low expression. (Bottom left) Cell type–A cells do not express gene B and cell type–B cells do not express gene A, resulting in a generally negative Pearson's correlation (R) between the two genes. (Bottom right) However, expression of gene A in cells of cell type A is highly correlated with the expression of gene B in spatially adjacent cells of cell type B, resulting in a positive iSCI. The significance of this iSCI is assessed by permutation.

Next, to identify genes with spatially heterogeneous expression, given a matrix of normalized gene expression magnitudes for the same set of spatially resolved cells, MERINGUE uses this adjacency weight matrix W in calculating Moran's I, a global measure of spatial autocorrelation popular in geospatial analysis, for each gene's expression magnitude (x) across the population of N cells (Moran 1950): When a gene's expression magnitude (x) between spatially adjacent cells (W_ij = 1) are positively correlated, Moran's I will be positive (Fig. 1B), indicative of spatial gene expression heterogeneity. Moran's I has a closed form, allowing P-values to be derived without reliance on computationally intensive permutations (Supplemental Fig. S2B; Moran 1950).

To further characterize the scale of significant spatial gene expression heterogeneity, using the same matrix of normalized gene expression magnitudes and adjacency weight matrix W, MERINGUE calculates a local indicator of spatial association (LISA) for each gene (Anselin 1995): When a gene's expression values (x) in a given cell (i) is positively correlated with that cell's spatially adjacent neighbors, the cell's LISA for the given gene will be highly positive. Again, LISA has a closed form, allowing P-values to be derived quickly. As such, MERINGUE defines the percent of cells with statistically significant LISAs as the percent of cells driving a spatially heterogeneous gene expression pattern. This use of LISA guards against the identification of spatially heterogeneous genes driven by small hotspots or outliers. Simulations suggest that false positives may be effectively eliminated by restricting to spatial heterogeneity driven by >5% of cells (Supplemental Fig. S2C).

Finally, to summarize genes into primary spatial patterns, MERINGUE calculates a spatial cross-correlation index between all pairs of genes identified with significant spatially heterogeneous expression driven by a sufficient percentage of cells: When one gene's expression magnitude (x) in a given cell (i) is positively correlated with another gene's expression magnitude (y) in the cell's spatially adjacent neighbors (j), the SCI for this gene pair will be positive. MERINGUE computes this spatial cross-correlation index for all gene pairs to derive a spatial cross-correlation matrix that is then used for hierarchical clustering and dynamic tree cutting to group these genes into primary spatial patterns (Fig. 1C; Langfelder et al. 2008).

In addition, MERINGUE further builds on this spatial cross-correlation index to identify spatially cross-correlated gene expression patterns that may be indicative of cell–cell communication. In particular, communicating cell types may express higher levels of particular ligand genes while being spatially adjacent to cells that express higher levels of corresponding receptor genes or vice versa. Thus, to identify such gene expression patterns that may be indicative of putative cell–cell communication, MERINGUE constructs an adjacency weight matrix W to only include adjacency relationships between cell types and calculates the spatial cross-correlation statistics for known receptor and ligand genes (Ramilowski et al. 2015). In this manner, when a receptor gene's expression magnitude (x) in a given cell (i) of cell type A is positively correlated with the corresponding ligand gene's expression magnitude (y) in cells of cell type B among the cell's spatially adjacent neighbors (j), the inter-cell-type SCI for this cell type pair will be highly positive. Statistical significance can then be assessed by permutation testing (Fig. 1D).

MERINGUE identifies genes with spatially heterogeneous expression patterns and is robust to changes in cellular densities

As a proof of principle, we first applied MERINGUE to spatial transcriptomics (ST) data of the mouse main olfactory bulb (MOB) and Slide-seq data of the mouse cerebellum (Ståhl et al. 2016; Rodriques et al. 2019). Briefly, for ST and Slide-seq, RNAs from tissue sections are captured onto an array of DNA barcoded spots or a monolayer of DNA barcoded beads, respectively. By resolving the DNA barcodes, both approaches enable matching of detected RNA abundances with their original spatially resolved spots or beads, resulting in RNA sequencing measurements with uniformly gridded 2D positional information. To validate MERINGUE, we expected that identified spatially heterogeneous genes in the MOB should mark transcriptionally distinct and spatially organized cell layers or combinations of cell layers (Fig. 2A; Supplemental Fig. S3A). Indeed, when we applied MERINGUE to analyze 7365 genes among 260 spots, of the 834 identified as significantly spatially heterogeneous genes (adjusted P-value < 0.05) driven by >5% of spots (Fig. 2B; Supplemental Fig. S3B; Supplemental Table S1), 90% (754/834) overlapped with genes that are significantly differentially expressed genes across cell layers (adjusted P-value < 0.05) as identified from ANOVA testing. Furthermore, these 834 spatially heterogeneous genes can be further partitioned into five primary spatial patterns that correspond to cell layers and combinations of cell layers as expected (Fig. 2C; Supplemental Fig. S3C). One well-characterized aspect of spatial organization in the MOB involves the convergence of axonal projections from olfactory receptor neurons expressing a given olfactory receptor (Olfr) into glomerular neuropils at fixed locations in the glomerular cell layer of the olfactory bulb (Ressler et al. 1994; Vassar et al. 1994; Mombaerts et al. 1996). Therefore, as an additional validation, we evaluated whether Olfr genes were spatially heterogeneous in a pattern that corresponds to the glomerular and surrounding cell layers. Although individual Olfr genes are very lowly expressed such that detection was generally limited to only a few copies in a few spots (Supplemental Fig. S3D), rendering assessment of spatial heterogeneity for individual Olfr genes infeasible, by aggregating the expression of all detected Olfr genes, we validate that MERINGUE was able to identify significant spatial heterogeneity (P-value = 0.0000283). The spatial expression pattern further corresponded approximately to the glomerular and surrounding cell layer as expected (Supplemental Fig. S3E). For Slide-seq data of the mouse cerebellum, we applied MERINGUE to analyze 9762 genes among 1589 beads previously annotated to correspond to the Purkinje layer (Supplemental Fig. S4A). We validate that Aldoc (also known as zebrin II) is identified as among the most significantly spatially heterogeneous genes (adjusted P-value < 0.05, >5% beads) (Supplemental Table S2), consistent with observations from the original publications (Rodriques et al. 2019).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 2.

Application of MERINGUE to 2D spatial transcriptomic data of the main olfactory bulb (MOB) and 3D-aligned ISH data of the Drosophila melanogaster embryo. (A) Spatially unaware single-cell clustering analysis identifies five transcriptionally distinct clusters corresponding to various known cell layers in the MOB. Spatial spots are colored based on their inferred cell layer annotation. (B) MERINGUE identifies genes with significantly spatially heterogeneous expression in the MOB. Select genes are shown. (C) MERINGUE groups genes with significantly spatially heterogenous expression in the MOB into five primary spatial patterns. Select patterns are shown. (D) MERINGUE's adjacency weight matrix visualized for aligned 3D in situ hybridization data of the D. melanogaster embryo. Each point is an aligned cell. Cells are connected with a red line if they are inferred to be adjacent. A top view and rotated side view are shown. (E) MERINGUE groups genes into spatial patterns in the D. melanogaster embryo. Representative genes from select identified patterns are shown.

We next compared MERINGUE to previously published computational methods for analyzing spatially resolved transcriptomic data, SpatialDE and SPARK (Svensson et al. 2018; Sun et al. 2020). We applied each method to analyze 7365 genes among 260 spots in the MOB to identify spatially heterogeneous genes (Supplemental Methods). We found the resulting significance of spatial heterogeneity in terms of −log₁₀(adjusted P-value) to be highly correlated across genes between all tested computational methods (R = 0.914 between MERINGUE and SpatialDE, R = 0.898 between MERINGUE and SPARK) (Supplemental Fig. S5A,B). The resulting set of significantly spatially heterogeneous genes identified by each tested computational method using a common significance threshold (adjusted P-value < 0.05) were also highly overlapping (Supplemental Fig. S5C). We further evaluated the computational efficiency of each method in terms of runtime and memory usage as a function of the number of genes and the number of cells in the data set (Supplemental Methods). We found that MERINGUE achieves improved computational efficiency compared to previously published computational methods (Supplemental Fig. S5D,E). Thus, MERINGUE is capable of identifying spatially heterogeneous genes consistent with previously published approaches in a scalable manner.

We developed MERINGUE to accommodate the nonuniform cellular densities inherent to tissues. Thus, we reasoned that changes in cellular densities should not substantially impact MERINGUE's ability to identify spatially heterogeneous genes. To assess MERINGUE's robustness to spatial variations in cellular densities, we artificially induced nonuniformity in the spatial distribution of ST spots by distorting their positional coordinates (Supplemental Methods; Supplemental Fig. S6A). Owing to its use of a distance-agnostic binary weight matrix, MERINGUE's resulting significance of spatial heterogeneity across genes was highly correlated between the uniform and nonuniform case as expected (Spearman's ρ = 0.862) (Supplemental Fig. S6B). Likewise, although MERINGUE was able to identify 834 significantly spatially heterogeneous genes (adjusted P-value < 0.05, >5% of spots) in the uniform density case, 544 (65%) of these genes were recovered in the nonuniform density case with the same adjusted P-value and spot percentage thresholds. The discrepancies between the uniform and nonuniform cases can be largely attributed to changes in the binary weight matrix (Supplemental Fig. S6C). Because SpatialDE and SPARK incorporate Euclidean distances between cells in their evaluation of spatial patterns, we reasoned that spatial variations in cellular density would impact their ability to identify spatially heterogeneous genes. We thus applied the same uniform and artificially induced nonuniform case comparison. As expected, the resulting significance of spatial heterogeneity across genes was less well correlated between the uniform and nonuniform density case for both SpatialDE (Spearman's ρ = 0.427) and SPARK (Spearman's ρ = 0.418) (Supplemental Fig. S6D). Likewise, although SpatialDE was able to identify 360 significantly (adjusted P-value < 0.05) spatially heterogeneous genes in the uniform density case, only 56 (16%) of these genes were recovered in the nonuniform case with the same adjusted P-value threshold. Similarly, SPARK was able to identify 664 significantly (adjusted combined P-value < 0.05) spatially heterogeneous genes in the uniform case, but only 66 (10%) of these genes were recovered in the nonuniform case with the same adjusted P-value threshold.

MERINGUE integrates 3D and multilayer tissue information

Although spatially resolved transcriptomics measurements generally provide positional information in the imaging (x–y) plane, z-direction information can be obtained through optically scanning through imaging planes or sequential tissue sections. To show integration of z-direction information, we first applied MERINGUE to 3D in situ hybridization (ISH) data, aligned across multiple stage 6 Drosophila melanogaster embryos for 84 selected marker genes (Fig. 2D; Supplemental Methods; Fowlkes et al. 2008; Karaiskos et al. 2017). The role of spatial patterning in shaping cellular identities has been well established in the D. melanogaster embryo and as such, the 84 marker genes were previously chosen for their known spatial patterning. Indeed, we validated that all 84 genes are identified by MERINGUE as significantly spatially heterogeneous (adjusted P-value < 0.05, >5% spots) as expected. We further validated that these genes could be grouped by MERINGUE into 14 primary spatial patterns that correspond to known regionally confined developmental fates and layers of the segmentation gene network (Fig. 2E; Supplemental Fig. S7; Supplemental Table S3; Ingham 1988; Karaiskos et al. 2017). For example, pattern 2 corresponds to the mesoderm and includes mesoderm determinant gene twist (twi), whereas pattern 13 corresponds to the thoracic segments and includes known gap gene Kruppel (Kr) (Preiss et al. 1985; Leptin 1991). Similarly, patterns 9 and 12 correspond to two spatially alternating striped patterns that include known pair-rule genes even skipped (eve) and odd skipped (odd), respectively (Macdonald et al. 1986; Coulter et al. 1990).

Alternatively, z-information may be derived through serial sections. We thus next applied MERINGUE to spatial transcriptomic data of four consecutive histological sections of a human breast cancer biopsy (Supplemental Methods; Ståhl et al. 2016). Analyzing each section independently, we identified 414 genes that show significant spatial variability (adjusted P-value < 0.05, >5% spots) in at least one section out of 6214 genes tested (Supplemental Table S4). Because the distance between cells across serial sections are greater than the distances between cells within sections, we sought to identify spatial patterns consistent across layers by testing for spatial correlation between mutual nearest neighbors in space across sections (Supplemental Fig. S8A). Such a multilayer integrated approach confirmed 242 significantly spatially heterogeneous genes as being consistent across sections (Supplemental Table S4; Supplemental Fig. S8B). Of the remaining 172 genes that were identified as significantly spatially heterogeneous in individual sections but not across sections, visual inspection showed that although these genes show spatial variability within sections, there was minimal correspondence across sections (Supplemental Fig. S8C). Such transcriptional patterns may be indicative of layer-specific subpopulations or transcriptional features. For structurally stereotypic tissues, consistency across tissue sections may be used as an additional criterion for identifying functionally relevant spatial patterns. Thus, MERINGUE is capable of accommodating 3D information to identify spatially heterogeneous genes in 3D as well as genes with spatial expression patterns consistent across serial sections.

MERINGUE identifies spatial patterns in the mouse hypothalamic preoptic region using spatially resolved single-cell gene expression data by MERFISH

Particularly in complex organs such as the mammalian brain, the ability to identify and interrogate the spatial organization of cell types may provide additional insights into potential functional roles underlying the spatial organization of neuronal populations (y Cajal 1911; Amaral and Witter 1989; Arber 2012). We applied MERINGUE to analyze spatially resolved single-cell transcriptomic data of the hypothalamic preoptic region obtained using multiplexed error-robust fluorescence in situ hybridization (MERFISH) (Moffitt et al. 2018). Briefly, MERFISH allows individual RNA molecules in cells to be imaged and identified by using a combinatorial labeling strategy that encodes RNA species with error-robust barcodes that can be read out bit-by-bit using sequential rounds of single-molecule fluorescence in situ hybridization (Chen et al. 2015). MERFISH has enabled simultaneous detection and identification of thousands of targeted RNA species, which can then be segmented into cells to provide spatially resolved single-cell transcriptome measurements (Chen et al. 2015; Xia et al. 2019). Moffitt et al. (2018) previously used a 155 gene panel to characterize the hypothalamic preoptic region (1.8 × 1.8 × 0.6 mm, Bregma +0.26 to −0.34) in adult mice to identify 31 excitatory and 39 inhibitory neuronal subtypes in addition to non-neuronal cell types using graph-based community-detection clustering analysis that relies solely on the gene expression of profiles of cells without considering the spatial information.

We applied MERINGUE to analyze the 155 genes along with five blank control barcodes, DAPI, and poly(dT) signals as negative controls within each cell type and subtype to identify additional aspects of spatial heterogeneity. Applying a rigorous approach to identify genes with spatially heterogeneous expression patterns that are consistent across tissue layers and reproducible across animals (Supplemental Methods), we were able to identify at least one such spatially heterogeneous gene in 34 of 83 cell types and subtypes analyzed (Fig. 3A; Supplemental Fig. S9; Supplemental Table S5). None of the blank control barcodes, DAPI, or poly(dT) signals were identified as consistently spatially variable. MERINGUE further identified significant spatial gene expression heterogeneity within neuronal subtypes in both the anterior and posterior of the preoptic region. Likewise, spatial gene expression heterogeneity was identified in both inhibitory and excitatory neuronal subtypes. These aspects of spatial heterogeneity were consistent with previous published spatially unaware variance and principal components-based analyses and visual inspection (Moffitt et al. 2018).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 3.

MERINGUE identifies spatial heterogeneity within cell types in the preoptic region of the mouse hypothalamus using MERFISH. (A, left) Expression of three sample genes in Inhibitory I-6 neurons in female naïve animal 7 (FN7). Each point is a cell. Cells are colored by expression with red denoting high expression and blue denoting low expression. Cells that are not I-6 cells are colored in gray. Gad1 is highly expressed in all I-6 neurons, whereas Necab1 and Sema3cshow significant spatial variation. (Right) Expression of Sema3c in I-6 in adjacent tissue sections in FN7 show similar spatial patterning. (B, left) Expression of three sample genes in Inhibitory I-11 neurons in FN7. Again, Gad1 is highly expressed in all I-11 cells, but Gabra1 and Nos1 show significant spatial variation. (Right) Expression of Nos1 in I-11 neurons in adjacent tissue sections in FN7 show similar spatial patterning. (C) Expression of Nos1 in I-11 neurons in representative male naïve animals MN5, MN8, and MN9 show similar spatial patterns to the female animal in B. (D) Fraction of Nos1⁺I-11 cells for male and female animals across tissue layers from the anterior to posterior preoptic region. Each dot represents one tissue layer in one animal. Lines represent fitted curves for males and female animals. (E) Expression of Tacr1 in E-8 neurons in female (top) and male (bottom) animals. (F) Scale spatial heterogeneity of Tacr1 in E-8 neurons for male and female animals across tissue layers from the anterior to posterior preoptic region. Each dot represents one tissue layer in one animal. Lines represent fitted curves for male and female animals.

By providing a quantitative framework to systematically identify and evaluate the statistical significance of spatial gene expression heterogeneity, MERINGUE identified that cells of inhibitory subtype I-6 in the anterior of the preoptic region can be partitioned into a superior and inferior spatial lobe marked by higher and lower expression of Sema3c and Necab1, respectively (Fig. 3A). These patterns are consistent across adjacent tissue sections. Likewise, cells of inhibitory subtype I-11 in the posterior preoptic region can be partitioned into a medial and lateral spatial group marked by lower expression of Gabra1, higher expression of Nos1, and higher expression of Gabra1 lower expression of Nos1, respectively, and this partition is consistent across adjacent tissue sections (Fig. 3B). Alternatively, Gad1, which marks inhibitory cells, is highly expressed among all cells and does not show significant spatial heterogeneity as expected. Inhibitory subtype I-11 was previously identified to be specifically activated by male mating based on the expression of immediate early gene Fos (Moffitt et al. 2018). Although inhibitory subtype I-11 showed significant spatial heterogeneity in both male and female animals (Fig. 3C), we found the fraction of Nos1⁺ I-11 neurons to be significantly higher in males than females (Student's t-test P-value = 0.03656) (Fig. 3D). Esr1 and Irs4 were also identified to be significantly spatially heterogeneous in I-11 neurons in a pattern similar to Nos1 (Supplemental Fig. S9). Esr1 and Irs4 have been previously shown to display sex-differences in their expression (Xu et al. 2012; Moffitt et al. 2018). These observations suggest the potential presence of a finer Nos1⁺ I-11 neuronal subpopulation that is sexually dimorphic. Furthermore, MERINGUE generally identified concordant spatial gene expression heterogeneity in both male and female animals (Supplemental Fig. S10), but tachykinin receptor 1 (Tacr1, also known as neurokinin receptor 1) was identified as significantly spatially heterogeneous in excitatory subtype E-8 neurons only in male and not female mice (Fig. 3E). No other tested neuronal subtype was identified to show such consistently statistically significant sexually dimorphic spatial heterogeneity. Previously, E-8 neurons were identified to be activated in male mice during mating based on expression of Fos (Moffitt et al. 2018). However, E-8 neurons did not show a significant difference in terms of their proportion to all cells between female and male mice (Student's t-test P-value = 0.268). Likewise, we confirmed that the fraction of cells expressing Tacr1 in E-8 neurons is not significantly different between male and female mice (Student's t-test P-value = 0.429). However, when we quantified the fraction of cells driving the spatial heterogeneity of Tacr1 expression based on LISA, we observe a significant difference between male and female mice (Fig. 3F) (Student's t-test P-value = 0.01316). Tacr1 knockout mice have been previously observed to show deficits in sexual behavior (Berger et al. 2012). The sexually dimorphic spatial organization of Tacr1 expression in E-8 neurons may thus suggest a sexually dimorphic difference in connectivity responsible for its sexually dimorphic activation in sexual behavior. In this manner, MERINGUE enables quantitative and systematic evaluation of spatial gene expression heterogeneity within transcriptionally distinct cell subtypes from single-cell spatially resolved transcriptomics data.

Spatially informed clustering identifies transcriptionally and spatially distinct subtypes of cells

Spatial organization may play an important role in shaping cellular identities. Likewise, we may expect unsupervised clustering based on transcriptional profiles alone to recover spatially organized cell populations. However, for the aligned ISH data of the D. melanogaster embryo, we find such gene expression clustering analysis to aggregate cells expressing different pair-rule genes into a single transcriptional cluster owing to these cells sharing many other commonly up-regulated and down-regulated genes despite their spatially distinct organization (Fig. 4A,B), consistent with previously published analyses (Karaiskos et al. 2017). However, as our spatial analysis was able to distinguish between the two alternating striped spatial patterns marked by expression of pair-rule genes eve and odd, respectively (Fig. 2E), we sought to incorporate spatial information to help distinguish these spatially distinct but transcriptionally similar groups of cells.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 4.

Spatially informed clustering distinguishes spatially distinct subpopulations of cells. (A) Expression-based clustering of 3035 stage 6 D. melanogaster embryo cells with 84 marker genes by aligned ISH identifies approximately five transcriptionally distinct clusters. (Top) UMAP embedding colored by identified cluster annotations. (Bottom) Spatial coordinates colored by identified cluster annotations. (B) Expression of select marker genes on the UMAP embedding with red denoting high expression and blue denoting low expression. (C) Spatially informed clustering splits expression-based clusters in a spatially coherent manner. (Top) Again, UMAP embedding colored by identified spatially informed cluster annotations. (Bottom) Spatial coordinates colored by identified spatially informed cluster annotations. (D) Correspondence between expression-based clusters in A and spatially informed clusters in C highlights high correspondence between most clusters with the exception of one cluster being split into two. (E) UMAP embedding of populous inhibitory neuronal subtypes in one posterior preoptic tissue section from one animal measured using MERFISH, where each point is a cell colored by the original subtype annotations. (F) Same UMAP embedding as E in which each point is a cell colored by the spatially informed clustering annotation. Black dashed lines highlight clusters that have now split. (G) Correspondence between expression-based clusters in E and spatially informed clusters in F highlights high correspondence between most clusters with the exception of cells originally annotated as I-2 and I-11 now being split into two. (H) Same UMAP embedding as E in which each point is a cell colored by Nos1expression for cells originally annotated as I-11. (I) Spatial location of cells within the tissue colored by their spatially informed cluster assignment for cells originally annotated as I-11. (J) Same UMAP embedding as E, in which each point is a cell colored by Cplx3 expression for cells originally annotated as I-2. (K) Spatial location of cells within the tissue colored by their spatially informed cluster assignment for cells originally annotated as I-2. Regions corresponding to the BNST and SDN are highlighted with blue and red dashed lines, respectively. Representative slice in representative animal shown. (L) Percentage of activated cells based on Fos expression during female (FS) and male (MS) sexual behavior for spatially informed clusters C3 and C8 originally annotated as I-2. Boxes in the box plot denote the median values and inner quartile ranges (IQR), and whiskers denote 1.5 × IQR with additional outliers represented as points.

Briefly, as in expression-based clustering, we constructed a neighbor graph in which nodes are cells and nodes are connected with an edge if the represented cells that are within the k-most transcriptionally similar cells for some user-selected resolution parameter k. We incorporated spatial information by weighing the edges of the network by the distance (d) between two neighborhoods (i, j) in the adjacency representation W (Methods): (1/d_ij + 1) + 1. Again, use of such a neighborhood representation can accommodate the nonhomogenous density of cells in tissues compared to a Euclidean distance-based measure of spatial distance. In this manner, if two cells are closer in space (d_ij is small), their transcriptional similarity will give greater weight in the graph-based clustering. Incorporating these spatial weights into our graph-based clustering with all other parameters kept constant, we were able to split the cluster of cells expressing either eve or odd into two subpopulations, as desired (Fig. 4C). Moreover, such spatially informed clustering generally preserved all other subpopulations and did not result in additional splitting for other subpopulations (Fig. 4D). Furthermore, we showed using simulated data how such incorporation of spatial information can be used to distinguish transcriptionally identical but spatially distinct clusters of cells (Supplemental Fig. S11). In a biological setting, however, cells from the same cell type may populate spatially distinct locations, but such distinct spatial locations alone would not necessarily indicate the presence of finer subtypes. We thus suggest that such spatially informed clustering to be complementary to differential expression analysis, whereby identified spatially distinct cell subpopulations should be analyzed for significantly differentially expressed genes to ensure the presence of significant, likely subtle, transcriptional differences. Therefore, by incorporating spatial information, in conjunction with differential expression analysis, we can identify finer, transcriptionally and spatially distinct subpopulations.

Having shown that incorporation of spatial information via graph weighting can be applied to identify finer transcriptionally and spatially distinct subpopulations of cells, we next sought to apply this approach to identify finer neuronal subtypes in the preoptic region profiled by MERFISH (Moffitt et al. 2018). Focusing on inhibitory neurons, we performed spatially informed clustering analysis on all inhibitory cells in the same animal and tissue layer and compared resulting clusters to previous annotations (Fig. 4E,F; Supplemental Methods). We found that among the eight most populous inhibitory neuronal subtypes (clusters with >100 cells each), our spatially informed clustering was able to produce comparable clusters with the exception of I-2 and I-11, which were each split into two subtypes (Fig. 4G). I-11 was split into two subtypes, cluster 10 (C10) and cluster 5 (C5), that significantly differentially expressed genes including Nos1 (Fig. 4H; Supplemental Fig. S12A), consistent with our observations of significant spatial heterogeneity in Nos1 expression among I-11 neurons. Indeed, the two I-11 subtypes appear to be spatially distinct with C10 positioned more medially and C5 more laterally in the posterior preoptic region (Fig. 4I). Likewise, I-2 was split into two subtypes, cluster 3 (C3) and cluster 8 (C8), that significantly differentially up-regulated genes including Cplx3 and Dgkk, respectively (Fig. 4J; Supplemental Fig. S12B,C). Previously, I-2 neurons were observed to overlap with both the sexually dimorphic nucleus of the preoptic area (SDN-POA) as well as other anatomical nuclei such as the bed nucleus of the stria terminalis (BNST) (Moffitt et al. 2018). By refining I-2 into two finer subtypes, C8 is observed to overlap more so with the BNST, while C3 comparably more so with the SDN-POA (Fig. 4K; Supplemental Fig. S12C). I-2 neurons were previously observed to show sexually dimorphic activation during mating and aggression based on the expression of immediate early gene Fos (Moffitt et al. 2018). When we compare activation of the two I-2 subtypes based on significant Fos expression, we observe comparatively greater activation during mating behavior in one subtype than the other (Fig. 4L). This suggests that the activation in I-2 neurons observed previously may be driven by one of the two I-2 subtypes. Although tuning parameters for regular graph-based clustering without spatial information can also achieve splitting of I-2 and I-11, other inhibitory neuronal clusters can become over split (Supplemental Fig. S12D). Therefore, by incorporating spatial information via graph weighting, MERINGUE provides an alternative approach to tease apart spatially distinct subpopulations without impacting other transcriptionally distinct subtypes.

MERINGUE identifies putative cell–cell communication between cell types

Spatially resolved transcriptomic data offers opportunity to identify gene expression patterns that may be indicative of putative cell–cell communication between spatially colocalized cell types. Previous computational approaches for inferring cell–cell communication from single-cell RNA-sequencing data have relied on correlations or coexpression of receptor genes in one cell type and corresponding expression of ligand genes in another cell type (Ramilowski et al. 2015; Vento-Tormo et al. 2018; Smillie et al. 2019; Fan et al. 2020). Spatially resolved transcriptomic data provide the opportunity to infer potential cell–cell communication by identifying spatially complementary expression patterns between genes corresponding to interacting surface proteins such as receptors and ligands on spatially neighboring cells. To enable such analyses, we further build on MERINGUE's spatial cross-correlation functionalities by developing an inter-cell-type spatial cross-correlation function to identify potential complementary spatial patterns of gene expression across spatially colocalized cell types (Fig. 1D; Supplemental Fig. S13A–D). However, unlike the spatial autocorrelation function, this inter-cell-type spatial cross-correlation function is not solvable; thus significance must be assessed using permutation to derive a null model. We enhance computational efficiency by implementing a parallelized, adaptive permutation testing approach and assess significance using a permutation-based random label null model. We confirm using simulations that such a permutation-based assessment produces the expected type-I error rate (Supplemental Fig. S13E).

We first applied our approach to identify gene expression patterns that may be indicative of putative cell–cell communication between cells on beads corresponding to the Purkinje layer with cells on spatially adjacent beads in Slide-seq data of the mouse cerebellum (Fig. 5A). We used a set of more than 2500 known receptor-ligand pairs previously supported by orthogonal biological validations (Ramilowski et al. 2015). Restricting to well-detected (CPM > 0 in more than 30 cells) receptor genes in the Purkinje layer beads and well-detected ligand genes in the spatially adjacent beads, we applied MERINGUE to test for significant spatial cross-correlation between all receptor and ligand gene pairs. We identified statistically significant inter-cell-type spatial cross-correlation between expression of protein tyrosine phosphatase, receptor type Z, polypeptide 1 (Ptprz1 [also known as PTPζ]) in beads corresponding to Purkinje layer and expression of its ligand Ptn (secreted growth factor pleiotrophin) in spatially adjacent beads (Fig. 5B,C). Ptprz1 has been previously identified to be expressed by Purkinje neurons, whereas Ptn has been previously identified to distribute along Bergmann glial fibers in postnatally developing cerebellum (Matsumoto et al. 1994). Although this Slide-seq data set does not provide single-cell resolution, we confirm significant coexpression of Ptprz1 with Purkinje cell specific promoter Pcp2 (Fisher's exact test P-value = 2.3 × 10⁻¹⁸), suggestive that the Ptprz1 expression may be attributed to Purkinje cells within the Purkinje layer beads. Likewise, we confirm significant coexpression of Ptn with Slc1a3 (glutamate aspartate transporter [also known as GLAST]), a glutamate transporter expressed by Bergmann glia (Fisher's exact test P-value = 4.5 × 10⁻³⁶). In contrast, restricting to well-detected ligand genes in the Purkinje layer beads and well-detected receptor genes in the spatially adjacent beads (Fig. 5D,E), we do not identify significant spatial cross-correlation between any receptor and ligand gene pairs, including between Ptn expression in beads corresponding to Purkinje layer and Ptprz1 expression in spatially adjacent beads, indicative of the cell type specificity of inferred receptor-ligand interactions. Previous studies with cerebellar slice culture systems have shown that PTN-PTPRZ1 signaling is involved in the morphogenesis of Purkinje dendrites (Tanaka et al. 2003). The identification of such putative cell–cell communication between Bergmann glia and Purkinje cells may be suggestive of the potential for glial signals to actively regulate neuronal function and contribute to sustained plasticity in adult brains (Barres 2008).

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 5.

MERINGUE identifies putative cell–cell communication in the cerebellum using Slide-seq data. (A) Adjacency relationship between Slide-seq beads. Orange dots correspond to beads previously annotated as corresponding to the Purkinje layer. Green dots correspond to beads that are spatially adjacent. Gray lines connect each bead with its spatial neighbors and is agnostic to bead density. (B, left) Expression of receptor Ptprz1 in beads annotated to correspond to the Purkinje layer. (Right) Expression of corresponding ligand Ptn in spatially adjacent beads. Same select region highlighted. (C) Bar plot of −log₁₀(adjusted P-value) for the inter-cell-type spatial cross-correlation statistic of all receptors in Purkinje layer beads versus ligands in the spatially adjacent beads. Red line indicates alpha = 0.2 multiple testing corrected significance threshold. (D, left) Expression of ligand Psap in beads annotated to correspond to the Purkinje layer. (Right) Expression of corresponding receptor Gpr37l1 in spatially adjacent beads. Same select region highlighted as D. (E) Bar plot of −log₁₀(adjusted P-value) for the inter-cell-type spatial cross-correlation statistic of all ligands in Purkinje layer beads versus receptors in the spatially adjacent beads.

We next sought to identify gene expression patterns that may be indicative of putative cell–cell communication between cell types using single-cell resolution MERFISH data of the preoptic region. Previously, Moffitt et al. (2018) visually noted that Cyp19a1 (also known as aromatase) enriched inhibitory I-2 neurons displayed substantial spatial overlap with estrogen receptor (Esr1) enriched neuronal subtypes. CYP19A1 is an enzyme that converts testosterone to estrogen, thereby modulating steroid hormone signaling in the preoptic region. The spatial organization of these CYP19A1-enriched neuronal subtypes with ESR1-enriched cells suggest that estrogen synthesized by these CYP19A1-expressing neurons may be interacting with estrogen receptors on spatially adjacent ESR1-expressing cells in a paracrine manner. To quantitatively assess to putative cell–cell communication between neuronal subtypes and spatially neighboring cells via such paracrine signaling, we applied MERINGUE to assess for significant spatial cross-correlation between Cyp19a1 expression in all neuronal subtypes and Esr1 expression in spatially adjacent neurons. Indeed, we are able to identify statistically significant spatial cross-correlation between Cyp19a1 expression in I-2 neurons and Esr1 expression in spatially adjacent neurons in a manner that is consistent across tissue layers and reproducible across animals (Fig. 6A). In contrast, neuronal subpopulations I-13 also express Cyp19a1 but do not show significant spatial cross-correlation with Esr1 in surrounding cells (Fig. 6B,C; Supplemental Fig. S14A). Furthermore, we also apply MERINGUE to test for spatial cross-correlation between Cyp19a1 expression in all neuronal subtypes and androgen receptor (Ar) expression in spatially adjacent neurons and do not identify consistently significant associations (Supplemental Fig. S14B), thereby highlighting the nonrandomness of the Esr1 juxtaposition. This thus highlights MERINGUE's potential to quantitatively and systematically identify complementary gene expression patterns that may be indicative of cell–cell communication.

View larger version:

• In this window
• In a new window

• Download as PowerPoint Slide

Figure 6.

MERINGUE systematically and quantitatively evaluates for putative cell–cell communication for neuronal subtypes in the preoptic region using MERFISH data. (A) Distribution of −log₁₀(P-values) for the spatial cross-correlation between Cyp19a1 expression in neuronal subtypes and Esr1 expression in adjacent cells across animals. Boxes in the box plot denote the median values and inner quartile ranges (IQR) and whiskers denote 1.5 × IQR with additional outliers represented as points. Red dotted line is the alpha = 0.05 significance threshold. Generally, inhibitory neuron subtype I-2 shows significant spatial cross-correlation between Cyp19a1 expression and Esr1 expression in adjacent cells in a manner that is consistent across animals. (B) Cyp19a1 expression in I-2 neurons in one tissue slice in one animal with red indicating high expression and blue indicating low expression. Representative slice and animal shown. Select areas are highlighted in the zoom-in. (C) Esr1 expression in cells neighboring I-2 neurons in one tissue layer in one animal, with red indicating high expression and blue indicating low expression. Representative slice and animal shown. The same select areas as B are highlighted in the zoom-in.