Multi-character whole-genome associated mixed linear model

Abstract

The invention discloses a multi-character whole-genome associated mixed linear model method MTMMade. The model is based on a multi-character mixed model, dominant and epistatic effects are innovatively introduced, associated variation sites of genotypes and phenotypes are comprehensively analyzed, and the detection capability of non-additive effect sites and the analysis precision of complex character heritable variation are improved. Through simulation and real rice data verification, the model significantly enhances the joint detection capability of the non-additive effect, effectively controls the false positive rate, and ensures the reliability of the result. The model is packaged into an R language statistical package, a user only needs to input genotype and phenotype data, analysis can be completed through a row of simplified codes, dependence on a professional background is greatly reduced, and non-professional users can conveniently carry out multi-character whole genome association analysis. The model and a software tool thereof provide important support for animal and plant genetic research and crop genetic improvement, and provide an accurate theoretical basis and a practical method for complex character analysis and breeding strategy optimization.

Description

Technical Field

[0001] The present invention belongs to the technical field of quantitative genetics, and particularly relates to a mixed linear model for multi-trait genome-wide association. Background Art

[0002] Genotype non-additive effects, that is, the influence of gene-gene interactions on traits, play a key role in the genetic regulation of many complex traits. In recent years, genome-wide association studies (GWAS) have gradually begun to focus on the importance of non-additive effects while revealing the additive effects of traits. [1] . Traditional GWAS models mainly focus on the additive effects of markers. However, more and more studies have shown that non-additive effects, such as dominance effects and epistatic effects, play an important role in complex traits. [2] . For example, Reynolds et al. [3] studied the lactation traits of dairy cows through non-additive GWAS and successfully identified 5 new quantitative trait loci (QTLs) with significant recessive effects on three milk production traits. Similarly, Xue et al. [4] identified 805 epistatic effect loci related to average daily gain, adjusted age, and loin muscle depth in 2360 Duroc pigs by simultaneously fitting additive and non-additive effect loci. These studies fully demonstrate the influence of non-additive effects on complex traits.

[0003] Non-additive effects are also considered an important factor in "missing heritability". [5] . In-depth study of non-additive effects not only helps the present application to more comprehensively understand the genetic basis of complex traits, but also provides new ideas for breeding applications. By integrating non-additive effect information, breeders can more effectively select and combine ideal genotypes, thereby improving breeding efficiency. [1] . The research of Li et al. [6] showed that compared with inbred rice, improving the grain quality of hybrid rice is more challenging, which is mainly attributed to the existence of non-additive effects such as dominance. They identified 44 additive effect loci, 97 dominance effect loci, and 13 loci with both additive and dominance effects in 12 rice grain quality traits through joint analysis of additive and non-additive effects, and these loci can explain more than 30% of the genetic variation of each trait and improve the accuracy of genomic prediction. These results all emphasize the importance of non-additive effects in the analysis of complex traits and genetic improvement.

[0004] In recent years, multi-trait GWAS methods have shown great potential compared with single-trait GWAS. For example, the PC-GWAS method can better reveal hidden multi-effect loci and significantly improve the detection ability by integrating the principal components of multiple related traits. [7]。The mvLMM method can simultaneously consider multiple traits and their correlations, providing more accurate association analysis results. [8] 。These methods have made significant progress in detection power and analysis speed and have been applied to the study of complex traits. However, many multi-trait GWAS methods often neglect the interpretability of multi-effect loci. For example, methods such as mvLMM, TATES, and PC-GWAS all face this problem. Although these methods improve the detection power, the increase in their model complexity often leads to a decrease in the interpretability of multi-effect loci. The MTMM method pays more attention to analyzing the effects of multi-effect loci, decomposing them into co-directional effects and counter-directional effects. However, its limitation is that it can only analyze two traits. The MTMM model has been successfully applied to the study of human blood lipid-related traits, successfully identifying multiple multi-effect loci that can simultaneously affect multiple blood lipid metabolism traits (such as triglycerides and low-density lipoproteins). [7] 。However, traditional MTMM models mainly focus on additive effects and neglect the important contributions of non-additive effects (such as dominance effects and epistatic effects) to trait performance. Although existing studies have shown that dominance effects play an important role in trait genetic variation, existing models (such as MTMM) have not fully exploited this information. In addition, epistatic effects (i.e., interactions between genes) have also been confirmed to be important factors affecting trait inheritance, but their analysis remains a difficult point in current research. Therefore, developing a model that can comprehensively consider additive effects, dominance effects, and epistatic effects is crucial for a more comprehensive understanding of the genetic mechanisms of complex traits and improving breeding efficiency. Summary of the Invention

[0005] The object of the present invention is to provide a mixed linear model for multi-trait genome-wide association, so that this method can simultaneously consider additive effects, dominance effects, and epistatic effects, thereby more comprehensively and accurately analyzing the genetic variation of complex traits.

[0006] To solve the above technical problems, the specific technical solutions of the present invention are as follows:

[0007] In some embodiments of the present application, a mixed linear model MTMMade for multi-trait genome-wide association is provided. This model adds a comprehensive consideration of dominance effects and epistatic effects on the basis of the MTMM model. At the same time, combined with a false positive filtering mechanism, through the verification of simulated data and real rice data, the MTMMade model significantly improves the detection accuracy and statistical power, specifically including the following steps:

[0008] Step 1: Data quality control, screening, checking, imputing, and quality control processing of phenotypic data and genotypic data;

[0009] Step 2: Calculate the genetic variance components associated with the phenotype and integrate and process the results by using a multi-trait mixed model;

[0010] Step 3: construct a multi-effect genotype matrix;

[0011] Step 4: Multi-effect multi-trait genome-wide association.

[0012] In some embodiments of the present application, step 1 includes the following steps:

[0013] Step 1.1) Sample name consistency filtering, that is, comparing the sample names in the genotype and phenotype data to ensure a complete match. Any samples with missing values ​​in the genotype or phenotype data will be deleted to generate a matched and non-missing genotype-phenotype dataset;

[0014] Step 1.2) Genotype data check and imputation: Check whether the genotype data is in PLINK additive coding format, where 0 is homozygous reference, 1 is heterozygous, and 2 is homozygous variant;

[0015] Missing genotype data will be filled with heterozygous genotypes, coded as 1;

[0016] Step 1.3) Minimum allele frequency quality control: traverse the variant sites in the genotype data, calculate the MAF, and eliminate variant sites with a MAF lower than 0.05.

[0017] In some embodiments of the present application, step 2 includes the following steps:

[0018] Step 2.1) Input data includes sample ID, two phenotypic values, denoted as Y1 and Y2, and the kinship matrix K between samples; first, Y1 and Y2 are merged into a mixed phenotype vector Y_ok, and an environmental factor vector Env is constructed to distinguish the data sources of Y1 and Y2; then, the kinship matrix K is standardized to obtain a standardized kinship matrix K_stand;

[0019] Step 2.2) Preliminary evaluation of heritability: Use linear mixed models to perform preliminary heritability evaluation on phenotypes Y1 and Y2 respectively; the model forms of Y1 and Y2 are as follows:

[0020] Y1=μ1+g1+e1

[0021] Y2=μ2+g2+e2

[0022] where μ1 and μ2 are the population means of Y1 and Y2 respectively, g1 and g2 are the genetic effects of Y1 and Y2 respectively, and e1 and e2 are the residuals of Y1 and Y2 respectively; the genetic effects are assumed to follow normal distributions of g1 ~ N(0, σ_g1^2 * K_stand) and g2 ~ N(0, σ_g2^2 * K_stand); through this model, the heritabilities of Y1 and Y2 are estimated, denoted as herit1 and herit2 respectively;

[0023] Step 2.3) Multi-trait mixed model: Use the multi-trait mixed model to simultaneously analyze the genetic and environmental effects of Y1 and Y2; the model form is as follows:

[0024] Y_ok = Xb + Zu + e

[0025] where Y_ok is the mixed phenotype vector, X is the design matrix containing the environmental factor Env, b is the fixed effect, Z is the design matrix, u is the random effect, including the genetic effect, which follows a normal distribution, and its variance-covariance structure is denotes the Kronecker product, varcov is the genetic variance-covariance matrix, e is the residual, and its variance-covariance structure is ve is the residual variance-covariance matrix; through this model, the genetic variance-covariance matrix varcov and the residual variance-covariance matrix ve are estimated;

[0026] Step 2.4) Correlation calculation: Based on the varcov and ve matrices, calculate the genetic correlation coefficient rho_g, environmental correlation coefficient rho_e, phenotypic correlation coefficient rho_p, and Pearson correlation coefficient pearson between phenotypes Y1 and Y2; at the same time, calculate the genetic covariance G, specific genetic variances G1 and G2 of single traits, environmental variance E, and environmental covariance EE; and calculate the proportions of each variance component, including the proportion of genetic variance var_G, proportion of specific genetic variance var_GE, proportion of environmental variance var_E, and proportion of environmental covariance var_EE;

[0027] Step 2.5) Significance assessment: Use the Delta method to calculate the standard errors of the genetic correlation coefficient rho_g and environmental correlation coefficient rho_e, denoted as se_g and se_e respectively; at the same time, calculate the corresponding p-values through the likelihood ratio test, denoted as p_g and p_e respectively;

[0028] Step 2.6) Result integration and data preparation: Integrate the results such as the genetic correlation coefficient rho_g, the environmental correlation coefficient rho_e, the standard errors se_g and se_e, the p-values p_g and p_e, the heritabilities herit1 and herit2, and the variance component ratios var_G, var_GE, var_E, and var_EE into the correlation list; also add the mixed model coefficient matrix M and the adjusted phenotypic data Y_t to the list; the correlation list is the genetic variance result of Y1 and Y2.

[0029] In some embodiments of the present application, step 3 includes the following steps:

[0030] Step 3.1) Construct the main effect gene matrix: Convert the quality-controlled genotype data into a main effect matrix; the genotype data adopts an additive coding method, where 0 represents the homozygous reference genotype, 1 represents the heterozygous genotype, and 2 represents the homozygous variant genotype; at the same time, perform dominant effect coding on the genotype data, coding both the homozygous reference and the homozygous variant as 0, and the heterozygous genotype as 1; combine the additive coding matrix and the dominant coding matrix to form a main effect gene matrix containing additive effects and dominant effects for subsequent identification of variant sites.

[0031] Step 3.2) Construct the epistatic effect gene matrix: Construct the epistatic effect gene matrix through non-repetitive permutation and combination; the types of epistatic effects include four types: additive × additive, additive × dominant, dominant × additive, and dominant × dominant.

[0032] Taking the epistatic effect type of additive × additive as an example, for each variant site in the genotype data, select any two additive sites for combination; for example, if there are n additive sites, the number of additive × additive combination methods is n(n - 1) / 2; for each pair of additive sites, multiply their additive genotype values to construct the corresponding additive × additive epistatic site; for the genotype values of epistatic effect types such as additive × dominant, dominant × additive, and dominant × dominant, use a similar method for construction; finally, the number of constructed epistatic sites will depend on different epistatic effect types and the number of sites in the genotype data.

[0033] Step 3.3) Quality control of the epistatic gene matrix: Since the dimension of the epistatic effect gene matrix grows exponentially with permutation and combination, to ensure the reliability of subsequent analysis, it is necessary to perform minor allele frequency filtering on the epistatic effect genotype matrix; the specific operation is the same as step 1.3, traverse the variant sites in the epistatic genotype matrix, calculate their minor allele frequencies, and filter out the variant sites with MAF less than 0.05.

[0034] In some embodiments of the present application, step 4 includes the following steps:

[0035] Step 4.1) Input data preparation: Use the multi-effect genotype matrix constructed in Step C, including additive effect, dominant effect, and epistatic effect, as the input. At the same time, input the adjusted phenotypic data Y_t, covariate matrix cof_t, mixed model coefficient matrix M, environmental factor vector Env, genetic variance-covariance matrix varcov, and sample ID obtained in Step 2. In addition, it also includes SNP_INFO containing variant site information.

[0036] Step 4.2) Effect scaling: Scale the multi-effect genotype matrix. Multiply the parts of the multi-effect genotype matrix related to phenotypes Y1 and Y2 by the corresponding genetic standard deviations in the varcov matrix, denoted as X_ok1 and X_ok2. Combine X_ok1 and X_ok2 row by row to obtain X_.

[0037] Step 4.3) Construct the transformed genotype matrix: Use the mixed model coefficient matrix M to transform the scaled genotype matrix X_ into X_t. X_t is a three-dimensional array, where X_t[,,1] is the matrix related to the total effect, and its calculation formula is X_t[,,1] = M %*% X_. X_t[,,2] is the matrix related to the genotype-by-environment interaction effect, and its calculation formula is X_t[,,2] = M %*% (Env * X_).

[0038] Step 4.4) Calculate the common effect: Use a linear model to analyze the common effect of each variant site. The model form is as follows:

[0039] Y_t = cof_t * β_cof + X_t[,,1] * β_com + ε

[0040] where Y_t is the adjusted phenotypic data, cof_t is the covariate matrix, X_t[,,1] is the part related to the total effect in the transformed genotype matrix, β_cof and β_com are the effect values of the covariate effect and the common effect respectively, and ε is the residual. Through linear regression, obtain the p-value pval_com, effect value beta_com, and standard error se_com of the common effect of each variant site.

[0041] Step 4.5) Calculate the interaction effect: Use a linear model to analyze the genotype-by-environment interaction effect of each variant site. The model form is as follows:

[0042] Y_t = cof_t * β_cof + X_t[,,1] * β_com + X_t[,,2] * β_inter

[0043] + ε

[0044] Among them, \(X_t[,,2]\) is the part related to the genotype-by-environment interaction effect in the transformed genotype matrix, and \(\beta_{inter}\) is the effect value of the interaction effect; through linear regression, the p-value \(pval_{inter}\), effect value \(\beta_{inter}\), and standard error \(se_{inter}\) of the interaction effect of each variant locus are obtained.

[0045] Step 4.6) Calculate the Full Model Effect: Using the likelihood ratio test method, compare the full model that includes all effects (main effects and interaction effects) with the model that only includes the covariate effects to evaluate the overall effect of each variant locus on the phenotype; use the least squares method to fit the model that includes covariates and the model that includes covariates and genotypes, calculate the residual sum of squares, and use the F-test to calculate the p-value; the calculation formula is as follows:

[0046] \(RSS_{full} = sum(lsfit(cbind(cof_t, x), Y_t, intercept =

[0047] FALSE)$residuals^2)\) is the residual sum of squares of the model that includes covariates and genotype effects, where \(x\) represents the current variant locus.

[0048] \(RSS_{env} = sum(lsfit(cof_t, Y_t, intercept =

[0049] FALSE)$residuals^2)\) is the residual sum of squares of the model that only includes covariates.

[0050] \(par_{env} = ncol(cof_t)\) is the degrees of freedom of the model that only includes covariates

[0051] \(par_{full} = par_{env}+2\) is the degrees of freedom of the model that includes covariates and genotypes

[0052] \(F_{full} = ((RSS_{env} / RSS_{full}-1)*(2*n -

[0053] \(par_{full}) / (par_{full}-par_{env}))\) calculates the F statistic, where \(n\) represents the sample size

[0054] \(pval_{full} = pf(F_{full}, par_{full}-par_{env}, 2*n-par_{full}, lower.tail = FALSE)\) calculates the p-value.

[0055] In some embodiments of the present application, the core mathematical expression of the MTMMade model is as follows:

[0056]

[0057] where y i is the observed value of the i-th trait; s i is a 0-1 vector, which is 1 for all values belonging to the i-th trait and 0 otherwise; μ i is the phenotypic mean of the i-th trait; A j represents the j-th additive effect SNP, D j represents the j-th dominant effect SNP, E j represents the j-th epistatic interaction effect SNP-SNP pair; here E j can be further divided into aa, ad, da, dd interactions; the model in the case of only considering interactions is as follows:

[0058]

[0059] where aa, ad, da, and dd represent the four epistatic effect loci of additive-additive, additive-dominant, dominant-additive, and dominant-dominant respectively; β is the respective effect value.

[0060] Compared with the prior art, the beneficial effects of the present invention are that

[0061] The MTMMade method conducts genome-wide association analysis by constructing a multi-trait mixed linear model that comprehensively considers additive effects, dominant effects, and epistatic effects, and can more comprehensively and accurately analyze the genetic variation mechanism of complex traits. This method can not only effectively identify the variant loci related to traits, but also deeply explore the roles of non-additive effects (such as dominant effects and epistatic effects) in the formation of complex traits, providing a strong theoretical basis and practical method for plant and animal genetic research and crop genetic improvement. Through this method, researchers and breeders can more accurately understand the complex relationship between genotypes and phenotypes, thereby providing key support for formulating more effective breeding strategies and achieving precise genetic improvement. In addition, encapsulating this method into a user-friendly R language statistical package significantly reduces the usage threshold for non-professional users, making multi-trait genome-wide association analysis easier to popularize and apply. BRIEF DESCRIPTION OF THE DRAWINGS

[0062] By reading the detailed description of the preferred embodiments below, various other advantages and benefits will become clear to those of ordinary skill in the art. The drawings are only for the purpose of showing the preferred embodiments and are not considered to be a limitation of the present invention. Moreover, throughout the drawings, the same reference numerals are used to represent the same components. In the drawings:

[0063] Figure 1 is a schematic diagram of the simulation experiment comparison between the model MTMMade and MTMM provided by the embodiment of the present invention;

[0064] Figure 2 Schematic diagram of the MTMMade and MTMM comparison results (Full) for 575 rice hybrid materials provided by the embodiments of the present invention;

[0065] Figure 3 Schematic diagram of the MTMMade and MTMM comparison results (Common) for 575 rice hybrid materials provided by the embodiments of the present invention. Detailed implementation manners

[0066] The following combines the accompanying drawings and embodiments to further describe in detail the specific implementation manners of the present invention. The following embodiments are used to illustrate the present invention, but do not limit the scope of the present invention.

[0067] To better understand the purpose, structure and function of the present invention, the following further describes the present invention in detail with reference to the accompanying drawings.

[0068] Embodiment 1

[0069] In this embodiment, on the MTMM model, a comprehensive consideration of the dominance effect and the epistatic effect is newly added. At the same time, combined with the false positive filtering mechanism, through the verification of simulated data and real rice data, the MTMMade model significantly improves the detection accuracy and statistical power. Among them, the core mathematical expression of the MTMMade model is as follows:

[0070]

[0071] Where y i is the observed value of the i-th trait; s i is a 0-1 vector, which is 1 for all values belonging to the i-th trait and 0 otherwise; μ i is the phenotypic mean of the i-th trait; A j represents the j-th additive effect SNP, D j represents the j-th dominance effect SNP, E j represents the j-th epistatic interaction effect SNP-SNP pair; here E j can be further divided into aa, ad, da, dd interactions; the model in the case of only considering interactions is as follows:

[0072]

[0073] Where aa, ad, da and dd represent the four epistatic effect loci of additive-additive, additive-dominance, dominance-additive, and dominance-dominance respectively; β is the respective effect value.

[0074] Specifically, it includes the following steps:

[0075] Step 1: Data quality control, screening, checking, filling and quality control of the comparison gene and genotype data, including the following steps:

[0076] Step 1.1) Sample name consistency filtering, that is, comparing the sample names in the genotype and phenotype data to ensure a complete match. Any samples with missing values ​​in the genotype or phenotype data will be deleted to generate a matched and non-missing genotype-phenotype dataset;

[0077] Step 1.2) Genotype data check and imputation: Check whether the genotype data is in PLINK additive coding format, where 0 is homozygous reference, 1 is heterozygous, and 2 is homozygous variant;

[0078] Missing genotype data will be filled with heterozygous genotypes, coded as 1;

[0079] Step 1.3) Minimum allele frequency quality control: traverse the variant sites in the genotype data, calculate the MAF, and eliminate variant sites with a MAF lower than 0.05.

[0080] Step 2: Calculate the genetic variance components associated with the phenotype and integrate and process the results by using a multi-trait mixed model, including the following steps:

[0081] Step 2.1) Input data includes sample ID, two phenotypic values, denoted as Y1 and Y2, and the kinship matrix K between samples; first, Y1 and Y2 are merged into a mixed phenotype vector Y_ok, and an environmental factor vector Env is constructed to distinguish the data sources of Y1 and Y2; then, the kinship matrix K is standardized to obtain a standardized kinship matrix K_stand;

[0082] Step 2.2) Preliminary evaluation of heritability: Use linear mixed models to perform preliminary heritability evaluation on phenotypes Y1 and Y2 respectively; the model forms of Y1 and Y2 are as follows:

[0083] Y1=μ1+g1+e1

[0084] Y2=μ2+g2+e2

[0085] Among them, μ1 and μ2 are the population means of Y1 and Y2, g1 and g2 are the genetic effects of Y1 and Y2, e1 and e2 are the residuals of Y1 and Y2; the genetic effects are assumed to obey the normal distribution of g1~N(0,σ_g1^2*K_stand) and g2~N(0,σ_g2^2*K_stand); through this model, the heritability of Y1 and Y2 is estimated, which are recorded as herit1 and herit2 respectively;

[0086] Step 2.3) Multi-trait mixed model: Use the multi-trait mixed model to simultaneously analyze the genetic and environmental effects of Y1 and Y2; the model form is as follows:

[0087] Y_ok = Xb + Zu + e

[0088] where Y_ok is the mixed phenotype vector, X is the design matrix containing the environmental factor Env, b is the fixed effect, Z is the design matrix, u is the random effect, including the genetic effect, which follows a normal distribution, and its variance-covariance structure is denotes the Kronecker product, varcov is the genetic variance-covariance matrix, e is the residual, and its variance-covariance structure is ve is the residual variance-covariance matrix; through this model, the genetic variance-covariance matrix varcov and the residual variance-covariance matrix ve are estimated;

[0089] Step 2.4) Correlation calculation: Based on the varcov and ve matrices, calculate the genetic correlation coefficient rho_g, environmental correlation coefficient rho_e, phenotypic correlation coefficient rho_p, and Pearson correlation coefficient pearson between phenotypes Y1 and Y2; at the same time, calculate the genetic covariance G, specific genetic variances G1 and G2 of single traits, environmental variance E, and environmental covariance EE; and calculate the proportions of each variance component, including the proportion of genetic variance var_G, proportion of specific genetic variance var_GE, proportion of environmental variance var_E, and proportion of environmental covariance var_EE;

[0090] Step 2.5) Significance evaluation: Use the Delta method to calculate the standard errors of the genetic correlation coefficient rho_g and environmental correlation coefficient rho_e, denoted as se_g and se_e respectively; at the same time, calculate the corresponding p-values through the likelihood ratio test, denoted as p_g and p_e respectively;

[0091] Step 2.6) Result integration and data preparation: Integrate the results such as the genetic correlation coefficient rho_g, environmental correlation coefficient rho_e, standard errors se_g and se_e, p-values p_g and p_e, heritabilities herit1 and herit2, and the proportions of variance components var_G, var_GE, var_E, and var_EE into the correlation list; also add the mixed model coefficient matrix M and the adjusted phenotypic data Y_t to the list; the correlation list is the genetic variance result of Y1 and Y2;

[0092] Step 3: Construct a multi-effect genotype matrix. Step 3 includes the following steps:

[0093] Step 3.1) Construct the main effect gene matrix: Transform the quality-controlled genotype data into a main effect matrix; the genotype data adopts an additive coding method, where 0 represents the homozygous reference genotype, 1 represents the heterozygous genotype, and 2 represents the homozygous variant genotype; at the same time, perform a dominant effect coding on the genotype data, coding both the homozygous reference and the homozygous variant as 0, and the heterozygous genotype as 1; combine the additive coding matrix and the dominant coding matrix to form a main effect gene matrix containing additive effects and dominant effects for subsequent identification of variant sites;

[0094] Step 3.2) Construct the epistatic effect gene matrix: Construct the epistatic effect gene matrix through non-repetitive permutation and combination; the types of epistatic effects include four types: additive×additive, additive×dominant, dominant×additive, and dominant×dominant;

[0095] Taking the epistatic effect type of additive×additive as an example, for each variant site in the genotype data, select any two additive sites for combination; for example, if there are n additive sites, the number of additive×additive combination methods is n(n - 1) / 2; for each pair of additive sites, multiply their additive effect values to construct the corresponding additive×additive epistatic site; for epistatic effect types such as additive×dominant, dominant×additive, and dominant×dominant, use a similar method for construction; finally, the number of constructed epistatic sites will depend on different epistatic effect types and the number of sites in the genotype data;

[0096] Step 3.3) Quality control of the epistatic gene matrix: Since the dimension of the epistatic effect gene matrix grows exponentially with permutation and combination, to ensure the reliability of subsequent analysis, it is necessary to perform a minor allele frequency filter on the epistatic effect gene matrix; the specific operation is the same as in Step 1.3, traverse the variant sites in the epistatic gene matrix, calculate their minor allele frequencies, and filter out the variant sites with MAF less than 0.05.;

[0097] Step 4: Genome-wide association of multiple effects and multiple traits. The following steps are included in Step 4:

[0098] Step 4.1) Input data preparation: Use the multi-effect genotype matrix constructed in Step C, including additive effects, dominant effects, and epistatic effects as input, and at the same time input the adjusted phenotype data Y_t, covariate matrix cof_t, mixed model coefficient matrix M, environmental factor vector Env, genetic variance-covariance matrix varcov, and sample ID obtained in Step 2; in addition, SNP_INFO containing variant site information;

[0099] Step 4.2) Effect scaling: Scale the multi-effect genotype matrix; Multiply the parts of the multi-effect genotype matrix related to phenotypes Y1 and Y2 by the corresponding genetic standard deviations in the varcov matrix, denoted as X_ok1 and X_ok2 respectively; Combine X_ok1 and X_ok2 row by row to obtain X_;

[0100] Step 4.3) Construct the transformed genotype matrix: Use the mixed model coefficient matrix M to transform the scaled genotype matrix X_ into X_t; X_t is a three-dimensional array, where X_t[,,1] is the matrix related to the total effect, and the calculation formula is X_t[,,1] = M %*% X_, and X_t[,,2] is the matrix related to the genotype-by-environment interaction effect, and the calculation formula is X_t[,,2] = M %*% (Env * X_);

[0101] Step 4.4) Calculate the common effect: Use a linear model to analyze the common effect of each variant site; The model form is as follows:

[0102] Y_t = cof_t * β_cof + X_t[,,1] * β_com + ε

[0103] where Y_t is the adjusted phenotype data, cof_t is the covariate matrix, X_t[,,1] is the part of the transformed genotype matrix related to the total effect, β_cof and β_com are the effect values of the covariate effect and the common effect respectively, and ε is the residual; Through linear regression, obtain the p-value pval_com, effect value beta_com, and standard error se_com of the common effect of each variant site;

[0104] Step 4.5) Calculate the interaction effect: Use a linear model to analyze the genotype-by-environment interaction effect of each variant site; The model form is as follows:

[0105] Y_t = cof_t * β_cof + X_t[,,1] * β_com + X_t[,,2] * β_inter

[0106] + ε

[0107] where X_t[,,2] is the part of the transformed genotype matrix related to the genotype-by-environment interaction effect, and β_inter is the effect value of the interaction effect; Through linear regression, obtain the p-value pval_inter, effect value beta_inter, and standard error se_inter of the interaction effect of each variant site;

[0108] Step 4.6) Calculate the Full Model Effect: Using the likelihood ratio test method, compare the full model that includes all effects (main effects and interaction effects) with the model that only includes the covariate effects to evaluate the overall effect of each variant site on the phenotype; use the least squares method to fit the model that includes covariates and the model that includes covariates and genotypes, calculate the residual sum of squares, and use the F-test to calculate the p-value; the calculation formula is as follows:

[0109] RSS_full = sum(lsfit(cbind(cof_t, x), Y_t, intercept =

[0110] FALSE)$residuals^2) is the residual sum of squares of the model that includes covariates and genotype effects, where x represents the current variant site;

[0111] RSS_env = sum(lsfit(cof_t, Y_t, intercept =

[0112] FALSE)$residuals^2) is the residual sum of squares of the model that only includes covariates;

[0113] par_env = ncol(cof_t) is the degrees of freedom of the model that only includes covariates

[0114] par_full = par_env + 2 is the degrees of freedom of the model that includes covariates and genotype

[0115] F_full = ((RSS_env / RSS_full - 1) * (2 * n -

[0116] par_full) / (par_full - par_env)) calculates the F statistic, where n represents the sample size

[0117] pval_full = pf(F_full, par_full - par_env, 2 * n - par_full, lower.tail = FALSE) calculates the p-value.

[0118] Through the above technical solutions, the technical effects generated in the embodiments of the present application are:

[0119] The MTMMade method can more comprehensively and accurately analyze the genetic variation mechanism of complex traits by constructing a multi-trait mixed linear model that comprehensively considers additive effects, dominance effects, and epistatic effects and combining it with genome-wide association analysis. This method can not only effectively identify the variant sites related to traits but also deeply explore the roles of non-additive effects (such as dominance effects and epistatic effects) in the formation of complex traits, providing a strong theoretical basis and practical methods for plant and animal genetic research and crop genetic improvement. Through this method, researchers and breeders can more accurately understand the complex relationship between genotypes and phenotypes, thus providing key support for formulating more effective breeding strategies and achieving precision genetic improvement.

[0120] Example 2

[0121] To facilitate the use of the multi-effect multi-trait genome-wide association analysis method by non-bioinformatics professionals and improve research efficiency, this application has encapsulated and developed an R language statistical package based on the MTMMade model. This R package aims to provide an easy-to-use, powerful, and highly flexible tool that enables users to perform multi-trait genome-wide association analysis without in-depth understanding of complex models and algorithms.

[0122] The main functions of the R package include:

[0123] 1. Data input and preprocessing:

[0124] 1.1 Sample name consistency check: Automatically check whether the sample names in the input genotype data and phenotype data are exactly matched, and remove the missing samples to ensure data consistency.

[0125] 1.2 Data format check:

[0126] a) Genotype data: Accept genotype data in PLINK additive raw format (numerical data encoded by 012), and support users to directly provide the raw file path or the loaded R variable.

[0127] b) Genotype coordinate file: Support users to provide a coordinate file containing variant site names, chromosomes, and base positions, or automatically generate this file according to the raw format.

[0128] c) Phenotype data: Accept a phenotype file containing sample names and multiple phenotype values.

[0129] 1.3 Quality control: Provide flexible quality control options, allowing users to customize the minimum allele frequency (MAF) filtering threshold and select whether to perform filtering.

[0130] 2. Kinship matrix calculation and standardization:

[0131] 2.1 Kinship matrix calculation: By default, the VanRaden method is used to calculate the kinship matrix, and multiple calculation methods such as EMMA and transposed dot product matrix are supported.

[0132] 2.2 Kinship matrix standardization: Allows users to choose whether to standardize the kinship matrix to meet different analysis requirements.

[0133] 3. Model establishment and parameter setting:

[0134] 3.1 Model type selection:

[0135] a) According to the user input data, a multi-effect multi-trait mixed linear model is automatically established, including additive effects and dominant effects.

[0136] b) Supports users to select analysis modes, including main effect analysis (main), epistatic effect analysis (epistasis), and simultaneous main effect and epistatic analysis (all).

[0137] 3.2 Flexible parameter setting: Users can adjust the analysis method through parameters, such as whether to standardize the kinship matrix, set the MAF filtering threshold, select the epistatic effect model (AA, AD, DA, DD), etc.

[0138] 3.3 Estimated parameter import: Allows users to import externally estimated variance components, or automatically calculate variance components through the model and save them for subsequent analysis.

[0139] 4. Result output:

[0140] 4.1 Output the significance p-value of each significant locus. The final result will be integrated with the genotype coordinate file (X.map) provided by the user. At the same time, it is possible to choose to sort by chromosome and position, or not to sort. By selecting different methods, additive, dominant, or epistatic effect results can be output.

[0141] 4.2 Result log record: Using verbose = TRUE can output the analysis log.

[0142] Main function of the R package (MTMMade):

[0143] The core function of this R package is MTMMade, which integrates functions such as kinship matrix calculation, variance component estimation, and multi-trait mixed model analysis, providing users with a flexible and efficient workflow. The main parameters of the MTMMade function are as follows:

[0144] · X or X.raw.path: Genotype data matrix or PLINK raw file path.

[0145] · Y: Phenotypic data, including sample IDs and multiple phenotypic values.

[0146] · K: Optional kinship matrix.

[0147] · X.map: Optional SNP information file.

[0148] · col1, col2: Specify the columns in the phenotypic data where two phenotypes are located.

[0149] ● method: Select the analysis mode, including "main" (main effect), "epistasis" (epistatic effect), or "all" (analyze both main and epistatic effects simultaneously).

[0150] ● kinship.method: Select the calculation method for the kinship matrix, including "VanRaden", "EMMA", or "Epistasis".

[0151] ● K.stand: Whether to standardize the kinship matrix.

[0152] ● epi.model: Select the epistatic effect model, including "AA", "AD", "DA", "DD".

[0153] · estimate: Optional pre-computed estimates (such as variance components).

[0154] · estimate.method: Variance component estimation method, including "default" and "errorcorrelation".

[0155] · maf.filter: Whether to filter SNPs based on MAF.

[0156] · maf.threshold: MAF filtering threshold.

[0157] · file.out: Whether to save the results to a file.

[0158] · outpfx: Prefix for the result file name.

[0159] · include.singleGWAS: Whether to include single-trait GWAS results.

[0160] · sort: Whether to sort the results by chromosome and position

[0161] · verbose: Whether to display detailed output information.

[0162] The MTMMade function returns a list containing GWAS output and variance component estimates (optional).

[0163] Example 3

[0164] See the appendix Figure 1 as shown

[0165] Based on the 1439 rice hybrid datasets published by Huang et al. (HUANG X, YANG S, GONG J, et al. Genomic analysis of hybrid rice varieties reveals numerous superior alleles that contribute to heterosis[J / OL]. Nature Communications, 2015, 6(1): 6258. DOI: 10.1038 / ncomms7258)(http: / / www.ncgr.ac.cn / RiceHap4), simulation experiments were conducted. A total of four simulated datasets were designed to simulate various scenarios of complex genetic structures with additive-dominant, epistatic effects, single QTN, and multiple QTNs.

[0166] Simulated datasets 1 - 4:

[0167] Simulated dataset 1: Based on the real 1439 samples and 934,911 SNP data (MAF > 0.05), 300 samples and 10,000 SNPs were randomly selected as the analysis basis. Four quantitative trait nucleotides (QTNs) were set in this dataset, of which 2 were additive QTNs and 2 were dominant QTNs. These 4 QTNs were alternately assigned common effects (Common) and interaction effects (Interaction), that is, 2 Common QTNs and 2 Interaction QTNs. The phenotypic variance explained rate (PVE) of each QTN was set to 5%. At the same time, 100 SNPs were selected from the remaining SNPs as background QTNs with minor effects, and their effect values followed a normal distribution with a mean of 0 and a variance of 0.005. The overall heritability of this dataset was set to 0.85. This dataset was designed to evaluate the detection ability of the MTMMade model for the common and interaction effects of additive and dominant QTNs in the presence of background effects.

[0168] Simulated dataset 2: It has the same basic data, samples, and SNP selection as dataset 1, but the number of additive QTNs and dominant QTNs is increased. A total of 12 QTNs are set, with half being additive and half being dominant. The PVE of each QTN varies between 5% and 7%, and the number and effect value distribution of background QTNs are the same as those in dataset 1. The overall heritability of this dataset is set to 0.85. This dataset is designed to evaluate the detection ability of the MTMMade model when the number of additive and dominant QTNs increases and the PVE changes.

[0169] Simulated dataset 3: Based on the same 1439 sample data as dataset 1, 300 samples are randomly selected, and 20,000 SNPs are randomly selected as the analysis basis. 20 QTNs are set in this dataset, with half being additive and half being dominant, and the effect values of the QTNs follow a geometric distribution. The overall heritability is set to 0.9. Instead of setting specific PVE values, the QTN effects are controlled by heritability and geometric distribution. This dataset is designed to evaluate the detection ability of the MTMMade model when the effects of multiple additive and dominant QTNs follow a geometric distribution.

[0170] Simulated dataset 4: Based on the same 1439 sample data as dataset 1, 400 samples are randomly selected, and 200 SNPs are randomly selected as the analysis basis. Epistatic interaction QTNs (QQIs) are designed in this dataset, only considering the additive × additive (AA) interaction type. 10 QQIs are set, and the PVE of each QQI is set to 5%. At the same time, 100 background QTNs with small effects are also set, and their effect values follow a normal distribution with a mean of 0 and a variance of 0.005. The overall heritability of this dataset is set to 0.85, and the simulation is repeated 20 times. This dataset is designed to evaluate the detection ability of the MTMMade model for epistatic interaction loci when there are epistatic interaction effects.

[0171] In the simulation experiment, this application uses classical GWAS model evaluation metrics to evaluate the performance of the MTMMade method, mainly including Power vs FDR, Power vs Type I error, and TPR. These evaluation metrics can help this application comprehensively understand the detection ability and stability of the MTMMade model under different genetic structures and effect types, so as to provide a reference for practical applications. The meanings of the metrics are as follows:

[0172] Power (statistical power): It is defined as the ability to correctly identify QTN loci related to phenotypes at a given significance level. The higher the Power value, the stronger the ability of the model to detect real effects.

[0173] FDR (False Discovery Rate): It is defined as the proportion of sites that are wrongly judged as significant among all the sites detected as significant. The calculation formula for FDR is: FDR = number of false positive sites / (number of false positive sites + number of true positive sites). The lower the FDR value, the stronger the ability of the model to control false positives.

[0174] Type I error: It is defined as the probability of wrongly detecting a significant association at non-QTN sites. In this application, the empirical null P-value distribution of all non-QTN sites is used to derive the Type I error, so as to evaluate the false positive control ability of the model.

[0175] TPR (True Positive Rate): It is defined as the proportion of QTN sites correctly detected by the model among all the set QTNs. The calculation formula for TPR is: TPR = number of correctly detected QTN sites / total number of set QTN sites. The higher the TPR value, the stronger the ability of the model to detect true QTNs.

[0176] For the analysis of the above four simulated datasets, this application used the R package MTMMade to perform multi-trait multi-effect genome-wide association analysis. The analysis results are as Figure 1 shown, clearly demonstrating the performance of MTMMade under different simulation scenarios. In the Power vs FDR evaluation ([[]] Figure 1 A-C), MTMMade showed significantly better performance than the traditional MTMM model on all simulated datasets, which means that MTMMade can better control the false discovery rate while maintaining high detection ability. In the Power vs Type I error evaluation ([[]] Figure 1 E-G), the advantage of MTMMade was also significant. At different detection thresholds, the Type I error of MTMMade was significantly lower than that of MTMM, further proving the reliability of MTMMade in controlling false positives. In addition, Figure 1 D shows the true positive rate (TPR) results of the MTMMade model under different effect types. The results show that for epistatic effects, the average TPR value of MTMMade is above 0.7, which indicates that MTMMade has a high ability to detect true QTN sites, especially in the detection of epistatic interaction effects, with a very significant advantage, because MTMM cannot detect epistatic sites. These results comprehensively show that MTMMade not only has higher detection ability under various simulation scenarios, but also can effectively control false positives, thus ensuring the reliability of the results, fully verifying the superiority of this method in the genetic analysis of complex traits.

[0177] Example 4

[0178] See the appendix Figures 2 - 3 As shown, to further verify the application effect of the MTMMade method in real data, this study used 575 rice hybrid materials and analyzed two important agronomic traits: yield per plant (YD) and filled grains per panicle (FGPP). The original genotype data contained 1,894,012 SNP loci. To control the influence of linkage disequilibrium (LD) on the analysis, this application used PLINK software to perform LD filtering on the genotype data, generating two datasets with 219,817 SNPs and 2,372 SNPs respectively. The LD filtering parameters were: plink--indep-pairwise 10000 1 0.99 and plink--indep-pairwise100000 1 0.3.

[0179] This application performed multi-trait genome-wide association analysis on the genotype data and phenotype data after LD filtering using MTMMade and the traditional MTMM method respectively. Figure 2 The full model (Full) analysis results of MTMMade and the MTMM method are shown. It can be clearly seen that the MTMM method failed to detect the significant dominant effect loci located on chromosome 11 and also could not identify many significant epistatic interaction loci. These results indicate that MTMM has obvious limitations in analyzing complex genetic effects. The MTMMade method can detect more significant loci in the full model results, including the dominant effect and epistatic effect loci that the MTMM method failed to detect.

[0180] Figure 3 The common effect (Common) analysis results of MTMMade and the MTMM method are shown. The results show that the MTMM method hardly detected any significant loci, while the MTMMade method can effectively identify multiple new significant association loci, fully demonstrating the superiority of the MTMMade method in multi-trait association analysis and complex effect analysis.

[0181] These results clearly show that the MTMMade method has stronger detection ability and higher resolution in real data analysis and can more comprehensively reveal the genetic mechanism of complex traits.

[0182] References:

[0183] [1] VARONA L, LEGARRA A, TORO M A, et al. Non-additive effects in genomic selection[J]. Frontiers in genetics, 2018, 9: 341678.

[0184] [2]NAGAI R, KINUKAWA M, WATANABE T, et al. Genome-wide detection of non-additive quantitative trait loci for semen production traits in beef and dairy bulls[J]. Animal, 2022, 16(3): 100472.

[0185] [3]REYNOLDS E G, LOPDELL T, WANG Y, et al. Non-additive QTL mapping of lactation traits in 124,000 cattle reveals novel recessive loci[J]. Genetics Selection Evolution, 2022, 54(1): 5.

[0186] [4]XUE Y, LIU S, LI W, et al. Genome-wide association study reveals additive and non-additive effects on growth traits in Duroc pigs[J]. Genes, 2022, 13(8): 1454.

[0187] [5]TSOURIS A, BRACH G, SCHACHERER J, et al. Non-additive genetic components contribute significantly to population-wide gene expression variation[J]. Cell Genomics, 2024, 4(1).

[0188] [6]LI L,ZHENG X,WANG J,et al.Joint analysis of phenotype-effect-generation identifies loci associated with grain quality traits in ricehybrids[J / OL].Nature Communications,2023,14(1):3930.DOI:10.1038 / s41467-023-39534-x.

[0189] [7]KORTE A, LMSSON B J,SEGURAV,et al.Amixed-model approach forgenome-wide association studies of correlated traits in structuredpopulations[J / OL].Nature Genetics,2012,44(9):1066-1071.DOI:10.1038 / ng.2376.

[0190] [8]ZHOU X,STEPHENS M.Efficient multivariate linear mixed modelalgorithms for genome-wide association studies[J / OL].Nature Methods,2014,11(4):407-409.DOI:10.1038 / nmeth.2848.

[0191] In the description of the present application, it should be understood that the orientation or positional relationship indicated by the terms "center", "upper", "lower", "front", "rear", "left", "right", "vertical", "horizontal", "top", "bottom", "inner", "outer", etc. is based on the orientation or positional relationship shown in the drawings. It is only for the convenience of describing the present application and simplifying the description, rather than indicating or implying that the device or element referred to must have a specific orientation, be constructed and operated in a specific orientation, and therefore should not be construed as a limitation to the present application.

[0192] The terms "first" and "second" are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the quantity of the indicated technical features. Thus, features defined with "first" and "second" may explicitly or implicitly include one or more of such features. In the description of this application, unless otherwise specified, the meaning of "a plurality of" is two or more.

[0193] In the description of this application, it should be noted that unless otherwise clearly specified and defined, the terms "installed", "connected", and "coupled" should be understood in a broad sense. For example, it may be a fixed connection, a detachable connection, or an integral connection; it may be a mechanical connection or an electrical connection; it may be directly connected or indirectly connected through an intermediate medium, and it may be the communication inside two elements. For those of ordinary skill in the art, the specific meanings of the above terms in this application can be understood according to specific circumstances.

[0194] The various embodiments in this specification are described in a progressive manner. Each embodiment focuses on the differences from other embodiments. The same or similar parts among the various embodiments can be referred to each other. For the devices disclosed in the embodiments, since they correspond to the methods disclosed in the embodiments, the description is relatively simple. For the relevant parts, reference can be made to the description in the method section.

[0195] The above description of the disclosed embodiments enables those skilled in the art to implement or use the present invention. Various modifications to these embodiments will be obvious to those skilled in the art, and the general principles defined herein can be implemented in other embodiments without departing from the spirit or scope of the present invention. Therefore, the present invention will not be limited to the embodiments shown herein, but rather to the broadest scope consistent with the principles and novel features disclosed herein.

Claims

1. A mixed linear model for genome-wide association of multiple traits, characterized in that: Based on the MTMM model, a comprehensive consideration of dominant effects and epistatic effects was added, and a false positive filtering mechanism was combined. Through the verification of simulated data and real rice data, the MTMMade model significantly improved the accuracy and statistical efficiency of detection. Specifically, the following steps are included: Step 1: Data quality control, which includes screening, checking, filling in and quality control of phenotypic and genotypic data; Step 2: Calculate the genetic variance components associated with the phenotype and integrate and process the results by using a multi-trait mixed model; Step 3: construct a multi-effect genotype matrix; Step 4: Multi-effect multi-trait genome-wide association.

2. A mixed linear model for genome-wide association of multiple traits according to claim 1, characterized in that: The step 1 includes the following steps: Step 1.1) Sample name consistency filtering, that is, comparing the sample names in the genotype and phenotype data to ensure a complete match. Any samples with missing values ​​in the genotype or phenotype data will be deleted to generate a matched and non-missing genotype-phenotype dataset; Step 1.2) Genotype data check and imputation: Check whether the genotype data is in PLINK additive coding format, where 0 is homozygous reference, 1 is heterozygous, and 2 is homozygous variant; Missing genotype data will be filled with heterozygous genotypes, coded as 1; Step 1.3) Minimum allele frequency quality control: traverse the variant sites in the genotype data, calculate the MAF, and eliminate variant sites with a MAF lower than 0.

05.

3. A mixed linear model for genome-wide association of multiple traits according to claim 1, characterized in that: The step 2 includes the following steps: Step 2.1) Input data includes sample ID, two phenotypic values ​​of all samples, denoted as Y1 and Y2, and the kinship matrix K between samples; first, Y1 and Y2 are merged into a mixed phenotype vector Y_ok, and an environmental factor vector Env is constructed to distinguish the data sources of Y1 and Y2; then, the kinship matrix K is standardized to obtain a standardized kinship matrix K_stand; Step 2.2) Preliminary evaluation of heritability: Use linear mixed models to perform preliminary heritability evaluation on phenotypes Y1 and Y2 respectively; the model forms of Y1 and Y2 are as follows: Y1=μ1+g1+e1 Y2=μ2+g2+e2 Among them, μ1 and μ2 are the population means of Y1 and Y2, g1 and g2 are the genetic effects of Y1 and Y2, e1 and e2 are the residuals of Y1 and Y2; the genetic effects are assumed to obey the normal distribution of g1~N(0,σ_g1^2*K_stand) and g2~N(0,σ_g2^2*K_stand); through this model, the heritability of Y1 and Y2 is estimated, which are recorded as herit1 and herit2 respectively; Step 2.3) Multi-trait mixed model: Use a multi-trait mixed model to analyze the genetic and environmental effects of Y1 and Y2 simultaneously; the model form is as follows: Y_ok=Xb+Zu+e Among them, Y_ok is the mixed phenotype vector, X is the design matrix containing the environmental factor Env, b is the fixed effect, Z is the design matrix, and u is the random effect, including the genetic effect, which obeys the normal distribution and has a variance-covariance structure of represents the Kronecker product, varcov is the genetic variance-covariance matrix, e is the residual, and its variance-covariance structure is ve is the residual variance-covariance matrix; through this model, the genetic variance-covariance matrix varcov and the residual variance-covariance matrix ve are estimated; Step 2.4) Correlation calculation: Based on the varcov and ve matrices, calculate the genetic correlation coefficient rho_g, environmental correlation coefficient rho_e, phenotypic correlation coefficient rho_p and Pearson correlation coefficient pearson between phenotypes Y1 and Y2; at the same time, calculate the genetic covariance G, the specific genetic variances G1 and G2 of the single trait, the environmental variance E and the environmental covariance EE; and calculate the proportion of each variance component, including the genetic variance proportion var_G, the specific genetic variance proportion var_GE, the environmental variance proportion var_E and the environmental covariance proportion var_EE; Step 2.5) Significance evaluation: The standard errors of the genetic correlation coefficient rho_g and the environmental correlation coefficient rho_e were calculated using the Delta method, denoted as se_g and se_e, respectively; at the same time, the corresponding p values ​​were calculated by the likelihood ratio test, denoted as p_g and p_e, respectively; Step 2.6) Calculation of mixed model coefficient matrix and related variables: For subsequent analysis, this step will calculate the mixed model coefficient matrix M, the adjusted phenotypic data Y_t and the covariate matrix cof_t. First, the genetic variance-covariance matrix varcov is Kronecker product with the kinship matrix K_stand to obtain the total genetic covariance matrix K_comb. Then, the residual variance-covariance matrix ve is Kronecker product with the identity matrix to obtain the residual covariance matrix I_comb. K_comb and I_comb are added to obtain the total covariance matrix bigK=K_comb+I_comb. Next, Cholesky decomposition is performed on bigK, and the inverse of the decomposed matrix is ​​calculated to obtain the mixed model coefficient matrix M=solve(chol(bigK)). On this basis, the adjusted phenotypic data Y_t=crossprod(M,Y_ok) is obtained by matrix multiplication of the mixed phenotypic vector Y_ok with the transposed matrix of M. At the same time, in order to control the fixed effects in the phenotypic data, we calculated the covariate matrix cof_t. The specific method is to merge the intercept vector Xo and the environmental factor vector Env in step 2 by column, and then multiply them with the mixed model coefficient matrix M, that is, cof_t = crossprod (M, cbind (Xo, Env)); Step 2.7) Data integration and result output: This step first integrates the calculated genetic and environmental correlations rho_g and rho_e, heritability herit1 and herit2, variance component proportions var_G, var_GE, var_E and var_EE into a list called correlation. Then, the correlation list, as well as the original phenotypic data Y, sample size n, covariate matrix cof_t, mixed model coefficient matrix M, environmental factor vector Env, genetic variance-covariance matrix varcov, adjusted phenotypic data Y_t, sample ID variable ecot_id, intercept vector Xo, and standardized kinship matrix K_stand are integrated into the estimate list. The estimate list is the final output of step 2.

4. A mixed linear model for genome-wide association of multiple traits according to claim 1, characterized in that: The step 3 includes the following steps: Step 3.1) Construct the main effect gene matrix: The quality-controlled genotype data are converted into a main effect matrix; the genotype data are coded in an additive manner, where 0 represents the homozygous reference genotype, 1 represents the heterozygous genotype, and 2 represents the homozygous variant genotype; at the same time, the genotype data are coded for dominant effects, with both the homozygous reference and the homozygous variant coded as 0, and the heterozygous genotype coded as 1; the additive coding matrix and the dominant coding matrix are merged to form a main effect gene matrix containing additive effects and dominant effects, which is used for the subsequent identification of variant sites; Step 3.2) Construct epistatic effect gene matrix: construct epistatic effect gene matrix through non-repeated permutation and combination; epistatic effect types include additive × additive, additive × dominant, dominant × additive and dominant × dominant; Step 3.3) Quality control of epistatic gene matrix: Since the dimension of epistatic effect gene matrix increases exponentially with the permutation combination, the epistatic effect gene matrix is ​​filtered for the minimum allele frequency; the specific operation is the same as step 1.3, traversing the variant sites in the epistatic gene matrix, calculating their minimum allele frequency, and filtering out variant sites with MAF less than 0.

05.

5. A mixed linear model for genome-wide association of multiple traits according to claim 1, characterized in that: The step 4 includes the following steps: Step 4.1) Input data preparation: Take the multi-effect genotype matrix constructed in step C as input, and input the adjusted phenotypic data Y_t, covariate matrix cof_t, mixed model coefficient matrix M, environmental factor vector Env, genetic variance-covariance matrix varcov, and sample ID obtained in step 2; in addition, SNP_INFO containing variant site information is also input; Step 4.2) Effect scaling: Scale the multi-effect genotype matrix; multiply the parts of the multi-effect genotype matrix related to phenotypes Y1 and Y2 by the corresponding genetic standard deviations in the varcov matrix, respectively, and record them as X_ok1 and X_ok2; merge X_ok1 and X_ok2 by row to obtain X_; Step 4.3) Construct the converted genotype matrix: Use the mixed model coefficient matrix M to convert the scaled genotype matrix X_ to the X_t array; in R language, X_t can be regarded as a set of multiple layers of two-dimensional matrices, each representing a different effect. The index [,,1] represents the first layer matrix of the X_t array, and the index [,,2] represents the second layer matrix, where X_t[,,1] is the matrix related to the total effect, and the calculation formula is X_t[,,1] = M%*%X_, and X_t[,,2] is the matrix related to the environmental interaction effect, and the calculation formula is X_t[,,2] = M%*%(Env*X_); Step 4.4) Calculate the common effect: Use a linear model to analyze the common effect of each variant site; the model form is as follows: Y_t=cof_t*β_cof+X_t[,,1]*β_com+ε Among them, Y_t is the adjusted phenotypic data, cof_t is the covariate matrix, X_t[,,1] is the part of the converted genotype matrix related to the total effect, β_cof and β_com are the effect values ​​of the covariate effect and the common effect, respectively, and ε is the residual; through linear regression, the p-value pval_com, effect value beta_com and standard error se_com of the common effect of each variant site are obtained; Step 4.5) Calculate the interaction effect: Use a linear model to analyze the environmental interaction effect of each variant site; the model form is as follows: Y_t=cof_t*β_cof+X_t[,,1]*β_com+X_t[,,2]*β_inter+ε Among them, X_t[,,2] is the part of the transformed genotype matrix related to the environmental interaction effect, and β_inter is the effect value of the interaction effect; through linear regression, the p-value pval_inter, effect value beta_inter and standard error se_inter of the interaction effect of each variant site are obtained; Step 4.6) Calculate the full model effect: Use the likelihood ratio test method to compare the full model containing all effects and the model containing only covariate effects to evaluate the overall effect of each variant site on the phenotype; use the least squares method to fit the model containing covariates and the model containing covariates and genotypes, calculate the residual sum of squares, and use the F test to calculate the p value; the calculation formula is as follows: RSS_full=sum(lsfit(cbind(cof_t,x),Y_t,intercept= FALSE)$residuals^2) is the residual sum of squares for the model including covariates and genotype effects, Where x represents the current mutation site; RSS_env=sum(lsfit(cof_t,Y_t,intercept= FALSE)$residuals^2) is the residual sum of squares of the model containing only covariates; par_env = ncol(cof_t) is the degrees of freedom for the model containing only covariates par_full = par_env + 2 is the degrees of freedom for the model including covariates and genotypes F_full=((RSS_env / RSS_full-1)*(2*n-par_full) / (par_full-par_env)) calculates the F statistic, where n represents the sample size pval_full= pf(F_full,par_full-par_env,2*n-par_full,lower.tail= FALSE) to calculate the p-value.

6. A mixed linear model for genome-wide association of multiple traits according to claim 1, characterized in that: The core mathematical expression of the MTMMade model is as follows: where y i is the observed value of the ith trait; s i is a 01 vector, which is 1 for all the values ​​belonging to the i-th trait and 0 otherwise; μ i is the phenotypic mean of the i-th trait; A j represents the jth additive effect SNP, D j represents the jth dominant effect SNP, E j represents the jth epistatic interaction SNP-SNP pair; here E j It can be further subdivided into four types of interactions: aa (additive x additive), ad (additive x dominant), da (dominant x additive), and dd (dominant x dominant); β represents the common effect of genotypes, which measures the average effect of SNP on all phenotypes, where β A , β D and β E Represent the common effects of additive, dominant and epistatic loci respectively; α represents the interaction effect between genotype and phenotype, which measures the difference in the influence of SNP on different phenotypes, among which (Aj×si)α, (Dj×si)α, (Ej×si)α represent the interaction terms between genotype and phenotype, which are only reflected in the terms related to the corresponding phenotype i; v represents the residual in the model, which contains the variation that cannot be explained by the model, including random errors and genetic random effects; when only considering the interaction, the simplified form of the MTMMade model is as follows: aa, ad, da and dd represent the four epistatic effect sites of additive × additive, additive × dominant, dominant × additive and dominant × dominant, respectively; β aa ,β ad ,β da and β dd represent the common effects of different epistatic effect sites, and α aa α ad ,α da ,α dd They represent the interaction effects between different epistatic effect sites and phenotypes.

Images