A gene-microorganism-phenotype ternary relationship mining method

Abstract

The application discloses a gene-microorganism-phenotype ternary relationship mining method based on multi-omics data integration and gradient boosting model, and relates to the technical field of bioinformatics and computational biology. The method integrates four types of omics information of genome-wide association study (GWAS) results, gene expression, microorganism abundance data and phenotype data, constructs a "gene-microorganism-phenotype" ternary group in combination with statistical genetics constraints and machine learning joint modeling, and systematically finds biological pathways of the host genome and the microorganism group cooperatively regulating phenotypes. The application solves the problems of the existing method, such as the fragmentation of genome and microorganism association analysis, the difficulty in capturing nonlinear interaction effects, and low computational efficiency, and provides efficient cross-omics relationship mining for livestock and poultry genetic breeding.

Description

Technical Field

[0001] This invention relates to the fields of bioinformatics and computational biology, specifically to a method for linking host genetics and microbial interactions, which can be applied to fields such as livestock and poultry genetic breeding and research on microbiome-host interaction mechanisms. Background Technology

[0002] With the development of high-throughput sequencing technology, genome-wide association studies (GWAS) have become a core tool for elucidating the genetic basis of complex traits. Meanwhile, advancements in metagenomic sequencing technology have made large-scale microbiome studies possible. A growing body of research indicates a complex interaction between the host genome and the symbiotic microbiome, jointly influencing quantitative traits in the host.

[0003] In the livestock and poultry industry, the interaction between the rumen microbiome and the host genome has been proven to be closely related to important economic traits such as feed conversion efficiency, carcass traits, and meat production traits. However, existing research methods have certain limitations:

[0004] (1) Fragmented association analysis: Traditional methods perform genome-phenotype association analysis (GWAS) and microbiome-phenotype association analysis (MeWAS) separately, making it difficult to examine the synergistic effects of genome and microbiome on phenotype within the same framework;

[0005] (2) Gene-level signal convergence bias: When GWAS signals at the SNP level are converged to the gene level, the difference in the number of SNPs contained in the gene region will introduce a systematic bias—genes containing more SNPs are more likely to show significant signals, thus leading to the overestimation of large genes.

[0006] (3) Ignoring nonlinear effects: Existing association analyses are mostly based on linear models, which make it difficult to capture the nonlinear interaction effects that may exist between gene expression and microbial abundance;

[0007] (4) Computational efficiency bottleneck: As the number of genes, microorganisms and phenotypes increases, the combination space of triples grows cubically, and the computational cost of exhaustive evaluation is extremely high.

[0008] (5) Lack of a systematic ternary relationship quantification framework: Although existing studies have reported the association between gene-microbe and microbe-phenotype, there is still a lack of a methodological framework that can integrate the two into a unified "gene-microbe-phenotype" pathway and quantify their contribution.

[0009] Therefore, there is an urgent need for a method that can efficiently integrate genomic, microbiome, and phenotypic data to systematically mine the "gene-microbe-phenotype" ternary relationship. Summary of the Invention

[0010] In order to overcome the shortcomings of the prior art, at least one technical problem raised in the background art is solved.

[0011] The technical solution adopted by this invention to solve its technical problem is: a gene-microorganism-phenotype ternary relationship mining method of this invention, comprising the following steps:

[0012] Step S1: Acquire and preprocess four types of omics data, including phenotypic matrix P, microbiome abundance matrix M, gene expression matrix G, and gene-microbe association matrix Z;

[0013] Step S2: Construct a gene-level GWAS association signal, pool the association P-values ​​at the SNP level to the gene level, obtain the gene-level association P-value for each gene-microbe pair, convert it to Z-score, and then perform linear regression correction with the number of SNPs in the gene region as the independent variable to obtain the corrected Z-score matrix that eliminates the bias of the number of SNPs.

[0014] Step S3: Apply threshold constraints to the corrected Z-score matrix to screen out gene-microbe candidate pairs with statistically significant genetic associations;

[0015] Step S4: Calculate the multi-level correlation matrix, including the gene-phenotype correlation coefficient matrix R, through vectorized matrix operations. GP Microbial-phenotypic correlation coefficient matrix R MP And the gene-microbe expression correlation coefficient matrix R GM ;

[0016] Step S5: For each gene-microbe candidate pair obtained through step S3, form a candidate ternary relationship with each phenotype combination, calculate a preliminary score according to the comprehensive scoring function, and retain the Top-N1 candidates according to the score; the comprehensive scoring function includes gene-microbe GWAS Z-score, R² aggregation term, square term of gene-microbe expression correlation coefficient, and directional consistency reward factor.

[0017] Step S6: For the candidate ternary relationships selected in Step S5, a gradient boosting decision tree model is used for joint modeling and evaluation. A gene single-factor model, a microbial single-factor model, and a joint model containing genes, microorganisms, and their interaction terms are constructed respectively. The joint gain and the importance of interaction features are calculated.

[0018] Step S7: Calculate the final comprehensive score based on the joint model performance, mediation gain, and interaction effect strength, and output the Top-N2 ternary relations in order of score.

[0019] Preferably, in step S2, the association P-values ​​at the SNP level are pooled to the gene level using the mBAT-combo method, and multivariate association analysis is performed using the association statistics and linkage disequilibrium structure of all SNPs within the gene region.

[0020] Preferably, the specific steps of linear regression correction in step S2 are as follows: converting the gene-level association P-value into the original Z-score using a standard normal quantile function, and then using the logarithmic number of SNPs as the basis for the Z-score. Using the original Z-score as the dependent variable and performing linear regression with the original Z-score as the independent variable, the regression residuals are taken and standardized with Z-score to obtain the corrected gene-level Z-score.

[0021] Preferably, the calculation of the multi-level correlation matrix in step S4 involves: performing column standardization on the input matrices X (dimension N×p) and Y (dimension N×q) respectively, to obtain... and The Pearson correlation coefficient matrix of dimension p×q is obtained by matrix multiplication in one step.

[0022] Where N is the number of samples.

[0023] Preferably, the comprehensive scoring function in step S5 is:

[0024] in, This is represented by the corrected gene level GWAS Z-score between gene i and microorganism j. Let be the R² aggregation function, and D be the directional consistency reward factor.

[0025] Preferably, the direction consistency reward factor D is calculated as follows: For a candidate ternary relation (i,j,k), the expected direction and the actual direction are defined:

[0026]

[0027]

[0028] in, As expected, For practical direction, Represented as a symbolic function, The Pearson correlation coefficient is expressed as the relationship between the abundance of microorganism j and phenotype k. It is expressed as the Pearson correlation coefficient between gene i expression level and phenotype k.

[0029] Preferably, the joint gain in step S6 is defined as:

[0030] ;

[0031] in, Let M be the predicted phenotypic value, where M is the predicted phenotypic value. G For a single-factor gene model, M M For microbial single-factor models, M GM For a joint model that includes genes, microorganisms, and their interactions, R² is the prediction determination coefficient.

[0032] Preferably, the formula for calculating the final comprehensive score in step S7 is:

[0033]

[0034] in, For gene-based GWAS statistics, For joint gain, This is an interaction importance factor.

[0035] Preferably, the method supports three calculation modes:

[0036] Ultra-fast mode: Only steps S1 to S5 are executed, using pure vectorized computation, and the R² aggregation function uses geometric mean. ;

[0037] Two-stage fast mode: Executes the complete process from steps S1 to S7, where step S6 uses simplified gradient boosting model parameters;

[0038] Full evaluation mode: Perform steps S1 to S7, where step S6 uses the full gradient boosting model parameters.

[0039] Preferably, the preset threshold θ of the threshold constraint in step S3 is set to 2.0, corresponding to a significance level of less than 0.05 for a two-sided test p-value.

[0040] The beneficial effects of this invention are as follows:

[0041] This invention provides a multi-omics integrated gene-microbe-phenotype ternary relationship mining method, which has the following advantages compared with the prior art: (1) Unified integration of multi-omics: For the first time, GWAS genetic association, gene expression, microbiome abundance and phenotype data are jointly analyzed within the same analytical framework; (2) SNP quantity bias correction: Eliminates systematic bias caused by differences in gene size; (3) Nonlinear relationship capture: Automatically identifies nonlinear effects and interaction effects using a gradient boosting model; (4) Mediation effect quantification: For the first time, a gene-microbe mediation effect quantification index based on machine learning is proposed; (5) Flexible multi-mode calculation: Provides three calculation modes to meet different needs from rapid exploration to in-depth analysis; (6) Directional consistency verification: Enhances the biological interpretability of the results. Attached Figure Description

[0042] The invention will now be further described with reference to the accompanying drawings.

[0043] Figure 1 The final result is a gene-microbe network diagram. Detailed Implementation

[0044] To make the technical means, creative features, objectives and effects of this invention easier to understand, the invention will be further described below in conjunction with specific embodiments.

[0045] Please see Figure 1 As shown in the figure, this invention provides a method for mining the gene-microorganism-phenotype ternary relationship, and the specific steps are as follows:

[0046] Step S1: Data Preparation and Preprocessing – Acquire and preprocess four types of omics data, including: phenotypic matrix P (dimension N×K, where N is the number of samples and K is the phenotype), microbiome abundance matrix M (dimension N×J, where J is the number of microorganisms), gene expression matrix G (dimension N×I, where I is the number of genes), and gene-microorganism association matrix Z (dimension I×J); where the element z_ij in the Z matrix represents the gene-based GWAS Z-score between the i-th gene and the j-th microorganism.

[0047] The above data undergoes sample matching, missing value imputation, and Z-score standardization. The specific process is as follows:

[0048] A common sample set is obtained by taking the intersection of the sample indices of the four types of data matrices; a common gene set is obtained by taking the intersection of the gene indices of the Z-score matrix and gene expression matrix based on gene-based GWAS; a common microbial set is obtained by taking the intersection of the microbial indices of the Z-score matrix and microbial abundance matrix based on gene-based GWAS; covariate columns such as age, sex, and batch are removed from the phenotypic data, and numerical phenotypes with at least 30 valid observations are retained; column Z-score standardization is performed on microbial abundance and gene expression data.

[0049] For example, a multi-omics dataset of a beef cattle population is used, containing the following four types of data:

[0050] (1) Phenotypic data: K phenotypic traits of N individuals are stored in a text file separated by tabs. The rows are individuals and the columns are traits. After excluding covariate columns such as age, sex and batch, valid numerical phenotypes are retained.

[0051] (2) Microbiome abundance data: relative abundance values ​​of J microbial taxa (species level) for N individuals;

[0052] (3) Gene expression data: the expression level of I genes in N individuals;

[0053] (4) Gene-microbe association matrix: Z-score value of the corrected association between I genes and J microorganisms.

[0054] Data preprocessing workflow: Find the intersection of common samples from the four data categories, ensuring the sample size is no less than 30; find the common gene set from the Z-score matrix and gene expression matrix; find the common microbial set from the Z-score matrix and microbial abundance matrix; for microbial abundance and gene expression data, first fill missing values ​​with the mean of each column, then perform column Z-score standardization (subtract the mean, divide by the standard deviation, and add 1 × 10⁻⁶ to the standard deviation). -8 (Avoid division by zero); missing values ​​in phenotypic data are replaced with NA.

[0055] Step S2: Construction of GWAS association signals at the gene level—First, based on the mixed linear model association analysis (LLM) method of GCTA software, each microbial abundance is used as a phenotype, and genome-wide association analysis is performed on each SNP; then, the mBAT-combo method is used to pool the SNP-level association P-values ​​to the gene level to obtain the gene-level association P-value for each gene-microbe pair; then, the P-values ​​are converted into Z-scores, and linear regression correction is performed with the logarithmic number of SNPs as the independent variable. The residuals are taken and standardized to obtain the corrected Z-score matrix that eliminates the bias of the number of SNPs.

[0056] The process for obtaining the gene-level association P-value for each gene-microbe pair is as follows:

[0057] The population structure and kinship were controlled by the genomic kinship matrix (GRM), and the abundance data of each microorganism was used as the phenotype to detect the association between each SNP site on the genome and the microbial abundance. The mBAT-combo method used the association statistics and linkage disequilibrium structure of all SNPs in the gene region to aggregate the SNP level signals into a gene-level comprehensive association P value through multivariate association analysis.

[0058] The methods for linear regression correction using the logarithmic number of SNPs as the independent variable specifically include:

[0059] The gene-level association P-values ​​output by mBAT-combo are converted to the original Z-score using a standard normal quantile function:

[0060]

[0061] in, The Z-score is the uncorrected original gene level, representing the statistical significance of the association between the gene and the microorganism. The inverse function (quantile function) of the cumulative distribution function of the standard normal distribution. P-value for gene-level association;

[0062] Logarithmic SNP count Using Zraw as the independent variable and Zraw as the dependent variable, perform ordinary least squares linear regression fitting:

[0063]

[0064] The residual ε is taken and Z-score standardized to obtain the corrected gene-level Z-score; where... This represents the total number of SNPs contained within a gene region. The intercept of the linear regression is... is the slope of the linear regression.

[0065] For example, GWAS association signal construction at the gene level:

[0066] SNP-level microbiome GWAS: For each microorganism (J species in total), its abundance value was used as the phenotype, and genome-wide association analysis was performed using the LMM method of GCTA software.

[0067]

[0068] in, It is the corrected phenotypic value vector. It is the correlation matrix of fixed effects. It is a fixed effects vector. It is the correlation matrix of the target set sites. Each conforms to ~ It follows a normal distribution.

[0069] Gene-level signal pooling: The mBAT-combo method in GCTA software is used as a basis to pool the association statistics of all SNPs in each gene region using a multivariate method. The running parameters include --mBAT-print-all-p to output the results for all genes.

[0070] Z-score transformation and correction: The P-value of the mBAT-combo output is first converted to the original Z-score using the upper tail quantile function of the standard normal distribution; for infinity values, qnorm(1×10) is used. -300 Replace with the value of ), and replace missing values ​​with 0. Then use For independent variable, Perform linear regression on the dependent variable, take the residuals and standardize them using the Z-score to obtain the corrected value. β0 and β1 are the intercept and slope of the linear regression model, and scale(·) is the standardization function. This correction eliminates the systematic bias of gene length (including the number of SNPs) on the association statistic:

[0071]

[0072] Step S3: GWAS constraint screening—Apply a threshold constraint |z_ij| ≥ θ (θ is a preset threshold, it is recommended to set it to 2, corresponding to a two-sided test with approximately P < 0.05) to the Z-score matrix Z, and screen out gene-microbe candidate pairs with statistically significant genetic associations;

[0073] It needs to be explained that a threshold constraint is applied to the corrected Z-score matrix. (In this example, θ=2.0) Significant gene-microbe pairs are preserved, greatly reducing the search space.

[0074] Step S4: Calculate the multi-level correlation matrix: Perform column standardization on the input matrices X (dimension N×p) and Y (dimension N×q) respectively to obtain... and The Pearson correlation coefficient matrix of dimension p×q is obtained by matrix multiplication in one step.

[0075]

[0076] When GPU acceleration is configured, the above operations are transferred to the GPU and performed using PyTorch tensor operations, which greatly accelerates the computational efficiency.

[0077] Specifically, vectorized matrix operations are used to calculate gene-phenotype correlation coefficient matrices in batches. (Dimension I×K), Microbial-Phenological Correlation Coefficient Matrix (Dimensions J×K), and gene-microbe expression correlation coefficient matrix (Dimension I×J); Simultaneously calculate the R² value matrix and directional consistency index corresponding to each correlation coefficient;

[0078] (a) Gene-phenotype correlation coefficient matrix (dimension I×K):

[0079]

[0080] (b) Microbial-phenotypic correlation coefficient matrix (dimension J×K):

[0081]

[0082] (c) Gene-microbe expression correlation coefficient matrix (dimension I×J):

[0083]

[0084] in, T represents the transpose of the standardized gene expression matrix. Represented as the transpose of the standardized microbial abundance matrix. It is represented as a standardized phenotypic matrix, where N represents the number of samples;

[0085] It needs to be explained that vectorized matrix operations are used to batch calculate the three types of correlation matrices and gene-phenotype correlation coefficient matrices. (Dimension I×K), Microbial-Phenological Correlation Coefficient Matrix (Dimensions J×K), and gene-microbe expression correlation coefficient matrix (Dimension I×J), completed in one step through standardized matrix multiplication. The Na values ​​in the result are replaced with 0, and the R² matrix is ​​calculated simultaneously.

[0086] Step S5: Candidate Triad Scoring and Preliminary Screening: For each gene-microbe pair that passes the constraints of Step S3... Candidate ternary relations are formed with all K phenotypic combinations. Calculate preliminary scores based on the comprehensive scoring function, and retain the Top-N1 candidates according to their scores; the scoring function is:

[0087]

[0088] in, This is represented by the corrected gene level GWAS Z-score between gene i and microorganism j. Let R² be the aggregation function, and D be the directional consistency reward factor;

[0089] The calculation method for the directional consistency reward factor D is as follows: For a candidate ternary relation (i,j,k), define the expected direction and the actual direction:

[0090]

[0091]

[0092] in, As expected, For practical direction, Represented as a symbolic function, The Pearson correlation coefficient is expressed as the relationship between the abundance of microorganism j and phenotype k. Expressed as the Pearson correlation coefficient between gene i expression level and phenotype k;

[0093] That is, if the two are the same, then D = 1.5; otherwise, D = 1.0.

[0094] In the candidate generation process, a batch processing and intermediate truncation strategy is adopted: significant gene-microbe pairs are divided into batches of fixed size and processed batch by batch. When the cumulative number of candidates exceeds a preset multiple of the target retention number, the top candidates are truncated according to their scores to control memory usage. Finally, after deduplication of all candidates, the Top-N1 candidates are retained to enter the fine evaluation in step S6.

[0095] It should be explained that a two-stage fast mode is used to calculate a comprehensive score for significant gene-microbe pairs and all phenotype combinations, retaining the top-5000 candidates. A directional consistency factor D is introduced into the score to verify whether the direction of the gene-microbe effect on the phenotype is consistent with the inferred direction of the "gene→microbe→phenotype" pathway. Batch processing and intermediate truncation strategies are used to control memory usage.

[0096] Step S6: Joint Modeling and Evaluation of Gradient Boosting Models – For the candidate ternary relations selected in Step S5, a gradient boosting decision tree model is used for detailed evaluation; for each candidate ternary relation… Three models were constructed separately, and the feature matrix of the joint model contained three columns: gene expression values. Microbial abundance value and interactive items Gene importance is obtained by averaging the feature importance of each fold model in cross-validation. Importance of microorganisms And the importance of interaction .

[0097] (a) Single-factor gene model:

[0098]

[0099] (b) Microbial single-factor model:

[0100]

[0101] (c) Gene-microbe joint model:

[0102]

[0103] The prediction R² of each model is evaluated using K-fold cross-validation, and the joint gain and feature importance are calculated; the joint gain is defined as:

[0104] ;

[0105] in, Let M be the predicted phenotypic value, where M is the predicted phenotypic value. G For a single-factor gene model, M M For microbial single-factor models, M GM For a joint model that includes genes, microorganisms, and their interactions, R² is the predictive determination coefficient;

[0106] For example, the top 5000 candidates were finely evaluated using XGboost. For each candidate, a single-factor gene model, a single-factor microbial model, and a joint model (including interaction terms) were constructed, and R² was evaluated using 3-fold cross-validation. Model hyperparameters included: the number of trees (n_estimators) set to 50-200, the maximum depth (max_depth) set to 3-5, and the learning rate (learning_rate) set to 0.05-0.2. Optimal parameters were automatically found using a Bayesian network.

[0107] Step S7: Final Overall Score and Result Output—Calculate the final overall score based on the joint model performance, mediation gain, and interaction effect strength. Output the Top-N2 ternary relationships and their detailed evaluation metrics, sorted by score. For prediction of interaction models , The importance of interaction model features; the final comprehensive score formula:

[0108]

[0109] in, For gene-based GWAS statistics, For joint gain, As an interaction importance factor;

[0110] Specifically, calculate the final comprehensive score, sort the scores in descending order, and output the top-500 ternary relations.

[0111]

[0112]

[0113]

[0114]

[0115] in, This is a gene-phenotype correlation coefficient matrix. This is a microbial-phenotypic correlation coefficient matrix. This is a gene-microbe expression correlation coefficient matrix. As a directional consistency factor, This is an interaction importance factor.

[0116] The results showed that the top-ranked ternary relationship, the combination of gene ENSBTAG00000008443 and the microorganism Ligilactobacillus_salivarius, significantly influenced the net meat yield, with a gene-based GWAS Z-score of 2.58, a combined model R² of 0.371, and a final composite score of 1.329. The top 10 ternary relationships involved several known microbial species associated with rumen fermentation and nutrient metabolism (such as Pseudodesulfovibrio_profundus, Bifidobacterium_scardovii, Staphylococcus_aureus, etc.), further validating the biological rationale of the method.

[0117] In summary, the gene-microbe-phenotype ternary relationship mining method of this invention supports three computation modes: (a) Ultra-fast mode: only steps S1 to S5 are executed, pure vectorized computation is used, and the geometric mean is used for R² aggregation in the scoring function.

[0118]

[0119] (a) Fast two-stage mode: Executes the complete process from S1 to S7, with simplified gradient boosting model parameters (50 trees, maximum depth of 3 layers, 3-fold cross-validation) in step S6, which can be completed in minutes; (b) Full evaluation mode: Executes steps from S1 to S7, with full gradient boosting model parameters (100 trees, maximum depth of 4 layers, learning rate of 0.1 and 5-fold cross-validation, including early stopping strategy) in step S6, which is suitable for scenarios that require the highest accuracy evaluation.

[0120] The above description is merely an embodiment of the present invention and does not limit the patent scope of the present invention. Any equivalent structural or procedural transformations made based on the content of the present invention specification, or direct or indirect applications in other related technical fields, are similarly included within the patent protection scope of the present invention.

[0121] The foregoing has shown and described the basic principles, main features, and advantages of the present invention. Those skilled in the art should understand that the present invention is not limited to the above embodiments. The embodiments and descriptions in the specification are merely illustrative of the principles of the invention. Various changes and modifications can be made to the invention without departing from its spirit and scope, and all such changes and modifications fall within the scope of the present invention as claimed. The scope of protection of the present invention is defined by the appended claims and their equivalents.

Claims

1. A method for mining the ternary relationship between genes, microorganisms, and phenotypes, characterized in that, Includes the following steps: Step S1: Acquire and preprocess four types of omics data, including phenotypic matrix P, microbiome abundance matrix M, gene expression matrix G, and gene-microbe association matrix Z; Step S2: Construct a gene-level GWAS association signal, pool the association P-values ​​at the SNP level to the gene level, obtain the gene-level association P-value for each gene-microbe pair, convert it to Z-score, and then perform linear regression correction with the number of SNPs in the gene region as the independent variable to obtain the corrected Z-score matrix that eliminates the bias of the number of SNPs. Step S3: Apply threshold constraints to the corrected Z-score matrix to screen out gene-microbe candidate pairs with statistically significant genetic associations; Step S4: Calculate the multi-level correlation matrix, including the gene-phenotype correlation coefficient matrix R, through vectorized matrix operations. GP Microbial-phenotypic correlation coefficient matrix R MP And the gene-microbe expression correlation coefficient matrix R GM ; Step S5: For each gene-microbe candidate pair obtained through step S3, form a candidate ternary relationship with each phenotype combination, calculate a preliminary score according to the comprehensive scoring function, and retain the Top-N1 candidates according to the score; the comprehensive scoring function includes gene-microbe GWAS Z-score, R² aggregation term, square term of gene-microbe expression correlation coefficient, and directional consistency reward factor. Step S6: For the candidate ternary relationships selected in Step S5, a gradient boosting decision tree model is used for joint modeling and evaluation. A gene single-factor model, a microbial single-factor model, and a joint model containing genes, microorganisms, and their interaction terms are constructed respectively. The joint gain and the importance of interaction features are calculated. Step S7: Calculate the final comprehensive score based on the joint model performance, mediation gain, and interaction effect strength, and output the Top-N2 ternary relations in order of score.

2. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, In step S2, the association P-values ​​at the SNP level are pooled to the gene level using the mBAT-combo method, and multivariate association analysis is performed using the association statistics and linkage disequilibrium structure of all SNPs in the gene region.

3. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, The specific steps for linear regression correction in step S2 are as follows: The gene-level association P-value is converted to the original Z-score using the standard normal quantile function, and the score is calculated using the logarithmic SNP count. Using the original Z-score as the dependent variable and performing linear regression with the original Z-score as the independent variable, the regression residuals are taken and standardized with Z-score to obtain the corrected gene-level Z-score.

4. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, Step S4 involves calculating the multi-level correlation matrix: The input matrix X (dimension N×p) and matrix Y (dimension N×q) are column-normalized to obtain... and The Pearson correlation coefficient matrix of dimension p×q is obtained by matrix multiplication in one step. Where N is the number of samples.

5. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, The comprehensive scoring function mentioned in step S5 is: in, This is represented by the corrected gene level GWAS Z-score between gene i and microorganism j. Let be the R² aggregation function, and D be the directional consistency reward factor.

6. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 5, characterized in that, The calculation method for the directional consistency reward factor D is as follows: For a candidate ternary relation (i,j,k), define the expected direction and the actual direction: ； in, As expected, For practical direction, Represented as a symbolic function, The Pearson correlation coefficient is expressed as the relationship between the abundance of microorganism j and phenotype k. It is expressed as the Pearson correlation coefficient between gene i expression level and phenotype k.

7. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, The joint gain described in step S6 is defined as follows: ； in, Let M be the predicted phenotypic value, where M is the predicted phenotypic value. G For a single-factor gene model, M M For microbial single-factor models, M GM For a joint model that includes genes, microorganisms, and their interactions, R² is the prediction determination coefficient.

8. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, The formula for calculating the final comprehensive score in step S7 is as follows: in, For gene-based GWAS statistics, For joint gain, This is an interaction importance factor.

9. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, The method supports three calculation modes: Ultra-fast mode: Only steps S1 to S5 are executed, using pure vectorized computation, and the R² aggregation function uses geometric mean. ; Two-stage fast mode: Executes the complete process from steps S1 to S7, where step S6 uses simplified gradient boosting model parameters; Full evaluation mode: Perform steps S1 to S7, where step S6 uses the full gradient boosting model parameters.

10. The method for mining the gene-microorganism-phenotype ternary relationship according to claim 1, characterized in that, In step S3, the preset threshold θ of the threshold constraint is set to 2.0, which corresponds to a significance level of less than 0.05 for a two-sided test p-value.

Images