PheWAS 实操：ieugwasr 包

安装依赖包

# 安装 ieugwasr 包
install.packages("ieugwasr")
# 或从 GitHub 安装最新版本
devtools::install_github("MRCIEU/ieugwasr")

# 安装其他依赖
install.packages(c("TwoSampleMR", "data.table", "dplyr", "openxlsx"))

先了解数据库中有哪些可用表型，便于后续筛选。

library(TwoSampleMR)

# 获取 IEU 数据库中所有可用 GWAS 汇总数据
ao <- available_outcomes()

# 查看前几行
head(ao)

# 保存完整表型列表
write.csv(ao, "available_outcomes.csv", row.names = FALSE)

# 筛选 UK Biobank 相关表型（可选）
# UK Biobank 数据通常以 "ukb" 开头
ukb_ao <- subset(ao, grepl("ukb", id) & author != "Pan-UKB team")

# 保存 UK Biobank 表型列表
write.csv(ukb_ao, "ukb_outcomes.csv", row.names = FALSE)

cat("总表型数:", nrow(ao), "\n")
cat("UK Biobank 表型数:", nrow(ukb_ao), "\n")

使用 ieugwasr 的 phewas 函数进行全表型关联分析。

library(ieugwasr)

# 设置工作目录
setwd("XXXXX/your-data-directory")

# 选择感兴趣的 SNP（示例：rs507666 - ABO 血型相关 SNP）
snp_id <- "rs507666"

# 执行 PheWAS 分析
# phewas 函数会自动扫描所有可用表型
# pval = 1 表示返回所有结果不过滤
phewas_result <- ieugwasr::phewas(snp_id, pval = 1)

# 查看结果维度
cat("分析表型数:", nrow(phewas_result), "\n")

# 保存完整结果
write.csv(phewas_result, "phewas_result_full.csv", row.names = FALSE)

phewas 结果中的 beta 是对数尺度，需要转换为 OR 值便于解读。

# 读取 phewas 结果（如果已保存）
phewas_result <- read.csv("phewas_result_full.csv")

# 读取可用表型信息（获取表型名称）
ao <- available_outcomes()

# 合并表型名称
phewas_result <- merge(phewas_result, ao[, c("id", "trait")], by = "id", all.x = TRUE)

# 计算 95% 置信区间
phewas_result$lo_ci <- phewas_result$beta - 1.96 * phewas_result$se
phewas_result$up_ci <- phewas_result$beta + 1.96 * phewas_result$se

# 转换为 OR 值
phewas_result$or <- exp(phewas_result$beta)
phewas_result$or_lci95 <- exp(phewas_result$lo_ci)
phewas_result$or_uci95 <- exp(phewas_result$up_ci)

# 保存带 OR 值的结果
write.csv(phewas_result, "phewas_result_or.csv", row.names = FALSE)

# 查看显著结果（p < 0.05）
significant <- phewas_result[phewas_result$p < 0.05, ]
cat("显著结果数 (p<0.05):", nrow(significant), "\n")

# FDR 校正
phewas_result$P_Adjust <- p.adjust(phewas_result$p, method = "fdr")
significant_fdr <- phewas_result[phewas_result$P_Adjust < 0.05, ]
cat("FDR 校正后显著 (P_Adjust<0.05):", nrow(significant_fdr), "\n")

如果只关注 UK Biobank 表型，可以筛选子集进行分析。

# 读取之前的 phewas 结果
phewas_result <- read.csv("phewas_result_or.csv")

# 读取 UK Biobank 表型列表
ukb_ao <- read.csv("ukb_outcomes.csv")

# 筛选 UK Biobank 相关表型
phewas_ukb <- subset(phewas_result, id %in% ukb_ao$id)

# 保存 UK Biobank 专属结果
write.csv(phewas_ukb, "phewas_ukb.csv", row.names = FALSE)

# 查看显著结果
sig_ukb <- phewas_ukb[phewas_ukb$p < 0.05, ]
cat("UK Biobank 显著结果数:", nrow(sig_ukb), "\n")

# 按 p 值排序查看 top 结果
sig_ukb <- sig_ukb[order(sig_ukb$p), ]
print(head(sig_ukb[, c("trait", "or", "or_lci95", "or_uci95", "p")], 10))

如果同一表型有多个数据源，可以用 Z 检验比较效应差异。

# Z 检验函数：比较两个来源的效应是否显著不同
Z_test <- function(data, trait_name) {
  # 筛选同一表型的两个来源
  subset_data <- data[data$trait == trait_name, ]
  
  if (nrow(subset_data) < 2) return(NULL)
  
  beta1 <- subset_data$beta[1]
  se1 <- subset_data$se[1]
  beta2 <- subset_data$beta[2]
  se2 <- subset_data$se[2]
  
  # 计算效应差异
  beta_diff <- beta1 - beta2
  se_diff <- sqrt(se1^2 + se2^2)
  z_score <- beta_diff / se_diff
  p_value <- 2 * (1 - pnorm(abs(z_score)))
  
  return(data.frame(
    trait = trait_name,
    z_score = z_score,
    p_value = p_value,
    significant = p_value < 0.05
  ))
}

# 示例：比较某个表型的两个来源
# 需要先筛选出有两个来源的表型
# trait_counts <- table(phewas_result$trait)
# multi_source_traits <- names(trait_counts[trait_counts > 1])

# 对每个多来源表型执行 Z 检验
# results <- lapply(multi_source_traits, function(t) Z_test(phewas_result, t))
# results <- do.call(rbind, results)