孟德爾隨機化（Mendelian Randomization, MR）是一種利用遺傳變異作為工具變量來探索潛在因果關系的統(tǒng)計方法浇坐。這種方法可以繞開傳統(tǒng)觀察性研究中的一些常見偏誤樊诺，如混雜因素、反向因果等問題。在R語言中遏弱，MendelianRandomization包提供了一系列功能剑按，以支持單變量和多變量MR分析的實施杰刽。

主要功能和方法

單變量MR方法

• IVW (Inverse Variance Weighted)：逆方差加權方法障涯，這是MR分析中的標準方法强窖，適用于估計一個單一遺傳工具變量與表型之間的關系。
  • 逆方差加權（IVW）方法之所以被這樣命名滓鸠，主要是因為它在估計過程中對每個工具變量（例如遺傳變異）的貢獻進行加權拾碌，而這種加權是基于各個工具變量影響暴露變量的估計方差的逆吐葱。簡單來說，每個工具變量的影響被賦予一個權重校翔，該權重是其方差的倒數(shù)弟跑。這樣做的目的是為了增強估計的精確度，特別是當不同的工具變量對于暴露變量的影響估計精度不一致時防症。
  • 如果你正在研究吸煙（暴露變量）和肺癌（結果變量）之間的關系孟辑，IVW方法會使用各種與吸煙相關的遺傳變異（如煙草使用量相關的SNP），并根據(jù)這些變異的精確度加權蔫敲，從而估計吸煙對肺癌風險的總體影響饲嗽。
• Median-based：中位數(shù)方法，當工具變量可能存在無效時奈嘿，提供一個更穩(wěn)健的因果估計貌虾。
  • 在研究維生素D水平（暴露變量）與骨折風險（結果變量）的關系時，如果部分SNP與其他疾病相關聯(lián)裙犹，中位數(shù)基方法可以通過專注于中間值來提供一個較為穩(wěn)健的估計酝惧。
• MR-Egger：可以提供對工具變量多效性的敏感性分析，同時也允許存在工具變量無效的情況伯诬。
  • 在探究膽固醇水平（暴露變量）對心血管疾餐泶健（結果變量）影響的研究中，如果一些膽固醇相關的SNP也影響其他病狀（如糖尿驳了啤）哩陕，MR-Egger可以幫助識別并糾正這種多效性的影響。
• Maximum Likelihood：極大似然估計法赫舒，考慮了SNP與暴露變量關聯(lián)的不確定性悍及，并可以處理樣本重疊的問題。
  • 研究體重（暴露變量）對心臟步影（結果變量）的影響時心赶，如果使用的SNP數(shù)據(jù)來源于部分相同的研究群體，極大似然方法能夠適當考慮這種樣本重疊缺猛，提供更準確的因果關系估計缨叫。

選擇使用孟德爾隨機化（MR）中的哪種單變量方法取決于數(shù)據(jù)的特點椭符、工具變量的質量以及所面對的具體統(tǒng)計挑戰(zhàn)。下面是每種方法的使用原則和推薦的使用順序：

1. 逆方差加權（IVW）

使用原則：

• 當所有的工具變量都認為是有效的耻姥，即它們僅通過影響暴露變量來影響結果變量销钝，沒有通過任何其他途徑影響結果。
• 沒有證據(jù)顯示工具變量存在多效性（即一個遺傳變異影響多個表型）琐簇。

優(yōu)先使用條件：

• 當你擁有多個信譽良好蒸健、與暴露變量強相關的SNPs時。
• 當遺傳工具變量與暴露之間的關聯(lián)非常明確且不存在顯著的遺傳異質性時婉商。

2. MR-Egger

使用原則：

• 用于檢測和校正多效性偏差似忧，即某些SNP可能通過不相關的途徑影響結果。
• 當存在工具變量無效性的疑慮時丈秩，可以用來估計這種偏差的大小和方向橡娄。

優(yōu)先使用條件：

• 當研究中的SNP可能不完全滿足排除限制（僅通過暴露影響結果）時。
• 當需要評估和校正多效性偏差時癣籽，作為分析的一部分進行敏感性測試挽唉。

3. 中位數(shù)基方法（Median-based）

使用原則：

• 當存在一部分工具變量可能無效（即它們可能同時影響其他表型）時，提供一種更為穩(wěn)健的因果估計方法筷狼。
• 當數(shù)據(jù)中存在異常值或者工具變量效應大小的分布極為不均時瓶籽。

優(yōu)先使用條件：

• 在對IVW結果的穩(wěn)健性有疑慮時，作為一種補充分析埂材。
• 當部分SNP的有效性受到質疑塑顺，或者存在顯著的遺傳異質性時。

4. 極大似然（Maximum Likelihood）

使用原則：

• 考慮到SNP與暴露之間關聯(lián)的不確定性俏险，同時可以處理樣本重疊的問題严拒。
• 當需要精細調整樣本重疊對因果估計的潛在影響時。

優(yōu)先使用條件：

• 當數(shù)據(jù)來自可能存在樣本重疊的不同研究時竖独，或者當單個研究中用于暴露和結果的樣本部分相同時裤唠。
• 當分析需要精確控制遺傳工具的不確定性和復雜性時。

總體策略

• 開始：通常以IVW方法作為首選莹痢，因為它直接种蘸、簡單且在工具變量全部有效時非常強大。
• 敏感性分析：使用MR-Egger檢查多效性和工具變量無效性的影響竞膳；同時航瞭，運用中位數(shù)方法來評估結果的穩(wěn)健性。
• 復雜情境處理：在面對樣本重疊或高度復雜的數(shù)據(jù)結構時坦辟，采用極大似然方法刊侯。

多變量MR方法

• MV-IVW：多變量逆方差加權，同時考慮多個遺傳工具變量锉走。
• MV-Egger：多變量Egger回歸滨彻，可以在多變量情境下評估多效性和工具變量無效的影響藕届。

里面還有多種方法，注意帶有mr_m為多變量的疮绷，不帶有的為單變量

image.png

應用示例

使用MendelianRandomization包進行MR分析的R代碼示例：

# 加載R包
library(MendelianRandomization)

# 輸入數(shù)據(jù)準備
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)

# 使用逆方差加權方法
IVWObject <- mr_ivw(MRInputObject)

# 使用中位數(shù)方法進行穩(wěn)健分析
MedianObject <- mr_median(MRInputObject)

# 使用MR-Egger方法分析潛在的多效性
EggerObject <- mr_egger(MRInputObject)

# 使用極大似然法考慮SNP-exposure關聯(lián)的不確定性
MaxLikObject <- mr_maxlik(MRInputObject)

# 多變量MR分析
MRMVInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig), bxse = cbind(ldlcse, hdlcse, trigse), by = chdlodds, byse = chdloddsse)
MRMVObject <- mr_mvivw(MRMVInputObject)

數(shù)據(jù)可視化

MendelianRandomization包還提供了數(shù)據(jù)可視化功能，可以用于生成散點圖嚣潜、森林圖等冬骚，幫助用戶評估MR分析的結果：

# 繪制基于IVW方法的MR分析結果圖
mr_plot(MRInputObject, line="ivw")

# 繪制多方法比較的圖形
mr_plot(mr_allmethods(MRInputObject, method="all"))

示例數(shù)據(jù)讀取

用了 MendelianRandomization 包中的 extract.pheno.csv 函數(shù)，該函數(shù)用于提取和整理特定的遺傳和表型數(shù)據(jù)懂算，以便進行孟德爾隨機化分析只冻。

# 定義指向包內(nèi)預設數(shù)據(jù)文件的路徑，這些文件包含PhenoScanner數(shù)據(jù)庫的SNP信息
path.noproxy <- system.file("extdata", "vitD_snps_PhenoScanner.csv",
                            package = "MendelianRandomization")
path.proxies <- system.file("extdata", "vitD_snps_PhenoScanner_proxies.csv",
                            package = "MendelianRandomization")


# 提取不使用代理SNPs的數(shù)據(jù)计技，用于暴露變量“l(fā)og(eGFR creatinine)”和結果變量“Tanner stage”
# 這些數(shù)據(jù)針對的是歐洲血統(tǒng)的研究
extract.pheno.csv(
  exposure = "log(eGFR creatinine)",  # 指定暴露變量
  pmidE = 26831199,  # 暴露數(shù)據(jù)的PubMed ID
  ancestryE = "European",  # 暴露數(shù)據(jù)的人種背景
  outcome = "Tanner stage",  # 指定結果變量
  pmidO = 24770850,  # 結果數(shù)據(jù)的PubMed ID
  ancestryO = "European",  # 結果數(shù)據(jù)的人種背景
  file = path.noproxy  # 指定使用不包含代理的數(shù)據(jù)文件
)

# 提取使用代理SNPs的數(shù)據(jù)喜德，同樣針對暴露變量“l(fā)og(eGFR creatinine)”和結果變量“Tanner stage”
extract.pheno.csv(
  exposure = "log(eGFR creatinine)", 
  pmidE = 26831199, 
  ancestryE = "European",
  outcome = "Tanner stage", 
  pmidO = 24770850, 
  ancestryO = "European",
  rsq.proxy = 0.6,  # 設置代理SNPs的R2閾值為0.6，表示中等關聯(lián)強度
  file = path.proxies  # 指定使用包含代理的數(shù)據(jù)文件
)

# 提取使用代理SNPs的數(shù)據(jù)垮媒，暴露變量“l(fā)og(eGFR creatinine)”和不同的結果變量“Asthma”
extract.pheno.csv(
  exposure = "log(eGFR creatinine)", 
  pmidE = 26831199, 
  ancestryE = "European",
  outcome = "Asthma",  # 指定不同的結果變量為“哮喘”
  pmidO = 20860503,  # 哮喘數(shù)據(jù)的PubMed ID
  ancestryO = "European",
  rsq.proxy = 0.6,  # 同樣設置代理SNPs的R2閾值為0.6
  file = path.proxies  # 使用包含代理的數(shù)據(jù)文件
)

從PhenoScanner提供的示例數(shù)據(jù)集中提取了用于孟德爾隨機化的必要遺傳信息舍悯。其中，rsq.proxy參數(shù)指定了接受的代理SNP與主SNP之間的最低相關性閾值睡雇，用于在主SNP數(shù)據(jù)不足時擴充可用的遺傳工具變量萌衬。這使得能夠在數(shù)據(jù)可用性受限的情況下，繼續(xù)進行有效的MR分析它抱。

單變量方法

IVW方法

逆方差加權（IVW）方法基于一個關鍵假設：所有的遺傳變異都是有效的工具變量秕豫，比如它們只通過影響暴露變量（如LDL cholesterol）來影響結果變量（如心血管疾病風險）。在這種假設下观蓄，可以將每個SNP的估計效應加權平均混移，以得到因果效應的整體估計。權重通常是每個估計的逆方差侮穿，以提高精確度較高的估計的貢獻度歌径。

# 加載MendelianRandomization包
library(MendelianRandomization)

# 創(chuàng)建輸入對象，包括每個SNP對暴露和結果的效應大小及其標準誤
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)
print(MRInputObject)

# 使用IVW方法進行估算亲茅，不使用穩(wěn)健回歸沮脖，不對SNP中的outliers進行懲罰
IVWObject1 <- mr_ivw(MRInputObject, model= "default", robust = FALSE, penalized = FALSE,
                     correl = FALSE, weights = "simple", psi = 0, distribution = "normal", alpha = 0.05)
print(IVWObject1)

# 使用穩(wěn)健回歸，并對SNP中的outliers進行懲罰芯急，以提高估計的穩(wěn)健性
IVWObject2 <- mr_ivw(MRInputObject, model= "default", robust = TRUE, penalized = TRUE,
                     correl = FALSE, weights = "simple", psi = 0, distribution = "normal", alpha = 0.05)
print(IVWObject2)

median-based方法

在孟德爾隨機化分析中勺届，中位數(shù)基方法（Median-based）是用來提供對個別SNP異常值或無效工具變量的穩(wěn)健因果估計。這種方法可以有效地減少異常數(shù)據(jù)點對總體估計結果的影響娶耍。下面的代碼片段展示了如何使用不同的加權策略來執(zhí)行中位數(shù)基方法免姿，并注釋解釋了每個函數(shù)調用的目的和參數(shù)的意義。

# 加載所需的MendelianRandomization包
library(MendelianRandomization)

# 創(chuàng)建輸入數(shù)據(jù)對象榕酒，包括SNP對LDL膽固醇和冠心病風險的影響及其標準誤
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)

# 使用加權的中位數(shù)方法進行因果估計胚膊，加權以反映每個估計的不確定性
WeightedMedianObject1 <- mr_median(MRInputObject, weighting = "weighted", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265)
print(WeightedMedianObject1)  # 打印結果故俐，顯示LDL升高與增加的CHD風險之間的顯著關聯(lián)

# 使用“penalized”加權法，這種方法在考慮權重時對潛在的異常值或強影響點進行懲罰
WeightedMedianObject2 <- mr_median(MRInputObject, weighting = "penalized", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265)
print(WeightedMedianObject2)  # 打印結果

# 使用簡單中位數(shù)方法紊婉，不進行任何加權
WeightedMedianObject3 <- mr_median(MRInputObject, weighting = "simple", distribution = "normal", alpha = 0.05, iterations = 10000, seed = 314159265)
print(WeightedMedianObject3)  # 打印結果

# 增加迭代次數(shù)到100,000药版，以提高估計的精度和穩(wěn)健性
WeightedMedianObject4 <- mr_median(MRInputObject, weighting = "weighted", distribution = "normal", alpha = 0.05, iterations = 100000, seed = 314159265)
print(WeightedMedianObject4)  # 打印結果

參數(shù)解釋

• weighting: 控制如何加權每個SNP的估計∮骼纾可選值包括"weighted"（考慮標準誤的逆）槽片，"penalized"（對潛在的異常值施加懲罰），和"simple"（不加權）肢础。
• distribution: 指定使用的分布類型还栓，這里使用正態(tài)分布（"normal"）。
• alpha: 顯著性水平传轰，常設置為0.05剩盒。
• iterations: 迭代次數(shù)，用于計算的迭代過程中提高估計的精度慨蛙。
• seed: 隨機數(shù)種子辽聊，確保結果的可重復性。

示例展示了如何在不同場景下應用中位數(shù)方法來評估LDL對冠心病風險的影響期贫，同時考慮到了數(shù)據(jù)中可能存在的異常值和統(tǒng)計的穩(wěn)健性身隐。通過調整加權策略和迭代次數(shù)，可以根據(jù)具體研究需求優(yōu)化因果估計的準確性和穩(wěn)健性唯灵。

MR-Egger方法

MR-Egger 方法在孟德爾隨機化分析中被用來檢測和調整工具變量的潛在多效性或偏倚贾铝。該方法類似于常規(guī)的回歸分析，但增加了對工具變量多效性偏倚的估計和調整埠帕。這可以幫助研究者評估工具變量是否僅通過暴露變量影響結果垢揩，還是也通過其他路徑影響結果，這種偏倚被稱為多效性偏倚或遺傳混雜敛瓷。

# 加載所需的MendelianRandomization包叁巨，這是執(zhí)行MR分析的必要工具包
library(MendelianRandomization)

# 創(chuàng)建輸入數(shù)據(jù)對象，包括SNP對LDL膽固醇和冠心病風險的影響及其標準誤
MRInputObject <- mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse)

# 執(zhí)行MR-Egger回歸分析呐籽，未應用穩(wěn)健回歸或懲罰方法
EggerObject1 <- mr_egger(MRInputObject, robust = FALSE, penalized = FALSE, correl = FALSE, distribution = "normal", alpha = 0.05)
print(EggerObject1)  # 打印結果锋勺，提供不考慮異常值調整的基本MR-Egger估計

# 執(zhí)行MR-Egger回歸分析，應用穩(wěn)健回歸和懲罰方法
EggerObject2 <- mr_egger(MRInputObject, robust = TRUE, penalized = TRUE, correl = FALSE, distribution = "normal", alpha = 0.05)
print(EggerObject2)  # 打印結果狡蝶，提供考慮異常值調整后的穩(wěn)健MR-Egger估計

# 分析和討論結果
# 使用穩(wěn)健回歸后庶橱，MR估計值的誤差變小了，但顯著性并未改變贪惹。
# 這表明LDL膽固醇水平的升高與增加的冠心病風險之間的因果關聯(lián)是穩(wěn)健的苏章，即使考慮潛在的多效性偏倚。

參數(shù)解釋

• robust: 是否使用穩(wěn)健回歸，幫助減少異常值和影響點對回歸分析結果的影響枫绅。
• penalized: 是否應用懲罰項泉孩，用于減少模型的過擬合和增強模型預測的穩(wěn)定性。
• correl: 這個參數(shù)用于指示是否考慮SNP之間的相關性并淋，通常在MR-Egger分析中設為FALSE寓搬。
• distribution: 指定模型估計的分布，通常為正態(tài)分布（"normal"）县耽。
• alpha: 顯著性檢驗的閾值句喷，常設置為0.05。

Maximum likelihood方法（極大似然估計法)

極大似然方法（Maximum Likelihood, ML）在孟德爾隨機化（Mendelian Randomization, MR）分析中提供了一個強有力的工具酬诀，尤其是在處理樣本重疊和SNP與暴露之間關聯(lián)不確定性的情況下脏嚷。這個方法可以通過調整psi參數(shù)來適應不同程度的樣本重疊骆撇，從而提供更為精確的因果推斷瞒御。

極大似然方法的優(yōu)勢

1. SNP-Exposure關聯(lián)的不確定性：極大似然估計法允許研究者考慮到每個SNP與暴露之間關聯(lián)強度的不確定性，這在簡單的逆方差加權方法（IVW）中通常是被忽略的神郊。
2. 樣本重疊情況：當使用來自相同樣本的遺傳和表型數(shù)據(jù)時肴裙，可能會出現(xiàn)樣本重疊，這會影響到估計的準確性涌乳。極大似然方法通過psi參數(shù)調整蜻懦，可以有效管理這一點。

參數(shù)psi的作用

• psi = 0：表示樣本之間完全不重疊夕晓，即傳統(tǒng)的獨立雙樣本MR研究場景宛乃。
• psi > 0：表示樣本之間有重疊，psi的大小代表暴露和結果變量之間的相關性蒸辆，即觀察性研究中暴露和結局的相關系數(shù)征炼。

這里展示了如何使用不同的psi值進行極大似然估計，以探索樣本重疊對因果估計的影響：

library(MendelianRandomization)

# 極大似然估計躬贡，樣本不重疊
MaxLikObject1 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0, distribution = "normal", alpha = 0.05)
print(MaxLikObject1)

# 極大似然估計谆奥，樣本輕微重疊，psi = 0.3
MaxLikObject2 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0.3, distribution= "normal", alpha = 0.05)
print(MaxLikObject2)

# 極大似然估計拂玻，樣本中等重疊酸些，psi = 0.6
MaxLikObject3 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0.6, distribution= "normal", alpha = 0.05)
print(MaxLikObject3)

# 極大似然估計，樣本高度重疊檐蚜，psi = 0.9
MaxLikObject4 <- mr_maxlik(MRInputObject, model = "default", correl = FALSE, psi = 0.9, distribution= "normal", alpha = 0.05)
print(MaxLikObject4)

通過這種方式魄懂，研究者可以評估樣本重疊在MR分析中的實際影響，確保因果估計的可靠性和精確性闯第。

多樣本

MRMVInputObject <- mr_mvinput(bx = cbind(ldlc, hdlc, trig),
                              bxse = cbind(ldlcse, hdlcse, trigse),
                              by = chdlodds, byse = chdloddsse) #MVMR的input格式會和單變量的有所不同
MRMVInputObject

MRMVObject1 <- mr_mvivw(MRMVInputObject,
                        model = "default",
                        correl = FALSE,
                        distribution = "normal",
                        alpha = 0.05)
MRMVObject2 <-  mr_mvegger(MRMVInputObject, 
                           orientate = 1,
                           correl = FALSE,
                           distribution = "normal", 
                           alpha = 0.05)

繪圖

結果取決于傳遞給 mr_plot 的對象類型逢渔。當對象是 MRInput對象，該函數(shù)使用 plot 命令（如果 interactive 設置為 FALSE）或 plotly語法（如果 interactive 設置為 TRUE）來繪制關聯(lián)估計值乡括。當肃廓。智厌。。的時候object 是一個 MRMVInput 對象盲赊，功能類似铣鹏，除了我們繪制估計的關聯(lián)結果在 y 軸上，關聯(lián)的擬合值與來自x軸上的逆方差加權方法哀蘑。如果 interactive 設置為 FALSE，則為靜態(tài)圖被生產(chǎn)绘迁。通過將標簽設置為 TRUE，基因變體的名稱出現(xiàn)在點上方缀台。這會產(chǎn)生一個視覺上不太吸引人的圖表，但更容易識別個人遺傳變異膛腐。如果 interactive 設置為 TRUE睛约，則繪圖是交互式的哲身，用戶可以懸停在各個點上查看相關遺傳變異的名稱及其關聯(lián)估計辩涝。當對象是 MRAll 對象時勘天，該函數(shù)生成一個 ggplot 來比較因果估計通過不同的方法提出怔揩。

mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
        line="egger", orientate = TRUE)

image.png

mr_plot(mr_input(bx = ldlc, bxse = ldlcse, by = chdlodds, byse = chdloddsse),
        line="ivw", interactive=FALSE) # produces a static graph

image.png

mr_plot(mr_allmethods(mr_input(bx = ldlc, bxse = ldlcse,
                               by = chdlodds, byse = chdloddsse), method="all", iterations = 50))

image.png

https://zhuanlan.zhihu.com/p/439009634

參考文獻和資源

• CRAN - Package MendelianRandomization

• MendelianRandomization.pdf

• Mendelian Randomization: Methods for Using Genetic Variants in Causal Estimation