Method for constructing molecular regulation network through multi-omics data and computer device

Abstract

The invention discloses a method for constructing a molecular regulation network through multi-omics data and a computer device. Correlation analysis is carried out on corn multi-omics data, a linear model (LM) and a mixed linear model (MLM) are used for quantifying causal effects of different omics levels in the multi-omics data, and MR analysis is carried out to obtain an MR analysis result; the edge between every two points in the causal network deduced by the causal network module represents the weight of the MR effect between every two points, and finally, a target phenotype molecular regulation and control network can be constructed, so that high-dimensional omics data can be efficiently and accurately analyzed. The computer device and the method can be used for mining regulatory genes of target character phenotypes and breeding or quality improvement of the target character phenotypes.

Description

Technical Field

[0001] This invention belongs to the field of bioinformatics technology, specifically relating to a method and computer device for constructing molecular regulatory networks using multi-omics data. Background Technology

[0002] Genome-wide association studies (GWAS) are a key genetic approach for identifying genomic regions associated with crop-specific traits. Typically, genomic variation initially produces differences at the molecular trait (mTrait) scale, affecting gene expression levels or metabolite composition, which in turn lead to observed differences in phenotypic traits (pTraits). However, traditional GWAS methods often face challenges for complex traits involving interactions of multiple quantitative trait loci (qTLs) with relatively minor influences and environmental factors, as they rely on simple correlation analyses between single-gene variations and phenotypes, thus limiting their ability to fully reveal complex genetic structures and infer causal relationships between genetic variation and phenotypic changes. The rapid development of multi-omics technologies such as transcriptomics and metabolomics in recent years has enabled researchers to integrate mTrait data from various sources, allowing for a more comprehensive exploration of the relationship between genetic variation and pTraits. Through multi-omics association analysis, complex pTraits are broken down into their basic mTrait components, thus enabling the analysis of causal relationships between genetic variation and phenotypic changes at the molecular level. Mendelian randomization (MR) analysis is a key technique in epidemiological causal inference, using genetic variation as an instrumental variable to estimate the causal impact of environmental exposures on phenotypic outcomes. This approach has gained significant appeal in recent years, particularly in elucidating the findings of multi-omics association studies. Therefore, combining multi-omics association analysis with MR methods helps identify key genes influencing specific crop traits, enhances the biological understanding of trait evoked responses, reveals the complexity of genetic regulatory networks, and offers considerable promise for crop genetic improvement research. Although tools such as GEMMA, GAPIT, and rMVP are widely used in GWAS, they primarily cater to traditional GWAS designs and have limited capabilities in multi-omics association analysis. Summary of the Invention

[0003] The technical problem to be solved by this invention is how to construct molecular regulatory networks based on multi-omics data and / or how to mine regulatory genes of target phenotypes based on multi-omics data.

[0004] To address the aforementioned technical problems, the present invention first provides a computer device, including a memory, a processor, and a computer program stored in the memory, wherein the processor can execute the computer program to perform the following steps:

[0005] S1. Data reception: Receiving multi-omics data; the multi-omics data includes genotype data related to the target trait, transcriptome data related to the target trait, metabolome data related to the target trait, and / or phenotypic data related to the target trait;

[0006] S2. Association Analysis: The genotype data and transcriptome data are analyzed using a mixed linear model, and the genotype data and metabolome data are analyzed using a mixed linear model to obtain high-quality genotype data associated with the target phenotype.

[0007] S3. Mendelian randomization analysis: Using the target phenotype-associated genotype data as instrumental variables, the transcriptomic data and / or the metabolomic data as exposure factors, the metabolomic data and / or the target trait-related phenotypic data as outcome variables, and the high-quality target phenotype-associated genotype data as intermediates, Mendelian randomization analysis is performed to estimate the causal relationship between the exposure factors and trait changes, and the MR analysis results are obtained.

[0008] S4. MR Network Analysis: Based on the MR analysis results, a causal network is constructed; the ClusterONE algorithm is applied to identify the functional regulatory modules within the causal network, and a molecular regulatory network related to the target phenotype is constructed.

[0009] Each point in the causal network represents a gene corresponding to a genotype associated with the target phenotype, and the edge between any two points in the causal network represents the weight of the MR effect between any two points; the weight of the MR effect is obtained based on the MR analysis results in S3.

[0010] In the aforementioned computer device, the weight of the MR effect between every two points in S4 can be calculated using the following method:

[0011] Each pair of points is named gene A and gene B. Based on the MR analysis results in S3, the MR effect MR-AB between gene A and gene B with the expression level of gene A as the exposure factor is obtained. Based on the MR analysis results in S3, the MR effect MR-BA between gene A and gene B with the expression level of gene B as the exposure factor is obtained. The average value of MR-AB and MR-BA is calculated using the following formula (6) to obtain the weight of the MR effect between gene A and gene B:

[0012]

[0013] In equation (6), W (A，B)The p-value represents the weight of the relationship between gene A and gene B, where P(A, B) represents the p-value of MR-AB and P(B, A) represents the p-value of MR-BA.

[0014] In the aforementioned computer device, the processor executing the computer program can also perform GO enrichment analysis, which includes the step of performing GO enrichment analysis on genes within the functional regulation module to obtain GO enrichment results of the genes.

[0015] In the aforementioned computer device, the correlation analysis described in S2 may include the following steps:

[0016] S2-1) Genotype dimensionality reduction: Based on the linkage disequilibrium principle, the genotype data is reduced in dimensionality using DBSCAN and / or PCA analysis to obtain pseudo-genotype files;

[0017] S2-2) Phenotypic association interval screening: The pseudogenotype files are filtered using a general linear model and a mixed linear model respectively to obtain candidate association intervals related to the target phenotype;

[0018] S2-3) Local association analysis: SNP information is obtained by SNP retrieval of the candidate association intervals. Local association analysis is performed on the SNP information and the transcriptome data using a mixed linear model. Local association analysis is performed on the SNP information and the metabolome data using a mixed linear model to obtain the genotype data associated with the target phenotype.

[0019] In the aforementioned computer device, the Mendelian randomization analysis described in S3 can employ linear models (LM) and mixed linear models (MLM) to quantify the causal effects at different omics levels in the multi-omics data, thereby obtaining the MR analysis results.

[0020] In the linear model (LM) described above, z represents the genotype of the leader SNP corresponding to the molecular trait (mTrait), x represents the expression of mTrait, and y represents the value of the phenotypic trait (pTrait). Subsequently, z is used in relation to y(b) zy ) and z for x(b zx The least squares estimate of the MR effect of x on y (using b) xy (represented). Then, the MR effect of x on y is expressed as follows: (1):

[0021] b xy =b zy / b zx Equation (1);

[0022] The χ² statistical test is performed using the following formula (2) to evaluate b. xy Statistical significance:

[0023]

[0024] In equation (2), T MR Represents the chi-square statistical test value, var(b) xy The result is obtained by calculation using the following formula (3):

[0025]

[0026] In equation (3), n is the number of samples. This represents the interpretable variance of x with respect to y; Let z represent the interpretable variance of z with respect to x.

[0027] In the above mixed linear model (MLM), b can be... zy and b zx These represent the effects of the leading SNP z on y and x, estimated through association analysis. The kinship of the samples can be directly considered in the MLM model.

[0028] The statistical significance of the MR effect of x on y(bxy) can also be tested using the following equation (4) with a χ² statistical test:

[0029]

[0030] In equation (4), var(b) xy The result is obtained by calculation using the following formula (5):

[0031]

[0032] In equation (5), (b zy ) and (b zx ) represent b respectively zy and b zx The standard error.

[0033] To address the aforementioned technical problems, the present invention also provides an apparatus for constructing a molecular regulatory network related to a target phenotype, the apparatus comprising the following modules:

[0034] A1. Data receiving module: used to receive multi-omics data; the multi-omics data includes genotype data related to the target trait, transcriptome data related to the target trait, metabolome data related to the target trait, and / or phenotypic data related to the target trait;

[0035] A2. Association Analysis Module: Used to perform association analysis on the genotype data and the transcriptome data using a mixed linear model, and to perform association analysis on the genotype data and the metabolome data using a mixed linear model to obtain high-quality genotype data associated with the target phenotype;

[0036] A3. Mendelian Randomization Analysis Module: This module is used to estimate the causal relationship between the exposure factors and trait changes by using the target phenotype-associated genotype data as instrumental variables, the transcriptomic data and / or the metabolomic data as exposure factors, the metabolomic data and / or the target trait-related phenotypic data as outcome variables, and the high-quality target phenotype-associated genotype data as intermediates. The result is obtained through Mendelian randomization analysis.

[0037] A4. MR Network Analysis Module: Used to construct a causal network based on the MR analysis results; apply the ClusterONE algorithm to identify the functional regulatory modules within the causal network, and construct a molecular regulatory network related to the target phenotype.

[0038] Each point in the causal network represents a gene corresponding to a genotype associated with the target phenotype, and the edge between any two points in the causal network represents the weight of the MR effect between any two points; the weight of the MR effect is obtained based on the MR analysis results in A3.

[0039] In the aforementioned apparatus, the weight of the MR effect between every two points in A4 can be calculated using the following method:

[0040] Each pair of points is named gene A and gene B. Based on the MR analysis results in A3, the MR effect MR-AB between gene A and gene B with the expression level of gene A as the exposure factor is obtained. Based on the MR analysis results in A3, the MR effect MR-BA between gene A and gene B with the expression level of gene B as the exposure factor is obtained. The average value of MR-AB and MR-BA is calculated using the following formula (6) to obtain the weight of the MR effect between gene A and gene B:

[0041]

[0042] In equation (6), W (A，B) The p-value represents the weight of the relationship between gene A and gene B, where P(A, B) represents the p-value of MR-AB and P(B, A) represents the p-value of MR-BA.

[0043] In the aforementioned apparatus, the Mendelian randomization analysis module described in A3 can employ linear models (LM) and mixed linear models (MLM) to quantify the causal effects at different omics levels in the multi-omics data, thereby obtaining the MR analysis results.

[0044] To address the aforementioned technical problems, the present invention also provides a method for constructing a molecular regulatory network related to a target phenotype, the method comprising the following steps:

[0045] S1. Data reception: Receiving multi-omics data; the multi-omics data includes genotype data related to the target trait, transcriptome data related to the target trait, metabolome data related to the target trait, and / or phenotypic data related to the target trait;

[0046] S2. Association Analysis: The genotype data and transcriptome data are analyzed using a mixed linear model, and the genotype data and metabolome data are analyzed using a mixed linear model to obtain high-quality genotype data associated with the target phenotype.

[0047] S3. Mendelian randomization analysis: Using the target phenotype-associated genotype data as instrumental variables, the transcriptomic data and / or the metabolomic data as exposure factors, the metabolomic data and / or the target trait-related phenotypic data as outcome variables, and the high-quality target phenotype-associated genotype data as intermediates, Mendelian randomization analysis is performed to estimate the causal relationship between the exposure factors and trait changes, and the MR analysis results are obtained.

[0048] S4. MR Network Analysis: Based on the MR analysis results, a causal network is constructed; the ClusterONE algorithm is applied to identify the functional regulatory modules within the causal network, and a molecular regulatory network related to the target phenotype is constructed.

[0049] Each point in the causal network represents a gene corresponding to a genotype associated with the target phenotype, and the edge between any two points in the causal network represents the weight of the MR effect between any two points; the weight of the MR effect is obtained based on the MR analysis results in S3.

[0050] In the above method, the weight of the MR effect between every two points in S4 can be calculated using the following method:

[0051] Each pair of points is named gene A and gene B. Based on the MR analysis results in S3, the MR effect MR-AB between gene A and gene B with the expression level of gene A as the exposure factor is obtained. Based on the MR analysis results in S3, the MR effect MR-BA between gene A and gene B with the expression level of gene B as the exposure factor is obtained. The average value of MR-AB and MR-BA is calculated using the following formula (6) to obtain the weight of the MR effect between gene A and gene B:

[0052]

[0053] In equation (6), W (A，B)The p-value represents the weight of the relationship between gene A and gene B, where P(A, B) represents the p-value of MR-AB and P(B, A) represents the p-value of MR-BA.

[0054] The above method may also include a GO enrichment analysis step, wherein the GO enrichment analysis includes the step of performing GO enrichment analysis on genes within the functional regulatory module to obtain GO enrichment results of the genes.

[0055] In the above method, the association analysis described in S2 may include the following steps:

[0056] S2-1) Genotype dimensionality reduction: Based on the linkage disequilibrium principle, the genotype data is reduced in dimensionality using DBSCAN and / or PCA analysis to obtain pseudo-genotype files;

[0057] S2-2) Phenotypic association interval screening: The pseudogenotype files are filtered using a general linear model and a mixed linear model respectively to obtain candidate association intervals related to the target phenotype;

[0058] S2-3) Local association analysis: SNP information is obtained by SNP retrieval of the candidate association intervals. Local association analysis is performed on the SNP information and the transcriptome data using a mixed linear model. Local association analysis is performed on the SNP information and the metabolome data using a mixed linear model to obtain the genotype data associated with the target phenotype.

[0059] In the above method, the Mendelian randomization analysis described in S3 can use linear models (LM) and mixed linear models (MLM) to quantify the causal effects at different omics levels in the multi-omics data, and obtain the MR analysis results.

[0060] In the linear model (LM) described above, z represents the genotype of the leader SNP corresponding to the molecular trait (mTrait), x represents the expression of mTrait, and y represents the value of the phenotypic trait (pTrait). Subsequently, z is used in relation to y(b) zy ) and z for x(b zx The least squares estimate of the MR effect of x on y (using b) xy (represented). Then, the MR effect of x on y is expressed as follows: (1):

[0061] b xy =b zy / b zx Equation (1);

[0062] The χ² statistical test is performed using the following formula (2) to evaluate b. xy Statistical significance:

[0063]

[0064] In equation (2), T MR Represents the chi-square statistical test value, var(b) xy The result is obtained by calculation using the following formula (3):

[0065]

[0066] In equation (3), n is the number of samples. This represents the interpretable variance of x with respect to y; Let z represent the interpretable variance of z with respect to x.

[0067] In the above mixed linear model (MLM), b can be... zy and b zx These represent the effects of the leading SNP z on y and x, estimated through association analysis. The kinship of the samples can be directly considered in the MLM model.

[0068] The statistical significance of the MR effect of x on y(bxy) can also be tested using the following equation (4) with a χ² statistical test:

[0069]

[0070] In equation (4), var(b) xy The result is obtained by calculation using the following formula (5):

[0071]

[0072] In equation (5), (b zy ) and (b zx ) represent b respectively zy and b zx The standard error.

[0073] The following applications of the computer device described above, the device described above, and / or the computer-readable storage medium described above are also within the scope of protection of this invention:

[0074] P1. Application in the development or preparation of products that uncover regulatory genes of the target phenotype;

[0075] P2. Application in the development or preparation of products that tap into breeding hotspots for the target phenotype;

[0076] P3. Application in breeding for the desired phenotype or quality improvement;

[0077] P4. Applications in the development or preparation of disease-related drugs.

[0078] The GO mentioned above can be interpreted as Gene Ontology. DBSCAN, also mentioned above, stands for Density-Based Spatial Clustering of Applications with Noise.

[0079] The computer program product described in this invention can be a software product that primarily implements its solution through a computer program.

[0080] The computer-readable storage medium refers to a carrier for storing data, which may be magnetic tape, disk, floppy disk, optical disk, magneto-optical disk, ROM, PROM, VCD, DVD, hard disk, flash memory, USB flash drive, CF card, SD card, MMC card, SM card, Memory Stick, or xD card, etc.

[0081] Due to the characteristics of multi-omics data, such as wide sources, high noise, and high redundancy, traditional single statistical models in association analysis are difficult to apply to the efficient and accurate analysis of high-dimensional omics data. This invention introduces dimensionality reduction and regional association analysis steps into genotype data; uses linear models (LM) and mixed linear models (MLM) to quantify the causal effects at different omics levels in the multi-omics data, and performs MR analysis to obtain MR analysis results; in the causal network module inference, the edge between each two points in the causal network represents the weight of the MR effect between each two points, ultimately achieving efficient and accurate analysis of high-dimensional omics data.

[0082] In this invention, the computer program is named MRBIGR (MR-based Inference of Genetic Regulation). In the embodiments, MRBIGR was used to analyze maize multi-omics data, constructing a molecular regulatory network for flavonoid metabolism. It was found that the gene p1, which regulates seed coat color, occupies an important position in the flavonoid metabolism network, and different genotypes of the p1 GWAS signaling site showed significant differences in the expression of flavonoid metabolism genes, suggesting that the regulatory network in which p1 is located may be a selection hotspot during the domestication process. Attached Figure Description

[0083] Figure 1 GWAS signaling results for p1 and three known structural flavonoid metabolites. A. GWAS analysis results of p1 expression level. B. GWAS analysis results of cosmos glycoside abundance; C. GWAS analysis results of senna flavin abundance; D. GWAS analysis results of hesperidin abundance. From left to right: Manhattan plot, QQ plot, and a local Manhattan plot within the 10Mb-100Mb region of chromosome 1.

[0084] Figure 2 The MR effects of p1 with known metabolites and target genes are shown. A shows the MR significance distribution of p1 with 102 known structural metabolites. The text highlighted in green boxes represents metabolites that passed the significance threshold (P = 4.9e-4). B shows the MR effect distribution of p1 with 102 known structural metabolites. C shows the MR significance distribution of p1 with 16 target genes.

[0085] Figure 3 Enrichment analysis of p1 target genes and construction of a causal network for flavonoid metabolism. A shows the GO enrichment results of p1 target genes. B shows the ClusterONE interaction network of p1, p1 target genes, and annotated flavonoid metabolism genes.

[0086] Figure 4 This invention presents a flowchart of constructing a molecular regulatory network (MRBIGR) based on multi-omics data. Detailed Implementation

[0087] The present invention will now be described in further detail with reference to specific embodiments. The given embodiments are merely illustrative of the invention and not intended to limit its scope. The embodiments provided below can serve as a guide for further improvements by those skilled in the art and do not constitute a limitation on the invention in any way.

[0088] Unless otherwise specified, the experimental methods used in the following examples are conventional methods, performed according to the techniques or conditions described in the literature in this field or according to the product instructions. Unless otherwise specified, the materials and reagents used in the following examples are commercially available.

[0089] Example 1: Flowchart of the method for constructing molecular regulatory networks based on multi-omics data (MRBIGR)

[0090] This invention obtains causal relationships between multi-omics data by performing association analysis, Mendelian randomization analysis, and MR network analysis on multi-omics data, and then constructing network modules and visualizing the results, thereby obtaining a molecular regulatory network. This method is called MRBIGR (MR-based Inference of Genetic Regulation). Figure 4 ).

[0091] 1. Multi-omics data download

[0092] The maize multi-omics dataset encompasses 527 associated populations representing global maize inbred line diversity, and omics data from these maize varieties were measured across different dimensions (Related literature: Yang X, Gao S, Xu S, et al. Characterization of a global germplasm collection and its potential utilization for analysis of complex quantitative traits in maize[J]. Molecular Breeding, 2011, 28(4): 511-526. DOI: 10.1007 / s11032-010-9500-7; Fu J, Cheng Y, Linghu J, et al. RNA sequencing reveals the complex regulatory network in the maize kernel. Nat Commun. 2013; 4:2832.), specifically including: 1) SNP genotype data of 527 maize inbred lines, and 2) transcriptome and metabolome data of 368 maize inbred lines 15 days after pollination (data source: http: / / www.maizego.org / Resources.html).

[0093] 2. Multi-omics data preprocessing

[0094] Material filtering was performed on the genotype data obtained in Step 1 from the multi-omics data, removing materials with a missing rate greater than 0.2. Then, SNP filtering was performed, removing sites with a missing rate greater than 0.1 and a minor allele frequency (MAF) less than 0.05. Genotype inference and phase reconstruction were performed on the missing genotype sites using Beagle software (version v3.3.2) (Browning and Browning 2007). If the posterior probability of the inferred genotype was less than 0.5, the SNP site containing that missing genotype was deleted to obtain highly accurate genotypes. The generated file was converted to HapMap format as a high-density genomic variation map (hapmapSNPs) for subsequent analysis. Considering that millions of hapmapSNPs may still contain significant redundancy, these SNP sites may be located within the same linkage disequilibrium (LD) region, reflecting the same genetic diversity. Therefore, this invention uses hapmapSNPs data as material and constructs tagged SNPs using the snp_clumping function of bigsnpr (Privé et al. 2018). The correlation between SNPs is set to 0.7, and the window size is 125Kb.

[0095] This invention uses RPKM (reads per kilobase per m The IPKM value was used to quantify the expression level of each gene in the transcriptome data from the multi-omics dataset obtained in step 1. To obtain a high-quality gene expression matrix, genes with low expression were first filtered out by a criterion that the RPKM value was ≥0.3 in at least a quarter of the samples. Then, the expression levels were normalized using normal quantiles, resulting in a high-quality expression matrix containing 25,558 genes.

[0096] Metabolomic data from the multi-omics dataset obtained in step 1 were analyzed to identify 735 metabolites, of which only 102 had known structures, originating from 15 metabolic pathways. Therefore, the analysis below focuses only on these 102 metabolites. Furthermore, to reduce the dispersion of metabolite abundance, the abundance of each metabolite was log2 transformed to obtain the transformed metabolomics data.

[0097] 3. Correlation Analysis

[0098] 3.1 Multi-omics association analysis

[0099] Association analysis was performed using MODAS (Multi-Omics Data Association Study, related literature: Liu S, Xu F, Xu Y, Wang Q, Yan J, Wang J, Wang X, Wang X. MODAS: exploring maizegermplasm with multi-omics data association studies. Sci Bull (Beijing). 2022 May 15; 67(9): 903-906. Doi: 10.1016 / j.scib.2022.01.021. Epub 2022 Jan 31. PMID: 36546022). Phenotypic data can be gene expression levels in transcriptome data (high-quality expression matrix data obtained in step 2), metabolite content in metabolome data (transformed metabolome data obtained in step 2), or phenotypic values. The main steps are as follows:

[0100] 3.1.1 Perform genotypic dimensionality reduction.

[0101] Based on the linkage disequilibrium principle, dimensionality reduction of the genotype data (high-quality genotype data obtained in step 2) was performed using DBSCAN (Density-Based Spatial Clustering of Applications with Noise, related literature: Martin Ester, Hans-Peter Kriegel, Joerg Sander, Xiaowei Xu (1996). A Density-Based Algorithm for Discovering Clusters in Large Spatial Databases with Noise. Institute for Computer Science, University of Munich. Proceedings of 2nd International Conference on Knowledge Discovery and Data Mining (KDD-96)) and PCA principal component analysis to obtain pseudo-genotype files. This step preserves genetic diversity while reducing the number of SNPs required for analysis, greatly improving the efficiency of GWAS analysis.

[0102] 3.1.2 Phenotypic association interval filtering.

[0103] Through two-step association analysis, the pseudo-genotype files obtained by dimensionality reduction were filtered using a general linear model and a mixed linear model (MLM) to obtain candidate association intervals related to the phenotype (containing phenotype-related genotype data, i.e., leader SNPs).

[0104] 3.1.3 Local correlation analysis.

[0105] Using the association interval information obtained in step 3.1.2, SNP retrieval is performed to obtain leader SNP (phenotypic association genotype data) information. Then, the MLM method is used to perform local association analysis on the leader SNP information and phenotypic data (transcriptomics data, metabolomics data, or phenotypic data). Finally, the results are merged and filtered to obtain high-quality phenotypic association intervals (high-quality phenotypic association genotype data).

[0106] 4. Mendelian randomization analysis

[0107] MRBIGR's Mendelian randomization analysis module employs linear (LM) and mixed linear (MLM) models to quantify causal effects at different omics levels. In this step, molecular traits (mTraits), such as gene expression, metabolite content, or protein content, are treated as exposure variables, while phenotypic traits (pTraits), such as metabolite content, protein content, or agronomic traits, are treated as outcome variables.

[0108] When performing Mendelian randomization (MR) analysis, an intermediate form is needed to estimate the causal relationship between the exposure variable and the trait change. This invention uses high-quality association intervals obtained from local association analysis as an intermediate form to obtain MR analysis results.

[0109] 4.1 Linear Model (LM) Analysis

[0110] In the linear model (LM), z represents the genotype of the leader SNP corresponding to the molecular trait (mTrait), x represents the expression of mTrait, and y represents the value of the phenotypic trait (pTrait). Subsequently, z is compared with y(b)... zy ) and z for x(b zx The least squares estimate of the MR effect of x on y (using b) xy (represented). Then, the MR effect of x on y is expressed as follows: (1):

[0111] b xy =b zy / b zx Equation (1);

[0112] In order to evaluate b xy The statistical significance of the χ² test was determined using the following formula (2):

[0113]

[0114] In equation (2), T MR Represents the chi-square statistical test value, var(b) xy The result is obtained by calculation using the following formula (3):

[0115]

[0116] In equation (3), n is the number of samples. This represents the interpretable variance of x with respect to y; Let z represent the interpretable variance of z with respect to x.

[0117] 4.2 Mixed Linear Model (MLM) Analysis

[0118] In a mixed linear model (MLM), b zy and b zx These represent the effects of the leading SNP z on y and x, estimated through association analysis. In this way, the kinship of samples is directly considered in the MLM model.

[0119] The statistical significance of the MR effect of x on y(bxy) can also be tested using the following equation (4) with a χ² statistical test:

[0120]

[0121] In equation (4), var(b) xy The result is obtained by calculation using the following formula (5):

[0122]

[0123] In equation (5), (b zy ) and (b zx ) represent b respectively zy and b zx The standard error.

[0124] The exposure variables and result files used here are flexible and can vary, including gene expression levels and metabolite abundance, or only gene expression levels for gene-pair analysis. Example 2 below specifically demonstrates causal relationship analysis between gene expression levels and metabolite abundance, causal regulatory relationship analysis between gene expression, and causal relationship analysis between metabolite abundance and agronomic traits.

[0125] 5. MR Network Analysis (net):

[0126] 5.1 Constructing a causal network

[0127] MRBIGR introduces a novel method for constructing MR-based networks with weights for the MR effects between gene pairs (i.e., each node in the network represents a gene, and the edge between two nodes represents the weight of the MR effect between gene pairs). The genes are those corresponding to the leading SNP genotypes. The effects of MR are inherently asymmetric; that is, the MR effect between gene A and gene B when gene A is expressed at exposure level differs from the MR effect when gene B is expressed at exposure level. Therefore, the weights assigned to the relationship between gene A and gene B are the average of these two MR effects. The calculation formula is as follows (6):

[0128]

[0129] In equation (6), W (A，B) The weights represent the relationship between gene A and gene B; when the expression level of gene A is exposure, P(A, B) represents the p-value of the MR between gene A and gene B, and when the expression level of gene B is exposure, P(B, A) represents the p-value of the MR between gene A and gene B. MRBIGR determines the significance threshold (cutoff in equation (6)) of the MR effect between gene pairs based on the data, and assigns the weight of the effect below the specified p-value threshold to 0. 5.2 Functional Module Identification

[0130] The ClusterONE (clustering with overlapping neighborhood expansion) algorithm (related literature: T. Nepusz, H. Yu, and A. Paccanaro. Detecting overlapping protein complexes in protein-protein interaction networks. Nature Methods, vol. 9, pp. 471-472, 2012) was applied to detect modules based on the MR-derived network (the network obtained in step 5.1) in MRBIGR. ClusterONE employs a greedy strategy to identify highly cohesive groups. Initially, it starts with randomly selected seed points and progressively adds points to enhance module cohesion until the cumulative edge weights within the module reach a sufficient level compared to the sum of edge weights connecting genes in the module. To identify modules in MRBIGR, the transformed MR effect is used as a weight metric between points. A module is considered biologically relevant when its calculated P-value is below a specified cutoff and the number of points within the module is at least 5. Subsequently, the hub_score function in the R-pack graph is applied to calculate the center score of points within each module. Points with a center score greater than or equal to 0.8 are defined as centers. Finally, the ggnet package is used for visualization of network modules.

[0131] 6. GO enrichment analysis

[0132] ClusterProfiler (version 3.14.0) was used to perform GO (Gene Ontology) functional enrichment analysis on the gene set to obtain the GO enrichment results of genes in the functional regulatory module (genes corresponding to the genotype data in the functional module). The significance threshold was set to pvalue ≤ 0.05. Since GO enrichment analysis requires specifying a background annotation set, and currently there is no complete maize annotation information in AnnotationHub, this invention provides the background annotation set in two ways: 1) Annotating the target sequence using InterProScan and extracting the GO number corresponding to each isoform; 2) Downloading maize GO annotation information from the Ensembl plant biomart database, version Ensembl plantGene 51 (https: / / plants.ensembl.org / biomart / martview).

[0133] Example 2: Constructing and analyzing a molecular regulatory network for flavonoid metabolism using the method of the present invention.

[0134] 1. p1 co-localization with GWAS signaling of flavonoid metabolites

[0135] The results of Liu et al.’s research indicate that there is an eQTL hotspot regulating the flavonoid metabolism pathway in a region of about 48 Mb on maize chromosome 1. This hotspot is located upstream of several Myb transcription factors, such as the genes p1 and p2 (pericarp color) that regulate seed coat color, and three Myb transcription factors including mybr102 (related literature: Liu H, Luo X, Niu L, Xiao Y, Chen L, Liu J, Wang X, Jin M, Li W, Zhang Q, Yan J (2017a) Distant eQTLs and Non-coding Sequences Play Critical Roles in Regulating Gene Expression and Quantitative Trait Variation in Maize. Mol Plant 10:414-426). This invention uses an eQTL hotspot regulating the flavonoid metabolism pathway located approximately 48 Mb on maize chromosome 1 as the genotype, and performs association analysis using the expression levels of p1, p2, and mybr102, as well as the abundance of 33 flavonoid metabolites (Table 1), as the phenotypes. The associations are visualized using MRBIGR from Example 1 to analyze whether these genes are directly related to the expression of flavonoid metabolites. The results show that only p1 among the three Mybr1 transcription factors exhibits a significant GWAS signal (…). Figure 1 (A). Furthermore, this site not only co-localizes with the GWAS signaling of nearly half of the flavonoid metabolites (14 / 33 = 42.4%), but also... Figure 1 The gene p1 (BD) also perfectly overlaps with the eQTL hotspot mentioned by Liu et al. (related literature: Liu H, Luo X, Niu L, Xiao Y, Chen L, Liu J, Wang X, Jin M, Li W, Zhang Q, Yan J (2017a) Distant eQTLs and Non-coding Sequences Play CriticalRoles in Regulating Gene Expression and Quantitative Trait Variation in Maize. Mol Plant 10:414-426), indicating that p1 may be a key gene regulating the flavonoid metabolism pathway.

[0136] Table 1. 33 Flavonoid Metabolites

[0137]

[0138]

[0139] 2. Causal role analysis of p1 in the flavonoid metabolic pathway

[0140] Through the cluster analysis and GWAS analysis described above, this invention identified p1 as a potentially key gene regulating the flavonoid metabolism pathway. Furthermore, this invention employed the Mendelian randomization (MR) analysis module of MRBIGR to perform MR causal analysis on p1 and 102 known metabolites to analyze the causal relationship between p1 and flavonoid metabolism.

[0141] Specifically, the peak SNP (Chr1.s_48424403) from the p1 GWAS analysis was used as the genetic variation affecting p1 expression (instrumental variable), p1 expression level as the exposure factor, and the content of 102 known structural metabolites as trait changes (outcome variables). Simultaneously, 4.9e-4 (0.05 / 102) was used as the significance threshold for MR analysis to reduce false positives in the results. Next, MRBIGR was used to visualize the MR analysis results. The results showed that 17 metabolites passed the significance test, and all of them were flavonoid metabolites. Figure 2 (A); In addition, a large number of flavonoid metabolites also showed a high MR effect distribution with p1 ( Figure 2 (B) Since a higher MR effect indicates a stronger regulatory relationship, the above results show that p1 plays a significant causal regulatory role in flavonoid metabolism.

[0142] Given the effectiveness of the MR method in inferring the causal relationship between p1 and flavonoid metabolites, this invention further infers the regulatory relationship between p1 and other genes. Specifically: First, GWAS analysis of 25,558 genes in the maize transcriptome data was performed using MRBIGR, identifying 9,520 eQTLs and corresponding 8,929 genes (including p1); then, the peak SNPs in the GWAS analysis corresponding to p1 and the 8,928 genes were used as genetic variations, and forward and reverse MR analyses were performed on these 8,929 genes based on expression level information. Only genes whose bidirectional P-values ​​with p1 both meet the significance threshold of 1.12e-4 (1 / 8929) were retained; finally, the results were visualized. After screening, a total of 16 genes significantly associated with p1 were obtained. Figure 2 (C). GO functional enrichment analysis of these genes using MRBIGR revealed that they were significantly enriched in pathways related to flavonoid and anthocyanin metabolism. Figure 3 (A)

[0143] 3. Construction of a causal molecular network of flavonoid metabolism genes

[0144] To further construct a causal network of flavonoid metabolism genes, this invention used MRBIGR to perform pairwise MR analysis on p1, 16 genes significantly associated with p1, and genes in the flavonoid metabolism pathway. Then, the MR effects between genes were converted into edge weights to construct a causal network. Finally, a significance threshold of 0.05 was used to filter the results, resulting in a causal network of flavonoid metabolism genes consisting of 21 genes. Figure 3 (B). The hub genes (hub_score≥0.8) in this network include p1, C2 (Zm00001d052673), ZmCGT1 (Zm00001d037382), ZmUGT1 (Zm00001d037383), and FNS1 (Zm00001d047452) mentioned above, further illustrating that p1 plays an important role in the flavonoid metabolism network.

[0145] The present invention has been described in detail above. Those skilled in the art will recognize that the invention can be practiced in a wide range of ways with equivalent parameters, concentrations, and conditions without departing from its spirit and scope, and without requiring unnecessary experiments. While specific embodiments have been provided, it should be understood that further modifications can be made to the invention. In summary, according to the principles of the invention, this application is intended to include any changes, uses, or improvements to the invention, including changes made using conventional techniques known in the art that depart from the scope disclosed herein.

Claims

1. A computer apparatus comprising a memory, a processor, and a computer program stored on the memory, wherein the computer program, when executed by the processor, causes the processor to perform steps comprising: The processor executes the computer program to implement the following steps: S1, data receiving: receiving multi-omics data; the multi-omics data includes genotype data related to a target trait, transcriptome data related to the target trait, metabolome data related to the target trait, and / or phenotype data related to the target trait; S2, association analysis: performing association analysis on the genotype data and the transcriptome data using a mixed linear model, and performing association analysis on the genotype data and the metabolome data using a mixed linear model to obtain high-quality target phenotype-associated genotype data; S3, Mendelian randomization analysis: using the target phenotype-associated genotype data as an instrumental variable, using the transcriptome data and / or the metabolome data as an exposure factor, using the metabolome data and / or the phenotype data related to the target trait as an outcome variable, and using the high-quality target phenotype-associated genotype data as an intermediate type to estimate the causal relationship between the exposure factor and trait change, performing Mendelian randomization analysis to obtain an MR analysis result; S4, MR network analysis: constructing a causal network based on the MR analysis result; applying a ClusterONE algorithm to identify functional regulation modules in the causal network, and constructing a molecular regulation network related to the target phenotype; Each point in the causal network represents a gene corresponding to the target phenotype-associated genotype, and the edge between each two points in the causal network represents the weight of the MR effect between the two points; the weight of the MR effect is obtained based on the MR analysis result in S3.

2. The computer device of claim 1, wherein: The weight of the MR effect between each two points in S4 is obtained by the following method: Name each two points as gene A and gene B, obtain the MR effect MR-AB between gene A and gene B when the expression level of gene A is used as an exposure factor based on the MR analysis result in S3, obtain the MR effect MR-BA between gene A and gene B when the expression level of gene B is used as an exposure factor based on the MR analysis result in S3, calculate the average value of MR-AB and MR-BA by the following formula (6) to obtain the weight of the MR effect between gene A and gene B: In formula (6), W (A，B) The weight representing the relationship between gene A and gene B, P(A, B) represents the p-value of the MR-AB, and P(B, A) represents the p-value of the MR-BA.

3. The computer apparatus according to claim 1 or 2, characterized in that: In the Mendelian randomization analysis in S3, linear models and mixed linear models are used to quantify the causal effects at different levels of multi-omics data to obtain the MR analysis result.

4. A device for constructing a molecular regulatory network related to a target phenotype, characterized in that: The device includes the following modules: A1, data receiving module: for receiving multi-omics data; the multi-omics data includes genotype data related to a target trait, transcriptome data related to the target trait, metabolome data related to the target trait, and / or phenotype data related to the target trait; A2, association analysis module: for performing association analysis on the genotype data and the transcriptome data using a mixed linear model, and performing association analysis on the genotype data and the metabolome data using a mixed linear model to obtain high-quality target phenotype-associated genotype data; A3, Mendelian randomization analysis module: used to estimate the causal relationship between the exposure factor and trait change by taking the purpose phenotype associated genotype data as the instrumental variable, taking the transcriptome data and / or the metabolome data as the exposure factor, taking the metabolome data and / or the purpose trait related phenotype data as the outcome variable, and taking the high quality purpose phenotype associated genotype data as the intermediate type, and performing Mendelian randomization analysis to obtain MR analysis results; A4, MR network analysis module: used to construct a causal network based on the MR analysis results; and identify functional regulation modules in the causal network by applying the ClusterONE algorithm to construct a purpose phenotype related molecular regulation network; Each point in the causal network represents a gene corresponding to the purpose phenotype associated genotype, and the edge between each two points in the causal network represents the weight of the MR effect between the two points; the weight of the MR effect is obtained based on the MR analysis results in A3.

5. The apparatus of claim 4, wherein: The weight of the MR effect between each two points in A4 is obtained by the following method: Name each two points as gene A and gene B, obtain the MR effect MR-AB between the gene A and the gene B when the expression level of the gene A is the exposure factor based on the MR analysis results in A3, obtain the MR effect MR-BA between the gene A and the gene B when the expression level of the gene B is the exposure factor based on the MR analysis results in A3, and calculate the average value of the MR-AB and the MR-BA by formula (6) to obtain the weight of the MR effect between the gene A and the gene B: In formula (6), W (A，B) The weight representing the relationship between gene A and gene B, P(A, B) represents the p-value of the MR-AB, and P(B, A) represents the p-value of the MR-BA.

6. The apparatus of claim 4 or 5, wherein: In the Mendelian randomization analysis module A3, linear models and mixed linear models are used to quantify the causal effects at different levels of the multi-omics data to obtain the MR analysis results.

7. A method of constructing a molecular regulatory network associated with a phenotype of interest, characterized by: The method comprises the following steps: S1, data receiving: receiving multi-omics data; the multi-omics data comprises purpose trait related genotype data, purpose trait related transcriptome data, purpose trait related metabolome data, and / or purpose trait related phenotype data; S2, association analysis: performing association analysis on the genotype data and the transcriptome data using a mixed linear model, and performing association analysis on the genotype data and the metabolome data using a mixed linear model to obtain high quality purpose phenotype associated genotype data; S3, Mendelian randomization analysis: used to estimate the causal relationship between the exposure factor and trait change by taking the purpose phenotype associated genotype data as the instrumental variable, taking the transcriptome data and / or the metabolome data as the exposure factor, taking the metabolome data and / or the purpose trait related phenotype data as the outcome variable, and taking the high quality purpose phenotype associated genotype data as the intermediate type, and performing Mendelian randomization analysis to obtain MR analysis results; S4, MR network analysis: constructing a causal network based on the MR analysis results; applying ClusterONE algorithm to identify functional regulatory modules in the causal network, and constructing a molecular regulatory network related to the target phenotype; Each point in the causal network represents a gene corresponding to the genotype associated with the target phenotype, and the edge between each two points in the causal network represents the weight of the MR effect between the two points; the weight of the MR effect is obtained based on the MR analysis results in S3.

8. The method of claim 7, wherein: The weight of the MR effect between each two points in S4 is obtained by the following method: Name each two points as gene A and gene B, obtain the MR effect MR-AB between gene A and gene B when the expression level of gene A is the exposure factor based on the MR analysis results in S3, obtain the MR effect MR-BA between gene A and gene B when the expression level of gene B is the exposure factor based on the MR analysis results in S3, calculate the average value of MR-AB and MR-BA by the following formula (6) to obtain the weight of the MR effect between gene A and gene B: In formula (6), W (A，B) The weight representing the relationship between gene A and gene B, P(A, B) represents the p-value of the MR-AB, and P(B, A) represents the p-value of the MR-BA.

9. The method according to claim 7 or 8, characterized in that: In the Mendelian randomization analysis in S3, linear model and mixed linear model are used to quantify the causal effect at different omics levels in the multi-omics data to obtain the MR analysis results.

10. The following any one application of the computer device of any one of claims 1-3, the device of any one of claims 4-6, and / or the computer readable storage medium of any one of claims 7-9: P1, in the development or preparation of products for mining regulatory genes of target phenotypes; P2, in the development or preparation of products for mining breeding hotspots of target phenotypes; P3, in the breeding or quality improvement of target phenotypes; P4, in the development or preparation of disease-related drugs.

Images