?詳情請聯(lián)系作者：
? Monocle2擬時分析及可視化我們整了好多呀，算是整明白了吧,。上一期我們做了Monocle3的擬時基因富集分析（Monocle3個性化分析作圖：擬時熱圖/擬時基因GO分析/擬時基因趨勢分析）,。Monocle2的還沒有做過，這里我們演示下,，其實(shí)方法是一樣的,。monocle2擬時熱圖修飾我們已經(jīng)做過了。這里的演示我們不再使用修飾,，用monocle2本來的擬時熱圖繪制函數(shù)進(jìn)行,，然后通過聚類不同的module提取基因進(jìn)行富集分析。本貼涉及到幾個修改函數(shù)和自寫函數(shù)的加載,，數(shù)據(jù)代碼已上傳QQ群,！
接下來我們就具體看看做法，首先還是做熱圖,。這里我們用R包里面的plot_pseudotime_heatmap函數(shù)即可,。#以下是monocle2做熱圖的代碼，很常規(guī)的操作library(monocle)library(Seurat)library(dplyr)library(viridis)library(pheatmap)library(grid)#======================================================================================
setwd("D:/KS項(xiàng)目/公眾號文章/monocle2擬時熱圖富集分析")mouse_monocle <- readRDS("D:/KS項(xiàng)目/公眾號文章/monocle2擬時結(jié)果個性化作圖/mouse_monocle.rds")peu_gene <- differentialGeneTest(mouse_monocle,fullModelFormulaStr = "~sm.ns(Pseudotime)",cores = 2)write.csv(peu_gene,file='peu_gene.csv')#保存好文件,，這個分析過程挺費(fèi)時間
peu_gene <- read.csv('peu_gene.csv', header = T, row.names = 1)peu_gene <- peu_gene[which(peu_gene$qval<0.01 & peu_gene$num_cells_expressed>100),]#篩選顯著是的peu_gene %>% arrange(qval) -> peu_gene#按照qval排個序peu_gene <- peu_gene[1:100,] #這里我們?nèi)∏?00個基因演示
#熱圖的做法就很多了,，可以用它默認(rèn)的函數(shù)，也可以用我們講過的ggplot或者pheatmap自己修飾p <-plot_pseudotime_heatmap(mouse_monocle[peu_gene$gene_short_name,], num_clusters = 4, cores = 2, show_rownames = T,return_heatmap =T, hmcols = viridis(256), use_gene_short_name = T)
p
class(p)# [1] "pheatmap"
 plot_pseudotime_heatmap函數(shù)好像并不能選擇性展示感興趣的基因,，然而提取數(shù)據(jù)用pheatmap或者ggplot作圖我們也做過了,，這里再去修飾一遍展示沒有必要，我就想用monocle2包熱圖函數(shù)完成。class一下plot_pseudotime_heatmap函數(shù)的產(chǎn)生的圖,，發(fā)現(xiàn)是pheatmap對象,，查看作者的原函數(shù)也發(fā)現(xiàn)這個熱圖作者是基于pheatmap做的，那就好辦了,，如果連續(xù)閱讀我們帖子的小伙伴可能會想到,，我們之前寫過一篇帖子：學(xué)習(xí)nature medicine函數(shù)：熱圖標(biāo)記任意基因，這里面我們引用了一個文章提供的函數(shù),，可以很方便的為pheatmap標(biāo)記基因,。。然而呢,，由于上圖是分裂的,，所以這個函數(shù)效果不好，那么辦法要么就是修改這個添加基因的函數(shù),，另一個方法就是修改plot_pseudotime_heatmap函數(shù),。權(quán)衡之下，偷懶的我選擇了修改plot_pseudotime_heatmap函數(shù),，因?yàn)樗眯薷模?/span>#改造熱圖函數(shù),，展示需要的基因source('./new_heatmap.R')p1 <-plot_pseudotime_heatmap(mouse_monocle[peu_gene$gene_short_name,], num_clusters = 4, cores = 2, show_rownames = T,return_heatmap =T, hmcols = viridis(256), use_gene_short_name = T)
p1
這樣就既有基因分組注釋，但是熱圖行不會分裂開,。下面使用add.flag.R函數(shù)標(biāo)記感興趣基因,。#只展示感興趣的基因genes <- c("Cxcr2","Celf2","Gsr","Sorl1","Ftl1","Ptprc", "Cst3","Il1b","Fgl2","Sirpb1c","Maml3","Peak1", "Orm1","H2afz","Cpa6","Ica1","Fmnl2","Cpne3", "Anxa3","Tbc1d8","Camp","Mmp8","S100a8","S100a9")
source('./add.flag.R')add.flag(p,kept.labels = genes,repel.degree = 0.2)
接下來就是提取module基因進(jìn)行富集分析了，我們直接包裝為一個函數(shù)Monocle2_gene_enrichment.R,，直接運(yùn)行即可,，當(dāng)然有個限制我們只寫了GO富集分析，感興趣可以自己添加KEGG或者其他的分析,。
#提取module基因,，進(jìn)行富集分析
library(clusterProfiler)library(ggplot2)source('./Monocle2_gene_enrichment.R')GOannlysis <- Monocle2_gene_enrichment(p,knum=4,species='org.Mm.eg.db', pvalueCutoff=0.05, qvalueCutoff=0.05)
table(GOannlysis$Cluster)#module1基因最少，沒有顯著性富集# 1 2 3 4 # 0 76 43 77 write.csv(GOannlysis, file = 'GOannlysis.csv')
最后,，將熱圖在AI中修飾,，旁邊展示富集通路，這個圖就完整了,。好了,，這就是monocle2的全部內(nèi)容了，覺得分享對您有幫助的,，點(diǎn)個贊分享下再走唄,！
|
