【原】使用clusterProfiler進(jìn)行GO富集分析

生信修煉手冊 2019-12-24

展開全文

clusterProfiler是一個(gè)功能強(qiáng)大的R包，同時(shí)支持GO和KEGG的富集分析，而且可視化功能非常的優(yōu)秀，本章主要介紹利用這個(gè)R包來進(jìn)行Gene Ontology的富集分析。

進(jìn)行GO分析時(shí)，需要考慮的一個(gè)基礎(chǔ)因素就是基因的GO注釋信息從何處獲取。Bioconductor上提供了以下19個(gè)物種的Org類型的包，包含了這些物種的GO注釋信息

packages	organism
org.Ag.eg.db	Anopheles
org.At.tair.db	Arabidopsis
org.Bt.eg.db	Bovine
org.Ce.eg.db	Worm
org.Cf.eg.db	Canine
org.Dm.eg.db	Fly
org.Dr.eg.db	Zebrafish
org.EcK12.eg.db	E coli strain K12
org.EcSakai.eg.db	E coli strain Sakai
org.Gg.eg.db	Chicken
org.Hs.eg.db	Human
org.Mm.eg.db	Mouse
org.Mmu.eg.db	Rhesus
org.Pf.plasmo.db	Malaria
org.Pt.eg.db	Chimp
org.Rn.eg.db	Rat
org.Sc.sgd.db	Yeast
org.Ss.eg.db	Pig
org.Xl.eg.db	Xenopus

clusterProfiler在做GO分析時(shí)，會讀取這些包中的GO信息，從這里也可以看出，Bioconductor的生態(tài)體系，不同R包的數(shù)據(jù)傳遞和依賴組成了其生態(tài)環(huán)境。

對于以上19個(gè)物種，只需要安裝對應(yīng)的org包，clusterProfile就會自動(dòng)從中獲取GO注釋信息，我們只需要差異基因的列表就可以了，使用起來非常方便。

在進(jìn)行富集分析時(shí)，支持以下兩種算法

1. Over-Representation Analysis

過表達(dá)分析其實(shí)就是費(fèi)舍爾精確檢驗(yàn)，對于以上19個(gè)物種，分析的代碼如下

ego <- enrichGO(
  gene          = gene,
  keyType = "ENTREZID",
  OrgDb         = org.Hs.eg.db,
  ont           = "CC",
  pAdjustMethod = "BH",
  pvalueCutoff  = 0.01,
  qvalueCutoff  = 0.05,
  readable      = TRUE)

gene就是差異基因?qū)?yīng)的向量，keyType指定基因ID的類型，默認(rèn)為ENTREZID, 該參數(shù)的取值可以參考keytypes(org.Hs.eg.db)的結(jié)果，建議采用ENTREZID, OrgDb指定該物種對應(yīng)的org包的名字，ont代表GO的3大類別，BP, CC, MF; pAdjustMethod指定多重假設(shè)檢驗(yàn)矯正的方法，cufoff指定對應(yīng)的閾值，readable=TRUE代表將基因ID轉(zhuǎn)換為gene symbol。

2. Gene Set Enrichment Analysis

ego <- gseGO(
  geneList     = geneList,
  OrgDb        = org.Hs.eg.db,
  ont          = "CC",
  nPerm        = 1000,
  minGSSize    = 100,
  maxGSSize    = 500,
  pvalueCutoff = 0.05,
  verbose      = FALSE)

GSEA分析通過置換檢驗(yàn)來計(jì)算p值，nPerm指定置換次數(shù)。

對于這19個(gè)物種之外的其他物種，也支持讀取基因的GO注釋文件，然后進(jìn)行分析，注釋文件的格式如下

GeneId	GO	Description
1	GO:0005819	spindle
2	GO:0072686	mitotic spindle
3	GO:0000776	kinetochore

只需要3列信息即可，第一列為geneID, 第二列為基因?qū)?yīng)的GO編號，第三列為GO的描述信息。這3列的順序是無所謂的，只要包含這3種信息就可以了。
讀取該文件，進(jìn)行分析的代碼如下

data <- read.table(
  "go_annotation.txt",
  header = T,
  sep = "\t")
go2gene <- data[, c(2, 1)]
go2name <- data[, c(2, 3)]
# 費(fèi)舍爾精確檢驗(yàn)
x <- enricher(
gene,
TERM2GENE = go2gene,
TERM2NAME = go2name)
# GSEA富集分析
x <- GSEA(
gene,
TERM2GENE = go2gene,
TERM2NAME = go2name)

對于GO富集分析的結(jié)果，clusterProfiler提供了以下幾種可視化策略

1. barplot

用散點(diǎn)圖展示富集到的GO terms,用法如下

barplot(ego, showCategory = 10)

生成的圖片如下

橫軸為該GO term下的差異基因個(gè)數(shù)，縱軸為富集到的GO Terms的描述信息， showCategory指定展示的GO Terms的個(gè)數(shù)，默認(rèn)展示顯著富集的top10個(gè)，即p.adjust最小的10個(gè)。注意的顏色對應(yīng)p.adjust值，從小到大，對應(yīng)藍(lán)色到紅色。

2. dotplot

用散點(diǎn)圖展示富集到的GO terms,用法如下

dotplot(ego, showCategory = 10)

生成的圖片如下

橫軸為GeneRatio, 代表該GO term下的差異基因個(gè)數(shù)占差異基因總數(shù)的比例，縱軸為富集到的GO Terms的描述信息， showCategory指定展示的GO Terms的個(gè)數(shù)，默認(rèn)展示顯著富集的top10個(gè)，即p.adjust最小的10個(gè)。

圖中點(diǎn)的顏色對應(yīng)p.adjust的值，從小到大，對應(yīng)藍(lán)色到紅色，大小對應(yīng)該GO terms下的差異基因個(gè)數(shù)，個(gè)數(shù)越多，點(diǎn)越大。

3. GO有向無環(huán)圖

調(diào)用topGO來實(shí)現(xiàn)GO有向無環(huán)圖的繪制，代碼如下

plotGOgraph(ego)

生成的圖片如下

矩形代表富集到的top10個(gè)GO terms, 顏色從黃色過濾到紅色，對應(yīng)p值從大到小。

4. goplot

和plotGOgraph展示的信息一樣，都是GO有向無環(huán)圖，只不過采用了igraph中的布局方式，用法如下

goplot(ego)

生成的圖片如下

5. emapplot

對于富集到的GO terms之間的基因重疊關(guān)系進(jìn)行展示，如果兩個(gè)GO terms系的差異基因存在重疊，說明這兩個(gè)節(jié)點(diǎn)存在overlap關(guān)系，在圖中用線條連接起來，用法如下

emapplot(ego, showCategory = 30)

生成的圖片如下

每個(gè)節(jié)點(diǎn)是一個(gè)富集到的GO term, 默認(rèn)畫top30個(gè)富集到的GO terms, 節(jié)點(diǎn)大小對應(yīng)該GO terms下富集到的差異基因個(gè)數(shù)，節(jié)點(diǎn)的顏色對應(yīng)p.adjust的值，從小到大，對應(yīng)藍(lán)色到紅色。