这是 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 存在因果影响。