【R>>Robust Rank Aggregation】多个数据

挖掘GEO芯片时会发现多个类似的研究，那么怎么把它们整合到一起呢？常见的分析方法：

• 差异分析→直接求交集

PMID: 35313857

• 去除批次效应

PMID:33557837

• Robust Rank Aggregation（RRA)

PMID: 33671634

下面就来学习下RRA方法，堪称多数据集合并利器。

RRA的核心思想：对排序好的基因求交集的时，还考虑基因的排序。讲人话就是，挑选那些在多个数据集都表现差异的基因，并且每次差异都排名靠前的那些。

下面就结合生信技能树健明大神的帖子演练一番：

#更换国内镜像
# options("repos"= c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")

#安装RobustRankAggreg
# if (!require("BiocManager", quietly = TRUE))
#   install.packages("BiocManager")
# BiocManager::install("RobustRankAggreg",ask = F,update=F)

上面是更换国内镜像及安装RobustRankAggreg，下面进入正题：

rm(list = ls())
library(RobustRankAggreg)
set.seed(1234)
glist <- list(sample(letters, 6), sample(letters, 17), sample(letters, 25))
freq=as.data.frame(table(unlist(glist)))

# 我们假设排在前面的基因fold change比较大
ag <- aggregateRanks(glist)
ag$Freq <- freq[match(ag$Name,freq$Var1),2]
ag

# 小试牛刀
library(clusterProfiler)
set.seed(1234)
data(geneList,package = "DOSE")
deg <- data.frame(gene=names(geneList),
                  logFC=as.numeric(geneList),
                  P.Value=0)

#模拟差异的基因列表
deg2=deg;deg2$gene=sample(deg$gene,length(deg2$gene))
deg3=deg;deg3$gene=sample(deg$gene,length(deg2$gene))
deg4=deg;deg4$gene=sample(deg$gene,length(deg2$gene))


#提取高表达基因
get_up <- function(df){
  df$g=ifelse(df$P.Value>0.05,'stable',
              ifelse( df$logFC >1,'up', 
                      ifelse( df$logFC < -1,'down','stable')))
  print( table(df$g))
  df=df[order(df$logFC,decreasing = T),]
  # rownames(df[df$g=='up',])
  df[df$g=='up','gene']
}

glist <- list(get_up(deg2),
               get_up(deg3),
               get_up(deg4))
ups <- aggregateRanks(glist,N=length(unique(unlist(glist))))
tmp <- as.data.frame(table(unlist(glist)))
ups$freq <- tmp[match(ups$Name,tmp[,1]),2]

#获得目标基因(gene significant)
gs <- ups[ups$Score<0.05,1]
updat=data.frame(deg2=deg2[deg2$gene %in% gs,'logFC'],
                 deg3=deg3[deg3$gene %in% gs,'logFC'],
                 deg4=deg4[deg4$gene %in% gs,'logFC'] )
rownames(updat) <- gs

#热图
library(pheatmap)
pheatmap(updat,display_numbers = T)

一顿操作猛入虎后，就得到了下面这张小热图，是不是有点内味儿啦。以后有类似的数据，可以直接替换deg2-deg4，就可以得到类似的结果啦，感觉还是蛮实用的。

声明：本文仅供学习交流用，禁止用于商业用途，如有侵权，请联系删除。
参考链接：
知乎：多个数据集整合神器-RobustRankAggreg包（曾建明）