1.多變量孟德爾隨機(jī)化（MVMR），內(nèi)容遴選自公眾號生信與臨床，相對于育種等其他公眾號更易懂好操作

這里有幾點(diǎn)注意事項(xiàng)：

?（1）選擇的幾個暴露彼此之間應(yīng)該存在一定的關(guān)聯(lián)，比如高密度脂蛋白（HDL），低密度脂蛋白（LDL）和甘油三酯（TG）等。

（2）一般來說，SNP只要與其中一個暴露強(qiáng)相關(guān)即可。

（3）應(yīng)該提取SNP在所有暴露和結(jié)局中的信息，不能有缺失。

（4）進(jìn)行MVMR分析需要關(guān)于暴露和結(jié)局的完整GWAS結(jié)果。

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其實(shí)很簡單，就是在回歸模型里多加了幾個自變量而已。

接下來，我將以“TwoSampleMR”包及其提供的數(shù)據(jù)進(jìn)行簡單實(shí)踐：

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

?

我們可以看一下暴露數(shù)據(jù)的格式，這里總共有144個SNP，每個SNP對應(yīng)于三個暴露的信息，因此會有上述的144*4 = 432行數(shù)據(jù)。從下圖，我們可以看出SNP rs10490626和HDL沒有關(guān)聯(lián)，但是它和LDL強(qiáng)相關(guān)，所以我們需要把它納入研究。

?

oucome_dat <- extract_outcome_data(exposure_dat$SNP, id_outcome) #提取結(jié)局?jǐn)?shù)據(jù)mvdat <- mv_harmonise_data(exposure_dat, ?outcome_dat) # 對數(shù)據(jù)進(jìn)行harmonisationres <- mv_multiple(mvdat) #計(jì)算結(jié)果，大家也可以試試mv_residual()函數(shù)res

結(jié)果如下圖所示:

從結(jié)果我們不難看出，在矯正了其它脂質(zhì)的影響之后，HDL和冠心病的發(fā)病風(fēng)險沒有關(guān)聯(lián)了。如果大家單獨(dú)算一對一的MR，會發(fā)現(xiàn)這三個脂質(zhì)均和冠心病相關(guān)，但在MVMR中HDL的顯著性消失了，這和一些臨床觀測的結(jié)果一致（有臨床研究觀察到HDL并沒有之前認(rèn)為的那種保護(hù)作用）。

?在進(jìn)行MVMR研究時，我們不建議使用很多的暴露，因?yàn)檫@會帶來比較嚴(yán)重的共線性問題，一般3~5個為宜。如果暴露間的共線性問題比較嚴(yán)重，建議使用“TwoSampleMR”包的mv_lasso_feature_selection()函數(shù)來幫助你去除不必要的暴露。

2.中介孟德爾隨機(jī)化

中介分析：在因果推斷中，我們不僅對暴露對結(jié)局的影響程度感興趣，而且對暴露影響該結(jié)局的機(jī)制或途徑感興趣。中介分析試圖確定暴露影響結(jié)局的因果途徑及其相對重要性。當(dāng)暴露難以或不可能干預(yù)時，中介分析可以確定哪些因素介導(dǎo)了該暴露與結(jié)局之間的關(guān)系，從而能夠?qū)@些中介進(jìn)行干預(yù)以減輕暴露的影響。暴露對感興趣結(jié)局的因果效應(yīng)（包括通過中介的潛在效應(yīng)），是暴露對結(jié)局的總效應(yīng)(Total effect, TE)。TE可以分解為兩部分，一部分是暴露對結(jié)局的直接效應(yīng)(Direct effect, DE)，不通過模型中包含的中介起作用，由下圖中從X到Y的路表示，標(biāo)記為β1。另一部分是間接效應(yīng)(Indirect effect, IE)，即暴露對結(jié)局僅通過模型中包含的中介產(chǎn)生影響，在下圖中是X通過M到Y(jié)表示的，這條邊的效應(yīng)大小為αβ2。暴露對結(jié)局的總效應(yīng)是這兩種效應(yīng)的和，即TE=DE+IE，中介介導(dǎo)的總效應(yīng)比例可以計(jì)算為“間接效應(yīng)/總效應(yīng)(IE/TE)”。

中介效應(yīng)的孟德爾隨機(jī)化也稱為Two-Step MR，操作還是挺簡單的：開始分析X與Y的因果關(guān)系，無論是否是陽性；再觀察X與中介M的因果，中介M與Y的因果（注意effect的方向不能搞反）。

接下來就是如何解釋。優(yōu)秀的結(jié)果是：X有causal effect on Y， 同時X對M， M對Y也有同樣的Effect，同正或同負(fù)，這樣是可以解釋通的；另一個種是并未觀察到X與Y直接的因果關(guān)系，但是觀察到了通過中介路徑的兩個顯著的因果關(guān)系，這也是可以很好的解釋的。那么我們在做雙樣本的時候，X相對限制了我們對結(jié)果的解釋，那么多變量孟德爾隨機(jī)化（MVMR）的引入就很有必要了，可以進(jìn)行一個綜合的解釋，文章也會亮眼很多很多。

基礎(chǔ)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")

年代久遠(yuǎn)，將就看一下（很亂），建議自己多鉆研R包的help部分，那里有詳細(xì)的解釋也能助你更好的理解代碼。