Method for constructing gene regulatory network based on time delay of single-cell data

Abstract

The present invention proposes a method for constructing a gene regulatory network with time delay based on single-cell data. This method calculates the variance of the expression values of genes based on single-cell data to obtain gene expression data, and calculates the maximum mutual information between regulatory genes and target genes based on time delay to construct a gene regulatory network. The advantages of the present invention are as follows: introducing the time delay factor into single-cell data, calculating the mutual information between the expression values of regulatory genes and target genes, and the mutual information between the expression values of regulatory genes and the rates of target genes, obtaining the maximum value of the mutual information under different time delays, and using the highest mutual information value as the interaction score between genes. The present invention improves the accuracy of gene regulatory network inference and provides a new perspective for exploring the connections between genes.

Description

Technical Field

[0001] The present invention belongs to the field of bioinformatics technology, and relates to a method for constructing a gene regulatory network. Specifically, it relates to a method for constructing a gene regulatory network with time delay based on single-cell data, which can be used for mining gene regulatory relationships. Background Art

[0002] In an organism, the gene expression level of a cell affects specific physiological activities. The expression level of a gene affects the expression of other genes, and at the same time is also affected by other genes. This complex gene regulatory relationship is the basis for the normal life activities of an organism. Therefore, mining gene regulatory relationships has become an extremely important goal in biology. Single-cell RNA sequencing technology can detect hundreds of thousands of cells, with a resolution accurate to each cell. This sequencing technology can not only measure the expression level of each cell, but also retain the heterogeneity of the cells, obtaining single-cell data. The gene expression data in single-cell data contains rich information. In an organism, the expression of one gene may affect the expression value or expression rate of another gene. Existing methods for constructing gene regulatory networks either only consider that one gene affects the expression value of another gene, or only consider that one gene affects the expression rate of another gene. In addition, in an organism, one gene controls the expression of another gene through its own products (RNA and proteins), and the processes of transcribing into RNA and processing proteins require a certain amount of time. Therefore, the phenomenon of time delay in gene regulation always exists in an organism. Currently, the vast majority of methods for constructing gene regulatory networks ignore the existence of time delay, resulting in inaccurate calculation of the regulatory intensity between genes and low accuracy in inferring gene regulatory networks.

[0003] For example, in the article "Gene regulation inference from single-cell RNA-seq data with linear differential equations and velocity inference" published by Aubin-Frankowski P C et al. in Bioinformatics in 2020, a method for constructing a gene regulatory network is disclosed. This method uses gene expression data, calculates the association strength between gene expression values and gene expression rates in the form of linear differential equations, and uses the feature selection algorithm TIGRESS to solve the equations to calculate the instantaneous regulatory relationships between genes and construct a gene regulatory network. There are two deficiencies in this method. First, when calculating the instantaneous regulatory relationships between gene pairs, it ignores the regulatory relationships between gene pairs based on time delay. Second, it only calculates the association values between gene expression values and gene expression rates, ignoring the association values between gene expression values, resulting in inaccurate inference of the finally constructed gene regulatory network. Summary of the Invention

[0004] The object of the present invention is to overcome the defects existing in the above-mentioned prior art, and a method for constructing a gene regulatory network based on time delay of single-cell data is proposed to solve the problem of low accuracy in inferring gene regulatory networks in the prior art.

[0005] To achieve the above object, the technical solution adopted by the present invention includes the following steps:

[0006] (1) Calculate the variance of gene expression values based on single-cell data:

[0007] Obtain single-cell data of T single cells, where each single cell contains G genes, and calculate the variance of expression values Val of each gene g , where 10 ≤ T ≤ 1000 and 100 ≤ G ≤ 20000;

[0008] (2) Obtain gene expression data:

[0009] Select N genes with variance of expression values greater than a preset threshold, infer the pseudotime of each single cell, and then sort the T single cells in ascending order of pseudotime to obtain the expression data of each gene under the T single cells. Among them, each gene can be regarded as a regulatory gene or a target gene, and the expression data of the nth gene is

[0010]

[0011] (3) Calculate the maximum mutual information between regulatory genes and target genes based on time delay:

[0012] Initialize M time delays D = {d1, d2,..., d m ,..., d M}, d M < T, and calculate the mutual information m of the expression data between the ith regulatory gene and the jth target gene according to the time delay d of the maximum value with the expression data of the ith regulatory gene and the rate of the jth target gene of the maximum value

[0013] (4) Construct a gene regulatory network:

[0014] Calculate the regulatory strength I between the ith regulatory gene and the jth target gene through the maximum value and calculated in step (3) i,j, and construct a gene regulatory network including N nodes and N×N edges with N genes as nodes and the regulatory strength between two genes as edges, and then remove the edges with regulatory strength less than the preset threshold to obtain the final gene regulatory network.

[0015] Compared with the prior art, the present invention has the following advantages:

[0016] The present invention infers N genes with expression value variances greater than a preset threshold in all single-cell data as the nodes of the gene regulatory network, and then calculates the regulatory strength between the regulatory gene and the target gene by calculating the maximum mutual information between the regulatory gene and the target gene based on time delay. When calculating the regulatory strength, the expression values of the regulatory gene and the target gene, the expression value of the regulatory gene, and the expression rate of the target gene are used for calculation. And the regulatory strength between every two genes is used as the edge of the gene regulatory network, fully considering the regulatory relationship between gene pairs based on time delay, that is, using the mutual information between the expression value of the regulatory gene and the expression rate of the target gene, and the mutual information between the expression value of the regulatory gene and the expression value of the target gene. Compared with the prior art, the problem of low accuracy in inferring the gene regulatory network is effectively improved. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] Figure 1 is a flowchart for the implementation of the present invention. DETAILED DESCRIPTION OF THE INVENTION

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

[0019] Referring to Figure 1 , the present invention includes the following steps:

[0020] Step 1) Calculate the expression value variance of genes based on single-cell data:

[0021] Obtain single-cell data of human embryonic stem cells hESC with T single cells and each single cell containing G genes. This data includes cells measured at 0, 12, 24, 36, 72, and 96 h, and calculate the expression value variance Val of each gene g , where 10≤T≤1000, 100≤G≤20000. In this example, T = 758 and G = 18385;

[0022] The expression value variance Val of each gene g The calculation formula is:

[0023]

[0024] Among them, Val g represents the expression value variance of the g-th gene, represents the expression value of the g-th gene under the t-th single cell, It represents the average expression value of the g-th gene in T single cells.

[0025] Step 2) Obtain gene expression data:

[0026] Select N genes with expression value variances greater than a pre-set threshold. Use the pseudo-trajectory inference method SlingShot to infer the pseudo-time of each single cell with 0h as the starting cluster and 96h as the ending cluster. Then sort the T single cells in ascending order of pseudo-time to obtain the expression data of each gene in the T single cells. Each gene can be regarded as a regulatory gene and a target gene. Among them, the expression data of the n-th gene is In this example, N = 100;

[0027] Step 3) Calculate the maximum mutual information between regulatory genes and target genes based on time delay:

[0028] Initialize M time delays D = {d1, d2,..., d m ,..., d M} where d M < T. The time delay represents the deviation of the starting data point in the expression data of the regulatory gene and the target gene, and the starting time point of the regulatory gene precedes that of the target gene. And calculate the mutual information of the expression data between the i-th regulatory gene and the j-th target gene according to the time delay d m and the maximum value as well as the mutual information between the expression data of the i-th regulatory gene and the rate of the j-th target gene and the maximum value In this example, D = {1, 2, 3, 4, 5, 6, 7, 8, 9, 10}, M = 10;

[0029] Calculate the maximum value of the mutual information of the expression data between the i-th regulatory gene and the j-th target gene and the maximum value of the mutual information between the expression data X i - of the i-th regulatory gene and the rate V j - of the j-th target gene The implementation steps are as follows: The implementation steps are as follows:

[0030] (3a) Initialize M time delays D = {d1, d2,..., d m ,..., d M} and, according to the time delay d m obtain from the expression data of the i-th regulatory gene and the j-th target gene Select the data of C pseudo-time points from it to form the expression data of the i-th regulatory gene The expression data of the j-th target gene Then calculate and The mutual information of 1 ≤ i, j ≤ N, d M < T, C = T - d M ;

[0031] (3b) Calculate the time delay d m below The maximum value of and the rate of the j-th target gene And according to the time delay d m Select the data of E pseudo-time points from V j to form the new rate of the j-th target gene Then calculate the expression data of the new i-th regulatory gene The rate of the j-th target gene The mutual information between The calculation uses the differential approximation of adjacent time points, and then calculates the time delay d m below The maximum value of where 2 ≤ M ≤ 10, E = T - d M -1;

[0032] Calculate and The mutual information of The calculation formula is:

[0033]

[0034]

[0035] where H represents The number of pseudo-times in, ψ(x) represents the Digamma function, k represents the parameter in the k-nearest neighbor estimation method, Represents the joint random variable The corresponding variable z below i =(x i , x j ) Use the k-nearest neighbor estimation method to calculate the distance z i The number of points with a distance not greater than ε i / 2, Represents the point k-nearest to x i , Represents the point k-nearest to x j , ||x|| represents the maximum norm of x, ε iDenote z i Twice the maximum distance under the maximum norm of the k-th nearest neighbor boundary of it.

[0036] Step 4) Construct a gene regulatory network:

[0037] According to the maximum values under M time delays and the maximum values under M time delays Calculate the regulatory strength I between every two genes, and construct a gene regulatory network with N genes as nodes and the regulatory strength between two genes as edges, including N nodes and N×N edges. Then, remove the edges with regulatory strength less than the preset threshold to obtain the final gene regulatory network.

[0038] Calculate the regulatory strength I of the i-th regulatory gene and the j-th target gene i,j , and the calculation formula is:

[0039]

[0040] Next, in combination with simulation experiments, the technical effects of the present invention will be further described:

[0041] 1. Simulation conditions and content:

[0042] The hardware used in the simulation is an AMD Ryzen 7 4800U CPU with a main frequency of 1.6 GHz and 16G of memory; the software is Python 3.7.7.

[0043] Perform a comparative simulation on the AUROC and AUPR of the present invention and the prior art, and the results are shown in Table 1. By comparing the area under the curve AUROC of the Receiver Operating Characteristic ROC curve and the area under the curve AUPR of the Precision-Recall PR curve, the accuracy of the constructed single-cell gene regulatory network is characterized. The higher the values of these two areas, the higher the prediction accuracy of the network. Among them, the abscissa of the ROC curve is the false positive rate, and the ordinate is the true positive rate. The false positive rate is defined as the ratio of the number of misclassified negative samples to the total number of negative samples, and the true positive rate is defined as the ratio of the number of correctly classified positive samples to the total number of positive samples; the abscissa of the PR curve is the recall rate, and the ordinate is the precision rate. The definition of the recall rate is the same as that of the true positive rate, and the precision rate is defined as the ratio of the number of correctly classified positive samples to the total number of samples classified as positive.

[0044] 2. Analysis of simulation results:

[0045] Table 1

[0046] AUROC AUPR Prior art 0.642 0.072 The present invention 0.714 0.139

[0047] As can be seen from Table 1, the present invention is higher than the prior art in both the AUROC value and the AUPR value, which proves that the number of correct edges in the network inferred by the present invention is higher than the number of correct edges inferred by the prior art, effectively improving the accuracy of gene regulatory network inference.

Claims

1. A method for constructing a gene regulatory network with time delay based on single-cell data, characterized in that Including the following steps: (1) Calculate the variance of gene expression values based on single-cell data: Obtain single-cell data of T single cells, where each single cell contains G genes, and calculate the variance Val of the expression values of each gene g , where 10 ≤ T ≤ 1000 and 100 ≤ G ≤ 20000; (2) Obtain gene expression data: Select N genes with expression value variances greater than a pre-set threshold, infer the pseudotime of each single cell, and then sort the T single cells in ascending order of pseudotime to obtain the expression data of each gene under the T single cells. Among them, each gene can be regarded as a regulatory gene and a target gene, and the expression data of the nth gene is (3) Calculate the maximum mutual information between regulatory genes and target genes based on time delay: Initialize M time delays D = {d1, d2,..., d m ,..., d M}, where d M < T, and calculate the maximum mutual information of the expression data between the i-th regulatory gene and the j-th target gene according to the time delay d m as well as the maximum mutual information between the expression data of the i-th regulatory gene and the rate of the j-th target gene . The implementation steps are as follows: ​ (3a) Initialize M time delays D = {d1, d2, ..., d m ,...,d M }, and according to the time delay d m From the expression data of the i-th regulated gene and the j-th target gene Select the data of C pseudo time points from the expression data of the new ith regulatory gene Expression data of the jth target gene Then calculate and Mutual information 1≤i,j≤N,d M <T, C = Td M ; (3b) Calculate the time delay d m lower maximum value of and the rate of the j-th target gene and according to the time delay d m select data of E pseudotime points from V j to form the new rate of the j-th target gene Then calculate the expression data of the new i-th regulatory gene the rate of the j-th target gene mutual information between Recalculate the time delay d m lower maximum value of where, 2 ≤ M ≤ 10, E = T - d M - 1; Calculation and mutual information The calculation formula is as follows: where H represents the number of pseudo-times in, ψ(x) represents the Digamma function, k represents the parameter in the k-nearest neighbor estimation method, represents the joint random variable the corresponding variable z under i =(x i ,x j ) calculates the distance z using the k-nearest neighbor estimation method i the number of points whose distance is no greater than ε i / 2, represents the k-th nearest point to x i , represents the k-th nearest point to x j , ||x|| represents the maximum norm of x, ε i represents twice the maximum value of the distance under the maximum norm from z i to the boundary of its k-th nearest point; (4) Construct a gene regulatory network: The maximum value calculated through step (3) and Calculate the regulatory strength I between the i-th regulatory gene and the j-th target gene i,j , and construct a gene regulatory network with genes as nodes and the regulatory strength between two genes as edges, including N nodes and N×N edges. Then, remove the edges with regulatory strength less than the preset threshold to obtain the final gene regulatory network.

2. The method for constructing a gene regulatory network based on time delay of single-cell data according to claim 1, wherein The variance Val of the expression value of each gene calculated in step (1) g , and the calculation formula is: Among them, represents the expression value of the g-th gene in the t-th single cell, represents the average value of the expression values of the g-th gene in T single cells.

3. The method for constructing a gene regulatory network based on time delay of single-cell data according to claim 1, wherein at the M time delays described in step (3b) the maximum value V j the rate of the j-th target gene in the t-th single cell the expression data of the i-th target gene the rate of the j-th target gene the mutual information I m between, at the M time delays the maximum value The calculation formulas are respectively:[[]] Among them, represents the expression value of the j-th gene in the t-th single cell.

4. The method for constructing a gene regulatory network based on time delay of single-cell data according to claim 1, wherein The calculation of the i-th regulatory gene and the j-th target gene I described in step (4) i,j , and the calculation formula is:

Images