現(xiàn)在的社會正處于大數(shù)據(jù)時(shí)代烦粒，生物學(xué)研究領(lǐng)域也迎來了大組學(xué)時(shí)代次绘。無論是植物。動物撒遣，還是微生物，都有許多物種被測序管跺。在后續(xù)的個(gè)性化分析過程中义黎，通常會得到龐大的基因數(shù)據(jù)，如何選擇一個(gè)高效的生物功能分析方法顯得尤其重要豁跑×椋基因的功能富集，可以讓你快速地了解所獲得的基因具有怎樣的生物功能艇拍，它們主要集中在哪些生物通路狐蜕。clusterProfiler、GOplot這兩個(gè)工具可以說是在眾多GO富集分析軟件中脫穎而出卸夕，它們不僅能夠支持多個(gè)數(shù)據(jù)庫（GO/KEGG/DOSE/MSigDb/Reactome）层释，還能繪出高大上的CNS級別富集圖。那么快集，今天主要就是介紹如何結(jié)合這兩個(gè)軟件贡羔，讓你的一組gene symbol秒變高富帥。
01. 準(zhǔn)備數(shù)據(jù)
經(jīng)過上游的生信分析我們會獲得許多具有生物學(xué)意義的gene set个初，可以是差異表達(dá)基因乖寒，也可是正選擇基因或者加速進(jìn)化基因。通常院溺，只要具有這些基因的gene symbol或者是geneid楣嘁，都可以利用該軟件進(jìn)行分析。
本次所使用的數(shù)據(jù)主要是gene symbol珍逸，如下圖所示
02. 讀入數(shù)據(jù)
setwd("C:/Users/lenovo/Desktop/RERconverge-master/R") # 設(shè)置工作路徑
geneset <- read.table(file = "genelist",header = F) #讀取數(shù)據(jù)
> head(geneset) #查看數(shù)據(jù)
V1
1 ADCY10
2 ADCY2
3 ADCY3
4 ADCY4
5 ADCY5
6 ADCY6
gene=as.character(geneset[,1]) #轉(zhuǎn)換字符格式
> head(gene)
[1] "ADCY10" "ADCY2" "ADCY3" "ADCY4"
[5] "ADCY5" "ADCY6"
03. 使用clusterProfiler進(jìn)行富集分析
library(org.Hs.eg.db)
library(clusterProfiler)
library(GOplot)
library(ggplot2)
library(stringr)
gene <- bitr(gene, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Hs.eg.db") # 將symbol轉(zhuǎn)換為entrzid
ego <- enrichGO(gene=gene$ENTREZID,
OrgDb='org.Hs.eg.db',
ont= "BP",
pAdjustMethod="BH",
pvalueCutoff = 0.01,
qvalueCutoff = 0.05,
readable = TRUE)
> dim(ego)
[1] 402 9
ego2 <- simplify(ego,cutoff=0.7,by="p.adjust",select_fun=min) #去除冗余逐虚，可以調(diào)整cutoff值
> dim(ego2)
[1] 144 9
heatlpot(ego2)
cnetplot(ego2,categorySize="pvalue",circular = TRUE, colorEdge = TRUE)
# KEGG富集
kk <- enrichKEGG(gene$ENTREZID, organism="hsa",keyType = "kegg",pvalueCutoff=0.01,pAdjustMethod="BH",qvalueCutoff=0.05)
barplot(kk, showCategory=20,title="Enrichment KEGG")
04. 整理GOplot輸入文件格式
###ego3 <- data.frame(ego2) ### 轉(zhuǎn)換為數(shù)據(jù)矩陣
GO <- ego[1:9,c(1,2,8,6)] ### 選取富集結(jié)果前9個(gè)，并提取GO id弄息，description痊班， adj.p以及基因ID列
GO$geneID <- str_replace_all(GO$geneID,"/",",") ### 修改geneID這一列的分隔符號
names(GO)=c("ID","Term","Genes","adj_pval")
GO$Category = "BP"
###構(gòu)建基因矩陣，并隨機(jī)生成logFC
genedata = data.frame(ID=geneset$V1,logFC=rnorm(length(geneset$V1),mean = 0,sd=2))
circ <- circle_dat(GO,genedata)
05. 使用GOplot繪圖
# 01. 和弦圖
chord <- chord_dat(data=circ, genes = genedata) # 生成帶有選定基因列表的矩陣
chord <- chord_dat(data=circ, process = GO$Term) #生成帶有選定GO term的列表矩陣
chord <- chord_dat(data=circ, genes=genedata,process = GO$Term) # 構(gòu)建數(shù)據(jù)
GOChord ( data=chord,
title = "GOChord plot",
space = 0.02, # go term處間隔大小
limit = c(3,5) #第一個(gè)值是至少分配給一個(gè)基因的go term數(shù)目摹量，第二個(gè)數(shù)值是至少分配給一個(gè)Go term的基因數(shù)
gene.order ='logFC',gene.space=0.25,gene.size=10,
lfc.col = c('firebrick3','white','royalblue3') #上下調(diào)基因顏色
ribbon.col=brewer.pal(length(GO$Term)),"set3" #GO term顏色
process.label = 18 # Go terms字體大小
)
# 02. 條形圖
GOBar(circ,display = 'multiple')
# 03. 氣泡圖
# 要添加標(biāo)題涤伐，請更改圓圈的顏色馒胆，對圖進(jìn)行構(gòu)圖，并更改標(biāo)簽閾值凝果，請使用以下參數(shù)：
GOBubble(circ, title = 'Bubble plot', colour = c('orange', 'darkred', 'gold'), display = 'multiple', labels = 3)
# 對于構(gòu)面圖祝迂，還可以通過將bg.col設(shè)置為TRUE，根據(jù)顯示的類別為面板的背景著色:
GOBubble(circ, title = 'Bubble plot with background colour', display = 'multiple', bg.col = T, labels = 3)
# 軟件包的更新版本中包含一個(gè)新函數(shù)reduce_overlap 器净，以減少冗余項(xiàng)的數(shù)量. 到目前為止型雳，已實(shí)現(xiàn)的方法非常簡單+緩慢，需要進(jìn)一步完善. 但是山害，通過減少冗余項(xiàng)的數(shù)量纠俭，可以像氣泡圖一樣顯著提高圖的可讀性. 該函數(shù)刪除基因重疊大于或等于設(shè)定閾值的所有術(shù)語. 該函數(shù)在不考慮GO層次結(jié)構(gòu)的情況下，每組代表一個(gè)術(shù)語:
# Reduce redundant terms with a gene overlap >= 0.75...
reduced_circ <- reduce_overlap(circ, overlap = 0.75)
# ...and plot it
GOBubble(reduced_circ, labels = 2.8)
# 04. 圈圖
GOCircle(circ)
# 外圈顯示了分配基因的logFC的每個(gè)項(xiàng)的散點(diǎn)圖. 默認(rèn)情況下浪慌，紅色圓圈顯示上調(diào)冤荆，藍(lán)色圓圈顯示下調(diào). 可以使用參數(shù)lfc.col更改顏色. 因此，更容易理解权纤，為什么在某些情況下钓简，高度有效的術(shù)語的z得分接近于零. Z分?jǐn)?shù)為零并不意味著該術(shù)語不重要. 至少沒有，只要其顯著豐富即可. 它只是表明z分?jǐn)?shù)是一個(gè)粗略的度量汹想，因?yàn)轱@然分?jǐn)?shù)并未考慮過程中單個(gè)基因的功能水平和激活依賴性. 您可以使用各種參數(shù)來更改圖的布局外邓，請參閱？ GOCirlce.nsub參數(shù)需要更多說明才能明智地使用. 首先古掏，它可以是數(shù)字或字符向量. 如果它是一個(gè)字符向量损话，則它包含要顯示的進(jìn)程的ID或術(shù)語說明（未顯示輸出）。
IDs <- c('GO:0007507', 'GO:0001568', 'GO:0001944', 'GO:0048729', 'GO:0048514', 'GO:0005886', 'GO:0008092', 'GO:0008047')
GOCircle(circ, nsub = IDs)
# 05. 熱圖
GOHeat(chord, nlfc = 1, fill.col = c('red', 'yellow', 'green'))
# 06. 聚類圖
GOCluster(circ, EC$process, clust.by = 'logFC', term.width = 2)
OCluster(circ, EC$process, clust.by = 'term', lfc.col = c('darkgoldenrod1', 'black', 'cyan1'))
# 07. 韋恩圖
l1 <- subset(circ, term == 'heart development', c(genes,logFC))
l2 <- subset(circ, term == 'plasma membrane', c(genes,logFC))
l3 <- subset(circ, term == 'tissue morphogenesis', c(genes,logFC))
GOVenn(l1,l2,l3, label = c('heart development', 'plasma membrane', 'tissue morphogenesis'))
參考鏈接：
- https://yulab-smu.github.io/clusterProfiler-book/index.html
- http://wencke.github.io/
- http://www.360doc.com/content/19/0416/12/52645714_829160785.shtml
- http://www.reibang.com/p/77246ff36214
- https://bioinfohome.com/index.php/2019/07/08/goplot-01/
- http://www.sohu.com/a/339298858_120098972
- http://wap.sciencenet.cn/blog-118204-1192033.html?mobile=1
- http://cran.r-project.org.icopy.site/web/packages/GOplot/vignettes/GOplot_vignette.html
- http://www.reibang.com/p/16fd75aff9d4?utm_campaign=hugo&utm_medium=reader_share&utm_content=note&utm_source=qq
