A multi-modal graph fusion driven spatial transcriptomic region identification method

Abstract

The application discloses a multi-modal graph fusion driven spatial transcriptome region identification method, and belongs to the field of biomedical informatics and computational biology. The method comprises the following steps: step 1, data acquisition; step 2, data processing; step 3, multi-modal adjacency matrix generation; step 4, differentiable multi-modal graph structure fusion; step 5, variational graph autoencoder training; and step 6, clustering. By constructing a spatial domain identification framework based on multi-modal graph fusion and variational autoencoder, the application significantly improves the accuracy, continuity and biological interpretability of spatial transcriptome clustering analysis, effectively solves the deficiencies of existing methods in multi-modal fusion, graph structure optimization and feature discriminability, and provides a more reliable computing tool for in-depth revealing of tissue spatial heterogeneity in the fields of brain science, tumor microenvironment analysis and other frontier biomedical researches.

Description

TECHNICAL FIELD

[0001] The present application belongs to the field of biomedical informatics and computational biology, and relates to a spatial transcriptome data analysis method, in particular to a spatial domain recognition method based on multi-modal graph fusion and dynamic structure optimization. BACKGROUND

[0002] The rapid development of spatial transcriptomics technology provides strong technical support for analyzing the microenvironment of tissues in situ. Through advanced technology platforms such as 10x Visium, MerFISH and seqFISH, researchers can obtain the whole genome expression profile while preserving the spatial location information, greatly advancing the depth of understanding of the structural heterogeneity of tissues. Accurate identification of the "spatial domain" with specific functions in tissues has become a key link to understanding the physiological functions and pathological states of tissues. However, spatial transcript data has high dimensionality, sparsity, high noise and multi-modal characteristics, making precise spatial domain recognition face serious technical challenges. Existing spatial domain recognition methods have obvious evolution process and development bottlenecks when dealing with complex spatial transcript data.

[0003] Early methods mainly borrowed clustering strategies from single-cell RNA sequencing, directly using algorithms such as K-means, Louvain and Leiden, and only relying on gene expression similarity for region division. This method completely ignores spatial location information, resulting in fragmented clustering results in complex tissues, which cannot reflect the true biological tissue structure.

[0004] To make up for the lack of spatial information, researchers began to try to introduce spatial constraints into the clustering process. BayesSpace builds a hierarchical Bayesian model and introduces a spatial smoothing prior, forcing adjacent points to tend to belong to the same spatial domain, which improves the spatial coherence of the clustering results to some extent. However, this method only realizes the shallow fusion of spatial information and fails to establish a deep connection between gene expression and spatial location, and is prone to over-smoothing in the boundary region of the tissue.

[0005] With the development of graph neural network technology, researchers began to explore deeper information fusion methods. By regarding each measurement point as a node in the graph and constructing adjacency relationships based on spatial coordinates, a leap from single-modal to dual-modal fusion was achieved. SpaGCN made an important attempt in this direction, integrating gene expression and spatial coordinates in the graph convolution network framework for joint learning. Subsequent methods such as STAGATE further introduced graph attention mechanisms, which can adaptively learn the importance between adjacent nodes, enhancing the modeling ability of spatial context.

[0006] However, the above methods still have obvious shortcomings in information integrity - they only integrate two modalities of spatial location and gene expression, but ignore the histological images that contain rich histological structure information. To overcome this limitation, researchers have begun to explore solutions for three-modal fusion. The EfNST method represents an important development in this direction, by using a pre-trained EfficientNet-B0 network to extract deep features of tissue sections, and combining graph attention mechanisms to achieve cross-modal feature interaction. Although this method has made progress in the utilization of morphological information, its feature extraction and graph learning processes are separate, and the pre-trained image features and gene expression features are separated in the semantic space, making it impossible to achieve true end-to-end collaborative optimization.

[0007] The common technical difficulties faced by current multi-modal fusion methods mainly manifest in three aspects: First, different modal data have heterogeneous characteristics, with spatial coordinates being continuous numerical values, gene expression being high-dimensional sparse matrices, and tissue images being two-dimensional pixel arrays. Traditional splicing or weighted fusion is difficult to fully exploit their deep connections. Second, existing methods mostly use static graph structures, which cannot dynamically adjust the fusion weights between modalities during training. Finally, the separation of feature extraction and fusion learning prevents deep interaction and complementary enhancement of modal information in a unified framework.

[0008] Existing spatial transcriptomic region identification methods have obvious limitations in fully utilizing multi-modal information and optimizing graph structure modeling. Specifically:

[0009] (1) Insufficient multi-modal information fusion: Most existing methods can only cooperatively utilize spatial location and gene expression modalities, and fail to effectively integrate tissue section images that contain rich morphological information, or only use simple feature splicing and other shallow fusion strategies, which cannot achieve deep interaction and adaptive collaboration of cross-modal information, limiting the model's comprehensive analysis ability of complex biological systems.

[0010] (2) Graph structure modeling is rigid: Graph neural network-based methods heavily rely on pre-defined, static adjacency relationships (such as fixed graph construction based on spatial K-nearest neighbors). This fixed graph structure cannot dynamically evolve and optimize based on learned features during model training, limiting the model's ability to capture local heterogeneity and complex interaction patterns in tissues.

[0011] (3) Weak discriminability of feature representation: Due to the limitations of insufficient fusion and rigid graph structure, the latent embedding learned by existing methods is not discriminative enough, making it difficult to clearly distinguish between functionally distinct but spatially adjacent fine tissue structures (such as different functional subregions of tumor margins), resulting in blurred spatial domain boundaries and reduced resolution.

[0012] (4) Sensitivity to noise and data sparsity: In high-dimensional, high-noise, and large zero-value spatial transcriptome data, the robustness of existing methods is insufficient, and they are easily disturbed by technical noise, leading to fragmentation of clustering results or regional identification errors. SUMMARY

[0013] In order to overcome the limitations of existing spatial transcriptome analysis techniques in multi-modal information fusion and static graph structure, the present application provides a multi-modal graph fusion driven spatial transcriptome region identification method. The method proposes a spatial domain identification framework based on multi-modal graph fusion and variational autoencoder joint optimization, which completely breaks through the limitations of traditional single or double modal analysis by constructing a three-modal unified graph structure of spatial position, gene expression and tissue morphology image. In particular, a differentiable edge-by-edge fusion mechanism is designed to dynamically learn the optimal weight ratio between spatial proximity, gene expression correlation and morphological similarity through a multi-layer perceptron (MLP), breaking the limitations of traditional fixed weight or manually set fusion coefficients. At the model architecture level, the present application innovatively integrates variational graph autoencoder and multi-modal graph fusion module in depth, realizing end-to-end joint training. The framework updates the graph structure parameters and node embedding representation simultaneously through a unified optimization process, so that the adjacency matrix changes from static preset to dynamic evolution, significantly improving the adaptability of the model to complex tissue structures. In particular, the model optimizes the reconstruction loss and KL divergence cooperatively, which not only guarantees the spatial separability of the latent embedding, but also maintains the biological consistency of the feature expression.

[0014] The purpose of the present application is realized by the following technical solutions:

[0015] A multi-modal graph fusion driven spatial transcriptome region identification method, comprising the following steps:

[0016] Step 1, data acquisition:

[0017] Download the spatial transcriptome data from the public database;

[0018] Step 2: data processing:

[0019] After data reading, perform preprocessing operations including gene expression data normalization and gene screening, the specific steps are as follows:

[0020] Step 2-1, gene expression data normalization: log-transform the original gene expression matrix, and use the SCANPY package for data normalization;

[0021] Step 2-2, gene screening: z-score standardization is performed on the normalized expression matrix, and the top highly variable genes are selected as input features to form the preprocessed gene expression matrix , For the number of genes, Number of spots;

[0022] Step 3: Generation of multimodal adjacency matrix:

[0023] Spatial adjacency matrices were constructed from three modalities: spatial location, gene expression, and histological image. Gene expression adjacency matrix Image adjacency matrix Spatial feature matrix Gene expression feature matrix and image feature matrix The specific steps are as follows:

[0024] Step 3-1, Spatial Adjacency Matrix Construction: Using the sine and cosine position coding method, the two-dimensional spatial coordinates are mapped into high-dimensional features to obtain the spatial feature matrix. and based on The spatial adjacency matrix is ​​constructed using the K-nearest neighbor algorithm. ;

[0025] Step 3-2, Gene Expression Adjacency Matrix Construction: First, the preprocessed expression matrix... Spatial smoothing enhancement is performed to obtain the spatially enhanced representation matrix. Subsequently Principal component analysis was performed to reduce the dimensionality and obtain the gene expression feature matrix. Based on the dimensionality-reduced features, the K-nearest neighbor algorithm is used to construct the gene expression adjacency matrix. ;

[0026] Step 3-3: Image Adjacency Matrix Construction: A pre-trained convolutional neural network is used to extract deep features of the organized image patches centered on each point. Principal component analysis is then used to reduce the dimensionality, resulting in a low-dimensional image feature matrix. ,based on The Euclidean distance is used to construct the image adjacency matrix using the radius nearest neighbor method. ;

[0027] Step 4: Differentiable multimodal graph structure fusion:

[0028] Step 4-1: Convert the spatial feature matrix Gene expression feature matrix Image feature matrix By concatenating the features along the feature dimension, a comprehensive feature representation matrix is ​​obtained. ;

[0029] Step 4-2: The adjacency matrices and their features of spatial, gene expression, and image modalities are fused using an edge-by-edge multilayer perceptron (MLP) to generate a multimodal adjacency matrix. ;

[0030] Step 5, variational graph autoencoder training:

[0031] Step 5-1, data preparation: the multi-modal adjacency matrix and the integrated feature expression matrix are input into the variational graph autoencoder (VAE) for training;

[0032] Step 5-2, encoder training: in the encoder part, first, the node features are processed through the graph convolution layer, and the multi-modal adjacency matrix is used to propagate the neighborhood features to obtain the latent features of each node; then, the encoder uses a Gaussian distribution to approximate the posterior distribution of the latent space, and calculates the mean vector and the standard deviation vector of each node;

[0033] Step 5-3, decoder training: the decoder receives the latent variable output by the encoder and maps it back to the input feature space; the decoder gradually expands the dimension of the latent embedding through multiple fully connected layers, and enhances the non-linear expression ability through the LeakyReLU activation function; finally, the decoder restores the latent features to the original input features;

[0034] Step 5-4, loss function calculation: to optimize the latent space and feature reconstruction, two parts of the loss are considered during the training process: reconstruction loss and KL divergence loss;

[0035] Step 5-5, optimization process: by minimizing the weighted sum of the reconstruction loss and the KL divergence loss, the parameters of the encoder and the decoder are optimized, and the Adam optimizer is used to update the model parameters jointly. In each round of training, the model performs forward propagation, calculates the loss, and updates the parameters through back propagation;

[0036] Step 6, clustering:

[0037] The latent variable is clustered using the K-Means algorithm to complete the spatial domain clustering division of the spatial transcriptome.

[0038] Compared with the prior art, the present application has the following advantages:

[0039] (1) Multi-modal information graph fusion: the present application integrates spatial coordinates, gene expression data and H&E stained tissue images to construct a unified framework of multi-modal information sources. On this basis, the adjacency relationship is effectively mapped to the same graph structure, so that the information of different modalities can be optimized in the same space.

[0040] (2) Dynamic adaptive graph fusion mechanism: The present application realizes a differentiable graph fusion mechanism through edge-by-edge multilayer perceptron, so that the graph structure can be adaptively adjusted during the training process according to the learned features, thereby optimizing the adjacency relationship and significantly improving the adaptability of the graph model to complex organizational structures.

[0041] (3) Deep multi-modal feature collaboration: The present application generates comprehensive node feature representation by standardizing and fusing spatial, gene expression and morphological data at the feature level through deep multi-modal feature collaboration. This method not only enhances the complementarity of each modality information, but also enables the model to achieve deep collaboration at the feature and structure levels, improving the discriminability of potential features and thus improving the recognition accuracy of the spatial domain.

[0042] (4) Variational graph autoencoder joint optimization to handle noise: The present application uses variational graph autoencoder joint optimization to ensure the separability and regularity of the latent space by simultaneously optimizing reconstruction loss and KL divergence. In this process, the model can effectively suppress the influence of noise and maintain high robustness and clustering quality under high noise and sparse data conditions.

[0043] Through the above technical solutions, the present application effectively solves the shortcomings of existing methods in multi-modal fusion, graph structure optimization and feature discriminability. These improvements significantly improve the accuracy, continuity and biological consistency of spatial domain recognition, providing a more reliable computational tool for in-depth revealing of tissue spatial heterogeneity in the fields of brain science, tumor microenvironment analysis and other frontier biomedical research. BRIEF DESCRIPTION OF DRAWINGS

[0044] Figure 1 is the overall framework diagram;

[0045] Figure 2 is data reading and multi-modal information processing;

[0046] Figure 3 is multi-modal feature and graph structure fusion;

[0047] Figure 4 is variational autoencoder training and spatial domain recognition;

[0048] Figure 5 is the visualization result of spatial domain recognition, (a) manually annotated cortical levels (L1-L6) and white matter (WM) regions, (b) clustering results of the present application method, (c) clustering results of other methods;

[0049] Figure 6 is clustering performance comparison (ARI score);

[0050] Figure 7 is UMAP and PAGA visualization;

[0051] Figure 8 The spatial domain visualization comparison chart of the method of the present application and eight comparison methods on Human_breast_cancer;

[0052] Figure 9 For spatial domain analysis on mouse forebrain data set, (a) pathologist annotated region result, (b) spatial clustering result of the method of the present application, (c) identification result of domain 17 and domain 36, (d) region identification result of other eight methods. DETAILED DESCRIPTION

[0053] The technical solutions of the present application will be further described below in conjunction with the drawings, but are not limited thereto, and any modification or equivalent replacement to the technical solutions of the present application without departing from the spirit and scope of the technical solutions of the present application shall be covered in the protection scope of the present application.

[0054] The present application provides a multi-modal graph fusion driven spatial transcriptome region identification method, as shown in Figure 1 and Figure 2 The method comprises the following steps:

[0055] Step 1, data acquisition:

[0056] Spatial transcriptome data is downloaded from a public database.

[0057] Step 2: data processing:

[0058] After data reading, pre-processing operation is performed:

[0059] Step 2-1, gene expression data normalization: log transformation is performed on the original gene expression matrix, and data normalization is performed using the SCANPY package.

[0060] Step 2-2, gene screening: z-score standardization is performed on the normalized expression matrix, and the top highly variable genes are selected as input features to form the pre-processed gene expression matrix ( is the number of genes, is the number of spots). For data sets with fewer highly variable genes (such as MerFish data set), all available genes are selected as input.

[0061] Step 3, multi-modal adjacency matrix generation:

[0062] Adjacency matrices and feature matrices are constructed from three modalities of spatial position, gene expression and histological image, respectively:

[0063] Step 3-1, Spatial Adjacency Matrix Construction: Utilize the sine-cosine position encoding method to map two-dimensional spatial coordinates into high-dimensional features, obtaining a spatial feature matrix , and based on this feature, use the K-Nearest Neighbor algorithm to construct the spatial adjacency matrix .

[0064] Step 3-2, Gene Expression Adjacency Matrix Construction: First, perform spatial smoothing enhancement on the preprocessed expression matrix to obtain the spatial enhanced expression matrix , then perform principal component analysis dimensionality reduction on the enhanced expression matrix to obtain the gene expression feature matrix , and based on the reduced features, use the K-Nearest Neighbor algorithm to construct the gene expression adjacency matrix .

[0065] Step 3-3, Image Adjacency Matrix Construction: Use a pre-trained convolutional neural network to extract the deep features of the tissue image block centered at each point, and perform dimensionality reduction through principal component analysis to obtain the low-dimensional image feature matrix , based on the Euclidean distance of this feature, use the radius neighbor method to construct the image adjacency matrix .

[0066] Step 4, Differentiable Multimodal Graph Structure Fusion:

[0067] Step 4-1, concatenate the spatial feature matrix , gene expression feature matrix , and image feature matrix in the feature dimension to obtain the comprehensive feature expression matrix .

[0068] Step 4-2, fuse the adjacency matrices of space, gene expression, and image modalities and their features through edge-by-edge multilayer perceptron (MLP) to generate a multimodal adjacency matrix . In this process, by weighting learning the edge features of different modalities, the relationships between different modalities are optimized, and a comprehensive graph structure is generated.

[0069] Step 5, Variational Graph Autoencoder Training:

[0070] Step 5-1, Data Preparation: Input the multimodal adjacency matrix and the comprehensive feature expression matrix into the Variational Graph Autoencoder (VAE) for training. The encoder part receives these inputs and propagates neighborhood information through the Graph Convolutional Network (GCN) to aggregate first-order neighborhood features.

[0071] Step 5-2, Encoder Training: In the encoder part, the node features are first processed through graph convolutional layers, and a multimodal adjacency matrix is ​​used. The neighborhood features are propagated to obtain the latent features of each node. Next, the encoder approximates the posterior distribution of the latent space using a Gaussian distribution and calculates the mean vector for each node. and standard deviation vector These parameters are learned by the neural network through training. A reparameterization technique is employed to make the latent embeddings differentiable, facilitating optimization via gradient backpropagation.

[0072] Step 5-3, Decoder Training: The decoder receives the latent variables output by the encoder. The latent embeddings are then mapped back to the input feature space. The decoder progressively expands the dimension of the latent embeddings through multiple fully connected layers and enhances the nonlinear expressive power through the LeakyReLU activation function. Ultimately, the decoder restores the latent features to the original input features, thereby minimizing the reconstruction error.

[0073] Step 5-4, Loss Function Calculation: To optimize the latent space and feature reconstruction, two parts of the loss are considered simultaneously during training: reconstruction loss and KL divergence loss. The KL divergence term constrains the distribution of the latent space, making the latent features approximate a standard normal distribution, thereby ensuring the separability and consistency of the latent space.

[0074] Step 5-5, Optimization Process: Optimize the encoder and decoder parameters by minimizing the weighted sum of the reconstruction loss and the KL divergence loss. The Adam optimizer is used to jointly update the model parameters. In each training round, the model performs forward propagation, calculates the loss, and updates the parameters via backpropagation.

[0075] Step 6, Clustering:

[0076] For latent variables The K-Means algorithm was used for clustering to complete the spatial domain clustering of the spatial transcriptome.

[0077] This invention significantly improves the accuracy, continuity, and biological interpretability of spatial transcriptome clustering analysis by constructing a spatial domain recognition framework based on multimodal graph fusion and variational autoencoders. Specifically:

[0078] (1) Dynamic graph structure optimization: Through a differentiable edge-by-edge MLP fusion mechanism, the optimal weight ratio between spatial proximity, gene expression similarity and morphological features is adaptively learned to generate a graph structure that continuously evolves during training, breaking through the limitations of traditional static graph modeling of complex organizational structures.

[0079] (2) Deep multimodal feature coordination: The spatial coordinates, gene expression, and tissue image three modal features are spliced and jointly learned in a unified framework, realizing the deep interaction and complementary enhancement of cross-modal information, and significantly improving the discriminability and biological consistency of the latent embedding.

[0080] (3) Variational autoencoder enhanced robustness: The variational graph autoencoder is used to train the fused graph structure end-to-end, and through the joint optimization of reconstruction loss and KL divergence, a regular and separable low-dimensional representation is learned while ensuring data fidelity, effectively improving the adaptability of the model to high noise and high sparse data.

[0081] Embodiments:

[0082] Step 1: Data acquisition.

[0083] Five spatial transcriptome datasets covering DLPFC (human dorsolateral prefrontal cortex), human breast cancer, and mouse brain samples were obtained from public databases (such as 10X Genomics, MerFish), and all datasets were accompanied by benchmark clustering labels. The dataset contains the following information: spatial position coordinates (two-dimensional coordinates of each spot), gene expression matrix (rows represent genes and columns represent spots), and H&E staining tissue section images.

[0084] Step 2: Data processing.

[0085] The original gene expression matrix was logarithmically transformed and normalized using the SCANPY tool; the normalized expression matrix was z-score standardized. The top 2000 high-variability genes were selected as input features to form the preprocessed expression matrix ( is the number of genes, is the number of spots).

[0086] Step 3: Multimodal adjacency matrix construction.

[0087] Step 3-1, spatial adjacency matrix construction: The two-dimensional spatial coordinates are mapped to high-dimensional features using the sine-cosine position encoding method, obtaining a spatial feature matrix with dimensions . Specifically, for each dimension, a one-dimensional frequency band is used, and the encoding of the point is:

[0088]

[0089] (1)

[0090] where represents the point the spatial coordinates of the space, , is the dimension of the spatial feature matrix . The band selection ensures that each coordinate dimension is effectively encoded at different frequencies.

[0091] Subsequently, according to the encoding space of the position, K- nearest neighbors are used for composition, generating a spatial adjacency matrix :

[0092] (2)

[0093] wherein, is an element of the spatial adjacency matrix, representing the relationship between point and point in the spatial coordinate mode, that is, and whether there is a spatial connection between them; represents the nearest neighbor set of point in the spatial position; is one of the spatial position points, representing the spatial relationship with other points .

[0094] Step 3-2, gene expression adjacency matrix construction: based on the preprocessed expression matrix , identify the Euclidean neighbors of each spot in the coordinate space. Then, calculate the average gene expression in the neighborhood of each point, and linearly superimpose it by weight α to generate a spatial enhanced expression matrix :

[0095] (3)

[0096] wherein, is the generated enhanced feature representation, enhanced by the features of the weighted neighborhood points, is the feature expression of at the gene expression level; is the neighborhood set of ; is the number of spots, is the index of a spot position.

[0097] Next, perform principal component analysis (PCA) on the rows and columns of respectively, and take the first principal components to obtain the gene expression feature matrix and . Then, perform dimensionality reduction on Knn is performed to obtain the gene expression adjacency matrix :

[0098] (4)

[0099] wherein, represents the relationship between the point and the point in the gene expression mode, represents the nearest neighbor set of the point in the gene expression data, is one of the position points in the space, and represents the spatial relationship with other points .

[0100] Step 3-3, image adjacency matrix construction: using the morphological information of the tissue section to depict the relationship between the space points. First, according to the coordinates of each spot in the spatial transcriptome data, a local image block centered on each spot is obtained. Then, by PCA dimension reduction, the image feature matrix :

[0101] (5)

[0102] wherein, is the number of spots; is the feature expression of the image layer; is the image feature dimension corresponding to each spot; is the image feature dimension corresponding to each spot; is input into the downstream model as the node feature of the image mode, and further obtains the image adjacency matrix :

[0103] (6)

[0104] wherein, represents the relationship between the point and the point in the image mode, i.e., whether the point and the point are adjacent in the image space, represents the image feature vector of , represents the image feature vector of , represents the Euclidean distance between and in the image feature space, is a radius threshold.

[0105] Step 4: differentiable multi-modal graph structure fusion. ​

[0106] Step 4-1, Adjacency matrix stacking: Stack the adjacency matrices in the edge dimension under three modalities to construct a fusion edge feature vector containing spatial modalities, gene expression modalities, and image modalities , which represents the similarity of nodes and nodes under different modalities. The specific representation is as follows:

[0107] (7)

[0108] wherein, represents the relationship between point and point in the gene expression modality, represents the relationship between point and point in the spatial coordinate modality, represents the relationship between point and point in the image modality.

[0109] Step 4-2, MLP fusion: is input into a multi-layer perceptron (MLP) to learn the fusion weight of cross-modal edges. The MLP contains two layers of fully connected mapping, each followed by a LeakyReLU activation function (negative half-axis slope ), and after the first layer, a Dropout regularization is added to alleviate overfitting. Finally, the output is compressed to the [0, 1] interval by the Sigmoid function, and the fusion edge weight is obtained:

[0110] (8)

[0111] (9)

[0112] (10)

[0113] (11)

[0114] wherein, , represent the intermediate representation processed by the LeakyReLU activation function, is the parameter set to be learned, represents the Sigmoid function, is an activation function with a negative slope of 0.2, is the bias term corresponding to the weight matrix, and the final output multi-modal adjacency matrix is The process is completely trained to measure the mode, so it can be optimized end-to-end by gradient back propagation during the training process.

[0115] Step 4-3, expression feature fusion: in the feature level fusion, the feature matrix of three modalities is respectively normalized (zero mean, unit variance) to eliminate the difference in scale of different modalities. Then, the features are spliced according to the feature dimension to obtain the comprehensive feature expression matrix

[0116] (12)

[0117] Step 5: training of variational graph autoencoder.

[0118] Step 5-1, encoder part training: the goal of the encoder part is to convert the input comprehensive feature expression matrix and the multi-modal adjacency matrix into a distribution in the latent space. The encoder propagates the neighborhood information of the nodes through the graph convolution layer (GCN), and approximates the posterior distribution of the latent space through the Gaussian distribution. The specific process is as follows:

[0119] (1) graph convolution layer: in the graph convolution layer of the encoder, each layer of graph convolution operation propagates the features of the nodes through the adjacency matrix

[0120] (13)

[0121] wherein is the multi-modal adjacency matrix; L=1,2, is the learnable weight matrix, is the feature dimension of the previous layer, is the feature dimension of the current layer, is the activation function with a negative slope of 0.2, is the bias term of the Lth layer, used to adjust the output.

[0122] (2) Gaussian distribution approximation: the encoder approximates the posterior distribution of the latent space through the Gaussian distribution. The latent feature representation of each spot is represented by the mean vector and the standard deviation vector, and these two parameters are learned by the encoder through the neural network.

[0123] (14)

[0124] wherein represents a Gaussian distribution;​​​​​ and It is calculated through the encoder's neural network, representing the nodes. The mean and standard deviation of the latent characteristics; Representing latent variables, at this point, It is the node feature matrix after graph convolution processing.

[0125] (3) Reparameterization Technique: To make the sampling process of the latent variables differentiable and to enable optimization through gradient backpropagation, a reparameterization technique is employed. This involves using the mean... and standard deviation and noise variables sampled from the standard normal distribution. Latent variables can be obtained. Sampling process:

[0126] (15)

[0127] in, Representing latent variables The Each sample value, Let the noise variable represent a standard normally distributed noise variable. This represents a standard normal distribution.

[0128] Step 5-2, Encoder Decoding: Decoding the latent variables in the latent space The input data is mapped back to the original data space and reconstructed. The decoder maps the latent space to the data space through layer-by-layer linear transformations and nonlinear activation functions.

[0129] Step 5-3, Loss Function: The training objective is to minimize the weighted sum of the reconstruction loss and the KL divergence loss. During training, the following loss function is minimized:

[0130] To measure the difference between the decoded features and the original input, we define the reconstruction loss:

[0131] (16)

[0132] in, The expression matrix reconstructed by the decoder, This represents the Frobenius norm, used to measure the element-wise reconstruction error of a matrix.

[0133] However, relying solely on reconstruction error is insufficient to constrain the distribution of the latent space. To avoid disordered distribution or overfitting in the latent space, a KL divergence term is further introduced to measure the approximate posterior. With standard Gaussian prior Differences:

[0134] (17)

[0135] where, is the embedding vector is the dimension of is the mean of the th latent variable, is the variance of the th latent variable. The KL divergence term can constrain the latent embedding to approximate the standard normal distribution, thus ensuring its distribution regularity and separability, enhancing the separability and consistency between latent variables.

[0136] Finally, the overall optimization objective function combines the reconstruction loss and the KL divergence:

[0137] (18)

[0138] where, is the weight parameter of the KL divergence. During training, the Adam optimizer is used to update the parameter set jointly. In each iteration, the model performs forward propagation, calculates the loss function and updates the parameters by backpropagation.

[0139] Step 6: Latent variable generation and clustering.

[0140] After training, in order to perform clustering analysis, the latent embedding is clustered using the K-Means algorithm to complete the spatial domain clustering division of the spatial transcriptome.

[0141] As can be seen from Figure 5 , the method of the present application clearly distinguishes each cortical level and white matter region, and is superior to other methods in level distinction.

[0142] As can be seen from Figure 6 , the method of the present application performs best and has strong stability. The closer the homogeneity and normalized mutual information (NMI) scores are to 1, the better the effect. The results show that the method of the present application performs excellently under multiple indicators, and the method of the present application can reach up to 0.62 under the ARI indicator, which is superior to other comparative methods.

[0143] As can be seen from Figure 7 , the method of the present application presents a clear trajectory from L1 to L6 to white matter, while the level boundaries of other methods are relatively blurred.

[0144] As can be seen from Figure 8It can be known that the method can better identify the boundary of the spatial domain, and the overall exhibits more coherent and clear spatial division. In comparison, EfNST and DeepST are divided into multiple regions in the IDC_5 region, and SpaGCN, Seurat and STAGATE are somewhat disordered in the overall division.

[0145] By Figure 9 It can be known from the visual comparison between the methods that the method is still relatively clear in the overall block separation.

Claims

1. A method for spatial transcriptome region identification driven by multimodal graph fusion, characterized in that... The method includes the following steps: Step 1, Data Acquisition: Download spatial transcriptome data from publicly available databases; Step 2: Data Processing After the data is read, preprocessing operations are performed, including gene expression data normalization and gene screening; Step 3: Generation of multimodal adjacency matrix: Spatial adjacency matrices were constructed from three modalities: spatial location, gene expression, and histological image. Gene expression adjacency matrix Image adjacency matrix Spatial feature matrix Gene expression feature matrix and image feature matrix ; Step 4: Differentiable multimodal graph structure fusion: Step 4-1: Convert the spatial feature matrix Gene expression feature matrix Image feature matrix By concatenating the features along the feature dimension, a comprehensive feature representation matrix is ​​obtained. ; Step 4-2: The adjacency matrices and their features of spatial, gene expression, and image modalities are fused using an edge-by-edge multilayer perceptron to generate a multimodal adjacency matrix. The specific steps are as follows: Adjacency Matrix Stacking: The adjacency matrices of the three modalities are stacked along the edge dimension to construct a fused edge feature vector that includes spatial modality, gene expression modality, and image modality. : ； in, Point and points The relationship between them under gene expression modalities Point and points The relationship between them in spatial coordinate modes Point and points The relationship between them in the image modality; MLP Fusion: The input is fed into an edge-by-edge multilayer perceptron to learn the fusion weights for cross-modal edges. The output is compressed to the [0,1] interval by a sigmoid function to obtain the fusion edge weights. ； ； ； ； in, , This represents the intermediate representation after processing with the LeakyReLU activation function. } represents the set of parameters to be learned. This represents the Sigmoid function. Using an activation function with a negative slope of 0.2, the final output multimodal adjacency matrix is: ; Step 5: Variational graph autoencoder training: Step 5-1, Data Preparation: Prepare the multimodal adjacency matrix. and comprehensive feature expression matrix The input variational graph autoencoder is used for training. Step 5-2, Encoder Training: In the encoder part, the node features are first processed through graph convolutional layers, and a multimodal adjacency matrix is ​​used. The encoder propagates neighborhood features to obtain the latent features of each node; then, it approximates the posterior distribution of the latent space using a Gaussian distribution and calculates the mean vector of each node. and standard deviation vector ; Step 5-3, Decoder Training: The decoder receives the latent variables output by the encoder. The latent embeddings are mapped back to the input feature space; the decoder progressively expands the dimension of the latent embeddings through multiple fully connected layers and enhances the nonlinear expressive power through the LeakyReLU activation function; finally, the decoder restores the latent features to the original input features. Step 5-4, Loss function calculation: To optimize the latent space and feature reconstruction, two parts of loss are considered simultaneously during training: reconstruction loss and KL divergence loss; Step 5-5, Optimization process: Optimize the parameters of the encoder and decoder by minimizing the weighted sum of the reconstruction loss and KL divergence loss, and jointly update the model parameters using the Adam optimizer. In each round of training, the model performs forward propagation, calculates the loss, and updates the parameters through backpropagation. Step 6, Clustering: For latent variables The K-Means algorithm was used for clustering to complete the spatial domain clustering of the spatial transcriptome.

2. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 1, characterized in that... The specific steps of step 2 are as follows: Step 2-1, Gene Expression Data Normalization: Perform a logarithmic transformation on the original gene expression matrix and normalize the data using the SCANPY package; Step 2-2, Gene Screening: Perform z-score normalization on the normalized expression matrix and screen for the top gene expression. Using highly variable genes as input features, a preprocessed gene expression matrix is ​​constructed. .

3. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 1, characterized in that... The specific steps of step 3 are as follows: Step 3-1: Spatial Adjacency Matrix Construction: Using the sine and cosine position coding method, the two-dimensional spatial coordinates are mapped into high-dimensional features to obtain the spatial feature matrix. and based on The spatial adjacency matrix is ​​constructed using the K-nearest neighbor algorithm. ; Step 3-2, Gene Expression Adjacency Matrix Construction: First, the preprocessed expression matrix... Spatial smoothing enhancement is performed to obtain the spatially enhanced representation matrix. Subsequently Principal component analysis was performed to reduce the dimensionality and obtain the gene expression feature matrix. Based on the dimensionality-reduced features, the K-nearest neighbor algorithm is used to construct the gene expression adjacency matrix. ; Step 3-3: Image Adjacency Matrix Construction: A pre-trained convolutional neural network is used to extract deep features of organized image patches centered at each point. Principal component analysis is then used to reduce the dimensionality, resulting in a low-dimensional image feature matrix. ,based on The Euclidean distance is used to construct the image adjacency matrix using the radius nearest neighbor method. .

4. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 3, characterized in that... In step 3-1, for each dimension, a one-dimensional frequency band is used, then the point is encoded. for: , ； in, Represented as a point spatial coordinates, , Spatial characteristic matrix The dimension; Subsequently, based on the encoding space of that location, a graph is constructed using K-nearest neighbors to generate a spatial adjacency matrix. : ； in, Point and points The relationship between them in spatial coordinate modes; Indicates the midpoint of a spatial location of A set of nearest neighbors; It is one of the points in space, indicating its position relative to other points. Spatial relationships.

5. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 3, characterized in that... The specific steps of step 3-2 are as follows: Based on preprocessed expression matrix Identify each spot in coordinate space A neighbor of the Euler family; The mean gene expression in the neighborhood of each location is calculated and linearly superimposed using weight α to generate a spatially enhanced expression matrix. : ； in, The enhanced feature representation is augmented using features from weighted neighborhood points. At the gene expression level Feature representation; for of A set of spatial neighborhoods; For the number of spots, Index for a specific spot location; To each Perform principal component analysis on rows and columns, and take the first row and column. One principal component is used to obtain the gene expression feature matrix. and ; right Perform k-nn to obtain the gene expression adjacency matrix. : ； in, Point and points The relationship between them under gene expression modalities Indicates the midpoint of gene expression data of A nearest neighbor set It is one of the points in space, indicating its position relative to other points. Spatial relationships.

6. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 3, characterized in that... The specific steps of step 3-3 are as follows: Based on the coordinates of each spot in the spatial transcriptome data, a local image patch centered on each spot is obtained; The image feature matrix is ​​obtained by PCA dimensionality reduction. : ； in, Number of spots; For image level Feature representation; The node features of the image modality are input into the downstream model, and the image adjacency matrix is ​​further obtained. : ； in, Point and points The relationship between them in the image modality express Image feature vectors, express Image feature vectors, express and Euclidean distance in image feature space This is the radius threshold.

7. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 1, characterized in that... The specific steps of step 5-2 are as follows: (1) Graph convolutional layer: In the graph convolutional layer of the encoder, each graph convolutional operation propagates the features of the nodes through the adjacency matrix. : ； in, This is a multimodal adjacency matrix; L=1,2, The weight matrix is ​​a learnable matrix. For activation functions with a negative slope of 0.2, It is the bias term of the Lth layer, used to adjust the output; (2) Gaussian distribution approximation: The encoder approximates the posterior distribution of the latent space using a Gaussian distribution. : ； in, Represent a Gaussian distribution; and Represents a node The mean and standard deviation of the latent characteristics; Indicate latent variables; (3) Reparameterization technique: using the mean and standard deviation and noise variables sampled from the standard normal distribution. to obtain latent variables Sampling process: ； in, Representing latent variables The Each sample value, Let the noise variable represent a standard normally distributed noise variable. This represents a standard normal distribution.

8. The multimodal graph fusion-driven spatial transcriptome region identification method according to claim 1, characterized in that... In step 5-4, the reconstruction loss Defined as: ； in, The expression matrix reconstructed by the decoder, Denotes the Frobenius norm; KL divergence loss Defined as: ； in, Embedded vector Dimensions For the first The mean of the 10 latent variables, For the first The variance of each latent variable.

Images