Spatial transcriptomics is a rapidly developing field of single-cell genomics that quantitatively measures gene expression while providing spatial information within tissues. A key challenge in spatial transcriptomics is identifying spatially structured domains, which involves analyzing transcriptomic data to find clusters of cells with similar expression patterns and their spatial distribution. To address these challenges, we propose a novel deep-learning method called DMGCN for domain identification. The process begins with preprocessing that constructs two types of graphs: a spatial graph based on Euclidean distance and a feature graph based on Cosine distance. These graphs represent spatial positions and gene expressions, respectively. The embeddings of both graphs are generated using a multi-view graph convolutional encoder with an attention mechanism, enabling separate and co-convolution of the graphs, as well as corrupted feature convolution for contrastive learning. Finally, a fully connected network (FCN) decoder is employed to generate domain labels and reconstruct gene expressions for downstream analysis. Experimental results demonstrate that DMGCN consistently outperforms state-of-the-art methods in various tasks, including spatial clustering, trajectory inference, and gene expression broadcasting.

Spatial transcriptomics technology not only provides high-resolution gene expression data but also completely preserves the spatial location information of each sequencing spot in tissue sections, offering an unprecedented multi-dimensional perspective for in-depth exploration of tissue heterogeneity and dissection of cellular microenvironments and functional partitions. Similar to the cell clustering task in single-cell data analysis, one of the core challenges in spatial transcriptomics data analysis is spatial domain identification—using algorithms to cluster sequencing spots with adjacent spatial locations and similar gene expression patterns into biologically meaningful functional regions. Although existing methods can achieve preliminary spatial domain partitioning, accurately capturing global semantic associations and balancing the modeling capabilities of local features and global structures remain critical scientific challenges when faced with non-linear gene expression patterns and spatial distribution relationships in complex tissues. Here, we propose a novel algorithmic framework that integrates multi-view graph convolutional networks and contrastive learning. Results across multiple datasets generated by different technologies demonstrate that our method significantly enhances the accuracy and biological interpretability of spatial domain identification.

(a) Dataset preprocessing. The preprocessing includes the construction of the spatial graph, based on the Euclidean distance, and the feature graph, based on the Cosine distance. (b) and (c). Multi-view Graph Convolutional (MGCN) encoder. Spatial, feature, and their combination embeddings of both graphs will be generated by the attention mechanism as well, and corrupted features will be trained and discriminated before attention. (d) The result of DMGCN can be used for downstream analysis, such as domain identification and gene expression broadcasting.

https://doi.org/10.1371/journal.pcbi.1013369.g001

Model performance of domain identification in DLPFC and downstream analysis

To evaluate the spatial domain recognition performance of DMGCN, the 10x Visium human dorsolateral prefrontal cortex (DLPFC) dataset [24] proves our model’s capacity in domain identification. Through the experiment on the DLPFC dataset with label information (Fig 2A), the quantitative evaluation, ARI or NMI, was applied to test the algorithm’s performance in handling the complex tissue structures on the specific dataset.

In this experiment, we compared DMGCN with five developed methods that excel in domain identification, allowing us to illustrate the model’s capabilities in comparison. Across all 12 slices, we found that DMGCN achieved the highest performance among all methods on DLPFC (Fig 2B), especially on slices 151509, 151671, and 151675 (Fig 2C).

In addition, compare the violin maps (Fig 2D) between the original gene expression and the reconstructed gene expression, the significant enhancement of these mark genes was found, and this result proves the DMGCN’s capacity to reconstruct the gene expression in the right way.

Performance of domain identification in HBC and downstream analysis

Following the outstanding performance of the DMGCN model on DLPFC samples, this paper applies the DMGCN model to the human breast cancer dataset for evaluation. In this experiment, we found that, compared to the manual annotation, the results of the DMGCN model, in which the ARI is 0.68 (Fig 3A), show that the DMGCN model has strong accuracy in spatial domain recognition for breast cancer pathology samples and performs better than other methods (Fig 3C), whatever for ARI or NMI (Fig 3B). Except that, the UMAP visualization (Fig 3D) shows that the DMGCN is better at processing the correlation and trajectory of different spatial clusters.

Based on the result of DMGCN on human breast cancer, we select domains 3, 7, and 11 that correspond to the IDC_2, Healthy, and IDC_8 in ground truth, and we found that the differentially marker genes [12], such CRISP3, AQP3, VTCN1, NUPR1, CPB1, FCGR3B, CGA, REPS2, can match corresponding pathology areas well (Fig 4A). Except for these marker genes, we still found that the marker genes, LINC00645 and SLITRK6, can be broadcast through the reconstructed gene expression (Fig 4B).

Performance on Stereo-seq ST data and downstream analysis

To further assess the generalizability of DMGCN in spatial domain identification, we performed additional validation on a mouse embryo dataset [25] acquired with Stereo-seq technology, complementing the earlier experiments conducted on the 10x Visium platform (Fig 5A). Our results demonstrate that DMGCN effectively captures the laminar structures within the mouse embryo dataset, with 0.37ARI, showing strong alignment with the annotated distributions (Fig 5B). In contrast, the outcomes from SCANPY [26] identified an extensive region enveloping the neural crest and brain, which exhibited limited correspondence with the reference annotations.

In comparison, the segmentation results produced by DMGCN demonstrated stronger agreement with anatomical annotations, both at tissue boundaries and within internal regions. Clusters corresponding to distinct anatomical structures, such as the heart, brain, and liver, were consistently associated with expected marker gene expression patterns (Fig 5C and 5D), thereby validating DMGCN’s localization accuracy.

Performance on MERFISH data and downstream analysis

We also implemented DMGCN on five adjacent slices of the mouse brain [27] acquired through MERFISH, which include 254 genes labeled across seven distinct structures. For this experiment, we focused on slices 1, 10, and 31 from sample 1, along with slice 112 from sample 3, to compare the performance of spatial domain identification.

We first examined the domains identified by various methods. Fig 6 illustrates that DMGCN surpasses most other approaches, closely aligning with the manually annotated results. It is important to note that the seven primary spatial domains appear in all four tissue sections, although their shapes and sizes differ across these sections. This pattern proves instrumental for benchmarking the continuity and smoothness of domains generated by different algorithms. Remarkably, DMGCN was able to identify more continuous regions, and when compared to other methods, it nearly achieves the highest performance.

Ablation Study

To explore how DMGCN operates, we performed a series of ablation studies on the DLPFC dataset. These experiments involved the systematic removal of key components, including the ZINB loss and the Co-convolution loss to evaluate their individual impacts on model performance.

• DMGCN-w/o-: DMGCN has no spatial ZINB loss;
• DMGCN-w/o-: DMGCN has no Co-convolution loss;

As illustrated in Fig 7, DMGCN consistently surpasses its variants. Notably, DMGCN-w/o- has the lowest median ARI at just 0.38, the lowest in our ablation tests, highlighting the importance of fitting the gene expression data to the ZINB distribution, which effectively addresses the sparsity of spatial transcriptome data. Conversely, while DMGCN-w/o- exhibits the highest variance among all variants, despite achieving a competitive median ARI. This instability implies that consistent and accurate spatial domain identification depends critically on a loss function capable of effectively fusing gene expression and spatial information.

The ablation study further confirms that the integration of multimodal views, spatial and feature, coupled with an attention mechanism, is indispensable to DMGCN’s design. By adaptively learning and combining these complementary information sources, the model achieves not only higher accuracy in domain identification but also significantly improved robustness and stability.

Dataset description

We validate the performance of our proposed DMGCN using SRT datasets from three different platforms, which can be downloaded from the public website.

1. [1]. Human dorsolateral pre-frontal cortex (DLPFC): The 10x Visium dataset comprises 12 slices, and all slices were manually annotated. The article presents the experimental results of slice 151507, with 4,221 spots and 33,538 genes, and 151674, with 3635 spots and 33,538 genes. The data links to https://www.10xgenomics.com/datasets/visium-hd-cytassist-gene-expression-libraries-of-human-crc.
2. [2]. Human breast cancer: The 10x Visium dataset comprises 3,798 spots and 36,601 genes, with 20 regions manually annotated. The data links to https://support.10xgenomics.com/spatial-gene-expression/datasets/1.1.0/V1_Breast_Cancer_Block_A_Section_1.
3. [3]. Mouse embryo: The Stereo-seq data comprises 19,527 spots and 27,106 genes, with 12 regions manually annotated for some significant spots. The data links to https://db.cngb.org/stomics/mosta/.
4. [4]. Mouse brain: The MERFISH data comprises six samples, with 254 genes and 5–7 regions manually annotated. The samples 1 and 3 were selected to test the performance. The data links to https://doi.brainimagelibrary.org/doi/10.35077/g.21.

Data preprocessing

During the data processing stage, we first normalize the raw data using the "normalize" function from the SCANPY toolkit. And it’s first to filter genes with 100 "min_cells", meaning retention of genes expressed in at least 100 cells, and especially, for the MERFISH-seq Mouse Brain dataset, the extra preprocessing is to filter the cells with 100 "min_genes", meaning retention of cells with at least 100 genes detected, and filter the genes with"min_counts"as 50, for only cells with a total UMI count of ≥1000 were retained.

And then, the highly variable genes will be selected for that model can focus on significant gene information and train effectively. Generally, the 3000 highly variable genes are always a choice for different domain identification methods, and except for the MERFISH-seq Mouse Brain data, it has 254 genes in total, and we set the number of highly variable genes as 30 for training.

Before training, we constructed the spatial proximity map as well as the feature correlation map needed for the multi-view neural network encoder. We used the KNN to obtain the spatial proximity map by calculating the K nearest neighbors of each spatial site based on the Euclidean distance, as well as the feature correlation map by calculating the K nearest neighbors of each gene based on the cosine distance.Finally, we obtained the features graph , and spatial graph .And the method of DMGCN to select the K value, and made some adjustments to fit our model, such as the K value of the spatial adjacent graph on the 10X Visium dataset was set as 25, and on the Stereo-seq dataset, we set the K value as 50. Other parameters that are not mentioned here can be found in the S3 Table.

Data augmentation

For constructing contrastive learning information of the features, the corrupted features is made by randomly shuffling the gene expression and will be trained by the co-convolution with the spatial graph and features graph. The still preserves the topological structure of the original neighborhood graph. Specifically, we generated the contrastive learning label for the calculation of the positive and negative loss after the discriminator.

Multi-view graph convolutional encoder

The initial integration of genomic and spatial information is achieved by transforming the spatial adjacency graph via convolution in the spectral domain, thereby aggregating features from adjacent nodes. Follow Spatial-MGCN [17], this enriched representation is subsequently passed to a multi-layer graph convolutional architecture, which propagates it according to the following hierarchical rule:

where is the weight matrix of layer l in the spatial convolution, initially , and ReLU is the activation function. Specifically, we have and is the degree matrix of .

Applying the same approach used for spatial convolution, we next process feature information by performing feature convolution on the feature association matrix and gene expression data:

where is the weight matrix of the lth layer in the feature convolution, initially , and is the degree matrix of .

There is an inherent interconnection between gene expression and spatial localization. When addressing this relationship, we need to account for both their unique attributes and their shared latent features. We introduce a parameter-sharing architecture to train the fused embedding of gene expression and coordinates as Spatial-MGCN [17]. The propagation rule is formulated as follows:

where is the weight matrix of the lth layer in the co-convolution, and and are the lth layer embeddings extracted from the spatial and feature maps, respectively, by co-convolution, with the initial ==.

For corrupted features , the same co-convolution method is used to extract contrastive information from to increase the model’s ability to separate positive and negative sample points:

where is the weight matrix of the lth layer in the feature convolution, initially , , and, is the degree matrix of and respectively.

Attention mechanism

In tasks involving spatial domain identification, the relative contributions of gene expression, spatial context, and their joint representation to the final outcome may vary [17]. To dynamically capture the importance of each latent embedding, we incorporate an attention mechanism that adaptively weights these representations. Specifically, for the joint embedding, we first apply a nonlinear transformation, followed by a Softmax activation to normalize the attention scores across features. The resulting coefficients are computed as follows:

where the transformation employs a tanh activation function , W denotes a trainable weight matrix and represents a bias vector. This attention operation is applied to the spatial embedding, feature embedding, and their combined representation to produce an integrated latent embedding:

Discriminator

To make the feature embedding of spatial and feature co-convolution more stable and better differentiate between positive and negative sample points, we construct a contrastive learning mechanism between the feature embedding of co-convolution and the feature embedding of corrupted features by using a discriminator to redistribute the weights according to the feature and spatial maps and to generate the contrastive learning labels of the spatial and feature.

Firstly, we processed the co-convolution embedding with the graphs by the average readout for aggregating information of neighboring nodes:

For the corrupted features, we readout it in the same way:

We define the discriminator to identify the positive and negative spot:

For the and , we discriminate them with respectively and use the attention mechanism to generate the new contrastive learning label:

where and mean the contrastive learning labels of corrupted features in spatial graph and features graph.

FCN decoder

The FCN decoder is used to generate the embedding from the result of the attention mechanism and reconstruct the gene expression matrix for the loss calculation and domain identification visualization:

where means embedding that used to visualization though K-means and reconstructed gene expression matrix that not only can fit the distribution of features to ZINB distribution.

Co-Convolution loss

To improve the capability that delineates the positive and the negative points, the discriminator (*) will be applied to , ,and , to better constrain the co-convolution that combines with spatial and feature information. The loss of positive and negative points will be calculated by consistency loss for spatial and features co-convolution:

and use the binary cross-entropy with logarithmic for the contrastive learning label to calculate the positive loss and negative loss:

where mean the positive and negative samples of the spot i, and means the contrastive learning label of spot i.

Finally, consistency loss and the summary of positive loss and negative loss will be multiplied with temperature parameter and , respectively:

ZINB loss

Based on the observed gene expression data in spatial transcriptomics, we model its distribution using a ZINB formulation as Spatial-MGCN [17]. Let and represent the output parameter matrices of the decoder corresponding to the zero-inflation probability, mean expression, dispersion and biased vectors respectively. The ZINB probability distribution is formally defined as:

The reconstruction loss for the original gene expression is defined as the negative log-likelihood under the ZINB model, formulated as:

Loss function of DMGCN

The training process jointly optimizes the multi-view GCN encoder, the decoder, and the attention mechanism by assigning the proportion of the loss function with certain weights. The final training objective of the model in this paper is defined as:

where α is a temperature parameter that weighs the effects of the ZINB loss.

Benchmark methods

In our research, we assessed the effectiveness of DMGCN by contrasting it with five prominent methods. These include DeepST [12] and SpaGCN [13], which incorporate gene expression, spatial information, and H&E images; CCST [14], GraphST [15], and STAGATE [16], which integrate gene expression data along with spatial location details; and SCANPY, which focuses exclusively on gene expression profiles. Each approach leverages data with different complexities, making them suitable for the task of domain identification using SRT data.

Code availability

The code, scripts, and parameter configuration of DMGCN have been shared through GitHub (https://github.com/Jinsl-lab/DMGCN).

S1 Text. Platform and model settings.

https://doi.org/10.1371/journal.pcbi.1013369.s001

(DOC)

S1 Table. Summary of the datasets used in this study.

https://doi.org/10.1371/journal.pcbi.1013369.s002

(XLSX)

S2 Table. Memory Usage (MB).

https://doi.org/10.1371/journal.pcbi.1013369.s003

(XLSX)

S3 Table. Parameter configuration of all datasets.

https://doi.org/10.1371/journal.pcbi.1013369.s004

(XLSX)

S4 Table. Performance of various loss terms on HBC.

https://doi.org/10.1371/journal.pcbi.1013369.s005

(XLSX)

S5 Table. Performance of various methods on MERFISH.

https://doi.org/10.1371/journal.pcbi.1013369.s006

(XLSX)

S6 Table. Performance of various methods on DLPFC.

https://doi.org/10.1371/journal.pcbi.1013369.s007

(XLSX)

S7 Table. Performance of various methods on HBC.

https://doi.org/10.1371/journal.pcbi.1013369.s008

(XLSX)

S1 Fig. The spatial domain identification performance of four different parameter combinations on the HBC dataset is shown, among which the combination of term=(1, 1, 10) achieves the highest accuracy.

https://doi.org/10.1371/journal.pcbi.1013369.s009

(TIF)

S2 Fig. Domain identification results of 151507, 151508, and 151510 on DLPFC.

https://doi.org/10.1371/journal.pcbi.1013369.s010

(TIF)

S3 Fig. Domain identification results of 151669, 151670, and 151672 on DLPFC.

https://doi.org/10.1371/journal.pcbi.1013369.s011

(TIF)

S4 Fig. Domain identification results of 151673, 151674, and 151676 on DLPFC.

https://doi.org/10.1371/journal.pcbi.1013369.s012

(TIF)

S5 Fig. Domain identification results of various methods on MOSTA.

https://doi.org/10.1371/journal.pcbi.1013369.s013

(TIF)