marker list的富集分析

前面獲得的來自于不同celltype的marker,或者不同cluster的marker豪娜,或者每個celltype中不同分組的差異基因餐胀,上述marker構(gòu)成的列表形式,可以直接輸入瘤载,獲得富集結(jié)果的list以及圖片否灾。


1. 加載包

library(Seurat)
library(ggplot2)
library(openxlsx)
library(homologene)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(stringr)

2. 富集分析,org指定物種"mmu"或者"hsa",markers_list為單細(xì)胞數(shù)據(jù)中獲得的maker的list

EnrichmentForMarkerlist <- function(markers_list,org,item.size=8,tag,path.out){
  if (org=="mmu") {
  db <- "org.Mm.eg.db"
}else if(org=="hsa"){
  db <- "org.Hs.eg.db"
}
GOenrich_result_list <- list()
KEGGenrich_result_list <- list()
for (i in names(markers_list)){
    print(paste("cluster",i,"enrichment",sep=" "))
    if (nrow(markers_list[[i]]) == 0){
        print(paste("nrow cluster",i,"is",nrow(markers_list[[i]]),sep=" "))
        next
    }
  #GO富集分析
  print(paste("cluster",i,"GO entichment",sep=" "))
    GOenrich_result <- enrichGO(rownames(markers_list[[i]]),db,keyType="SYMBOL", ont="all", pvalueCutoff=0.05)
#畫三種分類的圖
print(paste("cluster",i,"GO all entichment",sep=" "))
if (nrow(as.data.frame(GOenrich_result)) > 0){
    GO_bar_split <- barplot(GOenrich_result, split = "ONTOLOGY",showCategory=20) + #GOenrich_result不能是矩陣
    facet_grid(ONTOLOGY~., scale = "free") + 
    scale_y_discrete(labels=function(x) str_wrap(x, width=100)) +  #取消了iterm強(qiáng)制換行鸣奔,需要加載stringr包
    theme(axis.text.y=element_text(size=7)) 
    try(ggsave(paste0("cluster_",i,"_GO_bar_all_plot.png"),GO_bar_split,path=path.out,width=15,height=10))
#只畫BP的圖
 print(paste("cluster",i,"GO BP entichment",sep=" "))
    GOenrich_result_BP <- enrichGO(rownames(markers_list[[i]]),db,keyType="SYMBOL", ont="BP", pvalueCutoff=0.05)
    GO_bar_BP <- barplot(GOenrich_result_BP,showCategory=30,drop=T) + 
    theme(axis.text.y=element_text(size=item.size))  + 
    scale_y_discrete(labels=function(x) str_wrap(x, width=100)) #取消了iterm強(qiáng)制換行墨技,需要加載stringr包
    try(ggsave(paste0("cluster_",i,"_GO_bar_BP_plot.png"),GO_bar_BP,path=path.out,limitsize = FALSE))
}
#KEGG富集分析
    ID_transform_table <- bitr(rownames(markers_list[[i]]),fromType="SYMBOL",toType="ENTREZID",OrgDb=db)
    ID_Entrez_gene_set <- ID_transform_table$ENTREZID
    print("KEGGenrich start")
    KEGGenrich_result <- enrichKEGG(ID_Entrez_gene_set, organism = org, keyType = 'kegg', pvalueCutoff = 0.05,pAdjustMethod = 'BH', minGSSize = 10,maxGSSize = 500,use_internal_data = FALSE)
#繪制氣泡圖
print("KEGGenrich start map")
    try(KEGG_dot <- dotplot(KEGGenrich_result, showCategory=20) + theme(axis.text.y=element_text(size=item.size)) +
        scale_y_discrete(labels=function(x) str_wrap(x, width=100)))
    print("map saved")
    try(ggsave(paste0("cluster_",i,"_KEGG_dot_plot.png"),KEGG_dot,path=path.out,width=10,height=6))

try(GOenrich_result_list[[i]] <- GOenrich_result_BP)
try(KEGGenrich_result_list[[i]] <- KEGGenrich_result)
print("result saved")
}
print("start output")
  write.xlsx(GOenrich_result_list,file = paste0(path.out,"/",tag,'_GOenrich_result_list.xlsx'),row.names=T,sep='\t',overwrite=T)
  write.xlsx(KEGGenrich_result_list,file = paste0(path.out,"/",tag,'_KEGGenrich_result_list.xlsx'),row.names=T,sep='\t',overwrite=T)
  saveRDS(GOenrich_result_list,file=paste0(path.out,"/",tag,'_GO_EnrichResult.rds'))
  saveRDS(KEGGenrich_result_list,file=paste0(path.out,"/",tag,'_KEGG_EnrichResult.rds'))
}

3.加載數(shù)據(jù)

EnrichmentForMarkerlist(markers_list=all_cluster_markers,org="hsa",item.size=8,tag="clusters",path.out=path_out)

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 序言:七十年代末,一起剝皮案震驚了整個濱河市挎狸,隨后出現(xiàn)的幾起案子扣汪,更是在濱河造成了極大的恐慌,老刑警劉巖伟叛,帶你破解...
    沈念sama閱讀 217,185評論 6 503
  • 序言:濱河連續(xù)發(fā)生了三起死亡事件私痹,死亡現(xiàn)場離奇詭異,居然都是意外死亡统刮,警方通過查閱死者的電腦和手機(jī)紊遵,發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 92,652評論 3 393
  • 文/潘曉璐 我一進(jìn)店門,熙熙樓的掌柜王于貴愁眉苦臉地迎上來侥蒙,“玉大人暗膜,你說我怎么就攤上這事”揆茫” “怎么了学搜?”我有些...
    開封第一講書人閱讀 163,524評論 0 353
  • 文/不壞的土叔 我叫張陵,是天一觀的道長论衍。 經(jīng)常有香客問我瑞佩,道長,這世上最難降的妖魔是什么坯台? 我笑而不...
    開封第一講書人閱讀 58,339評論 1 293
  • 正文 為了忘掉前任炬丸,我火速辦了婚禮,結(jié)果婚禮上蜒蕾,老公的妹妹穿的比我還像新娘稠炬。我一直安慰自己焕阿,他們只是感情好,可當(dāng)我...
    茶點(diǎn)故事閱讀 67,387評論 6 391
  • 文/花漫 我一把揭開白布首启。 她就那樣靜靜地躺著暮屡,像睡著了一般。 火紅的嫁衣襯著肌膚如雪毅桃。 梳的紋絲不亂的頭發(fā)上褒纲,一...
    開封第一講書人閱讀 51,287評論 1 301
  • 那天,我揣著相機(jī)與錄音疾嗅,去河邊找鬼外厂。 笑死,一個胖子當(dāng)著我的面吹牛代承,可吹牛的內(nèi)容都是我干的汁蝶。 我是一名探鬼主播,決...
    沈念sama閱讀 40,130評論 3 418
  • 文/蒼蘭香墨 我猛地睜開眼论悴,長吁一口氣:“原來是場噩夢啊……” “哼掖棉!你這毒婦竟也來了?” 一聲冷哼從身側(cè)響起膀估,我...
    開封第一講書人閱讀 38,985評論 0 275
  • 序言:老撾萬榮一對情侶失蹤幔亥,失蹤者是張志新(化名)和其女友劉穎,沒想到半個月后察纯,有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體帕棉,經(jīng)...
    沈念sama閱讀 45,420評論 1 313
  • 正文 獨(dú)居荒郊野嶺守林人離奇死亡,尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點(diǎn)故事閱讀 37,617評論 3 334
  • 正文 我和宋清朗相戀三年饼记,在試婚紗的時候發(fā)現(xiàn)自己被綠了香伴。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點(diǎn)故事閱讀 39,779評論 1 348
  • 序言:一個原本活蹦亂跳的男人離奇死亡具则,死狀恐怖即纲,靈堂內(nèi)的尸體忽然破棺而出,到底是詐尸還是另有隱情博肋,我是刑警寧澤低斋,帶...
    沈念sama閱讀 35,477評論 5 345
  • 正文 年R本政府宣布,位于F島的核電站匪凡,受9級特大地震影響膊畴,放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜病游,卻給世界環(huán)境...
    茶點(diǎn)故事閱讀 41,088評論 3 328
  • 文/蒙蒙 一巴比、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧,春花似錦轻绞、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 31,716評論 0 22
  • 文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至兼砖,卻和暖如春奸远,著一層夾襖步出監(jiān)牢的瞬間,已是汗流浹背讽挟。 一陣腳步聲響...
    開封第一講書人閱讀 32,857評論 1 269
  • 我被黑心中介騙來泰國打工懒叛, 沒想到剛下飛機(jī)就差點(diǎn)兒被人妖公主榨干…… 1. 我叫王不留,地道東北人耽梅。 一個月前我還...
    沈念sama閱讀 47,876評論 2 370
  • 正文 我出身青樓薛窥,卻偏偏與公主長得像,于是被迫代替她去往敵國和親眼姐。 傳聞我的和親對象是個殘疾皇子诅迷,可洞房花燭夜當(dāng)晚...
    茶點(diǎn)故事閱讀 44,700評論 2 354

推薦閱讀更多精彩內(nèi)容