rm(list = ls())
library(clusterProfiler)
library(org.At.tair.db)
library(stringr)
library(ggplot2)
load("genes.Rdata")
head(genes)
## [1] "ATBCA3" "BCA3" "AtWSCP" "Kunitz-PI;1" "BCAT6"
## [6] "OSR1"
1.ID轉(zhuǎn)換
首先是要轉(zhuǎn)換ID,擬南芥是要TAIR id來做富集分析。怎么知道的呢诵盼？ 如果拿ENTRIZID來做九昧，會(huì)收到這樣式的報(bào)錯(cuò)信息：
那我們就知道應(yīng)該用AT開頭的ID了擅腰。那是啥ID呢，查查唄傲诵。
x = bitr(genes,fromType = "SYMBOL",toType = "TAIR",
OrgDb = "org.At.tair.db")
head(x)
## SYMBOL TAIR
## 1 ATBCA3 AT1G23730
## 2 BCA3 AT1G23730
## 3 AtWSCP AT1G72290
## 4 Kunitz-PI;1 AT1G72290
## 5 BCAT6 AT1G50110
## 6 OSR1 AT2G41230
2.做KEGG富集
ekk <- enrichKEGG(gene = x$TAIR,organism = 'ath')
ekk <- setReadable(ekk,OrgDb = org.At.tair.db,keyType = "TAIR")
#如果ekk是空的凯砍，這句就會(huì)報(bào)錯(cuò)，因?yàn)闆]富集到任何通路拴竹。
# 條帶圖畫一下
barplot(ekk)
發(fā)現(xiàn)縱坐標(biāo)有多余的東西悟衩，可以在enrichResult對象里面刪掉。
代碼里的兩個(gè)斜杠代表括號就是括號本身栓拜，不加兩個(gè)斜杠就會(huì)被當(dāng)作正則表達(dá)式的括號座泳，它另有含義惠昔。
ekk@result$Description = str_remove(ekk@result$Description," - Arabidopsis thaliana \\(thale cress\\)")
barplot(ekk)
改完再畫,這回好咯
3.GO富集分析
就比較簡單，指定一下keyType參數(shù)即可
ego <- enrichGO(gene = genes,OrgDb= org.At.tair.db,
keyType = "SYMBOL",
ont = "ALL")
barplot(ego, split = "ONTOLOGY") +
facet_grid(ONTOLOGY ~ ., space = "free_y",scales = "free_y")
