Method for jointly deducing dynamic cell communication and cell state transition rate

Abstract

The invention provides a method for jointly deducing dynamic cell communication and a cell state transition rate. Relates to the technical field of biological information. The method comprises the following steps: extracting candidate ligands, candidate receptors and characteristic genes from space transcriptome data to be processed; the method comprises the following steps: screening nodes which have an interaction relationship with candidate ligands, candidate receptors and characteristic genes from a pre-constructed prior database, and constructing a multi-layer signal network of which the structure is ligand-receptor-transcription factor-target genes; establishing a gene regulation kinetic model according to a ligand-receptor action relationship, a receptor-transcription factor action relationship and a transcription factor-target gene action relationship in the multilayer signal network; iteratively optimizing parameters to be estimated in the gene regulation and control kinetic model; and obtaining the change rate of the expression quantity of the target gene in the receiving cell based on the potential time after iterative optimization and the parameters to be estimated after iterative optimization. And combined analysis of dynamic cell communication and cell differentiation tracks can be realized.

Description

Technical Field

[0001] This invention relates to the field of bioinformatics, and in particular to a method for jointly inferring dynamic cell communication and cell state transition rates. Background Technology

[0002] The coordination of cell fate assignment and functional specialization in multicellular systems stems from the dynamic interaction between intrinsic cellular regulatory programs and microenvironment-mediated intercellular signal transduction. Understanding this signal transduction mechanism is key to comprehending the operation of multicellular systems.

[0003] To reveal the patterns of intercellular signal transduction, several cell-cell communication (CCC) inference methods have been developed to infer intercellular ligand-receptor (LR) interactions or intracellular ligand / receptor-target gene (TG) regulation from scRNA-seq data.

[0004] However, current methods for inferring intercellular communication are mainly designed to infer static signaling networks between different cell types or states. The limitations of this approach stem from the snapshot nature of single-cell and spatial transcriptomics data, which does not take into account the impact of cell communication on cell state transitions. As a result, it is impossible to capture the dynamic evolution of intercellular communication throughout the entire process of cell state transition (CST). Summary of the Invention

[0005] Therefore, it is necessary to provide a method for jointly inferring dynamic cell communication and cell state transition rate to address the above-mentioned technical problems, so as to solve the deficiency of existing technologies in being unable to capture the dynamic evolution of intercellular communication during cell state transition and to achieve synchronous inference of dynamic cell communication and cell state transition rate.

[0006] The following technical solution is adopted in this specification: This specification provides a method for jointly inferring dynamic cell communication and cell state transition rates, including: Obtain unprocessed spatial transcriptome data containing information on cell spatial location and gene expression, and extract candidate ligands, candidate receptors, and characteristic genes from the unprocessed spatial transcriptome data; Nodes that interact with candidate ligands, candidate receptors, and characteristic genes are screened from a pre-constructed prior database, and a multi-layer signaling network with a structure of ligand-receptor-transcription factor-target gene is constructed based on the nodes. Based on the ligand-receptor interaction, receptor-transcription factor interaction, and transcription factor-target gene interaction in multilayer signaling networks, a gene regulation kinetic model is established to describe the rate of change in target gene expression. The parameters to be estimated in the gene regulation kinetics model are iteratively optimized. In each iteration, potential time is allocated to cells in the spatial transcriptome data to be processed, so as to provide pseudo-time sequence of cells in the spatial transcriptome data on the developmental time axis to the iterative optimization process of the parameters to be estimated in the gene regulation kinetics model. Based on the iteratively optimized potential time and parameters to be estimated, the rate of change in target gene expression, which characterizes the rate of cell state transition in the receiving cells of the candidate receptor, is obtained.

[0007] Furthermore, the mathematical expression of the gene regulation kinetic model is: in, Indicates transcription factor At any moment The activity; ,when When =1, it means the first... Ligand-receptor signaling regulates transcription factors ;when When =0, it means the first... Ligand-receptor signaling does not regulate transcription factors ; Indicates the first Ligand-receptor signaling regulates transcription factors The maximum activation rate; Indicates in Time of the first The signal intensity of each ligand-receptor; Indicates the first Ligand-receptor signaling activates transcription factors Reaching 1 / 2 Required concentration; Indicates transcription factor The level of self-expression; Indicates transcription factor The degradation coefficient; This indicates the number of ligand-receptor pairs in a multilayer signaling network; Indicates target gene The amount of expression; , =1 indicates transcription factor Regulatory target genes ;when =0 indicates transcription factor Unregulated target genes ; Indicates transcription factor Target genes The maximum regulation rate; Indicates transcription factor target genes Reaching 1 / 2 Required concentration; Indicates target gene The activity of upstream transcription factors; This indicates the number of transcription factors in a multilayer signaling network.

[0008] Furthermore, the formula for calculating the ligand-receptor signal intensity is as follows: in, Indicates receiving cell The first received Each ligand-receptor signal intensity; Indicates sending cells With receiving cells Spatial distance between them; Indicates sending cells The expression level of the corresponding ligand; Indicates receiving cell The expression level of the corresponding receptor; Indicates the number of cells sent; Indicates the signal intensity of the diffuse ligand-receptor; This represents the neighborhood relationship of each receiving cell; This indicates the signal intensity of the contact ligand-receptor.

[0009] Furthermore, the iterative optimization of the parameters to be estimated in the gene regulation kinetics model specifically includes: I. Based on prior knowledge or by randomly setting pre-defined weight parameters in a deep neural network and parameters to be estimated in a gene regulation dynamics model; II. Potential time of cell allocation in spatial transcriptome data to be processed; III. Based on the latent time and the weight parameters in the preset deep neural network, the predicted value of target gene expression is obtained, and the Euclidean distance between the predicted value of target gene expression and the observed value of target gene expression in cells in the spatial transcriptome data to be processed is used as the data loss. IV. Based on the parameters to be estimated in the potential time and gene regulation dynamics model, obtain the model fitting value of the target gene expression rate change rate; use the error between the model fitting value of the target gene expression rate change rate and the target gene expression rate change rate obtained by automatic differentiation through a preset deep neural network as the physical loss; V. Obtain the total loss based on data loss and physical loss; VI. Obtain the parameter gradient of the total loss, and update the weight parameters in the preset deep neural network and the parameters to be estimated in the gene regulation dynamics model along the parameter gradient direction; VII. Repeat steps II to VI. When the total loss is less than the preset loss threshold or the preset maximum iteration step is reached, output the parameters to be estimated in the gene regulation kinetics model.

[0010] Furthermore, the potential time for cell allocation in the spatial transcriptome data to be processed specifically includes: Divide the interval [0, 1] into several discrete time points; Obtain observations of target gene expression levels in cells from the spatial transcriptome data to be processed; Several discrete time points and target gene expression level observations are input into a pre-set deep neural network, and the pre-set deep neural network outputs the target gene expression level prediction values ​​corresponding to several discrete time points. The discrete time point corresponding to the predicted target gene expression value with the smallest Euclidean distance from the observed target gene expression value is taken as the potential time of the cell corresponding to the observed target gene expression value.

[0011] Furthermore, the formula for calculating the data loss is as follows: in, Indicates the number of cells; Indicates the number of target genes; This indicates that the preset deep neural network is based on latent time. and the weight parameters in the preset deep neural network The output is the predicted value of target gene expression level; Indicates potential time The corresponding observed values ​​of target gene expression levels.

[0012] Furthermore, the formula for calculating the physical loss is as follows: in, This indicates that the preset deep neural network is based on latent time. The output is the predicted value of target gene expression level; This indicates differentiation with respect to time; , =1 indicates transcription factor Regulatory target genes ;when =0 indicates transcription factor Unregulated target genes ; Indicates transcription factor Target genes The maximum regulation rate; Indicates transcription factor target genes Reaching 1 / 2 Required concentration; Indicates target gene The activity of upstream transcription factors; This indicates the number of transcription factors in a multilayer signaling network.

[0013] Furthermore, the total loss is calculated as follows: in, and These are the first adjustment parameter and the second adjustment parameter, respectively. Total loss Minimization is the optimization objective, and the weight parameters in the preset deep neural network and the parameters to be estimated in the gene regulation dynamics model are jointly iteratively optimized.

[0014] Furthermore, the extraction of candidate ligands, candidate receptors, and characteristic genes from the spatial transcriptome data to be processed includes: In the spatial transcriptome data, all cells are clustered to obtain cell type clusters. The average gene expression level of the cell type cluster is calculated. Genes with an average gene expression level exceeding a preset threshold and an expression ratio exceeding a specified ratio in the sending cells are selected as candidate ligands. In the spatial transcriptome data, all cells are clustered to obtain cell type clusters. The average gene expression level of the cell type clusters is calculated. Genes with an average gene expression level exceeding a preset threshold and an expression ratio exceeding a specified ratio in the receiving cells are selected as candidate receptors. Differentially expressed genes of cell type clusters identified by the Seurat-based FindMarkers function are used as feature genes.

[0015] Furthermore, the node-based multilayer signaling network, which is a ligand-receptor-transcription factor-target gene structure, includes: Obtain pre-constructed prior databases, including: a ligand-receptor database storing the interaction relationship between ligands and receptors, a receptor-transcription factor database storing the interaction relationship between receptors and transcription factors, and a transcription factor-target gene database storing the interaction relationship between transcription factors and target genes. Search the ligand-receptor database for receptor nodes that match candidate ligands and candidate receptors, and construct a ligand-receptor subnetwork; Feature genes are mapped to target gene nodes in the transcription factor-target gene database, and a transcription factor-target gene subnetwork is constructed based on the target gene nodes and the transcription factor nodes associated with the target gene nodes. The receptor-transcription factor subnetwork was constructed by searching the receptor-transcription factor database for transcription factor nodes associated with target gene nodes and the receptor nodes linked to the transcription factor nodes. By sequentially linking the ligand-receptor subnetwork, the receptor-transcription factor subnetwork, and the transcription factor-target gene subnetwork, a multilayer signaling network with the structure of ligand-receptor-transcription factor-target gene is obtained.

[0016] This invention constructs a multi-layer signaling network to obtain complete cell communication links corresponding to the spatial transcriptome data to be processed, including intercellular ligand-receptor signal transduction, intracellular receptor-transcription factor transduction, and intracellular transcription factor-target gene regulation. Based on the ligand-receptor, receptor-transcription factor, and transcription factor-target gene interactions in the multi-layer signaling network, a gene regulation kinetic model is constructed. This model outputs the rate of change in target gene expression, which can quantify the regulatory effect of intercellular communication signal transduction on the dynamic evolution of target gene expression. During the iterative optimization of the model's parameters, potential time is allocated to cells in the data to be processed, providing pseudo-timelines of cells on the developmental timeline, injecting time dimension information into parameter optimization, and ensuring that the model can accurately fit the dynamic changes in gene expression with the developmental process. Attached Figure Description

[0017] The accompanying drawings, which are included to provide a further understanding of this application and form part of this application, illustrate exemplary embodiments and are used to explain this application, but do not constitute an undue limitation of this application. In the drawings:

[0018] Figure 1 This specification provides a flowchart illustrating a method for jointly inferring dynamic cell communication and cell state transition rates. Figure 2 This is a schematic diagram of a framework for constructing a multi-layer signal network provided by the present invention; Figure 3 This is a schematic diagram of a framework for constructing a gene regulation dynamics model provided by the present invention; Figure 4 This is a schematic diagram of a framework for iterative optimization of parameters to be estimated in a gene regulation dynamics model based on PINN-CELL, provided by the present invention. Figure 5 This is a schematic diagram of a framework for optimizing result output provided by the present invention; Figure 6 This is a schematic diagram of an iterative optimization process for parameters to be estimated in a gene regulation dynamics model, provided by the present invention. Detailed Implementation

[0019] To make the objectives, technical solutions, and advantages of this specification clearer, the technical solutions of this application will be clearly and completely described below in conjunction with specific embodiments and corresponding drawings. Obviously, the described embodiments are only a part of the embodiments of this application, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments in this specification without creative effort are within the scope of protection of this application.

[0020] The technical solution provided by this invention can be applied to disease mechanism research, drug development, developmental and regenerative biology, and other scenarios, enabling cell communication and inference of cell state transition rates in multicellular systems. In multicellular systems, the coordination of cell fate determination and functional specialization stems from the dynamic interaction between intrinsic cellular regulatory programs and microenvironment-mediated intercellular signals. Revealing the hierarchical regulatory mechanisms of spatiotemporal cell state transitions—encompassing molecular-scale signal transduction, cellular-scale state changes, and tissue-scale microenvironment patterns—remains a fundamental and challenging problem in biology. This challenge arises both from technical limitations in data acquisition and from the inherent complexity of the nonlinear coupling between intracellular and intercellular regulatory processes. The coordination of cell fate assignment and functional specialization in multicellular systems originates from the dynamic interaction between intrinsic cellular regulatory programs and microenvironment-mediated intercellular signal transduction. Revealing the hierarchical regulatory logic of spatiotemporal cell state transitions, encompassing molecular signal transduction, cellular state changes, and tissue-scale microenvironment patterns, remains a fundamental and challenging problem in biology. These challenges stem from both technical limitations in data acquisition and the inherent complexity of the nonlinear coupling between intracellular and intercellular regulatory processes.

[0021] Methods for inferring cell differentiation trajectories based on single-cell transcriptome data have emerged. These mainstream methods can be divided into two categories: trajectory inference and RNA velocity inference. Trajectory inference typically reconstructs the order of individual cells along differentiation pathways by calculating similarity or distance indices between cells in a low-dimensional manifold. RNA velocity inference, as an emerging method, infers the dynamic trends of gene expression based on the abundance and ratio between unspliced ​​immature mRNA and spliced ​​mature mRNA, thereby predicting the future state of cells. RNA velocity methods reveal the directionality of cell fate determination and improve trajectory inference. However, methods for inferring cell state transitions based on scRNA-seq data have limitations in resolving spatial developmental patterns of cells. With the advent of spatial transcriptomic (ST) technologies, which provide not only gene expression information but also spatial location information of cells or sites, some methods have been developed to infer spatial CST from ST data. In the aforementioned studies, neither machine learning-based trajectory inference methods nor kinetic modeling-based RNA velocity methods linked intercellular interactions with CST dynamics, resulting in a lack of mechanistic and functional insights into intercellular communication in cell fate determination in existing CST inference methods. Currently, numerous CCC inference methods have been developed to infer intercellular ligand-receptor interactions or intracellular ligand / receptor-target gene regulation from scRNA-seq data. Some CCC inference methods specifically for ST data have also been developed. However, current CCC inference methods are primarily designed to infer static signaling networks between different cell types or states. This limitation stems from the snapshot nature of single-cell and spatial transcriptomics data, thus failing to infer the dynamic changes in intercellular communication during cell state transitions. In other words, existing methods for inferring cell differentiation trajectories do not consider the impact of cell communication on cell state transitions, preventing them from inferring the dynamic changes in intercellular communication during these transitions. Therefore, integrating CCC mechanisms with CST dynamics within a cohesive framework is a crucial and fundamental challenge.

[0022] Based on this, this invention provides a method for jointly inferring dynamic cell communication and cell state transition rates, using the CCC-driven approach to infer cell differentiation trajectories from ST data. This method combines intercellular LR signaling activity with intracellular transcription factor (TF) activity cascades to construct a multi-scale nonlinear dynamic model to capture the dynamic changes in cell state within the tissue microenvironment. Furthermore, to accurately estimate the parameters in the model, this invention also designs PINN-CELL, a physically-informed neural network co-evolutionary learning algorithm that can simultaneously optimize model parameters and the latent time of cells. By combining CCC mechanism modeling with trajectory learning from ST data, this invention provides strong support for understanding the multi-scale dynamic regulation principles that determine cell fate.

[0023] The following is combined Figures 1-6 This invention describes a method for jointly inferring dynamic cell communication and cell state transition rates.

[0024] Figure 1 This is a flowchart illustrating a method for jointly inferring dynamic cell communication and cell state transition rates provided by the present invention, as shown below. Figure 1 As shown, the method includes: S100: Obtain unprocessed spatial transcriptome data containing information on cell spatial location and gene expression; extract candidate ligands, candidate receptors, and characteristic genes from the unprocessed spatial transcriptome data.

[0025] The steps for acquiring spatial transcriptome data to be processed may include: selecting target tissue sample slices and fixing the slices onto a dedicated chip with spatial barcode markers; using molecular biology steps such as nucleic acid probe hybridization, reverse transcription, and amplification to bind the mRNA of each cell to the spatial barcode on the chip; completing sequencing using a high-throughput sequencing platform to generate raw sequencing data; and then obtaining standardized spatial transcriptome data to be processed through bioinformatics preprocessing such as data quality control, sequence alignment, gene counting, and spatial coordinate calibration.

[0026] For example, in spatial transcriptome data, all cells are clustered to obtain cell type clusters. The average gene expression level of each cell type cluster is calculated. Genes with an average expression level exceeding a preset threshold and an expression ratio exceeding a specified proportion in the sending cells are selected as candidate ligands. Taking a preset receptor threshold min.expr=0.05 as an example, genes with an average expression level exceeding 0.05 and an expression ratio exceeding a specified proportion (e.g., 10%) in the sending cells can be defined as candidate ligands.

[0027] For example, in spatial transcriptome data, all cells are clustered to obtain cell type clusters. The average gene expression level of each cell type cluster is calculated. Genes with an average expression level exceeding a preset threshold and an expression ratio exceeding a specified proportion in receiving cells are identified as candidate receptors. Similarly, using a preset ligand threshold min.expr=0.05 as an example, genes with an average expression level exceeding 0.05 and an expression ratio exceeding a specified proportion (e.g., 10%) in receiving cells can be defined as candidate receptors.

[0028] For example, differentially expressed genes from cell type clusters identified by the FindMarkers function in Seurat are used as feature genes. Furthermore, the set of genes of interest can be included in the screening of candidate ligands, receptors, or feature genes, and together with the aforementioned candidate ligands, receptors, or feature genes, a multilayer signaling network can be constructed.

[0029] S200. Nodes that interact with candidate ligands, candidate receptors and characteristic genes are screened from a pre-constructed prior database, and a multi-layer signaling network with a structure of ligand-receptor-transcription factor-target gene is constructed based on the nodes.

[0030] For example, the pre-constructed prior databases include: LigRecDB, which stores the interaction between ligands and receptors; RecTFDB, which stores the interaction between receptors and transcription factors; and TFTGDB, which stores the interaction between transcription factors and target genes.

[0031] The pre-constructed prior database can be built by integrating multiple publicly available and literature-validated molecular interaction databases. Specifically, LigRecDB contains 920 ligands, 751 receptors, and 3659 non-redundant ligand-receptor interactions; RecTFDB contains 751 receptors, 525 transcription factors, and 17450 non-redundant receptor-transcription factor connections; and TFTGDB contains 525 transcription factors, 23021 target genes, and 373501 transcription factor-target gene regulatory relationships. During data integration, a restarted random walk algorithm is used to integrate these data sources to construct a multi-layered prior knowledge network with global connectivity.

[0032] For example, Figure 2 This is a schematic diagram of a framework for constructing a multi-layer signal network provided by the present invention, which can be seen in the following figure. Figure 2 As shown, spatial transcriptome data to be processed was obtained. This data includes: a gene-cell expression matrix to represent gene expression in different cells and cell coordinates: ( (), ) Then, a multilayer signaling network is constructed using a prior database. This network consists of three layers: the top layer is the ligand-receptor interaction layer, used to characterize the ligand-receptor (LR) signal transduction relationship. The middle layer represents the regulatory relationship between the receptor and the transcription factor TF, where the ligand-receptor regulatory relationship on the transcription factor is represented by matrix A=[ ] indicates that when When =1, it means the first... Ligand-receptor signaling regulates transcription factors ;when When =0, it means the first... Ligand-receptor signaling does not regulate transcription factors The lower layer is the regulatory layer between transcription factor TF and target gene (TG), used to characterize the regulatory relationship between transcription factor, target gene TF, and TG. This regulatory relationship is represented by matrix B=[ ]express, =1 indicates transcription factor Regulatory target genes ;when =0 indicates transcription factor Unregulated target genes The intermediate layer facilitates the transmission of ligand-receptor signals to transcription factor-target gene regulation.

[0033] For example, after identifying candidate ligands in the sending cell and candidate receptors and characteristic genes in the receiving cell from the input ST data to be processed, the Lig-Rec database is searched for receptor nodes matching the candidate ligands and ligand nodes matching the candidate receptors. Based on the results obtained through bidirectional searches of "candidate ligand → receptor node matching the candidate ligand" and "candidate receptor → ligand node matching the candidate receptor," an initial Lig-Rec subnetwork of ligands and receptors is constructed. Characteristic genes are mapped to target gene nodes in the TFTG database of transcription factors and target genes. Based on the target gene nodes and the transcription factor nodes linked to the target gene nodes, a TF-TG subnetwork of transcription factors and target genes is constructed. The Rec database of receptors and transcription factors is searched for... Transcription factor nodes were searched in TFDB. Based on the transcription factor nodes and the receptor nodes linked to them, a receptor-transcription factor Rec-TF subnetwork was constructed. Using the receptor nodes in the receptor-transcription factor Rec-TF subnetwork, the initial ligand-receptor Lig-Rec subnetwork was refined, such as by removing redundant spurious pairs, to obtain the ligand-receptor Lig-Rec subnetwork. The ligand-receptor Lig-Rec subnetwork, the receptor-transcription factor Rec-TF subnetwork, and the transcription factor-target gene TF-TG subnetwork were sequentially linked to obtain a multilayer signaling network with the structure of ligand-receptor-transcription factor-target gene. This multilayer signaling network can systematically and accurately characterize the complete link of intercellular signal transduction and intracellular signal transduction, and intuitively present the entire process of cell communication.

[0034] S300. Based on the ligand-receptor interaction, receptor-transcription factor interaction, and transcription factor-target gene interaction in multilayer signaling networks, establish a gene regulation kinetic model to describe the rate of change in target gene expression.

[0035] For example, Figure 3 This is a schematic diagram of a framework for constructing a gene regulation dynamics model provided by the present invention, as shown below. Figure 3 As shown, the steps for establishing a gene regulation kinetic model to describe the rate of change in target gene expression may include: firstly, obtaining the signal intensity matrix of intercellular ligand-receptor pairs (LRS) (the signal intensity matrix shows different intercellular LRS1, LRS2...). LRS J Signal intensity of ligand-receptor pairs Distribution), expression matrix of transcription factor TF (the expression matrix of transcription factor TF shows TF1, TF2 in different cells) TF M Expression levels of transcription factors Distribution of TG1 and TG2 in different cells) and expression matrix of target gene TG (expression matrix of target gene TG shows TG1, TG2 and TG2 in different cells) TG G Expression levels of target genes The distribution of transcription factors is then analyzed, and an ordinary differential equation model describing the transcriptional dynamics driven by intercellular communication—the gene regulation dynamics model—is constructed. This model includes equations for transcription factor activity. This is used to characterize the dynamic regulation of transcription factor activity by ligand-receptor signal intensity; equations targeting target gene expression levels. This model is used to characterize the regulation of target gene expression dynamics by transcription factor activity. The set of parameters to be estimated in the model is as follows: These equations enable a quantitative description of transcriptional dynamics driven by intercellular communication. Specifically:

[0036] The quantification process of intercellular ligand-receptor (LR) signal intensity is as follows: Based on the different mechanisms of ligand-receptor interaction, LR pairs are divided into two categories: diffusion-type and contact-type, and signal intensity calculation models are established for each. First, the diffusion-type LR signal attenuates with distance in space. Based on the law of mass action and the diffusion partial differential equation, the receiving cell signal is obtained. The first received The ligand-receptor signal intensities are: ;

[0037] ; in, Indicates sending cells With receiving cells Spatial distance between them Indicates sending cells The expression level of the corresponding ligand; Indicates receiving cell The expression level of the corresponding receptor.

[0038] Contact-type LR signals are transmitted only between physically nearest neighbor cells. The Delaunay triangulation algorithm is used to identify the neighborhood relationships of each cell. The signal strength is defined as: .

[0039] The combined definitions of the two types of signal strength are as follows: .

[0040] For example, the process of modeling the dynamics of TG expression driven by intercellular LR signaling and intracellular TF activity is as follows: First, assume that the first A TF at a moment The activity is Its variation is affected by the upstream LR signal strength With TF's own expression level Co-regulation. Using the Hill function to characterize the nonlinear relationship between LR signal intensity and TF activity, and introducing a linear TF self-regulation term, we can obtain:

[0041] ; in, , indicating the first Does ligand-receptor signaling regulate transcription factors? ,when When =1, it means the first... Ligand-receptor signaling regulates transcription factors ;when When =0, it means the first... Ligand-receptor signaling does not regulate transcription factors ; Indicates the first Ligand-receptor signaling regulates transcription factors The maximum activation rate; Indicates the first A ligand-receptor in Signal strength at any given time; Indicates the first Ligand-receptor signaling activates transcription factors Reaching 1 / 2 Required concentration; Indicates transcription factor The level of self-expression; Indicates transcription factor The degradation coefficient; This indicates the number of ligand-receptor pairs in a multilayer signaling network.

[0042] Secondly, assuming target gene The expression level is Its rate of change is determined by the upstream TF activity. The decision is influenced by its own degradation. Similarly, the Shell function is used to characterize the nonlinear regulatory relationship between them:

[0043] ; in, , representing transcription factors Does it regulate target genes? , =1 indicates transcription factor Regulatory target genes ;when =0 indicates transcription factor Unregulated target genes ; Indicates transcription factor Target genes The maximum regulation rate; Indicates transcription factor target genes Reaching 1 / 2 Required concentration; Indicates target gene The activity of upstream transcription factors; This indicates the number of transcription factors in a multilayer signaling network.

[0044] The constructed CCC driver gene regulatory dynamics model can be uniformly represented as: ; ; The set of parameters to be estimated is as follows: .

[0045] S400. Iteratively optimize the parameters to be estimated in the gene regulation kinetics model. In each iteration, by allocating potential time to cells in the spatial transcriptome data to be processed, pseudo-time sequence of cells in the spatial transcriptome data on the developmental time axis is provided to the iterative optimization process of the parameters to be estimated in the gene regulation kinetics model.

[0046] Optimizing gene regulation dynamics models based on ordinary differential equations (ODEs) is challenging because these models are dynamic, nonlinear, and high-dimensional, and ST data do not contain explicit cellular timeframes. Therefore, an iterative optimization algorithm based on Physics-Informed Neural Network-based CoEvolution Learning (PINN-CELL) is needed to jointly optimize model parameters and latent time variables. Other neural network algorithms can also be used to achieve this step.

[0047] The entire process utilizes the PINN-CELL algorithm to jointly estimate cell latent time and model parameters, deeply integrating data-driven approaches with the mathematical model (ODE) to characterize the dynamic evolution of cell states over time. In each iteration of PINN-CELL, a physics-informed neural network (PINN) is used to estimate the kinetic parameters in the model. For example, Figure 4 This is a schematic diagram of a framework for iterative optimization of parameters to be estimated in a gene regulation kinetics model based on PINN-CELL, as provided by this invention. Figure 4As shown, the steps in each iteration of PINN-CELL include: Step 1, allocating cell potential time: Starting with the root cell (defined as the cell that is at the earliest stage on the developmental timeline and can serve as the starting point for cell state transition, based on the gene expression patterns and spatial location information of cells in the spatial transcriptome), all cells are allocated cell-specific potential time, ranging from the initial time. Extending to the final moment In the above model, all learnable parameters can be divided into three groups, namely, cell-specific latent time. Gene-specific parameters and weight parameters of deep neural networks (DNN) Specifically, PINN-CELL first assigns potential time points to each cell by minimizing the Euclidean distance between the observed gene expression values ​​and the predicted values ​​estimated by the DNN. Step 2: Calculate data loss: [The text abruptly ends here, likely due to an incomplete sentence or missing information.] and the expression level of target genes Inputting the data into a deep neural network (DNN), the network calculates the predicted expression levels of target genes. and the observed expression levels of target genes The difference is used to obtain data loss; Step 3: Calculate physical loss: Calculate the model fitting value of the target gene expression rate change. The rate of change in target gene expression calculated automatically by a pre-set deep neural network Error between With error As the physical loss. Step 4: Integrate the data loss and physical loss (ODE loss) to obtain the final total loss. Minimize the total loss as the optimization objective, and adjust the weight parameters in the deep neural network based on this loss. Parameters to be estimated in gene regulation dynamics models Perform joint iterative optimization until the model converges and the gene-specific parameters are completed. and weight parameters of deep neural networks (DNN) The estimate.

[0048] S500: Based on the potential time and parameters to be estimated after iterative optimization, the rate of change in target gene expression, which characterizes the rate of cell state transition, is obtained in the receiving cells where the candidate receptor is located.

[0049] For example, by combining the optimized potential time and the parameters to be estimated into the gene regulation kinetics model, the rate of change in the expression level of the target gene in the receiving cell under the action of the signal can be accurately calculated, which can help in the subsequent analysis of the spatiotemporal differentiation trajectory of the cell. Figure 5 This is a schematic diagram of a framework for optimizing output provided by the present invention, such as... Figure 5 As shown, the final result can be:

[0050] Cellular state transition velocity (CST): In a two-dimensional coordinate system (x, y), the velocity distribution of cell state transition (CST) is displayed through a set of oriented points. Different points represent the velocity differences and directional trends of cells during the state transition process.

[0051] Cell pseudotime: In a two-dimensional coordinate system (x, y), the pseudotime distribution of cells on the developmental timeline is represented by points of different colors. The color gradient (from 0 to 1) reflects the order of the pseudotime sequence.

[0052] Gene / TF dynamics: Two heatmaps show the dynamic changes in the expression of target gene TG and transcription factor TF over pseudotime. The color intensity of the heatmap reflects the expression level. The horizontal axis is pseudotime and the vertical axis is gene / transcription factor.

[0053] Multilayer signaling network: Presents the topology of a multilayer signaling network, including ligands, receptors, transcription factors (TFs) and target genes (TGs), as well as the regulatory relationships (edges) between them, visually demonstrating the complete signaling pathway from ligand-receptor signal transduction to transcription factor-target gene regulation.

[0054] This invention constructs a multi-layered signaling network to obtain complete cell communication pathways corresponding to the spatial transcriptome data to be processed, including intercellular ligand-receptor signal transduction, intracellular receptor-transcription factor transduction, and intracellular transcription factor-target gene regulation. Based on the ligand-receptor, receptor-transcription factor, and transcription factor-target gene interactions within the multi-layered signaling network, a gene regulation kinetic model is constructed. This model can determine the rate of change in target gene expression, thereby inferring the dynamic changes in cell state within the tissue microenvironment. Furthermore, this invention can also infer the dynamic changes in cell communication during cell state transitions, achieving joint analysis of dynamic cell communication and cell differentiation trajectories.

[0055] Based on the above Figure 1 The illustrated embodiment, for example, in the above-described S400, when iteratively optimizing the parameters to be estimated in the gene regulation kinetics model, can be found in... Figure 6As shown, Figure 6 This specification provides a schematic diagram of a network framework for jointly inferring dynamic cell communication and cell state transition rates. The iterative optimization process includes: I. Based on prior knowledge or by randomly setting the weight parameters in a deep neural network and the parameters to be estimated in a gene regulation dynamics model.

[0056] II. Potential time for cell allocation in spatial transcriptome data to be processed.

[0057] Since the spatial transcriptome only contains "spatial location + gene expression" and lacks real developmental time information, the essence of allocating potential time is to construct a pseudo-time sequence to simulate the sequence of cell development. This pseudo-time sequence provides a constraint framework for the evolution of cell state over time for gene regulation dynamics models, allowing the model to iteratively optimize parameters based on temporally adjacent changes in cell expression, thus conforming to real regulatory laws.

[0058] For example, the step of allocating potential time may include: dividing the [0,1] interval into several discrete time points; obtaining the observed values ​​of target gene expression levels of cells in the spatial transcriptome data to be processed; inputting the several discrete time points and the observed values ​​of target gene expression levels into a preset deep neural network, and outputting the predicted values ​​of target gene expression levels corresponding to the several discrete time points through the preset deep neural network; and taking the discrete time point corresponding to the predicted value of target gene expression level with the smallest Euclidean distance from the observed value of target gene expression level as the potential time of the cell corresponding to the observed value of target gene expression level.

[0059] III. Based on the latent time and the weight parameters in the preset deep neural network, the predicted value of target gene expression is obtained, and the Euclidean distance between the predicted value of target gene expression and the observed value of target gene expression in cells in the spatial transcriptome data to be processed is used as the data loss.

[0060] For example, the formula for calculating data loss is: ; in, Indicates the number of cells; Indicates the number of target genes; This indicates that the preset deep neural network is based on latent time. and the weight parameters in the preset deep neural network The output is the predicted value of target gene expression level; Indicates potential time The corresponding observed values ​​of target gene expression levels.

[0061] IV. Based on the parameters to be estimated in the potential time and gene regulation dynamics model, obtain the model fitting value of the target gene expression rate change rate; use the error between the model fitting value of the target gene expression rate change rate and the target gene expression rate change rate obtained by automatic differentiation through a preset deep neural network as the physical loss.

[0062] For example, the formula for calculating physical loss is: ; in, This indicates that the preset deep neural network is based on latent time. The output is the predicted value of target gene expression level; This indicates differentiation with respect to time; , =1 indicates transcription factor Regulatory target genes ;when =0 indicates transcription factor Unregulated target genes ; Indicates transcription factor Target genes The maximum regulation rate; Indicates transcription factor target genes Reaching 1 / 2 Required concentration; Indicates target gene The activity of upstream transcription factors; This indicates the number of transcription factors in a multilayer signaling network.

[0063] V. Obtain the total loss based on data loss and physical loss.

[0064] For example, the formula for calculating the total loss is: ; in, and These are the first adjustment parameter and the second adjustment parameter, respectively.

[0065] VI. Obtain the parameter gradient of the total loss, and update the weight parameters in the preset deep neural network and the parameters to be estimated in the gene regulation dynamics model along the parameter gradient direction.

[0066] VII. Repeat steps II to VI. When the total loss is less than the preset loss threshold or the preset maximum iteration step is reached, output the parameters to be estimated in the gene regulation kinetics model.

[0067] Those skilled in the art will understand that all or part of the processes in the methods of the above embodiments can be implemented by a computer program instructing related hardware. The computer program can be stored in a non-volatile computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. Any references to memory, storage, databases, or other media used in the embodiments provided in this application can include at least one of non-volatile and volatile memory. Non-volatile memory can include read-only memory (ROM), magnetic tape, floppy disk, flash memory, or optical storage, etc. Volatile memory can include random access memory (RAM) or external cache memory. By way of illustration and not limitation, RAM can be in various forms, such as static random access memory (SRAM) or dynamic random access memory (DRAM), etc.

[0068] The technical features of the above embodiments can be combined in any way. For the sake of brevity, not all possible combinations of the technical features in the above embodiments are described. However, as long as there is no contradiction in the combination of these technical features, they should be considered to be within the scope of this specification.

Claims

1. A method for jointly inferring dynamic cell communication and cell state transition rates, characterized in that, include: Obtain unprocessed spatial transcriptome data containing information on cell spatial location and gene expression, and extract candidate ligands, candidate receptors, and characteristic genes from the unprocessed spatial transcriptome data; Nodes that interact with candidate ligands, candidate receptors, and characteristic genes are screened from a pre-constructed prior database, and a multi-layer signaling network with a structure of ligand-receptor-transcription factor-target gene is constructed based on the nodes. Based on the ligand-receptor interaction, receptor-transcription factor interaction, and transcription factor-target gene interaction in multilayer signaling networks, a gene regulation kinetic model is established to describe the rate of change in target gene expression. The parameters to be estimated in the gene regulation kinetics model are iteratively optimized. In each iteration, potential time is allocated to cells in the spatial transcriptome data to be processed, so as to provide pseudo-time sequence of cells in the spatial transcriptome data on the developmental time axis to the iterative optimization process of the parameters to be estimated in the gene regulation kinetics model. Based on the iteratively optimized potential time and parameters to be estimated, the rate of change in target gene expression, which characterizes the rate of cell state transition in the receiving cells of the candidate receptor, is obtained.

2. The method of jointly inferring dynamic cellular communication and cell state transition rates of claim 1, wherein, The mathematical expression for the gene regulation dynamics model is: wherein, represents a transcription factor at time ; , when =1, represents the th ligand-receptor signal regulated transcription factor ; when =0, represents the th ligand-receptor signal unregulated transcription factor ; represents the th ligand-receptor signal regulated transcription factor 's maximum activation rate; represents the th ligand-receptor's signal strength at time ; represents the th ligand-receptor signal regulated transcription factor 's concentration required to reach 1 / 2 ; represents a transcription factor 's self-expression level; represents a transcription factor 's degradation coefficient; represents the number of ligand-receptor pairs in a multi-layer signal network; represents the expression level of the target gene ; , =1 represents that the transcription factor regulates the target gene ; when =0 represents that the transcription factor does not regulate the target gene ; represents the maximum regulation rate of the transcription factor on the target gene ; represents the concentration required by the transcription factor to make the target gene reach 1 / 2 ; represents the activity of the upstream transcription factor of the target gene ; represents the number of transcription factors in the multi-layer signal network.

3. The method of jointly inferring dynamic cellular communication and cell state transition rates of claim 2, wherein, The formula for calculating the ligand-receptor signal intensity is as follows: in, Indicates receiving cell The first received Each ligand-receptor signal intensity; Indicates sending cells With receiving cells Spatial distance between them; Indicates sending cells The expression level of the corresponding ligand; Indicates receiving cell The expression level of the corresponding receptor; Indicates the number of cells sent; Indicates the signal intensity of the diffuse ligand-receptor; This represents the neighborhood relationship of each receiving cell; This indicates the signal intensity of the contact ligand-receptor.

4. The method of jointly inferring dynamic cellular communication and cell state transition rates of claim 1, wherein, The iterative optimization of the parameters to be estimated in the gene regulation kinetics model specifically includes: I. Based on prior knowledge or by randomly setting pre-defined weight parameters in a deep neural network and parameters to be estimated in a gene regulation dynamics model; II. Potential time of cell allocation in spatial transcriptome data to be processed; III. Based on the latent time and the weight parameters in the preset deep neural network, the predicted value of target gene expression is obtained, and the Euclidean distance between the predicted value of target gene expression and the observed value of target gene expression in cells in the spatial transcriptome data to be processed is used as the data loss. IV. Based on the parameters to be estimated in the potential time and gene regulation dynamics model, obtain the model fitting value of the target gene expression rate change rate; use the error between the model fitting value of the target gene expression rate change rate and the target gene expression rate change rate obtained by automatic differentiation through a preset deep neural network as the physical loss; V. Obtain the total loss based on data loss and physical loss; VI. Obtain the parameter gradient of the total loss, and update the weight parameters in the preset deep neural network and the parameters to be estimated in the gene regulation dynamics model along the parameter gradient direction; VII. Repeat steps II to VI. When the total loss is less than the preset loss threshold or the preset maximum iteration step is reached, output the parameters to be estimated in the gene regulation kinetics model.

5. The method of jointly inferring dynamic cellular communication and cell state transition rates of claim 4, wherein, The potential time of cell allocation in the spatial transcriptome data to be processed specifically includes: Divide the interval [0, 1] into several discrete time points; Obtain observations of target gene expression levels in cells from the spatial transcriptome data to be processed; Several discrete time points and target gene expression level observations are input into a pre-set deep neural network, and the pre-set deep neural network outputs the target gene expression level prediction values ​​corresponding to several discrete time points. The discrete time point corresponding to the predicted target gene expression value with the smallest Euclidean distance from the observed target gene expression value is taken as the potential time of the cell corresponding to the observed target gene expression value.

6. The method of jointly inferring dynamic cellular communication and cell state transition rates of claim 4, wherein, The formula for calculating the data loss is: in, Indicates the number of cells; Indicates the number of target genes; This indicates that the preset deep neural network is based on latent time. and the weight parameters in the preset deep neural network The output is the predicted value of target gene expression level; Indicates potential time The corresponding observed values ​​of target gene expression levels.

7. The method for jointly inferring dynamic cell communication and cell state transition rates as described in claim 6, characterized in that, The formula for calculating the physical loss is: in, This indicates that the preset deep neural network is based on latent time. The output is the predicted value of target gene expression level; This indicates differentiation with respect to time; , =1 indicates transcription factor Regulatory target genes ;when =0 indicates transcription factor Unregulated target genes ; Indicates transcription factor Target genes The maximum regulation rate; Indicates transcription factor target genes Reaching 1 / 2 Required concentration; Indicates target gene The activity of upstream transcription factors; This indicates the number of transcription factors in a multilayer signaling network.

8. The method for jointly inferring dynamic cell communication and cell state transition rates as described in claim 7, characterized in that, The total loss is calculated as follows: in, and These are the first adjustment parameter and the second adjustment parameter, respectively. Total loss Minimization is the optimization objective, and the weight parameters in the preset deep neural network and the parameters to be estimated in the gene regulation dynamics model are jointly iteratively optimized.

9. The method for jointly inferring dynamic cell communication and cell state transition rates as described in claim 1, characterized in that, The extraction of candidate ligands, candidate receptors, and characteristic genes from the spatial transcriptome data to be processed includes: In the spatial transcriptome data, all cells are clustered to obtain cell type clusters. The average gene expression level of the cell type cluster is calculated. Genes with an average gene expression level exceeding a preset threshold and an expression ratio exceeding a specified ratio in the sending cells are selected as candidate ligands. In the spatial transcriptome data, all cells are clustered to obtain cell type clusters. The average gene expression level of the cell type clusters is calculated. Genes with an average gene expression level exceeding a preset threshold and an expression ratio exceeding a specified ratio in the receiving cells are selected as candidate receptors. Differentially expressed genes of cell type clusters identified by the Seurat-based FindMarkers function are used as feature genes.

10. The method for jointly inferring dynamic cell communication and cell state transition rates as described in claim 1, characterized in that, The node-based multilayer signaling network, with a structure of ligand-receptor-transcription factor-target gene, includes: Obtain pre-constructed prior databases, including: a ligand-receptor database storing the interaction relationship between ligands and receptors, a receptor-transcription factor database storing the interaction relationship between receptors and transcription factors, and a transcription factor-target gene database storing the interaction relationship between transcription factors and target genes. Search the ligand-receptor database for nodes that match candidate ligands and candidate receptors, and construct a ligand-receptor subnetwork; Feature genes are mapped to target gene nodes in the transcription factor-target gene database, and a transcription factor-target gene subnetwork is constructed based on the target gene nodes and the transcription factor nodes associated with the target gene nodes. The receptor-transcription factor subnetwork was constructed by searching the receptor-transcription factor database for transcription factor nodes associated with target gene nodes and the receptor nodes linked to the transcription factor nodes. By sequentially linking the ligand-receptor subnetwork, the receptor-transcription factor subnetwork, and the transcription factor-target gene subnetwork, a multilayer signaling network with the structure of ligand-receptor-transcription factor-target gene is obtained.

Images