R语言实现孟德尔随机化研究
孟德尔随机化（Mendelian Randomization, MR）是近几年流行起来的用来进行因果推断的有效方法，它以遗传变异为工具变量来推导结局和暴露的因果关系，能有效避免传统流行病学研究的混杂偏倚。我们首先看下什么是工具变量? 孟德尔随机化的定义是“使用遗传变异进行工具变量分析”。在孟德尔随机中，遗传变异被用作工具变量(IV)评估暴露对结局的因果效应，遗传变异满足工具变量的基本条件总结为：
（1） 遗传变异与暴露有关。
（2） 该遗传变异与暴露-结果关联的任何混杂因素均不相关。
（3） 该遗传变异不会影响结果，除非可能通过与暴露的关联来实现。
具体的构架图如下：
对于上面的框架具体的阐述：：
1. 工具变量G与混杂因素无关联；
2. 工具变量G与暴露因素X有关联；
3. 工具变量G与结局变量Y无关联，G只能通过变量X与Y发生关联。
上述方程的使用必须满足条件：
1. 变量X与Y之间的关联一定会受到潜在混杂因素的影响，但工具变量G与变量X以及G与变量Y之间无潜在混杂因素影响；
2. 变量X与结局Y之间的关联无法直接观察获得，因为无法直接测量变量X，但是G是可测量的，并且G与X直接的关联是已知的或者可测量的，并独立于其他因素而存在。
通过上面的假设，便引入了孟德尔随机化研究的基础，面对不同的数据有以下几种设计策略：
1. 独立样本MR(One-sampleMR)：该方法利用单一研究样本，通过使用2阶段最小二乘法回归模型(2-stage least—squares regression，2SLS)，定量估计暴露因素X与Y之间的关联效应大小。第一步：建立G—X回归模型，获得暴露因素预测值(predicted value，P）；第二步：构建P—Y的回归模型，即获得暴露因素预测值P和结局变量Y之间的回归方程。
2. 两样本MR(Two-sampleMR)：两样本MR的设计策略是建立在G—X和G—Y的关联研究人群来自相同人群的两个独立样本(如GWAS与暴露，GWAS与结局的关联数据。)，要求两样本具有相似的年龄、性别和种族分布特征，因为样本量较大，该方法可以获得更大的把握度。目前，两样本MR因为全球大量GWAS合作组的公共数据而被广泛使用。
3. 两阶段MR(Two-stepMR)：与两样本MR不同的是，两阶段MR需要使用遗传工具变量来评价因果关联的可能中间变量M(Mediation)，来探讨环境暴露因素(X)是否通过表观遗传指标(M)而导致疾病(Y)改变。也就是需要在结局和工具之间构建一个中间变量。
第一阶段，遗传工具变量G1独立于混杂因素，指代暴露因素X与结局Y之间的关联，并且必须经过中间变量M才能实现；
第二阶段，另一独立遗传工具变量G2作为中间变量M的指代工具，分析中间变量M与结局Y之间的关联
4. 基因一暴露交互作用MR(Gene-exposureinteractions）interactions)：MR研究设计还可以用于探讨基因一暴露因素在疾病发生中的交互作用现象，同时要求基因与结局的关联必须取决于暴露因素的状态。这种方法可以区分基因直接作用于结局，还是基因通过暴露因素而作用于结局。
接下来我们看下MR的评估模型：
1．敏感度分析(sensitivityanalysis)
2．MR-Egger回归分析：以使用MR—Egger回归分析的方法来评价基因多效性带来的偏倚，MR-Egger回归直线的斜率可以估计定向多效性(directional pleiotropy)的大小。
3. Beavis效应：基于GWAS数据的MR研究可能会高估了遗传和暴露之间的关联，亦被称之为“胜利者的诅咒(the winner’s curse)
4. IVW中文叫做逆方差加权法，它的特点是回归时不考虑截距项的存在并且用结局方差（se的二次方）的倒数作为权重来进行拟合。与MR-eggr相比多了一个模型中的-1表示的就是去除截距项。
接下来我们就看下在R语言中实现MR的包TwoSampleMR。其主要基于两样本的设计策略，实现了IVW和MR-eggr的回归模型。我们首先看下包的安装：
devtools::install_github("MRCIEU/TwoSampleMR")然后我们直接基于实例来看下MR的实现：
## https://gwas.mrcieu.ac.uk/数据库中的实例
library(TwoSampleMR)
# 列出数据API内容
ao <- available_outcomes()
# 获得工具G的数据
exposure_dat <-extract_instruments("ieu-a-2")
# 获得G的影响的暴露因素
outcome_dat <-extract_outcome_data(snps=exposure_dat$SNP, outcomes="ieu-a-7")
# 连接暴露因素和结局因素
dat <- harmonise_data(exposure_dat,outcome_dat)
#查看MR的模型列表
mr_method_list()
# 获得结果，method_list=c("mr_egger_regression", "mr_ivw")
res <- mr(dat)
从上面我们可以看出所有的模型显示我们的G并不能通过BMI影响到心脏病。
mr_scatter_plot(res, dat)
从图中我们可以看出我们各模型的估计值均徘徊在原始的估计值周围并没有比原始估计值大。
##对比个模型和单个SNP估计值之间差异
res_single <- mr_singlesnp(dat)
p2 <- mr_forest_plot(res_single)
上图其实为森林图，我们可以看出和原始估计值具有很高的一致性。基本都是对结局有一定的影响
#总体评估的的森林图
res_loo <- mr_leaveoneout(dat)
p3 <- mr_leaveoneout_plot(res_loo)
#评估各模型之间的评估结果
res_single <- mr_singlesnp(dat)
p4 <- mr_funnel_plot(res_single)
##一对多的模型构建，目的对多个暴露同时进行分析
exp_dat <-extract_instruments(outcomes=c(2,100,1032,104,1,72,999))
table(exp_dat$exposure)
chd_out_dat <- extract_outcome_data(
snps= exp_dat$SNP,
outcomes= 7
)
dat2 <- harmonise_data(
exposure_dat= exp_dat,
outcome_dat= chd_out_dat
)
res<-mr(dat2)###可视化结果
res<-subset_on_method(res) #default isto subset on either the IVW method (>1 instrumental SNP) or Wald ratiomethod (1 instrumental SNP).
res<-sort_1_to_many(res,b="b",sort_action=4)#this sorts results by decreasing effect size (largest effect at top of theplot)
res<-split_exposure(res) # to keep the Yaxis label clean we exclude the exposure ID labels from the exposure column
res$weight<-1/res$se
min(exp(res$b-1.96*res$se)) #identify valuefor 'lo' in forest_plot_1_to_many
max(exp(res$b+1.96*res$se)) #identify valuefor 'up' in forest_plot_1_to_many
forest_plot_1_to_many(res,b="b",se="se",
exponentiate=T,ao_slc=F,lo=0.3,up=2.5,
TraitM="exposure",col1_width=2,by=NULL,
trans="log2",xlab="ORfor CHD per SD increase in risk factor (95% confidenceinterval)",weight="weight")
##为数据增加更多的信息
res$pval<-formatC(res$pval, format ="e", digits = 2)
forest_plot_1_to_many(res,b="b",se="se",
exponentiate=T,ao_slc=F,lo=0.3,up=2.5,
TraitM="exposure",by=NULL,
trans="log2",xlab="OR for CHD per SD increase in riskfactor (95% CI)",
weight="weight",subheading_size=11,
col1_title="Risk factor",
col1_width=2.5,
col_text_size=4,
addcols=c("nsnp","pval"),
addcol_widths=c(1.0,1.0),
addcol_titles=c("No. SNPs","P-val")
)
##对每个的暴露进行分组展示
res<-mr(dat2)
res<-split_exposure(res) # to keep the Yaxis label clean we exclude the exposure ID labels from the exposure column
res<-sort_1_to_many(res,group="exposure",sort_action=3,priority="Inversevariance weighted",trait_m="method")
forest_plot_1_to_many(res,b="b",se="se",
exponentiate=T,trans="log2",ao_slc=F,lo=0.03,
up=22,col1_width=2,by="exposure",TraitM="method",
xlab="ORfor CHD per SD increase in risk factor (95% confidence interval)",
subheading_size=12,col_text_size=4)
###获取GWAScatalog中的数据
library(MRInstruments)
data(gwas_catalog)##获取暴露数据
bmi_gwas <- subset(gwas_catalog,grepl("Speliotes", Author) & Phenotype == "Body massindex")
bmi_exp_dat <- format_data(bmi_gwas)#代谢数据获取
data(metab_qtls)
ala_exp_dat <-format_metab_qtls(subset(metab_qtls, phenotype=="Ala"))###从千人组计划获取LD值
bmi_exp_dat <- clump_data(bmi_exp_dat)###自定义的数据实例
bmi_file <-system.file("extdata", "bmi.txt",package="TwoSampleMR")
bmi_exp_dat <-read_exposure_data(bmi_file)
bmi2_file <-system.file("extdata/bmi.csv", package="TwoSampleMR")
outcome_dat <- read_outcome_data(
snps= bmi_exp_dat$SNP,
filename= bmi2_file ,
sep= ",",
snp_col= "rsid",
beta_col= "effect",
se_col= "SE",
effect_allele_col= "a1",
other_allele_col= "a2",
eaf_col= "a1_freq",
pval_col= "p-value",
units_col= "Units",
gene_col= "Gene",
samplesize_col= "n"
)
dat <- harmonise_data(bmi_exp_dat,outcome_dat)
# 获得结果，method_list=c("mr_egger_regression", "mr_ivw")
res <- mr(dat)