A spatial transcriptome data clustering method

Abstract

The present application relates to the field of spatial transcriptome data analysis, and particularly relates to a spatial transcriptome data clustering method, comprising obtaining original data and preprocessing; constructing a kernel density estimation function according to spatial density and gene expression similarity, dynamically adjusting the neighborhood range of each node based on the kernel density estimation; using graph wavelet transform to calculate low-frequency coefficients and high-frequency coefficients, combining the spatiotemporal activity score to obtain the node importance weight; weighting the neighborhood features according to the node importance weight, combining the feature distance and distribution difference to construct a cost matrix, and quantifying the transmission cost from the source domain to the target domain; using Sinkhorn algorithm to solve the regularized transmission plan, fusing the cross-domain information through Barycentric mapping, and generating the reconstructed features of the target domain; clustering using the reconstructed features to obtain the clustering results. The present application solves the problem of insufficient utilization of spatial information in traditional clustering methods, and provides a more accurate and adaptive solution for single-cell spatial omics analysis.

Description

Technical Field

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

[0002] Spatial transcriptome clustering divides tissues into groups with similar molecular characteristics based on the similarity of gene expression profiles, revealing the biological functions, cell type distribution, or state heterogeneity of different tissue regions. Non-spatial clustering methods use traditional clustering techniques, such as K-means and Louvain algorithms. These methods are limited to a small number of spots, and the clustering results may be discontinuous in tissue slices. While dynamic graph structure learning can adjust neighborhood ranges in recent spatial clustering methods, it lacks an optimization mechanism for information transfer weights in heterogeneous density regions, resulting in significant information loss in low-density areas. Optimal transport theory, although capable of modeling inter-node relationships, relies on a fixed spatial topology when applied alone, limiting its adaptability to spatially heterogeneous density distributions. Current technologies have not yet solved the problem of co-modeling both spatial and non-spatial clustering methods in spatial transcriptome data, leading to information redundancy or loss in heterogeneous regions. Summary of the Invention

[0003] To address the collaborative modeling problem of dynamic graph structure learning and optimal transport theory, this invention proposes a spatial transcriptome data clustering method, comprising the following steps:

[0004] Acquire raw data and perform preprocessing;

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

[0006] Low-frequency and high-frequency coefficients are calculated using graph wavelet transform, and node importance weights are obtained by combining spatiotemporal activity scores; wherein, the node mentioned in this invention refers to the spot in spatial transcriptome data, and each node includes two-dimensional spatial location information and a corresponding gene expression feature vector;

[0007] The neighborhood features are weighted according to the importance of the nodes, and a cost matrix is ​​constructed by combining the 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 transport plan, and cross-domain information is fused through Barycentric mapping to generate the reconstructed features of the target domain.

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

[0010] This invention achieves adaptive adjustment of the neighborhood range through kernel density estimation, and significantly improves the alignment accuracy of spatial and gene expression information by introducing probabilistic polarization optimal transport and Barycentric Wasserstein paths. Through these technological innovations, this invention solves the problem of insufficient spatial information utilization in traditional clustering methods, providing a more accurate and adaptive solution for single-cell spatial omics analysis. Attached Figure Description

[0011] Figure 1 This is a flowchart of a spatial transcriptome data clustering method according to the present invention;

[0012] Figure 2 This is a demonstration of the clustering results of one case of data from a human breast cancer dataset, based on the present invention. Detailed Implementation

[0013] The technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0014] This invention proposes a spatial transcriptome data clustering method, such as... Figure 1 This includes the following steps:

[0015] Acquire raw data and perform preprocessing;

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

[0017] Low-frequency and high-frequency coefficients are calculated using graph wavelet transform, and node importance weights are obtained by combining spatiotemporal activity scores.

[0018] The neighborhood features are weighted according to the importance of the nodes, and a cost matrix is ​​constructed by combining the 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 transport plan, and cross-domain information is fused through Barycentric mapping to generate the reconstructed features of the target domain.

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

[0021] To address the collaborative modeling challenge of dynamic graph structure learning and optimal transport theory, this invention proposes a neighborhood information fusion mechanism that combines dynamic graph structure learning and optimal transport theory. This mechanism generates a dynamic neighborhood range through kernel density estimation, expands the neighborhood range in low-density regions, and shrinks the connection radius in high-density regions. Simultaneously, optimal transport theory is introduced to quantify the topological changes of neighborhood expansion / shrinkage into a differentiable transport cost function, achieving joint optimization of graph structure and information transmission weights. The system demonstrates significant improvements in redundant connection suppression, information integrity preservation, and error control, verifying that the synergistic effect of dynamic graph learning and optimal transport theory can overcome the inherent limitations of single technologies and achieve more robust spatial topology modeling.

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

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

[0024] The second step is to 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 combining them with spatiotemporal activity scores to obtain the node importance weight.

[0025] Step 3: Optimal transmission-driven neighborhood information fusion, including: weighting neighborhood features according to the spatiotemporal importance of nodes, constructing a cost matrix by combining 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; using the Sinkhorn algorithm to solve the regularized transmission plan, fusing cross-domain information through Barycentric mapping, and generating reconstructed features of 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] A kernel density estimation function is constructed based on a joint measure of spatial coordinates and gene expression similarity:

[0029]

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

[0031] Furthermore, high / low density regions are dynamically divided using a density threshold τ. A node is considered a low-density region when its kernel density estimation function is less than the density threshold, and the node's neighborhood radius is set to N×r. base If N > 1, then it is a high-density region, and the neighborhood radius of the node is set to M × r. base 1 > M > 0, reference neighborhood radius r base The distance between any two nodes is the median of the distance set. In this invention, a node is a spot in the spatial transcriptome data. Each node includes two-dimensional spatial location information and a corresponding gene expression feature vector, with a reference neighborhood radius r. base The calculations include:

[0032] Based on the two-dimensional spatial location information of the nodes {x1,x2,...,x} n Find the Euclidean distance between any two nodes, and the two-dimensional spatial location information x of the nodes. i Two-dimensional spatial location 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, and n is the number of nodes;

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

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

[0035] Based on the set of two-dimensional spatial location information of the nodes {x1,x2,...,x} n} and the set of gene expression feature vectors of the nodes {g1, g2, ..., g n}, calculate the kernel density estimation function for each node;

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

[0037] As a preferred implementation, the neighborhood radius is extended to r in the low-density region.low =2r base The high-density region shrinks to r high =0.5r base This enables adaptive adjustment of the neighborhood radius.

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

[0039]

[0040] Where, Δa t,j Let be the rate of change of the low-frequency coefficient, and ∈ be the numerical stability constant, taken as 10. -7 .

[0041] Extracting high-frequency coefficients from feature vectors:

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

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

[0044] The time importance weights of computation nodes:

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

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

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

[0048] Neighborhood features are weighted according to their spatiotemporal importance:

[0049]

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

[0051] Construct the cost matrix:

[0052]

[0053] Among them, C ij This represents the transmission cost from source domain node i to target domain node j; The feature distance from source domain node i to target domain node j is used in this invention to measure the difference in gene expression between nodes. is the feature vector of source domain node i after standardization and dimensionality reduction; Let ||j| represent the gene expression vector features of node j in the target domain; ||·||2 represents calculating the L2 norm; KL(p) represents the distribution alignment term from source domain node i to target domain node j, used in this invention to measure the KL divergence difference in neighborhood distributions; 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 them.

[0054] Introducing probabilistic polarization optimal transmission constraints, we set intra-class / inter-class polarization thresholds to force the transmission strength between nodes of the same class to be no less than the threshold λ. inter Limit the transmission strength between heterogeneous nodes to no more than the threshold λ. intra :

[0055]

[0056] Where, N source N represents the total number of nodes in the source domain. target This indicates the total number of nodes in the target domain.

[0057] Construct regularization terms to ensure that transport plans are dense within classes and sparse between classes, avoiding negative cross-domain migration:

[0058]

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

[0060] The Sinkhorn algorithm is used to solve for 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 denoted as 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 μ, where Π1 is the product of the transmission plan matrix Π and the all-1 vector, denoted as the total output transmission amount of all nodes in the source domain, and μ is the expected output distribution of the source domain nodes. Expressing the request KL divergence between ν and ν ν is the transpose of Π1, representing the total input transmission volume of the target domain node; ν is the expected input distribution of the target domain node.

[0063] Using Barycentric mapping, cross-domain neighborhood information fusion is achieved:

[0064]

[0065] in, The reconstructed features of the target domain node j; For the elements of the optimized transmission plan matrix, represent the normalized transmission amount from source node i to destination node j; This represents the normalized transmission amount from source domain node k to destination domain node j.

[0066] like Figure 2 Figure (a) shows the spatial distribution of the features of the human breast cancer dataset reconstructed by this invention, clustered into 11 categories using k-means. The clustered regions are continuously distributed on the tissue slices, indicating that dynamic neighborhood adjustment effectively solves the clustering fragmentation problem caused by fixed neighborhoods in traditional methods. Figure (b) shows the tSNE distribution results, which show clear inter-class boundaries and tight intra-class clustering, indicating that the fused features can better distinguish cell types or functional states and are more discriminative.

[0067] In addition, Table 1 presents the ablation experiments constructed from the dynamic graph, where Baseline 1 has a fixed neighborhood radius and Baseline 2 has no node importance weights.

[0068] Table 1. Experimental results of the 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 presents the ablation experiments for optimal transport fusion, where Baseline 3 has no optimal transport fusion and Baseline 4 has no polarization constraints.

[0071] Table 2 Ablation Experiment Results of Optimal Transmission Module

[0072] Methods / Indicators Cross-domain feature correlation Intraclass tightness Inter-class separation Polarization constraint effectiveness 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] As can be seen from the above ablation experiment data, the method of this application has significant progress compared with the prior art.

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

Claims

1. A spatial transcriptome data clustering method, characterized in that, Includes the following steps: Acquire raw data and perform preprocessing; A kernel density estimation function is constructed based on spatial density and gene expression similarity. A density threshold τ is determined based on the kernel density estimation function, and a benchmark neighborhood radius r is also considered. base Dynamically adjust the neighborhood range of each node; Low-frequency and high-frequency coefficients are calculated using graph wavelet transform, and node importance weights are obtained by combining spatiotemporal activity scores. The neighborhood features are weighted according to the importance of the nodes, and a cost matrix is ​​constructed by combining the feature distance and distribution difference to quantify the transmission cost from the source domain to the target domain. Construct regularization terms to ensure that the transport plan is dense within classes and sparse between classes, avoiding negative cross-domain migration; The Sinkhorn algorithm is used to solve the regularized transport plan, and cross-domain information is fused through Barycentric mapping to generate the reconstructed features of the target domain. Clustering is performed using 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 If N > 1, then the neighborhood radius of the node is set to M × r. base 1 > M > 0, reference neighborhood radius r base Let τ be the median of the distance set between any two nodes, and let τ be 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 a joint measure of spatial coordinates and gene expression similarity, a kernel density estimation function is constructed, expressed as: Wherein, ρ(x) i ) represents the kernel density estimation function; h is the bandwidth parameter; n is the total number of nodes; d is the dimension of the spatial coordinates; x i Represents the two-dimensional spatial location information of the i-th node; g i Represents the gene expression feature vector of node i; || || denotes the norm; sim(g i ,g j ) represents the cosine similarity of gene expression, g i Let be the gene expression feature vector of the i-th node.

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

5. The method according to claim 1, characterized in that, The transmission cost from the source domain to the destination domain includes: Among them, C ij This represents the transmission cost from source domain node i to target domain node j; This represents the feature distance from source domain node i to target domain node j; is the feature vector of source domain node i after standardization and dimensionality reduction; Let ||j| represent the gene expression vector features of node j in the target domain; ||·||2 represents calculating the L2 norm; KL(p) represents the distribution alignment term from source domain node i to target domain node j; 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 them.

6. The method according to claim 1, characterized in that, The regularization term is represented as: in, For regularization terms of the transmission plan matrix; Π ij Let λ be the element in the i-th row and j-th column of the transmission plan matrix, representing the transmission amount from source node i to destination node j; intra λ represents the minimum threshold for transmission strength between nodes of the same type. inter This represents the highest threshold for transmission strength between heterogeneous nodes.

7. The method according to claim 6, characterized in that, Minimum threshold λ for transmission strength between similar nodes inter Represented as: Where, N source This represents the total number of nodes in the source domain.

8. The method according to claim 6, characterized in that, The highest threshold λ for transmission strength between heterogeneous nodes intra Represented as: Where, N target This indicates the total number of nodes in the target domain.

9. The method according to claim 1, characterized in that, Solving regularized transport plans using the Sinkhorn algorithm includes: Where Π represents the transmission plan matrix; Represents the transpose of the transmission plan matrix Π; Π ij Let C be the element in the i-th row and j-th column of the transmission plan matrix, representing the transmission amount from source node i to destination node j; ij denoted as 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-1 vector, representing the total output transmission amount of all nodes in the source domain; ν is the transpose of Π1, representing the total input transmission of the target domain node; μ is the expected output distribution of the source domain node; and ν is the expected input distribution of the target domain node.

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

Images