看完不会来揍我 | 孟德尔随机化（二）—— 代码实操

最近真的是超超超超超超超级多的小伙伴们在咨询孟德尔随机化相关的问题和课程，意想不到的那种多！那我怎么办嘞！整呗！主打的就是一个宠粉！

关于孟德尔随机化，我们之前就已经在孟德尔随机化（一）| 随处可见的孟德尔随机化到底是什么？中介绍过啦！对孟德尔随机化的理解还比较模糊的小伙伴可以去看看，说不定会有帮助哟！实操部分一直说要更新的，结果因为最近在忙毕业的事情，所以就拖到现在啦！抱歉抱歉啦！现在它终于要出场啦！咱们今天，就尽力给大家把孟德尔随机化讲透咯！让我们每个人都能够亲手做出自己的孟德尔随机化！

如果小伙伴们有需求的话，可以加入我们的交流群：一定要知道 | 永久免费的环境友好型生信学习交流群又双叒叕来啦！| 伴随不定期群友好物分享！在这里，你可以稍有克制地畅所欲言！
超级建议大家在入群前或入群后可以看一下这个：干货满满 | 给生信小白的入门小建议 | 掏心掏肺版！绝对干货满满！让你不虚此看！
如果有需要个性化定制分析服务的小伙伴，可以看看这里：你要的个性化生信分析服务今天正式开启啦！定制你的专属解决方案！全程1v1答疑！！绝对包你满意！

咱们直接开始边实操边讲解！冲冲冲！！！

分析流程

我们先大致了解一下进行**孟德尔随机化（Mendelian Randomization，MR）**的工作流程（里面有些词不理解没关系，咱们下面会详细介绍），如下所示：

为暴露选择工具变量（如有需要进行连锁不平衡（Linkage Disequilibrium，LD）clumping）；
从 IEU（https://gwas.mrcieu.ac.uk/）GWAS 数据库提取与感兴趣的结局相关的工具变量（除了这个数据库，还有其他数据库也可以，但只有这个数据库的数据可以直接用TwoSampleMR包获取，从其他数据库或者相关文献拿到的数据需要我们对其进行相应处理，咱们后面都会介绍的，大家不要担心！）；
协调或整合暴露数据和结局数据；
进行 MR 分析、敏感性分析、绘制图表、报告撰写等等等等。

下面是官网的示意图：

上面是对 MR 分析整体的概览，咱们现在开始掰开了揉碎了介绍喽！

工具变量的获取

从 GWAS 相关文献获取

示例文献：Novel loci affecting iron homeostasis and their effects in individuals at risk for hemochromatosis | Nature Communications（https://www.nature.com/articles/ncomms5926）

像这种得到的显著 SNP 位点不多的研究，它会把数据放在正文中，我们直接拿来用就好啦！别忘了引用！但是还需要注意一个问题，有些文章中展示的位点，虽然是显著的，但它们不一定相互独立（咱 MR 分析的基本条件嘛），如果这样的话，我们就需要再给它去除一下连锁不平衡。

从 MR 相关文献中获取

示例文献：Gout and the risk of Alzheimer′s disease: A Mendelian randomization study - Lee - 2019 - International Journal of Rheumatic Diseases - Wiley Online Library（https://onlinelibrary.wiley.com/doi/10.1111/1756-185X.13548）

不过，感觉它们多数也是从 IEU（https://gwas.mrcieu.ac.uk/）GWAS 数据库中拿到的。个人建议哈，最好还是咱们自己从原始数据出发开始处理，毕竟，万一个别文章分析过程有问题呢 🌚

利用 TwoSampleMR 包直接获取

咱们今天就主要介绍这种！请看下文！从文献中获取的数据的处理方法我们后面也会详细介绍滴！

TwoSampleMR 包的安装

那我们开始安装包！

############################ 孟德尔随机化 ############################

# 首先安装包
devtools::install_github("MRCIEU/TwoSampleMR")

# 或者也可以使用官方安装代码，其实都差不多的啦！
# 咱们在R包的下载与安装一文中也有详细介绍，有兴趣的小伙伴可以查看：https://mp.weixin.qq.com/s/jkCzmjTw_ILw4jTfaW7NpQ
# install.packages("remotes")
# remotes::install_github("MRCIEU/TwoSampleMR")

这里的包需要从 GitHub 进行下载安装，如果遇到 token 设置或者其他问题，大家可以查看：看完不会来揍我 | R包的下载与安装 | 再也没有一个包可以逃出你的手掌心啦，绝对保真！看完这篇，大概率是没有包能够逃出你的手掌心啦哈哈哈哈哈哈哈哈！

安装TwoSampleMR包的时候，它会帮助我们同时安装其他辅助包，有时候会提示我们是否需要更新其他包（提示可能不止一次，而且这个包安装需要一些时间，大家不要急，咱慢慢来！），咱们一般选3不更新就好啦！至于包是否需要经常更新，咱们主打一个能不更新就不更新！想了解详情的小伙伴也可以继续查看：看完不会来揍我 | R包的下载与安装 | 再也没有一个包可以逃出你的手掌心啦！

包安装好之后，我们就可以开始准备数据啦！

数据准备

使用 TwoSampleMR 包进行 MR 分析时，通常需要两个关键的数据集：一个用于估计遗传变异与暴露因子的关系（关联分析），另一个用于估计遗传变异与结局的关系（效应量估计）。接下来，咱们就来看看要怎么获取它们叭！

我有点好奇教育水平（educational attainment）和乳腺癌（breast cancer）之间是不是存在一些奇妙的关联，毕竟阿巴阿巴你懂的哈哈哈哈哈哈哈哈！（我就随脑子那么一想！大家随意发挥！）

注意注意：咱们这里先演示的是直接通过TwoSampleMR包下载数据并对其进行处理，如果是从文献中获取的数据，处理方式咱们后面演示好不好！

那么我们如何获取暴露因子（这里我们就指教育水平叭）和结局变量（这里指乳腺癌）的 SNP 数据呢？通常，这些数据可以从多个 GWAS 数据库或文献中获取。我们今天使用IEU 数据库（https://gwas.mrcieu.ac.uk/），这个数据库的优点是它允许我们直接通过TwoSampleMR包下载所需的数据，是不是超方便！

IEU数据库界面如下：

我们可以直接在中间的大框框里输入关键词，也可以点击右上角的datasets进入新的页面，在Trait contains的框框里输入关键词。比如我们这里就以educational attainment（教育水平）作为关键词进行输入，然后点击🔍或者Filter，就会看到有很多数据被筛选出来啦（虽然这里没有很多哈哈哈哈哈哈）！

这里我们选择ebi-a-GCST90029012这个数据。

点击上图的红色箭头，我们还可以进入数据详情界面，如下图所示：

接下来，我们再找乳腺癌（breast cancer）相关数据，以breast cancer为关键词进行输入，筛选得到下面的数据，我们选择ukb-b-16890这个数据进行后续分析。

所以嘞，咱们研究教育水平（educational attainment）对乳腺癌（breast cancer）的影响，就需要下载前面提到的ebi-a-GCST90029012和ukb-b-16890这两个数据，咱们用TwoSampleMR包进行下载。

获取暴露因子的 SNP 数据

# 加载包
library(TwoSampleMR)

# 下载暴露数据
exposure_data <- extract_instruments(
  outcomes = "ebi-a-GCST90029012",  # 指定 GWAS 数据集
  p1 = 5e-08,                       # 过滤 SNP 的 p 值阈值
  clump = TRUE,                     # 进行 SNP clumping
  r2 = 0.001,                       # clumping 的 r^2 阈值
  kb = 10000,                       # clumping 的窗口大小
  access_token = NULL               # 不使用访问令牌
)

# 其实后面成参数都是默认的，大家需要修改的时候再用即可（视情况调整），一般写成下面这样也是可以的！
# exposure_data <- extract_instruments(outcomes = "ebi-a-GCST90029012")

下面我们挨个介绍一下这几个参数，理解这些可以帮助大家更好地将其应用于自己的数据：

extract_instruments: 用于从 GWAS 数据中提取工具变量。
outcomes = 'ebi-a-GCST90029012': 这里指定了你希望从哪个 GWAS 数据集中提取工具变量。我们这里的数据集标识符是ebi-a-GCST90029012。
p1 = 5e-08: 用于过滤 SNP 的 p 值阈值，仅保留 p 值小于或等于 5e-08 的 SNP。p 值设置很小是为了找到与暴露强相关的 SNP，具体数值可视情况调整。
clump = TRUE: 这个参数告诉函数是否进行 SNP clumping。SNP clumping 是一个常用的步骤，用于减少 SNP的数量并确保它们是独立的。
r2 = 0.001: 这是 clumping 的阈值，它指定了 SNP 之间的连锁不平衡（linkage disequilibrium，LD）的 r^2 值的最小阈值。在这里表示，仅当 SNP 之间的 r^2 小于 0.001 时，它们才会被认为是独立的。
kb = 10000: 这是 clumping 的窗口大小，表示在进行 SNP clumping 时考虑的基因区域的大小。在这里，我们的窗口大小设为 10,000 kb。
access_token = NULL: 这个参数通常用于提供访问 GWAS 数据的访问令牌。这里我们设置为 NULL，表示不使用访问令牌。

这段代码从ebi-a-GCST90029012GWAS 数据集中提取工具变量，并进行 SNP clumping 以确保所提取的 SNP 是独立的。提取的工具变量数据将存储在 exposure_data 变量中，可用于后续的 MR 分析。

# 我们简单看一下数据长什么样子
head(exposure_data)[1:4, ]
#   chr.exposure pos.exposure beta.exposure se.exposure pval.exposure samplesize.exposure
# 1            1     41830086   -0.00759727 0.001156820   5.90065e-11              470941
# 2            1     44076469    0.00882349 0.001003490   3.19963e-19              470941
# 3            1     53740797   -0.00607082 0.000946858   1.70000e-10              470941
# 4            1    184668823   -0.00591639 0.000938774   2.10000e-10              470941
#          id.exposure        SNP effect_allele.exposure other_allele.exposure
# 1 ebi-a-GCST90029012 rs56044892                      T                     C
# 2 ebi-a-GCST90029012   rs549845                      A                     G
# 3 ebi-a-GCST90029012  rs6696068                      T                     G
# 4 ebi-a-GCST90029012 rs10911634                      G                     T
#   eaf.exposure                                                             exposure
# 1     0.211119 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 2     0.699305 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 3     0.383936 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 4     0.398718 Educational attainment (college completion) || id:ebi-a-GCST90029012
#   mr_keep.exposure pval_origin.exposure data_source.exposure
# 1             TRUE             reported                  igd
# 2             TRUE             reported                  igd
# 3             TRUE             reported                  igd
# 4             TRUE             reported                  igd

# 有点多不方便看哈，那咱们看看都有哪些列叭！顺便解释一下给大家！
colnames(exposure_data)
#  [1] "chr.exposure"           "pos.exposure"           "beta.exposure"         
#  [4] "se.exposure"            "pval.exposure"          "samplesize.exposure"   
#  [7] "id.exposure"            "SNP"                    "effect_allele.exposure"
# [10] "other_allele.exposure"  "eaf.exposure"           "exposure"              
# [13] "mr_keep.exposure"       "pval_origin.exposure"   "data_source.exposure"

咱们解释一下列名们的含义：

chr.exposure: 该列表示 SNP 所在的染色体。
pos.exposure: 表示 SNP 在染色体上的位置。
beta.exposure: 该列包含 SNP 与暴露因子之间的回归系数。这通常用于 MR 分析中作为工具变量的效应量。一定要有！

如果是分类变量的话，用 beta 值就不合适啦！我们会用到 OR 值，有的数据会给出，有的不会，那么就需要我们自行转换一下啦！使用generate_odds_ratios函数即可进行转换，懵懵哒不要慌！咱们在后面会进行演示的！
se.exposure: 标准误（standard error），用于估计 beta.exposure 的不确定性。一定要有！
pval.exposure: 表示 beta.exposure 是否显著不为 0 的 p 值，用于评估 SNP 与暴露因子之间关系的显著性。一定要有！
samplesize.exposure: 用于估计 beta.exposure 的样本大小。
id.exposure: GWAS 数据集的 GWAS ID。
SNP: SNP 的标识符或 ID。一定要有！
effect_allele.exposure: 效应等位基因（effect allele）。如果是从其他数据库或相关文献中拿到的数据，也可能写为alternative allele或者allele1。一定要有！
other_allele.exposure: 另一等位基因。如果是从其他数据库或相关文献中拿到的数据，也可能写为reference allele或者allele2。
eaf.exposure: 效应等位基因的频率（effect allele frequency）。
exposure: 暴露因子的名称或标识符。
mr_keep.exposure: 一个布尔值，表示是否在 MR 分析中保留该 SNP。
pval_origin.exposure: 暴露因子与 SNP 的原始关联分析的 p 值。
data_source.exposure: 提供 SNP 数据的数据源或来源。

获取结局变量的 SNP 数据

现在暴露数据导入完毕，我们开始获取结局数据。我们需要的是与上述经筛选的暴露数据中 SNP 匹配的结局数据，所以可以像下面一样这么干！

# 从指定 GWAS 数据集中提取与给定 SNP 关联的结局数据
outcome_data <- extract_outcome_data(
  snps = exposure_data$SNP,       # 指定要提取结局数据的 SNP 列表
  outcomes = 'ukb-b-16890',       # 指定 GWAS 数据集
  proxies = FALSE,                # 不使用代理 SNP
  maf_threshold = 0.01,           # 最小等位基因频率阈值
  access_token = NULL             # 不使用访问令牌
)

参数们的解释：

extract_outcome_data : 从指定的 GWAS 数据集中提取与给定 SNP 关联的结局数据。咱们这里是从ukb-b-16890GWAS 数据集中提取与 exposure_data$SNP 列表中的 SNP 关联的结局数据。
outcomes = 'ukb-b-16890’: 指定 GWAS 数据集，帮助我们确定结局变量是什么。在这里，ukb-b-16890就是我们要从中提取结局数据的 GWAS 数据集的标识符。也就是说，乳腺癌是我们的结局变量。
proxies = FALSE : 表示不使用代理 SNP（proxy SNP）来代替给定的 SNP，即只提取指定的 SNP 的结果（一般咱选不使用就好，除非数据不足，个人建议，仅供参考）。有时候暴露数据与结局数据的 SNP 可能是不一致的，也就是说有那么一部分暴露数据的 SNP 在结局数据中找不到，大家有时候就会考虑使用代理 SNP（我下面给大家简单介绍了一下代理 SNP，有兴趣的小伙伴们可以去瞅瞅！不小心又写长了，哎呀！控几不住我寄几）。
maf_threshold = 0.01 : 设置最小等位基因频率（MAF，minor allele frequency）阈值为 0.01，这意味着仅保留等位基因频率大于或等于 0.01 的结局数据。默认值是0.3，大样本 GWAS 可以适当调低。
access_token = NULL : 表示不使用访问令牌来访问 GWAS 数据，这与之前的代码段保持一致。

关于代理 SNP 的简单解释

代理 SNP（Proxy SNP）其实就是一个与我们实际关心的遗传变异（目标 SNP）高度关联的遗传变异。在 GWAS 或 MR 分析中，代理 SNP 通常用作目标 SNP 的替代来进行统计分析，因为有时目标 SNP 可能没有直接的测量数据或没有达到统计显著性。

为什么使用代理 SNP：

替代缺失数据：有时候，我们关心的目标 SNP 可能在某些数据集中没有测量数据或者没有足够的统计显著性。
增加数据的覆盖率：代理 SNP 可以提高数据的覆盖率，因为它们通常在更多的 GWAS 数据集中有可用的测量数据。
辅助分析：在进行复杂的 MR 分析或修复的环境相关性分析时，代理 SNP 可以作为辅助信息来增加统计效能。

如何确定代理 SNP：

确定代理 SNP 通常基于遗传变异与目标 SNP 之间的连锁不平衡（linkage disequilibrium, LD）关系。代理 SNP 通常是与目标 SNP 高度关联的遗传变异，这意味着当一个遗传变异发生变异时，另一个遗传变异也很可能会发生变异。以下是一些确定代理SNP的方法：

基于连锁不平衡: SNP 之间的 LD 通常用 r^2 值来衡量，r^2 值接近 1 表示两个 SNP 之间的高度关联。根据 r^2 值，我们可以选择与目标 SNP 高度关联（r^2 大于或等于某个预定阈值，例如 0.8 或 0.9，很多文章都会选择 0.8，更多文章选择不使用代理 SNP 🌚）的 SNP 作为代理 SNP。
考虑位置和功能: 除了 LD 关系外，还可以考虑代理 SNP 的位置（是否在目标基因或附近）和功能（是否已知为功能性变异或与疾病风险相关）。
验证代理 SNP: 在选择代理 SNP 后，建议进行进一步的验证，例如查看它是否在其他独立的数据集中与目标 SNP 高度关联，或者查看其与暴露和结局之间的关系。

注意事项：

LD 关系的不稳定性：由于不同人群之间的遗传结构差异，代理SNP在不同人群中的效果可能会有所不同。因此，最好在研究所考虑的目标人群中验证代理 SNP 的效果。
水果蔬菜偏误的可能性：使用代理 SNP 进行分析时，应注意水果蔬菜偏误（pleiotropy bias）的风险，并考虑使用修复的环境相关性分析或其他方法来评估和校正这种偏误。

是否必须使用代理 SNP：

不一定。在很多情况下，目标 SNP 有足够的测量数据和统计显著性，因此不需要使用代理 SNP。但在一些特定的分析场景中，代理 SNP 可能是有用的。

何时使用代理 SNP：

数据不足时：当目标 SNP 的数据不足或没有统计显著性时，可以考虑使用代理 SNP。
提高统计效能时：在进行复杂的 MR 分析或修复的环境相关性分析时，代理 SNP 可以作为辅助信息来提高统计效能。

何时不使用代理 SNP：

目标 SNP 有足够的数据和统计显著性：如果目标 SNP 有足够的测量数据和统计显著性，通常不需要使用代理 SNP。
分析目标明确：如果分析的目标是特定的目标 SNP，而不是它的代理，那么就不需要使用代理SNP。

统一效应等位与效应量（Harmonization）

将暴露数据和结局数据进行协调或整合，以确保它们在格式和结构上是一致的。这步必不可少！

# 协调或整合暴露数据和结局数据，以确保它们具有相同的格式和结构
data <- harmonise_data(
  exposure_dat = exposure_data,  # 指定暴露数据，这是之前从 GWAS 数据集中提取的工具变量数据
  outcome_dat = outcome_data      # 指定结局数据，这是从另一个 GWAS 数据集中提取的关联结局数据
# Harmonising Educational attainment (college completion) || id:ebi-a-GCST90029012 (ebi-a-GCST90029012) and Cancer code, self-reported: breast cancer || id:ukb-b-16890 (ukb-b-16890)
# Removing the following SNPs for being palindromic with intermediate allele frequencies:
#   rs10931821, rs10934957, rs12534506, rs1441098, rs1913145, rs2478208, rs2717609, rs4705755, rs62477728, rs7624274, rs7920624

# 这里可以看到它帮我们剔除掉了一些回文 SNP 和不兼容的 SNP，我们下面给大家细讲一下！

为什么要进行 Harmonization 这个步骤呢？

协同暴露数据SNP和结局数据SNP的等位基因方向：这是什么意思嘞！就是说呀，在不同的研究或数据源中，暴露和结局的 SNP 可能存在不一致的方向。就像下面这俩红绿框框中的 SNP，可以看到暴露数据与结局数据中 SNP 是不是一样的方向，太巧了！但是，有时候它们可能并不一致，刚好相反，有可能暴露是 G-A，结局是 A-G。那这种情况下我们该怎么办呢？Harmonization！这个步骤就可以帮助我们解决这个问题！
据 EFA 大小，剔除不能判断方向的回文 SNP：有些 SNP 的方向不明确，即它们与暴露或结局的关系不能确定。我们就可以通过判断效应等位基因频率（effect allele frequency，EFA）的大小，剔除这些无法判断方向的SNP。
剔除incompatible SNP：在数据协调或整合过程中，可能会出现不兼容的 SNP，即在暴露数据和结局数据中不一致或不能匹配的 SNP。识别并剔除这些不兼容的 SNP，可以确保数据的一致性和准确性。

总的来说，Harmonization 的主要目标是确保数据的一致性和准确性，通过协同 SNP 方向、剔除无法判断方向的回文 SNP 和不兼容的 SNP，从而提高孟德尔随机化分析的质量和可靠性。这些步骤都是为了减少由于数据不一致或不匹配而引入的偏见和误差，使得分析结果更加准确和可靠。

# 看看数据有哪些内容
colnames(data)
# [1] "SNP"                    "effect_allele.exposure" "other_allele.exposure" 
# [4] "effect_allele.outcome"  "other_allele.outcome"   "beta.exposure"         
# [7] "beta.outcome"           "eaf.exposure"           "eaf.outcome"           
# [10] "remove"                 "palindromic"            "ambiguous"             
# [13] "id.outcome"             "chr"                    "pos"                   
# [16] "se.outcome"             "samplesize.outcome"     "pval.outcome"          
# [19] "outcome"                "originalname.outcome"   "outcome.deprecated"    
# [22] "mr_keep.outcome"        "data_source.outcome"    "chr.exposure"          
# [25] "pos.exposure"           "se.exposure"            "pval.exposure"         
# [28] "samplesize.exposure"    "id.exposure"            "exposure"              
# [31] "mr_keep.exposure"       "pval_origin.exposure"   "data_source.exposure"  
# [34] "action"                 "SNP_index"              "mr_keep"   

# 我们可以在 Harmonization 之后检查一下数据中的 pval.outcome 列
# 看是否有 SNP 与 outcome 强相关

MR 分析正式开始

# 正式开始 MR 分析啦！
result <- mr(data)
result
#          id.exposure  id.outcome
# 1 ebi-a-GCST90029012 ukb-b-16890
# 2 ebi-a-GCST90029012 ukb-b-16890
# 3 ebi-a-GCST90029012 ukb-b-16890
# 4 ebi-a-GCST90029012 ukb-b-16890
# 5 ebi-a-GCST90029012 ukb-b-16890
#                                                       outcome
# 1 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 2 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 3 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 4 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 5 Cancer code, self-reported: breast cancer || id:ukb-b-16890
#                                                               exposure
# 1 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 2 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 3 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 4 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 5 Educational attainment (college completion) || id:ebi-a-GCST90029012
#                      method nsnp             b          se      pval
# 1                  MR Egger  198  0.0077366079 0.016666803 0.6430249
# 2           Weighted median  198  0.0003442573 0.005031674 0.9454529
# 3 Inverse variance weighted  198  0.0057960733 0.003880525 0.1352720
# 4               Simple mode  198 -0.0140984475 0.016365587 0.3900271
# 5             Weighted mode  198 -0.0167080531 0.014176057 0.2399748

我们解释一下 MR 输出结果的每列含义：

id.exposure 和 id.outcome：MR 分析中使用的暴露因子和结局的具体标识符或名称。
outcome：MR 分析中评估的结局变量，描述了研究的具体结果，例如在这里是“Cancer code, self-reported: breast cancer”。
exposure：MR 分析中评估的暴露变量，描述了研究的具体暴露，例如在这里是“Educational attainment (college completion)”。

method：使用的 MR 方法。在这里，默认使用了 MR Egger、Weighted median、Inverse variance weighted、Simple mode 和 Weighted mode 五种方法。MR 中最重要的是看 Inverse variance weighted 这个方法的 p 值。

MR 分析默认使用上述五种方法，我们也可以指定其他方法，使用mr_method_list()即可查看当前支持的统计方法，就像下面这样：

# 查看当前支持的统计方法
mr_method_list()
# obj                                                      name PubmedID                                                 Description use_by_default heterogeneity_test
# 1                  mr_wald_ratio                                                Wald ratio                                                                                TRUE              FALSE
# 2               mr_two_sample_ml                                        Maximum likelihood                                                                               FALSE               TRUE
# 3            mr_egger_regression                                                  MR Egger 26050253                                                                       TRUE               TRUE
# 4  mr_egger_regression_bootstrap                                      MR Egger (bootstrap) 26050253                                                                      FALSE              FALSE
# 5               mr_simple_median                                             Simple median                                                                               FALSE              FALSE
# 6             mr_weighted_median                                           Weighted median                                                                                TRUE              FALSE
# 7   mr_penalised_weighted_median                                 Penalised weighted median                                                                               FALSE              FALSE
# 8                         mr_ivw                                 Inverse variance weighted                                                                                TRUE               TRUE
# 9                  mr_ivw_radial                                                IVW radial                                                                               FALSE               TRUE
# 10                    mr_ivw_mre Inverse variance weighted (multiplicative random effects)                                                                               FALSE              FALSE
# 11                     mr_ivw_fe                 Inverse variance weighted (fixed effects)                                                                               FALSE              FALSE
# 12                mr_simple_mode                                               Simple mode                                                                                TRUE              FALSE
# 13              mr_weighted_mode                                             Weighted mode                                                                                TRUE              FALSE
# 14         mr_weighted_mode_nome                                      Weighted mode (NOME)                                                                               FALSE              FALSE
# 15           mr_simple_mode_nome                                        Simple mode (NOME)                                                                               FALSE              FALSE
# 16                       mr_raps                      Robust adjusted profile score (RAPS)                                                                               FALSE              FALSE
# 17                       mr_sign                                     Sign concordance test          Tests for concordance of signs between exposure and outcome          FALSE              FALSE
# 18                        mr_uwr                                     Unweighted regression                                              Doesn't use any weights          FALSE               TRUE

如果你想要指定其他方法，只需要使用method_list参数指定即可，就像下面这样：

# 指定其他方法
result_add_method <- mr(data, method_list = c("mr_ivw", "mr_egger_regression", "mr_weighted_median"))
result_add_method
#          id.exposure  id.outcome                                                     outcome
# 1 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 2 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 3 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890
#                                                               exposure                    method nsnp            b          se      pval
# 1 Educational attainment (college completion) || id:ebi-a-GCST90029012 Inverse variance weighted  198 0.0057960733 0.003880525 0.1352720
# 2 Educational attainment (college completion) || id:ebi-a-GCST90029012                  MR Egger  198 0.0077366079 0.016666803 0.6430249
# 3 Educational attainment (college completion) || id:ebi-a-GCST90029012           Weighted median  198 0.0003442573 0.005079968 0.9459706

nsnp：表示在 MR 分析中使用的 SNP 的数量。

b：这是 SNP 与暴露或结局之间关系的估计效应值，代表着两者之间的关联强度和方向。

注意注意：分类变量需要将 beta 转换 OR，连续变量就不需要啦！转换方式如下：

# beta 转换 OR 值
result_or <- generate_odds_ratios(result)
result_or
#          id.exposure  id.outcome                                                     outcome                                                             exposure
# 1 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 2 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 3 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 4 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890 Educational attainment (college completion) || id:ebi-a-GCST90029012
# 5 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890 Educational attainment (college completion) || id:ebi-a-GCST90029012
#                      method nsnp             b          se      pval        lo_ci      up_ci        or  or_lci95 or_uci95
# 1                  MR Egger  198  0.0077366079 0.016666803 0.6430249 -0.024930327 0.04040354 1.0077666 0.9753779 1.041231
# 2           Weighted median  198  0.0003442573 0.005031674 0.9454529 -0.009517824 0.01020634 1.0003443 0.9905273 1.010259
# 3 Inverse variance weighted  198  0.0057960733 0.003880525 0.1352720 -0.001809756 0.01340190 1.0058129 0.9981919 1.013492
# 4               Simple mode  198 -0.0140984475 0.016365587 0.3900271 -0.046174999 0.01797810 0.9860005 0.9548748 1.018141
# 5             Weighted mode  198 -0.0167080531 0.014176057 0.2399748 -0.044493124 0.01107702 0.9834308 0.9564822 1.011139

转换成功！

我们可以看到数据多了五列！不仅包括 OR 值，还包括其置信区间。那我们顺便简单介绍一下 β 值（效应量）转换为 OR 值（比率比）的原理、计算方法等等。

将 β 值转换为 OR 值是为了统一不同研究或数据集中的效应量，从而更容易进行比较和解释。转换的过程是基于 β 值和 OR 值之间的数学关系，通过这种方式，我们可以更清晰地理解暴露或解释变量对结果的影响，并与其他研究或数据集进行比较。

se：效应量估计的标准误，标准误是计算置信区间和 p 值的关键变量。
pval：统计显著性的 p 值。p 值用于评估 SNP 与暴露或结局之间关系的统计显著性。p 值越小，关系越可能是真实存在的。

那依据上述结果可以告诉我们什么呢？我在这里就尽力解读一下！大家可以多看相关文献，看看别人是怎么描述的，可以综合参考，然后咱们写文章的时候心里就有数啦！

咱们的 MR 分析结果显示：

MR Egger 的 p 值为 0.643，Weighted median 的 p 值为 0.945，Inverse variance weighted 的 p 值为 0.135，Simple mode 的 p 值为 0.390，Weighted mode 的 p 值为 0.240。这些 p 值都大于通常使用的显著性水平（例如 0.05），表明在这些 MR 方法中，没有找到 SNP 与教育水平或乳腺癌之间的显著关系。

**效应量（也就是 b 值）**也可以用于解释 SNP 与暴露或结局之间的关系的强度和方向。例如，正的 b 值表示 SNP 与暴露或结局之间的正相关，而负的 b 值表示负相关。

反正叭，总体来讲，这些结果表明在这些特定的 MR 方法中，并没有找到 SNP 与教育水平或乳腺癌之间的显著关系（要是真有关系就见鬼啦，啊，我是不是应该找个肯定相关的演示比较合适，啊，应该也还好，反正过程都一样对不啦，就是结果丑了点哈哈哈哈哈哈）！

敏感性分析

异质性检验

异质性（heterogeneity）检验主要用于评估所使用的 SNP 是否存在异质性。异质性指的是不同的 SNP 对暴露与结局的关系可能存在差异，这可能影响到 MR 分析的可靠性和解释。

# 异质性检验
mr_heterogeneity(data)
#          id.exposure  id.outcome                                                     outcome
# 1 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890
# 2 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890
#                                                               exposure                    method        Q Q_df      Q_pval
# 1 Educational attainment (college completion) || id:ebi-a-GCST90029012                  MR Egger 253.1165  196 0.003711509
# 2 Educational attainment (college completion) || id:ebi-a-GCST90029012 Inverse variance weighted 253.1351  197 0.004273611

# 也可以像 mr() 函数那样，指定方法
mr_heterogeneity(data, method_list = c("mr_egger_regression", "mr_ivw"))

异质性检验结果解读：

MR Egger: Q 值为 253.1165，Q_df 为 196，Q_pval 为 0.0037。这个p值（0.0037）显著小于常用的显著性水平（0.05），表明在 MR Egger 方法中，存在异质性，即不同的 SNP 对教育水平与乳腺癌之间的关系可能存在差异。
Inverse variance weighted: Q 值为 253.1351，Q_df 为 197，Q_pval 为 0.0043。这个 p 值（0.0043）同样显著小于常用的显著性水平，也表明在 Inverse variance weighted 方法（一般看这个）中，存在异质性。

其实有无异质性都不影响我们的分析，反正我们用的几乎都会使用随机效应模型（Inverse variance weighted，IVW）。当使用 IVM 进行 MR 分析时，异质性的存在并不像使用固定效应模型那样可能产生严重的影响。这是因为随机效应模型能够考虑到研究间的异质性，并为每个研究分配一个特定的权重，从而降低了异质性对分析结果的影响。我们可以像下面这样解释：

“虽然我们在分析中观察到了一些异质性，但由于我们使用的是随机效应模型（Inverse variance weighted, IVW），这些异质性对我们的分析结果影响较小。因此，我们得出的结论仍然是基于我们的主要分析方法，即IVW方法，它显示没有发现遗传变异与暴露或结果变量之间的显著关系。”

水平多效性检验

水平多效性（horizontal pleiotropy）检验在 MR 中是用于检测 SNP是否对多个结局产生效应的一种方法。多效性（pleiotropy）在这里指的是一个 SNP 对于多个表型或特征有影响的现象。水平多效性检验主要是为了检测和评估这种潜在的多效性。如果变量工具不通过暴露影响结局，这就违背了孟德尔的假说，也就说明其存在水平多效性。

# 水平多效性检验
mr_pleiotropy_test(data)
#          id.exposure  id.outcome                                                     outcome
# 1 ebi-a-GCST90029012 ukb-b-16890 Cancer code, self-reported: breast cancer || id:ukb-b-16890
#                                                               exposure egger_intercept           se      pval
# 1 Educational attainment (college completion) || id:ebi-a-GCST90029012   -1.407895e-05 0.0001175807 0.9048129

水平多效性检验结果解释：

Egger Intercept (-1.407895e-05)：这是 Egger 方法中的截距或偏见项，用于检测水平多效性。通常，如果该值显著不为 0，可能表明存在水平多效性。简单来说就是，当截距不为 0 的时候，也就是暴露在效应值为 0 的时候，对结局还有影响，那这个影响哪里来的嘞，这个效应值啊，可能是来自其他的混杂因素，也就表明存在水平多效性。在这里，Egger 截距接近 0，且 p 值为 0.905，说明没有发现显著的水平多效性。
se (0.0001175807)：效应量的标准误，用于计算 Egger 截距的置信区间和p值。
pval (0.9048129)：统计显著性的 p 值。这个 p 值很大，远大于通常的显著性水平（0.05），表明在这个水平多效性检验中没有找到显著的水平多效性。

水平多效性检验的结果表明，在这个特定的 MR 分析中，没有发现 SNP 对于教育水平或乳腺癌之间存在显著的水平多效性。这意味着 SNP 对这两个结局的影响主要是单一的，没有同时对多个结局产生显著的效应。

可视化可视化

散点图（scatter plot）

散点图用来展示遗传变异与暴露因子之间的关系。

# 开始可视化喽！

# 散点图
p1 <- mr_scatter_plot(result, data)
p1

在散点图中，每个点代表一个 SNP。横坐标通常表示 SNP 的效应量，纵坐标表示暴露因子的效应量。如果暴露因子与结局变量之间存在因果关系，我们期望看到 SNP 的效应量与暴露因子的效应量之间存在线性关系（像咱们这个结果就相关性一般般啦）。斜率通常代表因果效应的估计值。

关于散点图的详细解读，大家也可以移步：R高级绘图 | P1 | 带边缘分布散点图 | 代码注释 + 结果解读

# 要是只想展示某几种方法，在进行 MR 分析 的步骤时，指定方法即可
result_add_method <- mr(data, method_list = c("mr_ivw", "mr_egger_regression", "mr_weighted_median"))
p1_1 <- mr_scatter_plot(result_add_method, data)
p1_1

森林图（forest plot）

森林图用来展示每个 SNP 的因果效应估计值及其置信区间，主要看总效应情况。

# 森林图
result_single <- mr_singlesnp(data)
p2 <- mr_forest_plot(result_single)
p2

哎呀！非常抱歉我这 SPN 有点多，图图有点挤，大家凑合看看！不影响理解的！

在森林图中，每个条目表示一个 SNP，条目的横坐标对应值代表该 SNP 的效应量估计值，而线段代表置信区间。如果线段跨越了0点，表示效应量不显著。

对常规森林图的画法感兴趣的小伙伴，也可以移步：R语言绘图 | 高级森林图（Forest Plot）| 从数据处理到绘图细节 | 代码注释 + 结果解读

# 和上面一样，想展示指定方法的结果，在函数中指定即可
result_single_add_method <- mr_singlesnp(data, all_method = c("mr_ivw", "mr_two_sample_ml"))
p2_2 <- mr_forest_plot(result_single_add_method)
p2_2

留一图（leave-one-out plot）

我不确定中文是不是这么叫哈哈哈哈哈哈！

留一法（leave-one-out analysis）的思路是每次从数据集中剔除一个 SNP，然后使用剩余的 SNP 子集进行 MR 分析，从而得到一个没有该 SNP 的效应估计值。这个过程对于数据集中的每一个 SNP 都会执行一次。最后，我们可以比较留一法后的效应估计值与完整数据集的效应估计值，从而评估每个SNP对整体MR效应的贡献和稳定性。理想的结果是剔除后变化不大。

# 留一图
result_loo <- mr_leaveoneout(data)
p3 <- mr_leaveoneout_plot(result_loo)
p3

和上面的森林图解读类似，在留一图中，每个条目表示剔除该 SNP 的 SNP 子集进行 MR 分析结果，对应横坐标为其效应量估计值（也就是进行留一法后的效应量估计值）。判断每个 SNP 对 MR 分析结果的影响，如有离群值，则需要将其剔除后重新进行分析。

留一法交叉验证提供了一种评估每个 SNP 对整体 MR 效应稳定性和贡献的方法。通过留一图，我们可以直观地看到每个 SNP 对整体效应的影响，识别可能的“敏感” SNP，并评估整体 MR 效应的稳定性。

漏斗图（funnel plot）

漏斗图用于检查研究结果的发布偏见（publication bias）或其他偏见，其实就是对异质性的可视化。如果存在偏见，也就是有异质性，我们可能会看到效应量与其精度（通常用标准误或置信区间宽度表示）之间的关系。

# 漏斗图
result_single <- mr_singlesnp(data)
p4 <- mr_funnel_plot(result_single)
p4

在漏斗图中，每个点代表一个 SNP 的效应量，横坐标通常表示效应量，纵坐标表示精度（例如标准误或置信区间宽度）。在理想情况下，漏斗图应该是对称的，即效应量的精度（如标准误或置信区间宽度）与效应量之间无关。如果漏斗图呈现出一定程度的非对称性，这可能表明存在异质性或其他偏见，这时需要进一步考虑异质性的存在和可能的原因。此外，还要注意是否有明显偏离其他点的异常值，这些点可能需要特别关注。如果所有点都位于漏斗的底部，而没有位于顶部或中间，这可能表明存在发布偏见或其他系统性误差。

呼~