要聊起今天的話題就得從流行病學(xué)講起。人們通常采用隨機對照試驗（random control trial, RCT）研究暴露因素X與疾病結(jié)局Y之間的直接關(guān)聯(lián)證據(jù)，但該方法往往受限于人類醫(yī)學(xué)倫理與諸多的試驗設(shè)計，再加上近年來廣泛出現(xiàn)的大樣本GWAS數(shù)據(jù)，表觀遺傳學(xué)以及各種“組學(xué)”數(shù)據(jù)，使得在因果關(guān)系推斷的研究中面臨諸多挑戰(zhàn)。借助孟德爾隨機化（Mendelian Randomization, MR）的設(shè)計思想，將基因或者表觀遺傳標記作為待研究暴露因素的工具變量，為解決上述問題提供了有效的途徑。

本文以筆記的形式簡要介紹孟德爾隨機化的理論基礎(chǔ)和基于GWAS數(shù)據(jù)的R包實踐。

MR法研究設(shè)計原理

首先，通過回憶高中物理知識大家都知道孟德爾遺傳的基本思想遵循"親代等位基因隨機分配給子代"，在基因型決定表型的前提假設(shè)下，基因型通過表型與疾病建立聯(lián)系，因此，基因型可以作為工具變量來推斷表型與疾病的關(guān)聯(lián)。

在該模型中，遺傳變異可以直接準確測量，并且不受外界環(huán)境等因素影響，屬于長期而穩(wěn)定的暴露因素，因此，MR設(shè)計可以最大程度的降低偏倚的作用。

兩樣本MR（Two-sample MR）

單獨介紹兩樣本MR是因為其設(shè)計策略是建立在遺傳變異-暴露因素和遺傳變異-結(jié)局變量的關(guān)聯(lián)研究人群來自相同的人群的兩個獨立樣本，如暴露因素的GWAS研究和結(jié)局變量的GWAS研究。經(jīng)過改進的兩樣本MR方法具有一個顯著的優(yōu)點，即不依賴基因型數(shù)據(jù)，只需通過現(xiàn)有GWAS結(jié)果統(tǒng)計量即可估算暴露因素與結(jié)局變量之間的因果關(guān)聯(lián)。2018年發(fā)表在《eLife》上面的一篇文章（PMID：29846171）開發(fā)了一個數(shù)據(jù)庫和R包，專門針對現(xiàn)有GWAS數(shù)據(jù)進行兩樣本MR分析。

借助巨人的肩膀——TwoSampleMR

我們以實戰(zhàn)的形式介紹軟件的使用，具體的理論基礎(chǔ)和更多的細節(jié)請讀者參考原著論文和其他相關(guān)資料。

完成整個MR分析可以分為四個步驟：

• 選擇基因型數(shù)據(jù)作為工具變量（如果是SNP，需要進行LD clumping）

• 從GWAS summary中制作工具變量的exposure和outcome文件

• 將工具變量在暴露和結(jié)局變量間的效應(yīng)值進行協(xié)同統(tǒng)一，主要是針對兩個GWAS研究所用的參考基因型不同

• 進行MR分析，敏感性分析以及繪圖等

實戰(zhàn)代碼

#安裝
install.packages("devtools")
library(devtools)
install_github("MRCIEU/TwoSampleMR")

#獲取工具變量在暴露研究中的數(shù)據(jù)
bmi_file <- system.file("data/bmi.txt", package="TwoSampleMR")
bmi_exp_dat <- read_exposure_data(bmi_file)

#Clumping
bmi_exp_dat <- clump_data(bmi_exp_dat)

#從結(jié)局變量中取出工具變量所對應(yīng)的效應(yīng)值
outcome_dat <- read_outcome_data(
    snps = bmi_exp_dat$SNP,
    filename = "gwas_summary.csv",
    sep = ",",
    snp_col = "rsid",
    beta_col = "effect",
    se_col = "SE",
    effect_allele_col = "a1",
    other_allele_col = "a2",
    eaf_col = "a1_freq",
    pval_col = "p-value",
    units_col = "Units",
    gene_col = "Gene",
    samplesize_col = "n"
)
###在R包說明文檔中，昨天首先介紹的是直接從數(shù)據(jù)庫已有的GWAS中提取信息，這種方式比較快捷，但需要用戶自行注冊goole賬號，并且數(shù)據(jù)庫中還要有自己感興趣的GWAS數(shù)據(jù)###

#Harmonise data
dat <- harmonise_data(
    exposure_dat = bmi_exp_dat, 
    outcome_dat = chd_out_dat
)

#Perform MR
res <- mr(dat)
###默認情況下會執(zhí)行多種MR算法，并分別給出結(jié)果。也支持用戶指定算法###
mr_method_list()###獲取全部MR算法
mr(dat, method_list=c("mr_egger_regression", "mr_ivw"))###指定兩種MR算法

#敏感性分析
mr_heterogeneity(dat)

#繪制散點圖
p1 <- mr_scatter_plot(res, dat)
p1[[1]]

一些需要注意的事

熟悉遺傳學(xué)研究的朋友都知道基因是具有多效性的（pleiotropy），也就是說一個SNP可能不單單與目標暴露因素有關(guān)，也存在同時與其他暴露因素有關(guān)系的可能性，在這種情況下，需要進行敏感性分析（sensitivity analysis）來確定非特異SNP的存在對結(jié)果造成的影響。另一方面，如果多個SNPs共同作為工具變量，基因多效性帶來的偏倚也會存在，可以使用MR-Egger回歸分析的方法來評價偏倚大小。

參考文獻：

王莉娜, Zuofeng Z . 孟德爾隨機化法在因果推斷中的應(yīng)用[J]. 中華流行病學(xué)雜志, 2017, 38(4):547.

Gibran Hemani, Jie Zheng, Kaitlin H Wade, Charles Laurin, Benjamin Elsworth, Stephen Burgess, Jack Bowden, Ryan Langdon, Vanessa Tan, James Yarmolinsky, Hashem A. $The MR-Base platform supports systematic causal inference across the human phenome. eLife 2018.

轉(zhuǎn)自生信草堂公眾號，已授權(quán)

文獻請在公眾號獲取~