Correlation analysis method of biological phenotype data based on time or space dynamics

Abstract

The invention discloses a correlation analysis method for biological phenotype data based on time or space dynamics. By adopting the method provided by the invention, association analysis can be carried out on biological phenotype data based on different time or space dynamics, associated SNPs based on the time or space dynamics can be obtained through screening, important regulatory genes can be further explored through screening of the associated SNPs, and a gene fine positioning effect on biological complex characters is achieved.

Description

Technical Field

[0001] The invention relates to an association analysis method based on biological phenotype data under temporal or spatial dynamics. Background Art

[0002] Genome-wide association analysis (GWAS) is a commonly used analytical technique for studying complex biological traits. It uses the method of associating genome-wide single nucleotide polymorphism (SNP) data with population phenotypic data one by one to play a role in fine-grained positioning of genes for complex human diseases and biological agronomic traits. Since the birth of this technology in 2006, many GWAS software have been developed, among which the more classic ones are based on the "Q+K model", such as TASSEL and PLINK. In terms of Q, they regard factors such as age, gender, and population structure as covariates, and divide their effects into fixed effects; in terms of K, the model considers the impact of population kinship on individual random effects and constructs random effects based on kinship; the statistical principle of restricted maximum likelihood estimation is used to accurately separate fixed effects, random effects, and the effects of SNPs on biological traits. Through genome-wide SNP calculations, important SNPs with significant associations with phenotypes can be screened out to achieve the goal of genetic analysis of complex traits.

[0003] With the advent of the phenomics era, scientific researchers are paying more and more attention to collecting dynamically growing phenotypic data, including phenotypic data at different time points or in different environments. Based on this type of biological phenotype under temporal or spatial dynamics, in addition to dissecting the SNP effect, researchers also urgently need to calculate the interaction effect between SNP, time point and environment, which plays a key role in helping people understand the regulation of genes under different spatiotemporal states. The existing "Q+K" model can only use biological phenotypic data at a single time point; or treat multiple environments as covariates, which cannot reflect the changing laws of biological phenotypes under continuous environments; or handle multiple related traits. However, for biological phenotypic data with spatiotemporal characteristics, no association analysis method of the "Q+K model" has been developed. Based on the important biological significance of spatiotemporal phenotypes, the present invention discloses an association analysis method based on biological phenotypic data under temporal or spatial dynamics. Summary of the invention

[0004] In view of the above-mentioned defects, the present invention provides a method for association analysis based on biological phenotype data under temporal or spatial dynamics, the method comprising the following steps:

[0005] (1) Combine the biological phenotypic data and covariates of the same organism under temporal or spatial dynamics to construct a mixed model, where the covariates come from the population structure of the organism and the whole genome SNP data; the mixed model contains fixed effects, random effects and dynamic residuals;

[0006] (2) Growth curves were used to model fixed effects, Legendre curves were used to model random effects, and the prestructural dependence model was used to construct the variance structure of the dynamic residuals;

[0007] (3) Using the first-order Taylor expansion to linearize the mixed model, the mixed model is converted into a linear model;

[0008] (4) using the principle of restricted maximum likelihood estimation, deriving an iterative algorithm for combining the average information and maximum expectation of the linearized model of (3), optimizing the random effects and the variance structure after modeling using the iterative algorithm, and obtaining the variance parameters;

[0009] (5) The variance parameters of the whole genome SNP data of the organism are solved according to the methods (1)-(4), and the P value of the SNP data is obtained using the likelihood ratio test and the Wald test method to screen for significant SNPs based on temporal or spatial dynamics.

[0010] In one embodiment, the whole genome SNP data is used to construct the biological kinship matrix K, which is used to reflect the genetic similarity between biological individuals and further used to represent the proportional relationship between different individuals regarding random effects.

[0011] In one embodiment, the random effects include genetic random effects and general random effects.

[0012] In one embodiment, the covariate is regarded as a fixed effect affecting the biological phenotypic data; the biological kinship is regarded as a genetic random effect, and specifically, the biological kinship matrix is ​​regarded as the proportional relationship between different individuals with respect to the genetic random effect;

[0013] In one embodiment, the dynamic residual is the difference between the observed value of the biological phenotype data and the model predicted value; the dynamic residual has autocorrelation under temporal or spatial dynamics, and its law can be described by the structural forward model.

[0014] In one embodiment, step (5) of the analysis method includes substituting the variance parameters into the likelihood function of the mixed model, calculating the likelihood values ​​under the null hypothesis and the alternative hypothesis, respectively, to obtain the statistic of the likelihood ratio test; using the variance parameters to calculate the regression parameters of the fixed effects, and calculating the Wald statistic.

[0015] In one embodiment, the organism is one of an animal, a plant, a fungus, a bacterium, or a virus.

[0016] The method provided by the present invention can be used to carry out association analysis on biological phenotypic data based on different time or space dynamics, and screen for significant SNPs based on time or space dynamics. Through the screening of significant SNPs, important regulatory genes can be further discovered, which plays a role in the precise positioning of genes for complex biological traits. BRIEF DESCRIPTION OF THE DRAWINGS

[0017] The drawings described herein are used to provide a further understanding of the present application and constitute a part of the present application. The illustrative embodiments of the present application and their descriptions are used to explain the present application and do not constitute an improper limitation on the present application. In the drawings:

[0018] Figure 1 This is Example 1 of the present invention regarding the dynamic body weight phenotype trajectory of a mouse population at different time points (Data 1), where the gray and light gray lines represent the dynamic changes in body weight of male and female individuals.

[0019] Figure 2 This is a dynamic derivative diagram of the weight growth of the Gough Island mouse population from 1 to 16 weeks after birth in Example 2 of the present invention (Data 2). The gray and light gray lines represent the dynamic changes in body weight of male and female individuals of different families.

[0020] Figure 3 This is Example 3 of the present invention about the dynamic ground diameter phenotype trajectory of the one-year-old cuttings of the Tibetan Poplar population at different time points (data 3). The gray and light gray lines represent the dynamic changes in the ground diameter phenotype of individuals at high altitudes and low altitudes.

[0021] Figure 4 This is Example 4 of the present invention about the pseudo-dynamic plant height phenotypic trajectory of a nested associated population of corn in different environments; the light gray line represents the pseudo-dynamic changes of different individuals in the associated population with the environmental index, and the thick gray line represents the average dynamic changes of each family in the population.

[0022] Figure 5 The quantile graph (QQ graph) and the Manhattan graph of association analysis of Example 1 of the present invention;

[0023] Figure 6 The quantile map and the Manhattan map of the correlation analysis of the embodiment 2 of the present invention;

[0024] Figure 7 The quantile map and the Manhattan map of the correlation analysis of the third embodiment of the present invention;

[0025] Figure 8 It is the quantile map and the Manhattan map of correlation analysis of Example 4 of the present invention. DETAILED DESCRIPTION

[0026] The specific embodiments of the present invention are described in detail below. It should be understood that the specific embodiments described herein are only used to illustrate and explain the present invention, and are not used to limit the present invention. The present invention provides an association analysis method based on biological phenotype data under temporal or spatial dynamics, and the method comprises the following steps:

[0027] (1) A mixed model is constructed by combining biological phenotypic data and covariates of the same organism under temporal or spatial dynamics, where the covariates come from the population structure of the organism and the whole genome SNP data; the mixed model contains fixed effects, random effects and dynamic residuals.

[0028] Preferably, the whole genome SNP data is used to construct the biological kinship matrix K, which is used to reflect the genetic similarity between biological individuals and further used to represent the proportional relationship between different individuals regarding random effects.

[0029] Specifically, a biological population is obtained, biological phenotypic data of the same biological population under temporal or spatial dynamics is obtained, a population structure of all samples in the population is obtained, and whole-genome SNPs of the population are obtained, and n, d, and p are used to represent the number of samples contained in the population, the dynamic phenotypic dimension, and the number of SNPs, respectively;

[0030] The covariate is constructed as follows : Use the population structure and SNP of all samples in the population to construct the covariates of all individuals. Take individual i as an example, C i =(1,Q i ,SNP ij ) represents its covariate, where Q i represents the population structure of individual i, SNP ij represents the genotype of individual i at the jth SNP. If SNP ij Genotype is qq, Qa, QQ, then SNP ij They are coded as 0, 1, and 2 respectively;

[0031] The kinship matrix is ​​constructed as follows :Based on the SNP test results, M is the genotype matrix, whose dimension is n×p, filled with 0, 1, and 2, and p represents the number of SNPs. i Indicates the i-th column of the M matrix, representing SNP i , that is, all individuals in SNP i genotype on With v m i represents the population mean and variance of SNPi, and the kinship matrix K is constructed as follows:

[0032] (2) Growth curves were used to model fixed effects, Legendre curves were used to model random effects, and the prestructural dependence model was used to construct the variance structure of the dynamic residuals;

[0033] Preferably, the random effects include genetic random effects and general random effects. Covariates are considered as fixed effects that affect biological phenotypic data, and biological kinship is considered as genetic random effects. Specifically, the biological kinship matrix is ​​considered as the proportional relationship between different individuals regarding genetic random effects; the dynamic residual is the difference between the observed value of the phenotypic data and the predicted value of the model; the dynamic residual has autocorrelation under temporal or spatial dynamics, and its law can be described by the structural forward model.

[0034] The fixed effects are constructed as follows :Using growth curve g(C i β;t i ) represents the fixed effects under different time and space states, using β=(β 1 ,β Q ,β SNP ) represents the effect parameter of each component in the covariate, C i β can represent the integrated parameter result of the growth curve, t i Represents a time series. If faced with time-dynamic phenotypic data, Legendreorthogonal polynomials (LOPs), Logistic equation, Gompertz equation, etc. are used to construct g(C i β;t i ), if faced with spatially dynamic phenotypic data, calculate the average phenotypic value under each environment as t i , using the allometric equation to construct g(C i β;t i );

[0035] The genetic random effect was constructed as follows :Use r 1 (a i ;t i ) represents the heritable effect function under different temporal or spatial dynamics, a i is the model parameter, constructed using the LOPs equation, and the formula is as follows: Where nr1 represents the order of the LOPs equation, set nr1 = 4, according to t i Normalized, t min Represents t i The minimum value of t max Represents t i The maximum value of φ 0 (t i )=1, remember but

[0036] The general random effects method is as follows :Use r 2 (p i ;t i ) represents the effect function under different time or space dynamics, p i is the model parameter, and the LOPs equation is also used to construct r 2 (p i ;t i ), nr2 represents the order of the LOPs equation, and nr2 = 4 is set, then

[0037] The mixed model is constructed as follows : Use y i Represents the dynamic phenotype of individual i, using e i represents the dynamic residual of individual i, then y i It can be expressed as y i =g(C i β;t i )+r 1 (a i ;t i )+r 2 (p i ;t i )+e i ;

[0038] (3) The mixed model is linearized using the first-order Taylor expansion to transform the mixed model into a linear model; the random effects and residuals in the linear model have a mean of zero and variances of The structure of , where K is the kinship matrix, I n is the identity matrix, represents the Kronecker product, G 0 and P 0 The matrix is ​​constructed from the variance-covariance of the Legendre curves of different orders, R 0 The variance matrix representing the dynamic residual parameters is constructed from the structural pre-dependence model.

[0039] Specifically, the first-order Taylor expansion is used to linearize the hybrid model, and we can get Right now in, It is the form of the first-order derivative of each parameter in β for the fixed effect part. Considering all n individuals in the group, the first-order Taylor linearized mixed model can be further expressed as follows:

[0040]

[0041] Or it can be abbreviated as Where a, p, and e have a mean of zero and their variances are expressed as Where K is the correlation matrix between individuals, I n is the identity matrix, G 0 With P 0 are the variance matrices of genetic random effects and general random effects, respectively. The variance-covariance matrix of each order of LOPs is used to construct the R 0 The first-order structural pre-dependence model is used. Specifically:

[0042] G 0 =diag(Vg [1] ,Vg [2] ,…,Vg [nr1] ),

[0043] P 0 =diag(Vp [1] ,Vp [2] ,…,Vp [nr2] ),

[0044]

[0045] R 0 In structure represents the variance at each time point, and φ represents the correlation coefficient between different time points.

[0046] Further decomposition of K, we get K = U k D k U′ k , where D k is a diagonal matrix, and all elements on the diagonal are represented by (δ 1 ,δ 2 ,...,δ n ) Filling, U k is the eigenvector matrix, which is obtained by filling the eigenvectors into each column of the matrix in the order of the eigenvalues.

[0047] Further use Linearized Mixed Model Perform left multiplication, and we get remember

[0048] The above formula can be rewritten as:

[0049] Y=Xβ+Z 1 a+Z2 p+e (E1)At this time, the variance of Y can be written as:

[0050]

[0051] at this time The variance of

[0052] (4) using the principle of restricted maximum likelihood estimation, deriving an iterative algorithm for combining the average information and maximum expectation of the linearized model of (3), optimizing the random effects and the variance structure after modeling using the iterative algorithm, and obtaining the variance parameters;

[0053] In detail, based on the principle of restricted maximum likelihood estimation, the restricted maximum likelihood function is used to derive G 0 , P 0 , R 0 An iterative algorithm that combines the average information about the parameters in the matrix with the maximum expectation to calculate all variance parameters;

[0054] Restricted Maximum Likelihood Function Expressed as Where P = V -1 -V -1 X(X′V -1 X) -1 X′V -1 . Let θ represent the set of all unknown parameters, including G 0 ,P 0 and R 0 The unknown parameters in

[0055]

[0056] The formula for the iterative algorithm combining average information and maximum expectation is

[0057]

[0058] f is the current iteration number, f+1 represents the update of the next iteration, represents the set of first-order derivatives of all L with respect to all parameters, AE represents the combined information matrix, which is a weighted combination of the average information matrix (AI) and the expectation-maximum matrix (EM), that is, AE = (1-λ)AI + λEM, where λ is the weight coefficient. If λ = 0, the combined algorithm is completely switched to the AI ​​iterative algorithm; if λ = 1, the combined algorithm is completely switched to the EM iterative algorithm.

[0059] The iterative algorithm for combining average information and maximum expectation is described as follows: Includes the construction of first-order derivatives, the construction of AI matrices, the construction of EM matrices, the construction of AE matrices, and the iterative calculation steps using them.

[0060] The first-order derivative is constructed as follows :

[0061] The first-order derivative is expanded as follows:

[0062]

[0063] The formula (E3) with respect to Vg is given as follows: [j] , Vp [j] , The calculation method of the first-order derivative of and φ is:

[0064] Where W = (XZ 1 Z 2 ), The form is:

[0065]

[0066] The AI ​​matrix is ​​constructed as follows:

[0067] The number of rows and columns of the AI ​​matrix is ​​((nr 1 +1)+(nr 2 +1)+(d+1)), the element at the [i,j] position is calculated as follows:

[0068]

[0069] For most of the [i, j] elements, the last two terms of the above formula are equal to 0. Therefore, we can define F = [f 1 f 2 ... f J ], but

[0070]

[0071] but and is not equal to 0. Therefore, the (d+1)×(d+1) part of the matrix at the lower right corner of the AI ​​matrix needs to additionally calculate the second and third terms of formula (E9).

[0072] Construct f as follows j The vector is used to fill the jth column of the F matrix mentioned in formula (E10), and its form is as follows:

[0073]

[0074] For convenience, take d = 5 as an example, and H i (i=1,2,...,d) and J represents R 0 right and the first-order derivative of φ. That is The specific form is as follows:

[0075]

[0076]

[0077] For the lower right (d+1)×(d+1) part of the AI ​​matrix, we get and The second-order partial derivative form of , still taking d = 5 as an example,

[0078]

[0079]

[0080] Using equations (E17)-(E22), the second and third terms of equation (E9) can be calculated, and then combined with the AI ​​of equation (E10) [j,k] Add them together to get the complete AI matrix.

[0081] The EM matrix is ​​constructed as follows:

[0082]

[0083] The AE matrix is ​​constructed as follows :

[0084] Introduce the λ parameter and add the weights of AI and EM matrices to make

[0085] AE=(1-λ)AI+λEM(E28)

[0086] The iterative steps of the combination algorithm are as follows :

[0087] λ is tested between [0,1] by increasing the step size by 0.05 until the G constructed by the iterative parameter 0 ,P 0 and R 0 Until all the matrices are positive definite. The iteration formula of the AI+EM combination algorithm is as shown in E1. During the iteration process, the following formula is used to update

[0088]

[0089] To summarize, the iterative steps of the calculation are summarized as follows:

[0090] (I) Setting the initial value of θ;

[0091] (II) Construct a linearized mixed model according to (E1);

[0092] (III) Update according to formula (E29)-(E32)

[0093] (IV) Calculate the first-order derivative vector according to formulas (E4)-(E8);

[0094] (V) Construct the AI ​​matrix according to formulas (E9)-(E22);

[0095] (VI) Set λ=0 and try the AI ​​algorithm according to formula (E2); if the parameter falls into an unreasonable range during the iteration, increase the value of λ by a step of 0.05 and switch to the AE combination algorithm according to formula (E28);

[0096] Repeat steps (II)-(VI) until the iterations converge and the parameters are finally stable.

[0097] (5) The variance parameters of the whole genome SNP data of the organism are solved according to the methods (1)-(4), the P value of the SNP data is obtained by using the likelihood ratio test and the Wald test method, and the significant SNPs based on time or space dynamics are screened. Preferably, step (5) includes substituting the variance parameters into the likelihood function of the mixed model, calculating the likelihood values ​​under the null hypothesis and the alternative hypothesis, and obtaining the statistic of the likelihood ratio test; using the variance parameters to calculate the regression parameters of the fixed effect, and calculating the Wald statistic.

[0098] Specifically, for the whole genome SNPs of the population, the variance parameters of each SNP are calculated according to the above method, and the p-value is obtained by using the Z-wald and likelihood ratio test methods, and the significant SNPs based on time or space dynamics are screened out in the whole genome; in order to reasonably display the analysis results, the Manhattan plot and the quantile plot (Quantile-quantile plot, QQ plot) are used to display the p-value results of the whole genome SNPs. The Manhattan plot is a scatter plot drawn with the value of -log10(P), which shows the significance of different SNPs on different chromosomes. After the p-value on the QQ plot is transformed by the logarithm of -log10, the non-significant SNP is usually fitted to the y=x straight line in the QQ plot, while the significant site will show the phenomenon of the tail being upturned.

[0099] In order to make the purpose, technical methods and advantages of this application clearer, the technical solution of this application will be clearly and completely described below in combination with the specific embodiments of this application and the corresponding drawings. Obviously, the described embodiments are only part of the embodiments of this application, not all of them. Based on the embodiments in this application, all other embodiments obtained by ordinary technicians in this field without creative work are within the scope of protection of this application.

[0100] The following embodiments all adopt the above-mentioned association analysis method based on biological phenotype data under temporal or spatial dynamics provided by the present invention.

[0101] Example

[0102] Example 1 is used to illustrate the association analysis of biological phenotype data based on time dynamics

[0103] The data comes from http: / / phenome.jax.org / db / q? rtn=projects / projdet&reqprojid=539, involving the dynamic weight data of mice from Goff Island, South Atlantic, 1-16 weeks after birth. After SNP quality filtering, 11,833 SNPs were obtained. Considering the discarding of individuals with phenotypic missing greater than or equal to 4 times, the full data of 16 weights of 1212 strains of mice were obtained. The weaning time of mice was recorded during the construction and maintenance of the group. Therefore, in this embodiment, family information, gender, weaning time, etc. are considered as covariates, and the weight changes of 1212 strains of mice from 1-16 weeks after birth are integrated as phenotypic data, and 11,833 SNPs are used as genotypic data. The LOP curve is selected to construct a fixed effect and conduct GWAS analysis.

[0104] The results of Example 2 show that Figure 1 It was shown that the body weight dynamics of a mouse population does not follow a specific growth curve and can be modeled by a nonparametric approach using the Legendre formula. Figure 5 The tail of the QQ graph in A is upturned, and the Manhattan graph ( Figure 5 B) revealed that significant SNPs with high -logP values ​​were concentrated in the 16-21Mb range of chromosome 6 and the 118-129M range of chromosome 10, which may represent quantitative trait loci that affect the dynamic changes in mouse body weight.

[0105] Example 2 is used to illustrate the association analysis of biological phenotype data based on time dynamics

[0106] The data sources are the same as those in Example 1, using the same population and the same SNP data, but using the dynamic derivative of mouse body weight during the 1-16 week time series growth rate. Therefore, the dynamic changes in growth rate were used instead of the original body weight data in Example 1 to carry out subsequent data analysis, and the LOP curve was selected to construct the fixed effect.

[0107] The method of the present invention is used to analyze Figure 2 , Figure 2 It also shows that its derivative curve does not follow a specific growth curve law, but its shape is similar to Figure 1 The difference is large, and the non-parametric method of Legendre formula is also used to model it. Figure 6The tail of the QQ graph in A is upturned, and the Manhattan graph ( Figure 6 B) reveals that significant SNPs with high -logP values ​​are concentrated in the 115-129M interval of chromosome 10, which overlaps with the association analysis results of data 1, indicating that the results of the association analysis of weight traits based on the Legendre curve of weight derivative are equally reliable.

[0108] Example 3 is used to illustrate the association analysis based on biological phenotype data under spatial dynamics

[0109] The dynamic change data of the ground diameter of the one-year-old cuttings of Tibetan Poplar were used. The data came from the master's thesis of Wei Xueda of Beijing Forestry University in 2023, "Analysis of the Genetic Regulation Mechanism of Tibetan Poplar Growth Based on Systematic Mapping of Allometric Growth". The Tibetan Poplar population contains 200 lines. The ground diameter phenotype was investigated every 2 weeks 40 days after the branch cuttings. A total of 10 phenotypic data were obtained. The population has 24,949 SNP data. This example uses the high and low altitudes of the original sampling of Tibetan Poplar as covariates, and chooses to use the Gompertz curve to construct fixed effects.

[0110] The results of Example 3 show Figure 3 It shows that the ground diameter dynamic data of Tibetan Populus follows the growth curve of Gompertz. Figure 7 The tail of the QQ plot in A is upturned, and more than 95% of the SNPs have p-values ​​that match the expected values. Manhattan plot ( Figure 7 B) Nine significant SNPs were revealed, scattered on chromosomes 1, 4, 7, 9, 10, 11, and 19, which may represent quantitative trait loci affecting terrestrial diameter growth.

[0111] Example 4 is used to illustrate the association analysis based on biological phenotype data under spatial dynamics

[0112] The nested association mapping population (NAM) of maize containing 25 families constructed by Buckler et al. was used. The SNP data of the population was downloaded from https: / / cbsusrv04.tc.cornell.edu / users / panzea / download.aspx?filegroupid=5, and the phenotypic data of the population was downloaded from http: / / www.panzea.org / lit / data_sets.html / #phenos. 99,190 SNPs were obtained after low-quality filtering. Different locations and seasons were regarded as independent environments. 4558 individuals with both SNP sequencing and plant height phenotypic data at 8 different locations were screened out. In this example, families and SNPs were used to construct covariates, allometric growth curves were used to construct fixed effects, and the population phenotypic mean of each environment was calculated as the environmental index to replace t in model (1). i .

[0113] Example 4 The results show that: Figure 1-3 no the same, Figure 4 It shows the law of biological phenotype change with the environment: it grows exponentially with the environmental index, so the allometric curve is used to model this. Figure 8 The tail of the QQ graph in A is upturned, and the degree of upturn is greater than that of the first three examples, indicating that there are more significant SNPs. Manhattan graph ( Figure 8 B) revealed a significant SNP with higher -logP on chromosome 9, which is consistent with the results of the author's original analysis of the quantitative trait loci of plant height, indicating that the association analysis model proposed in the present invention has strong statistical power, and the significant SNP revealed can be used for further functional gene mining.

[0114] In summary, Figure 1-Figure 8 The analysis results of the four examples are shown. It can be seen that the growth dynamics of the four examples are quite different. Example 1 and Example 2 use the fourth-order Legendre model to draw their growth dynamics ( Figure 1 , Figure 2 ), Example 3 uses the Gompertz curve to draw the growth dynamics ( Figure 3 ), Example 4 uses allometric growth curves between different environments to draw its pseudo-dynamic growth ( Figure 4 ). In the QQ graphs of the four embodiments, most of the SNPs are close to the y=x line, indicating that the analysis method provided by the present invention has better control over random drift type SNPs, indicating a lower false positive, and the tail of the QQ graph is relatively satisfactory, indicating that the analysis method provided by the present invention has a stronger discovery rate for significant SNPs ( Figure 5 A, Figure 6A, Figure 7 A, Figure 8 A); Manhattan plots show both sparse and concentrated significant SNPs ( Figure 5 B, Figure 6 B, Figure 7 B, Figure 8 B). Therefore, the present invention has a good model effect for association analysis using dynamic growth data, and shows dynamic characteristics in the field of GWAS.

Claims

1. A method for association analysis based on biological phenotypic data under temporal or spatial dynamics, characterized in that: The method comprises the following steps: (1) constructing a mixed model by combining biological phenotypic data and covariates of the same organism under temporal or spatial dynamics, wherein the covariates are derived from the population structure and whole genome SNP data of the organism, and the mixed model comprises fixed effects, random effects and dynamic residuals; (2) using growth curves to model the fixed effects, using Legendre curves to model the random effects, and using a pre-structural dependency model to construct the variance structure of the dynamic residual; (3) linearizing the mixed model using a first-order Taylor expansion to convert the mixed model into a linear model; (4) using the principle of restricted maximum likelihood estimation, deriving an iterative algorithm for combining the average information and maximum expectation of the linearized model of (3), and using the iterative algorithm to optimize the random effects and the variance structure after modeling to obtain variance parameters; (5) The variance parameters of the whole genome SNP data of the organism are solved according to the methods (1)-(4), and the P value of the SNP data is obtained using the likelihood ratio test and the Wald test method, and the significant SNPs based on time or space dynamics are screened.

2. The method according to claim 1, characterized in that The whole genome SNP data is used to construct a kinship matrix K of the organism.

3. The method according to claim 1, characterized in that The random effects include genetic random effects and general random effects.

4. The method according to claim 2, characterized in that: The covariates are considered as fixed effects that affect the phenotypic data of the organisms, and the relatedness of the organisms is considered as a genetic random effect.

5. The method according to claim 1, characterized in that The dynamic residual is the difference between the observed value of the biological phenotype data and the model predicted value.

6. The method according to claim 1, characterized in that Step (5) includes substituting the variance parameter into the likelihood function of the mixed model, calculating the likelihood values ​​under the null hypothesis and the alternative hypothesis respectively, and obtaining the statistic of the likelihood ratio test; using the variance parameter to calculate the regression parameter of the fixed effect and calculate the Wald statistic.

7. The method according to claim 1, characterized in that The organism is one of animals, plants, fungi, bacteria and viruses.

Images