孟德尔随机化(Mendelian randomization，MR)是一种用于估计暴露因素(如生活方式或生物标志物)对结果/结局变量(如疾病或健康结果)的因果效应的统计方法。
注意一定不要搞颠倒了暴露因素和结局变量。
单变量MR:「目的」:评估一个特定暴露与一个特定结果之间的因果关系。
[如何工作」:使用一个或多个SNP作为工具变量来代表特定的暴露。
[应用场景]:当我们对一个明确的暴露和结果的因果关系感兴趣时。
[例子」:探讨BMI(暴露)是否会影响乳腺癌的风险(结果)。使用与BMI关联的SNP作为工具变量来评估BMI增加是否会增加乳腺癌的风险。
本实验设计的目的：分析BMI和乳腺癌之间是否有关系？（单变量MR分析）
暴露因素：BMI值
结局变量：患乳腺癌概率
实验结论：BMI和乳腺癌之间是有弱的负向的因果关系。
# 进行孟德尔随机化分析
###########################安装分析需要的R包##############
#install.packages("remotes")
library(remotes)
#remotes::install_github("MRCIEU/TwoSampleMR")
library("TwoSampleMR")
library(devtools)
#install_github("phenoscanner/phenoscanner") #读取访问PhenoScanner数据库
library(phenoscanner)
#remotes::install_github("MRCIEU/MRInstruments")
library(MRInstruments) ##该软件包包含许多 data.frames，其中每个都是 SNP 与性状关联的存储库。
library(ieugwasr) #一定要有这句话，否则会报错Error in r$status_code : $ operator is invalid for atomic vectors
library(tidyverse)
######获取暴露数据##############
##获取Speliotes研究的 BMI数据##此处使用BMI作为暴露因素
bmi_gwas <- subset(gwas_catalog,grepl("Speliotes", Author) & Phenotype == "Body mass index")
bmi_exp_dat <- format_data(bmi_gwas)
head(bmi_exp_dat)
#SNP chr.exposure gene.exposure effect_allele.exposure other_allele.exposure beta.exposure se.exposure pval.exposure units.exposure eaf.exposure exposure
#1 rs10150332 14 NRXN3 C T 0.13 0.03061224 3e-11 kg/m2 increase 0.21 Body mass index (kg/m2 increase)
#2 rs10767664 11 BDNF A T 0.19 0.03061224 5e-26 kg/m2 increase 0.78 Body mass index (kg/m2 increase)
##通过IEU gwas数据库找到了感兴趣的ID
bmi2014_exp_dat <- extract_instruments(outcomes = 'ieu-a-2')
#extract_instruments函数的参数
#p1 显著的snp的阈值，默认5e-08
#clump 是否对结果进行聚类,默认TRUE
#r2 SNP的R2过滤阈值，此处的值是允许的最大的R2. 默认0.001
#kb 计算LD R2的距离。默认是10000
##聚类分析
bmi_exp_dat <- clump_data(bmi2014_exp_dat)
#结局数据获取
#获取outcome数据，此处使用乳腺癌数据
ao <- available_outcomes()
breasr_cancer_list <- subset(ao,grepl("breast cancer",trait)) %>% arrange(year,desc(sample_size)) #也可以获得id
breasr_cancer_list
BMI_list <- subset(ao,grepl("BMI",trait)) %>% arrange(year,sample_size) #也可以获得id
BMI_list
#此处的编号是从 搜索breast cancer查询得到的结果
eur_out_dat <- extract_outcome_data(bmi_exp_dat$SNP,outcomes = "ukb-b-16890") #欧洲人的数据
##
#harmonise 协调暴露和结局数据的等位基因平衡
dat <- harmonise_data(
exposure_dat = bmi_exp_dat,
outcome_dat = eur_out_dat
)
#MR分析
#MR分析前，一定要有上面的协调过程，上面的结果可以得到暴露和结果性状中每个工具SNP的效应和标准误差。利用这些信息
#可以进行孟德尔随机化分析
res <- mr(dat)
res
###结果如下
mr_method_list()#查看mr支持的分析方法
#异质性分析
mr(dat,method_list = c("mr_egger_regression","mr_ivw")) #选择自己想要的方法
#上面的结果的Pval<0.05，说明研究中存在异质性。如果存在异质性，则需要想办法去调整样本。
#敏感性分析
mr_heterogeneity(dat) #如果p<0.05,同样说明敏感性很强，需要对样本进行处理。
#水平多效性
mr_pleiotropy_test(dat)
#Pval<0.05,则具有水平多效性显著。这个值如果显著的话，则可能做的MR就无法使用，需要更换暴露因素或结局，或者是调整样本。
#单SNP分析：分析每个SNP单独对结局变量的影响
#返回一个与mr()输出结果类似的data.frame,但它会对每个暴露-结果组合进行多次分析，每次使用不同的单SNP进行分析
res_single <- mr_singlesnp(dat)
#Leave-one-out analysis 留一法森林图
res_loo <- mr_leaveoneout(dat)
p3 <- mr_leaveoneout_plot(res_loo)
p3
#作图
##散点图
res <- mr(dat)
p1 <- mr_scatter_plot(res,dat)
p1
##森林图
res_single <- mr_singlesnp(dat)
p2 <- mr_forest_plot(res_single)
p2
#漏斗图
p4 <- mr_funnel_plot(res_single)
p4
散点图的解释
散点图中的每个点代表一个遗传变异，它向我们展示了每个遗传变异(SNP)与暴露(Exposure)和结果(Outcome)的关联。虚线表示使用不同模型拟合的线段，其中每个遗传(例如SNP)与暴露的关联可以直接预测其与结果的关联
该图中使用了IVW和MT Egger模式，可以看到斜率是负数，说明BMI指数和乳腺癌是负相关的关系，斜率绝对值越大，说明相关性越高。斜率大于0，代表暴露因素是结局的不利因素。斜率小于0，说明BMI是乳腺癌的有利因素。即BMI越大，越容易患乳腺癌。
一般以IVW的结果的Pval<0.05为准，说明显著。
森林图的解释
「标准的森林图」为每个SNP展示了效应大小和95%置信区间。从图中，我们可以看到每个SNP与结果之间的关联强度以及其不确定性
上图中可以看到rs10132280对整个的负向关联强度最大。
底部的红色的MR Egger和IVW也可以看到是负相关。
「留一法所指的森林图」，展示了在移除每个单独的SNP后的MR估计结果。每次移除一个SNP，都会进行MR估计，并将结果及其置信区间记录在图中。如果剔除了某一个SNP后，结果改变很大，说明存在某一个SNP对结果影响很大，这将大大降低分析结果的可靠性
从上图中可以看出，第一个snp rs1421085对结果的影响很大，应该剔除。
该方法是也是对敏感性显著的样本的处理方法
如果敏感性很强，或者异质性很强，这时候就需要去剔除离群的样本，通过上面的留一法森林图可以看出每个snp单独对结果的贡献率，离群值要丢掉，否则会影响MR分析的可靠性。
漏斗图
漏斗图用于检查各个遗传变异是否存在异质性，当没有异质性时，漏斗图呈现出对称的形状这意味着研究效应和它们的精确度之间没有系统性的关系(在进行孟德尔随机化分析时，异质性可能表明所选的遗传变异不仅仅与暴露有关，还可能与其他潜在的混淆因素有关，这可能会影响因果关系的解释)
MR分析步骤
- 读取exposure暴露数据，并过滤
- 获取outcome结局数据。
- harmonise 协调暴露和结局数据的等位基因平衡
- 孟德尔随机化分析MR
- 异质性分析
- 多效性分析
- 留一法敏感性分析
关键词解释
敏感性分析:统计学和模型分析中的一个关键概念，主要用于评估模型、估计或结果的稳定性和可靠性。在孟德尔随机化(MR)分析中，敏感性分析尤为重要，因为MR的效力和可靠性可能会受到多种偏见和混淆因素的影响。
MR中常见的敏感性分析有:
留一法(Leave-One-Out)分析(类似于留一交叉验证)
留一法 (Leave-One-Out Analysis)」:留一法用于评估单个SNP是否对整体MR估计产生了过度的影响。具体操作是一次移除一个SNP，然后使用其余的SNP重新进行MR估计。这一过程对每个SNP都进行一次。如果移除某个SNP后，整体估计发生显著变化，这可能意味着该SNP对MR估计有过度的影响或可能存在偏见。异质性测试
mr_heterogeneity()水平多效性检验（）
要直接从GitHub存储库安装TwoSampleMR，请执行以下操作:
library(remotes)
install_github(’MRCIEU/TwoSampleMR’)
TwoSampleMR提供了以下一些函数来执行相关分析:
available exposures() 查看可用的暴露数据，"from the lEU GWAS database thatrooted in the package”，下同。
extract instruments( 来提取特定暴露的数据。
clump_data()去除连锁不平衡的SNP，从而得到一组代表性的、不相关的SNP
available_outcomes()查看可用的结果数据。
extract_outcome_data(snp=exposure_dat$SNP,outcomes = "ukb-b-13348") 来提取特定结果的数据
harmonise_data(暴露变量，结局变量) 来协调暴露和结果数据。这一步确保暴露和结果数据在相同的基因位点上，并且效应大小都是针对相同的参考等位基因。mr0来执行MR分析，提供多种方法，如IVW、MR-Egger等: 其实就是统一正负链的一致性，让相同位置的SNP的参考碱基和主要变异的碱基统一。
mr_heterogeneity(dat)来进行异质性检验。
mr_pleiotropy_test(dat) 来进行水平多效性检验。
mr_leaveoneout()留一法检验，也被称为逐个剔除检验。
MR相关的概念解读来源于 http://www.360doc.com/content/23/1115/21/72552828_1104165784.shtml
TwoSampleMR官方文档 https://mrcieu.github.io/TwoSampleMR/articles/perform_mr.html
入门的关键
- 是要学会不同的数据读取方式，因为很多时候无法访问外网，所以需要把数据从 https://gwas.mrcieu.ac.uk/下载到本地之后，再读取处理。特别是下载原始的vcf的时候，本地电脑可能就处理不了这么大的文件。
- 学会理解输出的各种结果的意义
当我们计算出MR的结果后，接下来就要进行敏感性分析，这里我们主要从如下三方面进行检验:
(1)异质性检验(Heterogeneity test):主要是检验各个IV之间的差异，如果不同IV之间的差异大，那么这些I的异质性就大。
(2)多效性检验(Pleiotropy test):主要检验多个IV是否存在水平多效性，常用MR Egger法的截距项表示，如果该截距项与0差异很大，说明存在水平多效性。这个从散点图的散点在纵坐标的方向上分布可以看出。
(3)逐个剔除检验(Leave-one-out sensitivity test):主要是逐个剔除IV后计算剩下IV的MR结果，如果剔除某个IV后其它IV估计出来的MR结果和总结果差异很大，说明MR结果对该IV是敏感的。
异质性检验(Q pval远小于0.05具有异质性)，可剔除某些outcome的P值非常小的SNP，或可使用随机效应模型观察结局校正后是否还有显著性。
MR选择方法的三大原则
- 在没有异质性和多效性的时候，优先使用IVW模型。
- 只有异质性，没有多效性时，优先使用Weighted Median方法的结果(也可以使用IVW的随机效应模型)
- 有多效性时，优先使用MR-Egger方法计算出的结果。
最后提醒一下，SNP数目多的时候统计效力足，但是异质性和多效性可能会比较大，如果去除部分SNP后，可以消除异质性和多效性，但是会导致统计效力低，使结果变成阴性。因此，大家需要好好斟酌一下选择何种方法!
孟德尔分析词汇表查询 https://mr-dictionary.mrcieu.ac.uk/
