多個數(shù)據(jù)集整合神器RobustRankAggreg包的示例分析

RobustRankAggreg包在各種數(shù)據(jù)挖掘文章里面亮相的頻次之高，無需我多言，大家可以去查看一下引用它的文獻，基本上都是GEO數(shù)據(jù)庫挖掘文章：

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

作者把這4個數(shù)據(jù)集，分別獨立走差異分析，火山圖，熱圖等等標準流程，基本上讀一下我在生信技能樹的表達芯片的公共數(shù)據(jù)庫挖掘系列推文 就明白了；

• 解讀GEO數(shù)據(jù)存放規(guī)律及下載，一文就夠
• 解讀SRA數(shù)據(jù)庫規(guī)律一文就夠
• 從GEO數(shù)據(jù)庫下載得到表達矩陣 一文就夠
• GSEA分析一文就夠（單機版+R語言版）
• 根據(jù)分組信息做差異分析- 這個一文不夠的
• 差異分析得到的結(jié)果注釋一文就夠

你也可以很輕松的分析這幾個數(shù)據(jù)集：GSE7476, GSE13507, GSE37815 and GSE65635 ，然后作者就使用了RobustRankAggreg包對這4個數(shù)據(jù)集的差異分析結(jié)果進行整合，如下：

• 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的上調(diào)基因和下調(diào)基因的差異倍數(shù)進行熱圖可視化，如下：

當(dāng)然了，不僅僅是mRNA的表達芯片，其它，比如circRNA芯片也是如此，同樣是發(fā)表于2018的文章：A circRNA–miRNA–mRNA network identification for exploring underlying pathogenesis and therapy strategy of hepatocellular carcinoma

就是下載了3個GEO數(shù)據(jù)集，走差異分析，并且使用RobustRankAggreg包進行整合，最后僅僅是確定了6個circRNA。

幾百篇文章我們就不用一一解讀啦，反正都是獨立的數(shù)據(jù)集自己做自己的差異分析，然后把多個數(shù)據(jù)集的差異基因拿去使用RobustRankAggreg包進行整合。

RobustRankAggreg包說明書

這個RobustRankAggreg包超級簡單，有意思的是居然并不在bioconductor列表哦，可能是因為它最開始并不是為生物信息學(xué)領(lǐng)域的數(shù)據(jù)分析而創(chuàng)造的吧！因為不在bioconductor，所以它的示例教程一塌糊涂，需要一點背景才能理解。其重點就是aggregateRanks函數(shù)而已：

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包的函數(shù)，都是可以通過問號來調(diào)取其幫助文檔的，aggregateRanks函數(shù)也不例外。我們直接看一下示例代碼：

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個字母，假設(shè)是26個基因，然后做了3次隨機抽樣，假設(shè)是3個數(shù)據(jù)集的差異分析，拿到的上調(diào)基因，列表如下：

值得注意的是，每次抽樣，得到的字母列表的順序也是有意義的哦。我們的多次數(shù)據(jù)集差異分析結(jié)果，也制作成為這樣的表格即可哈！

然后直接使用aggregateRanks函數(shù)即可，得到的數(shù)據(jù)結(jié)果如下：

可以看到，a這個字母在3次隨機抽樣都抽到了，所以它的 exact p-value  非常小，就是統(tǒng)計學(xué)非常顯著啦！

然后，其余的出現(xiàn)了兩次的字母就比較多了，它們的得分之所以有區(qū)別，就在于它們的排序。

• n和g都是出現(xiàn)兩次，而且排名很靠前，所以p值是0.19，馬馬虎虎
• k出現(xiàn)了兩次，q出現(xiàn)一次，而且都有一個在各自的抽樣場合排名第一，k的另外一次在最后面所以權(quán)重很低，所以p值是0.33，很差了。
• 至于e和y，雖然也是出現(xiàn)了兩次，但是都排名超級靠后，所以p值也很辣雞，接近于1了。

 aggregateRanks函數(shù)其實就是對多個排好序的基因集，進行求交集的同時還考慮一下它們的排序情況?？傮w上來說，就是挑選那些在多個數(shù)據(jù)集都表現(xiàn)差異的基因，并且每次差異都排名靠前的那些。