這是 GRAPPLE(Pervasive PLEiotropy 下的全基因組 mR 分析)的 R 包,GRAPPLE 是孟德爾隨機化的綜合框架。
有關GRAPPLE框架的詳細信息,請參閱我們的手稿。
GRAPPLE 套件可以從這個 Github 儲存庫安裝在 R 中:
library(devtools)
install_github("jingshuw/grapple")
如果需要 GRAPPLE 進行資料預處理,可能還需要安裝 PLINK 進行 LD 聚類,並為 PLINK 準備合適的 LD 聚類參考資料集。歐洲人口聚類參考資料集的資源(也在GRAPPLE 論文中使用)是1000 個基因組歐洲參考面板,可在此下載:http://fileserve.mrcieu.ac.uk/ld/data_maf0.01_rs_ref. tgz ,來自 MRCIEU 引用。
使用 GRAPPLE 進行的分析從預處理原始 GWAS 摘要統計數據開始,其中包括選擇遺傳儀器、協調數據集以共享相同的效應和參考等位基因、計算 SNP 之間共享的估計邊緣關聯中的噪聲相關矩陣,並選擇用於模式標記檢測的候選SNP清單。
為此,GRAPPLE 需要原始 GWAS 摘要統計檔案清單。每個 GWAS 文件都是一個「.csv」或「.txt」文件,包含至少 6 個欄位的資料框,其列名稱為: SNP
、 effect_allele
、 other_allele
、 beta
、 se
、 pval
。 SNP
欄位包含每個 SNP 的 rsID。 effect_allele
和other_allele
列都需要大寫字母。 beta
列包含連續性狀的估計效果量和二元性狀的對數優勢比, se
列是對應beta
的標準差。
GRAPPLE 論文中所使用的上述格式的所有 GWAS 摘要統計檔案都可以從資料集中下載。有關每個文件的原始公開來源,請參閱手稿補充資訊的第 3 節。
GRAPPLE的Deta預處理需要一個疾病文件(結果文件),每個危險因子兩份文件,其中一個用於SNP選擇(選擇文件),另一個用於估計SNP對該危險因子的邊際效應(曝光文件) 。我們允許結果和暴露文件中的 GWAS 資料出現重疊個體,選擇文件和其他文件之間應該有盡可能少的重疊個體,以避免 SNP 選擇偏差。有時很難找到這樣的 GWAS 資料集來進行 SNP 選擇。一種選擇是使用來自其他祖先的 GWAS 數據。
例如,可以從資料集中下載 BMI 和 T2D 的摘要統計文件,也可以下載參考面板並將「.tgz」文件解壓縮到目前資料夾中。然後,我們可以從這些文件中選擇 SNP 作為遺傳工具。
library(GRAPPLE)
sel.file <- "BMI-ukb.csv"
exp.file <- "BMI-giant17eu.csv"
out.file <- "T2D-diagram12-M.csv"
plink_refdat <- "./data_maf0.01_rs_ref/data_maf0.01_rs_ref"
提取獨立 SNP 可能需要幾分鐘時間,具體取決於 GWAS 摘要統計文件的大小。
data.list <- getInput(sel.file, exp.file, out.file, plink_refdat, max.p.thres = 0.01, plink_exe = 'plink')
此函數提取所有選擇 p 值不超過 0.01 的獨立 SNP。它傳回三個元素。一個是“data”,它是所選 SNP 的總計統計的資料框,另一個是“marker.data”,它是所有候選標記 SNP 的資料框。 「data」和「marker.data」之間的差異在於,我們在 LD 聚集中使用比選擇候選標記 SNP (r2 = 0.05) 更嚴格的 r2 (0.001) 來選擇遺傳工具,以確保獨立性。第三個元素是 GWAS 隊列的估計相關矩陣。
我們使用 PLINK 來執行 LD 聚類。 plink 指令的預設名稱/路徑是「plink」(透過參數plink_exe
傳遞)。對於 Linux / Windows 系統的用戶,此命令不起作用,需要指定 exe 檔案的路徑,例如“./plink”,具體取決於他們安裝 plink 的位置。如果R無法執行plink命令,則會出現錯誤,指出未找到clumped檔案。
如果風險因素的數量為 1,GRAPPLE 可以透過在給定 p 值選擇閾值的情況下查找穩健輪廓可能性中的模式數量來檢測多效性途徑的數量。
我們可以嘗試上面的例子。
## Here we take the p-value threshold be 1e-4, but one can try a series of p-values and see if the result is consistent
diagnosis <- findModes(data.list$data, p.thres = 1e-4)
diagnosis$modes
diagnosis$p
可能性只有一種模式,因此我們沒有發現多種多效性途徑的任何證據。它為我們使用MR-RAP來估計風險因子的因果效應提供了一定的保證。
result <- grappleRobustEst(data.list$data, p.thres = 1e-4)
該函數傳回我們的估計結果。
這是我們可以檢測多種模式的另一個範例
library(GRAPPLE)
sel.file <- "CRP-Prins17.csv"
exp.file <- "CRP-Dehghan11.csv"
out.file <- "CAD-Nelson17.csv"
data.list <- getInput(sel.file, exp.file, out.file, plink_refdat, max.p.thres = 0.01)
diagnosis <- findModes(data.list$data, p.thres = 1e-4, marker.data = data.list$marker.data)
diagnosis$p
我們發現輪廓可能性的三種模式。檢查標記SNP、標記基因和映射的GWAS性狀,我們可以發現LDL-C可能是混雜性狀。然後我們可以執行 GRAPPLE 來調整 LDL-C 的影響
library(GRAPPLE)
sel.file <-c(sel.file, "LDL-gera18.csv")
exp.file <- c(exp.file, "LDL-glgc13.csv")
out.file <- "CAD-Nelson17.csv"
data.list <- getInput(sel.file, exp.file, out.file, plink_refdat, max.p.thres = 0.01)
result <- grappleRobustEst(data.list$data, p.thres = 1e-4)
在調整 LDL-C 後,我們找不到任何證據顯示 CRP 對 CAD 有因果影響。