MR联合Meta分析

将多个MR研究的OR值进行Meta合并，适用于汇总多个独立MR研究的结果。

# 1. 安装并加载meta包
# install.packages("meta")
library(meta)

# 2. 设置工作目录
setwd("your/path/to/data")

# 3. 读取数据 (包含OR, 95%CI, 样本量, p值)
data <- read.csv("MRMETA.csv", header = TRUE)

# 4. 数据整理 - 将OR转换为log值，计算SE
beta <- log(data[, "OR"])
loguci <- log(data[, "or_uci95"])
loglci <- log(data[, "or_lci95"])
se <- (loguci - loglci) / (2 * 1.96)

# 5. Meta分析 (随机效应模型)
metadata <- metagen(
  TE = beta,       # 效应值 (logOR)
  seTE = se,       # 标准误
  data = data,
  sm = "OR",       # 效应量类型
  n.e = ncase,      # 病例组样本量
  n.c = ncontrol,    # 对照组样本量
  pval = p,       # p值
  random = TRUE,     # 随机效应模型
  common = FALSE,
  studlab = paste(author, year, sep = ",")
)

# 查看结果
metadata

# 6. 绘制森林图
forest(metadata)

# 保存为PDF
pdf("forest_OR.pdf", width = 12, height = 8)
forest(metadata)
dev.off()

# 7. 亚组分析 (按年份分组)
metadata1 <- metagen(
  TE = beta, seTE = se, data = data, sm = "OR",
  n.e = ncase, n.c = ncontrol, pval = p,
  random = TRUE, common = FALSE,
  studlab = paste(author, year, sep = ","),
  subgroup = year   # 按年份进行亚组分析
)
forest(metadata1)

# 8. 敏感性分析 (Leave-one-out)
metainf(metadata, pooled = "random")
forest(metainf(metadata, pooled = "random"), comb.random = TRUE)

# 9. 漏斗图 - 检测发表偏倚
funnel(metadata)
pdf("funnel.pdf")
dev.off()

# 10. 剪补法 (Trim-and-fill)
tf0 <- trimfill(metadata)
funnel(tf0, pch = ifelse(tf0$trimfill, 1, 23), 
    level = 0.95, comb.random = TRUE)

# 异质性大时用MM法
metabias(metadata, method.bias = "mm")

# 异质性小时用Harbob法
metabias(metadata, method.bias = "score")

对METAL软件输出的GWAS Meta结果进行后处理，转换为TwoSampleMR可用的标准格式。

# 1. 读取METAL结果文件
library(vroom)
res <- vroom("METAANALYSIS1.TBL", col_names = TRUE)

# 2. 根据z-score和p-value计算beta和se
SE = 1 / sqrt((2 * res$`P-value`) * (1 - res$`P-value`) * 
       (res$Weight + res$Zscore^2))
Beta = SE * res$Zscore

# 3. 将beta、se加到结果数据中
res_betase <- as.data.frame(cbind(Beta, SE))
final <- cbind(res, res_betase)

# 4. 转化为TwoSampleMR标准格式
# 标准列名: SNP, beta, se, eaf, effect_allele, other_allele, p, N
GWAS <- final[, c(1, 2, 3, 4, 6, 8, 9)]
colnames(GWAS) <- c("SNP", "effect_allele", "other_allele", 
           "N", "p", "beta", "se")

# 转换为大写
GWAS$effect_allele <- toupper(GWAS$effect_allele)
GWAS$other_allele <- toupper(GWAS$other_allele)

# 5. 导出为CSV
write.csv(GWAS, "GWAS.csv", row.names = FALSE)

对多个研究的等位基因计数数据进行Meta分析，适用于SNP层面的关联研究。

# 1. 加载meta包
library(meta)

# 2. 读取SNP数据 (包含每个研究的2x2表数据)
mydata <- read.csv("snp_meta.csv")

# 3. 二分类变量Meta分析
# 数据格式: Case.CT.TT, Case.CC, Control.CT.TT, Control.CC
metadata <- metabin(
  Case.CT.TT,       # 病例组: 携带变异基因数
  Case.CT.TT + Case.CC,  # 病例组: 总数
  Control.CT.TT,      # 对照组: 携带变异基因数
  Control.CT.TT + Control.CC, # 对照组: 总数
  data = mydata,
  sm = "OR",
  comb.fixed = FALSE,   # 固定效应
  comb.random = TRUE,   # 随机效应
  studlab = paste(Author, Year, Country, sep = ",")
)

# 4. 森林图
forest(metadata)
pdf("forest_SNP.pdf", width = 12, height = 8)
dev.off()

# 5. 亚组分析 (按国家)
metadata1 <- metabin(
  Case.CT.TT, Case.CT.TT + Case.CC,
  Control.CT.TT, Control.CT.TT + Control.CC,
  data = mydata, sm = "OR",
  comb.fixed = FALSE, comb.random = TRUE,
  studlab = paste(Author, Year, Country, sep = ","),
  byvar = Country
)
forest(metadata1)