Spatial transcriptome data clustering method

Abstract

The invention relates to the field of spatial transcriptome data analysis, in particular to a spatial transcriptome data clustering method, which comprises the following steps: acquiring and preprocessing original data; constructing a kernel density estimation function according to the spatial density and the gene expression similarity, and dynamically adjusting the neighborhood range of each node based on kernel density estimation; calculating a low-frequency coefficient and a high-frequency coefficient by using graph wavelet transform, and obtaining a node importance weight in combination with a time-space activity score; weighting neighborhood features according to node importance weights, constructing a cost matrix in combination with feature distances and distribution differences, and quantifying transmission cost from a source domain to a target domain; a regularized transmission plan is solved by using a Sinkhorn algorithm, cross-domain information is fused through Barrier mapping, and reconstruction features of a target domain are generated; and clustering by using the reconstruction features to obtain a clustering result. According to the method, the problem of insufficient utilization of spatial information of a traditional clustering method is solved, and a more accurate and more adaptive solution is provided for single cell spatial omics analysis.

Description

Technical Field

[0001] The present invention relates to the field of spatial transcriptome data analysis, and in particular to a spatial transcriptome data clustering method. Background Art

[0002] Spatial transcriptome data clustering involves dividing different locations in a tissue into groups with similar molecular characteristics based on the similarity of gene expression profiles, thereby revealing the biological functions, cell type distribution, or state heterogeneity of different regions in the tissue. Non-spatial clustering methods use traditional clustering methods such as K-means and Louvain algorithms, which are limited to a small number of spots, and the clustering results may be discontinuous in tissue sections. Among the latest spatial clustering methods, dynamic graph structure learning can adjust the neighborhood range, but lacks a mechanism to optimize the information transfer weights in heterogeneous density regions, and still suffers from significant information loss in low-density areas. While optimal transmission theory can model inter-node associations, its sole application relies on a fixed spatial topological structure and has limited adaptability to spatially heterogeneous density distributions. Existing technologies have not yet addressed the problem of collaborative modeling of these two methods in spatial transcriptome data, resulting in information redundancy or loss in heterogeneous regions. Summary of the Invention

[0003] In order to solve the collaborative modeling problem of dynamic graph structure learning and optimal transmission theory, the present invention proposes a spatial transcriptome data clustering method, which includes the following steps:

[0004] Obtain raw data and perform preprocessing;

[0005] According to the spatial density and gene expression similarity, a kernel density estimation function is constructed, and the density threshold τ is determined based on the kernel density estimation function, and the combined benchmark neighborhood radius r base Dynamically adjust the neighborhood range of each node;

[0006] Graph wavelet transform is used to calculate low-frequency coefficients and high-frequency coefficients, and the node importance weight is obtained by combining the spatiotemporal activity score; wherein the node in the present invention refers to the spot in the spatial transcriptome data, and each node includes two-dimensional spatial position information and the corresponding gene expression feature vector;

[0007] The neighborhood features are weighted according to the node importance weights, and a cost matrix is constructed by combining feature distance and distribution difference to quantify the transmission cost from the source domain to the target domain.

[0008] The Sinkhorn algorithm is used to solve the regularized transmission plan, and the cross-domain information is fused through Barycentric mapping to generate the reconstructed features of the target domain.

[0009] Clustering is performed using the reconstructed features to obtain clustering results.

[0010] This paper uses kernel density estimation to achieve adaptive adjustment of neighborhood ranges. By introducing probabilistic polarization optimal transport and Barycentric Wasserstein paths, it significantly improves the alignment accuracy of spatial and gene expression information. Through these technological innovations, this paper addresses the underutilization of spatial information in traditional clustering methods and provides a more accurate and adaptive solution for single-cell spatial omics analysis. BRIEF DESCRIPTION OF THE DRAWINGS

[0011] Figure 1 This is a flow chart of a spatial transcriptome data clustering method of the present invention;

[0012] Figure 2 The figure shows the clustering results of one case of human breast cancer data set based on the present invention. DETAILED DESCRIPTION

[0013] The following will clearly and completely describe the technical solutions in the embodiments of the present invention in conjunction with the accompanying drawings. Obviously, the described embodiments are only part of the embodiments of the present invention, not all of the embodiments. Based on the embodiments of the present invention, all other embodiments obtained by ordinary technicians in this field without making creative efforts are within the scope of protection of the present invention.

[0014] The present invention proposes a spatial transcriptome data clustering method, such as Figure 1 , including the following steps:

[0015] Obtain raw data and perform preprocessing;

[0016] According to the spatial density and gene expression similarity, a kernel density estimation function is constructed, and the density threshold τ is determined based on the kernel density estimation function, and the combined benchmark neighborhood radius r base Dynamically adjust the neighborhood range of each node;

[0017] Graph wavelet transform is used to calculate low-frequency coefficients and high-frequency coefficients, and the node importance weight is obtained by combining the spatiotemporal activity score;

[0018] The neighborhood features are weighted according to the node importance weights, and a cost matrix is constructed by combining feature distance and distribution difference to quantify the transmission cost from the source domain to the target domain.

[0019] The Sinkhorn algorithm is used to solve the regularized transmission plan, and the cross-domain information is fused through Barycentric mapping to generate the reconstructed features of the target domain.

[0020] Clustering is performed using the reconstructed features to obtain clustering results.

[0021] In order to solve the collaborative modeling problem of dynamic graph structure learning and optimal transmission theory, the present invention proposes a neighborhood information fusion mechanism that combines dynamic graph structure learning with optimal transmission theory. This mechanism generates a dynamic neighborhood range through kernel density estimation, expands the neighborhood range in low-density areas, and shrinks the connection radius in high-density areas. The optimal transmission theory is simultaneously introduced to quantify the topological changes of neighborhood expansion / contraction into a differentiable transmission cost function, thereby realizing the joint optimization of graph structure and information transfer weights. The system shows significant improvements in redundant connection suppression, information integrity retention, and error control, verifying that the synergistic effect of dynamic graph learning and optimal transmission theory can overcome the inherent limitations of a single technology and achieve more robust spatial topological modeling.

[0022] As a preferred embodiment, the present invention provides a spatial transcriptome data clustering method for clustering spatial transcriptome data containing spatial location information and gene expression information, dividing different locations in a tissue into groups with similar molecular characteristics, thereby revealing the biological functions, cell type distribution, or state heterogeneity of different regions in the tissue. The specific implementation process includes the following steps:

[0023] Step 1: Preprocess the raw data, including selecting highly variable genes, normalization, and log normalization to correct the data sequencing depth;

[0024] Step 2: Quantify the spatiotemporal importance of nodes and construct a dynamic graph structure. Spatial importance is determined by constructing a kernel density estimation function based on spatial density and gene expression similarity. The neighborhood range of each node is dynamically adjusted based on the kernel density estimation to achieve adaptive adjustment of the neighborhood radius. Temporal importance is determined by calculating low-frequency and high-frequency coefficients using graph wavelet transform, and then combining the spatiotemporal activity scores to obtain the node importance weight.

[0025] Step 3: Neighborhood information fusion driven by optimal transmission. This includes weighting neighborhood features based on node spatiotemporal importance, constructing a cost matrix based on feature distance and distribution differences, quantifying the transmission cost from the source domain to the target domain, applying polarization constraints to enhance intra-class transmission and suppress inter-class interference, solving a regularized transmission plan using the Sinkhorn algorithm, and fusing cross-domain information through Barycentric mapping to generate reconstructed features for the target domain.

[0026] Step 4: Use the fused reconstructed features to perform traditional clustering methods such as K-means;

[0027] Step 5: Evaluate the clustering results and adjust parameters to optimize the model.

[0028] Based on the joint measurement of spatial coordinates and gene expression similarity, a kernel density estimation function is constructed:

[0029]

[0030] Among them, ρ(x i ) is the kernel density estimation function; h is the bandwidth parameter; n is the total number of nodes; d is the dimension of the spatial coordinate; x i Represents the two-dimensional spatial position information of the i-th node; |||| represents the norm; sim(g i ,g j ) is the cosine similarity of gene expression, g i is the gene expression feature vector of the i-th node.

[0031] Furthermore, the density threshold τ is used to dynamically divide high / low density areas. When the kernel density estimation function of a node is less than the density threshold, it is a low density area. The neighborhood radius of the node is set to N×r base , N>1; otherwise it is a high-density area, and the neighborhood radius of the node is set to M×r base , 1>M>0, reference neighborhood radius r base is the median of the distance set between any two nodes. In this invention, the node is the spot in the spatial transcriptome data. Each node includes two-dimensional spatial position information and the corresponding gene expression feature vector. The reference neighborhood radius r base The calculation includes:

[0032] According to the two-dimensional spatial position information of the node {x1,x2,...,x n}Find the Euclidean distance between any two nodes, the two-dimensional spatial position information x of the node i and the two-dimensional spatial position information x of the node j The distance between them is d ij , where i∈{1,2,…,n} and j∈{1,2,…,n} and i≠j, n is the number of nodes;

[0033] The median of all distances obtained is used as the base neighborhood radius r base .

[0034] Similarly, the density threshold τ is the number of bits in the kernel density estimation function set of the node, and its setting process includes:

[0035] According to the two-dimensional spatial position information set of the node {x1,x2,...,x n} and the node gene expression feature vector set {g1,g2,...,g n}, calculate the kernel density estimation function of each node;

[0036] The median of all kernel density estimation functions obtained is used as the density threshold τ.

[0037] As a preferred embodiment, the neighborhood radius is extended to r in low-density areas.low =2r base , the high-density area shrinks to r high =0.5r base , to achieve adaptive adjustment of the neighborhood radius.

[0038] Calculate the rate of change of the low-frequency coefficients:

[0039]

[0040] Where Δa t,j is the rate of change of the low-frequency coefficient, ∈ is the numerical stability constant, which is 10 -7 .

[0041] Extract high-frequency coefficients from the eigenvector:

[0042] b t,j =||W high ·g j ||2

[0043] Among them, W high is a high-pass wavelet filter.

[0044] Calculate the time importance weight of a node:

[0045] τ t,j =γ·Δa t,j +(1-γ)·b t,j

[0046] β t,j =δ·σ(τ t,j )

[0047] Among them, τ t,j is the spatiotemporal activity score of node j at time t; γ is the mixing coefficient; β t,j is the importance weight of node j at time t; δ is the global scaling factor; σ(·) is the activation function.

[0048] Weight neighborhood features according to their spatiotemporal importance weights:

[0049]

[0050] Among them, h i is the weighted neighborhood aggregation feature of node i, is the feature weight matrix.

[0051] Construct the cost matrix:

[0052]

[0053] Among them, C ij represents the transmission cost from source domain node i to target domain node j; is the characteristic distance from source domain node i to target domain node j, which is used to measure the gene expression difference between nodes in the present invention; is the feature vector of source domain node i after standardization and dimension reduction; is the gene expression vector feature of the target domain node j; ||·||2 means finding the L2 norm; is the distribution alignment item from source domain node i to target domain node j, which is used to measure the KL divergence difference of neighborhood distribution in this invention; KL(p i ||q j ) represents the neighborhood distribution p of source domain node i i The neighborhood distribution q of the target domain node j j The KL divergence between .

[0054] Introducing the probability polarization optimal transmission constraint, setting the intra-class / inter-class polarization threshold, and forcing the transmission strength between nodes of the same class to be no less than the threshold λ inter , limiting the transmission intensity between heterogeneous nodes to no more than the threshold λ intra :

[0055]

[0056] Among them, N source Indicates the total number of source domain nodes; N target Indicates the total number of nodes in the target domain.

[0057] Build regularization terms to ensure that the transfer plan is dense within the class and sparse between classes to avoid cross-domain negative transfer:

[0058]

[0059] Among them, Π ij is the element in the i-th row and j-th column of the transmission plan matrix Π, which represents the transmission amount from source node i to target node j.

[0060] Use the Sinkhorn algorithm to solve the optimal transmission plan with entropy regularization:

[0061]

[0062] Where π represents the transmission plan matrix; represents the transpose of the transmission plan matrix π; C ij represents the transmission cost from source domain node i to target domain node j; ∈1 is the regularization coefficient; KL(Π1||μ) represents the KL divergence between Π1 and μ, Π1 is the product of the transmission plan matrix Π and the all-one vector, represents the total output transmission of all nodes in the source domain, and μ is the expected output distribution of the source domain nodes; Express request The KL divergence between ν and is the transpose of Π1, which represents the total input transmission of the target domain nodes; ν is the expected input distribution of the target domain nodes.

[0063] Use Barycentric mapping to achieve cross-domain neighborhood information fusion:

[0064]

[0065] in, is the reconstructed feature of the target domain node j; is the element of the optimized transmission plan matrix, which represents the normalized transmission amount from source node i to target node j; represents the normalized transmission amount from source domain node k to target domain node j.

[0066] like Figure 2 Figure (a) shows the spatial distribution of the features of the human breast cancer dataset reconstructed using this method, clustered into 11 categories based on kmeans. The cluster regions are continuously distributed across the tissue sections, demonstrating that dynamic neighborhood adjustment effectively solves the clustering fragmentation problem caused by fixed neighborhoods in traditional methods. Figure (b) shows the tSNE distribution results, showing clear boundaries between classes and tight clustering within classes, indicating that the fused features can better distinguish cell types or functional states and are more discriminative.

[0067] In addition, Table 1 gives the ablation experiment of dynamic graph construction, where Baseline 1 has a fixed neighborhood radius and Baseline 2 has no node importance weight.

[0068] Table 1 Experimental results of dynamic graph ablation module

[0069] Methods\Indicators ARI NMI OT distance Neighborhood radius variance Baseline1 0.65 0.58 1.20 0.0 Baseline2 0.71 0.63 1.05 0.12 Method of the present invention 0.82 1.05 0.83 0.35

[0070] Table 2 shows the ablation experiment of optimal transmission fusion, where Baseline3 does not have optimal transmission fusion and Baseline4 has no polarization constraint.

[0071] Table 2 Ablation experiment results of optimal transmission module

[0072] Methods\Indicators Cross-domain feature correlation Intra-class closeness Inter-class separation Polarization constraint validity Baseline3 0.68 2.15 3.02 1.0 Baseline4 0.73 1.98 3.25 2.3 Method of the present invention 0.89 1.42 4.17 5.8

[0073] It can be seen from the above ablation experimental data that the method of the present application has made significant progress compared with the existing technology.

[0074] While embodiments of the present invention have been shown and described, it will be appreciated by those skilled in the art that various changes, modifications, substitutions, and variations may be made to these embodiments without departing from the principles and spirit of the invention, and that the scope of the invention is defined by the appended claims and their equivalents.

Claims

1. A spatial transcriptome data clustering method, characterized in that: The following steps are involved: Obtain raw data and perform preprocessing; According to the spatial density and gene expression similarity, a kernel density estimation function is constructed, and the density threshold τ is determined based on the kernel density estimation function, and the combined benchmark neighborhood radius r base Dynamically adjust the neighborhood range of each node; Graph wavelet transform is used to calculate low-frequency coefficients and high-frequency coefficients, and the node importance weight is obtained by combining the spatiotemporal activity score; The neighborhood features are weighted according to the node importance weights, and a cost matrix is constructed by combining feature distance and distribution difference to quantify the transmission cost from the source domain to the target domain. Build regularization terms to ensure that the transfer plan is dense within a class and sparse between classes to avoid cross-domain negative transfer; The Sinkhorn algorithm is used to solve the regularized transmission plan, and the cross-domain information is fused through Barycentric mapping to generate the reconstructed features of the target domain. Clustering is performed using the reconstructed features to obtain clustering results.

2. The method according to claim 1, characterized in that When the kernel density estimation function of a node is less than the density threshold, the neighborhood radius of the node is set to N×r base , N>1; otherwise the neighborhood radius of the node is set to M×r base , 1>M>0, reference neighborhood radius r base is the median of the distance set between any two nodes, and the density threshold τ is the median of the kernel density estimation function set of the nodes.

3. The method according to claim 1 or 2, characterized in that Based on the joint measurement of spatial coordinates and gene expression similarity, a kernel density estimation function is constructed, which is expressed as: Among them, ρ(x i ) is the kernel density estimation function; h is the bandwidth parameter; n is the total number of nodes; d is the dimension of the spatial coordinate; x i Indicates the two-dimensional spatial position information of the i-th node; g i represents the gene expression feature vector of node i; || || represents the norm; sim(g i ,g j ) is the cosine similarity of gene expression, g i is the gene expression feature vector of the i-th node.

4. The method according to claim 1, wherein The calculation of node importance weight includes: b t,j =δ·σ(τ t,j ) t t,j =γ·Δa t,j +(1-c)·b t,j Among them, β t,j is the importance weight of node j at time t; δ is the global scaling factor; σ(·) is the activation function; τ t,j is the spatiotemporal activity score of node j at time t; γ is the mixing coefficient; Δa t,j is the rate of change of the low-frequency coefficient of node j at time t; b t,j is the high frequency coefficient of node j at time t.

5. The method according to claim 1, wherein The cost of transferring data from the source domain to the target domain includes: Among them, C ij represents the transmission cost from source domain node i to target domain node j; represents the characteristic distance from source domain node i to target domain node j; is the feature vector of source domain node i after standardization and dimension reduction; is the gene expression vector feature of the target domain node j; ||·||2 means finding the L2 norm; represents the distribution alignment item from source domain node i to target domain node j; KL(p i ||q j ) represents the neighborhood distribution p of source domain node i i The neighborhood distribution q of the target domain node j j The KL divergence between .

6. The method according to claim 1, characterized in that The regularization term is expressed as: in, is the regularization term of the transmission plan matrix; ij is the element in the i-th row and j-th column of the transmission plan matrix, representing the transmission amount from source domain node i to target domain node j; intra is the minimum threshold of transmission strength between similar nodes; inter is the maximum threshold of transmission strength between heterogeneous nodes.

7. The method according to claim 6, characterized in that The minimum threshold λ of transmission strength between similar nodes inter Expressed as: Among them, N source Indicates the total number of source domain nodes.

8. The method according to claim 6, characterized in that The maximum threshold λ of transmission strength between heterogeneous nodes intra Expressed as: Among them, N target Indicates the total number of nodes in the target domain.

9. The method according to claim 1, characterized in that Solving the regularized transport plan using the Sinkhorn algorithm involves: Where π represents the transmission plan matrix; represents the transpose of the transmission plan matrix Π; ij is the element in the i-th row and j-th column of the transmission plan matrix, which represents the transmission amount from source node i to target node j; C ij represents the transmission cost from source domain node i to target domain node j; ∈1 is the regularization coefficient; Π1 is the product of the transmission plan matrix Π and the all-one vector, representing the total output transmission volume of all nodes in the source domain; is the transpose of Π1, representing the total input transmission of the target domain nodes; μ is the expected output distribution of the source domain nodes; ν is the expected input distribution of the target domain nodes.

10. The method according to claim 1, characterized in that Using Barycentric mapping to achieve cross-domain neighborhood information fusion includes: in, is the reconstructed feature of the target domain node j; N source Indicates the total number of source domain nodes; is the element of the optimized transmission plan matrix, which represents the normalized transmission amount from source node i to target node j; is the gene expression vector feature of source domain node i; represents the normalized transmission amount from source domain node k to target domain node j.

Images