A Mendelian randomization approach to causal inference based on summary statistics of genome-wide association studies

Abstract

The present invention discloses a Mendelian randomization causal inference method based on the summary statistical data of whole genome association study, comprising: S1, integrating the inference model, the polygenic effect model and the estimation error model, adjusting the LD effect, and obtaining the MR-MU model; S2, collecting the GWAS data of the studied phenotype and its risk factor candidate set; S3, screening the risk factor candidate set; S4, excluding the characteristics showing weak pleiotropy evidence; S5, selecting appropriate SNPs from the GWAS data as instrumental variables; S6, performing parameter estimation on the conditional likelihood function of the MR-MU model, respectively obtaining two sets of parameter estimation values ​​under the conditions of fixed and unfixed causal effects, and the estimation value of the causal effect under unfixed conditions, and calculating the function values ​​of the conditional likelihood function corresponding to the two sets of parameter estimation values; S7, calculating the likelihood ratio test statistic, identifying the risk factors of the studied phenotype, and statistically inferring whether there is a causal effect, so as to improve the accuracy and reliability of the causal inference.

Description

Technical Field

[0001] The present invention relates to a Mendelian randomization causal inference method based on the summary statistical data of whole genome association study, and belongs to the application field of combining biological genetic data analysis with computers. Background Art

[0002] Inferring the causal relationship between risk factors and related phenotypes (i.e., outcome variables) plays a vital role in biomedical, behavioral, and social sciences. A major challenge in causal inference is the influence of confounding factors, which can affect both risk factors and outcome variables and induce associations between them when they lack a causal relationship. In order to eliminate the influence of confounding factors, Mendelian randomization (MR) is currently introduced for causal inference in observational studies. In addition, genome-wide association studies (GWASs) have attracted much attention for their ability to reveal the relationship between genetic variation and phenotype. In particular, combined with GWASs data, MR can make more accurate causal inferences and promote progress in fields such as biomedicine.

[0003] Specifically, MR is a statistical method that uses instrumental variables (IVs), such as single nucleotide polymorphisms (SNPs), to infer the causal relationship between risk factors and outcome variables. Traditional MR relies on strong assumptions of valid IVs, including (i) IVs are associated with risk factors; (ii) IVs are independent of confounders; and (iii) IVs affect outcome variables only through risk factors. However, when applying GWAS summary statistics for MR, these IV assumptions are too strong to be met. In addition, when using MR for causal inference, it is important to distinguish between observed and unobserved confounders.

[0004] First, observed confounders are a set of variables that have been observed to be associated with both risk factors and outcome variables. If these variables are not properly considered in MR analysis, they may bias the estimated causal effect between risk factors and outcome variables.

[0005] Second, unobserved confounders refer to a set of unobserved variables that lead to spurious associations between risk factors and outcome variables. A common source is unknown genetic pathways shared between risk factors and outcome variable characteristics. This can also lead to genetic correlations between SNP effect sizes, even though there is no causal relationship between them.

[0006] like Figure 2As shown in the figure, the latest MR methods are divided into two categories: univariate MR (SVMR) and multivariate MR (MVMR). However, the latest MR methods still have some major limitations, which are introduced as follows: First, the SVMR method focuses on the instrumental variables of a single risk factor and relies heavily on the imposed model assumptions for causal inference. They have not yet fully utilized the available information from observed confounders to consider gene pleiotropy, resulting in unsatisfactory performance. Second, the MVMR method may be affected by more invalid instrumental variables than the SVMR method. The MVMR method usually selects instrumental variables for each risk factor and then takes their union. As the number of risk factor variables increases, the joint operation will bring too many instrumental variables. A considerable proportion of the instrumental variables may be associated with unobserved confounders. For example, they can affect risk factors and outcome variables through unknown shared pathways. Without considering genetic correlations, induced associations may be misunderstood as causal effects. Third, the latest MR methods largely ignore the effects of unobserved confounders hidden in the summary statistics of GWAS, such as population stratification, sample overlap, and so on.

[0007] Therefore, it is urgent to develop and design a new MR method based on GWASs summary statistics to improve the accuracy and reliability of causal inference by comprehensively considering confounding factors. Summary of the invention

[0008] In response to the above-mentioned existing technical problems, the present invention provides a Mendelian randomization causal inference method based on the summary statistical data of whole-genome association studies, which achieves the technical goal of improving the accuracy and reliability of causal inference by simultaneously considering both observed confounding factors and unobserved confounding factors.

[0009] To achieve the above technical objectives, the present invention provides a Mendelian randomization causal inference method based on the summary statistical data of whole genome association studies, comprising the following steps:

[0010] S1. Integrate the inference model for dealing with observed confounders, as well as the polygenic effect model and estimation error model for dealing with unobserved confounders, establish the MR-MU initial model for estimating the causal effect of the target risk factor on the outcome variable, and adjust the LD effect to obtain the MR-MU model;

[0011] S2. Collect GWAS data of the studied phenotype and its candidate risk factor set according to the studied phenotype;

[0012] S3, using the MR-APSS method to screen the risk factor candidate set and obtain potential risk factors;

[0013] S4. Eliminate characteristics showing weak pleiotropic effects from potential risk factors to obtain the final risk factors;

[0014] S5. Based on the final risk factors, select appropriate SNPs from GWAS data as instrumental variables;

[0015] S6. Based on the final risk factors and instrumental variables, the conditional likelihood function of the MR-MU model is estimated, and two sets of parameter estimates are obtained under the conditions of fixed and unfixed causal effects, as well as the estimated value of the causal effect under unfixed conditions, and the function values ​​of the conditional likelihood function corresponding to the two sets of parameter estimates are calculated;

[0016] S7. Based on the two sets of parameter estimates and the function values ​​of their corresponding likelihood functions, the likelihood ratio test statistic is calculated to identify the risk factors for the phenotype under study and to statistically infer whether there is a causal effect.

[0017] In the present invention, further, in S1, an inference model for processing observed confounders, and a polygenic effect model and an estimation error model for processing unobserved confounders are integrated to establish an MR-MU initial model for estimating the causal effect of the target risk factor on the outcome variable, comprising the following steps:

[0018] Assumptions represents the target risk factor; represents the outcome variable; express observed confounders; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and estimation errors; Indicates indivual.

[0019] Then the formula of the MR-MU initial model is as follows:

[0020]

[0021] in, represents a Bernoulli variable. If genotype instruments have non-zero instrumental variable intensities , then it is 1, otherwise it is 0; Indicates Genotype instrumental variable pairs residual effect; Indicates Genotype instrumental variable pairs indivual residual effect; express right causal effect; express right causal effect; Indicates Genotype instrumental variable pairs The residual direct effect of , and Respectively represent Genotype instrumental variable pairs , and polygenic effects; , and They represent the estimation errors respectively; , , Represent theoretical parameters , , The estimated value of .

[0022] It is assumed that the formulas of the polygenic effect model and the estimation error model are as follows:

[0023]

[0024] in, Indicates the mean and the variance is The multivariate normal density function of Is the length The zero vector of ; Indicates Genotype instrumental variable pairs and an estimate of the standard error of ; Indicates the mean and the variance is The multivariate normal density function of ; The coefficient matrix representing the variance term of the multivariate normal density function.

[0025] Assume that the formula of the inference model is as follows:

[0026]

[0027] in, It means that the mean is and the variance-covariance matrix The multivariate normal density function of ; Yes and indivual of The variance-covariance matrix of the variables; express The variance-covariance matrix of .

[0028] The present invention further comprises the following steps: in S1, adjusting the LD effect to obtain the MR-MU model:

[0029] set up The first SNPs and The correlation between the SNPs was calculated and the relevant variables in the initial MR-MU model were re-expressed to obtain the formula of the MR-MU model as follows:

[0030]

[0031] in, , , , , , ; Indicates the genome indivual and indivual The correlation between Indicates Genotype instrumental variable pairs indivual residual effect; , and Respectively represent Genotype instrumental variable pairs , and polygenic effects; Indicates Genotype instrumental variable pairs The remaining direct effect.

[0032] Combining formula (2) and (4), we get:

[0033]

[0034] in, represents the conditional probability function; represents the conditional multivariate normal density function; Indicates LD score of each SNP; , ; The dimension is , the identity matrix with all diagonal elements being 1; Indicates The probability that a valid instrument has a non-zero instrumental variable strength.

[0035] In the present invention, in S3, the screening criterion is: setting the threshold of the significance level to 0.05.

[0036] In the present invention, the exclusion criteria in S4 are: value For traits with less than five IVs, IVs are single nucleotide polymorphisms.

[0037] The present invention further comprises the following steps:

[0038] First use a moderate value threshold or equivalently use selection criteria , a set of candidate IVs are selected from the GWAS data; then the PLINK process is applied to ensure the independence between the candidate IVs to obtain the final IVs.

[0039] The summary statistics of the GWAS after IV selection are expressed as the following set:

[0040]

[0041] in, Indicates the use of selection criteria threshold The number of SNPs after IV selection when For the The probability that a valid IV has a non-zero IV strength; and As a collection of GWAS effect estimators after IV selection; As a collection of standard errors estimated by GWAS; As IV after selection The LD scores of the SNPs are summed.

[0042] In the present invention, further, in S6, the conditional likelihood function of the MR-MU model includes:

[0043] Will The elements in are considered random and are assigned a prior distribution as follows:

[0044]

[0045] in, express The variance coefficient parameter of , and obeys the Gamma distribution; , represents the hyperparameters of the Gamma distribution and specifies and ; represents the multivariate normal density function.

[0046] Under formulas (4), (5), (6) and (7), the conditional likelihood function formula of the MR-MU model is as follows:

[0047]

[0048] in, represents all latent variables;

[0049] and Respectively and The estimated value of As a collection of standard errors estimated by GWAS; and Respectively and The estimated value of As a collection of standard errors estimated by GWAS; As IV after selection The collection of LD scores of SNPs; Represents a set of parameters to be estimated.

[0050] The present invention further comprises that in S6, the conditional likelihood function of the MR-MU model is estimated to obtain two sets of parameter estimates under the conditions of fixed and unfixed causal effects, and an estimate of the causal effect under unfixed conditions, including:

[0051] S6-1. Estimation and , the steps are as follows:

[0052] Assuming that the GWAS summary statistics are calculated using standardized genotypes and phenotypes, and The formula is as follows:

[0053]

[0054] in, The variance matrix representing the multivariate normal density function of the estimation errors; Indicates The heritability of a phenotype, Indicates Phenotype and The co-inheritance rate of SNPs; M is the number of SNPs; represents the coefficient matrix of the variance term of the multivariate normal density function of the estimation error, express Middle diagonal elements, express Middle Line Elements of a column.

[0055] Use single-trait LD score regression to estimate and Let the diagonal elements in and For the diagonal elements, then we can calculate it using the following formula and :

[0056]

[0057] in, Indicates the expected value of the variable in brackets; Indicates The phenotype of SNPs were obtained through GWAS analysis ,and , which represents the ratio of the effect estimate to the estimated standard deviation; represents the sample size of the first phenotype, Indicates Sample size for each phenotype; Indicates LD score of each SNP; Indicates Sample size for each phenotype.

[0058] right and For the non-diagonal elements in and For the Row, No. The elements of the column are calculated according to the following formula and :

[0059]

[0060] in, , Respectively represent Phenotype, The phenotype of SNPs were obtained through GWAS analysis ,and , , both represent the ratio of the effect estimate to the estimated standard deviation; Indicates Phenotype and Public sample size for each phenotype; Indicates The phenotype and Genetic correlations of phenotypes; Indicates LD score of each SNP; express Middle Row, No. Elements of a column.

[0061] Will and The diagonal elements and off-diagonal elements are put together to form a matrix, and we get and The estimation formula is as follows:

[0062]

[0063] in, and Respectively and The estimated value of .

[0064] The present invention further comprises, in S6, performing parameter estimation on the conditional likelihood function of the MR-MU model, obtaining two sets of parameter estimation values ​​under the conditions of fixed and unfixed causal effects, and an estimation value of the causal effect under unfixed conditions, and further comprising:

[0065] S6-2, use the variational EM algorithm to estimate the parameters, and get are two sets of parameter estimates for zero and nonzero conditions, and The estimated value when it is non-zero is as follows:

[0066] Omit known or given ,set up For the selection criteria The conditional likelihood function in formula (8) is written as follows:

[0067]

[0068] in, represents KL divergence; log represents natural logarithm; Represents a collection of GWAS summary statistics after IV selection.

[0069] And assume the variational distribution is broken down into:

[0070]

[0071] in, Indicates that the variational distribution function is only related to the relevant item; Indicates that the variational distribution function is only related to the relevant item; Indicates the symbol for continuous multiplication; , , They represent the first Instrumental variables , and residual effect; is the conditional distribution term in the variational distribution function.

[0072] S6-2-1. Initialize parameter values.

[0073] S6-2-2. Repeat steps E and M in a loop until convergence conditions are reached to obtain parameter estimates.

[0074] The calculation method of step E is as follows:

[0075] E-1. Calculation , the formula is as follows:

[0076]

[0077]

[0078]

[0079] in, Express Take the expected value; ,so ;

[0080]

[0081]

[0082]

[0083]

[0084] E-2. Calculation , , the formula is as follows:

[0085]

[0086] in, It is a matrix operation function that sums the diagonal elements of the matrix.

[0087] E-3. Calculation , the formula is as follows:

[0088]

[0089]

[0090] in, is the matrix determinant calculation function; is the normal cumulative distribution function; Indicates the number of iterations in the EM algorithm; It is LD score of each SNP; express C The estimated value of the first element of the diagonal; express The element at row 1 and column 1 in .

[0091] The calculation method of the M step is as follows, and in M-2 is not executed when:

[0092] M-1, Update , , the formula is as follows:

[0093]

[0094]

[0095]

[0096]

[0097]

[0098] in, express The transpose of

[0099]

[0100]

[0101] M-2, Update , the formula is as follows:

[0102]

[0103] M-3, Update , the formula is as follows:

[0104]

[0105] M-4, Update , the formula is as follows:

[0106]

[0107]

[0108] in, represents a unit vector; Perform Cholesky decomposition to obtain the decomposed matrix ,satisfy .

[0109] And, calculate and , substitute The current value of and Go to update , the formula is as follows:

[0110]

[0111] M-5, Update , the formula is as follows:

[0112]

[0113] Repeat the E and M steps until the convergence condition is met, that is, .

[0114] The present invention further comprises: in S6, calculating the function value of the conditional likelihood function corresponding to the two sets of parameter estimation values ​​comprises:

[0115] S6-3. Calculation , the formula is as follows:

[0116]

[0117] in,

[0118]

[0119]

[0120]

[0121]

[0122]

[0123]

[0124]

[0125] and,

[0126]

[0127]

[0128]

[0129]

[0130] in, Represents the cumulative distribution function of the standard normal distribution; represents the density function of the standard normal distribution; Represents the identity matrix.

[0131] In summary, the present invention proposes the MR-MU method, which is an MR method that takes into account both observed and unobserved confounders. The MR-MU method integrates three components: an inference model, a polygenic effect model, and an estimation error model, thereby obtaining the MR-MU model. Among them, the inference model performs causal inference by explicitly considering observed confounders; the polygenic effect model considers the correlation of SNP effect sizes that may be caused by shared genetic pathways; and the estimation error model adjusts the impact of sample structure on estimation errors. In this way, by comprehensively considering confounders, the method of the present invention can improve the accuracy and reliability of causal inference, and has the following beneficial effects:

[0132] The present invention considers the influence of both observed and unobserved confounding factors, and comprehensively improves the reliability of causal inference. The present invention uses genome-wide summary statistics to estimate some parameters of the model, ensuring the identifiability of the model. The present invention takes into account the influence of selection bias and removes the influence of selection bias by constructing a conditional likelihood function. This is the key to ensuring the accuracy of causal effect estimation. The present invention allows for the inclusion of more instrumental variables of moderate strength, thereby improving statistical power (the statistical power here refers to the probability of the correct replacement hypothesis being accepted after the null hypothesis is rejected in the test hypothesis. Studies with insufficient statistical power often fail to identify major findings).

[0133] Compared with the prior art, the present invention mainly has the following technical advantages:

[0134] A. Unlike the SVMR method, the inference model of the MR-MU method of the present invention allows the inclusion of all measured confounding factors to consider the effects of genetic pleiotropy. Unlike the MVMR method, which simultaneously infers the causal effects of multiple risk factors, the MR-MU method of the present invention focuses on only one target risk factor at a time and selects instrumental variables based only on the risk factor of interest, rather than on the union of instrumental variables selected for each risk, which reduces the risk of invalid instrumental variables that may be introduced by the union operation in the MVMR method.

[0135] B. Unlike existing MR methods that only use a few strong instrumental variables, the MR-MU method of the present invention uses whole-genome summary data to estimate polygenic models and estimate error models.

[0136] C. The MR-MU method of the present invention takes into account unobserved confounding factors, such as population stratification and unknown shared genetic pathways, in its analysis. Under the assumption of LDSC, polygenic effects can be distinguished from confounding bias due to population structure, and the relevance of polygenic effects can be correctly considered.

[0137] D. Unlike most MR methods that only use instrumental variables with relatively high strength, the MR-MU method of the present invention allows the inclusion of more instrumental variables with moderate strength, thereby improving statistical power. BRIEF DESCRIPTION OF THE DRAWINGS

[0138] Figure 1 is a flow chart of the steps of the method of the present invention;

[0139] Figure 2 A comparison diagram of the model structure relationship between the MR-MU method provided by the present invention, the SVMR method and the MVMR method;

[0140] Figure 3 A schematic diagram of screening various features provided for an embodiment of the present invention;

[0141] Figure 4 A schematic diagram of model analysis results provided for an embodiment of the present invention. DETAILED DESCRIPTION

[0142] like Figure 1 , Figure 2 As shown, this embodiment proposes a Mendelian randomization causal inference method based on the summary statistical data of the whole genome association study, taking the identification of causal risk factors of coronary artery disease (CAD) as an example, the specific steps are introduced as follows.

[0143] S1. Integrate the inference model for dealing with observed confounders, as well as the polygenic effect model and estimation error model for dealing with unobserved confounders, establish the MR-MU initial model for estimating the causal effect of the target risk factor on the outcome variable, and adjust the LD effect to obtain the MR-MU model, so that the MR-MU method can be used to analyze and identify the risk factors of the phenotype under study. The steps are introduced as follows.

[0144] S1-1. Integrate the inference model, polygenic effect model and estimation error model to establish the initial MR-MU model.

[0145] Specifically, assuming represents the target risk factor, represents the outcome variable, express observed confounders that may be related to and All are related; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and estimation errors; Indicates indivual; , indicating the indivual.

[0146] The formula of the MR-MU initial model is as follows:

[0147]

[0148] in, represents a Bernoulli variable, if genotype instruments have non-zero instrumental variable intensities , then it is 1 and used for causal inference, otherwise it is 0; Indicates Genotype instrumental variables for target risk factors The residual effect of is used to measure the strength of the instrumental variable; Indicates Genotype instrumental variable pairs Observed confounders residual effect; Indicates target risk factors For outcome variables causal effect; express For outcome variables causal effect; Indicates Genotype instrumental variables for outcome variables The residual direct effect of , and Respectively represent Genotype instrumental variables for target risk factors , observed confounders and outcome variables polygenic effects; , and They represent the estimation errors respectively; , , Respectively represent the corresponding theoretical parameters , , The estimated value of .

[0149] S1-2. Based on the MR-MU initial model, the LD effect is adjusted to obtain the MR-MU model.

[0150] First, in order to remove the influence of unobserved confounding factors, the formulas of the polygenic effect model and the estimation error model are assumed to be as follows:

[0151]

[0152] in, Indicates the mean and the variance is The multivariate normal density function of Indicates the length is The zero vector of ; Indicates Genotype instrumental variables for target risk factors and outcome variables An estimate of the standard error of ; It means the mean and the variance is The multivariate normal density function of ; The coefficient matrix representing the variance term of the multivariate normal density function of the estimated errors.

[0153] It is noteworthy that the present invention allows The off-diagonal elements in are non-zero to take into account the correlation between gene effects caused by unobserved confounders, observed confounders or causal relationships. In addition, the present invention allows The diagonal elements in are greater than 1 to account for residual confounding factors such as population stratification that may not be fully adjusted in the GWAS summary statistics, and the off-diagonal elements can be non-zero, because population stratification, etc. may cause correlation between error terms. Different from the existing Mendelian randomization method that only uses a few strong instrumental variables, the present invention combines the target risk factors , outcome variables and observed confounders The summary statistics of the entire genome are used to estimate the parameters in the background model and , to ensure the identifiability of the model.

[0154] Secondly, after considering the unobserved confounders, the MR-MU method of the present invention begins to consider the observed confounders by using multiple regression to explain the observed confounders and using instrumental variables IV SNPs as target risk factors The inference model adopts the same strategy as the MVMR method.

[0155] Specifically, for Target risk factors The instrumental variable (IV) ), the inference model relies on the following assumptions: Target risk factors About, show ; Excluding observed confounders IV is independent of any other observed confounders other than The effect of IV strength ) and IV for observed confounders The impact of ) can be correlated with each other. These hypotheses suggest that: risk factor variables and The observed confounders in can be causally dependent on each other, The causal relationship between variables can be unidirectional or bidirectional; Direct effects on outcome variables Represents effects associated with unrelated pleiotropic effects and should be independent of these IV effects On the basis of these assumptions, the inference model is reasonable.

[0156] Then assume that the formula of the inference model is as follows:

[0157]

[0158] in, It means that the mean is and the variance-covariance matrix The multivariate normal density function of Indicates the length is The zero vector of represents the variance-covariance matrix;

[0159] Includes target risk factors and Observed confounders of The variance-covariance matrix of the variables, their variances are expressed as and ; express The variance-covariance matrix of .

[0160] It is worth noting that the variance ratio Can be used as an indicator of target risk factors The instrumental variable is an indicator of the strength of the pleiotropic signal in the corresponding measured variable. and The covariance between elements is set to zero because Should be independent of Any element in .

[0161] Furthermore, so far we have assumed that the SNPs are independent of each other. In reality, the SNPs are associated due to linkage disequilibrium (LD), which means that the GWAS summary statistics are actually estimates of the marginal true effect of the SNP on the phenotype.

[0162] set up The first SNPs and The correlation between the SNPs can be expressed by re-expressing the relevant variables in the initial MR-MU model, and the formula of the MR-MU model is as follows:

[0163]

[0164] in, represents a Bernoulli variable; , , , , , ; Indicates the genome SNPs and The correlation between SNPs; Indicates Genotype instrumental variable pairs Observed confounders residual effect; , and Respectively represent Genotype instrumental variables for target risk factors , observed confounders and outcome variables The polygenic effect of , and Synonymous, the reason why it is used here Instead of , is to In To make a distinction; Indicates Genotype instrumental variables for outcome variables The residual direct effect of Represents the weighted summation of the polygenic effects corresponding to all instrumental variables, which is equivalent to , and the same applies to other items.

[0165] And, combining formula (2) and (4), we get:

[0166]

[0167] in, represents the conditional probability function, As a conditional symbol, it means that the probability is calculated based on the parameter value following the conditional symbol (as known information); represents the conditional multivariate normal density function; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , , here Indicates the length is The zero vector of The dimension is , the identity matrix with all diagonal elements set to 1; Expressed as The probability that a valid instrumental variable has a non-zero instrumental variable strength; Indicates LD scores of SNPs.

[0168] The above MR-MU model conveys the key idea of ​​LD effect adjustment in the MR-MU method, that is, the contribution of gene effects to phenotypic variation is marked by the LD score, while the effects of residual confounding including population stratification and sample overlap are independent of the LD score. In addition, since the MR-MU model assumes that the contribution of polygenic effects to phenotypic variation is dominant, it can be estimated by using the LDSC of the whole genome summary statistics. and This estimation is crucial to ensure the identifiability of the MR-MU model.

[0169] S2. According to the phenotype under study, collect GWAS data of the phenotype under study and its risk factor candidate set, including the following sub-steps.

[0170] S2-1. According to the phenotype being studied, collect GWAS data of the phenotype being studied.

[0171] It should be noted that the GWAS data here mainly refers to the estimated genetic effect value or Z-score of the site, and the corresponding P-value.

[0172] In specific implementation, the GWAS data of the studied phenotype can be GWAS data from a single database or meta-analysis data of GWAS from multiple databases. In this embodiment, the studied phenotype is coronary artery disease (CAD), and the GWAS data set of CAD is selected from CAD (2015) of the CARDIoGRAMplusC4D Alliance. In addition, other CAD GWAS data sets can also be used for analysis, such as CAD from the UK Biobank (UKBB), and GWAS data sets that meta-analyze multiple GWAS data sets, etc., and are used to verify the accuracy and repeatability of the results.

[0173] S2-2. Based on the phenotype being studied, collect GWAS data for the candidate set of risk factors for the phenotype being studied.

[0174] It should be noted that the characteristics associated with the phenotype under study, in addition to the phenotype under study itself, also include a set of candidate characteristics that may potentially have a causal relationship with the phenotype under study, so this step collects GWAS data of the risk factor candidate set.

[0175] In the specific implementation, as shown in the table below, 46 different characteristics that may have a causal relationship with CAD are selected as candidate feature sets, including 24 blood and urine biomarkers, 12 blood cell composition characteristics, and nine other complex characteristics such as drinking, smoking, neuroticism, hypothyroidism, education level, systolic blood pressure, intelligence, height and BMI. These feature data are mainly derived from the GWAS summary statistics of the UK BioBank. In addition, potential risk factors for coronary artery disease (CAD) will be mined from these feature data in the future.

[0176]

[0177]

[0178] S3. Use the MR-APSS method to screen the risk factor candidate set obtained in S2 to obtain potential risk factors for the phenotype under study.

[0179] To improve the estimation efficiency, this step performs a pre-screening to determine which of the initial 46 features belong to the candidate feature set that shows a strong potential causal relationship with CAD, and obtains potential risk factors, that is, potential measured confounders. In addition, based on the effectiveness of considering unobserved confounders when evaluating causal relationships, the MR-APSS method is selected to evaluate the causal effect between each feature and CAD.

[0180] When implementing it, Figure 3 As shown in the figure, this step uses the CAD (2015) dataset from the CARDIoGRAMplusC4D consortium as a representative dataset for CAD, and successfully identifies 22 features by applying the MR-APSS method, including LDLdirect adjstantins, systolic blood pressure, glycosylated hemoglobin, body mass index, smoking, ALP enzyme, hypothyroidism, Blood HLSRC, neuroticism, HDL cholesterol, calcium, creatinine, education level, alanine aminotransferase, height, vitamin D, γ-glutamyl transferase, uric acid, triglycerides, sex hormone binding globulin, Blood NSCL, and aspartate aminotransferase. This process shows that at a significance level threshold of 0.05, these features have a significant causal effect on CAD (CAD (2015)) and are potential risk factors that should be focused on in the prevention and treatment of CAD. In addition, the final risk factors for CAD will be identified from these 22 features in the future.

[0181] S4. Eliminate characteristics showing evidence of weak pleiotropy from the potential risk factors obtained in S3 to obtain the final risk factors for the phenotype under study.

[0182] It should be noted that when measuring the impact of the potential risk factors obtained in S3 on CAD, the estimated results are easily affected by other confounding factors. For example, the same gene may affect multiple phenotypes, and the same phenotype may be affected by multiple genes. Therefore, it is necessary to measure which potential confounding factors each risk factor corresponds to. In particular, confounding factors can be divided into observed confounders and unobserved confounders based on whether they are observable. Among them, unobserved confounders have been processed in the MR-MU model through the model assumption method. In this embodiment, the potential risk factors obtained in S3 are also used as potential observed confounders. For example, when looking at the impact of LDL on CAD, its estimated value may be affected by the confounding effects of factors such as SBP and BMI. For each potential risk factor, the following method is used to screen the potential observed confounders it contains.

[0183] To minimize the risk of including relevant confounders, features showing evidence of weak pleiotropy were excluded from the set of potential observed confounders, and the criteria for exclusion were: failure to reach genomic significance in the GWAS data ( value ) with less than five IVs (single nucleotide polymorphisms (SNPs)). Since the method of the present invention intends to study the effects of potential risk factors on CAD, but these potential risk factors may also interact with each other, when using data to measure the effect of a potential risk factor on CAD, the measured effect may be mixed with the effects of other potential risk factors on CAD, so here we measure which other factors affect each potential risk factor.

[0184] In specific implementation, the target risk factors Among the selected instrumental variables (IV), those that failed to reach genomic significance in the GWAS data ( value ) with less than five IVs were excluded. Figure 3 As shown in the figure, the observed confounders set of each characteristic as a target risk factor in the MR-MU analysis is shown in detail, providing details on the treatment of confounders. Among them, the black grid represents the potential observed confounders included in the corresponding target risk factor. If it is the final risk factor for CAD, the characteristics corresponding to the black grid will be included in the observed confounders. In the example, the features corresponding to the white grids are not included in the observed confounding factors. This process helps ensure the accuracy and reliability of MR-MU analysis and makes the research results more precise and reliable by excluding irrelevant or weakly pleiotropic observed confounding factors.

[0185] S5. Based on the final risk factors for the studied phenotype obtained in S4, appropriate SNPs are selected as instrumental variables from the GWAS data obtained in S2.

[0186] It should be noted that in the MR-MU model, the instrumental variable (IV) is based on the target risk factor of interest. And, single nucleotide polymorphisms (SNPs) were used as IV.

[0187] When implementing it, first use a moderate The value threshold (usually ) or equivalently using selection criteria , a set of candidate IVs, i.e., a set of SNPs, were screened from the GWAS summary statistics. Then the PLINK procedure (using parameters: –clump-kb 1000, –clump-r2 0.001) was applied to ensure the independence between the candidate IVs, and the final IVs were obtained.

[0188] For ease of expression and to facilitate the derivation of the parameter estimation algorithm for the MR-MU model, the GWAS summary statistics after IV selection are expressed as a set, namely:

[0189]

[0190] in, Is to use the selection criteria threshold The number of SNPs after IV selection when For the The probability that an effective IV has a non-zero IV strength. And, for ease of expression, the following notation is introduced to represent the GWAS summary statistics after IV selection: , and As the collection of GWAS effect estimators after IV selection, as the collection of standard errors estimated by the GWAS, and As IV after selection The LD scores of the SNPs are summed.

[0191] In addition, Figure 3 shows the number of instrumental variables (IVs) corresponding to each risk factor in this example.

[0192] S6. Based on the final risk factors of the studied phenotype obtained in S4 and the instrumental variables obtained in S5, the conditional likelihood function of the MR-MU model is estimated, and two sets of parameter estimates are obtained under the conditions of fixed and non-fixed causal effects, as well as the estimated value of the causal effect under non-fixed conditions, and the function values ​​of the conditional likelihood function corresponding to the two sets of parameter estimates are calculated.

[0193] It should be noted that estimating GWAS summary statistics based on IV selection may be problematic in the MR-MU model. selected, they have an effect on observed confounders The impact may not be strong, which means Some elements of may be very small. If we try to estimate , which may lead to poor estimation efficiency and accuracy. To solve this problem, the method of the present invention does not directly estimate , but will The elements in are treated as random and a prior distribution is assigned to them as follows:

[0194]

[0195] in, express The variance coefficient parameter of is assumed to be Gamma distributed, and the two hyperparameters of Gamma distribution are , , and specify the hyperparameters and ; represents the multivariate normal density function.

[0196] Furthermore, another potential problem with instrumental variable (IV) selection involves the “winner’s curse”, in which IV selection based on statistical significance may bias the estimate of causal effects by affecting the probability density of the GWAS summary statistics of the selected SNPs. To correct for the “winner’s curse”, consider maximizing a conditional likelihood function based on instrumental variable selection conditions instead of maximizing the conventional likelihood, that is, under formulas (4), (5), (6) and (7), the conditional likelihood function of the MR-MU model is Written as the following formula:

[0197]

[0198] in, represents all latent variables; and Respectively and The estimated value of As a collection of standard errors estimated by GWAS; Represents a set of parameters to be estimated; represents the integration of z.

[0199] S6-1. Estimation and , the steps are as follows.

[0200] Assuming that the GWAS summary statistics are calculated using standardized genotypes and phenotypes, we can and The specific form is written as:

[0201]

[0202] in, Indicates The heritability of the shape, Indicates The shape and The co-inheritance rate of each property; M is the number of SNPs; Indicates the diagonal in C elements, express Middle Line Elements of a column.

[0203] Use single-trait LD score regression (LDSC) to estimate and diagonal elements in . LDSC assumes a random effects model to model the polygenic genetic architecture of complex phenotypes. Under the assumption of LDSC, polygenic effects can be marked by LD, but observed confounding factors (such as population stratification or hidden kinship) are not related to LD. Specifically, LD marks the effect, which means The expected square of Compared with its LD score is proportional to SNP with SNP The correlation between them.

[0204] set up and For the diagonal elements, the calculation formula is as follows: and :

[0205]

[0206] in, Indicates the expected value of the variable in brackets; Indicates The phenotype of SNPs were obtained through GWAS analysis , and the mathematical expression is , which is the ratio of the effect estimate to the estimated standard deviation; represents the sample size of the first phenotype; Indicates LD score of each SNP; Indicates Sample size for each phenotype.

[0207] right and For the non-diagonal elements in and For the Row, No. The elements of the column are calculated as follows: and :

[0208]

[0209] in, , Respectively represent Phenotype, The phenotype of SNPs were obtained through GWAS analysis ,and , , which represents the ratio of the effect estimate to the estimated standard deviation; Indicates Phenotype and Public sample size for each phenotype; Indicates LD scores of SNPs.

[0210] In specific implementation, and For example, this relationship is precisely given by the following formula:

[0211]

[0212] in, Indicates the first phenotype, SNPs were obtained through GWAS analysis ,and , which represents the ratio of the effect estimate to the estimated standard deviation; represents the heritability of the first trait; express Specifically, the target risk factor that includes the observed The square of the GWAS summary statistics Regressing to the LD score, we can get and The estimated values ​​of the diagonal elements in , Similarly, and The other diagonal elements in can also be obtained in this way.

[0213] Finally, and The matrix composed of the diagonal elements and off-diagonal elements is and C are estimated, and the estimated results are recorded as and , the formula is as follows:

[0214]

[0215] S6-2. Use the variational EM algorithm to estimate parameters and obtain the causal effect Two sets of parameter estimates under fixed and unfixed conditions, and causal effects Estimated value under non-fixed conditions.

[0216] It should be noted that the polygenic effect term and the error term are not considered as latent variables in the estimation algorithm, and the parameters in these two terms and C have already been estimated in advance by LDSC using genome-wide summary statistics, so there is no need to consider them as latent variables.

[0217] Since we can directly maximize the conditional likelihood function in formula (8) to obtain The optimal solution of is difficult, so the present invention develops a variational inference maximization algorithm (Variational EM) to estimate For simplicity, the known or given ,set up For the selection criteria The conditional likelihood function in formula (8) can be written as follows:

[0218]

[0219] in, represents KL divergence; It means taking the natural logarithm; represents conditional probability; Represents a collection of GWAS summary statistics after IV selection.

[0220] It can be seen that the above inequality is due to the fact that the Kullback-Leibler (KL) divergence KL(q(z|t)||Pr(z|Dt;θ)) is non-negative. Therefore, It provides a lower bound for the likelihood function, namely ELBO (evidence lowerbound). Unlike directly maximizing the marginal likelihood, the variational inference maximization algorithm maximizes the variational lower bound. .

[0221] The following introduces the derivation of the variational inference maximization algorithm (Variational EM). For the convenience of solving, it is assumed that the variational distribution can be broken down into:

[0222]

[0223] in, is the only variable in the variational distribution function that is The relevant items, is the only variable in the variational distribution function that is the relevant item; It is the multiplication symbol; , , They are the first Instrumental variables , and residual effect; is the conditional distribution term in the variational distribution function.

[0224] With this restricted form of variational distribution, becomes easy to calculate. Then we can apply the variational EM algorithm to maximize The algorithm is summarized as follows:

[0225] Step E: For fixed Optimizing ELBO .remember For a given The optimal variational distribution when ,but

[0226]

[0227] in, Indicates the parameter value that makes the function reach its maximum value.

[0228] Step M: For fixed optimization . Assume that the current parameters are , and Insert into In

[0229]

[0230] These two steps are repeated until convergence, and the parameter estimates are obtained.

[0231] It should be noted that the variational EM algorithm (Variational EM) used in parameter estimation in this step is to obtain parameter estimates after iterative convergence by running the E step and the M step multiple times. The estimated values ​​are then substituted into the conditional likelihood function of formula (9) to calculate the likelihood ratio test statistic for the likelihood ratio test.

[0232] Furthermore, this step will calculate two sets of parameter estimates, one set is fixed The parameter estimates under the condition that the other set is not fixed Parameter estimation under conditions where it is not fixed The parameter estimation under the condition can be obtained The difference between the two sets of estimates is that In this case, no update is required in the M step , the rest of the steps are the same, as follows.

[0233] S6-2-1. Initialize parameter values.

[0234] The estimated parameters include: and ,as well as As IV after selection The collection of LD scores of SNPs; It can be obtained from public LD reference data (e.g., the 1000 Genomes Project); correspond The length of the vector. And, set , , In addition, the unknown parameters involved in the calculation are initialized, and the values ​​of these unknown parameters will be updated in the subsequent iteration process, including: ; ; ; ; Initialize to , Put together The covariance of the matrix; for The element in the first row and first column of ; .

[0235] S6-2-2. After initializing the parameters, start executing steps E and M in sequence until the convergence condition is reached to obtain the estimated values ​​of each parameter, so as to further obtain the value of the conditional likelihood function.

[0236] The calculation method of the E step is as follows.

[0237] E-1. Calculation , The formula is as follows:

[0238]

[0239]

[0240]

[0241] in, Express Take the expected value, It means taking the expected value of all variables in the brackets. Indicates that Bring in In the expression, ,so .

[0242]

[0243]

[0244]

[0245]

[0246] E-2. Calculation , , the formula is as follows:

[0247]

[0248] in, Represents a matrix operation function, which sums the diagonal elements of the matrix.

[0249] E-3. Calculation , the formula is as follows:

[0250]

[0251]

[0252] in, Represents the matrix determinant calculation function; represents the normal cumulative distribution function; Indicates the number of iterations in the EM algorithm; Represents the estimated value of the first diagonal element in C; express The element at row 1 and column 1 in .

[0253] The calculation method of the M step is as follows, and No update required , M-2 is not executed.

[0254] M-1, Update , , the formula is as follows:

[0255]

[0256]

[0257]

[0258]

[0259]

[0260]

[0261] in, express The transpose of . The formula is as follows:

[0262]

[0263]

[0264] M-2, Update , the formula is as follows:

[0265]

[0266] M-3, Update , the formula is as follows:

[0267]

[0268] M-4, Update , the formula is as follows:

[0269]

[0270]

[0271] in, represents a unit vector. Perform Cholesky decomposition to obtain the decomposed matrix ,satisfy And, calculate and , substitute The current value of and Go to update , the formula is as follows:

[0272]

[0273] M-5, Update , the formula is as follows:

[0274]

[0275] Repeat the E and M steps until the convergence condition is met, that is, .

[0276] S6-3. Calculate the function values ​​of the conditional likelihood function corresponding to the two sets of parameter estimates respectively.

[0277] Specifically, calculate , the formula is as follows:

[0278]

[0279] in,

[0280]

[0281]

[0282]

[0283]

[0284]

[0285]

[0286]

[0287] and,

[0288]

[0289]

[0290]

[0291]

[0292] in, Represents the cumulative distribution function of the standard normal distribution; represents the density function of the standard normal distribution; Represents the identity matrix.

[0293] In specific implementation, Figure 3 shows The parameter estimation results under the assumption of , for different risk factors, different Estimated values, e.g. for LDL cholesterol The estimated value is 0.18, which is The estimated value is 0.211, and the corresponding And, in The parameter estimation under the assumption of is also done using the above estimation method, thus obtaining the corresponding The function value of .

[0294] S7. Based on the two sets of parameter estimates obtained in S6 and the function value of the corresponding conditional likelihood function, calculate the likelihood ratio test statistic to test the statistical significance of the causal effect of the risk factor, identify the risk factors for the phenotype under study, and statistically infer whether there is a causal effect.

[0295] It should be noted that in order to determine and Is there a causal effect between , by comparing a causal effect Fixed to zero (representing the null hypothesis H0: ) is similar to the model that allows for a causal effect Non-zero (represents the alternative hypothesis H1: ) to perform a likelihood ratio test. The two sets of parameter estimates obtained in S6 are Substitute into the following formula to obtain the likelihood ratio test statistic, thereby testing and Is there a causal effect between them?

[0296] set up Under the null hypothesis ( ) is the parameter estimate of Under the alternative hypothesis ( ), then the likelihood ratio test statistic is calculated as follows:

[0297]

[0298] Among them, under the null hypothesis, the test statistic , , are the likelihood values ​​after VEM converges under H0 and H1, respectively. If it is greater than the specified critical value, the null hypothesis is rejected. ; Otherwise, the null hypothesis is not rejected. The value can be obtained through the likelihood ratio test function in R language (a statistical analysis software). If the value is less than 0.05, it is considered to have a significant causal effect.

[0299] at last, Figure 4 The results were presented: From the 22 features screened, three major risk factors were successfully identified, including LDL cholesterol, systolic blood pressure, and body mass index. Among them, elevated LDL cholesterol (low-density cholesterol) and hypertension are two risk factors for CAD that have been supported by large-scale randomized controlled trials. Based on the Bonferroni correction Value threshold ( After considering observed and unobserved confounding factors, the present invention successfully reduced LDL cholesterol ( ,in express The estimated value of value and systolic blood pressure (SBP) , value ) were identified as risk factors for CAD (2015), such as Figure 3 The present invention also identified BMI as a risk factor for CAD (2015) ( , value ),like Figure 4 In addition, there is increasing evidence from large observational studies that patients with diabetes or obesity are more likely to die from cardiovascular disease. Consistent with this evidence, MR-MU methods have detected elevated glycated hemoglobin (HbA1c), which is commonly used in the diagnosis of diabetes, as a risk factor for CAD, such as Figure 4 shown.

[0300] In addition, it should be noted that the Mendelian randomization causal inference method based on the summary statistics of the whole genome association study of the present invention can also be applied to any other phenotypes besides CAD. For example, when studying the level of cognitive function of an individual, phenotypes such as depression, anxiety, and the level of cognitive function of parents can be used as potential risk factors, and the causal relationship between different risk factors and the level of cognitive function of an individual can be inferred by using the GWAS data of cognitive function for calculation.

Claims

1. A Mendelian randomization causal inference method based on summary statistics of genome-wide association studies, characterized in that: The steps include: S1. Integrate the inference model for dealing with observed confounders, as well as the polygenic effect model and estimation error model for dealing with unobserved confounders, establish the MR-MU initial model for estimating the causal effect of the target risk factor on the outcome variable, and adjust the LD effect to obtain the MR-MU model; S2. Collect GWAS data of the studied phenotype and its candidate risk factor set according to the studied phenotype; S3, using the MR-APSS method to screen the risk factor candidate set and obtain potential risk factors; S4. Eliminate characteristics showing weak pleiotropic effects from potential risk factors to obtain the final risk factors; S5. Based on the final risk factors, select appropriate SNPs from GWAS data as instrumental variables; S6. Based on the final risk factors and instrumental variables, the conditional likelihood function of the MR-MU model is estimated, and two sets of parameter estimates are obtained under the conditions of fixed and unfixed causal effects, and the function values ​​of the conditional likelihood function corresponding to the two sets of parameter estimates are calculated; S7. Based on the two sets of parameter estimates and the function values ​​of their corresponding likelihood functions, the likelihood ratio test statistic is calculated to identify the risk factors for the phenotype under study and to statistically infer whether there is a causal effect.

2. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 1, characterized in that: In S1, the inference model for dealing with observed confounders, as well as the polygenic effect model and the estimation error model for dealing with unobserved confounders are integrated to establish an initial MR-MU model for estimating the causal effect of the target risk factor on the outcome variable, including the following steps: Assumptions represents the target risk factor; represents the outcome variable; express observed confounders; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and standard errors; , Respectively represent Genotype instrumental variable pairs effect estimates and estimation errors; Indicates indivual; Then the formula of the MR-MU initial model is as follows: in, represents a Bernoulli variable. If genotype instruments have non-zero instrumental variable intensities , then it is 1, otherwise it is 0; Indicates Genotype instrumental variable pairs residual effect; Indicates Genotype instrumental variable pairs indivual residual effect; express right causal effect; express right causal effect; Indicates Genotype instrumental variable pairs The residual direct effect of , and Respectively represent Genotype instrumental variable pairs , and polygenic effects; , and They represent the estimation errors respectively; , , Represent theoretical parameters , , An estimated value of It is assumed that the formulas of the polygenic effect model and the estimation error model are as follows: in, It means the mean and the variance is The multivariate normal density function of Is the length The zero vector of ; Indicates Genotype instrumental variable pairs and An estimate of the standard error of ; It means the mean and the variance is The multivariate normal density function of ; The coefficient matrix of the variance term representing the multivariate normal density function; Assume that the formula of the inference model is as follows: in, It means that the mean is and the variance-covariance matrix The multivariate normal density function of ; Yes and indivual of The variance-covariance matrix of the variables; express The variance-covariance matrix of .

3. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 2, characterized in that: In S1, the LD effect is adjusted to obtain the MR-MU model, which includes the following steps: set up The first SNPs and The correlation between the SNPs was calculated and the relevant variables in the initial MR-MU model were re-expressed to obtain the formula of the MR-MU model as follows: in, , , , , , ; Indicates the genome indivual and indivual The correlation between Indicates Genotype instrumental variable pairs indivual residual effect; , and Respectively represent Genotype instrumental variable pairs , and polygenic effects; Indicates Genotype instrumental variable pairs The residual direct effect of Combining formula (2) and (4), we get: in, represents the conditional probability function; represents the conditional multivariate normal density function; Indicates LD score of each SNP; , ; The dimension is , the identity matrix with all diagonal elements set to 1; Indicates The probability that a valid instrument has a non-zero instrumental variable strength.

4. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 1, characterized in that: In S3, the screening criteria are: setting the threshold of the significance level to 0.

05.

5. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 1, characterized in that: In S4, the exclusion criteria are: failure to reach genomic significance in GWAS data value For traits with less than five IVs, IVs are single nucleotide polymorphisms.

6. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 3, characterized in that: The S5 comprises the following steps: First use a moderate value threshold or equivalently use selection criteria , a set of candidate IVs are selected from GWAS data, where IVs are single nucleotide polymorphisms; then the PLINK process is applied to ensure the independence between candidate IVs and obtain the final IVs; The summary statistics of the GWAS after IV selection are expressed as the following set: in, Indicates the use of selection criteria threshold The number of SNPs after IV selection when For the The probability that a valid IV has a non-zero IV strength; and As a collection of GWAS effect estimators after IV selection; As a collection of standard errors estimated by GWAS; As IV after selection The LD scores of the SNPs are summed.

7. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 6, characterized in that: In S6, the conditional likelihood function of the MR-MU model includes: Will The elements in are considered random and are assigned a prior distribution as follows: in, express The variance coefficient parameter of , and obeys the Gamma distribution; , represents the hyperparameters of the Gamma distribution and specifies and ; represents the multivariate normal density function; Under formulas (4), (5), (6) and (7), the conditional likelihood function formula of the MR-MU model is as follows: in, represents all latent variables; and Respectively and The estimated value of As a collection of standard errors estimated by GWAS; As IV after selection The collection of LD scores of SNPs; Represents a set of parameters to be estimated.

8. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 7, characterized in that: In S6, the conditional likelihood function of the MR-MU model is estimated to obtain two sets of parameter estimates under the conditions of fixed and unfixed causal effects, including: S6-1. Estimation and , the steps are as follows: Assuming that the GWAS summary statistics are calculated using standardized genotypes and phenotypes, and The formula is as follows: in, The variance matrix representing the multivariate normal density function of the estimation errors; Indicates The heritability of a phenotype, Indicates The phenotype and The co-inheritance rate of a phenotype; M is the number of SNPs; represents the variance term coefficient matrix of the multivariate normal density function of the estimation error, express Middle diagonal elements, express Middle Line Elements of a column; Use single-trait LD score regression to estimate and Let the diagonal elements in and For the diagonal elements, then we can calculate it using the following formula and : in, Indicates the expected value of the variable in brackets; Indicates The phenotype of SNPs were obtained through GWAS analysis ,and , which represents the ratio of the effect estimate to the estimated standard deviation; represents the sample size of the first phenotype, Indicates Sample size for each phenotype; Indicates LD score of each SNP; Indicates Sample size for each phenotype; right and For the non-diagonal elements in and For the Row, No. The elements of the column are calculated according to the following formula and : in, , Respectively represent Phenotype, The phenotype of SNPs were obtained through GWAS analysis ,and , , both represent the ratio of the effect estimate to the estimated standard deviation; Indicates The phenotype and Public sample size for each phenotype; Indicates The phenotype and Genetic correlations of phenotypes; Indicates LD score of each SNP; express Middle Row, No. Elements of a column; Will and The diagonal elements and off-diagonal elements are put together to form a matrix, and we get and The estimation formula is as follows: in, and Respectively and The estimated value of .

9. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 8, characterized in that: In S6, the conditional likelihood function of the MR-MU model is estimated to obtain two sets of parameter estimates under the conditions of fixed and unfixed causal effects, respectively, and further includes: S6-2, use the variational EM algorithm to estimate the parameters, and get For two sets of parameter estimates under zero and non-zero conditions, the steps are as follows: Omit known or given ,set up For the selection criteria The conditional likelihood function in formula (8) is written as follows: in, represents KL divergence; log represents natural logarithm; represents the set of GWAS summary statistics after IV selection; And assume the variational distribution is broken down into: in, Indicates that the variational distribution function is only related to the relevant item; Indicates that only the relevant item; Indicates the symbol for continuous multiplication; , , They represent the first Instrumental variables , and residual effect; is the conditional distribution term in the variational distribution function; S6-2-1, initialization parameter values; S6-2-2, repeatedly executing step E and step M until the convergence condition is reached to obtain the parameter estimation value; The calculation method of the E step is as follows: E-1. Calculation , the formula is as follows: in, Express Take the expected value; ,so ; and, E-2. Calculation , , the formula is as follows: in, It is a matrix operation function, which adds the diagonal elements of the matrix; E-3. Calculation , the formula is as follows: in, in, is the matrix determinant calculation function; is the normal cumulative distribution function; Indicates the number of iterations in the EM algorithm; It is LD score of each SNP; express C The estimated value of the first element of the diagonal; express The element in row 1 and column 1; The calculation method of the M step is as follows, and in M-2 is not executed when: M-1, Update , , the formula is as follows: in, in, express The transpose of and, M-2, Update , the formula is as follows: M-3, Update , the formula is as follows: M-4, Update , the formula is as follows: in, represents a unit vector; Perform Cholesky decomposition to obtain the decomposed matrix ,satisfy ; And, calculate and , substitute The current value of and Go to update , the formula is as follows: M-5, Update , the formula is as follows: Repeat the E and M steps until the convergence condition is met, that is, .

10. A Mendelian randomization causal inference method based on genome-wide association study summary statistics according to claim 9, characterized in that: In S6, calculating the function value of the conditional likelihood function corresponding to the two sets of parameter estimates includes: S6-3. Calculation , the formula is as follows: in, and, in, Represents the cumulative distribution function of the standard normal distribution; represents the density function of the standard normal distribution; Represents the identity matrix.

Images