本教程問題:
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)行獲得呢权逗?這個問題需要大家一起解決美尸,如果你知道可以分享一下,謝謝U遛薄师坎!*
library(org.Hs.eg.db)
hgGO <- godata('org.Hs.eg.db', ont="BP")
head(hgGO)
#save(hgGO, file="mmGO.rdata")
Q1,需要我們一起解決1伎选R倌汀!
#導(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í)資料!