Spatial Domain Identification Method Based on Feature Extraction of Spatial Transcriptomics Data

Abstract

The present invention discloses a spatial domain recognition method based on spatial transcriptomics data feature extraction, which mainly solves the problems of overfitting in the feature extraction of spatial transcriptomics data and low accuracy in spatial domain recognition in the prior art. The implementation scheme is as follows: preprocess the gene expression data and spatial information measured in spatial transcriptomics; construct a gene similarity network and a spatial neighborhood network based on the gene expression feature matrix and spatial information; perform data augmentation on the gene similarity network and the spatial neighborhood network; construct a feature extraction model, and input the augmented data into the model to calculate the contrastive loss and the reconstruction loss; train the model according to the calculated loss, input the non-augmented data into the trained model to obtain a low-dimensional embedding; perform clustering on the low-dimensional embedding to complete the spatial domain recognition. The present invention avoids overfitting in the feature extraction process, improves the accuracy of spatial domain recognition, and can be used to provide reference data for exploring biological development and disease treatment.

Description

Technical Field

[0001] The present invention belongs to the technical field of data mining, and particularly relates to a spatial domain recognition method, which can be used to provide reference data for exploring biological development and disease treatment. Background Art

[0002] In tissue sections, some regions have similar spatial gene expression profiles, forming specific structures or sub-structures in the tissue. These regions have different functional partitions due to the composition of cell types and gene expression, and thus form spatial domains with specific biological significance. The recognition of spatial domains is crucial for studying the influence of tissue structure and cell-cell interaction.

[0003] Single-cell transcriptome sequencing technology scRNA-seq can be used to provide high-resolution gene expression profiles. However, due to the inability to retain spatial location information during sample preparation, it brings limitations to downstream analysis. Spatial transcriptome sequencing technologies, including imaging techniques based on in situ hybridization and in situ sequencing technologies based on spatial barcodes, provide both gene expression profiles and spatial location information, and this spatial information is crucial for understanding the biological significance of healthy tissue development and the tumor microenvironment of diseases. Therefore, the introduction of spatial transcriptome data helps to better describe the spatial organization mode of cells. Mining regions with similar expression patterns from spatial transcriptome data through clustering to interpret the spatial organization mode of cells, that is, identifying spatial domains is one of the most important tasks in spatial transcriptomics.

[0004] Traditional clustering algorithms, such as Louvain and K-means, cannot effectively utilize the available spatial information. Therefore, the clustering results cannot continuously identify tissue regions with obvious hierarchical structures in tissue sections and cannot provide accurate references for downstream analysis. Therefore, there is a need to develop a spatial clustering method for spatial transcriptome data that simultaneously utilizes gene expression profiles and spatial position coordinates.

[0005] In 2021, Jian Hu et al. proposed a deep learning algorithm called SpaGCN in Nature Methods, which integrates gene expression, spatial position, and histological images through a graph convolutional network. The algorithm first constructs a graph representing the relationship between points by combining spatial position and histological images, and then uses a graph convolutional layer to aggregate gene expression information from neighboring points; then an unsupervised iterative clustering algorithm is adopted to cluster points using the aggregated expression matrix.

[0006] In 2021, Edward Zhao et al. proposed an algorithm named BayeSpace in Nature Biotechnology. The BayeSpace algorithm models the low-dimensional representation of the gene expression matrix and introduces a spatial neighbor structure in the prior algorithm through Bayesian statistical methods to encourage adjacent pixel points to belong to the same cluster, thus realizing spatial clustering.

[0007] In 2022, Shihua Zhang et al. proposed a new graph attention autoencoder-based framework STAGATE in Nature Communications. It uses a graph attention autoencoder to automatically learn the weights of the edges between nodes through an attention mechanism while embedding spatial information, considering the spatial similarity of the boundary pixel points in the spatial domain.

[0008] In 2022, Chang Xu et al. proposed a deep neural network framework DeepST in Nucleic Acids Research. It extracts histological image features with a neural network and creates a spatially enhanced gene expression matrix with gene expression and spatial positions, and jointly generates a latent representation of enhanced ST data using a graph convolutional network and a denoising autoencoder.

[0009] The above algorithms all have the following deficiencies:

[0010] First, due to the addition of histological features, while improving the clustering accuracy, it increases the complexity of the model, occupies a large amount of memory, and has a long running time.

[0011] Second, at the same time, some algorithms overcorrect the gene expression features due to overemphasis on spatial information, resulting in overfitting of clustering. Therefore, they cannot identify some fine regions and cannot accurately analyze the biological functions of the identified spatial domains.

[0012] Third, it is unstable, with large differences in the results of each reproduction, and it only performs well on datasets measured by spatial transcriptome sequencing methods based on in situ sequencing, but performs very poorly on imaging-based datasets, and cannot perform extensive analysis on spatial transcriptome datasets. Summary of the Invention

[0013] The purpose of the present invention is to address the above deficiencies of the prior art and propose a spatial domain recognition method based on the extraction of spatial transcriptomics data features, so as to extract the joint features of gene expression profiles and spatial position information in spatial transcriptomics data, improve the generalization ability on both sequencing methods of spatial transcriptomics based on sequencing and imaging, and complete the accurate analysis of the biological functions of spatial domains.

[0014] The technical solution of the present invention is as follows: preprocess the gene expression data and spatial information measured in the spatial transcriptome; construct a gene similarity network and a spatial neighborhood network based on the gene expression feature matrix and spatial information; perform data augmentation on the gene similarity network and the spatial neighborhood network; construct a feature extraction model, and input the augmented data into the model to calculate the contrastive loss and the reconstruction loss; train the model according to the calculated loss, input the non-augmented data into the trained model to obtain a low-dimensional embedding; perform clustering on the low-dimensional embedding to complete the recognition of the spatial domain. The implementation steps are as follows:

[0015] (1) Use spatial transcriptome sequencing technology to simultaneously measure the gene expression value and spatial position coordinates of each pixel point in the required tissue section, and obtain spatial transcriptome data including a pixel-gene expression matrix and the spatial position of each pixel point in the tissue section;

[0016] (2) Preprocess the gene expression matrix of the spatial transcriptome data:

[0017] (2a) Delete the expressed genes in the spatial transcriptome data with gene expression values less than three pixel points;

[0018] (2b) Perform numerical normalization on the deleted data so that the sum of the counts of each cell is the median of all cells, then perform logarithmic transformation on the normalized data, and standardize it to zero mean and unit variance;

[0019] (2c) Perform principal component analysis (PCA) on the standardized data, extract the first n principal components, and generate a feature matrix X of gene expression;

[0020] (3) Construct a spatial neighborhood network:

[0021] (3a) Calculate the Euclidean distance d between each pixel point in the tissue section based on the spatial coordinate information;

[0022] (3b) Select the first k nearest neighbors of each pixel point based on the Euclidean distance d calculated based on the spatial coordinates, and construct an adjacency matrix A representing spatial information;

[0023] (3c) Use the gene expression feature matrix X generated in step (2) as the node attribute feature matrix;

[0024] (3d) Based on the adjacency matrix A representing spatial information and the node attribute feature matrix X, construct a spatial neighborhood network G1(A, X);

[0025] (4) Construct a gene expression similarity network:

[0026] (4a) Based on the gene expression feature matrix X generated in step (2), calculate the Euclidean distance d' between the gene expression values of each pixel point in the tissue section;

[0027] (4b) Based on the Euclidean distance d′ between gene expression values, select the top k nearest neighbors for each pixel point to construct an adjacency matrix B representing gene expression similarity;

[0028] (4c) Based on the adjacency matrix B representing gene expression similarity and the node attribute feature matrix X, construct a gene expression similarity network G2(B, X);

[0029] (5) Data augmentation:

[0030] (5a) For the edges and node attribute features in the spatial neighbor network, perform masking according to a given edge masking probability p that follows a Bernoulli distribution r and a node feature masking probability p m to obtain the enhanced spatial neighbor network G1(A1, X1);

[0031] (5b) For the edges and node attribute features in the gene expression similarity network, perform masking according to a given edge masking probability p that follows a Bernoulli distribution r and a node feature masking probability p m to obtain the enhanced gene expression similarity network G2(B1, X2);

[0032] (6) Construct a feature extraction model for spatial transcriptome data composed of a cascaded encoder f(·), parallel decoders h(·), and projector g(·), and use the weighted sum of the contrastive loss L con and the reconstruction loss L recon as the loss function L;

[0033] (7) Train the feature extraction model for spatial transcriptome data:

[0034] (7a) Input the adjacency matrix A1 and node attribute feature matrix X1 of the enhanced spatial neighbor network G1(A1, X1) and the adjacency matrix B1 and node attribute feature matrix X2 of the gene expression similarity network G2(B1, X2) into the spatial transcriptome feature extraction model. The encoder generates low-dimensional embeddings Z1 and Z2, and the decoder generates reconstructed gene expression feature matrices and

[0035] (7b) Calculate the contrastive loss of the low-dimensional embeddings Z1 and Z2 and the reconstruction loss of the reconstructed gene expression feature matrices and with the node attribute feature matrix X. Update the network parameters of the encoder and decoder according to the calculated losses until the loss function L converges to obtain the trained spatial transcriptome feature extraction model;

[0036] (8) Input the adjacency matrix A and the node attribute feature matrix X of the spatial neighborhood network without data augmentation into the trained spatial transcriptome feature extraction model in (7b) to obtain the joint low-dimensional embedding Z containing spatial information and gene expression;

[0037] (9) Use the Leiden clustering algorithm to cluster the obtained joint low-dimensional embedding Z to obtain regions with consistent gene expression on the tissue section, that is, spatial domains.

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

[0039] 1) Since the present invention constructs a spatial neighborhood network and a gene expression similarity network by combining the spatial information and gene expression profiles of spatial transcriptome data, compared with existing methods, it can better balance spatial information and gene expression information, prevent overfitting of gene expression profiles, and improve the accuracy and robustness of spatial domain recognition.

[0040] 2) Since the present invention only uses gene expression profiles and spatial information and does not add features of histological images, the efficiency of the model is improved, which not only reduces the running time but also can cope with the challenges of larger datasets generated in the future.

[0041] 3) Since the present invention introduces a contrastive loss to train the low-dimensional embedding, making similar samples in the sample more similar and dissimilar samples farther away, it well fits the problem of spatial clustering. Compared with the prior art, the generalization ability on datasets generated by two sequencing methods of spatial transcriptome is improved.

[0042] 4) Since the present invention designs a cascaded model architecture of an encoder and a decoder and considers both contrastive loss and reconstruction loss, it will generate denoised data and better retain the biological significance in the original samples. Brief Description of the Drawings

[0043] Figure 1 is a flowchart for the implementation of the present invention;

[0044] Figure 2 is a schematic diagram of data augmentation in the present invention;

[0045] Figure 3 is a diagram of the feature extraction model of the spatial transcriptome data constructed in the present invention;

[0046] Figure 4 is a visualization diagram of the spatial clustering results using the present invention and the existing STAGATE and DeepST methods respectively. Detailed Embodiments

[0047] The following further elaborates on the embodiments and effects of the present invention with reference to the accompanying drawings.

[0048] Existing spatial transcriptome data includes in-situ hybridization-based imaging techniques and spatial barcoding-based in-situ sequencing techniques. Among them, imaging techniques include STAPmap and MERFISH, and in-situ sequencing techniques include spatialtranscriptomics, 10x Visium, and Slide-seq. In this embodiment, a spatial transcriptome dataset of 151673 slices of the human dorsolateral prefrontal cortex by 10x Visium spatial transcriptome sequencing is taken as an example. This dataset contains 3639 pixels, and each pixel has 33538 genes.

[0049] Refer to Figure 1 , the implementation steps of this example are as follows:

[0050] Step 1, preprocess the gene expression matrix of the spatial transcriptome data.

[0051] 1.1) Obtain the pixel-gene expression matrix data in the spatial transcriptome dataset of 151673 slices of the human dorsolateral prefrontal cortex by 10x Visium spatial transcriptome sequencing, and delete the genes with expression values in less than three pixels in the spatial transcriptome data to filter the data, obtaining the remaining 3639 pixels and 19151 genes after filtering;

[0052] 1.2) Perform median normalization on the filtered transcriptome data, that is, divide each column of data by the median of that column of data, then perform logarithmic transformation on the median-normalized data, and standardize it to zero mean and unit variance;

[0053] 1.3) Perform principal component analysis PCA on the standardized data, extract the first 300 principal components, and generate the feature matrix X of gene expression:

[0054] X = [x1; x2…; x i ;…; x n T

[0055] Among them, [;] represents the concatenation operation, and x i is the gene feature vector of the i-th pixel, i = 1....n, where n is the number of all pixels in the tissue section, and T represents transpose.

[0056] Step 2, construct a spatial neighborhood network.

[0057] 2.1) Obtain the spatial position coordinate data in the spatial transcriptome dataset of 151673 slices of the human dorsolateral prefrontal cortex by 10x Visium spatial transcriptome sequencing, and calculate the Euclidean distance d of each pixel in the spatial position based on the spatial coordinate information: ​

[0058]

[0059] Among them, (a i , b i ) and (a j , b j ) are the spatial coordinates of pixel point i and pixel point j on the tissue section, respectively;

[0060] 2.2) Select the first 5 nearest neighbors of each pixel point based on the Euclidean distance d calculated from the spatial coordinates, and construct an adjacency matrix A representing the spatial information:

[0061]

[0062] Among them, is the element in the i-th row and j-th column of the adjacency matrix A of the spatial neighborhood network. i and j represent two nodes in the spatial neighborhood network. If node i is included in the first 5 nearest neighbors of node j calculated based on the spatial coordinates, then i and j are adjacent; otherwise, i and j are not adjacent. i, j = 1....n, where n = 3639 represents the number of nodes included in the spatial neighborhood network;

[0063] 2.3) Use the gene expression feature matrix X generated in step 1 as the node attribute feature matrix;

[0064] 2.4) Based on the adjacency matrix A representing the spatial information and the node attribute feature matrix X, construct a spatial neighborhood network G1(A, X).

[0065] Step 3, construct a gene expression similarity network.

[0066] 3.1) Based on the gene expression feature matrix X generated in step 1, calculate the Euclidean distance d′ between the gene expression values of each pixel point in the tissue section:

[0067]

[0068] Among them, x jk and x ik are the values of the k-th dimension of the gene expression feature vectors of pixel point i and pixel point j, respectively, where k = 1....m and m = 300 is the dimension of each pixel point's gene expression feature vector;

[0069] 3.2) Based on the Euclidean distance d′ calculated from the gene expression values, select the first 5 nearest neighbors of each pixel point and construct an adjacency matrix B representing the gene expression similarity:

[0070]

[0071] Among them, is the element in the \(i\)-th row and \(j\)-th column of the adjacency matrix \(B\) of the gene expression similarity network. \(i\) and \(j\) represent two nodes in the gene expression similarity network respectively. If node \(i\) is included in the top 5 nearest neighbors of node \(j\) calculated based on the gene expression matrix, then \(i\) and \(j\) are adjacent; otherwise, \(i\) and \(j\) are not adjacent. \(i,j = 1,\cdots,n\), where \(n = 3639\) represents the number of nodes in the gene expression similarity network, which is the same as the number of nodes in the spatial neighborhood network;

[0072] 3.3) Based on the adjacency matrix \(B\) representing gene expression similarity and the node attribute feature matrix \(X\), construct the gene expression similarity network \(G_2(B,X)\).

[0073] Step 4, enhance the edges and node attribute features in the spatial neighbor network \(G_1(A,X)\) and the gene expression similarity network \(G_2(B,X)\).

[0074] To increase the training samples and improve the self-supervised ability of the model, it is necessary to perform data augmentation on the adjacency matrix \(A\) and the node attribute feature matrix \(X\) of the spatial neighbor network, as well as the adjacency matrix \(B\) and the node attribute feature matrix \(X\) of the gene expression similarity network.

[0075] Refer to Figure 2 , and the specific implementation of this step is as follows:

[0076] 4.1) According to each element \(a_{ij}\) in the adjacency matrix \(A\) of the spatial neighborhood network, sample an edge masking matrix ij according to the Bernoulli distribution

[0077]

[0078] where is the element in the \(i\)-th row and \(j\)-th column of the edge masking matrix . If \(a_{ij}\) ij = 1, then the value of is sampled from the Bernoulli distribution \(B(1 - p_{ij}\) r ). If \(a_{ij}\) ij = 0, then where \(p_{ij}\) r = 0.2 is the probability that each edge in the spatial neighborhood network is deleted. \(i\) and \(j\) represent two nodes in the spatial neighborhood network respectively. \(i,j = 1,\cdots,n\), where \(n = 3639\) represents the number of nodes in the spatial neighborhood network;

[0079] 4.2) Multiply the adjacency matrix \(A\) of the spatial neighborhood network with the sampling matrix generated in 4.1) element by element to obtain the enhanced adjacency matrix \(A_1\):

[0080]

[0081] In the formula, The operator represents multiplying the adjacency matrix A and the sampling matrix in the spatial neighborhood network element by element, where a ij is the element in the i-th row and j-th column of the adjacency matrix A of the spatial neighborhood network, and is the element in the i-th row and j-th column of the sampling matrix

[0082] 4.3) Sample a random vector according to the Bernoulli distribution B(1 - p m ), and generate a node feature masking vector with the same dimension as the gene feature vector where pm = 0.3 is the probability that the value in each node feature vector in the spatial neighborhood network is deleted;

[0083] 4.4) Multiply the node attribute feature matrix X of the spatial neighbor network and the node feature masking vector generated in 4.3) element by element to obtain the enhanced node attribute feature matrix X1:

[0084]

[0085] where [;] represents the concatenation operation, and x i is the i-th row of X1, representing the gene feature vector on node i in the spatial neighbor network, i = 1....n, and n = 3639 is the number of nodes in the spatial neighborhood network;

[0086] 4.5) According to each element b ij in the adjacency matrix B of the gene expression similarity network, sample an edge masking matrix R ∈ {0, 1} N×N :

[0087]

[0088] In the formula, is the element in the i-th row and j-th column of the edge masking matrix R. If b ij = 1, then the value of R ij is sampled from the Bernoulli distribution B(1 - p r ), and if b ij = 0, then R ij = 0, where p r ′ = 0.2 is the probability that each edge in the gene expression similarity network is deleted. i and j respectively represent two nodes in the gene expression similarity network, and i, j = 1....n, where n = 3639 represents the number of nodes included in the gene expression similarity network, which is the same as the number of nodes in the spatial neighborhood network;

[0089] 4.6) Multiply the adjacency matrix B of the gene expression similarity network element - by - element with the sampling matrix R generated in 4.5) to obtain the adjacency matrix B1 of the enhanced gene expression similarity network:

[0090]

[0091] In the formula, The operator represents multiplying the adjacency matrix B of the gene expression similarity network and the sampling matrix R element - by - element. b ij is the element in the i - th row and j - th column of the adjacency matrix B of the gene expression similarity network, and R ij is the element in the i - th row and j - th column of the sampling matrix R;

[0092] 4.7) Sample a random vector according to the Bernoulli distribution B(1 - p m ), and generate a node feature masking vector m with the same dimension as the gene feature vector, where p m ′ = 0.3 is the probability that the value in each node feature vector in the gene expression similarity network is deleted;

[0093] 4.8) Multiply the node attribute feature matrix X of the gene expression similarity network element - by - element with the node feature masking vector m generated in 4.7) to obtain the node attribute feature matrix X2 of the enhanced gene expression similarity network:

[0094]

[0095] Among them, [;] represents the concatenation operation, and x i is the i - th row in X2, representing the gene feature vector on node i in the gene expression similarity network.

[0096] Step 5, construct a feature extraction model for spatial transcriptome data.

[0097] Refer to Figure 3 The specific implementation of this step is as follows:

[0098] 5.1) Establish an encoder composed of an input GCN layer and two cascaded hidden GCN layers. Its input dimension is 300 - dimensional of the transcriptome gene feature dimension, the first hidden GCN layer is 256 - dimensional, and the second hidden GCN layer is 128 - dimensional. Use the PReLU function as the activation function between each GCN layer;

[0099] 5.2) Establish a decoder composed of an input fully - connected layer and three cascaded hidden fully - connected layers. Its input fully - connected layer is 128 - dimensional, the first hidden fully - connected layer is 128 - dimensional, the second hidden fully - connected layer is 256 - dimensional, and the third hidden fully - connected layer dimension is 300 - dimensional of the transcriptome gene feature dimension. Use the Relu function as the activation function between each fully - connected layer;

[0100] 5.3) Establish a projector composed of a cascaded input fully-connected layer and a hidden fully-connected layer. The input fully-connected layer is 128-dimensional, and the hidden fully-connected layer is 128-dimensional. No activation function is set between each layer;

[0101] 5.4) Cascade the encoder with the decoder and the projector respectively to form a feature extraction model for spatial transcriptome data;

[0102] 5.5) Let the loss function L of the feature extraction model be the weighted sum of the contrast loss and the reconstruction loss, expressed as follows:

[0103] L = λ con L con + λ recon L recon

[0104] where λ con = 1, λ recon = 0.01 are hyperparameters measuring the weights of the contrast loss and the reconstruction loss respectively. L recon represents the reconstruction loss, and L con represents the contrast loss.

[0105] Step 6, train the feature extraction model for spatial transcriptome data.

[0106] 6.1) Respectively extract the node low-dimensional embeddings Z1 of G1(A1, X1) and Z2 of G2(B1, X2) through the encoder f(·):

[0107] Z1 = f(X1, A1) = GC k+1 (GC k (X1, A1), A1)

[0108] Z2 = f(X2, B1) = GC k+1 (GC k (X2, B1), B1),

[0109] where GC k (·) represents the k-th layer network of the encoder. X1 and A1 are the node feature matrix and the adjacency matrix of the spatial neighborhood network respectively. X2 and B1 are the node feature matrix and the adjacency matrix in the gene expression similarity network respectively, and k = 1;

[0110] 6.2) Respectively use the two low-dimensional embeddings Z1 and Z2 obtained in step 6.1) as the inputs of the decoder h(·) to obtain the reconstructed gene expression feature matrix of G1(A1, X1) and the reconstructed gene expression feature matrix of G2(B1, X2)

[0111]

[0112]

[0113] 6.3) Input the two low-dimensional embeddings Z1 and Z2 generated in step 6.1) into the projector g(·) respectively to obtain the contrastive loss low-dimensional embedding Z′1 of Z1 and the contrastive loss low-dimensional embedding Z′2 of Z2 respectively:

[0114] Z′1 = g(Z1)

[0115] Z′2 = g (Z2);

[0116] 6.4) Calculate the contrastive loss l(z′ 1i , z′ 2i ) for each node i and other nodes k according to the results of step 6.3):

[0117]

[0118] In the formula, θ(·) represents the cosine similarity distance, and τ is a given hyperparameter; z′ 1i is the vector of the i-th row of Z′1, which represents the output of the projector when the pixel point i takes G1(A1, X1) as the input; z′ 2i is the vector of the i-th row of Z′2, which represents the output of the projector when the pixel point i takes G2(B1, X2) as the input, i, k = 1....n, and n = 3639 is the number of all pixel points in the tissue section;

[0119] 6.5) Calculate the contrastive loss L con of the entire model according to the contrastive loss of each node obtained in step 6.4):

[0120]

[0121] 6.6) Calculate the reconstruction loss L and according to the reconstructed gene feature matrices generated in step 6.2): recon :

[0122]

[0123] In the formula, x i is the vector of the i-th row of X, which represents the gene feature of the pixel point i; is the vector of the i-th row, which represents the gene feature of the pixel point i reconstructed by G1(A1, X1); is the vector of the i-th row, which represents the gene feature of the pixel point i reconstructed by G2(B1, X2);

[0124] 6.7) Calculate the loss function L according to the contrastive loss L con and the reconstruction loss L recon as follows:

[0125] L = λ con L con + λ recon L recon

[0126] 6.8) Update the network parameters of the encoder and decoder according to the loss function L obtained in step 6.7) until the loss function L converges, and obtain the trained spatial transcriptome feature extraction model.

[0127] Step 7: Input the adjacency matrix A and the node attribute feature matrix X of the spatial neighborhood network without data augmentation into the spatial transcriptome feature extraction model trained in step 6 to obtain the joint low-dimensional embedding Z containing spatial information and gene expression;

[0128] Step 8: Use the Leiden clustering algorithm to cluster the joint low-dimensional embedding obtained in step 7.

[0129] 8.1) Calculate the neighbors of each pixel point according to the joint low-dimensional embedding Z extracted in step 7, construct a neighborhood graph, and save the neighborhood label l';

[0130] 8.2) Reduce the dimension of the joint low-dimensional embedding Z through the UMAP algorithm to obtain the reduced-dimensional embedding Z';

[0131] 8.3) Obtain the clustering label l through the Leiden algorithm according to the neighborhood label l' obtained in step 8.1) and the low-dimensional embedding Z' obtained in step 8.2);

[0132] 8.4) Visualize the clustering label l and the low-dimensional embedding Z' through UMAP, stain each pixel point on the tissue section according to the clustering label l, and regard the pixel points of the same color as a domain, that is, realize the identification of the spatial domain.

[0133] The technical effects of the present invention are described below in combination with simulation experiments.

[0134] I. Simulation conditions:

[0135] The computer hardware CPU of the simulation experiment is Intel Core(TM)i7-8700, and the computer hardware memory is 32G;

[0136] Computer software: Python3.8 integrated development software on the WINDOWS10 system.

[0137] II. Simulation content:

[0138] Simulation 1: The present invention and six existing methods, namely SEDR, STAGATE, DeepST, scanpy, stlearn, and SpaGCN, were used to perform spatial clustering on datasets generated by two spatial transcriptome sequencing methods, namely the spatial transcriptome dataset of 12 sections of the human dorsolateral prefrontal cortex layer (DLPFC) of 10x Visium based on in situ sequencing and the spatial transcriptome dataset of the mouse visual cortex of STARmap based on imaging. The adjusted Rand index (ARI) was used as an evaluation index for the spatial clustering results of each method. The results are shown in Table 1:

[0139] Table 1 Evaluation of the present invention and six existing methods in the labeled dataset

[0140]

[0141] The six existing spatial domain recognition methods are sourced as follows:

[0142] SEDR, Ling S, Huazhu F, et al. Unsupervised Spatially Embedded Deep Representation of Spatial Transcriptomics[J]. bioRxiv, 2021.

[0143] STAGATE, Dong K, Zhang S. Deciphering spatial domains from spatially resolved transcriptomics with an adaptive graph attention auto-encoder[J]. Nature communications, 2022, 13(1): 1 - 12.

[0144] DeepST, Xu C, Jin X, Wei S, et al. DeepST: identifying spatial domains in spatial transcriptomics by deep learning[J]. Nucleic Acids Research, 2022.

[0145] Scanpy, Wolf F A, Angerer P, Theis F J. SCANPY: large-scale single-cell gene expression data analysis[J]. Genome Biology, 2018, 19(1): 1-5.

[0146] Stlearn, Pham D, Tan X, Xu J, et a1. stLearn: integrating spatial location, tissue morphology and gene expression to find cell types, cell-cell interactions and spatial trajectories within undissociated tissues[J]. bioRxiv, 2020.

[0147] SpaGCN, LiM, Hu J, Li X, et al. SpaGCN: Integrating gene expression, spatial location and histology to identify spatial domains and spatially variable genes by graph convolutional network[J]. Nature Methods, 2021, 18(10): 1342-1351.

[0148] As can be seen from Table 1, the present invention has better results than other methods on the 12 datasets of the DLPFC of 10x Visium, and the mean value is also higher than that of other methods. On the mouse visual cortex dataset of STAPmap, the performance of the present invention and STAGATE is significantly higher than that of other methods, but the present invention has higher accuracy than STAGATE. The simulation results show that the present invention maintains high accuracy and has good generalization ability in both in situ sequencing-based datasets and imaging-based datasets.

[0149] Simulation 2: The present invention and three existing methods, DeepST, SEDR, and STAGATE, were used to perform spatial clustering on the spatial transcriptome datasets of mouse brain slices and human breast cancer in 10x Visium. The Silhouette Coefficient score and Davies-Bouldin score were used as evaluation metrics for the spatial clustering results of each method. The results are shown in Table 2:

[0150] Table 2 Evaluation of the present invention and three existing methods in unlabeled datasets

[0151]

[0152] As can be seen from Table 2, on the spatial transcriptome dataset of the mouse brain in 10x Visium, the performance of the present invention and STAGATE is significantly higher than that of other methods, but the indicators of the present invention are slightly higher than those of STAGATE. On the spatial transcriptome dataset of human breast cancer in 10x Visium, the present invention has a greater advantage over other methods. The simulation results show that on some unlabeled datasets that require fine recognition, the clustering results of the present invention are better and better retain the biological significance in the original samples.

[0153] Simulation 3: The present invention and two existing methods, DeepST and STAGATE, were used to identify spatial domains through spatial clustering on the dataset of the coronal plane of the mouse brain in 10x Visium, and each pixel point on the histological section was stained with the clustering results. The results are as Figure 4 shown. Among them Figure 4 (a) represents the visualization diagram of the spatial clustering of the present invention, Figure 4 (b) represents the visualization diagram of the spatial clustering of STAGATE, Figure 4 (c) represents the visualization diagram of the spatial clustering of DeepST.

[0154] From Figure 4 it can be seen that the existing methods STAGATE and DeepST cannot accurately identify the spatial domains on the coronal section of the mouse brain and cannot clearly represent the differences between each domain. Especially in the hippocampal region of the dataset, it is not accurately divided into three clusters, while the spatial clustering results of the present invention are more in line with biological significance. The simulation results show that the features extracted by the present invention do not overfit the gene expression profile, improving the accuracy and robustness of spatial domain recognition.

Claims

1. A spatial domain recognition method based on feature extraction of spatial transcriptomics data, characterized in that It includes the following steps: (1) Use spatial transcriptome sequencing technology to simultaneously measure the gene expression value and spatial position coordinates of each pixel in the required tissue section, and obtain spatial transcriptome data including a pixel-gene expression matrix and the spatial position of each pixel in the tissue section; (2) Preprocess the gene expression matrix of the spatial transcriptome data: (2a) Delete the expressed genes with gene expression values less than three pixels in the spatial transcriptome data; (2b) Perform numerical normalization on the deleted data so that the sum of the counts of each cell is the median of all cells, then perform logarithmic transformation on the normalized data, and standardize it to zero mean and unit variance; (2c) Perform principal component analysis PCA on the standardized data, extract the first n principal components, and generate a feature matrix X of gene expression; (3) Construct a spatial neighborhood network: (3a) Calculate the Euclidean distance d in the spatial position between each pixel in the tissue section based on the spatial coordinate information; (3b) Select the first k nearest neighbors of each pixel based on the Euclidean distance d calculated based on the spatial coordinates, and construct an adjacency matrix A representing spatial information; (3c) Use the gene expression feature matrix X generated in step (2) as the node attribute feature matrix; (3d) Based on the adjacency matrix A representing spatial information and the node attribute feature matrix X, construct a spatial neighborhood network G1(A, X); (4) Construct a gene expression similarity network: (4a) Calculate the Euclidean distance d' between the gene expression values of each pixel in the tissue section based on the gene expression feature matrix X generated in step (2); (4b) Select the first k nearest neighbors of each pixel based on the Euclidean distance d' calculated based on the gene expression values, and construct an adjacency matrix B representing gene expression similarity; (4c) Based on the adjacency matrix B representing gene expression similarity and the node attribute feature matrix X, construct a gene expression similarity network G2(B, X); (5) Data augmentation: (5a) Mask the edge and node attribute features in the spatial neighbor network according to a given edge masking probability p that conforms to the Bernoulli distribution r and node feature masking probability p m to obtain the enhanced spatial neighbor network G1(A1, X1); (5b) Mask the edge and node attribute features in the gene expression similarity network according to the given edge masking probability p r ' and node feature masking probability p m ' to obtain the enhanced gene expression similarity network G2(B1, X2); (6)Construct a feature extraction model for spatial transcriptome data composed of the cascade of the encoder f(·) with the decoder h(·) and the projector g(·) respectively, and use the weighted sum of the contrastive loss L con and the reconstruction loss L recon as the loss function L; (7) Train the feature extraction model of the spatial transcriptome data: (7a) Input the adjacency matrix A1 and node attribute feature matrix X1 of the spatially neighboring network G1(A1, X1) after data augmentation, as well as the adjacency matrix B1 and node attribute feature matrix X2 of the gene expression similarity network G2(B1, X2) into the spatial transcriptome feature extraction model. The encoder generates low-dimensional embeddings Z1 and Z2, and the decoder generates a reconstructed gene expression feature matrix and (7b) Calculate the contrastive loss between the low-dimensional embeddings Z1 and Z2 and the reconstruction loss of the gene expression feature matrix, and update the network parameters of the encoder and decoder according to the calculated loss until the loss function L converges to obtain a trained spatial transcriptome feature extraction model; and the reconstruction loss with the node attribute feature matrix X, and update the network parameters of the encoder and decoder according to the calculated loss until the loss function L converges to obtain a trained spatial transcriptome feature extraction model; (8) Input the adjacency matrix A and the node attribute feature matrix X of the spatial neighborhood network without data augmentation into the trained spatial transcriptome feature extraction model in (7b) to obtain a joint low-dimensional embedding Z containing spatial information and gene expression; (9) Use the Leiden clustering algorithm to cluster the obtained joint low-dimensional embedding Z to obtain regions with consistent gene expression on the tissue section, that is, spatial domains.

2. The method according to claim 1, wherein The gene expression feature matrix X generated in step (2c) is expressed as follows: X = [x1; x2…; x i ; …; x n T ​ Among them, [;] represents the splicing operation, and x i is the gene feature vector of i pixel points, where i = 1....n, and n is the number of all pixel points in the tissue section, and T represents the transpose.

3. The method according to claim 1, wherein The Euclidean distance d in the spatial position between each pixel in the tissue section calculated in step (3a) is given by the following formula: Among them, (a i , b i ) and (a j , b j ) are the spatial coordinates of pixel point i and pixel point j on the tissue section, respectively.

4. The method according to claim 1, characterized in that, The adjacency matrix A of the spatial neighborhood network constructed in step (3b) is expressed as follows: Among them, is the element in the i-th row and j-th column of the adjacency matrix A of the spatial neighborhood network. i and j respectively represent two nodes in the spatial neighborhood network. If node i is included in the first k nearest neighbors calculated based on spatial coordinates of node j, then i and j are adjacent; otherwise, i and j are not adjacent. i, j = 1....n, where n represents the number of nodes included in the spatial neighborhood network.

5. The method according to claim 1, wherein The Euclidean distance d' between the gene expression values of each pixel in the tissue section calculated in step (4a) is given by the following formula: where x jk and x ik are the values of the k-th dimension of the gene expression feature vectors of pixel i and pixel j respectively, where k = 1....m and m is the dimension of the gene expression feature vector of each pixel.

6. The method according to claim 1, characterized in that The adjacency matrix B of the gene expression similarity network constructed in step (4b) is expressed as follows: wherein, is the element at the i-th row and j-th column in the adjacency matrix B of the gene expression similarity network. i and j respectively represent two nodes in the gene expression similarity network. If node i is included in the top k nearest neighbors calculated based on the gene expression matrix of node j, then i and j are adjacent; otherwise, i and j are not adjacent. i, j = 1....n, and n represents the number of nodes included in the gene expression similarity network.

7. The method according to claim 1, characterized in that, In step (5a), the edges and node attribute features in the spatial neighbor network G1(A, X) are masked according to probability, and it is achieved as follows: (5a1) For each element a in the adjacency matrix A of the spatial neighborhood network ij , sample an edge masking matrix according to the Bernoulli distribution which is represented as follows: In the formula, is the element in the \(i\)-th row and \(j\)-th column of the edge masking matrix . If \(a\) ij = 1, then the value of is sampled from the Bernoulli distribution \(B(1 - p\) r ). If \(a\) ij = 0, then . Among them, \(p\) r is the probability that each edge in the spatial neighborhood network is deleted. \(i\) and \(j\) respectively represent two nodes in the spatial neighborhood network, \(i, j = 1,\cdots,n\), and \(n\) represents the number of nodes contained in the spatial neighborhood network; (5a2) Multiply the adjacency matrix A of the spatial neighborhood network element - by - element with the sampling matrix generated in (5a1) to obtain the enhanced adjacency matrix A1: In the formula, The operator represents multiplying the adjacency matrix A in the spatial neighborhood network and the sampling matrix element by element, where a ij is the element at the i-th row and j-th column in the adjacency matrix A of the spatial neighborhood network, and is the element at the i-th row and j-th column in the sampling matrix (5a3) Sample a random vector according to the Bernoulli distribution B(1 - p m ) to generate a node feature masking vector with the same dimension as the gene feature vector where p m is the probability that the value in each node feature vector in the spatial neighborhood network is deleted; (5a4) Multiply the node attribute feature matrix X of the spatial neighbor network with the node feature masking vector generated in (5a3). Element-wise multiplication gives the enhanced node attribute feature matrix X1: Among them, [;] represents the splicing operation, x i is X1 的 The i-th row represents the gene feature vector on node i in the spatial neighbor network.

8. The method according to claim 1, wherein In step (5b), the edges and node attribute features in the gene expression similarity network G2(B, X) are masked according to probabilities, which is achieved as follows: (5b1) For each element \(b\) in the adjacency matrix \(B\) of the gene expression similarity network ij , sample an edge masking matrix \(R\in\{0, 1\}\) according to the Bernoulli distribution N×N , which is expressed as follows: Wherein, is the element in the i-th row and j-th column of the edge masking matrix R. If b ij = 1, then the value of R ij is sampled from the Bernoulli distribution B(1 - p r ). If b ij = 0, then R ij = 0. Wherein, p r ' is the probability that each edge in the gene expression similarity network is deleted. i and j respectively represent two nodes in the gene expression similarity network, i, j = 1....n, and n represents the number of nodes contained in the gene expression similarity network; (5b2) Multiply the adjacency matrix B of the gene expression similarity network element-wise with the sampling matrix R generated in (5b1) to obtain the adjacency matrix B1 of the enhanced gene expression similarity network: In the formula, The operator represents the element-wise multiplication of the adjacency matrix B of the gene expression similarity network and the sampling matrix R, and b ij is the element in the i-th row and j-th column of the adjacency matrix B of the gene expression similarity network, and R ij is the element in the i-th row and j-th column of the sampling matrix R; (5b3) Sample a random vector according to the Bernoulli distribution B(1 - p m ) to generate a node feature masking vector m with the same dimension as the gene feature vector, where p m ′ is the probability that the value in each node feature vector in the gene expression similarity network is deleted; (5b4) Multiply the node attribute feature matrix X of the gene expression similarity network element-wise with the node feature masking vector m generated in (5a3) to obtain the node attribute feature matrix X2 of the enhanced gene expression similarity network: Among them, [;] represents the splicing operation, and x i is the i-th row in X2, representing the gene feature vector on node i in the gene expression similarity network.

9. The method according to claim 1, wherein In the spatial transcriptome data feature extraction model constructed in step (6), its encoder, decoder, projector parameters, and loss function are as follows: The encoder f(·) is composed of a cascaded input graph neural network GCN layer and two hidden GCN layers. Its input dimension is the number of transcriptome gene features, the first hidden GCN layer is 256-dimensional, and the second hidden GCN layer is 128-dimensional. The PRelu function is used as the activation function between each layer of GCN; The decoder h(·) is composed of a cascaded input fully connected layer and three hidden fully connected layers. Its input fully connected layer is 128-dimensional, the first hidden fully connected layer is 128-dimensional, the second hidden fully connected layer is 256-dimensional, and the dimension of the third hidden fully connected layer is the number of transcriptome gene features. The Relu function is used as the activation function between each layer of fully connected layers; The projector g(·) is composed of a cascaded input fully connected layer and a hidden fully connected layer. Its input fully connected layer is 128-dimensional, and the hidden fully connected layer is 128-dimensional. No activation function is set between each layer; The loss function is the weighted sum of the contrastive loss and the reconstruction loss, which is expressed as: L = λ con L con + λ recon L recon Among them, λ con and λ recon are hyperparameters used to measure the weights of the contrast loss and the reconstruction loss respectively. L recon represents the reconstruction loss, and L con represents the contrast loss.

10. The method according to claim 1, wherein In step (7a), a low-dimensional embedding is generated by the encoder, and the gene expression feature matrix is reconstructed by the decoder, which is achieved as follows: (7a1) The encoder f(·) extracts the node low-dimensional embedding Z1 of G1(A1, X1) and the node low-dimensional embedding Z2 of G2(B1, X2): Z1 = f(X1, A1) = GC k+1 (GC k (X1, A1), A1) Z2 = f(X2, B1) = GC k+1 (GC k (X2, B1), B1), Among them, GC k (·) represents the k-th layer network of the encoder. X1 and A1 are the node feature matrix and adjacency matrix of the spatial neighborhood network respectively. X2 and B1 are the node feature matrix and adjacency matrix in the gene expression similarity network respectively, where k = 1; (7a2) Take the two low-dimensional embeddings Z1 and Z2 obtained in (7a1) as the inputs of the decoder h(·) respectively, and obtain the gene expression feature matrix reconstructed by G1(A1, X1) and the gene expression feature matrix reconstructed by G2(B1, X2) 11. The method according to claim 1, wherein In step (7b), the contrastive loss and the reconstruction loss are calculated, which is achieved as follows: (7b1) Input the low-dimensional embeddings generated in step (7a) into the projector g(·) respectively to obtain the low-dimensional embeddings Z′1 and Z′2 of the contrastive loss of Z1 and Z2 respectively: Z′1 = g(Z1) Z′2 = g(Z2); (7b2) Calculate the contrastive loss l(z′ 1i , z′ 2i ) between each node i and other nodes k according to the result of (7b1): where, θ(·) represents the cosine similarity distance, and τ is a given hyperparameter; z′ 1i is the vector of the i-th row of Z′1, which represents the output of the projector when the pixel point i takes G1(A1, X1) as the input; z′ 2i is the vector of the i-th row of Z′2, which represents the output of the projector when the pixel point i takes G2(B1, X2) as the input, where i, k = 1....n, and n is the number of all pixel points in the tissue section; Calculate the contrastive loss \(L\) of the entire model based on the contrastive loss of each node obtained in (7b2). con : (7b4) According to the reconstructed gene feature matrix generated in step (7a) and calculate the reconstruction loss L recon : where x i is the vector of the i-th row of X, representing the gene feature of pixel point i; is the vector of the i-th row, representing the gene feature of pixel point i reconstructed by G1(A1, X1); is the vector of the i-th row, representing the gene feature of pixel point i reconstructed by G2(B1, X2); ( 7b5) Calculate the loss function L according to the contrastive loss L con and the reconstruction loss L recon : L = λ con L con + λ recon L recon Among them, λ con and λ recon are hyperparameters used to measure the weights of the contrastive loss and the reconstruction loss, respectively.

12. The method according to claim 1, wherein In step (9), the Leiden clustering algorithm is used to cluster the jointly obtained low-dimensional embedding Z after training, which is achieved as follows: (9a) Calculate the neighbors of each pixel point according to the jointly obtained low-dimensional embedding Z extracted in step (8), construct a neighborhood graph, and save the neighborhood label l′; (9b) Reduce the dimension of the jointly obtained low-dimensional embedding Z through the UMAP algorithm to obtain the reduced-dimensional embedding Z′; (9c) Obtain the clustering label l through the Leiden algorithm according to the neighborhood label l′ in (9a) and the embedding Z′ in (9b); (9d) Perform UMAP visualization on the clustering label l and the low-dimensional embedding Z′, stain each pixel point on the tissue section according to the clustering label l, and regard the pixel points of the same color as a domain, that is, the recognition of the spatial domain is achieved.

Images