本教程問題：
1.mmGO.rdata文件的制作，是否不需要這一步，更改一下腳本即可實(shí)現(xiàn)呢？
2.這樣的圖形美化，目前的圖還不是很滿意，需要進(jìn)行美化。
這些問題，需要我們大家一起解決，切分享哦！
一、前言
我們在做組學(xué)時，常常遇到多多個處理進(jìn)行差異分析，獲得差異基因或差異代謝物。此后，是對這些差異結(jié)果進(jìn)行GO或KEGG富集分析。這是一個很普遍的結(jié)果，但是這樣一弄后，我們對多個組合差異基因進(jìn)行富集分析后，獲得多個GO或KEGG富集分析結(jié)果，圖很多，但是重復(fù)的term也是很多，那如何闡述就是個問題，以及很多圖列在論文中也不是很美觀，只能間部分圖放在附件中。
那么，我們是不是可以間多個GO或KEGG的term進(jìn)行合并繪圖呢？這個想法在很久以前就像做，但是一直沒找到相關(guān)的教程。今天，突發(fā)奇想，無意間看到這樣的教程，那就重復(fù)一下吧。
PS：在這個教程中，自己很多的參數(shù)還是沒有很好的了解，如果你看過或是已經(jīng)可以深入了解，希望你可以分享一下哦。
二、圖形效果
這是最終圖形的效果，總體來說，還是基本滿足我們的需求，但是很多細(xì)節(jié)還是需要調(diào)整。
三、開始繪圖
3.1 導(dǎo)入相關(guān)的包
## 導(dǎo)入R包
library(plyr)
library(stringr)
library(ape)
library(GOSemSim)
library(ggtree) ## 進(jìn)化樹
library(scales)
library(cowplot) ## 合圖
library(ggplot2)
設(shè)置路勁
setwd("D:\\小杜的生信筆記\\20220605GO富集聚類")
導(dǎo)入GO富集結(jié)果，這些結(jié)果在外面前期的教程中可獲得：clusterProfiler包 |GO、KEGG功能富集分析 | 值得收藏 ，GO、KEGG功能富集分析 | 功能富集網(wǎng)絡(luò)圖、熱圖繪制（代碼重新）
多組富集分析結(jié)果的文件名統(tǒng)一命名，便于導(dǎo)入
fnames <- Sys.glob("enrichGO*.csv")
fnames
# GO term的合集
ego.ID <- unique(ego.all[,c(2:4)])
head(ego.ID)
> head(ego.ID)
ID Description BgRatio Bg
1 GO:0046777 protein autophosphorylation 235/23239 235
2 GO:0034329 cell junction assembly 200/23239 200
3 GO:0034330 cell junction organization 248/23239 248
5 GO:0098742 cell-cell adhesion via plasma-membrane adhesion molecules 199/23239 199
11 GO:0034332 adherens junction organization 121/23239 121
12 GO:0043405 regulation of MAP kinase activity 293/23239 293
```
## 3.3 GO term過濾
```
#刪掉超過300個基因的GO term
#如果不想刪除，可以注釋掉下面這段
ego.ID$Bg <- as.numeric(str_split_fixed(ego.ID$BgRatio, "/",2)[,1]) #提取GO term包含的基因數(shù)量
ego.ID <- ego.ID[ego.ID$Bg < 300,]
dim(ego.ID)
#刪掉少于100個基因的GO term
#具體怎樣篩選，要根據(jù)你自己的數(shù)據(jù)靈活設(shè)置
ego.ID <- ego.ID[ego.ID$Bg > 100,]
dim(ego.ID)
head(ego.ID)
```
## 3.4 橫向合并多組富集結(jié)果，用于畫熱圖
```
## 橫向合并多組富集結(jié)果，用于畫熱圖
MyMerge <- function(x, y){
df <- merge(x, y, by= "ID", all.x= TRUE, all.y= TRUE)
return(df)
}
ego.m <- Reduce(MyMerge, fdataset)
head(ego.m)
```
```
#只保留GO term ID和各組的p.adjust
ego.m <- ego.m[,c(1,4,7,10)] #此處有三組，如果有更多組，要在后面繼續(xù)加
#提取篩選過的GO term
ego.m <- merge(ego.ID[,1:2], ego.m, by= "ID", all.x= TRUE)
rownames(ego.m) <- ego.m$Description
ego.m$ID <- NULL
ego.m$Description <- NULL
```
```
#把列名改為組名
#colnames(ego.m)<- str_remove(fnames, ".csv")
colnames(ego.m) <- paste0("G", seq(1:length(fnames)))
head(ego.m)
> head(ego.m)
G1 G2 G3
tissue homeostasis NA 0.008364684 0.004178739
regulation of cell-matrix adhesion NA 0.003506443 0.003381501
hematopoietic progenitor cell differentiation 0.007418708 0.002188681 NA
myeloid cell homeostasis 0.012630541 NA NA
myeloid leukocyte differentiation 0.003034455 NA NA
regulation of leukocyte migration 0.005619600 NA 0.012970792
3.4 按GO term的相似性聚類
mmGO.rdata文件是這樣獲得的，但是我在考慮如果我們沒有自己.db是要如何獲得，是否可以根據(jù)我們前面的分析進(jìn)行獲得呢？這個問題需要大家一起解決，如果你知道可以分享一下，謝謝??！*
library(org.Hs.eg.db)
hgGO <- godata('org.Hs.eg.db', ont="BP")
head(hgGO)
#save(hgGO, file="mmGO.rdata")
Q1，需要我們一起解決?。?！
#導(dǎo)入之前保存的文件
(load("mmGO.rdata"))
#計算GO term之間的相似性
ego.sim <- mgoSim(ego.ID$ID, ego.ID$ID, semData=mmGO, measure="Wang", combine=NULL)
ego.sim[1:3, 1:3]
> ego.sim[1:3, 1:3]
protein autophosphorylation cell junction assembly
protein autophosphorylation 1.000 0.061
cell junction assembly 0.061 1.000
cell junction organization 0.090 0.683
cell junction organization
protein autophosphorylation 0.090
cell junction assembly 0.683
cell junction organization 1.000
#用GO term作為行名、列名，便于查看和畫圖
rownames(ego.sim) <- ego.ID$Description
colnames(ego.sim) <- ego.ID$Description
ego.sim[1:3, 1:3]
繪制聚類圖
#給GO term分類
tree <- nj(as.dist(1-ego.sim))
ggtree(tree) + geom_tiplab() + #寫GO term
geom_text2(aes(subset=!isTip, label=node), hjust=-.3) + #寫node編號
coord_cartesian(xlim=c(-.1,1.3)) #左右兩側(cè)留出合適的空間
此圖節(jié)點(diǎn)就是我們后續(xù)分類的標(biāo)準(zhǔn)
3.5 繪圖
繪圖代碼在下面的教程中，如果獲得也可以在知乎、公眾號中進(jìn)行獲得代碼數(shù)據(jù)事例數(shù)據(jù)。
多組差異基因富集分析結(jié)果聚類 | 獲得共表達(dá)term的富集圖
@小杜的生信筆記，主要發(fā)表或收錄生物信息學(xué)的教程，以及基于R的分析和可視化（包括數(shù)據(jù)分析，圖形繪制等）；分享感興趣的文獻(xiàn)和學(xué)習(xí)資料!
