富集分析事後丸:過濾數據

作者：Y叔

來源: biobabble公眾號

（點以上文字會分別跳轉到有意思的原貼上，全體學霸都關注了Y叔哦~）

點本文左下方的閱讀原文，直接向Y叔提問。

有了事後丸之後，腰不酸了，背不疼了，腿啊也不抽筋了，走路也有勁兒了，我有個夢想，希望每個生信小白做分析之後，都能吃上一顆事後丸。

The cnetplot is a really attractive feature. I wonder if it is possible to subset the enrichment object that goes into the cnetplot function?

Case scenario: Run an enrichment analysis with pvalueCutoff = 1, to see all results. Plotting them all would be infeasible. How to subset the enrichment object to, say, first ten most significant terms, and then plot it with cnetplot?

github上這個問題問得挺好，其實啊，你用clusterProfiler系列包做富集分析的話，根本就不用卡p值、q值這些，因為你完全可以吃事後丸！拿到全部結果之後，想怎麼搞就怎麼搞。

基本上data.frame取子集那些方法，包括[, [[, $等操作符，都被我重新定義，也就是說，對於富集分析的結果，你完全可以像data.frame一樣的操作。

這裡用DOSE包做為例子，富集的結果存在x對象中：

> require(DOSE)
Loading required package: DOSE

DOSE v3.9.1  For help: https://guangchuangyu.github.io/DOSE

If you use DOSE in published research, please cite:
Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an R/Bioconductor package for Disease Ontology Semantic and Enrichment analysis. Bioinformatics 2015, 31(4):608-609

> data(geneList)
> de = names(geneList)[1:100]
> x = enrichDO(de, qvalueCutoff=1, pvalueCutoff=1)

我們可以用dim看看有多少結果，可以用head和tail偷窺一下結果：

> dim(x)
[1] 383   9
> head(x, 2)
                       ID            Description GeneRatio  BgRatio
DOID:0060071 DOID:0060071 pre-malignant neoplasm      5/77  22/8007
DOID:5295       DOID:5295     intestinal disease      9/77 157/8007
                   pvalue     p.adjust       qvalue
DOID:0060071 1.671524e-06 0.0006401937 0.0004609887
DOID:5295    1.759049e-05 0.0027885022 0.0020079362
                                                  geneID Count
DOID:0060071                    6280/6278/10232/332/4321     5
DOID:5295    4312/6279/3627/10563/4283/890/366/4902/3620     9

然後我們就可以像操作data.frame一樣來操作它，取個子集，但最後你會發現，輸出的是data.frame的對象y，也就是你實際上無法用enrichplot去對y作圖：

> y = x[x$pvalue < 0.05, ]
> dim(y)
[1] 121   9
> tail(y, 2)
                 ID                    Description GeneRatio  BgRatio
DOID:7474 DOID:7474 malignant pleural mesothelioma      2/77  36/8007
DOID:8692 DOID:8692               myeloid leukemia      4/77 142/8007
              pvalue  p.adjust    qvalue             geneID Count
DOID:7474 0.04659310 0.1487096 0.1070824          10232/332     2
DOID:8692 0.04748286 0.1502970 0.1082254 820/10232/332/3620     4
> class(y)
[1] "data.frame"

解決方法來了

於是我給[方法加了一個參數，叫asis，讓輸出保持為富集分析結果的對象，默認是FALSE，我們顯式設為TRUE即可。這功能要求DOSE >= 1.9.1。

> y = x[x$pvalue < 0.05, asis=T]
> class(y)
[1] "enrichResult"
attr(,"package")
[1] "DOSE"
> y
#
# over-representation test
#
#...@organism    Homo sapiens
#...@ontology    DO
#...@keytype     ENTREZID
#...@gene        chr [1:100] "4312" "8318" "10874" "55143" "55388" "991" "6280" "2305" ...
#...pvalues adjusted by \'BH\' with cutoff <1
#...121 enriched terms found
\'data.frame\':   121 obs. of  9 variables:
 $ ID         : chr  "DOID:0060071" "DOID:5295" "DOID:8719" "DOID:3007" ...
 $ Description: chr  "pre-malignant neoplasm" "intestinal disease" "in situ carcinoma" "breast ductal carcinoma" ...
 $ GeneRatio  : chr  "5/77" "9/77" "4/77" "4/77" ...
 $ BgRatio    : chr  "22/8007" "157/8007" "18/8007" "29/8007" ...
 $ pvalue     : num  1.67e-06 1.76e-05 2.18e-05 1.56e-04 2.08e-04 ...
 $ p.adjust   : num  0.00064 0.00279 0.00279 0.0136 0.0136 ...
 $ qvalue     : num  0.000461 0.002008 0.002008 0.009796 0.009796 ...
 $ geneID     : chr  "6280/6278/10232/332/4321" "4312/6279/3627/10563/4283/890/366/4902/3620" "6280/6278/10232/332" "6280/6279/4751/6286" ...
 $ Count      : int  5 9 4 4 13 6 13 5 5 6 ...
#...Citation
  Guangchuang Yu, Li-Gen Wang, Guang-Rong Yan, Qing-Yu He. DOSE: an
  R/Bioconductor package for Disease Ontology Semantic and Enrichment
  analysis. Bioinformatics 2015, 31(4):608-609

於是新的y，就還是富集分析生成的對象，這樣就可以用enrichplot來畫圖。目前只在開發版支持，超幾何分布和GSEA都全線支持。以後就有事後丸吃了！

這裡演示用的雖然只是pvalue，我知道很多人只會follow例子，所以要特別指出，任何column，你都可以拿來濾，比如geneID你可以濾出只有某些基因參與的通路出來，比如ID或Description你可以濾出你想要的通路來，而Count你可以限制參與的基因的數目等等，濾完之後，你就可以畫出相應的圖了，這簡直給了你玩壞它的可能性。

PS1: 微軟越來越硬了，一則消息炸開了鍋：

Unlimited free private Git repositories

從此github晉升網盤，我一直在用gitlab當網盤，現在可以用github了。

還記得《製作meme的通用方式，來了解一下》一文中的調侃嗎？雖然目前還是gitlab給得大方些。這次github可是要死掐gitlab了，畢竟同性社交嘛，社交的地方，還是要看人氣。其實gitlab的雲存儲，用的也是微軟家的，以我個人的經驗，gitlab當機的概率要高於github。

100年後看，蘋果 vs 微軟，微軟才是偉大的公司，後賈伯斯時候的蘋果只知道賺錢。100年後再看，Steve Jobs vs Bill Gates，Gates才是偉大的人，因為Gates身體力行在做慈善。捐錢容易，身體力行在推動一些東西，就難得了。

PS2: 月底要去南丹麥大學，有沒有丹麥的小夥伴要面基？

PS3: 2019年立個flag，我要寫一本書，掃碼圍觀：

往期精彩

讚賞

手機端點閱讀原文，在知識星球app向Y叔提問