STdGCN: spatial transcriptomic cell-type deconvolution using graph convolutional networks

Abstract

Spatially resolved transcriptomics integrates high-throughput transcriptome measurements with preserved spatial cellular organization information. However, many technologies cannot reach single-cell resolution. We present STdGCN, a graph model leveraging single-cell RNA sequencing (scRNA-seq) as reference for cell-type deconvolution in spatial transcriptomic (ST) data. STdGCN incorporates expression profiles from scRNA-seq and spatial localization from ST data for deconvolution. Extensive benchmarking on multiple datasets demonstrates that STdGCN outperforms 17 state-of-the-art models. In a human breast cancer Visium dataset, STdGCN delineates stroma, lymphocytes, and cancer cells, aiding tumor microenvironment analysis. In human heart ST data, STdGCN identifies changes in endothelial-cardiomyocyte communications during tissue development.

Background

The development of spatially resolved transcriptomic technologies has improved our understanding of the spatial organization and gene expression heterogeneity within tissues [1,2,3,4]. Spatial transcriptomics (ST) technologies enable the measurement of transcriptomes while maintaining the spatial context of the tissue, offering valuable insights into disease pathology research [5,6,7]. However, ST technologies still face certain limitations due to the trade-off between cell resolution and throughput. For instance, many fluorescence in situ hybridization (FISH) based technologies are unable to capture the full scope of the transcriptome within a single cell [8]. This poses computational challenges when attempting to infer the missing gene expression information from FISH data. In addition, many sequencing-based methods used in ST often lack sufficient resolution at the single cell level, raising computational challenges for deconvolving cell types detected within the same spatial spot/location. In contrast to ST technologies, single-cell RNA sequencing (scRNA-seq) is capable of detecting many thousands of gene expressions within individual cells. As a result, scRNA-seq can be effectively leveraged as a complementary tool to address certain limitations of ST technologies. For example, restoring cell type decomposition within a spot/location in ST data bears similarities to using scRNA-seq to deconvolve bulk RNA-seq data. Indeed, several models have been developed for bulk RNA-seq data deconvolution [9,10,11,12,13,14], and in theory, these models could be adapted to handle ST data. However, the cell numbers within a spot/location in ST data are considerably lower than those in bulk RNA-seq data. Consequently, applying a bulk RNA-seq deconvolution method to such a small sample size usually introduces noise from unrelated cell types and has been demonstrated by several benchmark analyses [15,16,17].

In recent years, several models have been developed for deconvolving ST data, employing either statistics-based or machine learning-based algorithms [3, 15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31]. Statistics-based algorithms typically assume a distribution, usually Poisson [18] or negative binomial [15, 17, 22, 26], for the gene transcript count. These methods first employ the scRNA-seq data to calculate the distribution coefficients for genes in each cell type. Cell-type proportions are then inferred using maximum-likelihood estimation or maximum a posteriori estimation. On the other hand, machine learning-based algorithms usually integrate the predicted ST data with scRNA-seq data to directly predict cell-type proportions through neural network systems. These models eliminate the need for an additional step to learn cell type gene signatures from scRNA-seq data [20, 21, 29]. However, a common limitation of these models is that they deconvolve the spots in a given ST dataset independently, without considering the neighborhood connections between spots. Indeed, it has been demonstrated that ST data can be divided into distinct spatial domains [32, 33]. By incorporating the spatial coordinate information between spots, the ST cell type decomposition learned from scRNA-seq can calibrated, leading to increased predictive accuracy.

In this study, we present the Spatial Transcriptomics deconvolution using Graph Convolutional Networks (STdGCN) as a novel approach for ST data cell-type deconvolution. Our model leverages graph convolutional networks (GCN), a widely used graph-based deep learning model. Notably, STdGCN integrates the expression profiles from scRNA-seq data with spatial localization information from the ST data for cell-type deconvolution. By integrating expression and spatial information, STdGCN achieves accurate predictions of cell-type proportions in ST data. To evaluate the performance of STdGCN, we conducted benchmarking against 17 state-of-the-art models. The results demonstrate that STdGCN consistently outperforms the benchmark models across various ST platforms. This highlights the superiority of STdGCN in cell-type deconvolution for ST data.

Results

STdGCN description

STdGCN is designed for cell-type deconvolution in ST data by using scRNA-seq data as a reference (Fig. 1). The underlying hypothesis is that both ST data and scRNA-seq data share common cell types, and their cell type gene transcript signatures exhibit similarities. In short, the normalized gene expression of a spot can be considered a combination of different cell types with varying proportions. The initial step of STdGCN involves identifying cell-type marker genes and generating pseudo-spots using the scRNA-seq data (Fig. 1A). Subsequently, it builds two link graphs to establish the GCN pipeline (Fig. 1B ~ C). The first link graph, known as the expression graph, is a hybrid graph comprising three sub-graphs, a pseudo-spot internal graph, a real-spot internal graph, and a real-to-pseudo-spot graph. Each sub-graph is constructed based on mutual nearest neighbors (MNN) using the expressed similarity between spots. The second link graph, the spatial graph, is built based on the Euclidean distance between real-spots in ST data. During the execution of STdGCN, the input feature matrix is forward propagated through the expression GCN layers and the spatial GCN layers, respectively. The outputs of the two GCN layers, namely Exp-feature and Spa-feature, are then concatenated column-wise into a single matrix. This concatenated matrix is subsequently fed into fully connected layers to predict the cell-type proportions for each spot. To train the model, we divide the pseudo-spots into a training dataset and a validation dataset. Only pseudo-spots in the training dataset are used for back propagation, whereas the validation dataset serves the purpose of early stopping. Using this approach, the cell-type proportions of real-spots can also be updated through the GCN pipeline, benefiting the learning of the pseudo-spots (Fig. 1D). The details are described in “Methods.”

Fig. 1

Schematic view of the STdGCN framework. STdGCN is a cell-type deconvolution method designed for spatial transcriptomics data, using scRNA-seq data as the reference. The workflow of STdGCN involves several key steps. Firstly, STdGCN employs the scRNA-seq reference data to identify cell-type marker genes and generate a pseudo-spot pool (A). It then builds two link graphs: a spatial graph (B) and an expression graph (C). The expression graph is a hybrid graph composed of three sub-graphs, a pseudo-spot internal graph, a real-spot internal graph, and a real-to-pseudo-spot graph. These sub-graphs are formed using mutual nearest neighbors (MNN) based on expression similarity. Based on the two link graphs, a GCN-based model is utilized to propagate information from both real- and pseudo-spots. The output of STdGCN is the predicted cell-type proportions for each spot (D). Additionally, STdGCN provides visualizations of the predicted results for the real-spot dataset in two different formats

Benchmarks of STdGCN in different datasets

To assess the performance of STdGCN, we conducted a comprehensive comparison with 17 state-of-the-art deconvolution models, including stereoscope [15], RCTD [18], SPOTlight [16], SpatialDWLS [19], Cell2location [26], DSTG [20], CellDART [21], DestVI [22], STRIDE [23], Tangram [24], BayesPrism [25], AdRoit [17], SpatialDecon [28], SD2 [29], NMFreg [3], CARD [30], and SONAR [34]. We selected three different ST platforms with single cell level cell type annotations: seqFISH [35], seqFISH + [33], and MERFISH [36]. Additionally, we conducted a methods comparison using slide-seq, a sequencing-based platform with near single-cell resolution (see Additional file 1: Supplementary Information) [37]. Synthetic multi-cellular ST datasets were generated using these platforms, allowing us to evaluate the performance of the models under controlled conditions (Fig. 2, “Methods”). Given that ST data and scRNA-seq data often originate from distinct datasets, we incorporated both internal and external single-cell reference of the seqFISH + dataset (“Methods”). To quantitatively assess accuracy of the 18 models, we employed three metrics: Jensen-Shannon divergence (JSD), root-mean-square error (RMSE), and Spearman correlation coefficient. Based on the three values, we ranked the models accordingly. Since each dataset consisted of multiple slices, the performance of a model for a given dataset was determined by calculating the mean rank across all slices within that dataset.

Fig. 2

Four representative multi-cellular synthetic benchmark ST slices corresponding to the four ST platforms. The scatter plots display the cell atlas of the mouse embryos sequenced with seqFISH (A), mouse preoptic region sequenced with MERFISH (C), mouse somatosensory (SS) region sequenced with seqFISH + (E), and mouse testis region sequenced with slide-seq (G). We resampled these datasets by dividing each slice into multiple square pixel areas (dash lines). Cells within each square pixel area were then merged into synthetic spots. Pie plots illustrate the cell-type proportions for each synthetic spot after resampling process (B, D, F, and H)

Table 1 illustrates the performance of the 18 models. STdGCN emerged the best outcomes, displaying the lowest average JSD in seqFISH, seqFISH + internal, and MERFISH datasets, as well as the lowest average RMSE in seqFISH and seqFISH + internal datasets. Specifically, in the seqFISH and seqFISH + datasets, STdGCN achieved the top ranking in both JSD and RMSE across all slices. Although STdGCN did not surpass other methods in Spearman’s correlation coefficient, its performance remained competitive. In the evaluating of seqFISH + using the external scRNA-seq dataset, STdGCN achieved 2nd place in both JSD and RMSE metrics, and 4th place in Spearman’s correlation coefficient. Despite not the best one, STdGCN consistently demonstrated stable performance across all benchmarked methods.

Table 1 The average ranks of the 18 benchmarked models across the three datasets

To further evaluate the performance of the models, we investigated their effectiveness in spots with different cell numbers. We divided the spots into two categories: a smaller group consisting of spots with five or fewer cells and a larger group comprising spots with more than five cells. We then benchmarked the models in each of the two groups across all datasets. For most models, the average JSD and RMSE of the larger group were lower than those of the smaller group (P value < 0.05, Wilcoxon test), and the Spearman’s correlation coefficient of the larger group were higher than those of the smaller group (P value < 0.05, Wilcoxon test) (Fig. 3A ~ D, Additional file 1: Figs. S1 ~ S2). The observations indicate that spots with smaller cell numbers pose greater challenges for accurate deconvolution. Despite this difficulty, when comparing the relative performance between two groups, the average ranks of most methods remained consistent (Additional file 2: Tables S1 ~ S3). Remarkably, STdGCN consistently outperformed other models in both groups.

Fig. 3

Assessment of STdGCN and benchmarked methods for estimation of cell-type proportions. A, B, and C The bar plots compare the performance of models with different numbers of cells within a spot. We divided the spots into two groups: a smaller group (spots with ≤ 5 cells) and a larger group (spots with > 5 cells). Bar plots display the Jensen-Shannon divergence (JSD) of the benchmarked methods for the two groups, as well as the total spots across the three datasets. D, E, and F Comparison of the models in spots with different number of cell type mixtures. We divided the spots into four groups based on the number of cell types within a spot: singleton (one cell type), doubleton (two cell types), tripleton (three cell types), and multiton (four or more cell types). Bar plots display the JSD of the benchmarked methods for the four groups, as well as the total spots across the three datasets. G, H, and I Bar plots display the JSD of STdGCN and benchmarked methods for distinguishing diverse cell types across the three datasets

Cell–cell communication plays a vital role in multi-cellular organisms, and it is influenced by physical location within the cellular microenvironment. Spatially adjacent cells are more likely to interact with each other. Thus, accurate deconvolution of cell types enables the identification of co-localization and spatial interactions. In light of this, we scrutinized the performance of the benchmarked models when dealing with different number of cell types within a single spot. We divided the spots into four groups based on the number of cell types: singleton (one cell type), doubleton (two cell types), tripleton (three cell types), and multiton (four or more cell types). Of the four groups, STdGCN consistently demonstrated strong performance (Fig. 3E ~ H, Additional file 1: Figs. S3 ~ S4, Additional file 2: Tables S4 ~ S6), and in most circumstances, STdGCN achieved the lowest JSD and RMSE scores.

Lastly, we assessed the models’ capability for predicting diverse cell types. For this analysis, we employed a binary JSD to measure the discrepancy between the predicted proportions and the actual proportions by summarizing the proportions of the remaining cell types within each spot (“Methods”). To measure the RMSE/Spearman’s correlation coefficient of a cell type, we exclusively calculated the RMSE/Spearman’s correlation coefficient between predicted proportions and the real proportions across all spots corresponding to that particular cell type within a given slice. In terms of the average ranks of the 18 models, STdGCN outperformed other models for most cell types (Fig. 3I, Additional file 1: Figs. S5 ~ S6). It is noteworthy that based on the average ranks of all cell types, STdGCN demonstrated a greater specialization in predicting cell types with higher frequencies (Additional file 2: Tables S7 ~ S9).

Developing human heart spatial transcriptomic dataset

In order to show the utility of STdGCN, we first employed STdGCN on a developing human heart [38] to assess its performance on a real-word multi-cellular ST dataset. This study used a combination of spatial transcriptomics [39, 40], scRNA-seq, and in situ sequencing (ISS), to capture transcriptomic profiles in human embryonic heart tissues (Fig. 4A ~ D). Leveraging the scRNA-seq data as the reference, we employed STdGCN to infer the cell-type composition within the spatial transcriptomics dataset. Remarkably, the cell-type distribution predicted by STdGCN exhibited a high degree of concordance with the ISS dataset (Fig. 4A, Additional file 1: Fig. S7). STdGCN accurately identified the spatial localization patterns of various cell types within the heart tissue: epicardial cells were predominantly located on the outer surface of the heart tissue, atrial cardiomyocytes were primarily situated within the atrium, and smooth muscle cells were located in the aorta and pulmonary artery area (Fig. 4E ~ G). These predicted cell types were consistent with existing medical knowledge [38], providing compelling evidence for the efficacy of STdGCN in capturing the spatial organization of cell types within complex tissue environments.

Fig. 4

Performance of STdGCN on the developing human heart spatial transcriptomic dataset. A The spatial structure and cell-type distribution of cells identified by in situ sequencing (ISS) in a PCW6.5 heart tissue section. B UMAP plot of the cells as in A. C The pie plot displays the predicted cell-type proportions in each spot of the PCW6.5–6 heart tissue section. D UMAP plot of the cells from the scRNA-seq dataset. E, F, and G The scatter plot displays the predicted proportions of epicardium cells (E), atrial cardiomyocytes (G), and smooth muscle cells from the PCW6.5–7 heart tissue section and their annotated results in ISS data. H and I The scatter plots display the annotated results of Myoz2-enriched cardiomyocytes (H) and ventricular cardiomyocytes (I) in ISS data and the predicted proportions from STdGCN. J The Venn diagram displays the overlapping of cell-type marker genes between ventricular cardiomyocytes, Myoz2-enriched cardiomyocytes, and atrial cardiomyocytes after filtering (see “Methods”). K The heatmap displays the co-localization score for all deconvolved cell types in the three stages (PCW4.5–5, PCW6.5, PCW9) of heart tissue sections. Red arrows are cardiomyocytes and blue arrows are endothelial cells. M The distributions of ventricular cardiomyocytes in the ISS dataset

In a previous review study [41], various models were compared in their ability to deconvolve the aforementioned dataset. The finding revealed that most models struggled to distinguish between Myoz2-enriched cardiomyocytes and ventricular cardiomyocytes (Fig. 3 in Chen et al. [41]) (Fig. 4H ~ I, Additional file 1: Figs. S8 ~ S9). It has been established that that the Myoz2-enriched cardiomyocytes are a subpopulation of cardiomyocytes in the healthy heart, characterized by the expression of MYOZ2 [38, 42]. Notably, both the UMAP plots of the ISS dataset and the scRNA-seq dataset demonstrated that the expression level of Myoz2-enriched cardiomyocytes closely resembled that of ventricular cardiomyocytes and atrial cardiomyocytes (Fig. 4B and D).

We investigated their respective cell-type marker genes (see “Methods”) and found that 41.6% of the Myoz2-enriched cardiomyocyte marker genes were also shared with atrial cardiomyocyte marker genes, while 100% of the Myoz2-enriched cardiomyocyte marker genes overlapped with ventricular cardiomyocyte marker genes (Fig. 4J). Notably, MYOZ2, which has been identified as the most significant marker gene for Myoz2-enriched cardiomyocytes (average log-fold change = 1.81, expressed in 99.0% Myoz2-enriched ventricular cardiomyocytes, FDR ≈ 1.08 × 10⁻⁸⁹) [38, 42], also serves as a marker gene for ventricular cardiomyocytes (average log-fold change = 1.30, expressed in 92.7% ventricular cardiomyocytes, FDR ≈ 0) (“Methods”). The considerable overlap of marker genes between ventricular cardiomyocyte and Myoz2-enriched cardiomyocyte elucidates why the majority of models fail to effectively differentiate between the two cell types. Despite the challenge, we compared STdGCN with eight previously well-performed models. The results showcased that STdGCN outperformed other models, successfully distinguishing between the two cell types (Fig. 4H ~ I, Additional file 1: Fig. S5).

Studies have emphasized the crucial role of communication between endothelial cells and cardiomyocytes in cardiac development. Endothelial cells not only facilitate the transmission of oxygenated blood supply to cardiomyocytes, but also provide local protective signals that promote cardiomyocyte organization and survival [43]. To investigate the extent of endothelial-cardiomyocyte communication within the analyzed dataset, we explored the cell type co-localization score (“Methods”) across three developmental stages of heart tissue sections (PCW4.5–5, PCW6.5, PCW9) (Fig. 4K). The background co-localization score, representing the mean pairwise co-localization score across all cell types, for the three stages are 0.0927, 0.0784, and 0.0845, respectively. The heatmap shows a higher co-localization score between endothelial cells (including “Endothelium/pericytes/adventitia” and “Capillary endothelium” in the dataset) and cardiomyocytes (including “Myoz2-enriched cardiomyocytes,” “ventricular cardiomyocytes,” “atrial cardiomyocytes” in the dataset) compared to the background co-localization score. This finding aligns with previous research and underscores the significance of endothelial-cardiomyocyte interactions in the context of this dataset.

Developing human breast cancer transcriptomic dataset

To assess the versatility of STdGCN, we then applied it to investigate the spatial organization of a human breast cancer ST dataset [44]. This dataset involved the analysis of breast tumor samples using both scRNA-seq (Fig. 5A ~ B) and 10X Genomics Visium ST sequencing. By utilizing the scRNA-seq data as a reference, we obtained the cell-type proportions corresponding to each spot in the 10X Genomics Visium dataset (Additional file 1: Figs. S10 ~ S15).

Fig. 5

Performance of STdGCN on the human breast cancer 10X Genomics Visium dataset. A and B UMAP plots of the cells from the scRNA-seq dataset divided by patient and cell type. C The heatmap displays the co-localization score (“Methods”) for normal cells, cancer cells, lymphocytes (integration of B-cells, T-cells, and plasmablasts), stroma (integration of cancer-associated fibroblasts, perivascular-like and endothelial cells), and myeloid. D Mean predicted cell-type proportions for patient CID4465 (TNBC) and CID4535 (ER +). E, F, and G The scatter plot displays the predicted result for patient CID4535. H The heatmap displays the co-localization score for patient CID4535. I, J, and K The scatter plot displays the predicted result for patient CID4465. L The heatmap displays the co-localization score for patient CID4465

By using STdGCN, we were able to quantitatively dissect the tumor microenvironment (TME) for these tissues in high resolution. For example, previous studies have highlighted the significance of tumor-infiltrating lymphocytes (TILs) as a strong prognostic marker for patient survival. In the case of locally developed breast cancer treated with neoadjuvant chemotherapy, the presence of TILs is a predictor of the response [45,46,47]. Leveraging the deconvolution results predicted by STdGCN, we compared the cell type co-localization score (“Methods”) between estrogen receptor-positive (ER +) samples and triple-negative breast cancer (TNBC) samples. The co-localization heatmap revealed higher TILs in ER + samples than TNBC samples, suggesting a more favorable prognosis for ER + samples (Fig. 5C). Additionally, the interactions between stromal cells, cancer cells, and lymphocytes have been reported to be important for maintaining the TME [48]. Cancer cells that interact with stromal cells can create favorable physical or molecular signaling environments that promote tumor progression [49], while certain stromal cells, such as cancer-associated fibroblasts (CAFs), can recruit immune cells to the tumor tissues and influence their behavior towards cancer cells [50]. Through STdGCN, we were able to analyze cell–cell interactions in spatial locations. For instance, we compared the TME of an ER + patient (CID4535) with that of a TNBC patient (CID4465) (Fig. 5D ~ L). Despite similar overall proportions of cell types between the two patients (Fig. 5D), the spatial structure and co-localization heatmap revealed substantial differences in their respective (Fig. 5E ~ L), indicating distinct microenvironmental characteristics. Both patients have three regions on their slices, and we analyzed the TME variations across these regions. The average proportions of cell types and the co-localization heatmaps reveal significant differences in the TMEs within individual patients (Additional file 1: Fig. S16). For patient CID4465, regions 1 and 3 have smaller proportions of T-cells compared to region 2 and all regions of patient CID4535. Furthermore, all regions of patient CID4465 exhibit lower TILs compared to those in patient CID4535.

Discussion

ST profiles have significantly enhanced our understanding of tissue architecture at the cellular level. However, many existing ST profiling technologies lack sufficient single-cell resolution. To this end, we introduce STdGCN, a graph-based deep learning framework that leverages cell type profiles learned from single-cell RNA-seq to effectively deconvolve the cell type mixtures in ST data. In order to evaluate the performance of STdGCN, we conducted a comprehensive benchmarking study involving four ST datasets and compared its performance against 17 state-of-the-art ST deconvolution methods [3, 15,16,17,18,19,20,21,22,23,24,25,26, 28,29,30]. The models were ranked based on their performance metrics. Across all four benchmarked datasets, STdGCN demonstrated the most superior performance in terms of JSD, RMSE, and Spearman’s correlation coefficient (Table 1, Additional file 2: Table S10).

To improve the performance of cell-type deconvolution, distinct designs were adapted based on the characteristics of ST datasets. First, STdGCN stands out from previous models by incorporating spatial structure information from ST data and expression profiles from scRNA-seq for cell-type deconvolution. This innovative approach capitalizes on the inherent similarity between spots within a spatial domain, allowing for precise calibration of predictions based on the expression profiles derived from scRNA-seq data. Although STdGCN is not the first study to utilize a GCN module for cell-type deconvolution, its novel architecture stands out by building the GCN pipeline on two link graphs—an expression graph and a spatial graph (Fig. 1B ~ C). This approach differentiates it from previously published studies. Next, to maximize the heterogeneity of the pseudo-spot pool, we developed a sampling-by-cell-type method to simulate pseudo-spots. This approach does not rely on cell-type distributions from scRNA-seq by assigning equal weights for each cell type. This effectively addresses the challenge of rare cell types being unlikely to be randomly selected using traditional methods. In addition, we imposed a maximum limit on the number of cell types within a single pseudo-spot to accurately mimic the characteristics observed in real-word ST data. To evaluate the extent of the impact of these approaches on the performance of STdGCN, we benchmarked the complete STdGCN model against STdGCN utilizing the conventional pseudo-spot generation approach and STdGCN without considering the spatial graph in deconvolving cell types of the seqFISH + dataset (Additional file 2: Table S11). Among the three methods, the complete STdGCN model achieved the best performance. It is noteworthy that, compared with the spatial structure information, a better pseudo-spot generation strategy offered a much greater improvement to STdGCN.

The modular design of STdGCN provides users with a high degree of flexibility. While default parameters are provided, STdGCN allows users to fine-tune the model according to the specific characteristics of their ST platforms. For example, in the preprocessing module, users have the option to input either a raw count matrix or a normalized expression matrix from both ST data and scRNA-seq data. Additionally, users can utilize their own genes of interest to customize the deconvolution process. The pseudo-spot simulation module offers users the flexibility to select the desired cell number and cell-type number within a spot, adapting to the specific characteristics of their ST platform. The number of simulated spots is also adjustable, allowing users to strike a balance between improved performance (i.e., more simulated spot number) and computational efficiency and memory usage. In the GCN module, users have the freedom to customize various aspects of the model, such as the weights of the link graphs, the depth of expression/spatial GCN layers, the depth of the fully connected neural network layers, and the dimensionality of the hidden layers. Lastly, STdGCN provides output in both csv table and the anndata format [51]. This enables users to seamlessly integrate the results with numerous downstream analysis tools, such as Scanpy, which are widely accessible and user-friendly.

In addition to evaluating model performance in deconvolution, we also consider computational efficiency as an important factor affecting user experience. In this study, we implemented both parallel computing and GPU acceleration to enhance the computational efficiency of STdGCN. We conducted a comprehensive comparison of running times across various benchmarked methods for deconvolving two datasets: the cortex seqFISH + slice using an external scRNA-seq reference and the developing human heart ST dataset (Additional file 2: Table S12). The cortex seqFISH + slice comprises only 109 spots and six cell types, while the developing human heart ST dataset consists of 3111 spots and 14 cell types. As anticipated, the running time for the developing human heart dataset is generally longer than that for the seqFISH + slice across most methods. Among these benchmarked methods, STdGCN exhibits consistent overall running times across both datasets. Although it may run slower than most methods on datasets with fewer spots and cell types, it significantly reduces processing time compared to many benchmarked methods when dealing with datasets containing more spots and cell types.

Despite the aforementioned merits, STdGCN still has limitations. One limitation of STdGCN is that it simply concatenates the real-spot internal graph, based on ST spots’ expression levels, with the spatial graph derived from spots’ coordinates during the deconvolution process. When a spot is located at the boundary of spatial domains, message passing from neighboring spots, especially those from different spatial domains, could potentially carry negative impacts. This can be mitigated by employing spatial clustering or spatial domain identification to refine the adjacency matrix, thereby diminishing inaccuracies or noise caused by boundaries. The second limitation is that the predictive performance of STdGCN varies across cell types. In particular, the relative performance of STdGCN in cell types with lower proportions is not as good as in the cell types with higher proportions (Additional file 2: Tables S4 ~ S6). This indicates that there is room for improvement in the current pseudo-spot generation method and the construction model of the link graph. Thus, our future work will focus on ameliorating the effectiveness of utilizing reference data.

In conclusion, STdGCN stands as the pioneering model that incorporates the spatial coordinates from ST data and cell type profiles learned from single-cell RNA-seq data for the purpose of deconvolving local cell-type proportions within complex tissues. By leveraging this innovative approach, STdGCN has the potential to greatly enhance our understanding of tissue spatial architecture and provide invaluable support for downstream analyses at the cellular level. With its ability to unravel the intricate composition of cell types within ST datasets, STdGCN opens new avenues for exploring the complexities of tissue organization and advancing our knowledge of biological systems.

Conclusions

Spatially resolved transcriptomics combines high-throughput measurement of transcriptomes with the preservation of spatial information about cellular organizations. In this study, we introduce STdGCN, a graph neural network model tailored for cell-type deconvolution of ST data, capable of leveraging abundant scRNA-seq data as a reference. STdGCN integrates expression profiles from single-cell data and spatial localization information from ST data to facilitate cell-type deconvolution. Extensive benchmarking experiments across multiple ST datasets demonstrated that STdGCN outperformed 17 cutting-edge models. In the analysis of a human breast cancer Visium dataset, STdGCN effectively distinguished spatial distributions between stroma, lymphocytes, and cancer cells, facilitating tumor microenvironment dissection. Furthermore, in a human heart ST dataset, STdGCN identified changes in potential endothelial-cardiomyocyte communications during tissue development. In conclusion, STdGCN represents a powerful tool for unraveling the intricate interplay between cellular composition and spatial organization within complex tissue environments. This framework holds great promise for advancing our understanding of spatially resolved transcriptomics and its implications for various biological processes and disease states.

Methods

Benchmark data selection and preprocess

To curate the datasets for benchmark, we conducted a thorough review of two published comparison studies on ST deconvolution [41, 52], as well as examined the models benchmarked in this study. After strict comparison process, we selected four datasets from four ST platforms, including seqFISH, MERFISH, seqFISH + , and slide-seq (Additional file 1: Table S13 [53,54,55]). These datasets were chosen because they offer single cell-level resolution. Moreover, the selected datasets come with cell type annotations provided by the authors, enabling us to obtain the ground truth after simulating low-resolution spots. In addition, the selected datasets include both sequencing-based technologies and imaging-based technologies, with gene sizes ranging from hundreds to genome wide (> 20,000) and the spot size ranging from tens to thousands. This diverse range of characteristics allows for a comprehensive evaluation of the models in various scenarios. Importantly, all four datasets contain multiple slices, which facilitates testing the robustness of the models across different samples. It is worth noting that the resolution of slide-seq approximates that of a single cell, but it is not at strictly single-cell resolution. Despite these limitations, it provides gene expression information at a genome-wide level, complementing other benchmark datasets. Consequently, we have chosen to utilize slide-seq data and regard it as a valuable resource for benchmarking spatial transcriptomic methodologies, and to avoid controversy, we put the benchmark of slide-seq data in Additional file 1: Supplementary Information.

To create synthetic spots for analysis, we resampled the datasets by dividing each slice into multiple square pixel areas. Cells within each square pixel area were merged to form a synthetic spot. As the ST platforms used provided single-cell-level data, we had access to the exact cell-type proportions for each synthetic spot. The single-cell-level data was also used as the internal single-cell reference in our study.

The seqFISH + dataset was obtained from Zhu et al. [33] (access via Giotto package [56]: getSpatialDataset(dataset = “seqfish_SS_cortex”, method = “wget”)), which encompasses the detection of mRNAs for 10,000 genes corresponding to a section of the mouse somatosensory (SS) region. We divided the SS region into two slices based on the author’s annotation, including a cortex slice and a sub-ventricular zone (SVZ) slice. For the cortex slice, we selected 400 × 400 square pixel area as a synthetic spot, whereas for the SVZ slice we chose 200 × 200 square pixel area. Any spot with cell number below two was discarded. There were 109 cortex synthetic spots and 59 SVZ synthetic spots after resampling, respectively. Given that ST data and scRNA-seq data often originate from distinct datasets, we opted to incorporate an external single-cell reference from the SpatialDWLS study [19] for benchmarking purposes. This approach enabled us to thoroughly evaluate the performance of the benchmarking models.

The seqFISH dataset was obtained from Lohoff et al. [35] (https://content.cruk.cam.ac.uk/jmlab/SpatialMouseAtlas2020/). This study detected mRNAs for 351 target genes in tissue sections of three mouse embryos at the 8–12 somite stage, with two slices obtained for each embryo. In terms of the cell types annotated by the authors from single-cell transcriptome atlases, we removed the cell types with very small proportions (< 3%) for cell-type deconvolution. The coordinates of the cells were stored in columns “x_global_affine” and “y_global_affine” within the file “metadata.Rds.” For each slice, we selected the 0.1 × 0.1 square “global_affine” area as a synthetic spot. Any spot with cell number less than two was discarded. The number of predictive synthetic spots obtained for each slice after applying the filtering criteria ranged from 1279 to 2350.

The MERFISH dataset was downloaded from Moffitt et al. [36, 57] (https://datadryad.org/stash/dataset/doi:https://doi.org/10.5061/dryad.8t8s248). This study detected mRNAs for a set of 155 target genes in the mouse preoptic region. Approximately 1.1 million cells were measured across 181 slices obtained from 36 mice. Cell types with proportions smaller than 3% were removed for cell-type deconvolution. The coordinates of the cells were stored in columns “xcoord” and “ycoord” in the file named “Moffitt_and_Bambah-Mukku_et_al_merfish_all_cells.csv.” For each slice, we selected a synthetic spot size of 50 × 50 square pixel area. Any spot with cell number less than two was discarded. The final number of predictive synthetic spots for each slice after filtering ranged from 938 to 1247.

Preprocess of the human heart dataset

The human heart spatial transcriptomics dataset was obtained from Asp et al. [38]. This dataset contains three sequencing platforms: spatial transcriptomics [39, 40], scRNA-seq, and ISS. The tissue sections were collected from human embryonic cardiac samples at three developmental stages in the first trimester: 4.5–5, 6.5, and 9 post-conception weeks (PCW). Herein, we employed STdGCN to infer the cell-type composition of the spatial transcriptomics data. This dataset consists of 19 slices, with four slices from 4.5–5 PCW, nine slices from 6.5 PCW, and six slices from 9 PCW. During data deconvolution, we considered only the spots located within the tissues. The coordinates of the spots can be obtained from their corresponding spot names. We used the scRNA-seq dataset provided in this study as the reference data. It is important to note that the spatial transcriptomics dataset lacks single-cell resolution, making it impossible to obtain realistic cell-type compositions for the spots. Additionally, this study offers the differentially expressed gene (DEG) analysis results based on the scRNA-Seq data (Table S3 in Asp et al. [38]). To filter cell-type marker genes for Myoz2-enriched cardiomyocytes, ventricular cardiomyocytes, atrial cardiomyocytes, we retained genes with FDR (marked as “'p_val_adj'”) < 0.01, average log-fold change (marked as “'avg_logFC'”) > 1, fraction of cells expressing the genes in the cell type (marked as “pct.1”) > 0.7, and fraction of cells expressing the genes in the remaining cell types (marked as “pct.2”) < 0.3 from the downloaded DEG analysis table. The results following the filtration of Myoz2-enriched cardiomyocytes, ventricular cardiomyocytes, and atrial cardiomyocytes are depicted in Fig. 4J.

Preprocess of the human breast cancer dataset

The human breast cancer 10X Genomics Visium dataset was obtained from Wu et al. [44]. The dataset consists of data from six breast cancer patients, including two ER⁺ samples and four TNBC samples. When performing data deconvolution, we considered only the spots located within the tissues. The coordinates of the spots can be found in the file “tissue_positions_list.csv.” Additionally, in the same study, the authors processed scRNA-seq from 26 breast tumor samples. For the purpose of deconvolution in our study, we employed this scRNA-seq data as the reference dataset.

Benchmark metrics

To assess the effectiveness of cell-type deconvolution algorithms across the synthetic real-word datasets, we selected three metrics: JSD, RMSE, and Spearman’s rank correlation coefficient. These metrics were chosen to quantify the disparities between the predicted cell-type proportions and the actual proportions within each spot.

The JSD is a symmetric form of the Kullback–Leibler (KL) divergence. Denote $P$ and $Q$ are the probability distributions representing the predicted and ground-truth results, respectively. The JSD can be expressed as:

$$\text{JSD}(P|\left|Q\right)=\frac{1}{2}{D}_{KL}(P||\frac{1}{2}(P+Q))+\frac{1}{2}{D}_{KL}(Q||\frac{1}{2}(P+Q))$$

(1)

where ${D}_{KL}$ is the Kullback–Leibler (KL) divergence, which can be expressed as:

$${D}_{KL}(P|\left|Q\right)=\sum_{x\in X}P\left(x\right)\text{log}\frac{P\left(x\right)}{Q\left(x\right)}$$

(2)

Typically, the JSD is a nonnegative number, and a smaller JSD represents a higher similarity between the predicted result and the ground-truth, indicating a better performance of the corresponding method in accurately estimating the cell-type proportions.

To assess the performance of the models specifically for cell type $A$, denote the predicted and ground-truth proportions of a spot are $P(A)$ and $Q(A)$. For the purpose of analysis, we considered the remaining cell types as a group, collectively referred to as the non-$A$ type. Their proportions are denoted as $1-P(A)$ and $1-Q(A)$, respectively. Equation (2) will be expressed as:

$${D}_{KL}(P|\left|Q\right)=P\left(A\right)\text{log}\frac{P\left(A\right)}{Q\left(A\right)}+(1-P(A))\text{log}\frac{\left(1-P\left(A\right)\right)}{\left(1-Q\left(A\right)\right)}$$

(3)

Based on Eq. (3), we can eventually get the binary JSD between predicted proportions compared to the ground truth of the cell type A.

We used the following equation to calculate RMSE:

$$\text{RMSE}(P|\left|Q\right)=\sqrt{\frac{\sum_{i=1}^{K}{(P\left({x}_{i}\right)-Q\left({x}_{i}\right))}^{2}}{K}}$$

(4)

When evaluating the performance of a model for deconvolving a spot, the variable $K$ represents the number of cell types. Conversely, when assessing the performance of the models for a specific cell type, $K$ represents the number of spots within a slice. RMSE is a nonnegative value, and a smaller RMSE represents a lower dissimilarity between the predicted result and the ground truth.

The Spearman’s rank correlation coefficient can be expressed as:

$${r}_{s}=1-\frac{6\sum {d}_{i}^{2}}{n({n}^{2}-1)}$$

(5)

where ${d}_{i}= P\left({x}_{i}\right)-Q\left({x}_{i}\right)$ represents the difference between the two ranks of each cell type, and n is the number of cell types. The value ${r}_{s}$ ranges from − 1 to 1, with values close to 1 indicating similarity in ranks for the predicted value and ground truth, and values close to − 1 indicating dissimilarity in ranks between the predicted value and the ground truth.

Cell type co-localization score

Cell–cell interactions often occur via receptor–ligand pairs. Thus, identifying cell type co-localization provides potential cell–cell interactions within a tissue. To quantify the degree of co-localization between two cell types within a ST slice, we have developed a cell type co-localization score denoted as ($\text{S}\left(\text{A},\text{ B}\right)$). $\text{S}\left(\text{A},\text{ B}\right)$ can be calculated using the following functions:

$$\text{S}\left(\text{A},\text{ B}\right)=\frac{{\sum }_{all spots}{S}_{A}=1 \text{and} {S}_{B}=1}{{\sum }_{all spots}{S}_{A}=1 \text{or} {S}_{B}=1}$$

where ${S}_{A}=1$ denotes that the location has cell type A. The value of $\text{S}\left(\text{A},\text{ B}\right)$ ranges between 0 and 1. A higher $\text{S}\left(\text{A},\text{ B}\right)$ a higher probability of co-localization between the two cell types within a spot. Considering that the proportions of all cell types within a spot predicted by STdGCN are greater than 0, to ensure meaningful co-localization analysis, we introduce a threshold for $\text{S}\left(\text{A},\text{ B}\right)$ calculations. When the predicted proportion of cell type $A$ within a spot is below the threshold, we set ${S}_{A}=0$. We used a threshold of 0.15 for the 10X Genomics Visium human breast cancer dataset, and 0.1 for the developing human heart spatial transcriptomic. This choice was made based on the characteristics of the dataset and the platform used.

Implementation of STdGCN

Single-cell data preprocessing and feature selection

The first step in the processing pipeline of STdGCN is to identify the cell-type marker genes from scRNA-seq data. The input of the single-cell expression matrix is flexible and can be used for either the raw count matrix or the normalized matrix. For the raw count matrix, we used the Scanpy package in Python to preprocess the data. We used the command “scanpy.pp.normalize_total” to normalize counts to a standardized value of 10,000 per cell. Subsequently, we logarithmized the data matrix using the command “sc.pp.log1p,” followed by “scanpy.pp.scale” to scale data to achieve unit variance and a mean of zero. During the preprocessing phase, we have incorporated options to allow users to make decisions regarding the filtering of highly variable genes and the regression of mitochondrial genes. These options provide flexibility based on specific analysis requirements. Following preprocessing step, we used the principal component analysis (PCA) (sc.tl.pca) for dimensionality reduction. To select cell-type marker genes for each cell type, we utilized logistic regression (sc.tl.rank_genes_groups(method = ‘logreg’)). The resulting gene features were obtained by aggregating the top “$n$” cell-type marker genes, where “$n$” represents the desired number of markers to be considered. To streamline the preprocessing and marker gene selection steps, we have integrated them into a single command named “find_marker_genes.” This integration enhances the ease of use and efficiency of the STdGCN pipeline.

Pseudo-spot generation

The generation of pseudo-spots involves randomly selecting a mixture of cells from scRNA-seq data. In our experience, the heterogeneity of the pseudo-spot pool is a key factor in constructing an effective GCN link graph. Conventional methods for pseudo-spot generation may result in an unbalanced pseudo-spot pool and limited heterogeneity, particularly due to the unbalanced proportions of cell types in reference single-cell data. To address this issue, we designed a sampling-by-cell-type method for obtaining a more balanced pseudo-spot generation. Simulating pseudo-spots requires specifying three user-defined hyper-parameters: ${N}_{up}$, the upper limit of cell numbers in a pseudo-spot; ${N}_{low}$, the lower limit of cell numbers in a pseudo-spot; and ${K}_{up}$, the upper limit of the number of cell types in a pseudo-spot. The use of ${K}_{up}$ is influenced by the fact that the number of cell types within a spot is usually limited in real-word ST datasets, unlike bulk sequencing data. For each pseudo-spot simulation, we first generated a random number between 1 and ${K}_{up}$ to determine the number of cell types ($k$) within a pseudo-spot. Subsequently, we randomly selected $k$ cell types from the entire cell types (${K}_{w}$) to form the cell-type set for the corresponding pseudo-spot. We next randomly selected a number (${N}_{s}$) between ${N}_{low}$ and ${N}_{up}$ to determine the cell number in this pseudo-spot. We then randomly picked up ${N}_{s}$ cells from the selected cell-type set to construct the pseudo-spot. The expression level of the pseudo-spot was determined by the summation of the counts of the randomly selected cells. To streamline the entire process of pseudo-spot generation, we have implemented a single command called “pseudo_spot_generation” within the STdGCN framework.

Real-spot and pseudo-spot preprocessing and integration

The preprocessing steps for real-spots and pseudo-spots were performed in a similar manner to the reference single-cell data. After preprocessing, we retained the previously selected marker genes for downstream analyses. To facilitate data integration, we employed PCA for dimensionality reduction, allowing the mapping of real- and pseudo-spots into a common latent space. The integrated dimensionality reduction matrix was used as the input for constructing the subsequent link graph.

Link graph construction

The STdGCN system includes two link graphs, an expression graph and a spatial graph. The expression graph captures the expression similarity between spots. It is a hybrid graph composed of three sub-graphs, a pseudo-spot internal graph, a real-spot internal graph, and a real-to-pseudo-spot graph. We employed the “cosine” distance to measure the distances between spots and the MNN [58] algorithm to establish links between spots. In MNN model, a link between two spots ($a$ and $b$) if and only if spot a belongs to one of the nearest neighbors of spot $b$, and spot $b$ is also one of the nearest neighbors of spot $a$. For the real-to-pseudo-spot link graph, we applied the integrated dimensionality reduction matrix of the two datasets as the input feature, with the parameter “nearest neighbors = 20” per cell. For the two internal link graphs, data integration was not required, we thus directly used the normalized expression matrix of the real-spot/pseudo-spot as the input feature. The parameter “nearest neighbors” was set to 10 for the real-spot internal graph and 20 for pseudo-spot internal graph. The adjacency matrix of the hybrid expression graph (${A}^{exp}$) is a binary matrix. ${A}_{ab}^{exp}=1$ represents a link between two spots and vice versa.

The spatial graph (${A}^{sp}$) represents the spatial distances between real-spots. We calculated the Euclidean distance between a given spot and all other spots (${D}_{ab}$) using their relative spatial coordinates. In our approach, we defined the distance between two neighboring spots as one unit (= 1). To construct the spatial graph, we established links between spots that have distances smaller than a predefined radius (${D}_{upper}$). The adjacency matrix of the spatial graph is defined as follows:

$$A_{ab}=\left\{\begin{array}{lc}0,&\text{when}\;D_{ab}\geq D_{upper}\\1/D_{ab},&\text{when}\;D_{ab}< D_{upper}\end{array}\right.$$

Unlike the expression link graph, which is represented as a binary matrix, the spatial graph is a weighted link graph based on the spatial distance between spots.

Model description

The STdGCN is built on the GCN [59], a semi-supervised learning algorithm designed for graph-structured data. GCN is based on an efficient variant of convolutional neural networks that directly operate on graphs. Denoting $\text{V}$ and $\text{E}$ are the sets of nodes and edges in the graph $\text{G}=(\text{V},\text{ E})$, the objective of GCN is to learn the graph representations by combining node embeddings with those of their neighboring nodes. The model takes two components as input: a feature matrix ($\text{X}$) representing the embeddings of all graph nodes, and an adjacency matrix ($\text{A}$) depicting the structural connections in the link graph. Each GCN layer can be written as a non-linear function:

$${H}^{(l+1)}=f\left({H}^{l},\widetilde{A}\right)=\sigma (\widetilde{A}{H}^{l}{W}^{l})$$

where ${H}^{l}$ is the output for the lth layer with ${H}^{l}=X$, ${W}^{l}$ is the weight matrix for the lth layer that needs to be learned during the training, $\sigma (\cdot )$ is the activation function, and $\widetilde{A}={D}^{-1/2}(A/\lambda +I){D}^{1/2}$ is the symmetric normalized adjacency matrix of $A$. The $\lambda =20$ is a stabilizing factor, as we have observed that incorporating this factor can increase the accuracy and robustness of STdGCN. Here, $I$ represents the identity matrix and $D$ is the diagonal degree matrix of $(A/\lambda +I)$.

The STdGCN constructs a comprehensive and diverse graph, containing both pseudo-spots and real-spots as individual nodes, enabling the adaptation of real-spots from pseudo-spots. The objective of STdGCN is to predict the cell-type proportions for the real-spots. Therefore, the label for each spot is a vector that includes proportions for all cell types (${Y}_{prop}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times {n}_{type}}$), where ${n}_{r}$ and ${n}_{p}$ denote the real- and pseudo-spot numbers and ${n}_{type}$ denotes the cell-type number. We used the expression level for the real- and pseudo-spots as the input matrix ($\text{X}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times g}$), where g represents the selected gene number. The input matrix was then forward propagated to two two-layer GCNs using the two previously constructed link graphs, an expression graph (${A}^{exp}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times ({n}_{r}+{n}_{p})}$) and a spatial graph (${A}^{sp}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times ({n}_{r}+{n}_{p})}$), with the function:

$$\left\{\begin{array}{c}{X}_{exp}=ELU({\widetilde{A}}^{exp}ELU({\widetilde{A}}^{exp}X{W}_{exp}^{(0)}){W}_{exp}^{(1)})\\ {X}_{sp}=ELU\left({\widetilde{A}}^{sp}ELU\left({\widetilde{A}}^{sp}X{W}_{sp}^{\left(0\right)}\right){W}_{sp}^{\left(1\right)}\right)\end{array}\right.$$

to obtain an expression embedding (${X}_{exp}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times m}$) and a spatial embedding (${X}_{sp}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times m}$), where m is the dimension of the embedding space and $ELU(\cdot )$ is the activation function that can be expressed as:

$$ELU\left(x\right)=\left\{\begin{array}{c}x, \text{if }x>0\\ {e}^{x}-1, \text{if} x\le 0\end{array}\right..$$

We next concatenated the two GCN output embedding matrix (${X}_{con}=[{X}_{exp}, {X}_{sp}]\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times 2m}$) as the input feature and propagated it through two additional two-layer fully connected neural networks using the following function:

$${\widetilde{Y}}_{prop}=softmax({b}^{(1)}+{a}^{1}ELU({b}^{(0)}+{a}^{(0)}{X}_{con}))$$

where ${a}^{(l)}$ and ${b}^{(l)}$ indicate the weight matrix and additive bias for the lth layer, respectively, ${\widetilde{Y}}_{prop}\in {\mathbb{R}}^{({n}_{r}+{n}_{p})\times {n}_{type}}$ is the prediction of cell-type proportions, and $softmax(\cdot )$ is the activation function that can be expressed as:

$$softmax\left({x}_{i}\right)=\frac{{x}_{i}}{{\sum }_{j}{x}_{j}}.$$

We used the KL divergence with L2 penalty as the loss function during the training of STdGCN.

Parameter settings

Herein, we used a two-layer GCN architecture for the implementation of STdGCN, as we found that utilizing two GCN layers yielded superior performance compared to employing a single layer or multiple layers. To train STdGCN, we adopted stochastic gradient descent (SGD) with a maximum of 3000 epochs (Python command: torch.optim.SGD (model.parameters(), lr = 0.2, momentum = 0.9, weight_decay = 0.0001, dampening = 0, nesterov = True)). We implemented an early stopping mechanism where training would cease if the validation loss did not decrease for 20 consecutive epochs. We used a self-adaption method called orch.optim.lr_scheduler.ReduceLROnPlateau to adjust the learning rate during training. The reduced learning rate factor was 0.1 and the learning rate was reduced after 5 epochs without improvement.

Running time benchmark and computational resource

To evaluate the computational efficiency of cell-type deconvolution algorithms, we conducted a comprehensive comparison of running time across the benchmarked methods for deconvolving the cortex seqFISH + slice using the external scRNA-seq reference. To enhance computational efficiency, we employed parallel computing and GPU acceleration if these methods offered such options. The workstation used to test all methods was two Intel(R) Xeon(R) CPU E5-2650 v4 @ 2.20GHz (30 MB cache size; 24 cores in total) and 516 GB of memory. The GPU used was the Tesla V100S-PCIE-32GB. The operating system used was CentOS Linux 7.