大家好,我是邓飞。
孟德尔随机化和GWAS一样,是比较常见的分析方法,这里用R语言包TwoSampleMR进行介绍,如何走通MR分析流程。
孟德尔随机化定义:
孟德尔随机化分析(Mendelian Randomization, MR)是一种流行的流行病学研究方法,主要用于探讨暴露因素(如生活习惯、环境因素等)与健康结果(如疾病、健康状况等)之间的因果关系。这种方法利用了遗传变异(通常是单核苷酸多态性,SNP)作为工具变量,以避免混杂因素和逆因果关系的影响。
孟德尔随机化和孟德尔的关系:
MR遵循“亲代等位基因随机分配给子代”的孟德尔遗传定律,选择合适的“基因变异”作为工具变量,指代无法测量的待研究暴露因素,通过测量遗传变异与暴露因素、遗传变异与疾病结局之间的关联,进而推断暴露因素与疾病结局之间的关联。
孟德尔随机化的由来:
1986年,Katan首次提出MR的遗传思想:由于配子形成时,遵循“亲代等位基因随机分配给子代”的孟德尔遗传规律,如果基因型决定表型,基因型通过表型而与疾病发生关联,因此可以使用基因型作为工具变量来推断表型与疾病的关联。
一、数据来源
示例数据使用官网的数据,进行了一点补充,对结果进行了可视化。(https://mrcieu.github.io/TwoSampleMR/articles/introduction.html)
整个步骤:
步骤1:提取暴露数据的GWAS
> ## 1, 安装TwoSampleMR,如果已安装,可以忽略
>
> # library(remotes)
> # install_github("MRCIEU/TwoSampleMR")
>
> ## 2, 载入TwoSampleMR包
> library(TwoSampleMR)
>
> ## 3,从数据库中提取暴露的GWAS summary数据
> exposure_dat = extract_instruments("ieu-a-2")
> dim(exposure_dat)
[1] 79 15
共有79行15列的暴露数据结果。
如果不会设置token,可以参考这篇博文:孟德尔随机化R包TwoSampleMR安装教程并设置token
步骤2:提取结局数据的GWAS
> ## 4,从数据库中提取结局变量的的GWAS summary数据,SNP用暴露数据的结果
> # Get effects of instruments on outcome
> outcome_dat = extract_outcome_data(snps=exposure_dat$SNP, outcomes = "ieu-a-7")
Extracting data for 79 SNP(s) from 1 GWAS(s)
> dim(outcome_dat)
[1] 79 16
共79行15列的结局数据,注意,这里直接使用暴露数据质控后的SNP,提取结局数据得到的结果,所以位点数是一样的。
步骤3:合并暴露数据和结局数据
> ## 5,将暴露数据和结局数据合并
> # Harmonise the exposure and outcome data
> dat = harmonise_data(exposure_dat, outcome_dat)
Harmonising Body mass index || id:ieu-a-2 (ieu-a-2) and Coronary heart disease || id:ieu-a-7 (ieu-a-7)
> dim(dat)
[1] 79 36
合并的数据,共79行,36列,这些数据可以用于孟德尔随机化的分析。
步骤4:孟德尔随机化分析及结果可视化
> ## 6,进行孟德尔随机化分析
> res = mr(dat)
Analysing 'ieu-a-2' on 'ieu-a-7'
> ## 7,异质化分析
> mr_heterogeneity(dat)
id.exposure id.outcome outcome exposure method Q
1 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7 Body mass index || id:ieu-a-2 MR Egger 143.3046
2 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7 Body mass index || id:ieu-a-2 Inverse variance weighted 143.6508
Q_df Q_pval
1 77 6.841585e-06
2 78 8.728420e-06
> ## 8,水平多效性分析
> mr_pleiotropy_test(dat)
id.exposure id.outcome outcome exposure egger_intercept se pval
1 ieu-a-2 ieu-a-7 Coronary heart disease || id:ieu-a-7 Body mass index || id:ieu-a-2 -0.001719304 0.003985962 0.6674266
> ## 9,留一法分析
> res_loo = mr_leaveoneout(dat)
> mr_leaveoneout_plot(res_loo)
$`ieu-a-2.ieu-a-7`
attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
id.exposure id.outcome
1 ieu-a-2 ieu-a-7
> ## 10,散点图
>
> p1 = mr_scatter_plot(res, dat)
> p1
$`ieu-a-2.ieu-a-7`
attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
id.exposure id.outcome
1 ieu-a-2 ieu-a-7
>
> ## 11,森林图
> res_single = mr_singlesnp(dat)
> mr_forest_plot(res_single)
$`ieu-a-2.ieu-a-7`
attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
id.exposure id.outcome
1 ieu-a-2 ieu-a-7
Warning messages:
1: Removed 1 row containing missing values or values outside the scale range (`geom_errorbarh()`).
2: Removed 1 row containing missing values or values outside the scale range (`geom_point()`).
>
> ## 12,漏斗图
> mr_funnel_plot(res_single)
$`ieu-a-2.ieu-a-7`
attr(,"split_type")
[1] "data.frame"
attr(,"split_labels")
id.exposure id.outcome
1 ieu-a-2 ieu-a-7
留一法的森林图:
孟德尔随机化的森林图:
孟德尔随机化的散点图:
孟德尔随机化的漏斗图:
还有哪些需要研究的?
如何读取自己的GWAS summary结果,并将格式整理为TwoSampleMR的格式?
如何对暴露数据GWAS结果进行质控,包括LD质控,F值质控,R2质控等?
怎么对已发表的文章进行结果图标的复现?
这些都是细枝末节,等我后续一一完成博客的文章,欢迎继续关注。
上面分析完整的代码汇总:
## 1, 安装TwoSampleMR,如果已安装,可以忽略
# library(remotes)
# install_github("MRCIEU/TwoSampleMR")
## 2, 载入TwoSampleMR包
library(TwoSampleMR)
## 3,从数据库中提取暴露的GWAS summary数据
exposure_dat = extract_instruments("ieu-a-2")
dim(exposure_dat)
## 4,从数据库中提取结局变量的的GWAS summary数据,SNP用暴露数据的结果
# Get effects of instruments on outcome
outcome_dat = extract_outcome_data(snps=exposure_dat$SNP, outcomes = "ieu-a-7")
dim(outcome_dat)
## 5,将暴露数据和结局数据合并
# Harmonise the exposure and outcome data
dat = harmonise_data(exposure_dat, outcome_dat)
dim(dat)
## 6,进行孟德尔随机化分析
res = mr(dat)
## 7,异质化分析
mr_heterogeneity(dat)
## 8,水平多效性分析
mr_pleiotropy_test(dat)
## 9,留一法分析
res_loo = mr_leaveoneout(dat)
mr_leaveoneout_plot(res_loo)
## 10,散点图
p1 = mr_scatter_plot(res, dat)
p1
## 11,森林图
res_single = mr_singlesnp(dat)
mr_forest_plot(res_single)
## 12,漏斗图
mr_funnel_plot(res_single)
想要更好的学习和交流,快来加入飞哥的知识星球,这是一个生物统计+数量遗传学+GWAS+GS的社区,在这里你可以向飞哥提问、帮你制定学习计划、跟着飞哥一起做实战项目,冲冲冲。点击这里加入吧:飞哥的学习圈子