多个数据集整合神器RobustRankAggreg包的示例分析

RobustRankAggreg包在各种数据挖掘文章里面亮相的频次之高，无需我多言，大家可以去查看一下引用它的文献，基本上都是GEO数据库挖掘文章：

比如发表在peerJ的BIOINFORMATICS AND GENOMICS的文章Identification of key candidate genes and biological pathways in bladder cancer 里面的：

作者把这4个数据集，分别独立走差异分析，火山图，热图等等标准流程，基本上读一下我在生信技能树的表达芯片的公共数据库挖掘系列推文 就明白了；

• 解读GEO数据存放规律及下载，一文就够
• 解读SRA数据库规律一文就够
• 从GEO数据库下载得到表达矩阵 一文就够
• GSEA分析一文就够（单机版+R语言版）
• 根据分组信息做差异分析- 这个一文不够的
• 差异分析得到的结果注释一文就够

你也可以很轻松的分析这几个数据集：GSE7476, GSE13507, GSE37815 and GSE65635 ，然后作者就使用了RobustRankAggreg包对这4个数据集的差异分析结果进行整合，如下：

• The integrated DEGs were screened using the RRA package (corrected     P < 0.05, logFC > 1 or −logFC < −1).
• The RRA method is based on the assumption that each gene in each dataset is randomly arranged.
• If the gene ranks high in all datasets, the associated     P-value is lower, the possibility of differential gene expression is greater.
• Through rank analysis, 343 integrated DEGs, consisting of 111 upregulated genes and 232 downregulated genes, were identified by the RRA method

并且把top20的上调基因和下调基因的差异倍数进行热图可视化，如下：

当然了，不仅仅是mRNA的表达芯片，其它，比如circRNA芯片也是如此，同样是发表于2018的文章：A circRNA–miRNA–mRNA network identification for exploring underlying pathogenesis and therapy strategy of hepatocellular carcinoma

就是下载了3个GEO数据集，走差异分析，并且使用RobustRankAggreg包进行整合，最后仅仅是确定了6个circRNA。

几百篇文章我们就不用一一解读啦，反正都是独立的数据集自己做自己的差异分析，然后把多个数据集的差异基因拿去使用RobustRankAggreg包进行整合。

RobustRankAggreg包说明书

这个RobustRankAggreg包超级简单，有意思的是居然并不在bioconductor列表哦，可能是因为它最开始并不是为生物信息学领域的数据分析而创造的吧！因为不在bioconductor，所以它的示例教程一塌糊涂，需要一点背景才能理解。其重点就是aggregateRanks函数而已：

options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
# https://bioconductor.org/packages/release/bioc/html/GEOquery.html
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("RobustRankAggreg",ask = F,update = F) 
library(RobustRankAggreg)
?aggregateRanks

一般来说，正常R包的函数，都是可以通过问号来调取其帮助文档的，aggregateRanks函数也不例外。我们直接看一下示例代码：

set.seed(1234567)
glist <- list(sample(letters, 4), sample(letters, 10), sample(letters, 12))
freq=as.data.frame(table(unlist(glist)))
# Aggregate the inputs
ag=aggregateRanks(glist = glist, N = length(letters))
ag$Freq=freq[match(ag$Name,freq$Var1),2]

的确是超级简单，可以看到，我们有26个字母，假设是26个基因，然后做了3次随机抽样，假设是3个数据集的差异分析，拿到的上调基因，列表如下：

值得注意的是，每次抽样，得到的字母列表的顺序也是有意义的哦。我们的多次数据集差异分析结果，也制作成为这样的表格即可哈！

然后直接使用aggregateRanks函数即可，得到的数据结果如下：

可以看到，a这个字母在3次随机抽样都抽到了，所以它的 exact p-value  非常小，就是统计学非常显著啦！

然后，其余的出现了两次的字母就比较多了，它们的得分之所以有区别，就在于它们的排序。

• n和g都是出现两次，而且排名很靠前，所以p值是0.19，马马虎虎
• k出现了两次，q出现一次，而且都有一个在各自的抽样场合排名第一，k的另外一次在最后面所以权重很低，所以p值是0.33，很差了。
• 至于e和y，虽然也是出现了两次，但是都排名超级靠后，所以p值也很辣鸡，接近于1了。

 aggregateRanks函数其实就是对多个排好序的基因集，进行求交集的同时还考虑一下它们的排序情况。总体上来说，就是挑选那些在多个数据集都表现差异的基因，并且每次差异都排名靠前的那些。