A method for constructing a gene regulatory network based on meta-analysis

Abstract

The application discloses a method for constructing a gene regulation network based on meta-analysis, and aims to solve the problem that a gene regulation network has a large error and a low accuracy because different research data are not completely homogeneous when the gene regulation network is expanded by combining data of multiple studies to enlarge a sample size, and the method comprises the following steps: performing meta-analysis on each transcriptomic gene expression dataset obtained to obtain a significant differential gene list; obtaining a transcription factor gene list according to a species to which the transcriptomic gene expression data belong, and generating a core gene list according to the transcription factor gene list and the significant differential gene; calculating a Pearson correlation coefficient of each core gene and each significant differential gene, and constructing a co-expression network according to the Pearson correlation coefficient; clustering the co-expression network to obtain a plurality of differential gene modules with high correlation of the core genes; and combining biological knowledge and a structural equation model to construct a corresponding gene regulation network according to each module. The application belongs to the field of gene regulation networks.

Description

TECHNICAL FIELD

[0001] The application relates to a method for constructing a gene regulatory network, in particular to a method for constructing a gene regulatory network based on meta-analysis, and belongs to the field of gene regulatory networks. BACKGROUND

[0002] With the rapid development of high-throughput sequencing technology, a large amount of second-generation sequencing data is generated and stored in public databases. By analyzing these data and constructing a gene regulatory network, it has become a very popular research topic in the field of biology. Due to the complex gene regulatory relationship widely existing in quantitative traits and diseases, by constructing a gene regulatory network, researchers can better determine the specific pathways involved in the disease and develop more effective treatment plans. In theory, the larger the sample size, the better the performance of the gene regulatory network, but due to the limited number of samples measured in a single study, there is often not enough data to accurately recover the gene causal relationship from gene expression levels, and in the process of constructing a gene regulatory network, the size of the sample size directly determines the effect of the gene regulatory network and the feasibility of practical application, but in a single biological study, increasing the sample size is costly, so how to construct an accurate gene regulatory network under the condition of small sample size is a problem that needs to be broken through and solved.

[0003] The common gene regulatory network construction method expands the sample size by simply merging the data of multiple studies, but in the real scene, due to noise and other reasons, the data between different studies is not completely homogeneous, which will lead to the fact that many classic gene regulatory network methods do not work well in the real scene, that is, the constructed gene regulatory network has large errors and low accuracy. SUMMARY

[0004] The present application is to solve the problem that the existing gene regulatory network expands the sample size by merging the data of multiple studies, but in the real scene, due to noise and other reasons, the data between different studies is not completely homogeneous, which leads to the fact that the constructed gene regulatory network has large errors and low accuracy, and further puts forward a method for constructing a gene regulatory network based on meta-analysis.

[0005] The technical scheme adopted by the present application is:

[0006] It comprises the following steps:

[0007] S1, obtaining a plurality of transcriptomic gene expression datasets and their corresponding sample category information in a transcriptome database, the sample category information being normal samples and abnormal samples, each transcriptomic gene expression dataset being represented as a gene expression matrix, each row of the gene expression matrix being a gene, and each column being a sample;

[0008] S2, meta-analysis is performed on each gene expression matrix to obtain a list of significantly different genes across datasets:

[0009] S21, the effect value of each gene in each gene expression matrix is calculated;

[0010] S22, according to the effect value of each gene in each gene expression matrix, the combined effect value of each gene is calculated by using the fixed effect model, the screening parameter and the screening standard of the gene are defined, the hypothesis test is performed on each gene, the corresponding screening parameter value is obtained, the gene is screened according to the screening standard, the screened gene and the corresponding multiple gene expression matrices of the screened gene are obtained;

[0011] S23, the multiple gene expression matrices obtained in S22 are merged according to the column to obtain one matrix, the self-help sampling of the samples in the merged matrix is performed to obtain bootstrap samples, the test statistic D+ is calculated according to the bootstrap samples, the screening parameter value of each gene is obtained according to the test statistic D+, and the multiple significantly different genes across datasets are obtained by using the screening standard;

[0012] S3, the corresponding transcription factor gene list is obtained according to the species to which the transcriptomic gene expression data belongs, and the core gene list is generated according to the transcription factor gene list and the significantly different genes;

[0013] S4, the Pearson correlation coefficient between each core gene in the core gene list and each significantly different gene is calculated according to the bootstrap samples obtained in S23, and the co-expression network is constructed according to the Pearson correlation coefficient;

[0014] S5, the co-expression network obtained in S4 is clustered to obtain multiple core gene high correlation difference gene modules;

[0015] S6, according to each core gene high correlation difference gene module, the corresponding gene regulation network is constructed by combining biological knowledge and structural equation model, that is, each core gene high correlation difference gene module corresponds to a gene regulation network.

[0016] Further, in S1, all normal samples are taken as the control group and all abnormal samples are taken as the experimental group.

[0017] Further, the specific process of S21 is:

[0018] S211, the effect value of a certain gene in the ith gene expression matrix is calculated:

[0019]

[0020] Wherein, T i represents the effect value, c(m) represents the factor for correcting positive bias, nE and n c respectively represent the sample number of the experimental group and the control group, and respectively represent the gene expression of the experimental group and the control group, and s represents the combined standard deviation of the experimental group and the control group, and respectively represent the variance of the experimental group and the control group;

[0021] S212, calculate the effect value variance of the gene in S211:

[0022]

[0023] wherein, v(T i ) represents the effect value variance;

[0024] Repeat the above process to obtain the effect value and variance of each gene in each gene expression matrix.

[0025] Further, the S22 specific process is:

[0026] S221, calculate the combined effect value of each gene:

[0027]

[0028]

[0029]

[0030] wherein, is the combined effect value, is the combined effect value variance, w i is the weight given to the effect value T i of each gene expression matrix, specifically the inverse of the variance of the effect value;

[0031] S222, define the screening parameter of the gene as the significance P value or the false discovery rate FDR, when the P value is used as the screening parameter, the screening standard is P<0.01 or P<0.05 or P<0.1; when the FDR is used as the screening parameter, the screening standard is FDR<0.01 or FDR<0.05 or FDR<0.1;

[0032] S223, convert each combined effect value into a Z statistic that meets the standard normal distribution, and then perform a two-tailed hypothesis test to obtain the corresponding screening parameter value:

[0033]

[0034] S224, screen the genes according to the screening criteria based on the screening parameter value, to obtain screened genes and a plurality of gene expression matrices corresponding to the screened genes.

[0035] Further, the S23 specific process is:

[0036] S231, merge the plurality of gene expression matrices obtained after S224 according to columns into one matrix, take each column in the merged matrix as a sampling unit to perform bootstrap sampling, define the number of bootstrap samplings and the sample size of each bootstrap sampling, obtain a plurality of sampling units through each bootstrap sampling, merge the plurality of sampling units according to columns into one gene expression matrix, take the gene expression matrix as a bootstrap sample corresponding to the bootstrap sampling round, and obtain a plurality of bootstrap samples;

[0037] It is assumed that the gene space contained in each sampling unit is consistent, that is, each sampling unit has the same gene expression data, so the same sampling unit repeatedly sampled in the above bootstrap sampling is taken as a new sampling unit, and the gene space represents a set of genes in a transcriptome data set or a set of rows of a gene expression matrix;

[0038] S232, calculate the weight a of each gene in each bootstrap sample:

[0039]

[0040] wherein a i is the weight of the i-th gene, y is the sample category information of the sampling unit, that is, the experimental group or the control group, F(i, y) is the F score between the i-th gene and y, and C(i, j) is the Pearson correlation coefficient between the i-th gene and the j-th gene in the corresponding transcriptome gene expression data set;

[0041] In each bootstrap sample, the genes are sorted according to the weight a of each gene from large to small to obtain a sorted gene list, and the above process is performed on all bootstrap samples respectively to obtain a plurality of gene lists, and the specific process is:

[0042]

[0043] wherein, is the ranking score of the i-th gene in the b-th gene list, p ib is the ranking of the i-th gene in the b-th gene list, and N is the number of genes in the b-th gene list;

[0044]

[0045] wherein, is used to calculate the deviation of the ranking score of the i-th gene in the b-th bootstrap sample from Q2, and Q2 is a constant;

[0046] S233, after b times of self-sampling, the following hypothesis test is performed on each gene in the gene space:

[0047] H0 hypothesis: the ith gene has no information on the given condition / characteristics, i.e. A i ≤ 0, A i is the median deviation expected ranking score of the ith gene in the bootstrap sample that accepts the hypothesis test, i.e.

[0048] H1 hypothesis: the ith gene has information on the given condition / characteristics, i.e. A i > 0;

[0049] After b times of self-sampling, b r b are obtained for each gene, the b r b are sorted in ascending order according to the size of r b , and the sorted r b list of each gene is obtained, i.e. r b is assigned a ranking of 1, 2,..., l, and it is remembered whether r b is positive or not, and the test statistic D+ of each gene is defined as the sum of the rankings of r b > 0; b .

[0050] S234, define an independent Bernoulli variable B, B is 1 when r b > 0; b , and B is 0 when r b < 0, then the mean and variance of the test statistic D+ are:

[0051]

[0052]

[0053] When A i = 0 under the H0 hypothesis, E(B) is , then the mean of the test statistic D+ is , and the variance is When the number of self-sampling is 200, D+ follows a standard normal distribution.

[0054] S235, perform hypothesis testing on the test statistic D+ of each gene with the standard normal distribution to obtain the significance P value of each gene, and then use the screening standard to obtain multiple significant difference genes across datasets. Combine all significant difference genes into a list to obtain a list of significant difference genes.

[0055] Further, the specific formula of the F score in S232 is:

[0056]

[0057] wherein, is the mean of expression level of the ith gene in the experimental group samples, is the mean of expression level of the ith gene in all samples, is the mean of expression level of the ith gene in the control group samples, is the variance of expression level of the ith gene in the experimental group samples, is the variance of expression level of the ith gene in the control group samples.

[0058] Further, the S3 specific process is:

[0059] S31, obtaining the corresponding transcription factor gene list from the transcription factor database according to the species to which the transcriptomic gene expression data belongs;

[0060] S32, extracting the same genes in the significant differential gene list obtained in S235 and the transcription factor gene list in S31 to generate a core gene list.

[0061] Further, the S4 specific process is:

[0062] performing Pearson correlation analysis between the core genes in the core gene list and the significant differential genes two by two to obtain a plurality of Pearson correlation coefficients, then obtaining the correlation value and the significance value between the expression level of each core gene and each significant differential gene, obtaining the co-expression relationship pair between the core gene and the significant differential gene according to the correlation value and the significance value, constructing a co-expression network according to the Pearson correlation coefficient and the co-expression relationship pair, the nodes of the co-expression network being all the core genes and the significant differential genes having the co-expression relationship pair with the core genes, and the edges being the co-expression relationship pairs between the core genes and the significant differential genes.

[0063] Further, the specific process of obtaining the co-expression relationship pair between the core gene and the significant differential gene according to the correlation value and the significance value is:

[0064] (1) when the absolute value of the correlation value of the core gene and the significant differential gene is ≥0.9 and the significance value is <0.01, the core gene and the significant differential gene are a co-expression relationship pair;

[0065] (2) when the absolute value of the correlation value of the core gene and the significant differential gene is ≥0.8 and the significance value is <0.01, the core gene and the significant differential gene are a co-expression relationship pair;

[0066] (3) When the absolute value of the correlation value of the core gene and the significantly different gene is ≥ 0.7, and the significance value is <0.01, the core gene and the significantly different gene are a co-expression relationship pair;

[0067] (4) When the absolute value of the correlation value of the core gene and the significantly different gene is ≥ 0.9, and the significance value is <0.05, the core gene and the significantly different gene are a co-expression relationship pair;

[0068] (5) When the absolute value of the correlation value of the core gene and the significantly different gene is ≥ 0.8, and the significance value is <0.05, the core gene and the significantly different gene are a co-expression relationship pair;

[0069] When the absolute value of the correlation value of the core gene and the significantly different gene is ≥ 0.7, and the significance value is <0.05, the core gene and the significantly different gene are a co-expression relationship pair.

[0070] Further, the S6 specific process is:

[0071] S61, using the pathways function of the R package graphtie, obtaining a plurality of known gene regulation relationships from a pathway database;

[0072] Using the R package igraph, converting the plurality of gene regulation relationships into a network form to obtain an a priori regulation network, one node in the a priori regulation network representing one gene, and one edge representing a regulation relationship between two genes, each edge being taken as one regulation relationship pair;

[0073] S62, according to the shortest path standard, searching for a regulation relationship pair between any two genes in each core gene high correlation difference gene module in the a priori regulation network, and screening the regulation relationship pairs in each core gene high correlation difference gene module according to the path length to obtain all sample-specific directed acyclic graphs corresponding to each core gene high correlation difference gene module in S1;

[0074] When the shortest path length = 1 or 2, the two genes are a regulation relationship pair;

[0075] S63, using the structural equation model and the bootstrap sample obtained in S23 to perform causal reasoning on the directed acyclic graph obtained in S62 to generate a gene regulation network corresponding to each core gene high correlation difference gene module:

[0076] Convert the directed acyclic graph into a structural equation model, calculate the significance and causal effect of the regulatory relationship pairs in the gene expression data according to the gene expression matrix obtained by S231, when the significance value of each regulatory relationship pair is <0.05, it is the final regulatory relationship pair, retain the regulatory relationship pair, otherwise, delete each regulatory relationship pair, generate a gene regulatory network for the retained regulatory relationship pair, repeat the above process to obtain the gene regulatory network corresponding to each core gene highly correlated difference gene module.

[0077] Beneficial effects:

[0078] The present application takes each transcriptomic gene expression dataset obtained in the transcriptome database as a gene expression matrix, and uses all existing data as sample size, which makes up for the problem of insufficient sample size when constructing a gene regulatory network. Meta-analysis is performed on the gene expression matrix to screen out significant difference genes across datasets from tens of thousands of genes, taking into account the heterogeneity between research data, which can better handle more and more genomic data under high-throughput sequencing technology. Specifically, the effect value of each gene in each gene expression matrix is calculated, as well as the combined effect value of each gene in all gene expression matrices, and the screening parameters and screening criteria of the gene are defined to screen the gene, obtaining multiple gene expression matrices corresponding to the screened genes, merging multiple gene expression matrices into one matrix, and performing bootstrap sampling (i.e. resampling) on the samples in the merged matrix to obtain bootstrap samples. According to the bootstrap samples, the test statistic D+ is calculated, the screening parameter value of each gene is obtained, and the screening criteria are used to obtain multiple significant difference genes across datasets. Next, the corresponding transcription factor gene list is obtained from the transcription factor database according to the species to which the transcriptomic gene expression data belongs, the same genes in the significant difference gene list and the transcription factor gene list are extracted, and a core gene list is generated. According to the bootstrap samples, the Pearson correlation coefficient of each core gene in the core gene list and each significant difference gene is calculated, and a co-expression network is constructed according to the Pearson correlation coefficient. The nodes of the co-expression network are all core genes and significant difference genes that have a co-expression relationship with the core genes, and the edges are the co-expression relationship pairs of the core genes and the significant difference genes. Cluster the co-expression network to obtain multiple core gene highly correlated difference gene modules. Combined with biological knowledge and structural equation model, construct the corresponding gene regulatory network according to each core gene highly correlated difference gene module, that is, each core gene highly correlated difference gene module corresponds to a gene regulatory network, which improves the accuracy and sensitivity of the gene regulatory network.

[0079] Compared with a conventional gene regulation network construction method, the application has better performance, is also beneficial to indirectly mine key disease or stress resistance genes, and reveals a complex regulation network related to the genes, thereby providing more candidate genes for wet experiments and effectively promoting development of a disease treatment plan and genetic improvement of crop resistance. BRIEF DESCRIPTION OF DRAWINGS

[0080] Figure 1 is a flowchart of the application;

[0081] Figure 2 is a co-expression network diagram of the embodiment;

[0082] Figure 3 is a gene regulation network diagram of the module of the embodiment; DETAILED DESCRIPTION

[0083] Detailed implementation one: combination Figures 1-3 In the embodiment, the method for constructing a gene regulation network based on meta-analysis comprises the following steps:

[0084] S1, obtaining a plurality of transcriptomic gene expression data sets and corresponding sample category information thereof in a transcriptome database, the sample category information being normal samples and abnormal samples, taking all normal samples as a control group and all abnormal samples as an experimental group, each transcriptomic gene expression data set being expressed as a gene expression matrix, each row of the gene expression matrix being a gene and each column being a sample.

[0085] The transcriptome database is NCBI, EBI, TCGA or the like. The transcriptome is RNA. According to the central dogma, the relationship between a gene and a transcriptome is that the gene guides synthesis of the transcriptome. Expression data of the gene and corresponding sample category information thereof are obtained in any one or more of the above transcriptome databases, the expression data being expressed as a plurality of gene expression matrices after quantification, wherein the normal samples are derived from normal tissues of a disease patient or tissues of a plant under normal development conditions, and the abnormal samples are derived from lesion tissues of the disease patient or tissues of the plant under stress conditions. The distribution of the normal samples and the abnormal samples in each gene expression matrix is disordered.

[0086] For example, the application obtains four transcriptomic gene expression data sets of Alzheimer's disease and corresponding sample category information thereof from the NCBI database. Each transcriptomic gene expression data set contains hippocampal tissues (the experimental group) from a patient with Alzheimer's disease after death and hippocampal tissues (the control group) from a normal person after death, and the positions of the normal tissues and the lesion tissues are consistent.

[0087] S2, meta-analysis is performed on each gene expression matrix to obtain a list of significantly different genes across datasets. Meta-analysis is a good statistical method that combines information from multiple studies and considers heterogeneity between studies to improve the power, stability and accuracy of the results. Therefore, a gene regulatory network can be constructed based on meta-analysis. The specific process is as follows:

[0088] S21, calculate the effect value of each gene in each gene expression matrix:

[0089] S211, calculate the effect value of a certain gene in the ith gene expression matrix:

[0090]

[0091] where T i represents the effect value, c(m) represents a factor for correcting positive bias, n E and n c represent the sample size of the experimental group and the control group, and represent the gene expression of the experimental group and the control group, and s represents the combined standard deviation of the experimental group and the control group, and represent the variance of the experimental group and the control group.

[0092] S212, calculate the variance of the effect value of the gene in S211:

[0093]

[0094] where v(T i ) represents the variance of the effect value.

[0095] Repeat the above process to obtain the effect value and variance of each gene in each gene expression matrix.

[0096] S22, according to the effect value of each gene in each gene expression matrix, calculate the combined effect value of each gene using a fixed effect model, define the screening parameter and screening standard of the gene, perform hypothesis testing on each gene, obtain the corresponding screening parameter value, screen the gene according to the screening standard, obtain the screened gene, and the screened gene corresponds to multiple gene expression matrices. The specific process is as follows:

[0097] S221, calculate the combined effect value of each gene:

[0098]

[0099]

[0100]

[0101] wherein, is the combined effect value, is the variance of the combined effect value, w i is the weight given to each gene expression matrix in the effect value T i , specifically the inverse of the variance of the effect value.

[0102] S222, define the screening parameter of the gene as the significance P value or the false discovery rate FDR (false discovery rate), when the P value is used as the screening parameter, the screening standard is P < 0.01 or P < 0.05 or P < 0.1; when the FDR is used as the screening parameter, the screening standard is FDR < 0.01 or FDR < 0.05 or FDR < 0.1. The standard threshold of the above screening parameter can be adjusted according to the characteristics of the transcriptome gene expression dataset and the requirements and purposes of the regulatory network analysis.

[0103] For example, when used for Alzheimer's disease, the screening standard of differential expression with P value as the screening parameter is P < 0.05, so that the accuracy of the regulatory network obtained finally can be improved.

[0104] S223, convert each combined effect value into a Z statistic value conforming to the standard normal distribution, and then perform a two-tailed hypothesis test to obtain the corresponding screening parameter value.

[0105]

[0106] P = 2 (1 (φ (| Z |)))

[0107] S224, based on the screening parameter value, screen the genes according to the screening standard to obtain screened genes and the multiple gene expression matrices corresponding to the screened genes.

[0108] S23, merge the multiple gene expression matrices obtained in S22 according to the columns into one matrix, and perform self-sampling (i.e. replacement resampling) on the experimental group samples and the control group samples in the merged matrix to obtain bootstrap samples, calculate the test statistic D+ according to the bootstrap samples, obtain the screening parameter value of each gene according to the test statistic D+, and obtain multiple significant differential genes across datasets by using the screening standard, the specific process being:

[0109] S231, merge the multiple gene expression matrices obtained after S224 according to columns into one matrix, take each column in the merged matrix as a sampling unit to perform bootstrap sampling, define the number of bootstrap sampling and the sample size of each bootstrap sampling, obtain multiple sampling units through each bootstrap sampling, merge the multiple sampling units according to columns into one gene expression matrix, take the gene expression matrix as a bootstrap sample corresponding to the bootstrap sampling round, and obtain multiple bootstrap samples, and the sample size of each bootstrap sampling needs to be less than the total sample size.

[0110] The application sets the gene space contained in each sampling unit (sample) to be consistent, that is, each sampling unit has the same gene expression data, so the same sampling unit repeatedly sampled by the bootstrap sampling is taken as a new sampling unit, and the gene space represents a set of genes in a transcriptome data set or a set of rows of a gene expression matrix.

[0111] For example, when the application is used for Alzheimer's disease, the total number of bootstrap samplings is 200, and the sample size of each bootstrap sampling is 90% of the original sample size.

[0112] S232, in each bootstrap sample, a sorted gene list is obtained by using a standard maximum correlation minimum redundancy algorithm (MRMR algorithm), and the specific process is as follows:

[0113] The weight a of each gene in each bootstrap sample is calculated:

[0114]

[0115] wherein, a i is the weight of the i-th gene, y is the sample category information of the sampling unit, that is, the experimental group or the control group, F(i, y) is the F score between the i-th gene and y, and C(i, j) is the Pearson correlation coefficient between the i-th gene and the j-th gene in the corresponding transcriptome gene expression data set. In this way, the weight a is used as the standard for gene ranking.

[0116]

[0117] wherein, is the mean of the expression level of the i-th gene in the experimental group sample, is the mean of the expression level of the i-th gene in all samples, is the mean of the expression level of the i-th gene in the control group sample, is the variance of the expression level of the i-th gene in the experimental group sample, is the variance of the expression level of the i-th gene in the control group sample. The application calculates F(i, y) and the weight a of each gene in each bootstrap sample.

[0118] In each bootstrap sample, genes are sorted according to the weight a of each gene from large to small, and a sorted gene list is obtained. According to the above process, all bootstrap samples are sorted respectively, and multiple gene lists are obtained. The specific process is as follows:

[0119]

[0120] wherein, is the ranking score of the ith gene in the bth gene list, p ib is the ranking of the ith gene in the bth gene list, and N is the total number of genes in Ω, that is, the number of genes in the bth gene list.

[0121] It can be noted that for the ith gene, is a random variable. Therefore, without loss of generality, another random variable can be defined as:

[0122]

[0123] wherein, is used to calculate the deviation of the ranking score of the ith gene in the bth bootstrap sample from Q2, Q2 is the second quartile (median) of the ranking score, and is a constant. Because the ranking score is a function of the gene ranking p ib , the ranking score of the gene in each gene list is symmetrically distributed about Q2 (Q2 = 0.5).

[0124] S233, in order to select information genes, after b times of bootstrap sampling, the following hypothesis test needs to be continuously performed on each gene in the gene space:

[0125] H0 hypothesis: the ith gene has no information on the given condition / character, that is, A i ≤ 0, A i is the median deviation expected ranking score of the ith gene in the bootstrap sample that accepts the hypothesis test, that is,

[0126] H1 hypothesis: the ith gene has information on the given condition / character, that is, A i > 0.

[0127] After b times of bootstrap sampling, b r b of each gene are obtained. According to the size of r b , the b r b are arranged in ascending order, and the sorted r b list of each gene is obtained, that is, r b is assigned a ranking of 1, 2,..., l, and it is remembered that r bwhether it is positive, define the test statistic D+ of each gene as the sum of the ranks of the corresponding genes, r b > 0, r b b b If the gene is in the 5th, 7th, 8th, and 15th positions in the list, then the test statistic D+ is 5 + 7 + 8 + 15 = 35.

[0128] S234, in order to calculate the mean and variance of the test statistic D+, define an independent Bernoulli variable B, B is 1 when r b > 0, and 0 when r b < 0, then the mean and variance of the test statistic D+ are:

[0129]

[0130]

[0131] Under the assumption of H0, A i = 0, E(B) is The mean of the test statistic D+ is The variance is When the number of bootstrap samples is large enough (such as b = 200), D+ follows a standard normal distribution.

[0132] S235, perform hypothesis testing on each gene's test statistic D+ and the standard normal distribution to obtain the significance P value of each gene, and then use the screening standard to obtain multiple significant difference genes across datasets. Combine all significant difference genes across datasets into a list to obtain a list of significant difference genes across datasets.

[0133] S3, according to the species to which the transcriptomic gene expression data belongs, obtain the corresponding transcription factor gene list from the transcription factor database, extract the same genes from the significant difference gene list obtained in S235 and the transcription factor gene list, and generate a core gene list.

[0134] S4, according to the gene expression matrix (bootstrap sample) obtained in S231, calculate the Pearson correlation coefficient between each core gene in the core gene list and each significant difference gene, and construct a co-expression network according to the Pearson correlation coefficient. The nodes of the co-expression network are all core genes and significant difference genes that have a co-expression relationship with the core genes, and the edges are the co-expression relationship pairs of the core genes and the significant difference genes. The specific process is as follows:

[0135] ​​performing Pearson correlation analysis between each core gene in the core gene list and each significantly different gene, obtaining a plurality of Pearson correlation coefficients, and then obtaining a correlation value and a significance value between the expression level of each core gene and each significantly different gene, and obtaining a co-expression relationship pair between the core gene and the significantly different gene according to the correlation value and the significance value, and constructing a co-expression network according to the plurality of Pearson correlation coefficients and the co-expression relationship pair:

[0136] (6) when the absolute value of the correlation value of the core gene and the significantly different gene is greater than or equal to 0.9 and the significance value is less than 0.01, the core gene and the significantly different gene are a co-expression relationship pair;

[0137] (7) when the absolute value of the correlation value of the core gene and the significantly different gene is greater than or equal to 0.8 and the significance value is less than 0.01, the core gene and the significantly different gene are a co-expression relationship pair;

[0138] (8) when the absolute value of the correlation value of the core gene and the significantly different gene is greater than or equal to 0.7 and the significance value is less than 0.01, the core gene and the significantly different gene are a co-expression relationship pair;

[0139] (9) when the absolute value of the correlation value of the core gene and the significantly different gene is greater than or equal to 0.9 and the significance value is less than 0.05, the core gene and the significantly different gene are a co-expression relationship pair;

[0140] (10) when the absolute value of the correlation value of the core gene and the significantly different gene is greater than or equal to 0.8 and the significance value is less than 0.05, the core gene and the significantly different gene are a co-expression relationship pair;

[0141] (11) when the absolute value of the correlation value of the core gene and the significantly different gene is greater than or equal to 0.7 and the significance value is less than 0.05, the core gene and the significantly different gene are a co-expression relationship pair;

[0142] The Pearson correlation analysis can be performed by using the cor.test function of R software. The screening threshold of the correlation value and the significance value can be adjusted according to the type of the disease, the requirement and the purpose of the regulatory network analysis. When used for Alzheimer's disease, when the absolute value of the correlation value of the relationship pair is greater than 0.8 and the significance value is less than 0.01, the relationship pair is a co-expression relationship pair, which can improve the accuracy and sensitivity of the final obtained regulatory network. Specifically, the co-expression network is composed of nodes and undirected edges between the nodes, wherein each node represents a gene, and each edge represents a co-expression relationship pair between each core gene and each significantly different gene.

[0143] S5, clustering the co-expression network obtained in S4 to obtain a plurality of core gene high correlation difference gene modules.

[0144] The core gene high correlation differential gene module is composed of core genes and significant differential genes in a co-expression relationship with the core genes. It is clustered from the co-expression network by the louvain algorithm. For example, when the present application is used for Alzheimer's disease, eight core gene high correlation differential gene modules are clustered.

[0145] S6, combine existing biological knowledge with structural equation modeling to construct a gene regulation network from gene expression data, and construct a corresponding gene regulation network according to each core gene high correlation differential gene module, i.e. each core gene high correlation differential gene module corresponds to a gene regulation network. The specific process is as follows:

[0146] S61, obtain prior knowledge of gene regulation relationships from a pathway database and convert it into a network form to obtain a prior regulatory network, wherein a node in the prior regulatory network is a gene and an edge between nodes is a causal relationship between genes. The specific process is as follows:

[0147] S611, use the pathways function of the R package graphtie to obtain multiple known gene regulation relationships from a pathway database (such as the KEGG database, WikiPathways database, Reactome database).

[0148] S612, use the R package igraph to convert multiple gene regulation relationships into a network form to obtain a prior regulatory network, wherein a node in the prior regulatory network represents a gene and an edge represents a regulation relationship between two genes. Each edge is taken as a regulation relationship pair.

[0149] S62, according to the shortest path standard, search for a regulation relationship pair between any two genes in each core gene high correlation differential gene module in the prior regulatory network, and filter the regulation relationship pairs in each core gene high correlation differential gene module according to the path length to obtain a directed acyclic graph specific to all samples in S1 corresponding to each core gene high correlation differential gene module.

[0150] Specifically, in the prior gene regulation network, a path is composed of a plurality of directed edges. A directed edge connecting two nodes can be regarded as a regulation relationship pair. When searching for the shortest path, one gene is taken as the starting point and another gene is taken as the end point. All paths between the two genes are searched in the prior gene regulation network and the length of each path (i.e. the number of directed edges included in the path) is calculated. The regulation relationship pair between the two genes is filtered according to the shortest path length between the two genes. The standard value of the shortest path can be adjusted according to the different types of diseases, the requirements and purposes of the regulation network analysis. When the shortest path length = 1 or 2, the two genes are a regulation relationship pair.

[0151] For example, when the present application is used for Alzheimer's disease, based on the analysis of a large number of complex disease transcriptome data, it is found that when the shortest path length is 1 or 2, the relationship between the genes is a causal relationship, which can improve the accuracy and sensitivity of the finally obtained regulatory network.

[0152] For example, when the present application is used for Alzheimer's disease, based on the analysis of a large number of complex disease transcriptome data, it is found that when the shortest path length is 1 or 2, the relationship between the genes is a causal relationship, which can improve the accuracy and sensitivity of the finally obtained regulatory network.

[0153] S63, using the structural equation model and the gene expression matrix (bootstrap sample) obtained in S231, the directed acyclic graph obtained in S62 is causally inferred, and the gene regulatory network corresponding to each core gene high correlation difference gene module is generated, and the specific process is:

[0154] The directed acyclic graph is converted into a structural equation model, and the significance and causal effect of the regulatory relationship pairs in the gene expression data are calculated according to the gene expression matrix obtained in S231, and when the significance value of each regulatory relationship pair is <0.05, it is the final regulatory relationship pair, the regulatory relationship pair is retained, otherwise, each regulatory relationship pair is deleted, the retained regulatory relationship pair is generated into a gene regulatory network, and the above process is repeatedly executed to obtain the gene regulatory network corresponding to each core gene high correlation difference gene module. The required gene regulatory network is obtained according to the actual situation.

[0155] Specifically, the significance and causal effect of the regulatory relationship pairs and genes in the gene expression data in S63 are calculated, and the structural equation model is used for calculation, and the SEMrun function of the SEMgraph R package is used for running.

[0156] The present application uses meta-analysis to construct a gene regulatory network from multiple transcriptome gene expression matrices, and automatically executes S2-S5. It makes up for the small sample size and long experimental period of traditional wet experiments, and can better process more and more transcriptome data under high-throughput sequencing technology. It is also beneficial to indirectly mine disease or stress-resistant genes, and reveal the complex regulatory network related to these genes, which indirectly promotes the development of disease treatment programs and crop resistance genetic improvement.

[0157] Embodiment:

[0158] This embodiment uses four people's transcriptome data sets to verify the disease-related key gene regulatory network construction method provided by the present application. The four data sets can be divided into control group and experimental group, and the control group is the brain tissue of healthy people, and the experimental group is the brain tissue of Alzheimer's disease patients. The flow chart of the method for constructing the gene regulatory network in the embodiment is as follows Figure 1As shown. First, 4 transcriptomic gene expression datasets were obtained from the transcriptome database, and a fixed effect model and a maximum relevance minimum redundancy feature selection algorithm were used to obtain 1106 significant differentially expressed genes with FDR < 0.05 as the screening standard. The human transcription factor list was obtained from the AnimalTFDB database, and 54 differentially expressed transcription factors were obtained. The differentially expressed transcription factors and the expression level of the genes were analyzed by Pearson correlation analysis, and the co-expression relationship pairs with an absolute value of the correlation value greater than 0.8 and a significance value less than 0.01 were screened, and 1206 gene-gene co-expression relationship pairs were obtained. The co-expression network as shown in Figure 2 was constructed. The co-expression network was clustered to obtain 8 Alzheimer's disease related modules. The largest module (module 8) was selected to construct the gene regulation network, that is, the causal relationship between the genes in module 8 with a shortest path less than or equal to 2 was obtained from the Reactome database and integrated to obtain a directed acyclic graph. The directed acyclic graph was used to perform causal reasoning using a structural equation model to obtain the module 8 gene regulation network in Figure 3 . Figure 3 The dashed line in the above figure represents a significantly down-regulated regulatory relationship, and the solid line represents a significantly up-regulated regulatory relationship.

Claims

1. A method for constructing a gene regulatory network based on meta-analysis, characterized in that: It comprises the following steps: S1, obtaining a plurality of transcriptomic gene expression datasets and their corresponding sample category information in a transcriptome database, the sample category information being normal samples and abnormal samples, each transcriptomic gene expression dataset being expressed as a gene expression matrix, each row of the gene expression matrix being a gene, and each column being a sample; S2, performing meta-analysis on each gene expression matrix to obtain a list of significantly different genes across datasets: S21, calculating the effect value of each gene in each gene expression matrix; S22, calculating the combined effect value of each gene using a fixed effect model according to the effect value of each gene in each gene expression matrix, defining a screening parameter and a screening criterion for the gene, performing hypothesis testing on each gene to obtain a corresponding screening parameter value, screening the gene according to the screening criterion to obtain a screened gene, and obtaining a plurality of gene expression matrices corresponding to the screened gene; S23, merging the plurality of gene expression matrices obtained in S22 according to the columns into one matrix, performing bootstrap sampling on the samples in the merged matrix to obtain bootstrap samples, calculating a test statistic D+ according to the bootstrap samples, obtaining a screening parameter value for each gene according to the test statistic D+, and obtaining a plurality of significantly different genes across datasets using the screening criterion; S3, obtaining a corresponding list of transcription factor genes according to the species to which the transcriptomic gene expression data belongs, and generating a core gene list according to the list of transcription factor genes and the significantly different genes; S4, calculating the Pearson correlation coefficient between each core gene in the core gene list and each significantly different gene according to the bootstrap samples obtained in S23, and constructing a co-expression network according to the Pearson correlation coefficient, the specific process being: performing Pearson correlation analysis between the core genes in the core gene list and the significantly different genes to obtain a plurality of Pearson correlation coefficients, then obtaining the correlation value and significance value between each core gene and each significantly different gene expression level, obtaining the co-expression relationship pair between the core gene and the significantly different gene according to the correlation value and the significance value, and constructing a co-expression network according to the Pearson correlation coefficient and the co-expression relationship pair, the nodes of the co-expression network being all core genes and significantly different genes having a co-expression relationship pair with the core genes, and the edges being the co-expression relationship pairs between the core genes and the significantly different genes; S5, clustering the co-expression network obtained in S4 to obtain a plurality of core gene highly correlated difference gene modules; S6, combining biological knowledge and structural equation models to construct a corresponding gene regulation network according to each core gene highly correlated difference gene module, i.e., each core gene highly correlated difference gene module corresponds to a gene regulation network, the specific process being: S61, using the pathways function of the graphtie R package to obtain a plurality of known gene regulation relationships from a pathway database; The multiple gene regulation relationships are converted into a network form using an R package igraph to obtain an initial regulation network. One node in the initial regulation network represents one gene, and one edge represents a regulation relationship between two genes. Each edge is regarded as a regulation relationship pair; S62, according to the shortest path standard, searching for a regulation relationship pair between any two genes in each core gene high correlation differential gene module in the initial regulation network, and screening the regulation relationship pairs in each core gene high correlation differential gene module according to the path length to obtain a directed acyclic graph corresponding to all sample specificities in S1 for each core gene high correlation differential gene module; When the shortest path length is 1 or 2, the two genes are a regulation relationship pair; S63, using a structural equation model and the bootstrap samples obtained in S23 to perform causal reasoning on the directed acyclic graph obtained in S62 to generate a gene regulation network corresponding to each core gene high correlation differential gene module: The directed acyclic graph is converted into a structural equation model. The significance and causal effect of the regulation relationship pairs in the gene expression data are calculated according to the gene expression matrix obtained in S231. When the significance value of each regulation relationship pair is less than 0.05, it is the final regulation relationship pair. The regulation relationship pairs are retained, otherwise, each regulation relationship pair is deleted. The retained regulation relationship pairs are used to generate a gene regulation network. The above process is repeatedly executed to obtain a gene regulation network corresponding to each core gene high correlation differential gene module.

2. The method for constructing a gene regulatory network based on meta-analysis according to claim 1, characterized in that: In S1, all normal samples are used as a control group, and all abnormal samples are used as an experimental group.

3. The method for constructing a gene regulatory network based on meta-analysis according to claim 2, characterized in that: The specific process of S21 is as follows: S211、calculating the effect value of a certain gene in the gene expression matrix S211、calculating the effect value of a certain gene in the gene expression matrix wherein, represents the effect size, represents the factor to correct for positive bias, and represents the sample size of the experimental group and the control group, respectively, and represents the gene expression of the experimental group and the control group, respectively, represents the pooled standard deviation of the experimental group and the control group, , and represents the variance of the experimental group and the control group, respectively; S212, calculating the effect value variance of the genes in S211: wherein, denotes the variance of the effect value; The above process is repeatedly executed to obtain the effect value and variance of each gene in each gene expression matrix.

4. The method for constructing a gene regulatory network based on meta-analysis according to claim 3, wherein: The specific process of S22 is as follows: S221, calculating the combined effect value of each gene: wherein, is the combined effect value, is the variance of the combined effect value, is the weight assigned to each gene expression matrix in the effect value , specifically the inverse of the variance of the effect value; S222 the screening parameter for a gene is significance value or false discovery rate when the value is the screening parameter, the screening criterion is or 5 or 1 ; when the is the screening parameter, the screening criterion is or or ; S223、convert each combination effect value to a standard normal distribution statistic, and then perform a two-tailed hypothesis test to obtain the corresponding screening parameter value: S224, based on the screening parameter value, screening the genes according to the screening standard to obtain screened genes and a plurality of gene expression matrices corresponding to the screened genes.

5. The method for constructing a gene regulatory network based on meta-analysis according to claim 4, wherein: The specific process of S23 is as follows: S231, merging the plurality of gene expression matrices obtained after screening in S224 according to columns into one matrix. Each column in the merged matrix is used as a sampling unit for bootstrapping. The number of bootstrapping times and the sample size of each bootstrapping are defined. A plurality of sampling units are obtained by bootstrapping each time. The plurality of sampling units are merged into one gene expression matrix according to columns. The gene expression matrix is used as a bootstrap sample corresponding to the bootstrap round to obtain a plurality of bootstrap samples; The gene space contained in each sampling unit is consistent, that is, each sampling unit has the same gene expression data. Therefore, the same sampling unit repeatedly sampled in the above bootstrapping is used as a new sampling unit. The gene space represents a set of genes in a transcriptome data set or a set of rows of a gene expression matrix; S232, compute weight of each gene in each bootstrap sample : in, It is the first The weight of each gene, This refers to the sample category information of the sampling unit, i.e., the experimental group or the control group. It is the first One gene and Between Score, It is the first in the corresponding transcriptomics gene expression dataset. The first gene and the first Pearson correlation coefficient between genes; In each bootstrap sample, the genes are sorted according to the weight a of each gene from large to small to obtain a sorted gene list. The above process is performed on all bootstrap samples respectively to obtain a plurality of gene lists. The specific process is as follows: in, For the first The first gene in the list Ranking score of each gene For the first The first gene in the list The ranking of individual genes For the first The number of genes in the gene list; in, Used to calculate the In the 1st bootstrap sample The ranking score of each gene and deviation, It is a constant; S233. After b self-sampling, the following hypothesis test is performed for each gene in S233. After b self-sampling, the following hypothesis test is performed for each gene in H0 hypothesis: the th gene has no information for the given condition / trait, i.e. , is the median expected rank score of the th gene in the bootstrap sample under the null hypothesis, i.e. ; H1 hypothesis: the first gene is informative for a given condition / trait, i.e. ; go through Secondary autosampler sampling was used to obtain each gene. indivual , According to Size pair indivual Sort in ascending order to obtain the sorted results for each gene. A list, that is Assign rankings 1, 2, ... And keep in mind Whether the result is positive or negative, the test statistic D+ for each gene is defined as the value of the corresponding gene. At 0 o'clock, The sum of their rankings; S234、define an independent Bernoulli variable as B, B is 1 when and 0 when the mean and variance of the test statistic D+ are: Under the H0 hypothesis In this case, for Then the mean of the test statistic D+ is The variance is When the number of self-sampling attempts is 200, D+ follows a standard normal distribution; S235, performing hypothesis test on the test statistic D+ of each gene with standard normal distribution to obtain the significance of each gene value, and obtaining multiple cross-dataset significant difference genes by using the screening standard, combining all the significant difference genes into a list to obtain a significant difference gene list.

6. The method for constructing a gene regulatory network based on meta-analysis according to claim 5, wherein: In the S232 The specific formula for the score is: wherein, is the mean of the expression level of the th gene in the experimental group samples, is the mean of the expression level of the th gene in all samples, is the mean of the expression level of the th gene in the control group samples, is the mean of the expression level of the th gene in the experimental group samples, is the mean of the expression level of the th gene in the control group samples.

7. The method for constructing a gene regulatory network based on meta-analysis according to claim 6, wherein: The S3 specific process is: S31, obtaining the corresponding transcription factor gene list from the transcription factor database according to the species to which the transcriptomic gene expression data belongs; S32, extracting the same genes in the significant difference gene list obtained in S235 and the transcription factor gene list in S31 to generate a core gene list.

8. The method for constructing a gene regulatory network based on meta-analysis according to claim 7, characterized in that: The specific process of obtaining the co-expression relationship pair between the core gene and the significant difference gene according to the correlation value and the significance value is: (1) when the absolute value of the correlation value of the core gene and the significant difference gene is ≥0.9, and the significance value is <0.01, the core gene and the significant difference gene are a co-expression relationship pair; (2) when the absolute value of the correlation value of the core gene and the significant difference gene is ≥0.8, and the significance value is <0.01, the core gene and the significant difference gene are a co-expression relationship pair; (3) when the absolute value of the correlation value of the core gene and the significant difference gene is ≥0.7, and the significance value is <0.01, the core gene and the significant difference gene are a co-expression relationship pair; (4) when the absolute value of the correlation value of the core gene and the significant difference gene is ≥0.9, and the significance value is <0.05, the core gene and the significant difference gene are a co-expression relationship pair; (5) when the absolute value of the correlation value of the core gene and the significant difference gene is ≥0.8, and the significance value is <0.05, the core gene and the significant difference gene are a co-expression relationship pair; (6) when the absolute value of the correlation value of the core gene and the significant difference gene is ≥0.7, and the significance value is <0.05, the core gene and the significant difference gene are a co-expression relationship pair.

Images