多变量与中介孟德尔随机化概念与操作思路（附基础双样本双向MR代码）
1.多变量孟德尔随机化（MVMR），内容遴选自公众号生信与临床，相对于育种等其他公众号更易懂好操作
这里有几点注意事项：
（1）选择的几个暴露彼此之间应该存在一定的关联，比如高密度脂蛋白（HDL），低密度脂蛋白（LDL）和甘油三酯（TG）等。
（2）一般来说，SNP只要与其中一个暴露强相关即可。
（3）应该提取SNP在所有暴露和结局中的信息，不能有缺失。
（4）进行MVMR分析需要关于暴露和结局的完整GWAS结果。
lm(beta.outcome ~ beta.exposure, weights = 1/se.outcome^2, data = data)
在MVMR研究中，R代码如下 （以三个暴露为例）：
lm(beta.outcome ~ beta.exposure1 + beta.exposure2 + beta.exposure3, weights = 1/se.outcome^2, data = data)
MVMR其实很简单，就是在回归模型里多加了几个自变量而已。
接下来，我将以“TwoSampleMR”包及其提供的数据进行简单实践：
library(TwoSampleMR)id_exposure <- c("ieu-a-299","ieu-a-300", "ieu-a-302") # 三个暴露分别是HDL cholesterol，LDL cholesterol和Triglyceridesid_outcome <- "ieu-a-7"exposure_dat <- mv_extract_exposures(id_exposure)dim(exposure_dat)[1] 432 9
我们可以看一下暴露数据的格式，这里总共有144个SNP，每个SNP对应于三个暴露的信息，因此会有上述的144*4 = 432行数据。从下图，我们可以看出SNP rs10490626和HDL没有关联，但是它和LDL强相关，所以我们需要把它纳入研究。
oucome_dat <- extract_outcome_data(exposure_dat$SNP, id_outcome) #提取结局数据mvdat <- mv_harmonise_data(exposure_dat, outcome_dat) # 对数据进行harmonisationres <- mv_multiple(mvdat) #计算结果，大家也可以试试mv_residual()函数res
结果如下图所示:
从结果我们不难看出，在矫正了其它脂质的影响之后，HDL和冠心病的发病风险没有关联了。如果大家单独算一对一的MR，会发现这三个脂质均和冠心病相关，但在MVMR中HDL的显著性消失了，这和一些临床观测的结果一致（有临床研究观察到HDL并没有之前认为的那种保护作用）。
在进行MVMR研究时，我们不建议使用很多的暴露，因为这会带来比较严重的共线性问题，一般3~5个为宜。如果暴露间的共线性问题比较严重，建议使用“TwoSampleMR”包的mv_lasso_feature_selection()函数来帮助你去除不必要的暴露。
2.中介孟德尔随机化
中介分析：在因果推断中，我们不仅对暴露对结局的影响程度感兴趣，而且对暴露影响该结局的机制或途径感兴趣。中介分析试图确定暴露影响结局的因果途径及其相对重要性。当暴露难以或不可能干预时，中介分析可以确定哪些因素介导了该暴露与结局之间的关系，从而能够对这些中介进行干预以减轻暴露的影响。暴露对感兴趣结局的因果效应（包括通过中介的潜在效应），是暴露对结局的总效应(Total effect, TE)。TE可以分解为两部分，一部分是暴露对结局的直接效应(Direct effect, DE)，不通过模型中包含的中介起作用，由下图中从X到Y的路表示，标记为β1。另一部分是间接效应(Indirect effect, IE)，即暴露对结局仅通过模型中包含的中介产生影响，在下图中是X通过M到Y表示的，这条边的效应大小为αβ2。暴露对结局的总效应是这两种效应的和，即TE=DE+IE，中介介导的总效应比例可以计算为“间接效应/总效应(IE/TE)”。
中介效应的孟德尔随机化也称为Two-Step MR，操作还是挺简单的：开始分析X与Y的因果关系，无论是否是阳性；再观察X与中介M的因果，中介M与Y的因果（注意effect的方向不能搞反）。
接下来就是如何解释。优秀的结果是：X有causal effect on Y， 同时X对M， M对Y也有同样的Effect，同正或同负，这样是可以解释通的；另一个种是并未观察到X与Y直接的因果关系，但是观察到了通过中介路径的两个显著的因果关系，这也是可以很好的解释的。那么我们在做双样本的时候，X相对限制了我们对结果的解释，那么多变量孟德尔随机化（MVMR）的引入就很有必要了，可以进行一个综合的解释，文章也会亮眼很多很多。
基础MR代码：
install.packages('remotes')
remotes::install_github('MRCIEU/TwoSampleMR', force = TRUE)
library(TwoSampleMR)
library(data.table)
b <- subset(a,a$'P-value'<5e-6)
library(epiDisplay)
library(tidyr)
library(rio)
Pe_exp_dat <- read_exposure_data( filename = "EUR.csv",
sep = ",",
snp_col = "SNP",
beta_col = "Effect",
se_col = "StdErr",
effect_allele_col = "Allele1",
other_allele_col = "Allele2",
eaf_col = "EAF",
pval_col = "P-value",
samplesize_col = "samplesize",
phenotype_col = "exposure",
clump=F)
Pe_exp_dat_clumped <- clump_data(Pe_exp_dat, clump_kb = 10000, clump_r2 = 0.01)
NAFLD_out_dat <- fread('NAFLD.txt', header = T, fill = TRUE)
write.csv(NAFLD_out_dat, 'NAFLD_out_dat.csv')
CP_out_dat <- fread("CP_out_dat.csv", header = T)
NAFLD_out_dat <- separate(DM_out_dat, 'INFO', sep = '=', c('AF', 'eaf'))
DM_out_dat$pval = -DM_out_dat$pval
write.csv(DM_out_dat, 'DM_out_dat.csv')
EUR_outcome_dat <- read_outcome_data(snps = exposure$SNP,
filename = "EUR.csv",
sep = ",",
snp_col = "SNP",
beta_col = "Effect",
se_col = "StdErr",
effect_allele_col = "Allele1",
other_allele_col = "Allele2",
eaf_col = "af_alt",
pval_col = "P-value",
samplesize_col = "samplesize",
phenotype_col = "outcome"
)
COPD_outcome_dat$id.outcome = 'outcome'
dat <- harmonise_data(exposure_dat = exposure, outcome_dat = EUR_outcome_dat)
write.csv(dat,"dat.csv")
res <- mr(dat)
write.csv(res,"res.csv")
or <- generate_odds_ratios(mr_res = res)
write.csv(or,"or.csv")
h <- mr_heterogeneity(dat)
write.csv(h,"h.csv")
p <- mr_pleiotropy_test(dat)
write.csv(p,"p.csv")
年代久远，将就看一下（很乱），建议自己多钻研R包的help部分，那里有详细的解释也能助你更好的理解代码。
