image
本次教程的figure仍然是讀者求助的圖端铛，算得上是kegg富集圖的新流派幢踏。據(jù)我的調(diào)查，該圖應(yīng)該是基迪奧云-平-臺首創(chuàng)(https://www.omicshare.com/tools/Home/Soft/enrich_circle)匿垄，之后小白魚的生統(tǒng)筆記進行了復(fù)現(xiàn)(仿一個網(wǎng)圖亿乳，使用circlize包繪制圈圖可視化基因集富集分析結(jié)果)。
最開始也是跟著上述的帖子學(xué)習(xí)共螺，之后自己對代碼進行了改寫该肴，重新安排圖形的布局，使之（在我看來）更有意義藐不。另一個改動是增加了kegg pathway的注釋信息匀哄，我在之前的帖子中也提到了如何獲取這個信息秦效，沒有這個信息是畫不了這張圖的。最后拱雏，將代碼模塊化棉安，一行代碼就能出圖底扳。
這張圖如何解讀：
- 不論是
組間比較還是cluster之間的比較铸抑，在得到某組/某cluster的差異基因后，都能進行KEGG富集分析衷模，選取20-35個term畫圖 - 圓圈
從外向內(nèi)看鹊汛，第1圈是通路編號和分類，具體編號對應(yīng)什么通路名稱可以在代碼輸出的Excel文件中查詢阱冶； -
第2圈表示這個通路有多少個基因刁憋； -
第3圈分為深綠和淺綠，二者加和始終是一樣的木蹬，表示高表達基因的數(shù)目至耻，深綠表示其中有多少基因?qū)儆谶@個通路，淺綠是不屬于這個通路的基因數(shù)目镊叁； -
第4圈是富集分析的顯著性尘颓，-log10(p.adjust)，有一條灰色的豎線表示-log10(0.01)晦譬； -
第5圈是富集因子疤苹，等于差異基因中落到這個通路的基因數(shù)除以這個通路的基因總數(shù)（第三圈深綠除以第二圈紫色）
代碼演示
library(Seurat)
library(tidyverse)
library(xlsx)
testseu=readRDS("testseu.rds")
Idents(testseu)="anno_new"
### 找差異基因 #########################################################################
marker_celltype=FindAllMarkers(testseu,logfc.threshold = 0.8,only.pos = T)
# 過濾
marker_celltype=marker_celltype%>%filter(p_val_adj < 0.01)
marker_celltype$d=marker_celltype$pct.1-marker_celltype$pct.2
marker_celltype=marker_celltype%>%filter(d > 0.2)
marker_celltype=marker_celltype%>%arrange(cluster,desc(avg_log2FC))
marker_celltype=as.data.frame(marker_celltype)
write.xlsx(marker_celltype,file = "markers_log2fc0.8_padj0.01_pctd0.2.xlsx",row.names = F,col.names = T)
### 富集分析 ###########################################################################
library(clusterProfiler)
library(org.Hs.eg.db)
R.utils::setOption("clusterProfiler.download.method","auto") #https://github.com/YuLab-SMU/clusterProfiler/issues/256
source("syEnrich.R")
syEnrich(marker_celltype,outpath = "markers_log2fc0.8_padj0.01_pctd0.2")
### 挑一類細胞來作為演示 #######################################################
kegg.res=read.xlsx("markers_log2fc0.8_padj0.01_pctd0.2.KEGG.xls",sheetIndex = 1,as.data.frame = T,header = T)
kegg.res=kegg.res%>%filter(p.adjust < 0.05)
kegg.res=kegg.res%>%filter(cluster == "Endothelial")
導(dǎo)入《KEGG通路的從屬/注釋信息如何獲取》這一講的文件
# 導(dǎo)入《KEGG通路的從屬/注釋信息如何獲取》這一講的文件
kegg_info=read.xlsx("kegg_info.xlsx",sheetIndex = 1,startRow = 3)
kegg_info=kegg_info[,c("ID","Pathway","big.annotion")]
# 合并兩個表格
kegg.res$ID=str_replace(kegg.res$ID,"hsa","")
kegg.res=kegg.res%>%inner_join(kegg_info,by = "ID")
# 畫圖展示的term控制在20到35個
kegg.res=kegg.res%>%arrange(p.adjust)
kegg.res=head(kegg.res,20)
write.table(kegg.res,file = "kegg.res.txt",quote = F,sep = "\t",row.names = F,col.names = T)
write.xlsx(kegg.res,file = "kegg.res.xlsx",col.names = T,row.names = F)
調(diào)用畫圖函數(shù)
### 調(diào)用畫圖函數(shù) ###############################################################
source("kegg_loop.R")
kegg_loop(enrich.res = kegg.res,filename = "Endothelial_kegg")
然后就能得到推文開頭那張圖了。
獲取代碼
后-臺回復(fù)2022A可知敛腌。
說點別的
這個系列到今天一共更新了10篇帖子(9篇2022A卧土，1篇2022B)，前后4個月的跨度像樊，算是更新很慢的了尤莺。一方面，學(xué)校的事情越來越多生棍，寫博客的時間越來越少颤霎；另一方面，這些原創(chuàng)度很高的繪圖帖子需要我花很多精力整理總結(jié)足绅。
