RRA整合多次差异分析结果

小洁忘了怎么分身

发布于 2026-02-04 11:21:28

1150

之前演示了如何使用meta分析整合多次差异分析结果，本文演示如何用RRA算法整合。

1. 数据下载和整理

示例数据是从GEO数据库下载3个胃癌与癌旁样本的基因表达芯片数据，实际使用可以是芯片与转录组或者单细胞混搭，只要都是按照p值排好顺序的基因即可，上下调基因要分开整合。

# 清理环境
rm(list = ls())
# 加载用于GEO数据下载和差异分析的包
library(tinyarray)

# 将目标GEO数据集编号存入向量
datasets = c("GSE19826", "GSE54129", "GSE79973")

# 依次下载3个GEO数据集
# geo_download函数会自动下载数据并解析为包含表达矩阵(exp)和样本信息(pd)的列表
geo1 = geo_download(datasets[1])
geo2 = geo_download(datasets[2])
geo3 = geo_download(datasets[3])

根据各数据集表达矩阵的数值范围，判断是否需要进行对数转换。通常，如果表达值范围较大（跨数量级，0~几万、几十万之间），则说明数据未经log2转换。我们以0-20的范围作为判断标准，对超出此范围的数据进行log2(x+1)转换，用于差异分析。

# 检查GSE19826的表达谱范围，发现数值较大，需进行log2转换
# range(geo1$exp) 
exp1 = log2(geo1$exp + 1)

# 其他两个数据集的表达值已在合理范围内，无需转换
exp2 = geo2$exp
exp3 = geo3$exp

2. 分组信息与探针注释

为每个数据集创建分组信息，明确区分“癌症”与“非癌”样本，并设置为因子类型，确保后续分析中“对照组在前，疾病组在后”。

library(stringr)


Group1 = ifelse(str_detect(geo1$pd$`tissue type:ch1`, "noncancer"), "noncancer", "cancer")
Group1 = factor(Group1, levels = c("noncancer", "cancer"))
table(Group1)

## Group1
## noncancer    cancer 
##        12        15

Group2 = ifelse(str_detect(geo2$pd$title, "normal"), "normal", "cancer")
Group2 = factor(Group2, levels = c("normal", "cancer"))
table(Group2)

## Group2
## normal cancer 
##     21    111

Group3 = ifelse(str_detect(geo3$pd$source_name_ch1, "mucosa"), "normal", "cancer")
Group3 = factor(Group3, levels = c("normal", "cancer"))
table(Group3)

## Group3
## normal cancer 
##     10     10

获取芯片平台的探针注释信息，建立探针ID与基因符号（Gene Symbol）的对应关系。

注意：本例中3个数据集均使用GPL570 [hgu133plus2]平台，因此可共用一套注释信息。在实际应用中，若数据集来自不同平台（GPL号不同），需要分别为每个平台获取注释，用在下面的get_deg函数中。

library(hgu133plus2.db)
ids <- toTable(hgu133plus2SYMBOL)
head(ids)

##    probe_id symbol
## 1 1007_s_at   DDR1
## 2   1053_at   RFC2
## 3    117_at  HSPA6
## 4    121_at   PAX8
## 5 1255_g_at GUCA1A
## 6   1294_at   UBA7

3. 差异基因筛选与差异基因提取

使用tinyarray包中的get_deg函数，对3个数据集分别进行差异表达分析。该函数整合了limma包的功能，能快速得到差异基因结果。

deg1 = get_deg(exp1, Group1, ids, entriz = F)
head(deg1)

##       logFC  AveExpr         t      P.Value adj.P.Val
## 1 -2.380063 3.865782 -5.603020 5.505999e-06 0.1548921
## 2 -1.025598 5.519599 -5.310656 1.216449e-05 0.2216978
## 3 -1.995034 4.485646 -4.780629 5.148657e-05 0.3332432
## 4  1.846215 5.252243  4.720927 6.057212e-05 0.3332432
## 5 -1.911447 3.387338 -4.633209 7.689642e-05 0.3332432
## 6 -2.101307 3.790645 -4.633087 7.692194e-05 0.3332432
##          B   probe_id       symbol change
## 1 3.044180 1555018_at        OR2C3   down
## 2 2.464743 1554793_at        UBE3C   down
## 3 1.393080  229094_at ATP6V0E2-AS1   down
## 4 1.271127 1554810_at      PLA2G4C     up
## 5 1.091645 1560557_at LOC105370549   down
## 6 1.091395 1562973_at LOC101928201   down

deg2 = get_deg(exp2, Group2, ids, entriz = F)
deg3 = get_deg(exp3, Group3, ids, entriz = F)

为使代码更简洁、可读性更强，我们自定义up和down函数，用于从差异分析结果中快速提取上调和下调的基因列表。

# 提取上调基因的函数
up = function(deg) {
  deg$symbol[deg$change == "up"]
}

# 提取下调基因的函数
down = function(deg) {
  deg$symbol[deg$change == "down"]
}

将各数据集的上调基因和下调基因分别整合到列表中（uplist和downlist）。get_deg的输出默认按p值排序，这恰好满足RRA算法对输入数据的要求（输入为排序好的基因列表）。

uplist = list(up(deg1), up(deg2), up(deg3))
downlist = list(down(deg1), down(deg2), down(deg3))

4. RRA整合

RRA算法背景介绍

在生物信息学分析中，单个研究的差异表达基因列表往往存在假阳性，且不同研究之间的结果重合度可能不高。这主要是由于样本异质性、批次效应、样本量不足或采用的分析平台不同等因素造成的。为了克服这些局限性，好的算法对多个数据集差异分析结果的整合至关重要。

稳健排序聚合（Robust Rank Aggregation, RRA） 是一种强大的Meta分析方法，专门用于整合多个排序后的基因列表（或其他元素列表）。它的核心思想是：如果一个基因在多个独立的分析中都 consistently（一致地）排在靠前的位置（即差异最显著），那么它是一个真正差异基因的可能性就远大于偶然。RRA算法会为每个基因计算一个Score值（类似于p值），这个值反映了其排序位置在多个列表中的一致性程度是否显著优于随机排序。通过这种方式，RRA能够筛选出在多个数据集中都表现稳健的候选基因，显著提升结果的可靠性。

执行RRA整合

我们使用RobustRankAggreg包，分别对上调和下调基因列表进行RRA整合。同时，我们统计每个基因在所有数据集中出现的频率，并将其合并到RRA结果中，作为额外筛选标准。

library(RobustRankAggreg)

# 对上调基因列表进行RRA整合
ups = aggregateRanks(glist = uplist, N = length(unique(unlist(uplist))))
# 计算每个基因出现的频率
freq1 = as.data.frame(table(unlist(uplist)))
# 将频率信息添加到RRA结果中
ups$Freq = freq1[match(ups$Name, freq1$Var1), 2]
head(ups)

##            Name        Score Freq
## COL10A1 COL10A1 0.0003280513    3
## SALL4     SALL4 0.0011799491    3
## THY1       THY1 0.0029836957    3
## BGN         BGN 0.0037554172    3
## IGFBP4   IGFBP4 0.0044134300    3
## COL5A2   COL5A2 0.0046833500    3

# 对下调基因列表进行RRA整合
downs = aggregateRanks(glist = downlist, N = length(unique(unlist(downlist))))
# 计算频率并添加
freq2 = as.data.frame(table(unlist(downlist)))
downs$Freq = freq2[match(downs$Name, freq2$Var1), 2]
head(downs)

##          Name        Score Freq
## VSTM2A VSTM2A 0.0001592915    3
## ERLN     ERLN 0.0006647707    3
## PSAPL1 PSAPL1 0.0006774730    3
## SMPD3   SMPD3 0.0013483127    3
## ADH7     ADH7 0.0018878245    3
## DISP1   DISP1 0.0026651509    3

基于RRA整合结果筛选核心差异基因。这里我们设定筛选条件为：Score < 0.05（统计上显著）且 Freq > 1（至少在2个数据集中都被识别为差异基因）。这些阈值可根据实际需求进行调整。

up_gene = ups$Name[ups$Score < 0.05 & ups$Freq > 1]
length(up_gene)

## [1] 35

down_gene = downs$Name[downs$Score < 0.05 & downs$Freq > 1]
length(down_gene)

## [1] 86

5. 核心差异基因的可视化

为了直观验证这些核心差异基因在各个数据集中的表达变化趋势是否一致，我们提取它们在3个数据集中的logFC值，构建矩阵并通过热图进行可视化。

# 构建核心上调基因的logFC矩阵
updat = data.frame(
  deg1$logFC[match(up_gene, deg1$symbol)],
  deg2$logFC[match(up_gene, deg2$symbol)],
  deg3$logFC[match(up_gene, deg3$symbol)],
  row.names = up_gene
)
colnames(updat) = datasets

# 构建核心下调基因的logFC矩阵
downdat = data.frame(
  deg1$logFC[match(down_gene, deg1$symbol)],
  deg2$logFC[match(down_gene, deg2$symbol)],
  deg3$logFC[match(down_gene, deg3$symbol)],
  row.names = down_gene
)
colnames(downdat) = datasets

为了使热图清晰，我们选取RRA分析结果中最显著（排序最靠前）的10个上调基因和10个下调基因进行展示。

# 合并logFC矩阵
dat = rbind(updat[1:10, ], downdat[1:10, ])
dat

##           GSE19826    GSE54129  GSE79973
## COL10A1   3.071534  2.15271981  4.462640
## SALL4     2.607545  1.18278826  3.422764
## THY1      1.342529  2.89722251  2.563030
## BGN       1.709515  3.00616305  2.370870
## IGFBP4    1.106299  3.39094842  1.004140
## COL5A2    1.098713  1.55004533  1.951849
## SPARC     1.186296  2.02695018  2.062874
## CDH3      1.905570  1.17492855  3.227229
## COL6A3    1.196825  2.59742060  2.058583
## NID2      1.287962  1.85254663  2.278844
## VSTM2A   -1.490487 -1.99546404 -3.848826
## ERLN     -1.111827 -3.80918861 -2.682419
## PSAPL1   -2.216413 -3.11882447 -4.934091
## SMPD3    -1.018454 -2.53184664 -1.478866
## ADH7     -1.448377 -3.21195574 -4.268211
## DISP1    -1.027092 -1.77041055 -1.424939
## ACER2    -1.792120 -4.42151077 -2.999481
## SULT1B1  -1.246034 -2.12849222 -2.754854
## DCAF12L1 -2.277390  0.03808734 -3.329189
## RDH12    -2.105713 -2.51700831 -3.622408

# 使用pheatmap绘制热图
pheatmap::pheatmap(
  dat,
  cluster_cols = F,                   # 不对列（数据集）进行聚类
  display_numbers = TRUE,             # 在热图上显示logFC数值
  number_format = "%.2f",             # 数值格式化为两位小数
  color = colorRampPalette(c("navy", "white", "firebrick3"))(50) # 设置颜色
)