MR Base Overview

MR Base是一个用于孟德尔随机化的在线数据库以及分析平台，由Bristol大学的MRC Integrative Epidemiology Unit开发。它包含了一个在线网页应用（http://mrbase.org）， MRInstruments包以及TwoSampleMR包。

安装R包

访问NHGRI-EBI GWAS目录和必要的数据格式化功能需要安装MRInstruments和TwoSampleMR 包。这些将由 MR_Practicals R包自动安装。
要直接从 GitHub 安装MRPracticals，请首先确保您安装了 R 3.5 或更高版本、R Studio 和 R 工具 3.5 版。

接下来，打开 R Studio 并安装devtools、knitr和rmarkdown包:

install.packages(c("devtools", "knitr", "rmarkdown"))

接下来，我们需要从 Github 安装TwoSampleMR和MRInstrumentsR 包。

library(devtools)
install_github(c("MRCIEU/TwoSampleMR","MRCIEU/MRInstruments"))

然后MRPracticals可以使用以下命令安装软件包：

library(devtools)
install_github("WSpiller/MRPracticals",build_opts = c("--no-resave-data", "--no-manual"),build_vignettes = TRUE)

要更新包，只需install_github("WSpiller/MRPracticals", build_opts = c("--no-resave-data", "--no-manual"))再次运行命令。

MR Workflow

Two-sample summary MR 的流程如下：

选择感兴趣的暴露变量的instruments，根据连锁不平衡（Linkage disequilibrium，LD）清洗数据。
根据感兴趣的结局从MR Base GWAS数据库提取instruments。
协调（Harmonise）暴露和结果数据集，每个样本使用相同的参考等位基因。
执行MR分析，敏感性分析，汇集报告。

Step 1: 获取暴露摘要估计

GWAS 目录

NHGRI-EBI GWAS目录包含从GWAS获得的重要关联目录。数据的MR Base版本包含具有执行MR所需数据的协调研究，确保用于报告特定研究效应大小的单位都相同。
要使用GWAS目录，我们可以首先使用以下命令访问包含所有可用暴露信息的参考数据集：

data(gwas_catalog)

加载此参考数据集后，我们可以使用 grepl() 函数搜索感兴趣的表型。例如，如果我们想使用关键词“Blood”查找表型，则可以使用以下命令：

exposure_gwas <- subset(gwas_catalog, grepl("Blood", Phenotype_simple))

我们还可以按作者搜索具体研究：

exposure_gwas <- subset(gwas_catalog, grepl("Neale", Author))

这些命令允许我们缩小搜索条件，直到获得所需的数据。在后面这个例子中，我们将使用身体质量指数（BMI）作为暴露，使用Locke等人（2015）确定的genetic instruments。

exposure_gwas <- subset(gwas_catalog, grepl("Locke", Author) &
                     Phenotype_simple == "Body mass index")

head(exposure_gwas[,c(7:12,18:21)])

为暴露变量手动选择特定的GWAS estimates

每种instrument都有多个条目，还提供特定于性别的GWAS数据。在这个例子中，我们有兴趣使用对所有具有欧洲血统（European ancestry，EA）的参与者进行的GWAS，可以通过运行以下命令来实现的：

exposure_gwas <- subset(exposure_gwas, grepl("EA", Phenotype_info))
exposure_gwas <-exposure_gwas[!grepl("women", exposure_gwas$Phenotype_info),]
exposure_gwas <-exposure_gwas[!grepl("men", exposure_gwas$Phenotype_info),]
head(exposure_gwas[,c(7:12,18:21)])

##          Phenotype Phenotype_info PubmedID   Author Year        SNP
## 2  Body mass index           (EA) 25673413 Locke AE 2015  rs1000940
## 4  Body mass index           (EA) 25673413 Locke AE 2015 rs10132280
## 8  Body mass index           (EA) 25673413 Locke AE 2015  rs1016287
## 11 Body mass index           (EA) 25673413 Locke AE 2015 rs10167079
## 13 Body mass index           (EA) 25673413 Locke AE 2015 rs10182181
## 16 Body mass index           (EA) 25673413 Locke AE 2015 rs10269783
##    effect_allele other_allele  beta          se
## 2              G            A 0.019 0.003316327
## 4              C            A 0.023 0.003571429
## 8              T         <NA> 0.023 0.003571429
## 11             G            A 0.024 0.005102041
## 13             G            A 0.031 0.003061224
## 16             A            G 0.014 0.003112245

在这里，我们简单地删除了所有性别特定的estimates。为了确保我们的instruments足够strong，我们必须将 p 值小于 5×10⁻⁸ 的关联纳入：

exposure_gwas<-exposure_gwas[exposure_gwas$pval<5*10^-8,]

然后，我们使用format_data（）函数创建一个数据集。如果任何instruments缺少进行MR分析的必要信息，该函数还将显示警告消息。

exposure_data<-format_data(exposure_gwas)

最后，根据LD清洗SNP数据集很重要，因为使用相关联的SNP会导致重复计数，从而导致有偏差的因果效应估计。可以使用clump_data函数来完成：

exposure_data<-clump_data(exposure_data, clump_r2 = 0.001)

识别相关联的SNP的阈值已设置为默认值（0.001），但如果合理，可以更改此阈值。执行上述操作后，总共有62个独立的instruments适合后续MR分析。

自动获取instruments

如果感兴趣的暴露变量的ID号已知，则可以使用extract_instruments（）函数提取instruments。例如，Locke 等人(2015) 的研究ID号是2，因此可以通过运行以下内容来获得该instrument：

quick_extract<-extract_instruments(2)

但请注意，由于此函数使用默认值，因此与手动选择instruments相比，instruments数量可能会有所不同。在这种情况下，我们有 79 种instruments，而手动获得的有 62 种。其中一个原因是自动选择包括从混合人群而不是欧洲血统的个体中获得估计值。这两种方法都有其优点，尽管手动选择可以更好地控制用于instruments选择的标准。

Step 2 获取结局摘要估计

一旦确定了暴露变量的instruments，就需要从结局变量中提取这些SNP。MR Base包含来自大量研究的完整GWAS摘要统计数据。要获取有关可用 GWAS 的详细信息，我们可以运行以下命令：

ao<-available_outcomes()
head(ao)[,c(3,4,6,8,9,20)]

##      category consortium   id ncase ncontrol    nsnp sample_size
## 1 Risk factor   ADIPOGen    1    NA       NA 2675209       39883
## 2     Disease     IIBDGC   10 14763    15977   13898       30740
## 3 Risk factor      GIANT  100    NA       NA 2725796       60586
## 4 Risk factor      SSGAC 1000    NA       NA 6524475      161460
## 5 Risk factor      SSGAC 1001    NA       NA 8146841      293723
## 6 Risk factor       <NA> 1002    NA       NA 2474010       32161
##                 trait
## 1         Adiponectin
## 2     Crohn's disease
## 3   Hip circumference
## 4 Depressive symptoms
## 5  Years of schooling
## 6              Leptin

请注意，当访问MR Base中的数据时，将弹出一个窗口，要求您登录Google帐户进行授权。请登录，之后可以继续分析。
head 函数仅显示可用于查找感兴趣的GWAS结局变量的有用列的摘要。与手动选择暴露数据集的情况一样，我们可以使用类似的命令来搜索感兴趣的结果。例如，要搜索收缩压（SBP），请执行以下操作：

outcome_gwas <- subset(ao, grepl("Systolic", trait))
head(outcome_gwas)[,c(3,4,6,8,9,11,16,20)]

##          category consortium          id ncase ncontrol     nsnp sample_size
## 971   Risk factor       ICBP         827    NA       NA  2673126       74064
## 1436         <NA>  Neale Lab   UKB-a:360     0        0 10894596      317754
## 12815  Continuous    MRC-IEU UKB-b:20175     0        0  9851867      436419
## 18087  Continuous    MRC-IEU  UKB-b:6503     0        0  9851867       39749
##                                            trait
## 971                      Systolic blood pressure
## 1436  Systolic blood pressure  automated reading
## 12815 Systolic blood pressure, automated reading
## 18087    Systolic blood pressure, manual reading

在本例中，我们将选择来自UK Biobank（Neale Lab）的数据，其中包含最多数量的遗传变异。本研究的ID号为UKB-a:360。
现在已知所需结局变量研究的 ID，我们需要提取与步骤 1 中获得的instruments相对应的 SNP 集。

outcome_data <- extract_outcome_data(
    snps = exposure_data$SNP, outcomes = "UKB-a:360")

上面的内容中，我们指定了格式化暴露数据框中的SNP rsid编号列，以及包含结局数据的研究的ID号。

代理LD （LD proxys）的说明

默认情况下，如果结局变量GWAS中不存在特定请求的SNP，则将搜索LD中具有请求的SNP（目标）的SNP（代理）。LD 代理使用 1000 个基因组欧洲样本数据定义。返回代理SNP对结局变量的影响，以及代理SNP，代理SNP的效果等位基因以及目标SNP的相应等位基因（相位）。
数据清洗：现在获得了暴露数据和结局数据，但对数据进行清洗很重要。这意味着SNP对暴露的影响以及该SNP对结局的影响必须分别对应于相同的等位基因。为了清洗暴露和结局数据，我们可以运行以下内容：

H_data <- harmonise_data(
    exposure_dat = exposure_data, 
    outcome_dat = outcome_data
)

重复条目

数据清洗后，可能会发现他们的数据集包含重复的暴露-结局摘要集。例如，当 GWAS 针对同一性状发布了来自单独 GWAS 分析的多个结果时，可能会出现这种情况。我们建议清洗其数据集，以便仅保留具有最高预期功效的暴露-结局组合。这可以通过使用power.prune函数选择具有最大样本量的暴露结局摘要集来完成：

H_data<-power_prune(H_data)

Step 3 进行MR分析

获取效应估计值

一旦暴露和结局数据协调一致，我们就有了每个SNP可用于暴露和结局特征的效果估计和标准误。我们可以使用此信息来执行MR。为此，只需运行：

mr_results<-mr(H_data)
mr_results

##   id.exposure id.outcome
## 1      4tEuYO  UKB-a:360
## 2      4tEuYO  UKB-a:360
## 3      4tEuYO  UKB-a:360
## 4      4tEuYO  UKB-a:360
## 5      4tEuYO  UKB-a:360
##                                                      outcome        exposure
## 1 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
## 2 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
## 3 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
## 4 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
## 5 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
##                      method nsnp          b         se         pval
## 1                  MR Egger   62 0.03778843 0.06393780 5.567279e-01
## 2           Weighted median   62 0.10794749 0.02261103 1.805066e-06
## 3 Inverse variance weighted   62 0.07046435 0.02636638 7.528678e-03
## 4               Simple mode   62 0.10996751 0.06425331 9.207540e-02
## 5             Weighted mode   62 0.10996751 0.03070508 6.782338e-04

在这种情况下，由于没有明确具体方法，因此已应用并提出了一系列常见的双样本汇总MR方法。可以使用以下内容缩小方法集的范围：

mr(H_data, method_list=c("mr_egger_regression", "mr_ivw"))

##   id.exposure id.outcome
## 1      4tEuYO  UKB-a:360
## 2      4tEuYO  UKB-a:360
##                                                      outcome        exposure
## 1 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
## 2 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
##                      method nsnp          b         se        pval
## 1                  MR Egger   62 0.03778843 0.06393780 0.556727910
## 2 Inverse variance weighted   62 0.07046435 0.02636638 0.007528678

可以通过运行以下命令获得可用方法的完整列表：

mr_method_list()

head(mr_method_list())[,1:2]

##                             obj                 name
## 1                 mr_wald_ratio           Wald ratio
## 2              mr_two_sample_ml   Maximum likelihood
## 3           mr_egger_regression             MR Egger
## 4 mr_egger_regression_bootstrap MR Egger (bootstrap)
## 5              mr_simple_median        Simple median
## 6            mr_weighted_median      Weighted median

生成具有 95% 置信区间的OR值

发表文章时，经常需要列出OR的95%CI。分析二元结局变量时，通过运行以下命令log-OR转变为OR：

generate_odds_ratios(mr_results)

请注意，本例中的分析使用的是连续结局变量

敏感性分析

虽然使用mr()函数可以获得一组使用一系列方法的效应估计，但评估异质性作为可能的多效性偏差的度量是有用的。为了测试平均多效性效应，我们可以评估MR-Egger截距不为零的程度：

mr_pleiotropy_test(H_data)

##   id.exposure id.outcome
## 1      4tEuYO  UKB-a:360
##                                                      outcome        exposure
## 1 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
##   egger_intercept          se      pval
## 1      0.00106525 0.001896711 0.5764598

此外，我们可以通过运行以下命令获得 IVW 和 MR-Egger 异质性的 Q 统计量：

mr_heterogeneity(H_data, method_list=c("mr_egger_regression", "mr_ivw"))

##   id.exposure id.outcome
## 1      4tEuYO  UKB-a:360
## 2      4tEuYO  UKB-a:360
##                                                      outcome        exposure
## 1 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
## 2 Systolic blood pressure  automated reading || id:UKB-a:360 Body mass index
##                      method        Q Q_df       Q_pval
## 1                  MR Egger 291.6671   60 3.681426e-32
## 2 Inverse variance weighted 293.2004   61 4.432881e-32

在这里，我们看到使用95%阈值的MR-Egger截距并不显着，但在全局水平上，单个SNP估计之间似乎存在异质性。这由 Q 统计量和异质性的相应 p 值表示。对这些结果的一种可能的解释是，一些SNP是多效性的，但平均多效性效应接近于零（因此是平衡的）。在后续分析中，我们将使用 RadialMR包根据它们对总体异质性的贡献来研究潜在的异常值。

生成MR结果的图示

我们可以使用散点图描述 SNP 效应对暴露的影响与 SNP 效应对结局的关系：

plot1 <- mr_scatter_plot(mr_results, H_data)
plot1

plot1

我们还可以使用从单个 SNP 获得的估计值生成森林图：

res_single <- mr_singlesnp(H_data)
plot2 <- mr_forest_plot(res_single)
plot2

plot2

以及显示leave-one-out分析结果的图：

res_loo <- mr_leaveoneout(H_data)
plot3 <- mr_leaveoneout_plot(res_loo)
plot3

plot3

最后，我们可以生成一个漏斗图来评估异质性：

plot4 <- mr_funnel_plot(res_single)
plot4

plot4

孟德尔随机化之MRBase和TwoSample MR Practical