Mendelian randomization(MR)

2023-12-19 21:42:24 浏览数 (1)

孟德尔随机化：根据孟德尔遗传规律，亲代的等位基因随机分配给子代，此过程相当于随机对照研究(RCT)的随机分组过程：不受混杂因素（社会地位、行为等）的影响；满足时间顺序合理性（遗传变异继承于父母，且保持不变）

基本思想：利用与暴露因素具有强相关的遗传变异作为工作变量，借此推断暴露因素与研究结局之间的因果效应。

一.基础知识

1.数据来源：

现有全基因组关联研究（GWAS）数据大部分来源于欧洲，少量来源于亚洲、非洲、混合人种（暴露与结局最好来自同一个地方，消除人群异质性）

2.两样本MR的一个关键要点：

暴露与结局数据不存在或存在少量可接受范围内的样本重叠（最好不用同一库的）

3.下载的GWAS数据通常含以下变量

Rs号（pos,chr可以转化为rs号）

effect_allele(A1;效应等位基因)

other_allele(A2)

beta（OR;效应值)：表示SNP对表型起到的作用（结局是连续型变量用beta，大于0则增加风险；结局是二分类变量用OR，OR大于1则增加风险）

se

P

EAF(Freq;MAF；效应等位基因频率)

注：beta与OR可以相互转化：beta <- ln(OR)

P、beta、se知道其中两个可以计算另一个:se <- abs(beta/qnorm(p/2,lower.tail=F)

4.工具变量的获取

  • 从GWAS原文中获取
  • 从已发表的MR文章中获取（注意引用）
  • R program提取

5.核心假设（若筛选出的SNP较少，可忽略假设2和3；独立性假设和排他性假设具有相似性）

  • 相关性假设（SNP与暴露强相关）：

p<5e-08(视条件而定);F>10代表强相关; 注：p最低可放宽到5e-06（可能存在弱工具变量） ； P值越小，F越大

Clumping:LD r2<0.001（r2越小，SNP越少）,kb=10000(视条件而定)#连锁不平衡r2小于0.001(越小SNP之间越独立)，kb(指连锁不平衡的区域长度，在遗传学上我们认为在染色体上距离很近的遗传位点通常是“捆绑”在一起遗传给后代的)

代码语言：text复制
R2 <- 2*(1-MAF)*MAF*（beta）2 #计算每个SNP的R2将其相加，带入下面公式
F=(R2/(1-R2))((n-k-1)/k)# n为sample size , k为最终的SNP个数
  • 独立性假设（SNP与混淆因素无关）
  • 排他性假设（SNP只能通过暴露对结局发生作用，而不能通过其他途径）：p>5e-5 或p>5e-8或P结局>P暴露

PhenoScanner (cam.ac.uk)#假设2和假设3（将假设1筛选出的SNP逐个输入该网站验证）

6.核心假设的违背

  • 弱工具变量：当遗传变异与暴露因素不具有强相关关系（P），或者遗传变异仅能解释小部分的表型变异时（F值），称为弱工具变量

解决方案：1.增加样本量 2.增加表型解释度：相对于单个遗传变异，多个遗传变异能解释更大比例的表型变异

  • 多效性：当遗传变异可通过“遗传变异-暴露因素-结局”以外的其他通路影响结局发生时，该遗传变异具有多效性，这种多效性可导致独立性和排他性假设不成立

解决方案：使用生物学功能明确的遗传变异作为工具变量

  • 连锁不平衡：基因组位置相近的遗传变异倾向于共同遗传

解决方案：使用生物学功能明确的遗传变异作为工具变量

  • 人群分层：遗传变异频率在不同遗传背景的人群间存在差异，导致遗传变异与结局之间出现虚假关联

解决方案：纳入具有相同遗传背景的人群开展遗传关联研究

7.

  • 逆方差加权法（IVW）：假设所有遗传变异SNPs都是有效的工具变量，总体偏差为零。但是即使只有一个基因变异是无效的，它也会有偏倚
  • 加权中位数法（weighted median）:即使多达一半的权重来自无效的工具，该估计量也是一致的。相比于IVW，该方法具有更小的I类错误，且是对MR-Egger的补充
  • MR-egger回归：斜率系数：估计因果效应；截距：检验基因多效性，当截距接近于0则认为遗传多效性影响较小

8.森林图中置信区间跨过0或1表示结果不显著（p>0.05）

6.暴露数据和结局数据的SNP对应的效应等位基因不一致：Harmonization

代码语言：text复制
dat <- harmonise_data(exposure_dat = RI_exp_dat_clumped,outcome_dat = outcome_dat)

7.一部分暴露数据的SNP在结局数据中找不到

8.代理SNP如何确定？

9.选择结局数据注意要点

  • SNP量至少四五百万

10.mr_keep=FALSE的观测分析时会被删掉

11.统计效能power值的计算

mRnd: Power calculations for Mendelian Randomization (cnsgenomics.com)

二.实操

1.实操第一步：暴露数据

将数据下载到桌面，打开RStudio，导入暴露数据exposure.txt：

代码语言：text复制
setwd("C:/Users/76325/Desktop")
exposure1 <- read.table("exposure.txt",header=T)

相关性设置：将exposure中P<5e-08的观测筛选出来

代码语言：text复制
exposure2 <- (exposure1,p<5e-08)#p为表中的变量名
  • 将数据导出，进行数据预处理（保留分析所需栏(SNP/bata/se/effect_allele/other_allele/eaf/p)，对每一栏进行命名）
代码语言：text复制
write.csv(exposure2,file="exposure_RI.csv")
  • 将数据重新读回R，名称为exposure_RI
  • 重新整理数据
代码语言：text复制
RI_exp_dat <- read_exposure_data(filename=exposure_RI,sep=",",snp_col="SNP",bata_col="bata",
se_col="se",effect_allele_col=""effect_allele,other_allele_col="other_allele",eaf_col=“eaf”,
pval_col="p")

独立性设置

代码语言：text复制
RI_exp_dat_clumped <- clump_data(RI_exp_dat,clump_kb= 10000,clump_r2=0.01,clump_p1=1,
clump_p2=1,pop="EUR")#EUR表示样本来自欧洲
  • 将数据导出，计算F值，大于10的数据留下

2..实操第二步：结局数据

  • 将结局数据导入R，命名为PCOS_outcome
  • 将暴露与结局数据进行整合
代码语言：text复制
exposure_outcome <- merge(RI_exp_dat_clumped,PCOS_outcome,by.x="SNP",by.y="MarkerName")
  • 导出整合数据，删除exposure部分，剩余outcome部分，名称为outcome.csv
  • 将数据重新读回R
  • 重新整理数据
代码语言：javascript复制
outcome_dat <- read_outcome_data(snps=RI_exp_dat_clumped,filename=outcome.csv,sep=",",snp_col="SNP",bata_col="bata",
se_col="se",effect_allele_col=""effect_allele,other_allele_col="other_allele",eaf_col=“eaf”,
pval_col="p")

3、实操第三步：Harmonization（协同）

代码语言：text复制
dat <- harmonise_data(exposure_dat = RI_exp_dat_clumped,outcome_dat = outcome_dat)

4、实操第四步： MR estimation

代码语言：text复制
> mr(dat)#默认用五种方法分析
mr_method_list()#查看总共有多少种方法
#mr(dat,method_list=c("mr_ivw","mr_raps"))#指定方法分析
Analysing 'ieu-a-2' on 'ieu-a-7'
  id.exposure id.outcome                              outcome
1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
3     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
4     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
5     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
                       exposure                    method nsnp         b         se
1 Body mass index || id:ieu-a-2                  MR Egger   79 0.5024935 0.14396056
2 Body mass index || id:ieu-a-2           Weighted median   79 0.3870065 0.07201409
3 Body mass index || id:ieu-a-2 Inverse variance weighted   9 79 0.4459091 0.05898302
4 Body mass index || id:ieu-a-2               Simple mode   79 0.3401554 0.15562464
5 Body mass index || id:ieu-a-2             Weighted mode   70.3888249 0.09681892
          pval
1 8.012590e-04 #在MR-Egger的intercept分析时，P必须大于0.05，否则IVW结果不可靠
2 7.699254e-08
3 4.032020e-14  #IVW的P值<0.05，其他的大于0.05也可以当作阳性结果，但其他方法的beta值与IVW方向必须一致
4 3.183437e-02   #若方向不一致，将p值逐渐缩小 5e-9
5 1.352008e-04
  • 将beta值转换为OR值
代码语言：text复制
generate_odds_ratios(mr_res=mr(dat))
  • 数据可视化：散点图
代码语言：text复制
mr_scatter_plot(mr_results = mr(dat,method_list = c("mr_ivw","mr_egger_regression",
"mr_weighted_median")),dat)

5.实操第五步：敏感性分析（我们希望P>0.05）

5.1 异质性检测

代码语言：text复制
> mr_heterogeneity(dat)#Q值小于0.05说明存在异质性
  id.exposure id.outcome                              outcome
1     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
2     ieu-a-2    ieu-a-7 Coronary heart disease || id:ieu-a-7
                       exposure                    method        Q Q_df       Q_pval
1 Body mass index || id:ieu-a-2                  MR Egger 143.3046   77 6.841585e-06
2 Body mass index || id:ieu-a-2 Inverse variance weighted 143.6508   78 8.728420e-06
  • 若有异质性（有离群值），则
代码语言：text复制
run_mr_presso(dat,NbDistribution = 5000)#数字越大，精度越高
  • 异质性可视化（漏斗图：SNP是否对称分布在IVW线两侧）
代码语言：text复制
mr_funnel_plot(singlesnp_results =mr_singlesnp(dat))

5.2 多效性检测（SNP通过其他暴露因素影响到结局）：若存在多效性，则通过看文献了解存在哪些暴露因素，用phenoscanner网站剔除与暴露因素有关的SNP

代码语言：text复制
mr_pleiotropy_test(dat)#若p<0.05说明存在多效性（有截距），则不做了

5.3 Leave-one-out analysis（对结果有显著影响的SNP要剔除）

代码语言：text复制
mr_leaveoneout(dat)
  • 可视化（置信区间：beta /-1.96se）
代码语言：text复制
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))

三、在线数据库：IEU

代码语言：text复制
library(devtools)
devtools::install_github("MRCIEU/TwoSampleMR")
library(TwoSampleMR)
exposure <- extract_instruments(outcomes = "ieu-a-2",p1=5e-08)
outcome <- extract_outcome_data(snps = exposure$SNP,outcomes = "ieu-a-7")
dat <- harmonise_data(exposure_dat = exposure,outcome_dat = outcome)
mr(dat)

代码语言：text复制
#异质性检测
mr_heterogeneity(dat)#Q值大于0.05说明不存在异质性
run_mr_presso(dat,NbDistribution = 5000)

#异质性可视化
mr_funnel_plot(singlesnp_results =mr_singlesnp(dat))

#多效性检测
mr_pleiotropy_test(dat)#P<0.05则别做了

#leaveoneout analysis
mr_leaveoneout_plot(leaveoneout_results = mr_leaveoneout(dat))

#可视化

代码语言：text复制
mr_scatter_plot(mr_results = mr(dat,method_list = c("mr_ivw","mr_egger_regression","mr_weighted_median")),dat)

四、写作技巧

1.Introduction（注意暴露与结局应有侧重点）

background(first paragraph)

  • Epidemiology(发病率高),diagnosis（诊断困难，发现即晚期）,treatment（覆盖率低因此预防很重要）,prognosis（预后不好）,disease burden...Then taken together,point out the importance and urgency for disease prevention ; Previous studies on the association between exposure and outcome（general to specific,experimental and clinical）

Gap(second paragraph)

  • What is unknown based on the existing evidence?(Causality)
  • The inherent defects for the conventional observational studies.
  • Why filling gap(undetermined causality) is important?
  • How will we fill the gap?(introduce MR)

Strategy(third paragraph)

  • Introduce the basic principles of MR design
  • Introduce previous topic-related MR studies

Conclusion(third/fourth paragraph)

  • Briefly describe what you will do in current study,and what implications will have

2.Method

study overview:

Briefly describe the study design and the corresponding fundamental assumption

Data source:

Exposure & Outcome

Genetic variant selection criteria:

Association threshold,LD cutoff,exclude SNP associated with outcome,Harmonization to exclude palindromic and incompatible SNPs, F>10 ,and so on

statistical analysis:

MR estimations like IVW,weighted median, MR-Egger,RAPS,and so on; Sensitivity analysis like Cochran Q test（异质性检测）,Egger-intercept test（多效性检测）, leave-one-out analysis

Risk factor analysis

Power analysis（推荐达到80%）

3.Result（呈现方法的结果）

4.Discussion

  • findings(开门见山，发现了什么)
  • How finding fit in with existing literature
  • Interpretation（mechanism）:干预exposure降低outcome的风险；
  • strengths
  • limitations
  • Conclusion(implications of findings)

5.Abstract(concise:250~350 words)

  • Background(1~2 sentences)
  • Method
  • Results(show the most important results)
  • Conclusions(1~2 sentences)
r语言

0 人点赞