meta分析整合多次差异分析结果

小洁忘了怎么分身

发布于 2026-01-27 14:16:46

1270

今天是生信星球陪你的第1065天

公众号里的文章大多数需要编程基础，如果因为代码看不懂，而跟不上正文的节奏，可以来找我学习，相当于给自己一个新手保护期。我的课程都是循环开课，点进去咨询微信，随时可以报名↓ 生信分析直播课程(每月初开一期) 生信新手保护学习小组长期报名中（每月一期）单细胞陪伴学习小组（每月一期）

背景知识

在数据挖掘中，多个数据差异分析之后的联合，有多种方法：

简单取交集 (Venn Diagram)

这是最直观的方法。分别对两个数据集做差异分析，然后画 Venn 图取重叠部分。

优点：逻辑简单，读者一看就懂。缺点：存在严重的**“漏斗效应”**。很多基因可能在一个数据集中显著（P<0.05），在另一个数据集中处于边缘显著（比如 P=0.06），直接取交集会无情地丢弃这些其实很重要的基因。

效应量 Meta 分析 (Effect Size Meta-analysis) 这是目前的“黄金标准”。它不是简单的数值平均，而是根据每个研究的标准误 (SE) 进行加权。数据越准、样本量越大的研究，话语权越重。

优点：能综合考虑变化方向、幅度和数据质量，结果最稳健可信，且支持跨平台（芯片+测序）整合。缺点：需要完整的 logFC 和 SE 数据，门槛稍高。

秩聚合算法 (RRA) 这是一种基于“排名”的整合方法。它不看具体的表达量数值，而是看基因在各个列表中是否都排在前列。

优点：容错率极高，能整合量纲完全不同的数据（如转录组+蛋白质组+甲基化）。缺点：丢失了具体的“倍数变化”信息，无法量化差异幅度。

直接合并去批次 (ComBat/SVA) 将两个表达矩阵强行拼合，并使用算法去除批次效应。

适用场景：通常仅限于同平台的数据（如都是 Affymetrix 芯片）。风险：如果强行合并芯片和测序数据，极易引入非线性偏差；且容易“过度校正”，误把生物学差异当成批次效应抹除。

Fisher's Method (P值合并) 这是 Meta 分析的简易版，仅针对 P 值进行合并（将多个 P 值融合成一个统计量）。

适用场景：当你只有 Gene Symbol 和 P-value，没有 logFC 时的救命稻草。硬伤：它分不清方向。如果某基因在一个数据集中上调，在另一个中下调，Fisher 法可能会误判为“显著”，导致假阳性。

共识加权基因共表达网络 (Consensus WGCNA) 这是一种进阶策略。它不局限于寻找单个差异基因，而是挖掘在两个数据集中“共存的基因模块”。

价值：如果 TCGA 和 GEO 都能构建出同一个共表达模块，说明该通路在生物学上非常保守且核心，适合讲述更高维度的调控网络故事。

今天介绍这里面的方法2(效应量 Meta 分析)的代码。

1.载入数据

rm(list = ls())
library(dplyr)
library(metafor)
library(ggplot2)

load("degs.Rdata")
head(deg1)

##             logFC  AveExpr         t      P.Value    adj.P.Val
## Mug2    -3.836029 6.883903 -27.90973 2.676041e-12 2.764409e-08
## Orm2     3.054059 6.040447  24.29455 1.380244e-11 6.244226e-08
## H2-K1   -3.346620 8.882681 -23.16197 2.422945e-11 8.221052e-08
## Gm10260  5.463509 4.118297  26.93697 4.073695e-12 2.764409e-08
## Rps18   -2.921302 7.197548 -21.07808 7.337186e-11 1.492161e-07
## Cd36     4.472500 6.481416  20.99240 7.696087e-11 1.492161e-07
##                B change  Symbol
## Mug2    18.54363   DOWN    Mug2
## Orm2    17.01484     UP    Orm2
## H2-K1   16.56112   DOWN   H2-K1
## Gm10260 16.45028     UP Gm10260
## Rps18   15.48881   DOWN   Rps18
## Cd36    15.39094     UP    Cd36

head(deg2)

##                   logFC   AveExpr         t      P.Value
## 1600002H07Rik -2.520145  6.292221 -18.32099 5.872667e-10
## Cyp3a11        2.321189 11.282366  14.21200 1.008935e-08
## Slc37a4       -1.857103  6.303268 -12.39953 4.520266e-08
## Cdip1         -1.658108  6.411937 -12.35854 4.686372e-08
## Pfkfb1        -1.851370  5.545345 -12.36736 4.650084e-08
## Mgll          -1.495570  6.969573 -12.01666 6.358669e-08
##                  adj.P.Val         B change        Symbol
## 1600002H07Rik 8.027349e-06 12.798728   DOWN 1600002H07Rik
## Cyp3a11       6.895568e-05 10.406196     UP       Cyp3a11
## Slc37a4       9.446013e-05  9.069660   DOWN       Slc37a4
## Cdip1         9.446013e-05  9.042104   DOWN         Cdip1
## Pfkfb1        9.446013e-05  8.997119   DOWN        Pfkfb1
## Mgll          1.045003e-04  8.765187   DOWN          Mgll

我的两个数据实验设计相同，物种相同，都用了limma包做差异分析，deg1和deg2是两个数据各自差异分析的结果表格。如果是多次差异分析，也同样支持，准备多次差异分析的结果，然后后续代码rma的method参数从FE改为REML即可。

2.数据清洗

1.两个差异分析结果都需要有Symbol、logFC和SE(标准误)列。

差异分析常用的三大R包中，deseq2直接提供了SE值，即lfcSE 列。

limma的差异分析结果中没有提供SE值，可以使用SE = |logFC| / t计算。edgeR或者t-test也没有提供SE值，但提供了p值和logFC，也可以计算SE。

edgeR的SE计算示例代码(不运行):

DEG2 <- DEG2 %>%
  mutate(Symbol = rownames(.)) %>%
  mutate(
    # 1. 计算 Z 分数 
    Z_score = qnorm(1 - (PValue / 2)),
    # 2. 反推 SE
    SE = abs(logFC) / Z_score
  )

2.列名修改，避免多个数据相同的信息列名混乱。

3.剔除基因名NA、重复以及无穷值，防止因异常数据而报错。

4.两个差异分析结果数据取交集。

deg1 <- deg1 %>%
  transmute(
    Symbol,
    logFC_dat1 = logFC,
    SE_dat1    = abs(logFC / t),
    P_dat1     = P.Value
  ) %>%
  filter(!is.na(Symbol)) %>%
  distinct(Symbol, .keep_all = TRUE) %>%
  filter(is.finite(logFC_dat1), is.finite(SE_dat1), SE_dat1 > 0)

deg2 <- deg2 %>%
  transmute(
    Symbol,
    logFC_dat2 = logFC,
    SE_dat2    = abs(logFC / t),
    P_dat2     = P.Value
  ) %>%
  filter(!is.na(Symbol)) %>%
  distinct(Symbol, .keep_all = TRUE) %>%
  filter(is.finite(logFC_dat2), is.finite(SE_dat2), SE_dat2 > 0)

dat <- inner_join(deg1, deg2, by = "Symbol")
head(dat)

##    Symbol logFC_dat1   SE_dat1       P_dat1  logFC_dat2   SE_dat2
## 1    Mug2  -3.836029 0.1374442 2.676041e-12  0.49740075 0.2591405
## 2    Orm2   3.054059 0.1257096 1.380244e-11  1.17220273 0.4348370
## 3   H2-K1  -3.346620 0.1444877 2.422945e-11 -0.09312900 0.1379991
## 4 Gm10260   5.463509 0.2028257 4.073695e-12 -0.02046381 0.4529364
## 5   Rps18  -2.921302 0.1385943 7.337186e-11  0.31219117 0.1859874
## 6    Cd36   4.472500 0.2130533 7.696087e-11  0.40889796 0.3730486
##       P_dat2
## 1 0.07971635
## 2 0.01988533
## 3 0.51292278
## 4 0.96472774
## 5 0.11979833
## 6 0.29514775

3.执行meta分析

在 Meta 分析中，通常有固定效应（FE）和随机效应（REML）两种模型。两研究使用固定效应 FE，三个及以上研究选择随机效应（REML）。因为只有两个数据时，异质性估计极不稳定，固定效应（FE）模型反而是更稳健、更敏感的选择。

针对上述数据的每一行，定义计算函数。

run_meta_per_gene <- function(row_data) {
  yi  <- c(as.numeric(row_data["logFC_dat1"]), as.numeric(row_data["logFC_dat2"]))
  sei <- c(as.numeric(row_data["SE_dat1"]),    as.numeric(row_data["SE_dat2"]))

# 排除异常值
if (any(!is.finite(yi)) || any(!is.finite(sei)) || any(sei <= 0)) {
    return(c(NA_real_, NA_real_))
  }

  vi <- sei^2

  res <- tryCatch(
    rma(yi = yi, vi = vi, method = "FE"),   # 固定效应
    error = function(e) NULL
  )

if (is.null(res)) return(c(NA_real_, NA_real_))

  c(as.numeric(res$beta), as.numeric(res$pval))
}

# 批量运行
res <- apply(dat, 1, run_meta_per_gene)
res <- as.data.frame(t(res))
colnames(res) <- c("Pooled_ES", "Meta_Pval")
head(res)

##   Pooled_ES     Meta_Pval
## 1 -2.884635 9.304021e-125
## 2  2.908911 3.380353e-128
## 3 -1.645183  4.656462e-61
## 4  4.547510 2.909540e-133
## 5 -1.766831  6.490362e-57
## 6  3.473059  1.266056e-78

Pooled_ES即Pooled Effect Size,合并效应值，是加权后的共识 logFC，它是两个数据集中 logFC 的综合结论。

Pooled Effect Size不是简单的 (dat1 + dat2) / 2。它是基于数据质量的“加权平均值”。

Meta 分析算法会根据每个研究的标准误 (SE) 自动分配话语权： SE 越小，说明该研究数据越稳、精度越高，它在合并时的权重（Weight）就越大。 SE 越大，说明数据波动大、噪音多，它的权重就越小。

举个例子： dat1 (SE极小) 显示：某基因下调 4倍。 dat2 (SE较大) 显示：该基因下调 2倍。最终结果：Pooled ES 可能是 3.5倍。 (结果会强烈向更靠谱的 dat1 倾斜，而不是简单的中间值 3)。

4.meta分析结果整理

res <- bind_cols(dat, res) %>%
  mutate(
    FDR = p.adjust(Meta_Pval, method = "BH")
  ) %>%
  filter(is.finite(Pooled_ES), is.finite(Meta_Pval), is.finite(FDR)) %>%
  # 防止FDR极小被数值舍入到0，导致 -log10(FDR) = Inf
  mutate(FDR_plot = pmax(FDR, 1e-300))

head(res)

##    Symbol logFC_dat1   SE_dat1       P_dat1  logFC_dat2   SE_dat2
## 1    Mug2  -3.836029 0.1374442 2.676041e-12  0.49740075 0.2591405
## 2    Orm2   3.054059 0.1257096 1.380244e-11  1.17220273 0.4348370
## 3   H2-K1  -3.346620 0.1444877 2.422945e-11 -0.09312900 0.1379991
## 4 Gm10260   5.463509 0.2028257 4.073695e-12 -0.02046381 0.4529364
## 5   Rps18  -2.921302 0.1385943 7.337186e-11  0.31219117 0.1859874
## 6    Cd36   4.472500 0.2130533 7.696087e-11  0.40889796 0.3730486
##       P_dat2 Pooled_ES     Meta_Pval           FDR      FDR_plot
## 1 0.07971635 -2.884635 9.304021e-125 2.998918e-121 2.998918e-121
## 2 0.01988533  2.908911 3.380353e-128 1.452763e-124 1.452763e-124
## 3 0.51292278 -1.645183  4.656462e-61  6.670641e-58  6.670641e-58
## 4 0.96472774  4.547510 2.909540e-133 1.875635e-129 1.875635e-129
## 5 0.11979833 -1.766831  6.490362e-57  7.607294e-54  7.607294e-54
## 6 0.29514775  3.473059  1.266056e-78  2.720544e-75  2.720544e-75

5.绘制火山图

cutoff_fdr <- 0.05
cutoff_es  <- 0.585  # 约等于1.5倍

res <- res %>%
  mutate(
    group = case_when(
      FDR < cutoff_fdr & Pooled_ES >  cutoff_es ~ "Up",
      FDR < cutoff_fdr & Pooled_ES < -cutoff_es ~ "Down",
      TRUE ~ "Stable"
    )
  )

table(res$group)

## 
##   Down Stable     Up 
##   1043  10737   1113

my_colors <- c("Up" = "#D53E4F", "Down" = "#3288BD", "Stable" = "grey80")

p <- ggplot(res, aes(x = Pooled_ES, y = -log10(FDR_plot))) +
  geom_point(
    data = subset(res, group == "Stable"),
    color = "grey80", size = 1, alpha = 0.5
  ) +
  geom_point(
    data = subset(res, group != "Stable"),
    aes(color = group, size = -log10(FDR_plot)),
    alpha = 0.8
  ) +
  scale_color_manual(values = my_colors) +
  scale_size_continuous(range = c(1, 4), guide = "none") +
  geom_vline(xintercept = c(-cutoff_es, cutoff_es), linetype = "dashed") +
  geom_hline(yintercept = -log10(cutoff_fdr), linetype = "dashed") +
  labs(
    x = "Pooled Effect Size",
    y = "-log10(FDR)",
    title = "Meta-analysis Volcano Plot"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    legend.position = "none",
    axis.title = element_text(size = 14),
    axis.text  = element_text(size = 12),
    plot.title = element_text(hjust = 0.5)
  )

p

在 Meta 分析的火山图中，横坐标不是单一研究的 logFC，而是综合计算出的 Pooled ES：

Pooled ES > 0 (右侧)：综合所有证据后，判定该基因在处理组整体上调。例如 Pooled ES = 2，意味着综合来看表达量约升高了 2^2=4 倍。 Pooled ES < 0 (左侧)：综合判定该基因整体下调。 Pooled ES ≈ 0 (中间)：基因无明显变化，或者两个数据集结论严重冲突（一个极高、一个极低，导致互相抵消）。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2026-01-26，如有侵权请联系 cloudcommunity@tencent.com 删除

算法

芯片