How it all started
操作步骤
# ———————- 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_genes <- c(“PHB1”, “PHB2”)
all_gene_positions <- lapply(target_genes, get_gene_position)
# ———————- 2. 提取 eQTL 数据并筛选目标区域———————-
window <- 1e100 # 1MB = 1,000,000 bp
eqtl_study_ids <- c(“eqtl-a-ENSG00000167085”, “eqtl-a-ENSG00000215021”)
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) {
warning(“eQTL 数据获取失败，返回空数据框”)
return(data.frame())
}
# 使用实际列名 chr.exposure 和 pos.exposure 筛选
filtered_eqtl <- eqtl_data %>%
filter(
chr.exposure == !!target_chr, # 染色体列名是 chr.exposure
pos.exposure >= !!region_start, # SNP 位置 ≥ 上游边界（起始-1MB）
pos.exposure <= !!region_end # SNP 位置 ≤ 下游边界（结束+1MB）
)
return(filtered_eqtl)
}
all_exposure_eqtl <- list()
for (i in seq_along(target_genes)) {
gene_position <- all_gene_positions[[i]]
region_start <- gene_position$start – window
region_end <- gene_position$end + window
exposure_eqtl <- fetch_eqtl_in_region(
eqtl_study_id = eqtl_study_ids[i],
target_chr = gene_position$chr,
region_start = region_start,
region_end = region_end
)
all_exposure_eqtl[[i]] <- exposure_eqtl
cat(“目标区域内筛选到”, nrow(exposure_eqtl), “个 eQTL SNP\n”)
}
# 合并所有 eQTL 数据
combined_exposure_eqtl <- bind_rows(all_exposure_eqtl)
cat(“合并后筛选到”, nrow(combined_exposure_eqtl), “个 eQTL SNP\n”)
# ———————- 3. 提取结局数据（疾病/表型 GWAS） ———————
# 函数：提取与 eQTL SNP 匹配的结局 GWAS 数据（支持合并后 SNP 列表）
fetch_outcome_for_eqtl <- function(combined_eqtl_snps, outcome_study_id) {
# 提取结局数据（仅保留与合并后 eQTL SNP 匹配的 SNP）
outcome_data <- extract_outcome_data(
snps = combined_eqtl_snps$SNP, # 使用合并后的 SNP 列表
outcomes = outcome_study_id # 结局的 GWAS ID（如 “ukb-d-HEARTFAIL”）
)
# 检查匹配的 SNP 数量（改为警告而非停止）
if (is.null(outcome_data) || nrow(outcome_data) == 0) {
warning(paste(“警告：结局数据（ID=”, outcome_study_id, “）中未匹配到任何 SNP！”))
return(data.frame()) # 返回空数据框
} else if (nrow(outcome_data) < 2) {
warning(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(combined_exposure_eqtl, outcome_study_id) # 使用合并后的 eQTL 数据
if (nrow(outcome_gwas) > 0) {
cat(“结局数据中匹配到”, nrow(outcome_gwas), “个 SNP\n”)
}
# ———————- 4. 数据协调（统一 SNP 方向与格式） ———————-
harmonise_data_safe <- function(exposure_dat, outcome_dat) {
if (nrow(exposure_dat) == 0 || nrow(outcome_dat) == 0) {
warning(“数据协调失败：暴露或结局数据为空！”)
return(data.frame())
}
harmonised_data <- harmonise_data(
exposure_dat = exposure_dat, # 合并后的 eQTL 数据（暴露）
outcome_dat = outcome_dat # GWAS 数据（结局）
)
# 检查协调后 SNP 数量（改为警告而非停止）
if (nrow(harmonised_data) < 2) {
warning(paste(“警告：协调后仅剩余”, nrow(harmonised_data), “个有效 SNP，MR 分析结果可能不可靠！”))
}
return(harmonised_data)
}
harmonised_data <- harmonise_data_safe(combined_exposure_eqtl, outcome_gwas)
# ———————- 5. MR 分析（输出所有方法结果） ———————-
run_mr_analysis <- function(harmonised_data) {
if (nrow(harmonised_data) == 0) {
warning(“无法进行 MR 分析：协调后无有效 SNP！”)
return(NULL)
}
# 运行默认 MR 方法（IVW、MR-Egger、加权中位数等，取决于 SNP 数量）
mr_result <- mr(harmonised_data)
return(list(
all_methods = mr_result, # 包含所有运行的 MR 方法结果
ivw = mr_result[mr_result$method == “Inverse variance weighted”, ] # 单独提取 IVW 结果（若有）
))
}
# 执行 MR 分析
mr_result <- run_mr_analysis(harmonised_data)
# 输出所有方法结果（优化显示格式）
if (!is.null(mr_result)) {
if (nrow(mr_result$all_methods) > 0) {
cat(“\n==================== 所有 MR 方法分析结果 ====================\n”)
# 打印表头
cat(sprintf(“%-30s | %-15s | %-6s | %-10s | %-10s | %-10s\n”,
“方法”, “SNP数量”, “效应值”, “标准误”, “P值”, “方法类型”))
# 逐行打印结果
for (i in 1:nrow(mr_result$all_methods)) {
row <- mr_result$all_methods[i, ]
cat(sprintf(“%-30s | %-15d | %-10.6f | %-10.6f | %-10.6f | %s\n”,
row$method, row$nsnp, row$b, row$se, row$pval,
ifelse(row$method == “Inverse variance weighted”, “主要方法（IVW）”, “辅助方法”)))
}
cat(“=================================================================\n”)
} else {
cat(“无有效 MR 方法结果\n”)
}
# 若需要，仍可单独打印 IVW 结果（可选）
if (nrow(mr_result$ivw) > 0) {
cat(“\n—————— 逆方差加权（IVW）详细结果 ——————\n”)
print(mr_result$ivw, row.names = FALSE)
}
}
# ———————- 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 留一法分析”))。
