A spatial cell-cell communication inference method based on regulation sub-guided neural stochastic differential equation

Abstract

The present application relates to a kind of space intercellular communication inference method based on the regulation of sub-guide neural stochastic differential equation. It solves the problem that the existing technology is strongly dependent on static space neighborhood, and the evidence of downstream regulation function is not fully utilized, and expression-driven false positive is prone to occur. It includes learning dynamic diffusion adjacency matrix of tissue-dependent through neural stochastic differential equation, and combining the transcription factor regulon activity inferred by SCENIC to weight the functional support of candidate ligand-receptor interaction, so as to improve the accuracy, spatial interpretability and biological credibility of spatial intercellular communication event identification. The present application has the advantages that it can describe the tissue-dependent communication accessibility beyond fixed Euclidean distance or static k-neighbor graph, make high-ranking communication events have stronger functional interpretability, and facilitate verification combined with tissue structure and biological background.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics technology, specifically to a spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations. Background Technology

[0002] Cell-cell communication (CCC) is fundamental to tissue structure maintenance, immune regulation, developmental pattern formation, and disease-related tissue remodeling. In most biological systems, extracellular ligand-receptor interactions are typically associated with intracellular signal transduction, transcription factor activation, and changes in downstream target gene expression. Therefore, elucidating CCC requires not only identifying potential ligand-receptor pairs but also determining whether the interaction occurs at an appropriate tissue location and further assessing whether it is accompanied by functional transcriptional regulatory responses in the receiving cells.

[0003] With the development of single-cell transcriptome sequencing and spatial transcriptomics technologies, researchers have proposed various computational methods for inferring intercellular communication. Traditional methods, such as CellPhoneDB and CellChat, mainly rely on ligand-receptor databases and cell population expression co-occurrence patterns to identify candidate communication events. Some recent spatial methods further incorporate spatial distance, neighborhood graphs, tissue topology, or attention mechanisms to improve the spatial plausibility of predicted interactions. These methods provide important tools for communication analysis in complex tissue microenvironments.

[0004] However, existing technologies still have the following shortcomings. First, many spatial intercellular communication methods still use fixed Euclidean distance, fixed k-nearest neighbor graphs, or pre-defined spatial neighborhoods to describe communication accessibility, which makes it difficult to fully characterize the effects of cell density, tissue structure, local permeability, chemotactic cues, and random diffusion perturbations on signal propagation in real tissues. Second, spatial proximity and ligand-receptor co-expression alone cannot prove that communication events have downstream functional effects. Without the support of transcription factor regulatory activity or target gene programming, candidate interactions may only reflect expression abundance or database co-occurrence. Third, ligand-receptor pairs with strong expression support but spatial inaccessibility or lack of functional evidence are easily misranked, leading to spatial false positives and expression-driven false positives.

[0005] Therefore, it is necessary to provide a new method for inferring spatial intercellular communication, which can further learn tissue-dependent dynamic spatial reachability based on ligand-receptor expression, and introduce downstream transcription factor regulator activity as evidence for functional verification, thereby prioritizing the screening of highly reliable communication events that are both spatially reachable and regulatoryally supported. Summary of the Invention

[0006] The purpose of this invention is to address the technical problems of existing spatial intercellular communication inference methods, such as strong dependence on static spatial neighborhoods, insufficient utilization of evidence of downstream regulatory functions, and susceptibility to expression-driven false positives. This invention provides a spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations. This method learns a tissue-dependent dynamic diffusion adjacency matrix through neural stochastic differential equations and combines the transcription factor regulon activity inferred by SCENIC to perform functional support weighting on candidate ligand-receptor interactions, thereby improving the accuracy, spatial interpretability, and biological credibility of spatial intercellular communication event identification.

[0007] To achieve the above objectives, the present invention employs the following technical solution: a spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations, comprising the following steps:

[0008] S1. Data Input and Preprocessing: Input spatial transcriptome expression matrix, spatial coordinates, candidate ligand-receptor database, and prior knowledge of transcription factor regulation;

[0009] S2. Construction of local ligand-receptor activity: A k-nearest neighbor graph is constructed for spatial coordinates, and local neighborhood expression is integrated according to the distance decay kernel; for candidate ligand-receptor pairs, the product of the ligand expression after neighborhood integration and the receptor expression is used as the local ligand-receptor activity to obtain the spatial activity matrix;

[0010] S3 and SCENIC regulator function verification: Transcription factors associated with candidate ligand-receptor pairs were extracted from prior transcription factors. The regulatory relationship between transcription factors and target genes was inferred using GRNBoost2. A co-expression module was constructed. Motif-supported regulons were obtained through cisTargetmotif pruning. The activity of each regulon in cells or spots was calculated using AUCell.

[0011] S4. Dynamic diffusion modeling of neural stochastic differential equations: using standardized spatial coordinates as the initial state, the continuous-time diffusion process of organization dependence is learned through drift network and diffusion network;

[0012] S5. Calculate the dynamic adjacency matrix: Transform the diffusion endpoint obtained from the neural stochastic differential equation back to the original spatial scale, and calculate the edge weights based on the distance between the diffused transmitter position and the original receiver position, thereby obtaining the dynamic spatial communication topology and forming a dynamic adjacency matrix.

[0013] S6. Calculate the spatial specificity score: The ligand-receptor activity distribution is propagated on the dynamic adjacency matrix, the consistency of the activity distribution before and after propagation is compared, and the penalty term is calculated by randomly permuting the adjacency background to obtain the ligand-receptor spatial specificity score corrected by random spatial background.

[0014] S7. Integration of transcription factor functional support: For candidate ligand-receptor pairs with downstream transcription factor annotations and corresponding regulons that can be recognized by SCENIC, their regulon activity is included as a functional enhancer in the final score.

[0015] S8. Output Analysis Results: Based on the final score, high-confidence ligand-receptor pairs are selected, and the communication network between the sending and receiving cell types is further inferred by combining cell type annotations, while generating a spatially resolved communication field map.

[0016] In step S1, let the spatial coordinates of the i-th cell or spot be c. i The expression level of the g-th gene is X. ig The initial spatial distance and initial adjacency weight are defined as follows:

[0017] ;

[0018] Where, d ij N represents the Euclidean distance between cell or spot i and j. k(i) Let A represent the set of k nearest neighbors of i. 0 ij A represents the initial spatial adjacency weight, and σ0 represents the distance decay scale parameter; 0 ij It is retained only if j belongs to the k nearest neighbor set of i, otherwise it is set to 0; σ0 is the scale parameter that controls the initial spatial distance decay.

[0019] In step S2, the expression levels of each gene are integrated into their local spatial neighborhood to obtain a local microenvironment-sensing expression matrix:

[0020] ;

[0021] Among them, X ig This represents the gene expression level after integration within a local spatial neighborhood; for the k-th candidate ligand-receptor pair (l k ,r k Its local ligand-receptor communication activity is defined as:

[0022] ;

[0023] Among them, Z ik This indicates the local communication activity of the k-th ligand-receptor pair on the cell or spoti.

[0024] I (·) This is an indicator function.

[0025] In step S3, let Regulon tAUCell represents a regulator centered on transcription factor t and supported by motifs. i (Regulon t The fraction () represents the activity of the regulator in the i-th cell or spot. For candidate ligand-receptor pair k, its regulator support fraction is defined as:

[0026] ;

[0027] Among them, TF k Let TF represent the downstream regulator support strength of the k-th candidate ligand-receptor pair; if candidate ligand-receptor pair k has no matching downstream transcription factor, or the corresponding transcription factor does not form an effective regimen, then let TF k =0.

[0028] In step S4, the neural stochastic differential equations, with standardized spatial coordinates as the initial state, learn the continuous-time diffusion process:

[0029] ;

[0030] Among them, y t This represents the diffusion state initialized by spatial coordinates;

[0031] f θ For drift networks;

[0032] g φ For diffusion networks;

[0033] W t For a Wiener stochastic process, for the k nearest neighbor edge i→j, the self-supervised objective of ligand-receptor affinity is:

[0034] ;

[0035] ;

[0036] q ij Q represents the edge affinity self-supervised target obtained from ligand-receptor expression. 0.99 L represents the 99th percentile normalization factor. SDE Represents the training loss of the neural stochastic differential equation; by minimizing L SDE This ensures that the spatial reachability after diffusion remains consistent with the ligand-receptor expression affinity.

[0037] In step S5, the diffusion endpoint y T (i) Denormalize to the original scale coordinates z i And calculate the dynamic adjacency matrix:

[0038] ;

[0039] z i A represents the transmitter location mapped back to the original coordinate scale after the diffusion endpoint is reached. dyn ij Represents the communication reachability weights in the dynamic adjacency matrix; the dynamic adjacency matrix preserves the local spatial topological skeleton, but its edge weights are determined by the diffusion reachability learned from neural stochastic differential equations.

[0040] In step S6, for the k-th ligand-receptor pair, let G k This represents its spatial activity vector across all cells or spots. The activity vector after dynamic topological propagation is:

[0041] ;

[0042] The observation spatial consistency score and the random background score are defined as follows:

[0043] ;

[0044] ;

[0045] Among them, G k Let A represent the spatial activity vector of the k-th ligand-receptor pair. rand b Let b be the background adjacency matrix obtained by the b-th random permutation, and let the ligand-receptor spatial specificity score be defined as:

[0046] ;

[0047] Among them, LRSS k denoted as the ligand-receptor spatial specificity score after random background correction, and μ represents the random background penalty coefficient.

[0048] In step S7, the space-specific score and regulator support score are integrated into the final communication score:

[0049] ;

[0050] Among them, Score k λ represents the final communication score, and λ represents the transcription factor functional support weight.

[0051] In step S8, let a be the sending cell type and b be the receiving cell type. Then, the communication strength at the cell type level can be obtained by aggregation of high-confidence candidate ligand-receptor pairs:

[0052] ;

[0053] Among them, Spec a (lk ) indicates ligand l k Expression specificity in cell type a

[0054] Spec b (r k ) represents receptor r k Expression specificity in receiving cell type b

[0055] Coloc ab (k) represents the colocalization support of ligand-receptor pair k in the sending-receiving cell type spatial neighborhood.

[0056] In the aforementioned spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations, the spatial transcriptome expression matrix is ​​X∈R (n×p) This indicates that the set of spatial coordinates is represented by C={c i ∈R d} indicates that the candidate ligand-receptor and its downstream transcription factor annotation set is represented by P={(l k ,r k ,T k )}express.

[0057] Compared with the prior art, the advantages of the present invention are as follows:

[0058] 1. This invention learns a dynamic spatial diffusion adjacency matrix through neural stochastic differential equations, which can describe the reachability of organization-dependent communication beyond fixed Euclidean distance or static k-nearest neighbor graphs.

[0059] 2. This invention introduces SCENICregulon activity as evidence of downstream transcriptional regulation, making high-ranking communication events more functionally explanatory.

[0060] 3. This invention reduces false positives caused by expression abundance and random topology through random spatial background correction;

[0061] 4. This invention can simultaneously output high-confidence ligand-receptor pairs, spatial communication fields, cell type communication networks, and transcription factor support signals, which facilitates verification in conjunction with tissue structure and biological background. Attached Figure Description

[0062] Figure 1 This is a flowchart illustrating the overall workflow of the present invention; wherein:

[0063] Figure a shows the inputs required for this invention, including spatial transcriptome data, spatial location, ligand-receptor database, and transcription factor database;

[0064] Figure b illustrates the SCENIC-based regulator validation process, including gene regulatory network inference, co-expression module construction, topic trimming, regulatory factor activity calculation, and ligand-receptor extraction.

[0065] Figure c illustrates the use of neural stochastic differential equations to learn dynamic diffusion trajectories from initial spatial cells and form a dynamic adjacency matrix based on diffusion effects;

[0066] Figure d shows spatial specificity scoring, random background correction, ligand-receptor sequencing, transcription factor-supported signaling pathway output, and cell type communication inference.

[0067] Figure 2 This is a schematic diagram illustrating the benchmark performance evaluation results of the present invention on nine simulated datasets; wherein:

[0068] Figure a shows the performance of each method on average AUC, SpatialFPR, and DecoyFPR;

[0069] Figure b shows the AUC distribution across the datasets;

[0070] Figure c shows the overall score AS;

[0071] Figure d shows the DecoyFPR heatmaps on different datasets;

[0072] Figure e illustrates the trade-off between AUC and Spatial FPR;

[0073] Figure 3 This is a schematic diagram illustrating the cross-method consistency, transcription factor support, and representative ligand-receptor spatial field validation results of the present invention in real tissue data; wherein:

[0074] Figure a shows the number of global top 300 ligand-receptor pairs shared by the present invention and the baseline method;

[0075] Figure b shows the joint distribution of communication score and TFregulon activity;

[0076] Figures c and d illustrate the ligand expression, receptor expression, interaction field, and continuous interaction landscape of TGFBR2::TGFBR1;

[0077] Figures e and f show the spatial field verification results of SRGN::CD44;

[0078] Figures g and h show weak evidence of spatial colocalization of baseline-specific candidates TNFSF13B::TNFRSF17 and CD80::CD28. Detailed Implementation

[0079] The present invention will now be described in further detail with reference to the accompanying drawings and specific embodiments.

[0080] like Figure 1-3 As shown, a spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations includes the following steps:

[0081] S1. Data Input and Preprocessing: Input spatial transcriptome expression matrix, spatial coordinates, candidate ligand-receptor database, and transcription factor regulatory priors; the expression matrix provides the expression levels of ligands, receptors, and transcription factors; spatial coordinates describe the tissue location of cells or spots; the candidate ligand-receptor database defines the search space; the transcription factor regulatory priors are used to establish the correspondence between candidate communication events and downstream regulatory procedures.

[0082] S2. Construction of local ligand-receptor activity: A k-nearest neighbor graph is constructed for spatial coordinates, and local neighborhood expression is integrated according to the distance decay kernel; for candidate ligand-receptor pairs, the product of the ligand expression and receptor expression after neighborhood integration is used as the local ligand-receptor activity to obtain the spatial activity matrix; in order to reduce the influence of isolated expression noise on communication inference, this step makes communication activity not only determined by the expression of isolated cells or spots, but also affected by the local spatial microenvironment.

[0083] S3 and SCENIC regulator function verification: Transcription factors associated with candidate ligand-receptor pairs were extracted from prior transcription factors. The regulatory relationship between transcription factors and target genes was inferred using GRNBoost2. A co-expression module was constructed, and motif-supported regimenls were obtained through cisTargetmotif pruning. The activity of each regimenlon in cells or spots was calculated using AUCell. The regimenlon activity serves as downstream functional support for candidate ligand-receptor pairs.

[0084] S4. Dynamic diffusion modeling of neural stochastic differential equations: using standardized spatial coordinates as the initial state, the continuous-time diffusion process of organization dependence is learned through drift network and diffusion network;

[0085] The training objective consists of ligand-receptor affinity on the k-nearest neighbor edges to ensure that the learned diffusion reachability is consistent with the candidate communication relationships supported by expression.

[0086] Among them, the drift network is used to learn the deterministic direction of change of spatial reachability, and the diffusion network is used to learn the intensity of local random diffusion.

[0087] During training, the ligand-receptor affinity on the initial k nearest neighbor edges is used as the self-supervised target. The affinity is obtained by the aggregation result of the sending end ligand expression and the receiving end receptor expression on the candidate ligand-receptor set, and is processed by log1p transformation and quantile normalization.

[0088] S5. Calculate the dynamic adjacency matrix: Transform the diffusion endpoint obtained from the neural stochastic differential equation back to the original spatial scale, and calculate the edge weights based on the distance between the diffused transmitter position and the original receiver position, thereby obtaining the dynamic spatial communication topology and forming a dynamic adjacency matrix.

[0089] The dynamic adjacency matrix preserves the initial local topological skeleton, but the edge weights are determined by the learned diffusion reachability.

[0090] S6. Calculate the spatial specificity score: The ligand-receptor activity distribution is propagated on the dynamic adjacency matrix, the consistency of the activity distribution before and after propagation is compared, and the penalty term is calculated by randomly permuting the adjacency background to obtain the ligand-receptor spatial specificity score corrected by random spatial background.

[0091] To control for false positives caused by expression abundance, random adjacency, or global co-expression, adjacency relationships are randomly permuted while maintaining edge weight distribution. Random spatial backgrounds are repeatedly generated, and random background similarity is calculated. The ligand-receptor spatial specificity score is defined as the observation spatial consistency score minus the random background penalty term.

[0092] S7. Integration of transcription factor functional support: For candidate ligand-receptor pairs with downstream transcription factor annotations and corresponding regulons that can be recognized by SCENIC, their regulon activity is included as a functional enhancer in the final score.

[0093] The ligand-receptor spatial specificity score was integrated with the SCENIC regiment functional support score. If a candidate ligand-receptor pair has a matching downstream transcription factor, and this transcription factor is recognized by SCENIC as a valid regiment, then the candidate pair receives enhanced functional support; if no matching downstream transcription factor exists or the corresponding regiment is not recognized, no functional enhancement occurs. This yields a final communication score that combines spatial constraints and evidence of downstream regulation. A higher final communication score indicates that the candidate interaction is more likely to simultaneously satisfy spatial accessibility, expression support, and downstream regulatory response.

[0094] S8. Output Analysis Results: Based on the final score, high-confidence ligand-receptor pairs are selected, and the communication network between the sending and receiving cell types is further inferred by combining cell type annotations, while generating a spatially resolved communication field map.

[0095] Candidate ligand-receptor pairs are sorted in descending order based on their final communication scores, and rare local outliers are excluded by combining non-zero expression cell or spot number thresholds. For cell type communication inference, the expression specificity of the ligand in the sending cell type, the expression specificity of the receptor in the receiving cell type, the spatial distribution strength of the ligand-receptor pair, and the neighborhood co-localization support of the sending and receiving cell types are calculated separately to obtain the communication strength at the cell type level.

[0096] In step S1, let the spatial coordinates of the i-th cell or spot be c. i The expression level of the g-th gene is X. ig The initial spatial distance and initial adjacency weight are defined as follows:

[0097] ;

[0098] Where, d ij N represents the Euclidean distance between cell or spot i and j. k(i) Let A represent the set of k nearest neighbors of i. 0 ij A represents the initial spatial adjacency weight, and σ0 represents the distance decay scale parameter; 0 ij It is retained only if j belongs to the k nearest neighbor set of i, otherwise it is set to 0; σ0 is the scale parameter that controls the initial spatial distance decay.

[0099] In step S2, the expression levels of each gene are integrated into their local spatial neighborhood to obtain a local microenvironment-sensing expression matrix:

[0100] ;

[0101] Among them, X ig This represents the gene expression level after integration within a local spatial neighborhood; for the k-th candidate ligand-receptor pair (l k ,r k Its local ligand-receptor communication activity is defined as:

[0102] ;

[0103] Among them, Z ik This indicates the local communication activity of the k-th ligand-receptor pair on the cell or spoti.

[0104] I (·) This is an indicator function used to filter background activity that completely lacks evidence of the original expression.

[0105] In step S3, let Regulon t AUCell represents a regulator centered on transcription factor t and supported by motifs.i (Regulon t The fraction () represents the activity of the regulator in the i-th cell or spot. For candidate ligand-receptor pair k, its regulator support fraction is defined as:

[0106] ;

[0107] Among them, TF k Let TF represent the downstream regulator support strength of the k-th candidate ligand-receptor pair; if candidate ligand-receptor pair k has no matching downstream transcription factor, or the corresponding transcription factor does not form an effective regimen, then let TF k =0.

[0108] In step S4, the neural stochastic differential equations, with standardized spatial coordinates as the initial state, learn the continuous-time diffusion process:

[0109] ;

[0110] Among them, y t This represents the diffusion state initialized by spatial coordinates;

[0111] f θ This is a drift network used to learn deterministic diffusion directions;

[0112] g φ It is a diffusion network used to learn the intensity of random diffusion;

[0113] W t For a Wiener stochastic process, for the k nearest neighbor edge i→j, the self-supervised objective of ligand-receptor affinity is:

[0114] ;

[0115] ;

[0116] q ij Q represents the edge affinity self-supervised target obtained from ligand-receptor expression. 0.99 L represents the 99th percentile normalization factor. SDE Represents the training loss of the neural stochastic differential equation; by minimizing L SDE This ensures that the spatial reachability after diffusion remains consistent with the ligand-receptor expression affinity.

[0117] In step S5, the diffusion endpoint y T (i) Denormalize to the original scale coordinates z i And calculate the dynamic adjacency matrix:

[0118] ;

[0119] z i A represents the transmitter location mapped back to the original coordinate scale after the diffusion endpoint is reached. dyn ij Represents the communication reachability weights in the dynamic adjacency matrix; the dynamic adjacency matrix preserves the local spatial topological skeleton, but its edge weights are determined by the diffusion reachability learned from neural stochastic differential equations.

[0120] In step S6, for the k-th ligand-receptor pair, let G k This represents its spatial activity vector across all cells or spots. The activity vector after dynamic topological propagation is:

[0121] ;

[0122] The observation spatial consistency score and the random background score are defined as follows:

[0123] ;

[0124] ;

[0125] Among them, G k Let A represent the spatial activity vector of the k-th ligand-receptor pair. rand b Let b be the background adjacency matrix obtained by the b-th random permutation, and let the ligand-receptor spatial specificity score be defined as:

[0126] ;

[0127] Among them, LRSS k denoted as the ligand-receptor spatial specificity score after random background correction, and μ represents the random background penalty coefficient.

[0128] In step S7, the space-specific score and regulator support score are integrated into the final communication score:

[0129] ;

[0130] Among them, Score k The final communication score is represented by λ, where λ represents the transcription factor functional support weight. The final score increases accordingly when a candidate ligand-receptor pair exhibits both high spatial specificity and strong downstream regulator activity.

[0131] In step S8, let a be the sending cell type and b be the receiving cell type. Then, the communication strength at the cell type level can be obtained by aggregation of high-confidence candidate ligand-receptor pairs:

[0132] ;

[0133] Among them, Spec a (l k ) indicates ligand l k Expression specificity in cell type a

[0134] Spec b (r k ) represents receptor r k Expression specificity in receiving cell type b

[0135] Coloc ab (k) represents the colocalization support of ligand-receptor pair k in the sending-receiving cell type spatial neighborhood.

[0136] Among them, the spatial transcriptome expression matrix is ​​X∈R (n×p) This indicates that the set of spatial coordinates is represented by C={c i ∈R d} indicates that the candidate ligand-receptor and its downstream transcription factor annotation set is represented by P={(l k ,r k ,T k )}express.

[0137] Performance Validation on Simulated Data: The proposed method was evaluated on nine simulated datasets with biological backgrounds. These datasets included data constructed from human tonsil slide-tags snRNA-seq, mouse brain development Stereo-seq, 10xGenomics breast cancer spatial transcriptome, and breast cancer simulation scenarios with varying signal intensities / sparseness. Evaluation metrics included AUC, Spatial FPR, Decoy FPR, and the overall ranking score AS. Results showed that the proposed method achieved an average AUC of 0.794, an average Spatial FPR of 0.0185, and the highest overall score AS=24, indicating its superior performance in recovering real-world communication events and controlling for unreachable false positives.

[0138] Real-world case validation: In human tonsil spatial transcriptome data, this method was used to evaluate the B_germinal_center related communication axis and compared with methods such as CellChatSpatial, CellPhoneDBv5, HoloNet, DeepTalk, and FastCCC. Results showed that this method identified FDC_LZDZ and T_follicular_helper as the main receiving cell types of the B_germinal_center, consistent with the histological and immunological relationships between B cells and follicular dendritic cells and Tfh cells in the germinal center.

[0139] Spatial field and transcriptional regulation support validation: This method overlaps to some extent with the baseline method's global top 300 ligand-receptor pairs, while also identifying specific candidate pairs with dynamic spatial diffusion and regulon support. For example, TGFBR2::TGFBR1 and SRGN::CD44 form a locally enriched continuous communication field in tissues, while the baseline-specific candidates TNFSF13B::TNFRSF17 and CD80::CD28 show weak evidence of spatial co-localization. Joint analysis of communication scores and TF regulon activity for 819 SCENIC-supported ligand-receptor pairs revealed that the high-score / high-TF-activity quadrant contained 16 candidate pairs, indicating that this method can link high-ranking communication events to downstream transcriptional regulatory states.

[0140] Example 1: As Figure 2 As shown, to verify the effectiveness of the method of the present invention on quantifiable benchmarks, nine simulated datasets with real biological backgrounds were used for evaluation, and compared with methods such as CellChat spatial mode, CellChat non-spatial mode, CellPhoneDB, SingleCellSignalR, scSeqComm, Connectome, NATMI, LogFC, and geometric mean. The nine simulated datasets include human amygdala slide-tags snRNA-seq data, mouse brain development Stereo-seq data, 10xGenomics breast cancer spatial transcriptome data, and simulated breast cancer scenarios under different signal intensities, sparsity, and interference settings.

[0141] The evaluation metrics include AUC, SpatialFPR, and DecoyFPR. AUC measures the ability to recover the order of real communication events; SpatialFPR measures the proportion of spatially unreachable sender-receiver pairs incorrectly ranked in high-ranking results; and DecoyFPR measures the degree to which interfering ligand-receptor pairs with expression support but lacking pre-defined functional communication tags are preferentially ranked. The overall score AS is obtained by summing the rankings of these metrics. The results show that the method of this invention achieves an average AUC of 0.794, an average SpatialFPR of 0.0185, and the highest overall score AS=24, indicating that this invention achieves a good balance between recovering real communication events and controlling spatial false positives.

[0142] Example 2: Figure 3As shown, to verify the biological interpretability of the method of this invention in real tissues, it was applied to human tonsil spatial transcriptome data. First, the overlap between the method of this invention and the baseline method in the global top 300 ligand-receptor pairs was compared to determine whether the method could retain some consensus communication signals. Simultaneously, dynamic spatial diffusion and the regulon function supported the reordering of candidate interactions. The results showed that the method of this invention shared 72, 62, 58, and 57 top ligand-receptor pairs with LIANA-SCA, CellChat non-spatial mode, LIANA-logFC, and CellChat spatial mode, respectively, indicating that the invention can retain a certain degree of consensus signal.

[0143] Furthermore, this invention identified RegSDE-Comm-specific candidate ligand-receptor pairs that were not present in the top 300 results of the baseline method. Among them, TGFBR2::TGFBR1 and SRGN::CD44 exhibited high spatial colocalization and continuous spatial communication fields. The ligand and receptor signals of TGFBR2::TGFBR1 are not uniformly distributed in the tissue; after dynamic spatial accessibility constraints, their interaction fields are concentrated in locally compatible regions. Although SRGN::CD44 involves the widely expressed receptor CD44, this invention can still identify local communication enrichment in its spatially aligned region with SRGN. In contrast, the baseline-specific candidates TNFSF13B::TNFRSF17 and CD80::CD28 had lower spatial colocalization scores and sparser spatial communication fields, indicating that although some baseline candidates have expression or database support, the evidence for spatial colocalization is weak.

[0144] Furthermore, this invention jointly analyzed the communication scores and TF regulon activities of 819 SCENIC-supported ligand-receptor pairs, and categorized candidates according to their 85th percentiles. The quadrant with high communication scores and high TF activity contained 16 ligand-receptor pairs, representing high-confidence candidates supported by both spatial communication signals and downstream regulon function. This result demonstrates that this invention can link high-ranking ligand-receptor events to the downstream transcriptional regulatory state in receiver cells.

[0145] The principle of this embodiment is as follows: This invention learns the dynamic spatial diffusion adjacency matrix through neural stochastic differential equations and introduces SCENICregulon activity into ligand-receptor sequencing, which can identify spatially reachable, expression-supported, and downstream-regulated intercellular communication events in spatial transcriptome data.

[0146] The specific embodiments described herein are merely illustrative of the spirit of the invention. Those skilled in the art to which this invention pertains may make various modifications or additions to the described specific embodiments or use similar methods to substitute them, without departing from the spirit of the invention or exceeding the scope defined by the appended claims.

Claims

1. A method for inferring spatial intercellular communication based on regulator-guided neural stochastic differential equations, characterized in that, Includes the following steps: S1. Data Input and Preprocessing: Input spatial transcriptome expression matrix, spatial coordinates, candidate ligand-receptor database, and prior knowledge of transcription factor regulation; S2. Construction of local ligand-receptor activity: A k-nearest neighbor graph is constructed for spatial coordinates, and local neighborhood expression is integrated according to the distance decay kernel; for candidate ligand-receptor pairs, the product of the ligand expression after neighborhood integration and the receptor expression is used as the local ligand-receptor activity to obtain the spatial activity matrix; S3 and SCENIC regulator function verification: Transcription factors associated with candidate ligand-receptor pairs were extracted from prior transcription factors. The regulatory relationship between transcription factors and target genes was inferred using GRNBoost2. A co-expression module was constructed. Motif-supported regulons were obtained through cisTargetmotif pruning. The activity of each regulon in cells or spots was calculated using AUCell. S4. Dynamic diffusion modeling of neural stochastic differential equations: using standardized spatial coordinates as the initial state, the continuous-time diffusion process of organization dependence is learned through drift network and diffusion network; S5. Calculate the dynamic adjacency matrix: Transform the diffusion endpoint obtained from the neural stochastic differential equation back to the original spatial scale, and calculate the edge weights based on the distance between the diffused transmitter position and the original receiver position, thereby obtaining the dynamic spatial communication topology and forming a dynamic adjacency matrix. S6. Calculate the spatial specificity score: The ligand-receptor activity distribution is propagated on the dynamic adjacency matrix, the consistency of the activity distribution before and after propagation is compared, and the penalty term is calculated by randomly permuting the adjacency background to obtain the ligand-receptor spatial specificity score corrected by random spatial background. S7. Integration of transcription factor functional support: For candidate ligand-receptor pairs with downstream transcription factor annotations and corresponding regulons that can be recognized by SCENIC, their regulon activity is included as a functional enhancer in the final score. S8. Output Analysis Results: Based on the final score, high-confidence ligand-receptor pairs are selected, and the communication network between the sending and receiving cell types is further inferred by combining cell type annotations, while generating a spatially resolved communication field map.

2. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 1, characterized in that, In step S1, let the spatial coordinates of the i-th cell or spot be c. i The expression level of the g-th gene is X. ig The initial spatial distance and initial adjacency weight are defined as follows: ； Where, d ij N represents the Euclidean distance between cell or spot i and j. k(i) Let A represent the set of k nearest neighbors of i. 0 ij A represents the initial spatial adjacency weight, and σ0 represents the distance decay scale parameter; 0 ij It is retained only if j belongs to the k nearest neighbor set of i, otherwise it is set to 0; σ0 is the scale parameter that controls the initial spatial distance decay.

3. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 2, characterized in that, In step S2, the expression levels of each gene are integrated into their local spatial neighborhood to obtain a local microenvironment-sensing expression matrix: ； Among them, X ig This represents the gene expression level after integration within a local spatial neighborhood; for the k-th candidate ligand-receptor pair (l k ,r k Its local ligand-receptor communication activity is defined as: ； Among them, Z ik This indicates the local communication activity of the k-th ligand-receptor pair on the cell or spoti. I (·) This is an indicator function.

4. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 3, characterized in that, In step S3, let Regulon t AUCell represents a regulator centered on transcription factor t and supported by motifs. i (Regulon t The fraction () represents the activity of the regulator in the i-th cell or spot. For candidate ligand-receptor pair k, its regulator support fraction is defined as: ； Among them, TF k Let TF represent the downstream regulator support strength of the k-th candidate ligand-receptor pair; if candidate ligand-receptor pair k has no matching downstream transcription factor, or the corresponding transcription factor does not form an effective regulon, then let TF k =0.

5. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 4, characterized in that, In step S4, the neural stochastic differential equation, with standardized spatial coordinates as its initial state, learns a continuous-time diffusion process: ； Among them, y t This represents the diffusion state initialized by spatial coordinates; f θ For drift networks; g φ For diffusion networks; W t For a Wiener stochastic process, for the k nearest neighbor edge i→j, the self-supervised objective of ligand-receptor affinity is: ； ； q ij Q represents the edge affinity self-supervised target obtained from ligand-receptor expression. 0.99 L represents the 99th percentile normalization factor. SDE Represents the training loss of the neural stochastic differential equation; by minimizing L SDE This ensures that the spatial reachability after diffusion remains consistent with the ligand-receptor expression affinity.

6. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 5, characterized in that, In step S5, the diffusion endpoint y T (i) Denormalize to the original scale coordinates z i And calculate the dynamic adjacency matrix: ； z i A represents the transmitter location mapped back to the original coordinate scale after the diffusion endpoint is reached. dyn ij Represents the communication reachability weights in the dynamic adjacency matrix; the dynamic adjacency matrix preserves the local spatial topological skeleton, but its edge weights are determined by the diffusion reachability learned from neural stochastic differential equations.

7. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 6, characterized in that, In step S6, for the k-th ligand-receptor pair, let G k This represents its spatial activity vector across all cells or spots. The activity vector after dynamic topological propagation is: ； The observation spatial consistency score and the random background score are defined as follows: ； ； Among them, G k Let A represent the spatial activity vector of the k-th ligand-receptor pair. rand b Let b be the background adjacency matrix obtained by the b-th random permutation, and let the ligand-receptor spatial specificity score be defined as: ； Among them, LRSS k denoted as the ligand-receptor spatial specificity score after random background correction, and μ represents the random background penalty coefficient.

8. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 7, characterized in that, In step S7, the space-specific score and regulator support score are integrated into the final communication score: ； Among them, Score k λ represents the final communication score, and λ represents the transcription factor functional support weight.

9. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 8, characterized in that, In step S8, let a be the sending cell type and b be the receiving cell type. Then, the communication strength at the cell type level can be obtained by aggregation of high-confidence candidate ligand-receptor pairs: ； Among them, Spec a (l k ) indicates ligand l k Expression specificity in cell type a Spec b (r k ) represents receptor r k Expression specificity in receiving cell type b Coloc ab (k) represents the colocalization support of ligand-receptor pair k in the sending-receiving cell type spatial neighborhood.

10. The spatial intercellular communication inference method based on regulator-guided neural stochastic differential equations according to claim 1, characterized in that, The spatial transcriptome expression matrix is ​​defined as X∈R (n×p) This indicates that the set of spatial coordinates is represented by C={c i ∈R d } indicates that the candidate ligand-receptor and its downstream transcription factor annotation set is represented by P={(l k ,r k ,T k )}express.

Images