A method for gene-environment interaction analysis that can correct population structure stratification

Abstract

The application relates to the technical field of biostatistics, in particular to a gene-environment interaction analysis method capable of correcting population structure stratification, obtaining phenotype data, genotype data, environmental factor data and confounding factor data of samples to form sample data; constructing a regression model based on the sample data, fitting a constraint model with a marginal genetic effect and a marginal gene-environment interaction both being 0 and calculating corresponding residual errors; estimating a site genetic variation rate at an individual level based on a genetic principal component for a genetic locus to be tested; adopting a normal distribution approximation method to calculate a p value based on a score statistic to test the significance of a marginal genetic effect of the site; and selecting and calculating a test statistic for testing the significance of a marginal gene-environment interaction of the site according to the significance of the marginal genetic effect; and the method can avoid over-high false positive rates or false negative rates in analysis results caused by population confounding.

Description

Technical Field

[0001] This invention relates to the field of biostatistics, and more specifically to a gene-environment interaction analysis method that can correct for population structure stratification. Background Technology

[0002] With the development of high-throughput sequencing technology and electronic health record systems, health characteristics such as whole-genome genetic data, clinical trial indicators, complex disease diagnoses, and personal exercise records are entering an era of information and digitalization. Currently, many large biobanks have emerged globally, and their biomedical big data provide abundant research resources for systems biology research and genome-wide association studies.

[0003] Genome-wide association studies (GWAS) have identified tens of thousands of genetic variations associated with complex human diseases. Complex traits or diseases are typically influenced by both genetic and environmental factors, including a wide range of environmental factors such as living environment, biological characteristics, lifestyle, medical practices, and biomarkers. Studying the impact of gene-environment interactions on complex traits or diseases helps to understand the mechanisms underlying these traits or diseases, and has significant theoretical implications and broad application prospects for fields such as precision medicine, personalized prevention, and intelligent diagnosis and treatment of complex diseases.

[0004] In recent years, phenotypic and genetic data in large biobanks have exhibited new characteristics. For phenotypic data, the prevalence of most diseases in large biobanks is less than 1%, resulting in a high degree of imbalance in phenotypic distribution. With the improvement of information technology, phenotypic data is no longer limited to traditional continuous or binary data; increasingly complex phenotypic data are being applied in genome-wide association studies. However, effective analytical methods for complex phenotypic data remain relatively scarce.

[0005] Population stratification or demographic structure is a significant confounding factor in genome-wide association studies (GWAS). Large biobanks may contain a large number of individuals from heterogeneous or mixed populations. The genetic variation rates of genetic loci and the distribution of phenotypes often differ among different populations. Due to the confounding effect of population structure, individuals from mixed populations are typically excluded from the analysis.

[0006] Traditional normal distribution approximation methods can lead to a high Type I error rate in imbalanced phenotypic distributions and loci with low genetic variation rates, resulting in an excessively high false positive rate in statistical inference. To address the inaccuracy of the normal distribution approximation, several genome-wide association analysis (GWAS) algorithms based on the saddlepoint approximation (SPA) method have been designed and developed. The core idea of ​​the saddlepoint approximation method is to use the cumulant generating function to estimate the cumulative distribution function of the distribution. Because it uses more higher-order moment information, it is more accurate than the traditional normal distribution approximation method. With the continuous deepening of theoretical research on the saddlepoint approximation method, several fast and accurate GWAS analysis algorithms based on the saddlepoint approximation method have been designed and developed.

[0007] Existing methods for analyzing gene-environment interactions still have limitations due to the challenges posed by large datasets, rare genetic variations, and imbalanced phenotypic distributions. Specifically, conventional analysis methods (such as score tests, Wald tests, and likelihood ratio tests) are computationally inefficient, requiring significant computation time for genome-wide gene-environment interaction analysis. Currently, there is a lack of fast and accurate algorithms for analyzing gene-environment interactions in multi-ethnic or mixed populations within large-scale data analytics.

[0008] In summary, big data analysis of gene-environment interactions urgently requires a fast, accurate analysis method that is applicable to complex phenotypes and mixed populations, in order to avoid excessively high false positive or false negative rates in the analysis results due to the mixed population. Summary of the Invention

[0009] The purpose of this invention is to provide a gene-environment interaction analysis method that can correct for population structure stratification. It is suitable for the analysis of mixed populations, can control the Type I error rate to not exceed a given significance level, and has good statistical power, avoiding excessively high false positive or false negative rates in the analysis results.

[0010] The objective of this invention is achieved through the following technical solution:

[0011] A gene-environment interaction analysis method that can correct for population structure stratification is proposed. This method can be applied to genome-wide association analysis of complex phenotypes such as continuous phenotypes, binary phenotypes, survival data phenotypes, multi-category phenotypes, and longitudinal data phenotypes in mixed populations. The method includes the following steps:

[0012] Step S1: Obtain phenotypic data, genotypic data, environmental factor data, and confounding factor data; confounding factor data includes age, sex, and principal genetic components, etc.

[0013] Step S2: Construct a regression model based on the sample data, fit a constrained model where both marginal genetic effect and marginal gene-environment interaction are 0, and calculate the corresponding residuals;

[0014] Step S3: For the genetic locus to be tested, estimate the locus genetic variation rate at the individual level based on the principal components of the genetic data;

[0015] Step S4: Calculate the p-value using the normal distribution approximation method based on the score statistic to test the significance of the marginal genetic effect at the locus;

[0016] Step S5: Select and calculate the test statistic for the significance of marginal gene-environment interactions at the loci based on the significance of the marginal genetic effect;

[0017] Step S6: For the test statistic of marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of marginal gene-environment interaction.

[0018] Step S7: Calculate the p-value using the normal distribution approximation method based on the ancestry-specific score statistic to test the significance of the ancestry-specific marginal genetic effect of the locus;

[0019] Step S8: Select and calculate the test statistic for the significance of ancestry-specific marginal genetic effects for testing loci based on the significance of the marginal genetic effects.

[0020] Step S9: For the test statistic of the marginal gene-environment interaction of ancestry specificity, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of the gene-environment interaction of ancestry specificity.

[0021] Step S10: Using the Cauchy combination method, the statistical p-value for testing the significance of gene-environment interactions and the statistical p-value for testing the significance of ancestry-specific gene-environment interactions across all ancestry types are combined to obtain the final statistical p-value, thereby achieving the analysis of gene-environment interactions across the entire genome.

[0022] Furthermore, the construction of a regression model based on sample data, fitting a constrained model where both marginal genetic effects and marginal gene-environment interactions are zero, and calculating the corresponding residuals specifically includes:

[0023] Assume the sample size of the participants in the study is n, and these n individuals can come from different ethnic groups or subgroups. The individuals are not related by blood. For individual i, 1 ≤ i ≤ n, let X i Let G represent an m-dimensional confounding factor vector. iIndicates the genotype of the locus, with values ​​of 0, 1, or 2, E i Y represents environmental factors. i To represent a certain phenotype, for different complex traits, corresponding regression models are used for analysis. The linear predictor term in the regression model is... Where, β X Indicates confounding factor X i The m-dimensional coefficient vector, β E β G and β G×E These represent environmental factors E. i Genotype G i and gene-environment interaction term G i E i The corresponding coefficient;

[0024] To test the significance of the marginal gene-environment interaction at the locus, i.e., to test the null hypothesis H0: β G×E =0, under the assumption H c :β G =β G×E =0, meaning the marginal genetic effect and marginal gene-environment interaction at the locus are both zero. The model is fitted under this constraint to obtain the maximum likelihood estimate of the model parameters under this constraint, and the n-dimensional residual vector R = (R1,…,R) under the constraint model is calculated. n ) T ;

[0025] In step S3, for the locus to be tested, the individual-level locus genetic variation rate is estimated based on the principal components of the genetic data, specifically including:

[0026] For the locus to be tested, assume that the individual-level genetic variation vector of n individuals is q = (q1, ..., qn). n ) T , where q i This represents the rate of genetic variation of individual i at that locus;

[0027] Let G = (G1, ..., G n ) T The genotype vector representing the locus, let It represents an n×(D+1) dimensional genetic principal component matrix containing column vectors with all elements equal to 1 and the first D genetic principal component vectors containing all population structure information;

[0028] Based on the genetic principal components, the response variable and explanatory variable were fitted to G and the genetic principal component X, respectively. PC The linear regression model calculates an estimate of q. This yields an estimate of q. in, I(.) denotes the indicator function;

[0029] Based on the genetic principal components, the response variable and explanatory variables were fitted as (I(G1≥0.5),…,I(G1≥0.5), respectively. n ≥0.5)) T and principal component X PC The logistic regression model can yield an estimate of q. in, The function σ(x) = exp(x) / (1 + exp(x)) This represents the linear predictor term in the logistic regression model. Maximum likelihood estimation;

[0030] To balance the accuracy of the estimation with computational efficiency, if If the proportion of elements in the interval [0,1] exceeds a pre-defined positive number between 0 and 1 (the default value is 0.9), then use... As an estimate of q, otherwise use As an estimate of q.

[0031] In step S4, the calculation of the p-value based on the score statistic using the normal distribution approximation method to test the significance of the marginal genetic effect at the locus specifically includes:

[0032] Calculate the score test statistic used to test for marginal genetic effects.

[0033]

[0034] And it is believed that the test statistic Assuming H G :β G When = 0 is true, it follows the expected value. variance is The normal distribution of, let As a method for testing hypothesis H using the normal distribution approximation G :β G =0 yields the two-sided p-values, where for The observed values ​​are Φ(.), which represents the cumulative distribution function of the standard normal distribution.

[0035] In step S5, the step of selecting and calculating the test statistic for testing the significance of marginal gene-environment interactions at a test site based on the significance of marginal genetic effects specifically includes:

[0036] Testing the significance of marginal gene-environment interactions is equivalent to testing the null hypothesis H0: β G×E Is the expression = 0 true?

[0037] Let ∈ be a given positive number between 0 and 1;

[0038] If we consider hypothesis H G :β G If the p-value obtained by testing = 0 is less than or equal to ∈, then use As a test of the null hypothesis H0:β G×E The test statistic for = 0, where

[0039] If we consider hypothesis H G :β G If the p-value obtained by testing =0 is greater than ∈, then the Wald test (or likelihood ratio test) and the test statistic are used respectively. Regarding the null hypothesis H0:β G×E =0 is used for testing, where To adapt to the residual vector of genotype, I n Let W be an n-order identity matrix, then W = (1 n G) is an n×2 dimensional matrix, 1 n =(1,1,…,1) T For an n-dimensional column vector containing only 1s, we will use The two p-values ​​obtained from the Wald test (or likelihood ratio test) are combined using the Cauchy combination (CCT) method to obtain the final p-value.

[0040] In step S6, the test statistic for testing marginal gene-environment interactions is calculated using a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method to determine the significance of the marginal gene-environment interaction. Specifically, this includes:

[0041] Assume the genotype vector of the sample at this locus is G = (G1, ..., G...). n ) T Fitting β G =β G×E The residual vector obtained from the constraint model with 0 is R = (R1, ..., R2). n ) T ;

[0042] Under the Hardy-Weinberg equilibrium hypothesis in genetics, the locus genotype G of each individual is... i ,1≤i≤n are considered to be mutually independent and follow a binomial distribution Binom(2,q) i A random variable q iThis represents the allele frequency of the genetic variation of the i-th individual at that locus.

[0043] Based on the locus genetic variation rate at the individual level, an approximate test statistic is used to evaluate the null hypothesis H0:β. G×E The statistical p-value is calculated using the conditional distribution when = 0 is true.

[0044] When using the normal distribution approximation method to calculate the statistical p-value, the specific steps are as follows:

[0045] If the test statistic is S G×E Given (R, E, λ), the normal distribution approximation method uses an expectation of variance is The normal distribution of the null hypothesis H0:β G×E When = 0 holds true, the statistic S G×E The distribution is approximated, and the two-sided p-value is calculated using the normal distribution approximation method. Where E = (E1, ..., E n ) T s represents the environmental factor vector. G×E For S G×E The observed values, where Φ(.) represents the cumulative distribution function of the standard normal distribution;

[0046] If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation and variance as follows: and The normal distribution of the null hypothesis H0:β G×E =0 statistic The distribution of is approximated, and the two-sided p-value obtained using the normal distribution approximation method is: in for The observed values.

[0047] When using the saddle point approximation method to calculate the statistical p-value, the specific steps are as follows:

[0048] If the test statistic is S G×E Then, under the null hypothesis H0:β G×E When = 0 is true, G i The i = 1, ..., n are considered to be independent and follow a binomial distribution Binom(2, q i Given a sequence of random variables (R, E, λ), in order to estimate S... G×E Given the cumulant generating function when the null hypothesis H0 holds, first estimate G. i Moment generating function:

[0049]

[0050] the first-order derivative and second-order derivative are respectively:

[0051]

[0052] G i the estimate of the cumulant generating function is the first-order derivative and second-order derivative thereof are respectively:

[0053]

[0054] Therefore, when H0 holds, S G×E the estimate of the cumulant generating function is:

[0055]

[0056] the first-order derivative and second-order derivative thereof are respectively:

[0057]

[0058] For any given real number s0, first calculate the value that satisfies ζ, then calculate to obtain and To approximate the distribution of S G×E under the null hypothesis H0, according to the Barndorff-Nielsen saddlepoint approximation formula, is used to approximately estimate the probability Pr(S G×E <s0|(R,E,λ)). For a given observed value s G×E of the test statistic S G×E , the two-sided p-value calculated by the above saddlepoint approximation algorithm is p l +p r , where the left p-value p l and the right p-value p r are respectively:

[0059]

[0060] wherein and respectively represent the probabilities obtained by the above saddlepoint approximation algorithm and estimates;

[0061] If the test statistic is then when the null hypothesis H0:β G×E = 0 holds, substitute G iThe i = 1, ..., n are considered to be independent and follow a binomial distribution Binom(2, q i A sequence of random variables, in Given the conditions, The estimate of the cumulant generating function when the null hypothesis H0 holds is:

[0062]

[0063] Its first and second derivatives are as follows:

[0064]

[0065] For any given real number First calculate to make Established Then the calculation yields... and use To approximate the probability For a given Observations The bilateral p-value calculated using the saddle point approximation algorithm described above is p. l +p r Let p be the value on the left side. l and the p-value on the right. r They are respectively:

[0066]

[0067] in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate.

[0068] The statistical p-value is calculated using a hybrid test strategy that combines the normal distribution approximation method with the saddle point approximation method, as detailed below:

[0069] Let r be a pre-selected positive number;

[0070] If we test the null hypothesis H0: β G×E The statistic for = 0 is S G×E s G×E For S G×E The observed values, if The normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value.

[0071] If we test the null hypothesis H0: β G×EThe statistic for =0 is for The observed values, if If the p-value is calculated using the normal distribution approximation method, then the p-value is calculated using the saddle point approximation method; otherwise, the p-value is calculated using the saddle point approximation method.

[0072] In step S7, the calculation of the p-value based on the ancestry-specific score statistic using a normal distribution approximation method to test the significance of the ancestry-specific marginal genetic effect of the locus specifically includes:

[0073] Suppose that the n individuals in the research cohort come from K ancestral populations, let G = (G1, ..., G2) n ) T Genotype vector representing a locus. Represents the ancestry-specific genotype vector of the kth ancestry (ancestral population);

[0074] For individual i, i≤n, let This represents the number of haplotypes in the k-th ancestry (ancestral population) at that locus, i.e., the local ancestry count. Let the vector representing the corresponding local lineage count be... An estimate of the ancestry-specific allele frequencies of the k-th ancestry (ancestral population);

[0075] make Let represent the ancestry-specific marginal genetic effect of the k-th ancestry (ancestral population), and let Let the ancestry-specific marginal gene-environment interaction of the k-th ancestry (ancestral population) be represented. Then, the linear prediction term in the regression model can be rewritten as follows:

[0076] Testing the significance of the marginal genetic effect of ancestry specificity is equivalent to testing the hypothesis. Is it true?

[0077] Calculate the score test statistic for testing the marginal genetic effect specific to ancestry.

[0078]

[0079] And it is believed that the test statistic In the assumption When established, it obeys the expected value. variance is The normal distribution of, let The two-sided p-value obtained by using the normal distribution approximation method to test the significance of the marginal genetic effect of lineage specificity, where for The observed values, where Φ(.) represents the cumulative distribution function of the standard normal distribution;

[0080] In step S8, selecting and calculating the test statistic for the significance of ancestry-specific marginal genetic effects at the test loci based on the significance of ancestry-specific marginal genetic effects specifically includes:

[0081] Testing the significance of ancestry-specific marginal gene-environment interactions is equivalent to testing the hypothesis. Is it valid?

[0082] Let ∈ be a given positive number between 0 and 1;

[0083] If we consider the hypothesis If the p-value obtained from the test is less than or equal to ∈, then use As a hypothesis to be tested The test statistic, of which

[0084] If we consider the hypothesis If the p-value obtained from the test is greater than ∈, then the Wald test (or likelihood ratio test) and the test statistic are used respectively. For the hypothesis The test was conducted, among which To adapt to the genotype-specific residual vector, I n W is an n-order identity matrix. (k) =(1 n G (k) ) is an n×2 dimensional matrix, 1 n =(1,1,…,1) T For an n-dimensional column vector containing only 1s, we will use Wald test (or likelihood ratio test) hypothesis The two obtained p-values ​​are combined using the Cauchy combination method to obtain the final p-value;

[0085] In step S9, the test statistic for testing the marginal gene-environment interaction specific to ancestry is calculated using a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method to determine the significance of the gene-environment interaction specific to ancestry. Specifically, this includes:

[0086] When using the normal distribution approximation method to calculate the statistical p-value, the specific steps are as follows:

[0087] If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation of variance is The normal distribution is related to the hypothesis. Statistics at the time of establishment If the distribution is approximated, then the two-sided p-value calculated using the normal distribution approximation method is... in for The observed values, where Φ(.) represents the cumulative distribution function of the standard normal distribution;

[0088] If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation and variance as follows: and The normal distribution is related to the hypothesis. Statistics at the time of establishment If the distribution is approximated, then the two-sided p-value obtained using the normal distribution approximation method is: in for The observed values.

[0089] When using the saddle point approximation method to calculate the statistical p-value, the specific steps are as follows:

[0090] If the test statistic is Then in the assumption When it was established Consider them as independent entities that follow a binomial distribution, Binom(2,q) (k) A sequence of random variables, in Given the conditions, in order to estimate In the assumption The cumulative generating function at the time of its establishment is first estimated. Moment generating function:

[0091]

[0092] The estimate of the cumulant generating function is: Its first and second derivatives are as follows:

[0093]

[0094] Therefore, in When it was established The estimate of the cumulant generating function is:

[0095]

[0096] Its first and second derivatives are as follows:

[0097]

[0098] For any given real number s0, first calculate such that The established ζ (k) Then calculate to get and To approximate Under the original hypothesis The distribution below, according to the Barndorff-Nielsen saddle point approximation formula, uses... To approximate the probability Then for a given test statistic Observations The bilateral p-value calculated using the saddle point approximation algorithm described above is p. l +p r Let p be the value on the left side. l and the p-value on the right. r They are respectively:

[0099]

[0100] in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate;

[0101] If the test statistic is Then in the assumption When it was established Consider them as independent entities that follow a binomial distribution, Binom(2,q) (k) A sequence of random variables, in Given the conditions, In the assumption The estimate of the cumulant generating function when it holds true is:

[0102]

[0103] Its first and second derivatives are as follows:

[0104]

[0105] For any given real number First calculate to make Established Then the calculation yields... and use To approximate the probability For a given Observations The bilateral p-value calculated using the saddle point approximation algorithm described above is p. l +p r Let p be the value on the left side. l and the p-value on the right. r They are respectively:

[0106]

[0107] in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate.

[0108] The statistical p-value is calculated using a hybrid test strategy that combines the normal distribution approximation method with the saddle point approximation method, as detailed below:

[0109] Let r be a pre-selected positive number;

[0110] If we test the null hypothesis The statistic is for The observed values, if The normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value.

[0111] If we test the null hypothesis The statistic is for The observed values, if The normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value.

[0112] In step S10, the Cauchy combination method combines the statistical p-value for testing the significance of gene-environment interactions with the statistical p-value for testing the significance of ancestry-specific gene-environment interactions across all ancestry levels to obtain the final statistical p-value, thereby achieving genome-wide gene-environment interaction analysis. Specifically, this includes:

[0113] The significance of the marginal gene-environment interaction at the locus will be tested, i.e., the null hypothesis H0: β will be tested. G×E The p-value obtained when p=0 is used to test the significance of ancestry-specific marginal gene-environment interactions at the K ancestral populations, i.e., to test the hypothesis. The obtained K statistical p-values ​​are combined using the Cauchy combination (CCT) method to obtain the final statistical p-value, so as to realize the analysis of gene-environment interactions across the entire genome.

[0114] After obtaining the final statistical p-value, a genome-wide gene-environment interaction analysis was performed at a given significance level.

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

[0116] This method is suitable for analyzing research cohorts including individuals from mixed populations or multiple ethnic groups. For genome-wide gene-environment interaction association analysis, it can control the Type I error rate to within a given significance level, exhibiting good statistical power and avoiding excessively high false positive or false negative rates due to population contamination. Specifically, this method constructs a regression model based on sample data, fits a constrained model where both marginal genetic effects and marginal gene-environment interactions are zero, and calculates the corresponding residuals. For the genetic locus to be tested, the locus genetic variation rate at the individual level is estimated based on the principal components of the genetic data. Based on the score statistic, a normal distribution approximation method is used to calculate the p-value to test the significance of the locus's marginal genetic effect. Based on the significance of the marginal genetic effect, a test statistic for testing the significance of the marginal gene-environment interaction at the locus is selected and calculated. For the test statistic testing the marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of the marginal gene-environment interaction. Based on ancestry... The specific score statistic uses a normal distribution approximation method to calculate the p-value to test the significance of the ancestry-specific marginal genetic effect at the locus. Based on the significance of the ancestry-specific marginal genetic effect, a test statistic for testing the significance of the ancestry-specific marginal gene-environment interaction at the locus is selected and calculated. For the test statistic testing the ancestry-specific marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of the ancestry-specific gene-environment interaction. The Cauchy combination method is used to combine the statistical p-value testing the significance of the gene-environment interaction with the statistical p-value testing the significance of the ancestry-specific gene-environment interaction across all ancestry groups to obtain the final statistical p-value, thereby achieving genome-wide gene-environment interaction analysis. This invention has advantages such as wide applicability, fast analysis speed, and high accuracy in gene-environment interaction analysis for mixed populations or multi-ethnic groups, avoiding excessively high false positive or false negative rates in the analysis results due to population mixing or stratified population structure. Attached Figure Description

[0117] The present invention will now be described in further detail with reference to the accompanying drawings and specific implementation methods.

[0118] Figure 1 This is a flowchart illustrating the gene-environment interaction analysis method for correctable population structure stratification according to the present invention.

[0119] Figure 2 This is a schematic diagram of the hybrid testing strategy that combines the normal distribution approximation method and the saddle point approximation method of the present invention;

[0120] Figure 3 This is a numerical simulation result of the Type I error rate in the heterogeneous population analysis of survival data phenotypes according to the present invention. Detailed Implementation

[0121] The present invention will now be described in further detail with reference to the accompanying drawings.

[0122] To further illustrate the features of the present invention, the present invention will be described in detail below with reference to the accompanying drawings and specific embodiments;

[0123] like Figure 3 As shown, the horizontal axis represents the grouping of genetic variation rates at loci, the vertical axis represents the empirical Type I error rate, ER: Event Rate, and MAF: Minor Allele Frequency, representing the genetic variation rate. The difference in genetic variation rate (MAF) between European (EUR) and East Asian (EAS) populations is represented by the Difference. MAF =q EUR -q EAS The genetic variation sites were divided into five groups: Diff MAF <<0(Diff MAF <-0.05), Diff MAF <0(-0.05≤Diff MAF <-0.01), Diff MAF ~0(-0.01≤Diff MAF ≤0.01), Diff MAF >0(0.01 <Diff MAF ≤0.05), and Diff MAF >>0(Diff MAF >0.05); based on the minimum genetic variation rate (MAF) in European (EUR) and East Asian (EAS) populations, i.e., min(q EUR ,q EAS The genetic variation sites are divided into three groups: minMAF low (min(q EUR ,q EAS)≤0.01),minMAF mod (0.01 <min(q EUR ,q EAS )≤0.05),minMAF high (min(q EUR ,q EAS ()>0.05); Based on the above two grouping indicators, all loci were divided into 15 (5×3) groups; Considering the event incidence rates of three pairs of European and East Asian populations: (ρ EUR ,ρ EAS ) = (0.1, 0.01) (low event incidence, ER) low ), (ρ EUR ,ρ EAS ) = (0.3, 0.05) (event occurrence rate, ER) mod ) and (ρ EUR ,ρ EAS ) = (0.5, 0.2) (High event incidence, ER high For each pair of genetic variation rate (MAF) and event occurrence rate (ER), a total of 10 were performed. 7 This is the second test for the significance of the locus's genetic effect, with a significance level of α = 5 × 10⁻⁶. -6 ;

[0124] like Figure 1 As shown, this embodiment discloses a gene-environment interaction analysis method that can correct for population structure stratification, and applies it to genome-wide association analysis of survival data phenotypes in mixed populations, including the following steps S1 to S10:

[0125] Step S1: Obtain phenotypic data, genotypic data, environmental factor data, and confounding factor data; confounding factor data includes age, sex, and principal genetic components, etc.

[0126] Step S2: Construct a regression model based on the sample data, fit a constrained model where both marginal genetic effect and marginal gene-environment interaction are 0, and calculate the corresponding residuals;

[0127] Step S3: For the genetic locus to be tested, estimate the locus genetic variation rate at the individual level based on the principal components of the genetic data;

[0128] Step S4: Calculate the p-value using the normal distribution approximation method based on the score statistic to test the significance of the marginal genetic effect at the locus;

[0129] Step S5: Select and calculate the test statistic for the significance of marginal gene-environment interactions at the loci based on the significance of the marginal genetic effect;

[0130] Step S6: For the test statistic of marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of marginal gene-environment interaction.

[0131] Step S7: Calculate the p-value using the normal distribution approximation method based on the ancestry-specific score statistic to test the significance of the ancestry-specific marginal genetic effect of the locus;

[0132] Step S8: Select and calculate the test statistic for the significance of ancestry-specific marginal genetic effects for testing loci based on the significance of the marginal genetic effects.

[0133] Step S9: For the test statistic of the marginal gene-environment interaction of ancestry specificity, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of the gene-environment interaction of ancestry specificity.

[0134] Step S10: Using the Cauchy combination method, the statistical p-value for testing the significance of gene-environment interactions and the statistical p-value for testing the significance of ancestry-specific gene-environment interactions across all ancestry types are combined to obtain the final statistical p-value, thereby achieving the analysis of gene-environment interactions across the entire genome.

[0135] Specifically, in step S1 above, it is assumed that the number of individual samples participating in the study is n, and the n individuals can come from different ethnic groups or subgroups. The individuals are not related by blood. For individual i, 1≤i≤n, let X i Let G represent a k-dimensional confounding vector. Confounding factors may include age, sex, and principal genetic components that contain information about the population structure of the sample. i Indicates the genotype of the locus, with values ​​of 0, 1, or 2, E i Y represents environmental factors. i It represents a certain phenotype.

[0136] Specifically, in step S2 above, for complex traits or phenotypes with different data formats, corresponding regression models are used for analysis (for example, using a linear regression model to analyze continuous phenotypes; using a Cox proportional hazards regression model to analyze survival data phenotypes; using a logistic regression model to analyze binary phenotypes; using a proportional dominance logistic regression model to analyze multi-class phenotypes; and using a mixed-effects multi-location scale model to analyze longitudinal data phenotypes; the linear predictor term in the regression model used is...). Where, β X Indicates confounding factor X i The m-dimensional coefficient vector, β E βG and β G×E These represent environmental factors E. i Genotype G i and gene-environment interaction term G i E i The corresponding coefficient.

[0137] To test the significance of the marginal gene-environment interaction at the locus, i.e., to test the null hypothesis H0: β G×E =0, under the assumption H c :β G =β G×E =0, meaning the marginal genetic effect and marginal gene-environment interaction at the locus are both zero. The model is fitted under this constraint to obtain the maximum likelihood estimate of the model parameters under this constraint, and the n-dimensional residual vector R = (R1,…,R) under the constraint model is calculated. n ) T .

[0138] The following three examples illustrate the process of constructing regression models based on sample data, fitting constrained models where the marginal genetic effects and marginal gene-environment interactions at loci are both zero, and calculating residuals in general genome-wide association analyses of confounding populations with continuous, binary, and survival data phenotypes, respectively:

[0139] Continuous phenotype:

[0140] Assuming the phenotypes are continuous data and the number of individuals in the sample is n, a linear regression model is used for analysis:

[0141]

[0142] Among them, Y i X represents a continuous phenotype. i Let E represent a k-dimensional confounding factor vector. i G represents environmental factors. i Indicates the genotype of the locus, ε i Let β represent the random error term that follows a normal distribution. X Indicates confounding factor X i a k-dimensional coefficient vector, with coefficient β E Indicates environmental factor E i The effect, coefficient β G Indicates the locus genotype G i The marginal genetic effect, coefficient β G×E This indicates the interaction between marginal genes and the environment.

[0143] To test the significance of the marginal gene-environment interaction at the locus, under the hypothesis H c :β G =βG×E =0, meaning the fitting model is under the constraint that both the marginal genetic effect of the locus and the marginal gene-environment interaction are 0:

[0144]

[0145] Obtain model parameters (β) X ,β E Under constraint β G =β G×E Maximum likelihood estimation at 0 And calculate Used to represent the linear prediction term η i The estimation under the constrained model, let Indicates in β G =β G×E The residuals obtained after fitting the constrained model under the constraint condition of 0 are given by R = (R1, ..., R2). n ) T This represents the corresponding residual vector.

[0146] Binary phenotype:

[0147] Assuming the phenotype is binary and the number of individuals in the sample is n, a logistic regression model is used for analysis:

[0148]

[0149] Among them, Y i Indicates a binary phenotype (values ​​0 or 1), X i Let E represent a k-dimensional confounding factor vector. i G represents environmental factors. i Indicates the genotype of the locus, ε i Let β represent the random error term that follows a normal distribution. X Indicates confounding factor X i a k-dimensional coefficient vector, with coefficient β E Indicates environmental factor E i The effect, coefficient β G Indicates genotype G i The marginal genetic effect, coefficient β G×E This indicates the marginal gene-environment interaction. Indicates that given X i E i and G i Under the condition of Y i The probability that = 1.

[0150] To test the significance of the marginal gene-environment interaction at the locus, under the hypothesis H c :β G=β G×E =0, meaning the fitting model is under the constraint that both the marginal genetic effect of the locus and the marginal gene-environment interaction are 0:

[0151]

[0152] Obtain the model parameters (β) X ,β E Under constraint β G =β G×E Maximum likelihood estimation at 0 And calculate Used to represent the linear prediction term η i The estimation under the constrained model, let Representing π i The estimation under the constrained model, let Indicates in β G =β G×E The residuals obtained after fitting the constrained model under the constraint condition of 0 are given by R = (R1, ..., R2). n ) T This represents the corresponding residual vector.

[0153] Survival data typology:

[0154] Assuming the phenotypes are survival data and the number of individuals in the sample is n, the following Cox proportional hazards regression model is used for analysis:

[0155]

[0156] Wherein, λ(t; X i E i G i Let λ0(t) represent the risk function, λ0(t) represent the baseline risk function, and X0(t) represent the risk function. i Let E represent a k-dimensional confounding factor vector. i G represents environmental factors. i Indicates the genotype of the locus, β X Indicates confounding factor X i a k-dimensional coefficient vector, with coefficient β E Indicates environmental factor E i The effect, coefficient β G Indicates genotype G i The marginal genetic effect, coefficient β G×E This indicates the marginal gene-environment interaction, causing C i Indicates the time the event was deleted. Indicates the time when the event occurred. This represents the observed survival data. Used to indicate whether an event has occurred, where the function I(.) represents the indicative function.

[0157] To test the significance of the marginal gene-environment interaction at the locus, under the hypothesis H c :β G =β G×E =0, meaning the fitting model is under the constraint that both the marginal genetic effect of the locus and the marginal gene-environment interaction are 0:

[0158]

[0159] Obtain model parameters (β) X ,β E Under constraint β G =β G×E Maximum likelihood estimation at 0 And calculate Used to represent the linear prediction term η i The estimation under the constrained model, let Indicates at time t i The risk set, making Represents individual i up to time t. i The cumulative risks make Represents Λ i In β G =β G×E The estimate obtained after fitting the constrained model under the constraint condition of 0 is let Indicates in β G =β G×E The martingale residuals calculated after fitting the constrained model under the constraint condition of 0 are given by R = (R1, ..., R). n ) T This represents the corresponding residual vector.

[0160] Specifically, in step S3 above, for the locus to be tested, it is assumed that the individual-level genetic variation rate vector of n individuals is q = (q1, ..., q n ) T , where q i q represents the genetic variation rate of individual i at this locus. Since n individuals can come from different populations, and the genetic variation rate of a locus is not necessarily the same in different populations, the elements of q can be unequal, and the genetic variation rate is individual-specific.

[0161] Let G = (G1, ..., G n ) T The genotype vector representing a locus, based on the Hardy-Weinberg law in genetics, considers G... i It follows a binomial distribution Binom(2,q) iGiven a random variable G, since it is assumed that the n individuals are not related, the elements in G are independent of each other. The following describes how to estimate the individual-level genetic variation rate of the n individuals.

[0162] If n individuals come from the same population, then they can be considered to have the same locus genetic variation rate, i.e., q1 = q2 = ... = q n At this point, each element in the genetic variation rate vector q of the locus can be derived from... Make an estimate.

[0163] If n individuals come from multiple ethnic groups, mixed populations, or multi-ethnic groups, the locus genetic variation rate q of different individuals is... i It is highly likely that they will be different. In genetic research, principal component analysis is often used to obtain the genetic principal components that contain population structure information of the sample. The genetic principal components that contain the structure information of the entire population can be used to accurately estimate the genetic variation rate at the individual level.

[0164] make It represents an n×(D+1) dimensional genetic principal component matrix consisting of column vectors containing all elements equal to 1 and the first D genetic principal component vectors containing all population structure information of the sample individuals.

[0165] Based on the genetic principal component, the response variable and explanatory variable can be fitted to G and the genetic principal component X, respectively. PC The linear regression model calculates an estimate of q. However, for sites with low genetic variation rates... It may deviate significantly from the true value. Therefore, let q c =0, can Based on this, another estimate of q is obtained. in, I(.) denotes the characteristic function. All elements are between 0 and 1.

[0166] Based on genetic principal components, the response variable and explanatory variables can be fitted as (I(G1≥0.5),…,I(G1≥0.5), respectively. n ≥0.5)) T and principal component X PC The logistic regression model is used to obtain an estimate of q. in, The function σ(x) = 1 / (1 + exp(-x)) This represents the linear predictor term in the logistic regression model. The maximum likelihood estimate.

[0167] Compared to the estimates obtained using a linear regression model Estimation obtained using logistic regression model The estimation of the individual-level genetic variation rate vector q is more accurate, but the computational cost is higher.

[0168] To balance the accuracy and computational efficiency of estimating the individual-level genetic variation rate vector q, a hybrid estimation strategy utilizing linear regression and logistic regression models is proposed. Specifically: If If the proportion of elements in the interval [0,1] exceeds a pre-defined positive number between 0 and 1 (default value 0.9), then use... As an estimate of q; otherwise, let As an estimate of q.

[0169] Specifically, in step S4 above, the score test statistic used to test the marginal genetic effect is calculated.

[0170]

[0171] And it is believed that the test statistic Assuming H G :β G When = 0 is true, it follows the expected value. variance is The normal distribution of, let As a method for testing hypothesis H using the normal distribution approximation G :β G =0 yields the two-sided p-values, where for The observed values ​​are Φ(.), which represents the cumulative distribution function of the standard normal distribution.

[0172] Specifically, in step S5 above, testing the significance of the marginal gene-environment interaction is equivalent to testing the null hypothesis H0: β G×E Is the expression = 0 true?

[0173] Let ∈ be a given positive number between 0 and 1.

[0174] If we consider hypothesis H G :β G If the p-value obtained by testing = 0 is less than or equal to ∈, then use As a test of the null hypothesis H0:β G×E The test statistic for = 0, where

[0175] If we consider hypothesis H G :β GIf the p-value obtained by testing =0 is greater than ∈, then the Wald test (or likelihood ratio test) and the test statistic are used respectively. Regarding the null hypothesis H0:β G×E =0 is used for testing, where To adapt to the residual vector of genotype, I n Let W be an n-order identity matrix, then W = (1 n G) is an n×2 dimensional matrix, 1 n =(1,1,…,1) T For an n-dimensional column vector containing only 1s, we will use The two p-values ​​obtained from the Wald test (or likelihood ratio test) are combined using the Cauchy combination (CCT) method to obtain the final p-value.

[0176] Specifically, in step S6 above, it is assumed that the sample genotype vector of this locus is G = (G1, ..., G...). n ) T Fitting β G =β G×E The residual vector obtained from the constraint model with 0 is R = (R1, ..., R2). n ) T .

[0177] Under the Hardy-Weinberg equilibrium hypothesis in genetics, the locus genotype G of each individual is... i ,1≤i≤n are considered to be mutually independent and follow a binomial distribution Binom(2,q) i A random variable q i This represents the allele frequency of the genetic variation of the i-th individual at that locus. A hybrid test strategy combining the normal distribution approximation and the saddle point approximation is used to approximate the test statistic in relation to the null hypothesis H0: β. G×E The conditional distribution when =0 is true will be explained in detail below.

[0178] When using the normal distribution approximation method to calculate the statistical p-value, the specific steps are as follows:

[0179] If the test statistic is S G×E Given (R, E, λ), the normal distribution approximation method uses an expectation of variance is The normal distribution of the null hypothesis H0:β G×E When = 0 holds true, the statistic S G×E The distribution is approximated, and the two-sided p-value is calculated using the normal distribution approximation method. Where E = (E1, ..., E n ) Ts represents the environmental factor vector. G×E For S G×E The observed values ​​are Φ(.), which represents the cumulative distribution function of the standard normal distribution.

[0180] If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation and variance as follows: and The normal distribution of the null hypothesis H0:β G×E =0 statistic The distribution of is approximated, and the two-sided p-value obtained using the normal distribution approximation method is: in for The observed values.

[0181] When using the saddle point approximation method to calculate the statistical p-value, the specific steps are as follows:

[0182] The saddle point approximation method approximates the distribution by utilizing all higher-order moments contained in the cumulant generating function. Because it incorporates more information about the distribution, its approximation accuracy is higher than that of the normal distribution approximation method, which only uses the first two moments. In genome-wide association studies of confounded populations, to obtain accurate statistical inferences when testing the significance of locus genetic effects, the saddle point approximation method is introduced below to approximate the test statistic at the null hypothesis H0:β. G×E The specific process of calculating the p-value from the conditional distribution when =0 is true.

[0183] If the test statistic is S G×E Then, under the null hypothesis H0:β G×E When = 0 is true, G i The i = 1, ..., n are considered to be independent and follow a binomial distribution Binom(2, q i Given a sequence of random variables (R, E, λ), in order to estimate S... G×E Given the cumulant generating function when the null hypothesis H0 holds, first estimate G. i Moment generating function:

[0184]

[0185] The first and second derivatives are as follows:

[0186]

[0187] G i The estimate of the cumulant generating function is: Its first and second derivatives are as follows:

[0188]

[0189] Therefore, when H₀ holds, the G×E estimated cumulant generating function is:

[0190]

[0191] Its first-order derivative and second-order derivative are respectively:

[0192]

[0193] For any given real number s₀, first calculate the ζ that satisfies, then calculate and To approximate the distribution of S G×E under the null hypothesis H₀, according to the Barndorff-Nielsen saddlepoint approximation formula, is used to approximately estimate the probability Pr(S G×E <s₀|(R,E,λ)). For a given test statistic S G×E observed value s G×E , the two-sided p-value calculated by the above saddlepoint approximation algorithm is p l +p r , the left-sided p-value p l and the right-sided p-value p r are respectively:

[0194]

[0195] wherein and respectively represent the probabilities obtained by the above saddlepoint approximation algorithm and estimations.

[0196] If the test statistic is then when the null hypothesis H₀:β G×E = 0 holds, let G i , i = 1, …, n be regarded as a sequence of mutually independent random variables following the binomial distribution Binom(2,q i ), under given conditions, when the null hypothesis H₀ holds, the estimated cumulant generating function is:

[0197]

[0198] Its first-order derivative and second-order derivative are respectively:

[0199]

[0200] For any given real number First calculate to make Established Then the calculation yields... and use To approximate the probability For a given Observations The bilateral p-value calculated using the saddle point approximation algorithm described above is p. l +p r Let p be the value on the left side. l and the p-value on the right. r They are respectively:

[0201]

[0202] in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate.

[0203] The statistical p-value is calculated using a hybrid test strategy that combines the normal distribution approximation method with the saddle point approximation method, as detailed below:

[0204] Because the saddle point approximation method is less computationally efficient than the normal distribution approximation method when calculating the statistical p-value, in order to balance accuracy and computational efficiency, the test statistic S is used... G×E or S G×E When calculating the statistical p-value, the following method is used: Figure 2 The hybrid testing strategy shown combines the normal distribution approximation method with the saddle point approximation method.

[0205] Let r be a pre-selected positive number.

[0206] If we test the null hypothesis H0: β G×E The statistic for = 0 is S G×E s G×E For S G×E The observed values, if If the p-value is calculated using the normal distribution approximation method, then the p-value is calculated using the saddle point approximation method; otherwise, the p-value is calculated using the saddle point approximation method.

[0207] If we test the null hypothesis H0: β G×E The statistic for =0 is for The observed values, if If the p-value is calculated using the normal distribution approximation method, then the p-value is calculated using the saddle point approximation method; otherwise, the p-value is calculated using the saddle point approximation method.

[0208] Specifically, in step S7 above, assuming that the n sample individuals in the research cohort come from K ancestral populations, let G = (G1, ..., G... n ) T Genotype vector representing a locus. This represents the genotype vector specific to the kth ancestry (ancestral population).

[0209] For individual i, i≤n, let This represents the number of haplotypes in the k-th ancestry (ancestral population) at that locus, i.e., the local ancestry count. Let the vector representing the corresponding local lineage count be... This represents an estimate of the ancestry-specific allele frequencies of the k-th ancestry (ancestral population).

[0210] make Let represent the ancestry-specific marginal genetic effect of the k-th ancestry (ancestral population), and let Let the ancestry-specific marginal gene-environment interaction of the k-th ancestry (ancestral population) be represented. Then, the linear prediction term in the regression model can be rewritten as follows:

[0211] Testing the significance of the marginal genetic effect of ancestry specificity is equivalent to testing the hypothesis. Is it true?

[0212] Calculate the score test statistic for testing the marginal genetic effect specific to ancestry.

[0213]

[0214] And it is believed that the test statistic In the assumption When established, it obeys the expected value. variance is The normal distribution of, let The two-sided p-value obtained by using the normal distribution approximation method to test the significance of the marginal genetic effect of lineage specificity, where for The observed values ​​are Φ(.), which represents the cumulative distribution function of the standard normal distribution.

[0215] Specifically, in step S8 above, testing the significance of ancestry-specific marginal gene-environment interactions is equivalent to testing the hypothesis. Whether it is valid or not.

[0216] Let ∈ be a given positive number between 0 and 1.

[0217] If we consider the hypothesis If the p-value obtained from the test is less than or equal to ∈, then use As a hypothesis to be tested The test statistic, of which

[0218] If we consider the hypothesis If the p-value obtained from the test is greater than ∈, then the Wald test (or likelihood ratio test) and the test statistic are used respectively. For the hypothesis The test was conducted, among which To adapt to the genotype-specific residual vector, I n W is an n-order identity matrix. (k) =(1 n G (k) ) is an n×2 dimensional matrix, 1 n =(1,1,…,1) T For an n-dimensional column vector containing only 1s, we will use Wald test (or likelihood ratio test) hypothesis The two obtained p values ​​are combined using the Cauchy combination method to obtain the final p value.

[0219] Specifically, in step S9 above, the statistical p-value for testing the marginal gene-environment interaction of ancestry specificity is calculated using a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method based on the locus genetic variation rate at the individual level, in order to test the significance of the gene-environment interaction of ancestry specificity.

[0220] A hybrid testing strategy combining the normal distribution approximation method and the saddle point approximation method is used to approximate the test statistic under the hypothesis. The conditional distribution when true is explained in detail below.

[0221] When using the normal distribution approximation method to calculate the statistical p-value, the specific steps are as follows:

[0222] If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation of variance is The normal distribution is related to the hypothesis. Statistics at the time of establishment If the distribution is approximated, then the two-sided p-value calculated using the normal distribution approximation method is... in for The observed values ​​are Φ(.), which represents the cumulative distribution function of the standard normal distribution.

[0223] If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation and variance as follows: and The normal distribution is related to the hypothesis. Statistics at the time of establishment If the distribution is approximated, then the two-sided p-value obtained using the normal distribution approximation method is: in for The observed values.

[0224] When using the saddle point approximation method to calculate the statistical p-value, the specific steps are as follows:

[0225] If the test statistic is Then in the assumption When it was established Consider them as independent entities that follow a binomial distribution, Binom(2,q) (k) A sequence of random variables, in Given the conditions, in order to estimate In the assumption The cumulative generating function at the time of its establishment is first estimated. Moment generating function:

[0226]

[0227] The estimate of the cumulant generating function is: Its first and second derivatives are as follows:

[0228]

[0229] Therefore, in When it was established The estimate of the cumulant generating function is:

[0230]

[0231] Its first and second derivatives are as follows:

[0232]

[0233] For any given real number s0, first calculate such that The established ζ (k) Then calculate to get and To approximate Under the original hypothesis The distribution below, according to the Barndorff-Nielsen saddle point approximation formula, uses... To approximate the probability Then for a given test statistic Observations The bilateral p-value calculated using the saddle point approximation algorithm described above is p. l +p r Let p be the value on the left side. l and the p-value on the right. r They are respectively:

[0234]

[0235] in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate.

[0236] If the test statistic is Then in the assumption When it was established Consider them as independent entities that follow a binomial distribution, Binom(2,q) (k) A sequence of random variables, in Given the conditions, In the assumption The estimate of the cumulant generating function when it holds true is:

[0237]

[0238] Its first and second derivatives are as follows:

[0239]

[0240] For any given real number First calculate to make Established Then the calculation yields... and use To approximate the probability For a given Observations The bilateral p-value calculated using the saddle point approximation algorithm described above is p.l +p r Let p be the value on the left side. l and the p-value on the right. r They are respectively:

[0241]

[0242] in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate.

[0243] The statistical p-value is calculated using a hybrid test strategy that combines the normal distribution approximation method with the saddle point approximation method, as detailed below:

[0244] Because the saddle point approximation method is less computationally efficient than the normal distribution approximation method in calculating the statistical p-value, in order to balance accuracy and computational efficiency, the test statistic is used... or When calculating the statistical p-value, the following method is used: Figure 2 The hybrid testing strategy shown combines the normal distribution approximation method with the saddle point approximation method.

[0245] Let r be a pre-selected positive number.

[0246] If we test the null hypothesis The statistic is for The observed values, if If the p-value is calculated using the normal distribution approximation method, then the p-value is calculated using the saddle point approximation method; otherwise, the p-value is calculated using the saddle point approximation method.

[0247] If we test the null hypothesis The statistic is for The observed values, if If the p-value is calculated using the normal distribution approximation method, then the p-value is calculated using the saddle point approximation method; otherwise, the p-value is calculated using the saddle point approximation method.

[0248] Specifically, in step S10 above, the significance of the marginal gene-environment interaction at the locus will be tested, that is, the null hypothesis H0:β will be tested. G×E The p-value obtained when p=0 is used to test the significance of ancestry-specific marginal gene-environment interactions at the K ancestral populations, i.e., to test the hypothesis. The obtained K statistical p-values ​​are combined using the Cauchy combination (CCT) method to obtain the final statistical p-value.

[0249] After obtaining the final statistical p-value, a genome-wide gene-environment interaction analysis was performed at a given significance level.

[0250] According to steps S1 to S10, this gene-environment interaction analysis method that can correct for population structure stratification, when applied, constructs a regression model based on sample data, fits a constrained model where both marginal genetic effects and marginal gene-environment interactions are zero, and calculates the corresponding residuals; for the genetic locus to be tested, the locus genetic variation rate at the individual level is estimated based on the principal components of the genetic model; the p-value is calculated using the normal distribution approximation method based on the score statistic to test the significance of the marginal genetic effect at the locus; based on the significance of the marginal genetic effect, a test statistic for testing the significance of the marginal gene-environment interaction at the locus is selected and calculated; for the test statistic for testing the marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of the marginal gene-environment interaction; based on blood... The p-value of the pedigree-specific score was calculated using the normal distribution approximation method to test the significance of the pedigree-specific marginal genetic effect. Based on the significance of the pedigree-specific marginal genetic effect, a test statistic for testing the significance of the pedigree-specific marginal gene-environment interaction was selected and calculated. For the test statistic testing the pedigree-specific marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method was used to calculate the statistical p-value to test the significance of the pedigree-specific gene-environment interaction. The Cauchy combination method was used to combine the statistical p-value testing the significance of the gene-environment interaction with the statistical p-value testing the significance of the pedigree-specific gene-environment interaction across all pedigrees to obtain the final statistical p-value, thus achieving the analysis of gene-environment interactions across the entire genome. Therefore, this gene-environment interaction analysis method that can correct for population structure stratification has the advantages of wide applicability, fast analysis speed and high accuracy in large-scale genome-wide association analysis. It can avoid a large number of false negative or false positive results, and is suitable for association analysis of gene-environment interactions in samples from mixed populations or multiple ethnic groups. It is also suitable for the analysis of various data phenotypes (including but not limited to continuous phenotypes, binary phenotypes, survival data phenotypes, multi-category phenotypes and longitudinal data phenotypes and other complex phenotypes).

[0251] The following is a specific numerical simulation example:

[0252] This example provides a numerical simulation example to verify the effectiveness of the method proposed in this paper. In this example, the proposed algorithm is applied to gene-environment interaction analysis for survival data phenotypes in admixed populations, and numerical simulations are performed to evaluate the type I error rate of the proposed method.

[0253] The specific implementation steps are as follows:

[0254] (1) Generation of sample population structure

[0255] In the numerical simulation, the individual sample size is n=10000, and it is assumed that n individuals come from an admixed population of European (EUR) and East Asian (EAS) populations, and there is no genetic correlation between individuals. For the i-th sample individual, let denote its ancestry component vector, where and represent the ancestry proportions from European (EUR) ancestry and East Asian (EAS) ancestry respectively, and satisfy It is assumed that the first n / 2=5000 individuals come from a population dominated by European ancestry, and the ancestry component vector is generated by the Dirichlet distribution Dirichlet(9,1); the latter n / 2=5000 individuals come from a population dominated by East Asian ancestry, and the ancestry component vector is generated by the Dirichlet distribution Dirichlet(1,9).

[0256] (2) Generation of minor allele frequencies and genetic variant loci

[0257] Based on the difference in minor allele frequency (MAF) of loci between European population (EUR) and East Asian population (EAS) in the 1000 Genomes Project, that is DiffMAF=q EUR -q EAS , genetic variant loci are divided into five groups: DiffMAF<<0(DiffMAF<-0.05),DiffMAF<0(-0.05≤DiffMAF<-0.01),DiffMAF~0(-0.01≤DiffMAF≤0.01),DiffMAF>0(0.01<DiffMAF≤0.05), and DiffMAF>>0(DiffMAF>0.05); based on the minimum value of minor allele frequency (MAF) in European population (EUR) and East Asian population (EAS), that is min(q EUR ,q EAS ), genetic variant loci are divided into three groups: minMAFlow(min(q EUR ,q EAS)≤0.01),minMAFmod(0.01 <min(q EUR ,q EAS )≤0.05),minMAFhigh(min(q EUR ,q EAS Based on the above two grouping indicators, all genetic variation loci were divided into 15 (5×3) groups. Within each group, the genetic variation rate (q) of 1000 randomly selected loci was analyzed. EUR ,q EAS ), and generate 1000 corresponding genetic variation loci based on 1000 pairs of genetic variation rates. Specifically, for each locus, individual i has a locus genotype G. i By binomial distribution Binom(2,q) i ) generated, where In addition, 100,000 mutually independent common genetic variation sites (q) were generated. EUR +q EAS >0.1) and use it to generate the principal components of the genetic sequence of n individuals.

[0258] (III) Generation of Survival Data Tables

[0259] The survival data phylogenetics for n individuals are generated through two steps. For the i-th sample individual, the occurrence time of the event is first generated. and deletion time C i Then, the survival data phylogenetic model was calculated. and indication quantity The function I(.) is the indicator function. Censoring time C i Generated by a Weibull distribution with a scale parameter of 0.15 and a shape parameter of 1. Event occurrence time. Generated from a Cox proportional hazards model with a Weibull baseline risk function:

[0260]

[0261] Among them, U i Generated by a uniform distribution U(0,1), the linear prediction term η i =β1X i1 +0.5X i2 +0.5X i3 +0.5E i +β G G i +β G×E G i E i Confounding factors Indicates the proportion of East Asian ancestry, with confounding factor X.i2 Generated by Bernoulli distribution (0.5), with confounding factor X. i3 Generated by a standard normal distribution, environmental factor E i Generated by a uniform distribution U(0,1), G i The genotype of the locus is represented by the parameter β. G The genetic effect of a locus is represented by the parameter β. G×E This represents gene-environment interactions. Appropriate selection of parameters λ and β1 is needed to obtain the expected event incidence rates (ER) in European and East Asian populations. EUR and ER EAS .

[0262] In this embodiment, the numerical simulation of the Type I error rate considers the event occurrence rates of three pairs of European and East Asian populations: (ER EUR ,ER EAS = (0.1, 0.01) (Low event rate, ERlow), (ER EUR ,ER EAS ) = (0.3, 0.05) (event occurrence rate, ERmod) and (ER EUR ,ER EAS = (0.5, 0.2) (High event rate, ERhigh).

[0263] (iv) Numerical simulation of Type I error rate in survival data phenotypic analysis:

[0264] To evaluate the Type I error rate, in β G×E =β G Survival phenotypes were generated under a constraint model with β=0. For each case, 1000 datasets containing phenotypes and confounding factors were generated. Therefore, for each pair of genetic variation rate (MAF) and event occurrence rate (ER), a total of 1,000,000 tests were performed to assess the significance of marginal gene-environment interactions at each locus. For all methods, the genetic effect β was fitted... G×E =β G When the constraint model is 0, the covariates in the model include X. i2 X i3 Environmental factors E i And the first four principal genetic components.

[0265] At a significance level of α = 5 × 10 -5 Below this, the marginal genetic effect and marginal gene-environment interaction at the locus are both 0, i.e., β G×E =β G In the case where = 0, evaluate and compare the empirical Type I error rates of the proposed new algorithm with those of the normal distribution approximation method.

[0266] (V) Numerical simulation results of the Type I error rate in survival data phenotypic analysis:

[0267] Figure 3 This demonstrates the significance level at α = 5 × 10⁻⁶. -5 The new algorithm effectively controls the Type I error rate in all cases, based on empirical data. The normal distribution approximation method cannot accurately control the Type I error rate, and false positives may occur in the analysis results for cases with low genetic variation rates.

[0268] The numerical simulation results above show that the gene-environment interaction analysis method with correctable population structure stratification proposed in this invention can effectively control the Type I error rate when analyzing samples from mixed populations or multi-ethnic groups, and can avoid excessive false positive or false negative results caused by population mixing or population structure stratification.

[0269] This document uses specific examples to illustrate the principles and implementation methods of the present invention. The descriptions of the above embodiments are only for the purpose of helping to understand the method and core ideas of the present invention. Furthermore, those skilled in the art will recognize that, based on the ideas of the present invention, there will be changes in the specific implementation methods and application scope. Therefore, the content of this specification should not be construed as a limitation of the present invention.

Claims

1. A gene-environment interaction analysis method that can correct for population structure stratification, applicable to genome-wide association analysis of gene-environment interactions in complex phenotypes of mixed populations, characterized by: The method includes the following steps: Step S1: Obtain phenotypic data, genotypic data, environmental factor data, and confounding factor data; confounding factor data includes age, sex, and principal genetic components; Step S2: Construct a regression model based on the sample data, fit a constrained model where both marginal genetic effect and marginal gene-environment interaction are 0, and calculate the corresponding residuals; Step S3: For the genetic locus to be tested, estimate the locus genetic variation rate at the individual level based on the principal components of the genetic data; Step S4: Calculate the p-value using the normal distribution approximation method based on the score statistic to test the significance of the marginal genetic effect at the locus; Step S5: Select and calculate the test statistic for the significance of marginal gene-environment interactions at the loci based on the significance of the marginal genetic effect; Step S6: For the test statistic of marginal gene-environment interaction, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of marginal gene-environment interaction. Step S7: Calculate the p-value using the normal distribution approximation method based on the ancestry-specific score statistic to test the significance of the ancestry-specific marginal genetic effect of the locus; Step S8: Select and calculate the test statistic for the significance of ancestry-specific marginal genetic effects for testing loci based on the significance of the marginal genetic effects. Step S9: For the test statistic of the marginal gene-environment interaction of ancestry specificity, a hybrid test strategy combining the normal distribution approximation method and the saddle point approximation method is used to calculate the statistical p-value to test the significance of the gene-environment interaction of ancestry specificity. Step S10: Using the Cauchy combination method, the statistical p-value for testing the significance of gene-environment interactions and the statistical p-value for testing the significance of ancestry-specific gene-environment interactions across all ancestry types are combined to obtain the final statistical p-value, thereby achieving the analysis of gene-environment interactions across the entire genome.

2. The gene-environment interaction analysis method for correcting population structure stratification according to claim 1, characterized in that: In step S1, let the number of sample individuals participating in the study be... , Individuals can come from different ethnic groups or subgroups, and there is no kinship between individuals. ,make Represent a The confounding factor vector contains age, sex, and principal genetic components that contain information about the population structure of the sample. Represents the genotype of the locus, with values ​​ranging from 0 to 1. , Indicates environmental factors, It represents a certain phenotype.

3. The gene-environment interaction analysis method for correcting population structure stratification according to claim 2, characterized in that: In step S2, for complex traits or phenotypes with different data formats, corresponding regression models are used for analysis. The linear prediction term in the regression model is... , in, Indicating confounding factors of dimensional coefficient vector, , and These represent environmental factors. ,genotype and gene-environment interactions The corresponding coefficients; in order to test the significance of the marginal gene-environment interaction at the locus, i.e., to test the null hypothesis. , in the assumption This involves fitting a model under the constraint that both the marginal genetic effect and the marginal gene-environment interaction at a locus are zero, obtaining the maximum likelihood estimate of the model parameters under this constraint, and calculating the model under the constraint. 3D residual vector .

4. The gene-environment interaction analysis method for correcting population structure stratification according to claim 3, characterized in that: In step S3, let The vector of individual-level genetic variation rate for each individual is: ,in Represents an individual The rate of genetic variation at that locus; make The genotype vector representing the locus, let Represents a column vector containing all elements equal to 1 and a front vector containing all population structure information. Genetic principal component vectors Dimensional genetic principal component matrix; Based on genetic principal components, the response variable and explanatory variables were fitted respectively. and principal components of genetics The linear regression model was calculated to obtain An estimate ; Therefore, we obtain An estimate ,in, , , Indicates characteristic functions; Based on genetic principal components, the response variable and explanatory variables were fitted respectively. and principal components of genetics The logistic regression model can be obtained An estimate ,in, ,function , This represents the linear predictor term in the logistic regression model. Maximum likelihood estimation; To balance the accuracy of the estimation with computational efficiency, if If the proportion of elements in the interval [0,1] exceeds 0.9, then let As The estimate, otherwise, let As The estimate.

5. The gene-environment interaction analysis method for correcting population structure stratification according to claim 4, characterized in that: In step S4, the score test statistic used to test the marginal genetic effect is calculated: And it is believed that the test statistic In the assumption When established, it obeys the expected value. variance is The normal distribution of, let As a method for testing hypotheses using the normal distribution approximation. The obtained two-sided p-values, where for The observed values, The cumulative distribution function represents the standard normal distribution.

6. The gene-environment interaction analysis method for correcting population structure stratification according to claim 5, characterized in that: In step S5, the significance of the marginal gene-environment interaction is tested, i.e., the null hypothesis is tested. Is it valid? make Given a positive number between 0 and 1; If we consider the hypothesis The p-value obtained from the test is less than or equal to Then use As a test of the null hypothesis The test statistic, of which ; If we consider the hypothesis The p-value obtained from the test is greater than Then use the Wald test and the test statistic respectively. Regarding the null hypothesis The test was conducted, among which To adapt to the residual vector of genotype, for An identity matrix of order 1. for 3D matrix For elements all equal to 1 3D column vectors will use The two p-values ​​obtained from the Wald test are combined using the Cauchy combination method to obtain the final p-value.

7. The gene-environment interaction analysis method for correcting population structure stratification according to claim 6, characterized in that: In step S6, it is assumed that the sample genotype vector of this site is , fitting The residual vector obtained from the constraint model is ; Under the Hardy-Weinberg equilibrium hypothesis in genetics, the locus genotype of each individual is... Treating them as independent entities that follow a binomial distribution a random variable, where Indicates the first The rate of genetic variation at this locus for each individual sample; Based on the locus genetic variation rate at the individual level, an approximate test statistic is used to verify the null hypothesis. The statistical p-value is calculated using the conditional distribution when the condition is true. When using the normal distribution approximation method to calculate the statistical p-value, the specific steps are as follows: If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation of The variance is The normal distribution of the null hypothesis Statistics at the time of establishment The distribution is approximated, and the two-sided p-value is calculated using the normal distribution approximation method. ,in Represents an environmental factor vector. for The observed values, The cumulative distribution function representing the standard normal distribution; If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation and variance as follows: and The normal distribution of the null hypothesis Statistics at the time of establishment The distribution of is approximated, and the two-sided p-value obtained using the normal distribution approximation method is: ,in for Observed values; When using the saddle point approximation method to calculate the statistical p-value, the specific steps are as follows: If the test statistic is Then, under the null hypothesis When it was established They are considered to be independent and follow a binomial distribution. A sequence of random variables, in Given the conditions, in order to estimate Under the original hypothesis The cumulative generating function at the time of its establishment is first estimated. Moment generating function: The first and second derivatives are as follows: The estimate of the cumulant generating function is: Its first and second derivatives are as follows: Therefore, in When it was established The estimate of the cumulant generating function is: Its first and second derivatives are as follows: For any given real number First calculate to make Established Then calculate to get and In order to approximate Under the original hypothesis The distribution below, according to the Barndorff-Nielsen saddle point approximation formula, uses... To approximate the probability For a given test statistic Observations The bilateral p-values ​​calculated using the saddle point approximation algorithm described above are... Record the p-value on the left side. and the p-value on the right They are respectively: in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate; If the test statistic is Then, under the null hypothesis When it was established They are considered to be independent and follow a binomial distribution. A sequence of random variables, in Given the conditions, Under the original hypothesis The estimate of the cumulant generating function when it holds true is: Its first and second derivatives are as follows: For any given real number First calculate to make Established Then calculate to get and ;use To approximate the probability For a given Observations The bilateral p-values ​​calculated using the saddle point approximation algorithm described above are... Record the p-value on the left side. and the p-value on the right They are respectively: in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate; The statistical p-value is calculated using a hybrid test strategy that combines the normal distribution approximation method with the saddle point approximation method, as detailed below: make It is a pre-selected positive number; If we test the null hypothesis The statistic is , for The observed values, if If the p-value is positive, the normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value. If we test the null hypothesis The statistic is , for The observed values, if If the p-value is positive, the normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value.

8. The gene-environment interaction analysis method for correcting population structure stratification according to claim 7, characterized in that: In step S7, it is assumed that the research queue contains The sample individuals came from A bloodline, making Genotype vector representing a locus; Indicates the first A lineage-specific genotype vector; For individuals ,make Indicates the number at that site The number of haplotypes in a pedigree, i.e., the local pedigree count. Let the vector representing the corresponding local lineage count be... Indicates the first An estimate of the ancestry-specific allele frequencies of a ancestry. make Indicates the first The ancestry-specific marginal genetic effect of an individual bloodline makes Indicates the first For ancestry-specific marginal gene-environment interactions, the linear predictor term in the regression model can be rewritten as follows: ; Testing the significance of the marginal genetic effect of ancestry specificity is equivalent to testing the hypothesis. Is it true? Calculate the score test statistic for testing the marginal genetic effect specific to ancestry. And it is believed that the test statistic In the assumption When established, it obeys the expected value. variance is The normal distribution of, let The two-sided p-value obtained by using the normal distribution approximation method to test the significance of the marginal genetic effect of lineage specificity, where for The observed values, The cumulative distribution function represents the standard normal distribution.

9. The gene-environment interaction analysis method for correcting population structure stratification according to claim 8, characterized in that: In step S8, the significance of the ancestry-specific marginal gene-environment interaction is tested, i.e., the hypothesis is tested. Is it true? make Given a positive number between 0 and 1; If we consider the hypothesis The p-value obtained from the test is less than or equal to Then use As a hypothesis to be tested The test statistic, of which ; If we consider the hypothesis The p-value obtained from the test is greater than Then, use the Wald test (or likelihood ratio test) and the test statistic respectively. For the hypothesis The test was conducted, among which To adapt to the genotype-specific residual vector, for An identity matrix of order 1. for 3D matrix For elements all equal to 1 3D column vectors will use Wald test hypothesis The two obtained p values ​​are combined using the Cauchy combination method to obtain the final p value.

10. The gene-environment interaction analysis method for correcting population structure stratification according to claim 9, characterized in that: In step S9, when calculating the statistical p-value using the normal distribution approximation method, the specific steps are as follows: If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation of The variance is The normal distribution is related to the hypothesis. Statistics at the time of establishment If the distribution is approximated, then the two-sided p-value calculated using the normal distribution approximation method is... ,in for The observed values, The cumulative distribution function representing the standard normal distribution; If the test statistic is Then in Given the conditions, the normal distribution approximation method uses the expectation and variance as follows: and The normal distribution is related to the hypothesis. Statistics at the time of establishment If the distribution is approximated, then the two-sided p-value obtained using the normal distribution approximation method is: ,in for Observed values; When calculating the statistical p-value using the saddle point approximation method, the specific steps are as follows: If the test statistic is Then, under the assumption When it was established They are considered to be independent and follow a binomial distribution. A sequence of random variables, in Given the conditions, in order to estimate In the assumption The cumulative generating function at the time of its establishment is first estimated. Moment generating function: The estimate of the cumulant generating function is: Its first and second derivatives are as follows: Therefore, in When it was established The estimate of the cumulant generating function is: Its first and second derivatives are as follows: For any given real number First calculate to make Established Then calculate to get and In order to approximate Under the original hypothesis The distribution below, according to the Barndorff-Nielsen saddle point approximation formula, uses... To approximate the probability Then for a given test statistic Observations The bilateral p-values ​​calculated using the saddle point approximation algorithm described above are... Record the p-value on the left side. and the p-value on the right They are respectively: in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate; If the test statistic is Then, under the assumption When it was established They are considered to be independent and follow a binomial distribution. A sequence of random variables, in Given the conditions, In the assumption The estimate of the cumulant generating function when it holds true is: Its first and second derivatives are as follows: For any given real number First calculate to make Established Then calculate to get and ,use To approximate the probability For a given Observations The bilateral p-values ​​calculated using the saddle point approximation algorithm described above are... Record the p-value on the left side. and the p-value on the right They are respectively: in and They represent the probabilities obtained by the above saddle point approximation algorithm, respectively. and The estimate; The statistical p-value is calculated using a hybrid test strategy that combines the normal distribution approximation method with the saddle point approximation method, as detailed below: make It is a pre-selected positive number; If we test the null hypothesis The statistic is , for The observed values, if If the p-value is positive, the normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value. If we test the null hypothesis The statistic is , for The observed values, if If the p-value is positive, the normal distribution approximation method is used to calculate the p-value; otherwise, the saddle point approximation method is used to calculate the p-value. In step S10, the significance of the marginal gene-environment interaction at the locus will be tested, i.e., the null hypothesis will be tested. The obtained statistical p-value is related to the test site at... The significance of ancestry-specific marginal gene-environment interactions within a ancestry, i.e., testing the hypothesis. Received The statistical p-values ​​are combined using the Cauchy combination method to obtain the final statistical p-value, so as to achieve genome-wide gene-environment interaction analysis; After obtaining the final statistical p-value, a genome-wide gene-environment interaction analysis was performed at a given significance level.

Images