這一節(jié)畫的圖是比較新的，圖中我用紅色箭頭標(biāo)出的是pathway
一級注釋信息(big annotation，自己想的，非專有名詞)，縱軸花花綠綠的標(biāo)注是pathway的二級注釋(small annotation)。如何獲取注釋算一個(gè)難點(diǎn)，我上一講也已經(jīng)講過：KEGG通路的從屬/注釋信息如何獲取。
整個(gè)圖反映的是有多少基因落到了對應(yīng)的分類里面。
辯證地看，整張圖都是pathway注釋，沒有具體的pathway名稱，跟平常做的富集分析很不一樣。把圖里面的二級注釋換成具體的pathway會(huì)更好。另外，這個(gè)圖中橫坐標(biāo)是基因數(shù)，但我覺得在富集分析中這個(gè)數(shù)值并不重要（基因ratio比單個(gè)數(shù)值重要），我們可以換成p值。
在開始畫圖之前，需要整理一個(gè)表格，表格中至少包含：pathway ID、pathway description、pathway注釋/big annotation、幾個(gè)富集指標(biāo)、被富集到的基因。
- 前三個(gè)信息，上一講的表格已經(jīng)整理好了，可以拿來用；
-
幾個(gè)富集指標(biāo)在
clusterprofiler的輸出結(jié)果中也有； - 被富集到的基因貌似clusterprofiler做kegg富集不能直接顯示基因symbol，所以返回的結(jié)果需要稍微加工一下。
整理表格的代碼名稱為run.R，如下所示：
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")
### 挑一類細(xì)胞來作為演示 #######################################################
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_info=read.xlsx("kegg_info.xlsx",sheetIndex = 1,startRow = 3)
kegg_info=kegg_info[,c("ID","Pathway","big.annotion")]
# 合并兩個(gè)表格
kegg.res$ID=str_replace(kegg.res$ID,"hsa","")
kegg.res=kegg.res%>%inner_join(kegg_info,by = "ID")
write.table(kegg.res,file = "kegg.res.txt",quote = F,sep = "\t",row.names = F,col.names = T)
我以單細(xì)胞分析中的kegg富集分析作為演示，只取其中一個(gè)cluster的富集結(jié)果來畫圖。
上述代碼中間用到的富集代碼叫syEnrich.R，這個(gè)文件只需要輸入單細(xì)胞seurat對象運(yùn)行FindAllMarkers得到的差異基因，就可以返回GO/KEGG富集結(jié)果，同時(shí)被富集到某個(gè)通路的基因symbol也會(huì)被列出。
運(yùn)行run.R之后，最終的表格如下圖所示：
然后開始畫圖，代碼名稱為3.plot.R，這里就不演示了，最終可以得到的圖如下：
獲取代碼
包含這張圖會(huì)用到的所有代碼，數(shù)據(jù)整理以及畫圖，超貼心有沒有！
這個(gè)系列都會(huì)采取限時(shí)公開的方式共享代碼，24小時(shí)內(nèi)是免費(fèi)的。超過這個(gè)時(shí)間如何獲取，公粽號發(fā)送2022A可知。
鏈接：https://pan.baidu.com/s/1hebbeQH4DgYA8JqFWXO4dg
提取碼：abo5
