Livestock and poultry dynamic character genome association analysis method and system based on stochastic regression model

Abstract

The embodiment of the invention discloses a livestock and poultry dynamic character genome association analysis method and system based on a stochastic regression model, and relates to the technical field of data processing. The method comprises the following steps: constructing a target stochastic regression model to fit an individual phenotype trajectory, and estimating a regression coefficient of an individual phenotype curve; taking the regression coefficient as a phenotypic value, and positioning a major gene of a dynamic character growth trajectory by adopting a multivariable mixed model association analysis method; wherein the target stochastic regression model comprises a first processing layer and a second processing layer, the first processing layer is used for optimizing an individual phenotype curve according to a Bayesian information criterion so as to reduce the number of dimensions of repeated measurement, and the second processing layer is used for eliminating genetic correlation among individuals by using spectral decomposition so as to reduce the number of dimensions of repeated measurement. And the multivariate hybrid model is converted into a plurality of independent univariate hybrid models through typical transformation so as to optimize the solving method of the multi-gene variance component.

Description

Technical Field

[0001] The present application relates to the field of data processing technology, and specifically to a method and system for livestock and poultry dynamic trait genome association analysis based on a random regression model. Background Art

[0002] Dynamic traits of livestock and poultry refer to traits that change with time (life period) or other quantitative factors (physiological state, production level, metabolic rate and environmental conditions, etc.), also known as "longitudinal traits" or "repeated records". Its characteristic is that the phenotypic record of the trait can reflect the dynamic changes in the growth and development process of the organism. However, due to the large number of measurement time points, the genetic correlation of the phenotypic records of the traits at each measurement time point, and the imbalance of data, especially in the context of high-throughput sequencing technology, a large number of high-throughput genetic markers are used, making the genomic association analysis of dynamic traits extremely complex, such as occupying a large amount of computer memory and exponentially increasing calculation time.

[0003] Based on this, it is necessary to study a more scientific method for genomic association analysis of dynamic traits of livestock and poultry to improve the detection accuracy and efficiency of major effect genes. Summary of the invention

[0004] One aspect of the embodiments of this specification provides a method for livestock and poultry dynamic trait genome association analysis based on a random regression model, the method comprising:

[0005] Construct a target random regression model to fit the individual phenotypic trajectory and estimate the regression coefficient of the individual phenotypic curve;

[0006] The regression coefficient is used as the phenotypic value, and the multivariate mixed model association analysis method is used to locate the main effect gene of the dynamic trait growth trajectory; wherein, the target random regression model includes a first processing layer and a second processing layer, the first processing layer is used to optimize the individual phenotypic curve according to the Bayesian information criterion to reduce the dimensionality of repeated measurements, and the second processing layer is used to use spectral decomposition to eliminate genetic correlation between individuals, and convert the multivariate mixed model into multiple independent univariate mixed models through typical transformation to optimize the solution method of multiple gene variance components.

[0007] In some embodiments, constructing a target random regression model to fit an individual phenotypic trajectory comprises:

[0008] Using the same linear or nonlinear biological model to fit multiple factors that affect the phenotypic value of a dynamic trait, the multiple factors include fixed factor regression effects, polygenic regression effects, and permanent environmental regression effects;

[0009] Determining a first model based on the dynamic changes of various factors affecting the dynamic trait phenotypic values ​​and the target markers;

[0010] The first model is decomposed into a first processing layer and a second processing layer to obtain the target random regression model, wherein the first processing layer is a random regression model for individual phenotypic trajectories, and the second processing layer is a multivariate mixed model for phenotypic regression.

[0011] In some embodiments, the first model is represented as follows:

[0012] y=f(xF+za+G+Pe,t)+e

[0013] Where y is the column vector of repeated observation phenotypic values ​​corrected for time-independent factors; F, a, G, and Pe represent the regression effects of fixed factors, target markers, polygenes, and permanent environment, respectively; x and z represent the indicator variables of fixed factor regression effects and target marker regression effects, respectively; e represents the independent residual error between individuals;

[0014] The first processing layer and the second processing layer are respectively:

[0015]

[0016] Among them, ① represents the first processing layer, ② represents the second processing layer, and P represents the phenotypic value obtained by fitting multiple factors.

[0017] In some embodiments, the random regression model of the first processing layer estimates the phenotypic regression covariance matrix by using the restricted maximum likelihood method, estimates the individual phenotypic regression coefficients based on the phenotypic regression covariance matrix, and then selects the optimal linear or nonlinear sub-model according to the Bayesian information criterion to improve the accuracy of the analysis;

[0018] The processing process of the second processing layer includes:

[0019] Firstly, the influence of fixed factors on phenotypic regression coefficients is corrected, so that the multivariate mixed model is converted into a second model that removes the regression effect of fixed factors;

[0020] The GEMMA method is used to estimate the genome regression covariance matrix and the permanent environment regression covariance matrix of the second model;

[0021] The genome regression covariance matrix and the permanent environment regression covariance matrix are further processed by spectral decomposition to obtain corresponding eigenvalue diagonal matrices and eigenvector matrices.

[0022] In some embodiments, the processing of the second processing layer further includes:

[0023] Each variable in the multivariate mixed model is transformed according to the characteristic root diagonal matrix, the characteristic vector matrix, the genome regression covariance matrix, and the permanent environment regression covariance matrix to obtain a first transformed variable corresponding to each variable;

[0024] Substituting the corresponding variables in the second model based on the first conversion variables to obtain a third model;

[0025] The third model is decomposed to convert the multivariate mixture model into multiple independent univariate mixture models.

[0026] In some embodiments, decomposing the third model to convert the multivariate mixture model into a plurality of independent univariate mixture models comprises:

[0027] Based on the genomic kinship correlation matrix estimated using high-throughput markers and the target eigenvalues ​​in the eigenvalue diagonal matrix, the corresponding regression effect in each univariate mixed model is determined.

[0028] In some embodiments, the method further comprises:

[0029] Performing spectral decomposition on the genome kinship correlation matrix to obtain a matrix of eigenvalues ​​and eigenvectors of the genome kinship correlation matrix;

[0030] Transforming each variable in the univariate mixed model based on the eigenvalues ​​and eigenvector matrix of the genomic kinship correlation matrix to obtain a second transformed variable corresponding to each variable;

[0031] According to the second transformation variable, each univariate mixed model is transformed into a linear regression model consisting of a target label and a residual error;

[0032] The method of using multivariate mixed model association analysis to locate the main effect gene of the dynamic trait growth trajectory includes:

[0033] Inferring marker genetic effects based on the least squares method and a linear regression model corresponding to each of the univariate mixed models;

[0034] Add the test statistics corresponding to the d marker genetic effects to obtain the one-factor, multiple-effect QTN analysis result, where the statistic follows a chi-square distribution with d degrees of freedom;

[0035] Based on the analysis results of the one-cause-multiple-effect QTN, the major effect gene is located from multiple candidate QTNs.

[0036] In some embodiments, locating the major effect gene from multiple candidate QTNs based on the one-cause pleiotropic QTN analysis results includes: further improving the detection accuracy and efficiency of the major effect gene by considering the linkage disequilibrium between the candidate QTNs.

[0037] Another aspect of the embodiments of this specification also provides a system for livestock and poultry dynamic trait genome association analysis based on a random regression model, the system comprising:

[0038] Model building module, used to build the target random regression model;

[0039] The genomic association analysis module is used to fit the individual phenotypic trajectory based on the target random regression model and estimate the regression coefficient of the individual phenotypic curve, and then use the regression coefficient as the phenotypic value to locate the main effect gene of the dynamic trait growth trajectory by adopting the multivariate mixed model association analysis method; wherein, the target random regression model includes a first processing layer and a second processing layer, the first processing layer is used to optimize the individual phenotypic curve according to the Bayesian information criterion to reduce the dimensionality of repeated measurements, and the second processing layer is used to eliminate genetic correlation between individuals by using spectral decomposition, and convert the multivariate mixed model into multiple independent univariate mixed models by typical transformation to optimize the solution method of multiple gene variance components.

[0040] Another aspect of the embodiments of this specification also provides a computer-readable storage medium, which stores computer instructions. When the computer instructions are executed by a processor, the method for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model as described in any one of the above is implemented.

[0041] The beneficial effects that the method and system for livestock and poultry dynamic trait genome association analysis based on the random regression model provided in the embodiments of this specification may bring include at least:

[0042] (1) By decomposing the first model determined by the target marker and multiple factors affecting the phenotypic value of the dynamic trait into a first processing layer and a second processing layer, and then optimizing the individual phenotypic curve in the first processing layer according to the Bayesian information criterion, the dimension of repeated measurements can be reduced, and spectral decomposition is used in the second processing layer to eliminate the genetic correlation between individuals, and the multivariate mixed model is converted into multiple independent univariate mixed models through typical transformation, which can optimize the solution method of the multi-gene variance component. On the basis of marker-by-marker testing, multiple candidate main effect genes are jointly analyzed to improve the statistical detection efficiency;

[0043] (2) The Bayesian information criterion is used to select the optimal linear or nonlinear sub-model, which can balance the model complexity and goodness of fit, and avoid overfitting through penalty terms, thereby selecting a simpler model with higher explanatory power and improving the accuracy of subsequent genetic analysis;

[0044] (3) By considering linkage disequilibrium among candidate QTNs, the detection accuracy and efficiency of major effect genes can be further improved.

[0045] Additional features will be described in part in the following description. They will become apparent to those skilled in the art by reviewing the following and accompanying drawings, or may be understood by the production or operation of the examples. The features of this specification may be implemented and obtained by practicing or using various aspects of the methods, tools, and combinations described in the following detailed examples. BRIEF DESCRIPTION OF THE DRAWINGS

[0046] This specification will be further described in the form of exemplary embodiments, which will be described in detail by the accompanying drawings. These embodiments are not restrictive, and in these embodiments, the same number represents the same structure, wherein:

[0047] Figure 1 is an exemplary flow chart of a method for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model according to some embodiments of this specification;

[0048] Figure 2 It is an exemplary module diagram of a livestock and poultry dynamic trait genomic association analysis system based on a random regression model according to some embodiments of this specification. DETAILED DESCRIPTION

[0049] In order to more clearly illustrate the technical solutions of the embodiments of this specification, the following is a brief introduction to the drawings required for the description of the embodiments. Obviously, the drawings described below are only some examples or embodiments of this specification. For ordinary technicians in this field, without paying creative work, this specification can also be applied to other similar scenarios based on these drawings. Unless it is obvious from the language environment or otherwise explained, the same reference numerals in the figures represent the same structure or operation.

[0050] It should be understood that the "system", "device", "unit" and / or "module" used in this specification is a method for distinguishing different components, elements, parts, portions or assemblies at different levels. However, if other words can achieve the same purpose, the words can be replaced by other expressions.

[0051] As shown in this specification and claims, unless the context clearly indicates an exception, the words "a", "an", "an" and / or "the" do not refer to the singular and may also include the plural. Generally speaking, the terms "comprise" and "include" only indicate the inclusion of the steps and elements that have been clearly identified, and these steps and elements do not constitute an exclusive list. The method or device may also include other steps or elements.

[0052] Flowcharts are used in this specification to illustrate the operations performed by the system according to the embodiments of this specification. It should be understood that the preceding or following operations are not necessarily performed precisely in order. Instead, the steps may be processed in reverse order or simultaneously. At the same time, other operations may also be added to these processes, or one or more operations may be removed from these processes.

[0053] Dynamic traits of livestock and poultry refer to the changes in traits exhibited by livestock and poultry (such as pigs, cattle, sheep, chickens, ducks, etc.) at different life stages or under the influence of other quantitative factors, such as physiological state, production level, metabolic rate and environmental conditions. Such traits are often referred to as "longitudinal traits" or "repeated records". Its notable feature is that the phenotypic record of the trait can reflect the dynamic changes of the organism during growth and development. However, due to the large number of measured time points, there are genetic correlations between the phenotypic records of the traits at each time point, and the data distribution is uneven. Especially in the context of the widespread application of high-throughput sequencing technology, the use of a large number of high-throughput genetic markers makes the genomic association analysis of dynamic traits extremely complicated. This complexity is mainly reflected in two aspects: first, a large amount of computer memory is required during the analysis process, and second, the time required for calculation will increase exponentially, which brings huge challenges to the research work.

[0054] In view of the above problems, the present application provides a fast algorithm for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model to improve the efficiency of detecting major effect genes of important economic traits of livestock and poultry and the accuracy of genome selection. The method and system for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model provided in the embodiments of this specification are described in detail below in conjunction with the accompanying drawings.

[0055] Figure 1 This is an exemplary flow chart of the method for livestock and poultry dynamic trait genomic association analysis based on a random regression model according to some embodiments of this specification. In some embodiments, the livestock and poultry dynamic trait genomic association analysis method 100 based on a random regression model can be executed by processing logic, which can include hardware (e.g., circuits, dedicated logic, programmable logic, microcode, etc.), software (instructions running on a processing device to perform hardware simulation), etc. or any combination thereof. In some embodiments, Figure 1 One or more operations in the flowchart of the method 100 for analyzing the association between dynamic traits of livestock and poultry based on a random regression model can be performed by Figure 2 The livestock and poultry dynamic trait genome association analysis system 200 based on the random regression model is implemented. For example, the livestock and poultry dynamic trait genome association analysis method 100 based on the random regression model can be stored in a storage device in the form of instructions and used by Figure 2The model building module 210 and / or the genomic association analysis module 220 are shown calling and / or executing.

[0056] Reference Figure 1 In some embodiments, the method 100 for analyzing the association between dynamic traits of livestock and poultry based on a random regression model may include:

[0057] Step S110, constructing a target random regression model. In some embodiments, step S110 may be performed by Figure 2 The model building module 210 is shown executing.

[0058] In the embodiment of the present application, the same linear or nonlinear biological model f(b, t) can be used to fit various factors that affect the dynamic trait phenotypic value, which may include the dynamic changes of the fixed factor regression effect F, the polygenic regression effect G, and the permanent environmental regression effect Pe. Among them, the biological model refers to a theoretical framework or mathematical expression used to describe, explain and predict biological phenomena. In the embodiment of the present application, the biological model can be understood as a mathematical model (including linear and / or nonlinear models) that uses mathematical formulas and equations to describe the dynamic traits of livestock and poultry. Fixed factor regression effects refer to those variables that remain unchanged in the study, such as sex, breed or treatment group. These factors usually appear in the form of categorical variables in the model, and their regression coefficients represent the average effect of the fixed factors on the phenotypic value of the trait. The polygenic regression effect involves the effect of genetic background on the trait. In the embodiment of the present application, the polygenic effect can refer to those genetic characteristics that are jointly acted by multiple genes, but the effect of a single gene is small. The effects of these genes can be estimated by genomic selection technology and reflected in the model as continuous variables, thereby revealing how the effect of genetic background on the phenotypic value of the trait changes under different growth stages or environmental conditions. Permanent environmental regression effects refer to those environmental factors that are closely related to individuals and remain relatively stable throughout the life cycle of the individual. For example, for dairy cows, their feeding environment in a specific pasture may have a lasting impact on their milk production. In the embodiments of the present application, the dynamic changes in the permanent environmental effect can reflect how these environmental factors affect the phenotypic values ​​of the traits over time.

[0059] Further, in the embodiment of the present application, a first model can be determined based on the dynamic changes of various factors affecting the phenotypic value of the dynamic trait and the target marker, and the first model can be understood as a genomic random regression model of a dynamic trait. Among them, the target marker refers to a specific gene or DNA sequence related to a specific dynamic trait, and these dynamic traits may include growth rate, feed conversion efficiency, meat production performance, reproductive capacity, etc., and these traits will change over time.

[0060] In the embodiment of the present application, the aforementioned first model can be expressed as:

[0061] y=f(xF+za+G+Pe,t)+e (1)

[0062] Where y is the column vector of repeated observed phenotypic values ​​corrected for the time-independent factor (t); F, a, G, and Pe represent the regression effects of fixed factors, target markers, polygenes, and permanent environment, respectively; x and z represent the indicator variables of the fixed factor regression effect and the target marker regression effect, respectively; and e represents the independent residual error between individuals.

[0063] In the present embodiment, the phenotype value can be Then the first model is decomposed into a first processing layer and a second processing layer to obtain a target random regression model. The first processing layer and the second processing layer can be expressed as:

[0064]

[0065] Among them, ① represents the first processing layer, which is a random regression model for individual phenotypic trajectories; ② represents the second processing layer, which is a multivariate mixed model for phenotypic regression; P represents the phenotypic value obtained by fitting multiple factors.

[0066] In the embodiment of the present application, the first processing layer can be used to optimize the individual phenotypic curve according to the Bayesian information criterion to reduce the dimensionality of repeated measurements, and the second processing layer can be used to eliminate genetic correlation between individuals using spectral decomposition, and convert the multivariate mixed model into multiple independent univariate mixed models through typical transformation to optimize the solution method of the multi-gene variance component, and on the basis of marker-by-marker testing, jointly analyze multiple candidate main effect genes to improve statistical detection efficiency. For more information about the first processing layer and the second processing layer, please refer to the following text, which will not be discussed in detail here.

[0067] Step S120 , fitting the individual phenotype trajectory based on the target random regression model, and estimating the regression coefficient of the individual phenotype curve. In some embodiments, step S120 may be performed by the genome association analysis module 220 .

[0068] In an embodiment of the present application, after the target random regression model is constructed in the above manner, the individual phenotypic trajectory can be fitted based on the target random regression model, and the regression coefficient of the individual phenotypic curve is estimated. Specifically, the SNP test data can be obtained, and then the SNP test data is processed by the target random regression model, the individual phenotypic trajectory is fitted, and the regression coefficient of the individual phenotypic curve is estimated. After completing the processing of the SNP test data, researchers can further use these regression coefficients for in-depth analysis. For example, by comparing the regression coefficient differences between different individuals, the influence of genes on phenotypic changes can be revealed. In addition, in combination with the genotype data of the individual, the correlation between specific genes and phenotypes can be further explored, thereby providing a scientific basis for the detection of major effect genes.

[0069] In some embodiments, the random regression model of the first processing layer can use the REML (Restricted Maximum Likelihood) method to estimate the phenotypic regression covariance matrix, and estimate the individual phenotypic regression coefficients based on the phenotypic regression covariance matrix, and then select the optimal linear or nonlinear sub-model according to the BIC criterion (Bayesian Information Criterion) to improve the accuracy of the analysis.

[0070] Specifically, the REML method is a statistical method that can be used to estimate the parameters in the variance component model, and it focuses on the variance components of random effects, rather than all parameters. In this way, REML can provide unbiased fixed effect estimation, and can also estimate the variance of random effects. In an embodiment of the present application, by adopting the REML method to estimate the phenotypic regression covariance matrix, and based on the phenotypic regression covariance matrix to estimate the individual phenotypic regression coefficient, the variance components of these random effects can be effectively separated, thereby more accurately estimating the genetic effects of an individual.

[0071] In some embodiments, the process of estimating individual phenotypic regression coefficients based on the phenotypic regression covariance matrix can be implemented based on BLUP (Best Linear Unbiased Prediction).

[0072] Further, in an embodiment of the present application, the optimal linear or nonlinear submodel can be selected according to the BIC criterion to improve the accuracy of the analysis. Specifically, the BIC criterion is a statistical model selection criterion, which can be weighed between model complexity and goodness of fit, and overfitting is avoided by a penalty term, thereby selecting a more concise model with higher explanatory power. When selecting the optimal linear or nonlinear submodel, the BIC criterion can consider the number of model parameters and the likelihood of data. In an embodiment of the present application, by calculating the BIC value of different submodels, the optimal model (for example, the model with the minimum BIC value is selected as the optimal model) can be effectively screened out, thereby improving the accuracy of genetic analysis.

[0073] Step S130 , using the regression coefficient as the phenotypic value, and using a multivariate mixed model association analysis method to locate the main effect gene of the dynamic trait growth trajectory. In some embodiments, step S130 can be performed by the genome association analysis module 220 .

[0074] After obtaining the regression coefficient of the individual phenotypic curve through the above steps, the regression coefficient can be used as the phenotypic value, and then the main effect gene of the dynamic trait growth trajectory is located by using the multivariate mixed model association analysis method through the second processing layer.

[0075] Specifically, in the embodiment of the present application, the effect of fixed factors on the phenotypic regression coefficient can be corrected first, so that the multivariate mixed model is converted into a second model that removes the fixed factor regression effect. The second model can be expressed as:

[0076]

[0077] Furthermore, the genomic regression covariance matrix V of the second model can be estimated using the GEMMA method (Zhou, et al., 2012). G and the permanent environmental regression covariance matrix V Pe , and then the genome regression covariance matrix and the permanent environment regression covariance matrix are further processed by spectral decomposition to obtain the corresponding eigenroot diagonal matrix and eigenvector matrix.

[0078] It should be noted that GEMMA (Genome-wide Efficient Mixed Model Association) is an efficient software package for genome-wide association analysis (GWAS), which is mainly used to detect the correlation between genotype and phenotype to identify genetic variations that affect specific traits. GEMMA supports linear and generalized linear mixed models and can play a role in single nucleotide polymorphism (SNP) and other genetic data analysis. About using GEMMA method to estimate the genomic regression covariance matrix V of the second modelG and the permanent environmental regression covariance matrix V Pe The specific implementation process can be regarded as the prior art and will not be discussed in detail in this specification.

[0079] Spectral decomposition, also known as eigendecomposition, or similar standard form decomposition, is a method of decomposing a matrix into a product of matrices represented by its eigenvalues ​​(or eigenroots) and eigenvectors. In the embodiment of the present application, the above spectral decomposition can be expressed as follows:

[0080]

[0081] Among them, S Pe and S λ is the permanent environmental regression covariance matrix V Pe and The characteristic root diagonal matrix, U Pe and U λ is the corresponding eigenvector matrix.

[0082] In an embodiment of the present application, in order to improve the efficiency of the association analysis of multiple regression coefficients, the multivariate mixed model of the second processing layer can be converted into multiple independent univariate mixed models through typical transformations. Specifically, each variable in the multivariate mixed model can be converted according to the characteristic root diagonal matrix, the characteristic vector matrix, the genome regression covariance matrix, and the permanent environment regression covariance matrix obtained in the aforementioned process to obtain the first conversion variable corresponding to each variable. The process can be expressed as follows:

[0083]

[0084] Among them, P * represents the first conversion variable corresponding to the phenotype value P; a * represents the first conversion variable corresponding to the target label a; G * Represents the first transformation variable corresponding to the polygenic regression effect G; Pe * Represents the first transformation variable corresponding to the permanent environmental regression effect Pe.

[0085] Furthermore, the corresponding variables in the second model are replaced based on the first conversion variables to obtain a third model. The third model can be expressed as follows:

[0086]

[0087] It is understandable that since V(P * )=V(G * +Pe * )=S λ +I d(Where V(*) represents the variance component corresponding to variable *, d is the number of regression coefficients, I d represents the product of the residual variance corresponding to the dth regression coefficient (or the dth trait) and the identity matrix. The third model above is decomposed into d independent univariate mixed models:

[0088]

[0089] Among them, l = 1, 2, ..., d, which represents the lth single variable (which can be understood as the lth trait), represents the phenotypic value corresponding to the lth single variable, represents the target label corresponding to the lth single variable, represents the polygenic regression effect corresponding to the lth univariate, Represents the permanent environmental regression effect corresponding to the lth single variable.

[0090] In the embodiment of the present application, the corresponding regression effect in each univariate mixed model can be determined based on the genomic kinship correlation matrix estimated by high-throughput markers and the target characteristic root in the aforementioned characteristic root diagonal matrix. Specifically, in the process of solving d independent univariate mixed models respectively: in, Denotes the aforementioned characteristic root diagonal matrix S λ The lth characteristic root in , K is the genomic kinship matrix estimated using high-throughput markers, so (I n It represents the product of the residual variance corresponding to n traits and the identity matrix, and is used to represent the variance caused by non-genetic or environmental factors).

[0091] Furthermore, in order to eliminate the correlation between individuals, the aforementioned genomic kinship correlation matrix K can be spectrally decomposed to obtain the characteristic roots and eigenvector matrix of the genomic kinship correlation matrix, and then each variable in the univariate mixed model is transformed based on the characteristic roots and eigenvector matrix of the genomic kinship correlation matrix to obtain the second transformed variable corresponding to each variable. Specifically, the spectral decomposition of the genomic kinship correlation matrix K can be expressed as:

[0092]

[0093] Among them, S k and U k is the matrix of characteristic roots and eigenvectors of the genomic kinship correlation matrix K.

[0094] Furthermore, each variable in the univariate mixed model (4) can be transformed based on the eigenvalues ​​and eigenvector matrix of the genomic kinship matrix to obtain a second transformed variable corresponding to each variable. Specifically, in some embodiments, it can be set as follows:

[0095]

[0096] Among them, p' l express The corresponding second conversion variable, z' represents the second conversion variable corresponding to z, g' l express The corresponding second conversion variable, pe' l express The corresponding second conversion variable.

[0097] Through the above transformation, the above polygenic regression effect (g' l ) and the regression effect of the permanent environment (pe' l ) is regarded as a fixed quantity, then according to the above second conversion variable, each univariate mixed model (4) can be converted into a model consisting of the target marker (a' l ) and residual error (e' l ). The linear regression model can be expressed as:

[0098] ”'

[0099] p l =za l +e l (5)

[0100] Furthermore, the marker genetic effect can be inferred based on the least squares method and the linear regression model (5) corresponding to each of the univariate mixed models. This process can be expressed as:

[0101]

[0102] in, represents the estimated value corresponding to the lth target label, express The variance of It represents the test statistic corresponding to the genetic effect of the lth marker, which reflects the degree of association between the target marker and the trait and its reliability.

[0103] To further verify the significance of these markers, the test statistic corresponding to the genetic effect of d markers can be Add together to obtain the result of a one-factor, multi-effect QTN (Quantitative Trait Nucleotide, quantitative trait locus) analysis, where the statistic is obtained based on the genetic effect of the target marker and obeys a chi-square distribution with d degrees of freedom. Specifically, when a marker is significantly correlated with a trait, its corresponding test statistic will be larger, indicating that the marker is statistically significant. In the embodiment of the present application, the one-factor, multi-effect QTN analysis result can be expressed as: It can reflect the influence of the same quantitative trait locus on multiple traits.

[0104] Furthermore, the main effect gene can be located from multiple candidate QTNs based on the analysis results of the one-cause-multiple-effect QTN. Specifically, in some embodiments, multiple candidate QTNs can be obtained from the SNP test data one by one at a significance level lower than the Bonferroni correction threshold, and then the univariate mixed model of the multiple candidate QTNs can be jointly analyzed by the stepwise elimination method to determine the main effect gene. Exemplarily, in some embodiments, each candidate QTN can first be subjected to a univariate analysis to evaluate its contribution to the trait variation. In this way, those QTNs that have a significant effect on the trait can be identified. Subsequently, the stepwise elimination method can be used, that is, the QTN with the smallest effect on the trait is excluded in each iteration until the main effect gene that contributes the most to the trait variation is found.

[0105] In the process of determining the main effect gene, the interaction between QTNs can also be considered. For example, some QTNs may show significant effects only under a specific genetic background, a phenomenon called epistasis. Therefore, a multivariate model including interaction terms between QTNs can be constructed to more comprehensively evaluate the effect of QTN on traits.

[0106] In some embodiments of the present application, the detection accuracy and efficiency of the major effect gene can be further improved by considering the linkage disequilibrium between the candidate QTNs. For example, the linkage disequilibrium values ​​(such as D' and r) between different QTNs can be calculated. 2 ), highly linked QTN pairs can be identified. These highly linked QTN pairs often have similar genetic backgrounds, so they can be treated as a whole in the analysis, thereby reducing redundant information and improving analysis efficiency. In addition, the use of linkage disequilibrium information can also help us better understand the distribution pattern of QTN on chromosomes. For example, some QTNs may show higher linkage disequilibrium values ​​in specific chromosomal regions, which may suggest the presence of important regulatory elements or gene clusters in these regions. By further analyzing these regions, candidate genes associated with the traits can be revealed, thereby providing a more targeted basis for functional verification and genomic association analysis.

[0107] Figure 2 This is a schematic diagram of a module of a livestock and poultry dynamic trait genome association analysis system based on a random regression model according to some embodiments of this specification. In some embodiments, Figure 2 The livestock and poultry dynamic trait genomic association analysis system 200 based on the random regression model shown can be implemented in the form of software and / or hardware. For example, it can be configured in the form of software and / or hardware to a processing device and / or a terminal device for processing SNP (Single Nucleotide Polymorphism) test data and performing genomic association analysis based on the SNP test data to detect major effect genes.

[0108] Reference Figure 2 In some embodiments, the livestock and poultry dynamic trait genomic association analysis system 200 based on the random regression model may include a model construction module 210 and a genomic association analysis module 220. The model construction module 210 may be used to construct a target random regression model; the genomic association analysis module 220 may be used to fit the individual phenotypic trajectory based on the target random regression model, and estimate the regression coefficient of the individual phenotypic curve, and then use the regression coefficient as the phenotypic value to locate the main effect gene of the dynamic trait growth trajectory by using the multivariate mixed model association analysis method; wherein the target random regression model includes a first processing layer and a second processing layer, the first processing layer is used to optimize the individual phenotypic curve according to the Bayesian information criterion to reduce the dimensionality of repeated measurements, and the second processing layer is used to use spectral decomposition to eliminate genetic correlation between individuals, and convert the multivariate mixed model into multiple independent univariate mixed models through typical transformation to optimize the solution method of polygenic variance components.

[0109] For more details about the above modules, please refer to other places in this manual (for example Figure 1 part and its related description), which will not be repeated here.

[0110] It should be understood that Figure 2The livestock and poultry dynamic trait genome association analysis system 200 based on the random regression model and its modules can be implemented in various ways. For example, in some embodiments, the system and its modules can be implemented by hardware, software, or a combination of software and hardware. Among them, the hardware part can be implemented using dedicated logic; the software part can be stored in a memory and executed by an appropriate instruction execution system, such as a microprocessor or a dedicated design hardware. Those skilled in the art will understand that the above methods and systems can be implemented using computer executable instructions and / or included in a processor control code, such as a carrier medium such as a disk, CD or DVD-ROM, a programmable memory such as a read-only memory (firmware), or a data carrier such as an optical or electronic signal carrier. Such code is provided. The system and its modules of this specification can not only be implemented by hardware circuits such as ultra-large-scale integrated circuits or gate arrays, semiconductors such as logic chips, transistors, etc., or programmable hardware devices such as field programmable gate arrays, programmable logic devices, etc., but can also be implemented by software executed by various types of processors, and can also be implemented by a combination of the above hardware circuits and software (e.g., firmware).

[0111] It should be noted that the above description of the livestock and poultry dynamic trait genome association analysis system 200 based on the random regression model is provided for illustrative purposes only and is not intended to limit the scope of this specification. It is understandable that for those skilled in the art, according to the description of this specification, without departing from this principle, the various modules can be arbitrarily combined, or a subsystem can be formed to connect with other modules. For example, Figure 2 The model building module 210 and the genome association analysis module 220 described in the above can be different modules in one system, or one module can realize the functions of two or more modules mentioned above. For another example, the livestock and poultry dynamic trait genome association analysis system 200 based on the random regression model can also include other modules. Such variations are all within the protection scope of this specification.

[0112] The embodiment of the present application also provides a computer-readable storage medium, which stores computer instructions. When the computer instructions are executed by a processor, the method for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model described in any embodiment of this specification can be implemented. The specific content of the method can be referred to above and will not be repeated here.

[0113] In summary, the beneficial effects that may be brought about by the embodiments of the present specification include but are not limited to: (1) In the method and system for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model provided in some embodiments of the present specification, the first model determined based on the target marker and multiple factors affecting the phenotypic value of the dynamic trait is decomposed into a first processing layer and a second processing layer, and then the individual phenotypic curve is optimized according to the Bayesian information criterion in the first processing layer, so that the dimension of repeated measurements can be reduced, and spectral decomposition is used in the second processing layer to eliminate genetic correlation between individuals, and the multivariate mixed model is converted into multiple independent univariate mixed models through typical transformation, so that the solution method of multi-gene variance components can be optimized, and on the basis of marker-by-marker testing, the combined Analyze multiple candidate major effect genes to improve statistical detection efficiency; (2) In the livestock and poultry dynamic trait genomic association analysis method and system based on random regression model provided in some embodiments of the present specification, the optimal linear or nonlinear sub-model is selected by Bayesian information criterion, and a trade-off can be made between model complexity and goodness of fit. Overfitting can be avoided by penalty terms, thereby selecting a simpler model with higher explanatory power, thereby improving the accuracy of subsequent genetic analysis; (3) In the livestock and poultry dynamic trait genomic association analysis method and system based on random regression model provided in some embodiments of the present specification, the detection accuracy and efficiency of major effect genes can be further improved by considering the linkage disequilibrium between candidate QTNs.

[0114] It should be noted that different embodiments may produce different beneficial effects. In different embodiments, the beneficial effects that may be produced may be any one or a combination of the above, or may be any other beneficial effects that may be obtained.

[0115] The basic concepts have been described above. Obviously, for those skilled in the art, the above detailed disclosure is only for example and does not constitute a limitation of this specification. Although not explicitly stated here, those skilled in the art may make various modifications, improvements and corrections to this specification. Such modifications, improvements and corrections are suggested in this specification, so such modifications, improvements and corrections still belong to the spirit and scope of the exemplary embodiments of this specification.

[0116] At the same time, this specification uses specific words to describe the embodiments of this specification. For example, "one embodiment", "an embodiment", and / or "some embodiments" refer to a certain feature, structure or characteristic related to at least one embodiment of this specification. Therefore, it should be emphasized and noted that "one embodiment" or "an embodiment" or "an alternative embodiment" mentioned twice or more in different positions in this specification does not necessarily refer to the same embodiment. In addition, certain features, structures or characteristics in one or more embodiments of this specification can be appropriately combined.

[0117] In addition, it will be understood by those skilled in the art that various aspects of this specification may be illustrated and described by a number of patentable categories or situations, including any new and useful process, machine, product or combination of substances, or any new and useful improvements thereto. Accordingly, various aspects of this specification may be performed entirely by hardware, entirely by software (including firmware, resident software, microcode, etc.), or by a combination of hardware and software. The above hardware or software may be referred to as "data blocks", "modules", "engines", "units", "components" or "systems". In addition, various aspects of this specification may be represented as a computer product located in one or more computer-readable media, which includes computer-readable program code.

[0118] A computer storage medium may include a propagated data signal containing computer program code, for example, in baseband or as part of a carrier wave. The propagated signal may be in a variety of forms, including electromagnetic, optical, etc., or a suitable combination. A computer storage medium may be any computer-readable medium other than a computer-readable storage medium, which can be connected to an instruction execution system, device or apparatus to communicate, propagate or transmit the program for use. The program code on the computer storage medium may be transmitted via any suitable medium, including radio, cable, fiber optic cable, RF, or similar media, or any combination of the above media.

[0119] The computer program code required for the operation of each part of this specification can be written in any one or more programming languages, including object-oriented programming languages ​​such as Java, Scala, Smalltalk, Eiffel, JADE, Emerald, C++, C#, VB.NET, Python, etc., conventional procedural programming languages ​​such as C language, Visual Basic, Fortran2003, Perl, COBOL2002, PHP, ABAP, dynamic programming languages ​​such as Python, Ruby and Groovy, or other programming languages, etc. The program code can be run entirely on the user's computer, or run on the user's computer as an independent software package, or run partially on the user's computer and partially on a remote computer, or run entirely on a remote computer or processing device. In the latter case, the remote computer can be connected to the user's computer through any network form, such as a local area network (LAN) or a wide area network (WAN), or connected to an external computer (e.g., via the Internet), or in a cloud computing environment, or used as a service such as software as a service (SaaS).

[0120] In addition, unless explicitly stated in the claims, the order of the processing elements and sequences described in this specification, the use of alphanumeric characters, or the use of other names are not intended to limit the order of the processes and methods of this specification. Although the above disclosure discusses some invention embodiments that are currently considered useful through various examples, it should be understood that such details are only for illustrative purposes, and the attached claims are not limited to the disclosed embodiments. On the contrary, the claims are intended to cover all modifications and equivalent combinations that are consistent with the essence and scope of the embodiments of this specification. For example, although the system components described above can be implemented by hardware devices, they can also be implemented only by software solutions, such as installing the described system on an existing processing device or mobile device.

[0121] Similarly, it should be noted that in order to simplify the description disclosed in this specification and thus help understand one or more embodiments of the invention, in the above description of the embodiments of this specification, multiple features are sometimes combined into one embodiment, figure or description thereof. However, this disclosure method does not mean that the features required by the subject matter of this specification are more than the features mentioned in the claims. In fact, the features of the embodiments are less than all the features of the single embodiment disclosed above.

[0122] In some embodiments, numbers describing the number of components and attributes are used. It should be understood that such numbers used in the description of the embodiments are modified by the modifiers "about", "approximately" or "substantially" in some examples. Unless otherwise specified, "about", "approximately" or "substantially" indicate that the numbers are allowed to vary by ±20%. Accordingly, in some embodiments, the numerical parameters used in the specification and claims are approximate values, which may change according to the required features of individual embodiments. In some embodiments, the numerical parameters should take into account the specified significant digits and adopt the general method of retaining digits. Although the numerical domains and parameters used to confirm the breadth of their range in some embodiments of this specification are approximate values, in specific embodiments, the setting of such numerical values ​​is as accurate as possible within the feasible range.

[0123] Each patent, patent application, patent application publication, and other materials, such as articles, books, specifications, publications, documents, etc., cited in this specification are hereby incorporated by reference in their entirety. Except for application history documents that are inconsistent with or conflicting with the contents of this specification, documents that limit the broadest scope of the claims of this specification (currently or later attached to this specification) are also excluded. It should be noted that if the descriptions, definitions, and / or use of terms in the materials attached to this specification are inconsistent or conflicting with the contents described in this specification, the descriptions, definitions, and / or use of terms in this specification shall prevail.

[0124] Finally, it should be understood that the embodiments described in this specification are only used to illustrate the principles of the embodiments of this specification. Other variations may also fall within the scope of this specification. Therefore, as an example and not a limitation, alternative configurations of the embodiments of this specification may be considered consistent with the teachings of this specification. Accordingly, the embodiments of this specification are not limited to the embodiments explicitly introduced and described in this specification.

Claims

1. A method for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model, characterized in that: include: Construct a target random regression model to fit the individual phenotypic trajectory and estimate the regression coefficient of the individual phenotypic curve; The regression coefficient is used as the phenotypic value, and the multivariate mixed model association analysis method is used to locate the main effect gene of the dynamic trait growth trajectory; wherein, the target random regression model includes a first processing layer and a second processing layer, the first processing layer is used to optimize the individual phenotypic curve according to the Bayesian information criterion to reduce the dimensionality of repeated measurements, and the second processing layer is used to use spectral decomposition to eliminate genetic correlation between individuals, and convert the multivariate mixed model into multiple independent univariate mixed models through typical transformation to optimize the solution method of multiple gene variance components.

2. The method according to claim 1, characterized in that The method of constructing a target random regression model to fit individual phenotypic trajectories includes: Using the same linear or nonlinear biological model to fit multiple factors that affect the phenotypic value of a dynamic trait, the multiple factors include fixed factor regression effects, polygenic regression effects, and permanent environmental regression effects; Determining a first model based on the dynamic changes of various factors affecting the dynamic trait phenotypic values ​​and the target markers; The first model is decomposed into a first processing layer and a second processing layer to obtain the target random regression model, wherein the first processing layer is a random regression model for individual phenotypic trajectories, and the second processing layer is a multivariate mixed model for phenotypic regression.

3. The method according to claim 2, characterized in that The first model is expressed as follows: y=f(xF+za+G+Pe,t)+e Where y is the column vector of repeated observation phenotypic values ​​corrected for time-independent factors; F, a, G, and Pe represent the regression effects of fixed factors, target markers, polygenes, and permanent environment, respectively; x and z represent the indicator variables of fixed factor regression effects and target marker regression effects, respectively; e represents the independent residual error between individuals; The first processing layer and the second processing layer are respectively: Among them, ① represents the first processing layer, ② represents the second processing layer, and P represents the phenotypic value obtained by fitting multiple factors.

4. The method according to claim 2, characterized in that The random regression model of the first processing layer estimates the phenotypic regression covariance matrix by using the restricted maximum likelihood method, estimates the individual phenotypic regression coefficients based on the phenotypic regression covariance matrix, and then selects the optimal linear or nonlinear sub-model according to the Bayesian information criterion to improve the accuracy of the analysis; The processing process of the second processing layer includes: Firstly, the influence of fixed factors on phenotypic regression coefficients is corrected, so that the multivariate mixed model is converted into a second model that removes the regression effect of fixed factors; The GEMMA method is used to estimate the genome regression covariance matrix and the permanent environment regression covariance matrix of the second model; The genome regression covariance matrix and the permanent environment regression covariance matrix are further processed by spectral decomposition to obtain corresponding eigenvalue diagonal matrices and eigenvector matrices.

5. The method according to claim 4, characterized in that The processing process of the second processing layer also includes: Transform each variable in the multivariate mixed model according to the characteristic root diagonal matrix, the characteristic vector matrix, the genome regression covariance matrix, and the permanent environment regression covariance matrix to obtain a first transformed variable corresponding to each variable; Substituting the corresponding variables in the second model based on the first conversion variables to obtain a third model; The third model is decomposed to convert the multivariate mixture model into multiple independent univariate mixture models.

6. The method according to claim 5, characterized in that Decomposing the third model to convert the multivariate mixed model into a plurality of independent univariate mixed models includes: Based on the genomic kinship correlation matrix estimated using high-throughput markers and the target eigenvalues ​​in the eigenvalue diagonal matrix, the corresponding regression effect in each univariate mixed model is determined.

7. The method according to claim 6, characterized in that The method further comprises: Performing spectral decomposition on the genome kinship correlation matrix to obtain a matrix of eigenvalues ​​and eigenvectors of the genome kinship correlation matrix; Transforming each variable in the univariate mixed model based on the eigenvalues ​​and eigenvector matrix of the genomic kinship correlation matrix to obtain a second transformed variable corresponding to each variable; According to the second transformation variable, each univariate mixed model is transformed into a linear regression model consisting of a target label and a residual error; The method of using multivariate mixed model association analysis to locate the main effect gene of the dynamic trait growth trajectory includes: Inferring marker genetic effects based on the least squares method and a linear regression model corresponding to each of the univariate mixed models; Add the test statistics corresponding to the d marker genetic effects to obtain the one-factor, multiple-effect QTN analysis result, where the statistic follows a chi-square distribution with d degrees of freedom; Based on the analysis results of the one-cause-multiple-effect QTN, the major effect gene is located from multiple candidate QTNs.

8. The method according to claim 7, characterized in that The method of locating the major effect gene from multiple candidate QTNs based on the one-cause-multiple-effect QTN analysis results includes: further improving the detection accuracy and efficiency of the major effect gene by considering the linkage disequilibrium between the candidate QTNs.

9. A system for genome association analysis of dynamic traits of livestock and poultry based on a random regression model, characterized in that: include: Model building module, used to build the target random regression model; The genomic association analysis module is used to fit the individual phenotypic trajectory based on the target random regression model and estimate the regression coefficient of the individual phenotypic curve, and then use the regression coefficient as the phenotypic value to locate the main effect gene of the dynamic trait growth trajectory by adopting the multivariate mixed model association analysis method; wherein, the target random regression model includes a first processing layer and a second processing layer, the first processing layer is used to optimize the individual phenotypic curve according to the Bayesian information criterion to reduce the dimensionality of repeated measurements, and the second processing layer is used to eliminate genetic correlation between individuals by using spectral decomposition, and convert the multivariate mixed model into multiple independent univariate mixed models by typical transformation to optimize the solution method of multiple gene variance components.

10. A computer-readable storage medium storing computer instructions, which, when executed by a processor, implement the method for genomic association analysis of dynamic traits of livestock and poultry based on a random regression model as described in any one of claims 1 to 8.

Images