How it all started
操作步骤
一、准备工作
1.先到NCBI中查询基因，获取ENSG号
2.打开IEU OpenGWAS project网站，输入ENSG号，获取GWAS ID
3.获取并配置token
3.1 登录opengwas，注册一下
登录注册，用github账号，https://api.opengwas.io/profile/
3.2 创建token
把token复制一下。
3.3 在R语言中检测是否有token
如果之前没有设置过，肯定是没有的。那就需要在文档文件夹中，新建一个.Renviron文件
3.4 把token放到新建的.Renviron文件中
OPENGWAS_JWT=”这里粘贴你的token”
3.5 重启R语言（必须）
重启R语言，然后键入下面代码，测试token是否设置成功：
## 测试token是否有效
library(ieugwasr)
ieugwasr::get_opengwas_jwt() #返回秘钥则提示成功
二、代码运行
# ———————- 0. 环境初始化与依赖加载 ———————-
# 清空环境变量
rm(list = ls())
options(scipen = 999) # 禁止科学计数法，提升结果可读性
# 检查并安装必要工具包（pacman 用于批量管理包）
if (!require(“pacman”)) install.packages(“pacman”, update = FALSE, ask = FALSE)
library(pacman)
# 设置 Bioconductor 镜像（加速下载）
options(BioC_mirror = “https://mirrors.ustc.edu.cn/bioc/”)
if (!require(“pacman”)) install.packages(“pacman”)
pacman::p_load(
biomaRt, # 获取基因在基因组中的位置信息（染色体、起始/结束坐标）
dplyr, # 数据清洗与筛选（如筛选上下游 SNP）
tidyr # 数据整理（可选）
)
# 批量加载/安装依赖包（核心包：TwoSampleMR 用于 MR 分析；ieuwasr 用于获取 GWAS 数据）
library(devtools)
library(TwoSampleMR)
library(ieugwasr)
library(biomaRt)
library(dplyr)
library(tidyr)
# ———————- 1. 获取目标基因的基因组位置信息 ———————-
get_gene_position <- function(gene_symbol) {
ensembl <- useEnsembl(biomart = “genes”, dataset = “hsapiens_gene_ensembl”)
gene_info <- getBM(
attributes = c(“hgnc_symbol”, “chromosome_name”, “start_position”, “end_position”),
filters = “hgnc_symbol”,
values = gene_symbol,
mart = ensembl
)
if (nrow(gene_info) == 0) stop(“基因位置获取失败”)
list(chr = as.character(gene_info$chromosome_name[1]),
start = as.integer(gene_info$start_position[1]),
end = as.integer(gene_info$end_position[1]))
}
target_gene <- “PHB1”
gene_position <- get_gene_position(target_gene)
# ———————- 2. 提取 eQTL 数据并筛选目标区域———————-
window <- 1e6
region_start <- gene_position$start – window
region_end <- gene_position$end + window
# 正确的函数定义（参数名 target_chr）
fetch_eqtl_in_region <- function(eqtl_study_id, target_chr, region_start, region_end) {
eqtl_data <- extract_instruments(eqtl_study_id)
if (is.null(eqtl_data) || nrow(eqtl_data) == 0) stop(“eQTL 数据获取失败”)
# 使用实际列名 chr.exposure 和 pos.exposure 筛选
filtered_eqtl <- eqtl_data %>%
filter(
chr.exposure == !!target_chr, # 染色体列名是 chr.exposure
pos.exposure >= !!region_start, # 位置列名是 pos.exposure
pos.exposure <= !!region_end
)
if (nrow(filtered_eqtl) < 2) stop(“目标区域 SNP 数量不足”)
return(filtered_eqtl)
}
# 正确的函数调用（参数名 target_chr）
eqtl_study_id <- “eqtl-a-ENSG00000167085”
exposure_eqtl <- fetch_eqtl_in_region(
eqtl_study_id = eqtl_study_id,
target_chr = gene_position$chr,
region_start = region_start,
region_end = region_end
)
cat(“目标区域内筛选到”, nrow(exposure_eqtl), “个 eQTL SNP\n”)
# ———————- 3. 提取结局数据（疾病/表型 GWAS） ———————-
# 函数：提取与 eQTL SNP 匹配的结局 GWAS 数据
fetch_outcome_for_eqtl <- function(eqtl_snps, outcome_study_id) {
# 提取结局数据（仅保留与 eQTL SNP 匹配的 SNP）
outcome_data <- extract_outcome_data(
snps = eqtl_snps$SNP, # 使用 eQTL 数据中的 SNP 列表
outcomes = outcome_study_id # 结局的 GWAS ID（如 “ukb-d-HEARTFAIL” 表示心力衰竭）
)
# 检查匹配的 SNP 数量
if (is.null(outcome_data) || nrow(outcome_data) < 2) {
stop(paste(“错误：结局数据（ID=”, outcome_study_id, “）中仅匹配到”,
nrow(outcome_data), “个 SNP，无法进行 MR 分析！”))
}
return(outcome_data)
}
# 示例：获取心力衰竭的 GWAS 数据（替换为你的结局 ID）
outcome_study_id <- “ukb-d-HEARTFAIL”
outcome_gwas <- fetch_outcome_for_eqtl(exposure_eqtl, outcome_study_id)
cat(“结局数据中匹配到”, nrow(outcome_gwas), “个 SNP\n”)
# ———————- 4. 数据协调（统一 SNP 方向与格式） ———————-
harmonised_data <- harmonise_data(
exposure_dat = exposure_eqtl, # eQTL 数据（暴露）
outcome_dat = outcome_gwas # GWAS 数据（结局）
)
# 关键检查：协调后 SNP 数量（需 ≥2）
if (nrow(harmonised_data) < 2) {
stop(“错误：协调后有效 SNP 数量为”, nrow(harmonised_data), “，无法进行多 SNP MR 分析！”)
}
# ———————- 5. MR 分析（默认多种方法） ———————-
# 运行默认 MR 方法（IVW、MR-Egger、加权中位数等）
mr_result <- mr(harmonised_data)
# 查看核心结果（以逆方差加权法 IVW 为例）
ivw_result <- mr_result[mr_result$method == “Inverse variance weighted”, ]
print(“逆方差加权（IVW）分析结果：”)
print(ivw_result)
# ———————- 6. 可视化与结果保存 ———————-
# 森林图：对比不同 MR 方法的效应值
mr_forest_plot(mr_result, title = paste(target_gene, “eQTL-MR 分析结果”))
# 散点图：展示 SNP 的暴露-结局效应关系
mr_scatter_plot(mr_result, dat = harmonised_data, title = paste(target_gene, “eQTL 与结局关联散点图”))
# 保存结果到 CSV
write.csv(mr_result, paste0(target_gene, “_eqtl_mr_results.csv”), row.names = FALSE)
# ———————- 7. 敏感性分析（评估结果稳健性） ———————-
# 水平多效性检验（MR-Egger 截距）
pleiotropy_test <- mr_pleiotropy_test(harmonised_data)
print(“水平多效性检验结果（P < 0.05 提示存在多效性）：”)
print(pleiotropy_test)
# 异质性检验（Cochran’s Q）
heterogeneity_test <- mr_heterogeneity_test(harmonised_data)
print(“异质性检验结果（P < 0.05 提示 SNP 间效应不一致）：”)
print(heterogeneity_test)
# 留一法分析（排除单个 SNP 后的结果变化）
leaveoneout_result <- mr_leaveoneout(harmonised_data)
mr_leaveoneout_plot(leaveoneout_result, title = paste(target_gene, “eQTL-MR 留一法分析”))
结局和暴露只有1个SNP
# 分析前准备
# 使用“ls()”函数获取当前环境中的所有变量、函数和对象的名称，然后传递给”list”参数
# 将“list”参数传递给“rm()”函数，表示删除这些变量、函数和对象
# 因此，“rm(list=ls())”表示删除当前环境中的所有变量、函数和对象，相当于清空当前环境
rm(list = ls())
# 判断是否已经安装”pacman”，如果没有则安装它
if(!require(“pacman”)) install.packages(“pacman”, update=F, ask=F)
# 设置Bioconductor镜像地址为中国科技大学的镜像
options(BioC_mirror=”https://mirrors.ustc.edu.cn/bioc/”)
# 加载”pacman”包，方便加载其他得到R包
library(“pacman”)
# 下载并打开devtools
install.packages(“devtools”)
library(devtools)
# 下载并打开TwoSampleMR
devtools::install_github(“MRCIEU/TwoSampleMR”,force = TRUE)
library(TwoSampleMR) # 引用包
# 下载并打开ieuwasr，用于直接访问IEU GWAS数据库
install.packages(“ieuwasr”)
library(ieugwasr)
# 2024.5.1后都要token才可以继续访问gwas服务器，否则提示Error in if (nrow(d) == 0) return(NULL) : 参数长度为零
# Sys.getenv(“R_USER”) # 获取.Renviron文件
# file.edit(‘~/.Renviron’) # 如果不存在则创建配置文件.Renvion，将秘钥复制到其中
# ieugwasr::get_opengwas_jwt() # 返回秘钥则提示成功
user() # 检查返回数据
# 暴露数据，获取靶基因的eQTL数据
# 示例：从IEU OpenGWAS获取PHB1的eQTL数据（ID为eqtl-a-ENSG00000167085）
phb1_eqtl <- extract_instruments(“eqtl-a-ENSG00000167085”)
# 看一下数据
head(phb1_eqtl)
# 从指定 GWAS 数据集中提取与给定 SNP 关联的结局数据
hf_outcome <- extract_outcome_data(phb1_eqtl$SNP,outcomes = “ukb-d-HEARTFAIL”)
# rm(hf_outcome) # 如果获取失败或获取结果为NULL，运行此行，删除hf_outcome，更换GWAS ID以便再次进行获取
# 协调暴露与结局数据
dat_hf <- harmonise_data(phb1_eqtl, hf_outcome)
# 运行Wald比率法MR分析
if(nrow(dat_hf) == 1) {
# 计算Wald比率
beta_exposure <- dat_hf$beta.exposure
beta_outcome <- dat_hf$beta.outcome
se_exposure <- dat_hf$se.exposure
se_outcome <- dat_hf$se.outcome
mr_estimate <- beta_outcome / beta_exposure
mr_se <- sqrt((se_outcome^2 / beta_exposure^2) + (beta_outcome^2 * se_exposure^2 / beta_exposure^4))
mr_p <- 2 * pnorm(-abs(mr_estimate / mr_se))
res_hf <- data.frame(
method = “Wald ratio”,
b = mr_estimate,
se = mr_se,
p = mr_p
)
} else {
res_hf <- mr(dat_hf, method_list = “mr_wald_ratio”)
}
# 查看分析结果
print(res_hf)
#结果解读
该结果表格包含三列，各列含义及结果解读如下：
method：分析方法为 “wald ratio”（Wald 比率法）。
b：效应估计值，为 -0.0001447013，接近 0，表明暴露（基因表达Prohibitin1）对结局（心力衰竭，Heart Failure）的影响幅度极小。
se：标准误，为 0.001005473，反映效应估计值的离散程度。
p：假设检验的 P 值，为 0.8855687，大于常用的显著性阈值（如 0.05），说明在当前分析中，未发现暴露与心力衰竭结局之间存在统计学意义的关联。
综上，基于 Wald 比率法的分析结果，未显示靶基因（通过 eQTL 数据体现）与心力衰竭结局之间有显著的因果关联，效应估计值本身也极小，进一步支持无实质关联的结论。
# 由于只有一个SNP，故无法进行敏感性分析。
