A gene regulatory network inference method based on linear mixed model

Abstract

The present invention discloses a gene regulatory network inference method based on a linear mixed model. The input is single-cell gene expression data, which can be represented by mean, random effect and noise. The random effect and noise are both matrix random variables and obey matrix normal distribution. Then, a row covariance matrix is ​​obtained by parameter estimation to represent the correlation between genes. The present invention greatly improves the computational efficiency and can obtain more reliable, scalable and robust downstream analysis.

Description

Technical Field

[0001] The present invention belongs to the field of gene technology, and in particular relates to a gene regulation network inference method. Background Art

[0002] Gene regulatory networks (GRNs) play a close regulatory role in different tissue types, developmental stages or cell states of biological systems. Therefore, accurate identification of gene interactions is crucial for understanding and analyzing major human diseases or biological processes. For example, transcriptional dysregulation revealed by disease-related gene interactions has been reported in various diseases, including cancer, neurological diseases and psychiatric diseases, which has attracted attention to the functions of corresponding genes in the development of diseases. Single-cell sequencing (scRNA-seq) technology is a high-resolution technology that can elucidate gene expression levels from the resolution of single cells, which makes it possible to construct cell type-specific regulatory networks and analyze changes in gene activity in specific cell types. Single-cell sequencing technology has also brought new topics to computational biology, such as cell differentiation, embryonic development, trajectory inference and immune system response, to help people understand cellular heterogeneity and biological processes.

[0003] In recent years, many methods have been developed to infer gene regulatory networks from single-cell sequencing data. Although they have advantages under certain conditions, they still have some limitations. For example, the ordinary differential equation-based methods SCODE and GRISLI can infer GRNs from data sets with time information to describe the dynamics of gene expression; the regression model-based methods SINCERITIES and GENIE3 aim to find suitable prediction functions to explain the underlying network. The above two types of methods have a common feature, that is, they need to make assumptions about the pattern of gene expression to meet the needs of the algorithm. For example, SCODE assumes that the expression change rate of each transcription factor (TF) depends linearly on their expression profile. Similarly, GENIE3 assumes that the expression level of the target gene is the weighted sum of the expression values ​​of all driver genes. These assumptions are statistically significant, but may not conform to the real biological background. Correlation-based models, such as LEAP, scLink, and PIDC, use correlations calculated by graph models or statistical learning methods to characterize the interactions between genes. The randomness of gene expression between cells and the heterogeneity of tissues make it inevitable to bring errors when constructing networks through correlation coefficients, which may eventually lead to inaccurate networks. Models based on Boolean networks use the numbers "1" or "0" to indicate the state of a node as "on" or "off", such as the SCNS method. Boolean networks can only determine whether there is an edge between genes, but cannot determine the strength of this relationship, which is not conducive to distinguishing the importance of genes in different cell states or biological processes. Summary of the invention

[0004] In order to overcome the shortcomings of the prior art, the present invention provides a gene regulatory network inference method based on a linear mixed model, wherein the input is single-cell gene expression data, which can be represented by mean, random effect and noise, wherein the random effect and noise are both matrix random variables and obey matrix normal distribution; then, the row covariance matrix is ​​obtained by parameter estimation to represent the correlation between genes; the present invention greatly improves the computational efficiency and can obtain more reliable, scalable and robust downstream analysis.

[0005] The technical solution adopted by the present invention to solve the technical problem includes the following steps:

[0006] Step 1: Use a linear mixed model to represent gene expression data:

[0007] Y=M+G+E;G~MN p×n (0,V g , K), E~MN p×n (0,V e , I n×n ), (1)

[0008] Where Y is the p×n dimensional preprocessed gene expression matrix, p is the number of genes, and n is the number of cells; the matrix M represents the mean of Y, and the value M in the matrix M is ij represents the average value of the expression of the i-th gene in each cell, so the elements of each column of M are the same; G is a p×n-dimensional random effect matrix variable that obeys a matrix normal distribution with a mean of 0, and the row covariance matrix is ​​represented by a p×p-dimensional V g , V g represents the correlation between genes. The column covariance matrix is ​​a known n×n dimensional matrix K, which represents the correlation between cell-cell expression and can be calculated through prior information. E is a p×n dimensional noise matrix with a mean of 0. The row covariance matrix is ​​represented by a p×p dimensional V e , the column covariance matrix is ​​the identity matrix;

[0009] Step 2: Subtract the mean from Y, that is, The variance part of Y

[0010] Step 3: Perform eigendecomposition on the matrix K, that is Among them U k is an n×n dimensional orthogonal matrix of eigenvectors, D k is an n×n-dimensional diagonal matrix filled with corresponding eigenvalues;

[0011] Step 4: Transform and vectorize matrices Y, G, and E:

[0012] g=vec(GUk ), e=vec(EU k ), vec means vectorization, that is, connecting each column in sequence into a vector;

[0013] Step 5: The linear mixed model becomes:

[0014]

[0015] Among them, MVN represents multivariate normal distribution, represents the Kronecker product;

[0016] For each cell i, we have:

[0017] y i =g i +e i ; g i ~MVN(0,δ i V g ), e i ~MVN(0,V e ), (3)

[0018] Among them, V g and V e is an unknown parameter;

[0019] Step 6: Derivation of likelihood function; in order to solve the unknown parameter V g and V e , first derive the log-likelihood function, and then estimate it using the PX-EM algorithm;

[0020] Step 6-1: Likelihood function:

[0021]

[0022]

[0023] Where l(y i |g i , V g , V e ) is the likelihood function of the ith cell in y, l(g i |V g , V e ) is the likelihood function of the ith cell of g;

[0024] Step 6-2: Multiply equation (4) and (5) to obtain:

[0025]

[0026] Among them, l(y i , g i |Vg , V e ) is the joint likelihood function;

[0027] Perform log transformation on equation (6):

[0028]

[0029] Step 6-3: Integrate the obtained joint likelihood function of each cell:

[0030]

[0031] Finally, the PX-EM algorithm is used to estimate the parameter V g until convergence condition is reached.

[0032] Preferably, the preprocessing in step 1 is specifically:

[0033] Calculate the variance of the expression counts of each gene in all cells, and then screen 100-500 highly variable genes. After obtaining the gene set through the above operations, standardize the count matrix to obtain CPM, recorded as C; then perform log10 transformation on C to obtain the gene expression matrix

[0034] The beneficial effects of the present invention are as follows:

[0035] 1. The present invention provides additional information to the algorithm by modeling the randomness of gene expression between cells to reduce the adverse effects of cell heterogeneity on the inferred network;

[0036] 2. The present invention uses the row covariance matrix of the random effect term of the model to represent the interaction relationship between genes, avoiding the prior assumption of the expression pattern of the gene;

[0037] 3. The present invention greatly improves computational efficiency by using the PX-EM algorithm for parameter optimization, and can obtain more reliable, scalable and robust downstream analysis;

[0038] 4. Through the analysis of multiple real data sets, the GRNLMM method of the present invention has certain advantages in accurately predicting GRN compared with other popular methods. BRIEF DESCRIPTION OF THE DRAWINGS

[0039] Figure 1 These are illustrations related to the method of the present invention, A: input and model; B: correlation matrix heat map; C: gene regulatory network; D: five functional module diagrams.

[0040] Figure 2This is the result diagram (1) of Example 1 of the present invention, A: ROC curve diagram; B: PRC curve diagram; C: true positive edge number bar graph; D: correlation heat map of five classes; E: functional module diagram of five classes.

[0041] Figure 3 Result diagram (2) of Example 1 of the present invention, A: heat map of gene activity of cells at different stages; B: graph of changes in NANOG, SOX2 and ZFX expression levels over time; C: functional heat map of five class enrichments.

[0042] Figure 4 Result diagram of Example 2 of the present invention, A: ROC curve diagram; B: true positive edge number bar graph; C: correlation heat map of five classes; D: functional module diagram of five classes: E: functional heat map of enrichment of five classes; F: gene activity heat map of cells at different stages.

[0043] Figure 5 Result diagram of Example 3 of the present invention, A: motif sequence comparison of known transcription factors and potential transcription factors of FOS, JUN and TFDP1; B: functional heat map of enrichment of five classes; C: gene activity heat map of cells at different stages; D: correlation heat map of five classes; E: violin plot of expression levels of ATF5 and NFKB1 in tumor and normal cells; F: functional module diagram of five classes. DETAILED DESCRIPTION

[0044] The present invention is further described below in conjunction with the accompanying drawings and embodiments.

[0045] In response to the existing problems, the present invention proposes a gene regulatory network inference algorithm GRNLMM based on a linear mixed model. This method mainly realizes the accurate calculation of the correlation between each pair of genes by inputting scRNA-seq data, which can facilitate downstream analysis, such as clustering, detecting functional modules, and gaining insight into the changes in gene function between two different states.

[0046] like Figure 1 As shown, a gene regulatory network inference method based on a linear mixed model includes the following steps:

[0047] Step 1: Use a linear mixed model to represent gene expression data:

[0048] Y=M+G+E;G~MN p×n (0,V g , K), E~MN p×n (0,V e , I n×n ), (1)

[0049] Where Y is the p×n dimensional preprocessed gene expression matrix, p is the number of genes, and n is the number of cells; the matrix M represents the mean of Y, and the value M in the matrix M is ij represents the average value of the expression of the i-th gene in each cell, so the elements of each column of M are the same; G is a p×n-dimensional random effect matrix variable that obeys a matrix normal distribution with a mean of 0, and the row covariance matrix is ​​represented by a p×p-dimensional V g , V g represents the correlation between genes. The column covariance matrix is ​​a known n×n dimensional matrix K, which represents the correlation between cell-cell expression and can be calculated through prior information. E is a p×n dimensional noise matrix with a mean of 0. The row covariance matrix is ​​represented by a p×p dimensional V e , the column covariance matrix is ​​the identity matrix;

[0050] The above preprocessing means: first calculate the variance of the expression counts of each gene in all cells, and then screen 100-500 highly variable genes (different for each data set). After obtaining the gene set through the above operations, the count matrix is ​​standardized to obtain CPM (count per million), recorded as C. Finally, C is log10 transformed to reduce the impact of some larger observations, and the gene expression matrix is ​​obtained.

[0051] Step 2: Subtract the mean from Y, that is, The variance part of Y

[0052] Step 3: Perform eigendecomposition on the matrix K, that is Among them U k is an n×n dimensional orthogonal matrix of eigenvectors, D k is an n×n-dimensional diagonal matrix filled with corresponding eigenvalues;

[0053] Step 4: Transform and vectorize matrices Y, G, and E:

[0054] g=vec(GU k ), e=vec(EU k ), vec means vectorization, that is, connecting each column in sequence into a vector;

[0055] Step 5: The linear mixed model becomes:

[0056]

[0057] Among them, MVN represents multivariate normal distribution, represents the Kronecker product;

[0058] For each cell i, we have:

[0059] y i =g i +e i ; g i ~MVN(0,δ i V g ), e i ~MVN(0,V e ), (3)

[0060] Among them, V g and V e is an unknown parameter;

[0061] Step 6: Derivation of likelihood function; in order to solve the unknown parameter V g and V e , first derive the log-likelihood function, and then estimate it using the PX-EM algorithm;

[0062] Step 6-1: Likelihood function:

[0063]

[0064]

[0065] Where l(y i |g i , V g , V e ) is the likelihood function of the ith cell in y, l(g i |V g , V e ) is the likelihood function of the ith cell of g;

[0066] Step 6-2: Multiply equation (4) and (5) to obtain:

[0067]

[0068] Among them, l(y i , g i |V g , V e ) is the joint likelihood function;

[0069] Perform log transformation on equation (6):

[0070]

[0071] Step 6-3: Integrate the obtained joint likelihood function of each cell:

[0072]

[0073] Finally, the PX-EM algorithm is used to estimate V g and V e The two parameters are adjusted until convergence conditions are reached. Specific embodiment:

[0075] Example 1: Inferring gene regulatory networks from human embryonic stem cell datasets

[0076] The gene expression data of the final endoderm (DE) cells differentiated from human embryonic stem cells were used, which contained 758 cells (reference: Hirotaka Matsumoto, Hisanori Kiryu, Chikara Furusawa, Minoru SHKo, Shigeru BH Ko, Norio Gouda, Tetsutaro Hayashi, and Itoshi Nikaido. Scode: an efficient regulatory network inference algorithm from single-cell rna-seq during differentiation. Bioinformatics, 33(15): 2314-2321, 2017.). In order to illustrate the performance of GRNLMM in inferring gene regulatory networks, it was compared with three other GRN inference algorithms, GENIE3, SCODE and scLink. The results are shown in Figure 2. Figure 2 As shown in (A) and (B), GRNLMM has the highest auruc value. In addition, the number of real regulatory edges among the top 1000 edges predicted by each method is calculated, as shown in Figure 2 As shown in (C), GRNLMM predicted the largest number of true regulatory relationships. Then, in order to better illustrate that the predicted correlation coefficient has good properties, the correlation coefficient matrix obtained by GRNLMM was clustered and 5 clusters were output (K-means). The heat map of the correlation coefficients between genes in these 5 clusters is shown, as shown in Figure 2 (D). The figure shows that there is a high correlation between the first and fifth categories. NANOG, SOX2, and POU5F1 in these two categories have the highest correlation. Previous studies have pointed out that the POU5F1 gene plays an important role in the development of embryonic stem cells. Then, gene enrichment analysis was performed on these five categories, such as Figure 2 As shown in (E), five gene subsets were extracted based on the correlation between genes. The results showed that their enriched functions were exactly related to embryonic stem cell development, indicating that the correlation coefficient obtained by GRNLMM has good properties for analyzing the functions between genes.

[0077] Subsequently, the dynamics of cellular gene regulation were analyzed based on the time stage of each cell. Figure 3 As shown in (A), a GRN is predicted for each time stage, and then a heat map is sorted according to the calculated gene regulatory activity. The calculation of gene regulatory activity is expressed as follows:

[0078]

[0079] Where p is the number of genes, and the score ranges from 0 to 1. The larger the score, the stronger the regulatory activity of the gene. Figure 3 It can be seen that the regulatory activity of genes such as SP6 and ZFX changes significantly over time, and previous studies have shown that these genes have a certain role in the development of embryonic stem cells. In addition, several genes SOX2, NANOG and ZFX were selected, and the regulatory intensity between them changed significantly over time, so their expression level graphs over time were drawn, as shown in Figure 1. Figure 3 As shown in (B), the expression levels of NANOG and SOX2 decreased over time, while the expression level of ZFX remained almost unchanged over time, which was consistent with the predicted correlation. Finally, a heat map analysis was performed on the enrichment analysis results of the five previously clustered classes to illustrate that the functions enriched by these gene sets are mutually exclusive, such as Figure 3 As shown in (C), the functions enriched in these five classes are basically non-overlapping, indicating that clustering using the correlation matrix predicted by GRNLMM can better distinguish different functional modules.

[0080] Example 2: Capturing the response of mouse dendritic cells stimulated with lipopolysaccharide

[0081] The dataset includes more than 1,700 mouse dendritic cells obtained from bone marrow (reference: Alex K Shalek, Rahul Satija, Joe Shuga, John J Trombetta, Dave Gennert, Diana Lu, Peilin Chen, Rona SGertner, Jellert T Gaublomme, Nir Yosef, et al. Single-cell RNA-seq reveals dynamic paracrine control of cellular variation. Nature, 510(7505): 363–369, 2014.). Specifically, the dataset is divided into cell expression levels measured at 1h, 2h, 4h, and 6h after lipopolysaccharide stimulation. According to the ground truth for evaluation reference, 414 genes (degree greater than or equal to 7) were selected from large to small. First, continue to verify the accuracy of GRNLMM in predicting regulatory relationships. Figure 4 As shown in (A) and (B), GRNLMM has the highest auroc value and the largest number of true regulatory edges (in the top 5000 regulatory strengths) compared to other methods. Subsequently, clustering was performed based on the predicted correlation matrix, and the correlation heat map is shown in Figure 4 (C) Compared with the results in the previous application, the correlation in the figure is relatively small, which may be because the selection of these genes has nothing to do with the correlation between genes, but is related to their degree in the ground truth. Then, gene enrichment analysis was performed on these five classes separately, as shown in Figure 4 As shown in (D), based on the gene regulation intensity and gene enrichment analysis results, a gene subset was obtained. Their functions include regulating T cell differentiation and vesicle fusion, etc., which are consistent with the view that cells secrete cytokines after being stimulated by lipopolysaccharide. Similarly, GRNLMM also achieved good results in distinguishing different functional modules, such as Figure 4 (E) As shown. Finally, the GRNs were predicted for the cells at these four different time stages, and the changes in gene regulatory activity between different stages were analyzed. Figure 4As shown in (F), the regulatory activity obtained in this data set has a clear gradient over time. The activity of some genes such as Gbp2, Gbp6, Cd40 and Ccl22 is significantly upregulated over time. Gbp2 and Gbp6 play a role in many life processes, such as the response of cells to lipopolysaccharide, and the method of the present invention accurately captures this phenomenon. In addition, stimulation of lipopolysaccharide promotes the secretion of various cytokines by cells, which then cause specific immune responses, and Cd40 can enhance antigen binding activity. Ccl22 can participate in many processes, such as the response of cells to cytokine stimulation. These methods can accurately identify the changes in the regulatory intensity of genes at different cell stages.

[0082] Example 3: Discovery of potential transcription factors in diffuse large B-cell lymphoma cells

[0083] The dataset includes 1568 diffuse large B-cell lymphoma (DLBCL) cells and 82 normal cells from the corresponding patients. 100 highly variable genes were selected for subsequent experimental analysis (data can be downloaded from https: / / www.ncbi.nlm.nih.gov / geo / query / acc.cgi?acc=GSE182434).

[0084] First, the gene regulatory networks were predicted for tumor cells and normal cells respectively. Since this dataset does not have ground truth for evaluation, we attempted to analyze the accuracy of the method prediction from the binding site sequences of transcription factors, such as Figure 5 As shown in (A). Previous studies have shown that IRF1 is a transcription factor for FOS, NFKB1 is a transcription factor for JUN, and MYC is a transcription factor for TFDP1. Here, we find the genes predicted by GRNLMM to be most correlated with FOS, JUN, and TFDP1, and compare the binding sites of these three genes with IRF1, NFKB1, and MYC. Figure 5 (B) shows that these three pairs of binding sequences have similar patterns, indicating that STAT1, KLF2 and TCF12 may be potential transcription factors for FOS, JUN and TFDP1, respectively. Then, five clusters were obtained based on the correlation matrix predicted by tumor cells, and their correlation heat maps were drawn. Figure 5As shown in (D), the first and fourth clusters have obvious correlations. Gene enrichment analysis of these five clusters revealed that the two clusters have similar functions. For example, some genes in the first cluster have the function of regulating the response to DNA damage stimuli, and in DLBCL, inhibition of DNA repair can effectively promote tumor cell apoptosis. In addition, some genes in the fourth cluster positively regulate the activity of cyclin-dependent protein serine / threonine kinases, while serine or threonine kinases play an important role in cell cycle regulation and are regulated by some oncogenes. Therefore, the genes contained in these two clusters have special effects on immune response or cell apoptosis after tumorigenesis.

[0085] Similarly, GRNs were predicted for two different types of cells, and the differences in gene regulatory activity between cancer cells and normal cells were analyzed. Figure 5 As shown in (C), the gene with the most significant change in regulatory activity is TSC22D3, which can be used to encode the anti-inflammatory protein glucocorticoid (GC)-induced leucine zipper and play a role in immunosuppression. In addition, the AP-1 complex plays a variety of roles in inflammation and tumor development, and FOSB is a member of the AP-1 family and is involved in the induced transcription of multiple genes containing AP-1 sites. Therefore, the calculated changes in regulatory activity are consistent with its biological process. Next, we searched for differentially expressed genes between normal cells and tumor cells, and analyzed the reasons for their differential expression, such as Figure 5 (E) As shown. Two genes, NFKB1 and ATF5, were found. Abnormal activation of NFKB1 is associated with inflammatory diseases, while sustained inhibition of NFKB1 can lead to improper development of immune cells or delayed cell growth. ATF5 is involved in regulating cell cycle processes and transcriptional regulation. The results showed that the expression of these two genes was low in normal cells, but their expression began to increase in tumor cells. At the same time, it was also found that the correlation between the two genes predicted by GRNLMM also changed in the two cell types. Coincidentally, it was also verified through the transcription factor database that NFKB1 happens to be the transcription factor of ATF5. This shows that the predicted gene correlation is accurate and consistent with the biological background. Finally, the five clustered categories were analyzed for functional modules, as shown in the following figure. Figure 5 As shown in (F), the functions of the five gene subsets obtained after GRNLMM prediction correlation screening are all involved in immune response from different angles. For example, the third cluster can positively regulate IkappaB kinase and NFkappaB signaling. The derivatives of NF-kappaB signaling pathway can lead to debilitating or even fatal inflammation and support some forms of cancer. This shows that GRNLMM can be used to identify corresponding functional modules.

Claims

1. A gene regulatory network inference method based on a linear mixed model, characterized in that: The steps include: Step 1: Use a linear mixed model to represent gene expression data: Y=M+G+E;G~MN p×n (0,V g ,K),E~MN p×n (0,V e ,AND n×n ), (1) Where Y is the p×n dimensional preprocessed gene expression matrix, p is the number of genes, and n is the number of cells; the matrix M represents the mean of Y, and the value M in the matrix M is ij represents the average value of the expression of the i-th gene in each cell, so the elements of each column of M are the same; G is a p×n-dimensional random effect matrix variable that obeys a matrix normal distribution with a mean of 0, and the row covariance matrix is ​​represented by a p×p-dimensional V g , V g represents the correlation between genes. The column covariance matrix is ​​a known n×n dimensional matrix K, which represents the correlation between cell-cell expression and can be calculated through prior information. E is a p×n dimensional noise matrix with a mean of 0. The row covariance matrix is ​​represented by a p×p dimensional V e , the column covariance matrix is ​​the identity matrix; Step 2: Subtract the mean from Y, that is, The variance part of Y Step 3: Perform eigendecomposition on the matrix K, that is Among them U k is an n×n dimensional orthogonal matrix of eigenvectors, D k is an n×n-dimensional diagonal matrix filled with corresponding eigenvalues; Step 4: Transform and vectorize matrices Y, G, and E: g=vec(GU k ), e=vec(EU k ), vec means vectorization, that is, connecting each column in sequence into a vector; Step 5: The linear mixed model becomes: Among them, MVN represents multivariate normal distribution, represents the Kronecker product; For each cell i, we have: y i =g i +e i ;g i ~MVN(0,δ i V g ),e i ~MVN(0,V e ), (3) Among them, V g and V e is an unknown parameter; Step 6: Derivation of likelihood function; in order to solve the unknown parameter V g and V e , first derive the log-likelihood function, and then estimate it using the PX-EM algorithm; Step 6-1: Likelihood function: Where l(y i |g i , V g , V e ) is the likelihood function of the ith cell in y, l(g i |V g , V e ) is the likelihood function of the ith cell of g; Step 6-2: Multiply equation (4) and (5) to obtain: Among them, l(y i , g i |V g , V e ) is the joint likelihood function; Perform log transformation on equation (6): Step 6-3: Integrate the obtained joint likelihood function of each cell: Finally, the PX-EM algorithm is used to estimate the parameter V g until convergence condition is reached.

2. A gene regulatory network inference method based on a linear mixed model according to claim 1, characterized in that: The pre-processing in step 1 is specifically as follows: Calculate the variance of the expression counts of each gene in all cells, and then screen 100-500 highly variable genes. After obtaining the gene set through the above operations, standardize the count matrix to obtain CPM, recorded as C; then perform log10 transformation on C to obtain the gene expression matrix

Images