Gene regulatory network inference method and system based on single-cell multi-omics data

Abstract

The invention provides a gene regulatory network inference method and system based on single-cell multi-omics data, a model based on a discrete evolution dynamical system is constructed by integrating paired scRNA-seq and scATAC-seq data, the change of gene expression is modeled into a dynamic process regulated by transcription factor efficacy and chromatin accessibility by using a pseudo time sequence, and the gene regulatory network inference method and system based on the single-cell multi-omics data are obtained. Therefore, the regulation intensity of causality is quantitatively deduced in the context of cell state evolution, and false correlation is remarkably reduced. Besides, the scheme also alleviates the sparsity and noise problems of single cell data through strategies such as adaptive cell grouping, realizes high-precision and high-biological-consistency cell type specific dynamic GRN inference in combination with priori knowledge such as transcription factor binding site matching, reduces the dependence on high-time-sequence-resolution data, and has the advantages of high sensitivity, high accuracy and high biological consistency. And the applicability and robustness in dynamic process analysis such as cell differentiation and state transition are improved.

Description

Technical Field

[0001] This invention relates to the fields of biomedicine, cell gene data processing, and gene screening, and in particular to a method and system for accurately and dynamically inferring cell gene regulatory networks based on single-cell multi-omics data. Background Technology

[0002] The emergence of high-throughput technologies such as single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq) has provided unprecedented insights into transcriptional regulatory mechanisms from the perspectives of transcriptome status and chromatin accessibility, respectively. Gene regulatory networks (GRNs) consist of the regulatory coupling relationships between transcription factors and target genes. Constructing GRNs at single-cell resolution is crucial for understanding the regulatory mechanisms behind cellular heterogeneity and connecting molecular maps with phenotypic expression. However, elucidating complex gene regulatory mechanisms urgently requires effective methods for integrating and analyzing multi-omics data. Existing GRN inference methods can be mainly categorized based on data type and modeling strategy. Methods based on single-omics data have inherent limitations. For example, methods relying solely on scRNA-seq data infer gene expression by calculating statistical dependencies (such as mutual information and regression analysis), but their performance is limited by data sparsity and incomplete modeling of regulatory mechanisms, resulting in lower accuracy in benchmark tests. Methods relying solely on scATAC-seq data construct networks by identifying transcription factor binding sites and enhancers in accessible regions, but because the correlation between chromatin accessibility and gene expression is not a simple linear one, their inference results differ from expression-based GRN inferences. Therefore, integrating multimodal data has become the mainstream strategy for improving the accuracy of GRN inference.

[0003] Against this backdrop, various methods have emerged for integrating scRNA-seq and scATAC-seq data. These integration methods encompass strategies for handling paired data (cells that are identical or in one-to-one correspondence), unpaired data, or metacell-based strategies, and their mathematical frameworks are diverse, including regression models, correlation analysis, probabilistic models, and deep learning. Some tools further integrate pseudo-temporal analysis, aiming to capture more biologically meaningful regulatory relationships. While these integration methods excel in constructing static GRNs to identify potential regulators, their core limitation lies in the fact that they often provide static network snapshots of population average levels or specific states, making it difficult to capture the dynamic evolution of regulatory relationships between transcription factors and target genes during cellular state transitions. In summary, existing technologies have the following limitations: (1) Static models are difficult to capture dynamic regulatory processes: Most current GRN inference methods are based on population averages or state-specific static "snapshots", which cannot effectively characterize the continuous dynamic evolution of the regulatory relationship between transcription factors (TF) and target genes (TG) in processes such as cell differentiation or development, resulting in insufficient analytical capabilities.

[0004] (2) Limited data integration and noise handling capabilities: Although existing methods have attempted to integrate scRNA-seq and scATAC-seq data, when faced with the inherent high sparsity and significant noise of single-cell data, they often rely on linear regression or simplified models (such as linear or idealized dynamic assumptions), making it difficult to accurately quantify the nonlinear regulatory effect of chromatin accessibility on gene expression and easily generating false associations.

[0005] (3) Insufficient ability to infer causality: Existing methods are mostly based on statistical correlation or co-expression relationship between genes, lacking clear event logic to distinguish the true causal regulatory direction, and cannot effectively identify the intensity of TF’s directional regulation of TG.

[0006] Therefore, how to reconstruct GRNs that reflect dynamic and causal regulatory relationships from sparse and noisy paired single-cell multi-omics data remains a problem that urgently needs to be improved in this field. Summary of the Invention

[0007] To address the aforementioned problems in existing technologies, this invention aims to provide a novel solution to overcome the limitations of the static, linear, and non-causal inference models, enabling the reconstruction of GRNs that reflect dynamic, causal regulatory relationships from sparse, noisy paired single-cell multi-omics data. Specifically, this invention provides the following technical solution: On the one hand, this invention provides a method for inferring gene regulatory networks based on single-cell multi-omics data, the method comprising: S1. Preprocess the paired single-cell RNA sequencing data and single-cell ATAC sequencing data to obtain the basic feature matrix (with the same cell name). S2. Based on the basic feature matrix, cells are pseudo-temporally sorted and developmental trajectories are constructed, wherein the developmental trajectory contains multiple branches; S3. Prior knowledge for DEDS model training based on multiple branches; S4. Constructing a DEDS model, including: constructing a discrete evolutionary dynamic system to describe the dynamic evolution of GRN along pseudo-time, the discrete evolutionary dynamic system including the evolution of TF expression, the evolution of TG activity and the evolution of TG expression; and constructing a regression prediction model based on the discrete evolutionary dynamic system to reconstruct GRN. S5. Train the DEDS model; S6. Infer the cell gene regulatory network based on the trained DEDS model.

[0008] Preferably, in step S4, the method for constructing the discrete evolutionary dynamic system is as follows: S41. The evolution of TF expression is represented by the TF expression evolution equation (for potential regulatory relationships). Based on the first In the group Construction of observations Evolutionary mapping , ): in, To influence conversion rate, This is a correction item; Indicates the first In each group Express For the In each group Express The influence of the saturation function; express right The scale of the impact; Used to screen for effective regulatory relationships; The evolution of S42 and TG activity is expressed by the TG activity evolution equation (for potential regulatory relationships). Based on the first In the group Construction of observations Evolutionary mapping , ): in, To influence conversion rate, Transform the regulatory effect The change in the Kth group, This is a correction item; Indicates the first In each group Express For the In each group active The influence of the saturation function; Indicating the K-1th group active; express right The scale of the impact; The evolution of S43 and TG expression is represented by the TG expression evolution equation as follows (for...) Based on the first In the group Construction of observations Evolutionary mapping , ): in, Indicates the first In each group production rate Indicates the first In each group The natural rate of order reduction, Indicates the first In each group The use of induced order reduction rate; This is a correction item; This represents the pseudo-time of cell group K. Represents the (K-1)th group Background requirements.

[0009] Preferably, The calculation method is as follows: in, Control of the interval Attention levels (negatively correlated) are usually set ; Regulation on the range The increase in attention (positively correlated); Control separately In the critical interval The degree and slope of the steep increase in the vicinity; Adjusting the effective regulatory relationship ( The efficiency value of ) Normalization factor; Indicates the intensity of regulation; Preferably, The calculation method is as follows: Model parameters Among them, setting , The calculation method is as follows (for example, here will be used) Abbreviated as ): Among them, model parameters .

[0010] Preferably, The calculation method is as follows: Model parameters ;set up ; The calculation method is as follows (for example, here will be used) Abbreviated as ): Among them, model parameters .

[0011] More preferably: Where model parameters , , All are from and The constituent factors are used as model parameters during training.

[0012] Preferably, the fitting loss Set to: in, This indicates the number of groups in which cells are grouped.

[0013] On the other hand, the present invention also provides a gene regulatory network inference system based on single-cell multi-omics data. This system is used to execute the gene regulatory network inference method based on single-cell multi-omics data as described above. The system includes: The preprocessing module preprocesses paired single-cell RNA sequencing data and single-cell ATAC sequencing data to obtain the basic feature matrix; The temporal and grouping module performs pseudo-temporal sorting of cells based on the basic feature matrix and constructs developmental trajectories, which contain multiple branches. The DEDS model module is based on prior knowledge formed by multiple branches for DEDS model training. The DEDS model is constructed by: building a discrete evolutionary dynamic system to describe the dynamic evolution of the GRN along pseudo-time, the discrete evolutionary dynamic system including the evolution of TF expression, the evolution of TG activity, and the evolution of TG expression; constructing a regression prediction model based on the discrete evolutionary dynamic system to reconstruct the GRN; and training the DEDS model.

[0014] Compared to existing technologies, this invention offers the following significant advantages: By integrating paired single-cell multi-omics data and constructing a discrete evolutionary dynamic system based on pseudo-time series, scDDS effectively overcomes the limitations of traditional methods in dealing with high data noise, sparsity, and nonlinear regulation. This method models regulatory relationships as dynamic processes evolving over time, enabling more accurate identification of causal interactions between transcription factors and target genes and significantly reducing spurious associations. Simultaneously, utilizing chromatin accessibility information enhances the reliability of inferring regulatory directions, exhibiting higher quantitative accuracy and biological consistency in simulating perturbation effects and identifying key regulatory loops compared to relying solely on gene expression data or other multi-omics integration methods. Furthermore, by replacing real-time sampling with pseudo-time analysis, this invention reduces reliance on high-temporal-resolution data, improving its applicability and robustness in analyzing dynamic processes such as cell differentiation and state transitions. Attached Figure Description

[0015] To more clearly illustrate the technical solutions in the embodiments of the present invention or the prior art, the drawings used in the description of the embodiments or the prior art will be briefly introduced below. Obviously, the drawings described below are only some embodiments of the present invention. For those skilled in the art, other drawings can be obtained based on these drawings without creative effort.

[0016] Figure 1 This is a schematic diagram of the single-cell gene regulatory network inference method according to an embodiment of the present invention; Figure 2 A schematic diagram of the transcriptional regulation mechanism; Figure 3 This is a diagram illustrating situations where cell data is missing. Figure 4 This is a schematic diagram illustrating the filtering and screening of cells and genes. Figure 5 Schematic diagram for generating gene expression matrix and gene activity matrix; Figure 6 This is a schematic diagram of the pseudo-time analysis and cell grouping algorithm in the embodiment; Figure 7 This is a schematic diagram illustrating the relationship between TF expression, TG activity, and TG expression. Figure 8 This is a diagram showing the prediction results of the regulation intensity based on the evolutionary process and a binary representation of these results. Figure 9 A schematic diagram for calculating the matching degree of transcription factor binding to promoter (dependent on TFDB and promoter position weight matrix). Figure 10 A schematic diagram showing the TF-TG matching degree calculation results based on TFDB and starter PWM; Figure 11 This is a schematic diagram of the binarized TF-TG matching degree calculation results based on TFDB and starter PWM; Figure 12 This is a schematic diagram illustrating the training of a model based on ibGRN and pbGRN, using a genetic algorithm and gradient descent. Detailed Implementation

[0017] The embodiments of the present invention will now be described in detail with reference to the accompanying drawings. It should be understood that the described embodiments are only a part of the embodiments of the present invention, and not all of them. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without creative effort are within the scope of protection of the present invention.

[0018] Those skilled in the art should understand that the following specific embodiments or implementation methods are a series of optimized configurations listed to further explain the specific content of the invention. These configuration methods can be combined or used in conjunction with each other, unless the invention explicitly states that some or a specific embodiment or implementation method cannot be associated with or used in conjunction with other embodiments or implementation methods. Furthermore, the following specific embodiments or implementation methods are merely optimized configurations and are not intended to limit the scope of protection of the invention.

[0019] To address at least some of the shortcomings of existing technologies, the present invention aims to provide a new technical solution that overcomes the deficiencies of existing static or simplified dynamic models. The primary objective of this invention is to achieve a leap from static association to causal inference, namely, to explicitly define directional transcription factor-target gene interactions based on transcription factor binding logic. The core objective of this invention is to provide a scheme capable of accurately capturing the dynamic evolutionary characteristics of GRNs by integrating paired scRNA-seq and scATAC-seq data to construct a mathematical model based on a discrete evolutionary dynamic system. This model aims to utilize pseudo-time series to model changes in gene expression as a dynamic process regulated by transcription factor efficacy and chromatin accessibility, thereby quantitatively inferring the intensity of causal regulation within the context of cell state evolution. Furthermore, this scheme alleviates the sparsity and noise problems of single-cell data through strategies such as adaptive cell grouping, and combines prior knowledge such as transcription factor binding site matching to ultimately achieve high-precision, highly biologically consistent cell type-specific dynamic GRN inference, providing a powerful tool for revealing the regulatory logic of key biological processes such as cell fate determination. The specific details of this scheme are described below with reference to the accompanying drawings.

[0020] This invention relates to a method for inferring gene regulatory networks (GRNs) based on paired single-cell RNA sequencing (scRNA-seq) and single-cell ATAC sequencing (scATAC-seq) data, and a Discrete Evolutionary Dynamic System (DEDS) model. This method is named scDEDS. The following explains the meanings of the abbreviations used in this method: DEDS – Discrete Evolutionary Dynamic System; scRNA-seq – Single-cell RNA sequencing; scATAC-seq – Single-cell ATAC sequencing; TF – Transcription factor; TG – Target gene; GRN – Gene regulatory network; TSS – Transcription start site; iGRN – Initial gene regulatory network (initial GRN); pGRN – Predicted gene regulatory network (predicted GRN); ibGRN – Initial binary GRN; pbGRN – Predicted binary GRN.

[0021] The core process of this invention is as follows: Figure 1 As shown, the main steps include: data preprocessing, pseudo-temporal analysis and cell grouping, DEDS model construction, and GRN regression prediction. Overall, this approach starts with paired scRNA-seq and scATAC-seq data, uses pseudo-temporal analysis to group cells into metacells along their developmental trajectories, and then establishes discrete dynamic equations to describe the evolution of transcription factor (TF) expression, target gene (TG) activity, and TG expression. Finally, a regression model predicts the regulatory intensity and outputs GRNs. This approach innovatively combines dynamic modeling with machine learning optimization, significantly improving the accuracy and interpretability of GRN inference. The following sections explain each of the four key steps in detail.

[0022] I. Data Preprocessing The data preprocessing stage aims to extract high-quality feature matrices from raw sequencing data and define key regulatory elements, namely transcription factors (TFs) and target genes (TGs).

[0023] First, define the promoter region: using the transcription start site (TSS) as a reference, delineate the promoter region as 0 to 0.5 meters upstream and downstream of the TSS. Within the range of base pairs ( (This is a preset distance parameter, in bp). This definition serves as a benchmark for chromatin accessibility analysis, such as... Figure 2 As shown.

[0024] Acquire paired single-cell RNA sequencing data (scRNA-seq) and single-cell ATAC sequencing data (scATAC-seq) (preferably with more than 500 cells in this embodiment) and convert them into, for example, Seurat objects (e.g., using the Seurat v5.3.0 tool). If necessary, batch effect removal can be performed on the scRNA-seq data.

[0025] Annotate chromatin fragments in scATAC-seq data (for example, you can use the annotatePeak function of the ChIPseeker tool), considering only genes on the same strand as the fragment (set the parameter sameStrand=True), and ignoring the overlapping region between the fragment and the TSS (set the parameter ignoreOverlap=True) to ensure annotation specificity.

[0026] TG was defined as genes expressed in scRNA-seq data and annotated as promoter-accessible in scATAC-seq data. TF was retrieved from the JASPAR2024 database (e.g., using the getMatrixSet function in JASPAR2024 v0.99.6) for all TF genes. Furthermore, the promoter base sequence for each TG was obtained (using, for example, the getSeq function in Biostrings tools), and these sequences were annotated using the JASPAR2024 database. Genes of interest (including TG and TF) were used to screen scRNA-seq and scATAC-seq data, such as... Figure 3 As shown.

[0027] Subsequently, cells and genes were filtered: scRNA-seq data with a missing value frequency less than [value missing] were retained. Cells with a missing value rate of less than Genes to remove low-quality entries, such as Figure 4 As shown, where and This is a threshold parameter, with a value ranging from 0 to 1.

[0028] Finally, a gene expression matrix is ​​generated based on the selected scRNA-seq data, and a gene activity matrix is ​​generated based on the selected scATAC-seq data (using tools such as the GeneActivity function of Signac). Figure 5 As shown. All data can be stored as Seurat objects, forming the foundational feature matrix for subsequent analysis.

[0029] II. Pseudo-time analysis and cell grouping To capture the continuous evolution of cell states, this embodiment utilizes the Monocle tool to perform pseudo-temporal sequencing of cells and constructs developmental trajectories based on gene expression similarity. Through machine learning algorithms, individual cells are arranged on one or more continuous trajectories, thereby inferring the relative position of each cell in the process; this relative position is the pseudo-temporal sequence. The specific process is shown in Table 1 below.

[0030] Table 1 Pseudo-time sorting Steps, Objectives, Core Functions, and Key Parameters: Brief Explanation 1. Object Construction: Creates a container for Monocle2 analysis: `newCellDataSet(expression_matrix, phenoData, featureData, expressionFamily=negbinomial.size())`. This constructs the core object from the expression matrix, cell phenotype information, and gene information. For UMI-based count data, the `negbinomial.size()` family of functions is typically used. 2. Data preprocessing, standardization, and estimation of variability estimateSizeFactors(cds) estimateDispersions(cds) Similar to normalization in other processes, this prepares for subsequent differential analysis. 3. Selecting Ordering Genes: Select the gene set used to construct the trajectory. `setOrderingFilter(cds, ordering_genes)` tells Monocle2 which genes (`ordering_genes`) to use to learn the trajectory. This gene set selection strategy directly affects the trajectory morphology. 4. Data Dimensionality Reduction: Reduce dimensionality to construct the trajectory skeleton. `reduceDimension(cds, max_components=2, method='DDRTree')`: By default, the DDRTree method is used to reduce the data to 2 dimensions for visualization of the trajectory structure. 5. Cell Sort: Calculate the pseudotime value for each cell. `orderCells(cds, root_state=NULL)`: Sort cells according to gene expression patterns, assigning each cell a State (state / branch) and a Pseudotime (pseudotime value, starting from 0).

[0031] The trajectory branches can be manually divided based on the dimensionality reduction results generated by the monocle. The starting point is at pseudo-time point 0, and branches are divided by combining each bifurcation point. Each trajectory branch represents a cell fate direction.

[0032] In this embodiment, the trajectory branches are divided based on the dimensionality reduction result, and the division method is as follows: For single-cell RNA sequencing data, pseudo-temporal sorting was performed using the whole genome, and the data was manually partitioned based on dimensionality reduction visualization. Each branch yields the pseudo-time for each cell in each branch. .

[0033] Each branch reflects multiple possible fate paths during differentiation, corresponding to different biological state transitions; for the current branch, it includes... Individual cells, single cells according to their pseudo-time ( Sort in ascending order. expression level, The activity and expression level are denoted as follows: , and . This indicates the i-th TF contained in the current branch. This indicates the j-th TG contained in the current branch.

[0034] For branches ( ), which includes TF and Each TG is denoted as TG ( )and ( (When performing branch-specific analysis, superscripts will be omitted for simplicity.) , recorded as , , , In the unambiguous case, and It can be further simplified to and .

[0035] It should be further explained here that, in this embodiment, after the cells are grouped, each cell group synthesizes an independent metacell. The Kth cell group is the Kth metacell. In other words, in some expressions, "cell group" and "metacell" can be regarded as equivalent.

[0036] Cells are sorted by pseudo-time values ​​and then grouped into cell groups to alleviate data sparsity issues. In this embodiment, the grouping algorithm employs an adaptive strategy: iteratively grouping cells to ensure that the number of non-missing values ​​for each representative feature in each group is higher than a given threshold, such as... Figure 6 As shown. Specifically, based on the total number of cells. Determine the minimum group size Its calculation model is as follows: Where [⋅] represents the floor function. Parameters a, b, and c need to be determined based on three sets of... The data points are fitted (points are manually selected based on data characteristics and analysis needs). In this embodiment, a fitting accuracy of at least 10 is required. −3 Record the first In the group , and The number of non-missing values ​​are respectively , and ,satisfy The adaptive grouping algorithm divides cells through the following iterative process: starting from initialization, each allocation... Each cell is transferred to a new group, and the new group is tested to see if it meets the criteria. If the conditions are not met, the new group size is expanded when there are sufficient remaining cells; otherwise, it is merged with an existing group. This iterative process continues, updating the remaining cell pool after each grouping until all cells have been divided. In this embodiment, the ideal state is... ;like It is necessary to adjust the fitting points or increase them. To mitigate the problem of excessive zero values. If the results are still unsatisfactory or the target gene is not detected, the data should be discarded. In this embodiment, representative features include TF expression, TG activity, and TG expression. The parameters mentioned above include: total cell number. minimum group size , TF number TG number Cell group number Number of cells per group ( ).

[0037] Ultimately, for primordial cells (i.e., cell group K), its pseudo-time is defined as (Supplementary definition of initial) This means taking the maximum pseudo-time of all individual cells in cell group K as the pseudo-time of that cell group. Its representative characteristics (Transcription factor expression levels) (Target gene activity) and (Target gene expression level) is defined as the corresponding non-zero value in cells within the group. , and The average value, This represents the pseudo-time of a single cell within a cell group. This step discretizes the continuous cell states, providing time-series data for dynamic modeling, wherein preferably: .

[0038] III. DEDS Model Construction The DEDS model is a crucial core component of this scheme, describing the pseudo-temporal evolution of TF expression (TFE), TG activity (TGA), and TG expression (TGE) using discrete dynamic equations. This model is based on a negative feedback mechanism of transcriptional regulation, namely: each... The TF and the first TG regulation relationship By regulation intensity Quantification, its evolution is driven by time-delayed feedback between TF expression (TFE), TG activity (TGA), and TG expression (TGE).

[0039] After completing the trajectory branching, prior knowledge is prepared for DEDS model training. In this embodiment, the prior knowledge based on the branches is formed by independently analyzing each branch. The position weight matrix (PWM) of the transcription factors is obtained (this can be achieved using, for example, the getMatrixSet function in JASPAR2024). Then, the matching degree between each transcription factor and the target gene promoter is calculated (this can be achieved using, for example, the matchPWM function in Biostrings v2.74.1, with the parameter set to min_score=80%). The regulatory strength is then defined. This is the ratio of the matching degree to 1.05 times the maximum matching degree. This normalization process avoids the value from being too close to 1 during training. Using genes as nodes... The initial gene regulatory network (iGRN) was constructed using these as edge weights. A threshold was then used to... Binarization is performed to obtain This leads to the initial binary GRN (ibGRN): when hour Otherwise, it is 0. Positive weighted edges are considered positive samples, denoted as . The rest are negative samples, denoted as To construct a holistic cell type network, branch-specific iGRNs and ibGRNs were merged with the number of branched cells as the weight, such as... Figure 9-10 As shown, missing values ​​due to genetic differences are filled with zeros.

[0040] In addition, the constructed DEDS model requires prerequisite settings: for the regulation relationship Its regulatory intensity is denoted as A higher value indicates a stronger regulatory effect. The physical quantities involved... ( (expression level) ( (Active) and ( All expression levels are positive, and there is a time-delayed interaction. Under the assumption of non-extreme, background steady state, The background demand value is a positive constraint with gentle fluctuations. Similarly, Its demand value Regulation is subject to each The positive impact, which is caused by Parameterization. Assume that within a branch, each pair... and All have regulatory strength Its regulatory mechanism includes two feedback loops: specifically, feedback... By adjusting To offset and Differences - When Insufficient supply leads to an upward trend, while excessive supply leads to a downward trend. This trend follows... Scaling is performed. The total change integrates all Contributions. Similar feedback. Based on Regulation and integrate all The impact of transcriptional cooperativity and chromatin accessibility. Furthermore, considering the roles of transcriptional synergy and chromatin accessibility, we hypothesize the existence of a positively bounded regulatory mechanism. All assumptions apply to the relationship with the metacell. Within the pseudo-time interval (and its corresponding cell group), each variable is denoted as... , , , , .

[0041] In this embodiment, the definition of saturation in the evolution equation is based on the form of the Hill function, which is used to describe the saturation effect of regulation: in, This typically represents the strength of upstream regulatory factors. This represents the normalized biological response level. D is the apparent dissociation constant, which biologically represents the concentration of x required for the response y to reach half of its maximum value (i.e., y=0.5). A smaller D value indicates higher affinity or a more sensitive response. The synergistic parameter determines the shape of the dose-response curve; when n>1, it indicates a positive synergistic effect, the curve is S-shaped, and the modulation process has switching characteristics; when n=1, it degenerates into a simple hyperbolic non-synergistic combination. Here, the effect on saturation is set based on the Hill function form.

[0042] In this embodiment, the evolution equations in the model are as follows: Figure 7 As shown: First, the evolution of TF expression is represented by the TF expression evolution equation (for potential regulatory relationships). Based on the first In the group Construction of observations Evolutionary mapping , ): Among them, factors affecting conversion rate Transform the regulatory effect The change in the Kth group, This is a correction item; Indicates the first In each group Express For the In each group Express Influence saturation function, model parameters To facilitate formula expression, we set... , The calculation method is as follows (for example, here will be used) Abbreviated as ): The above formula represents right Impact scale, model parameters ; Used to screen for effective regulatory relationships: in, Control of the interval Attention levels (negatively correlated) are usually set ; Regulation on the range The increase in attention (positively correlated); joint control In the critical interval The degree of sharp increase in the vicinity (and) Positive correlation) and steep slope (with (positive correlation), general settings ; Adjusting the effective regulatory relationship ( The efficiency value of ( ). Normalization factor make sure And make It is a relatively small positive value. Within the extended interval. Above, function Must meet The constraints are as follows. Among them, model parameters , , , , All are from and The constituent factors are used as model parameters during training.

[0043] Secondly, similarly, the evolution of TG activity is represented by the TG activity evolution equation (for potential regulatory relationships). Based on the first In the group Construction of observations Evolutionary mapping , ): Among them, factors affecting conversion rate Transform the regulatory effect The change in the Kth group, This is a correction item; Indicates the first In each group Express For the In each group active Influence saturation function, model parameters To facilitate formula expression, we set... , The calculation method is as follows (for example, here will be used) Abbreviated as ): The above formula represents right Impact scale, model parameters .set up Where model parameters , , All are from and The constituent factors are used as model parameters during training.

[0044] Finally, the evolution of TG expression is represented by the TG expression evolution equation (for...). Based on the first In the group Construction of observations Evolutionary mapping , ): in, Indicates the first In each group ; production rate; Indicates the first In each group The natural degradation rate; Indicates the first In each group The use of induced degradation rate, They represent The maximum production rate, maximum natural degradation rate, and maximum use-induced degradation rate were used as model parameters during training. for Production saturation; for The natural degradation saturation; for The use of induced degradation saturation; This is a correction item. It is trained as model parameters.

[0045] These equations together constitute the discrete evolutionary dynamics system of this embodiment, which describes the evolutionary dynamics of GRN along pseudo-time.

[0046] Based on the discrete evolutionary dynamics system established above, this embodiment constructs a regression prediction model to reconstruct the GRN, specifically as follows: Figure 8 , Figures 11-12 As shown. The regression prediction model is trained as follows: First, the edges in iGRN are divided into a training set (70%), a validation set (15%), and a test set (15%) to ensure a balance between positive and negative samples. The training objective is to minimize the loss function. In this embodiment, the loss function is set as follows: Among them, the fitting loss Set to: In the above formula, Indicates the number of cell groups. This indicates the calculation of variance.

[0047] In this embodiment, during model training, a genetic algorithm combined with gradient descent is used to optimize parameters until the validation set performance is stable. Ultimately, the predicted regulatory strength is determined. Used to construct branch-specific GRNs, abbreviated as pGRNs, and then weighted by cell count to merge them into overall GRNs, such as Figure 12 As shown. The specific training and reconstruction evaluation process is as follows: The first step is to divide the edges in iGRN into training, validation, and test sets. In this embodiment, the default division ratio is 70%, 15%, and 15%, but this can be changed according to requirements to ensure a balanced ratio of positive and negative samples in each set. This division scheme is then assigned to each branch for independent processing.

[0048] The second step is to compare the current parameters with the positive samples in the training set. and its corresponding Input loss function L. First, run a multi-generation genetic algorithm, for example, 30 generations. In this embodiment, the population size can be set to 20. The crossover probability decreases from 0.95 to 0.8 as the loss decreases, and the mutation probability decreases from 0.75 to 0.3. Then, perform multi-step gradient descent with an adaptive learning rate (for example, set to 3 steps). For the current parameters, if it is the first iteration, the initial value can be used; otherwise, the result of the previous iteration is inherited.

[0049] The third step is to perform a test on each training sample, based on the current parameters and their corresponding parameters. Minimize computation of , recorded as All Values ​​sorted in ascending order , To match the size of the training set, pad with 0s and 0s at the beginning and end. form There are several intervals. The midpoint of each interval is used as the binary classification threshold. The candidate values ​​are evaluated based on the following criteria: if The prediction is then a negative sample. Conversely, positive samples are considered positive samples, and the confusion matrix is ​​then calculated. Statistics. This will produce the optimal... The interval is denoted as .

[0050] The fourth step is to obtain the validation set samples using the method described in step three. and However, only those that fall into within The remaining After sorting, use and Use the endpoints to fill in the intervals, and take the midpoint of the interval as... Candidate values. Calculated according to the method in step three. Record the best value and its corresponding interval .

[0051] Fifth, obtain the test set samples using the same procedure as in step five. .

[0052] Step 6, if If there is no improvement after multiple iterations, the number of iterations can be set manually according to the needs (e.g., 4 times). If so, the early stopping mechanism is activated; otherwise, the process returns to the second step and enters the next round of training.

[0053] Step 7: After training is completed, As a branch pGRN, the branch pGRN is merged with cell number as the weight to obtain cell type pGRN and predicted binary GRN, i.e., pbGRN, whose non-zero value is set to 1.

[0054] Step 8: Finally, calculate the initial binary gene regulatory network ibGRN ( ) and pbGRN ( The confusion matrix is ​​calculated using accuracy, recall, precision, F1 score, AUC, and... The coefficients are used to evaluate the network reconstruction performance.

[0055] Combining the discrete evolutionary dynamics system and regression prediction mentioned above, a DEDS model is constructed, thereby building a trainable DEDS model to complete the inference of gene regulatory networks.

[0056] Furthermore, this solution can also be implemented systematically. In a preferred embodiment, the structure of this system can be set as follows: The preprocessing module preprocesses paired single-cell RNA sequencing data and single-cell ATAC sequencing data (with the same cell name) to obtain the basic feature matrix. The temporal and grouping module performs pseudo-temporal sorting of cells based on the basic feature matrix and constructs developmental trajectories, which contain multiple branches. The DEDS model module is based on prior knowledge formed by multiple branches for DEDS model training. The DEDS model is constructed by: building a discrete evolutionary dynamic system to describe the dynamic evolution of the GRN along pseudo-time, the discrete evolutionary dynamic system including the evolution of TF expression, the evolution of TG activity, and the evolution of TG expression; constructing a regression prediction model based on the discrete evolutionary dynamic system to reconstruct the GRN; and training the DEDS model.

[0057] Furthermore, the system also includes an output module for outputting the inference results of the cell gene regulatory network for user use.

[0058] Furthermore, the configuration of the system described above in this solution can execute the gene regulatory network inference method based on single-cell multi-omics data as given in the above embodiments.

[0059] In another embodiment, this solution can be implemented using a device, which may include corresponding modules that perform one or more steps in the various embodiments described above. Therefore, each or more steps in the various embodiments can be performed by a corresponding module, and the electronic device may include one or more of these modules. A module may be one or more hardware modules specifically configured to perform a corresponding step, or implemented by a processor configured to perform a corresponding step, or stored in a computer-readable medium for implementation by a processor, or implemented through some combination thereof.

[0060] This device can be implemented using a bus architecture. The bus architecture can include any number of interconnect buses and bridges, depending on the specific application of the hardware and overall design constraints.

[0061] Those skilled in the art will understand that all or part of the processes in the above embodiments can be implemented by a computer program instructing related hardware. The program can be stored in a computer-readable storage medium, and when executed, it can include the processes of the embodiments of the above methods. The storage medium can be a magnetic disk, optical disk, read-only memory (ROM), or random access memory (RAM), etc.

[0062] The above description is merely a specific embodiment of the present invention, but the scope of protection of the present invention is not limited thereto. Any variations or substitutions that can be easily conceived by those skilled in the art within the technical scope disclosed in the present invention should be included within the scope of protection of the present invention. Therefore, the scope of protection of the present invention should be determined by the scope of the claims.

Claims

1. A method for inferring gene regulatory network based on single-cell multi-omics data, characterized in that, The method comprises: S1, preprocessing paired single-cell RNA sequencing data and single-cell ATAC sequencing data to obtain a basic feature matrix; S2, based on the basic feature matrix, the cells are pseudo-time sorted, and a development trajectory is constructed, the development trajectory comprising a plurality of branches; S3, forming a DEDS model training priori knowledge based on the plurality of branches; S4, constructing a DEDS model, comprising: constructing a discrete evolution dynamic system for describing the kinetic evolution of GRN along the pseudo-time, the discrete evolution dynamic system comprising the evolution of TF expression, the evolution of TG activity and the evolution of TG expression; constructing a regression prediction model based on the discrete evolution dynamic system to reconstruct the GRN; S5, training the DEDS model; S6, inferring the cell gene regulation network based on the trained DEDS model.

2. The method of claim 1, wherein, In S1, the preprocessing comprises: S11, obtaining paired single-cell RNA sequencing data and single-cell ATAC sequencing data with the same cell label; S12, demarcate the promoter upstream and downstream of the transcription start site chromatin fragments in single-cell ATAC sequencing data, only considering genes on the same strand as the chromatin fragments, ignoring the overlapping region of the chromatin fragments and the transcription start site; the gene existing in the single-cell RNA sequencing data and annotated as the promoter in the single-cell ATAC sequencing data is TG; is a preset distance parameter; S13, obtaining the base sequence of the promoter of each TG, obtaining all TFs; annotating the base sequence of the promoter; S14, based on TG and TF, screening single-cell RNA sequencing data and single-cell ATAC sequencing data; generating a basic feature matrix based on the screened single-cell ATAC sequencing data.

3. The method of claim 1, wherein, In S2, the branches are divided according to the dimension reduction result, and the division method is: For single-cell RNA sequencing data, perform pseudo-time ordering with whole gene, divide into branches, get the pseudo-time of each cell in each branch ; For the current branch, contains cells, individual cells are arranged in ascending order of their pseudo-time , expression amount, activity and expression amount respectively recorded as , and ; , indicates the i-th TF contained in the current branch, indicates the j-th TG contained in the current branch.

4. The method of claim 1, wherein, The S2 further comprises, after the cells are ordered by pseudo-time, the cells are further grouped into groups of cells, wherein the first group contains cells, and satisfies , represents the total number of cells, and the grouping method is: S21, according to the total number of cells determining a minimum group size : where [•] denotes the ceiling function, and parameters a, b, c are determined according to multiple sets of Data points were fitted to obtain; S22, record the number of non-missing values of groups , and are respectively , and ; based on the following iterative process to divide the cells: each time allocate cells to a new group, and check whether the new group meets the decision condition ; if the condition is not met, then expand the size of the new group when the remaining cells are sufficient, otherwise merge with the existing group; the iterative process continues until all cells are divided; , and respectively represent the expression of , activity and expression, represent the i-th TF contained in the current branch, represent the j-th TG contained in the current branch; S23, for the cell group with pseudo-time defined as .

5. The method of claim 1, wherein, In S3, the priori knowledge is obtained in the following way: S31, independently analyzing each branch: obtaining the position weight matrix of the transcription factor, and then calculating the matching degree of each TF and TG promoter; S32, set control intensity is the ratio of the matching degree and 1.05 times the maximum matching degree; TG and TF are taken as nodes, and the gene regulation network is constructed as the edge weight S33, passing through a threshold value obtained by performing a binarization process , and further obtaining a binary GRN, wherein is an element of the binary GRN; S34, constructing an overall cell type network, taking the number of branch cells as the weight, merging the initial gene regulation network and the initial binary GRN with branch specificity, wherein the missing values caused by gene differences are filled with zeros.

6. The method of claim 1, wherein, In S4, the premise of constructing the DEDS model is set as: For regulatory relationship , the regulatory strength is denoted as ; represents the ith TF contained in the current branch, represents the jth TG contained in the current branch; Within one branch, and contains two feedback loops: feedback loop is adjusted to counteract the difference from , which creates an upregulation tendency when is insufficient and a downregulation tendency when in excess, scaled by ; feedback loop regulates based on and integrates the effects of all . The above prerequisites apply to the pseudo-time intervals associated with the cell and its corresponding cell group; denotes the expression amount, denotes the expression amount, denotes the activity, denotes the background demand value of denotes the background demand value of 7. The method of claim 1, wherein, In S4, the method for constructing the discrete evolution dynamic system is: S41, the evolution of TF expression is represented by a TF expression evolution equation: wherein, is an influence on conversion rate, is a correction term; represents the influence of the expression in the group on the expression in the group ; represents the scale of the influence ; for screening effective regulatory relationships; S42, the evolution of TG activity is represented by a TG activity evolution equation: wherein, is an influence on conversion, is a correction term; denotes the activity in the Kth group of K groups of expression of the activity in the Kth group of K groups of activity a saturation function; denotes the activity in the K-1th group of K groups denotes the influence size on S43, the evolution of TG expression is represented by a TG expression evolution equation: in, Indicates the first In each group production rate Indicates the first In each group The natural rate of order reduction, Indicates the first In each group The use of induced order reduction rate; This is a correction item; This represents the pseudo-time of cell group K. Represents the (K-1)th group Background requirements.

8. The method of claim 1, wherein, In S5, training the DEDS model includes training the regression prediction model, and in the training of the regression prediction model, the training target is to minimize the loss function, and the loss function is: wherein, represents the fitting loss, represents the variance.

9. The method of claim 8, wherein, In S5, the training and evaluation method of the regression prediction model is: S51, dividing the edges in the gene regulation network into a training set, a validation set and a test set, and assigning the division scheme to each branch for independent processing; S52, compare the current parameters with the positive samples of the training set and their corresponding input the loss function, and train by using the genetic algorithm and gradient descent method; S53、for each training sample, based on the current parameters and the corresponding compute the fitting loss of , denoted as ; arrange all the fitting loss values in ascending order as , where is the training set size, and the first and last are padded with 0 and to form intervals; take the midpoint of each interval as a candidate value of the binary classification threshold , and evaluate the standard: if , predict as negative sample, otherwise as positive sample; then compute the confusion matrix and statistic, and record the interval that produces the optimal statistic as ; represents the regulation intensity; S54, obtaining the verification set samples according to the method in step S53 and only retaining the within ; filling the interval with and as endpoints after sorting the remaining , taking the midpoint of the interval as the candidate value, calculating the statistical quantity, and recording the optimal value and the corresponding interval ; S55, obtaining the operation process of step S54 on the test set samples ; S56, if If the successive iterations are not improved, the early stopping mechanism is initiated, otherwise, return to S52 for the next round of training. S57. After training is completed, The branch prediction GRN is used as a weight; the branch prediction GRNs are merged with cell number as the weight to obtain cell type prediction GRN and prediction binary GRN. S58, finally calculating the confusion matrix of the initial binary GRN and the predicted binary GRN, and evaluating the network reconstruction performance.

10. A gene regulatory network inference system based on single-cell multi-omics data, characterized in that, The system is used to execute the method of any one of claims 1-9, and the system comprises: a preprocessing module, which preprocesses paired single-cell RNA sequencing data and single-cell ATAC sequencing data to obtain a basic feature matrix; a time sequence and grouping module, which sorts cells based on the basic feature matrix, and constructs a development trajectory, the development trajectory comprising a plurality of branches; The DEDS model module is based on a plurality of branches to form prior knowledge for training of the DEDS model; the DEDS model is constructed, including: constructing a discrete evolution dynamic system for describing dynamics evolution of the GRN along pseudo-time, the discrete evolution dynamic system including evolution of TF expression, evolution of TG activity and evolution of TG expression; constructing a regression prediction model based on the discrete evolution dynamic system to reconstruct the GRN; and training the DEDS model.

Images