单细胞数据使用ClusterGVis来可视化不同细胞类型的marker基因的表达情况
只需要提供Seurat对象和Seurat::FindAllMarkers的计算结果即可：
devtools::install_github("junjunlab/ClusterGVis")
library(ClusterGVis)
library(Seurat)
load("./data_harmony_final.RData")
Idents(data_harmony) <- data_harmony$celltype
pbmc <- data_harmony
pbmc.markers.all <- Seurat::FindAllMarkers(pbmc,only.pos = TRUE,min.pct=0.25,logfc.threshold=0.25)
# get top 20 genes
pbmc.markers <- pbmc.markers.all
%>% dplyr::group_by(cluster)
%>% dplyr::top_n(n = 20, wt = avg_log2FC)
# check
head(pbmc.markers)
# prepare data from seurat
#以细胞群的均值作为输入,且默认对基因进行去重
objectst.data <- prepareDataFromscRNA(
object = pbmc,
diffData = pbmc.markers,
showAverage = TRUE)
# check
str(objectst.data)
# [2] 可选，不对基因进行去重处理，添加后缀
# objectst.data2 <- prepareDataFromscRNA(object = scRNA_harmony,
# diffData = pbmc.markers,
# showAverage = TRUE,
# keep.uniqGene = FALSE,
# sep = "_") #"gene_1","gene_2",....
# [3] 可选，以每个细胞单独的表达作为输入
# objectst.data3 <- prepareDataFromscRNA(object = scRNA_harmony,
# diffData = pbmc.markers,
# showAverage = FALSE)
富集分析
library(org.Hs.eg.db)
# enrich for clusters
enrich.go <- enrichCluster(object = objectst.data,
OrgDb = org.Hs.eg.db,
type = "BP", # c("BP", "MF", "CC", "KEGG", "ownSet")
organism = "hsa",
pvalueCutoff = 0.5,
topn = 5,seed = 5201314)
# check
head(enrich)
# add gene name
markGenes = unique(pbmc.markers$gene)[sample(1:length(unique(pbmc.markers$gene)),40, replace = F)]
# line plot
visCluster(object = objectst.data,plot.type = "line")
# heatmap
pdf('sc1.pdf',height = 10,width = 6,onefile = F)
visCluster(object = objectst.data,
plot.type = "heatmap",
column_names_rot = 35,
markGenes = markGenes,cluster.order = c(1:7))
dev.off()
# add bar
pdf('sc2.pdf',height = 10,width = 14,onefile = F)
visCluster(object = objectst.data,
plot.type = "both",
column_names_rot = 35,
show_row_dend = F,
markGenes = markGenes,
markGenes.side = "left",
annoTerm.data = enrich.go,
line.side = "left",
cluster.order = c(1:7), go.col = rep( ggsci::pal_d3()(7), each = 5),
textbar.pos = c(0.8,0.2),add.bar = T)
dev.off()
enrich.go.BP
同时执行 GO 和 KEGG
enrich.kegg <- enrichCluster(object = objectst.data,
OrgDb = org.Hs.eg.db,
type = "KEGG",
organism = "hsa",
pvalueCutoff = 0.05,
topn = 5,
seed = 5201314)
head(enrich.kegg)
## group Description pvalue ratio
## hsa05340...1 C1 Primary immunodeficiency 1.909178e-04 40.0
## hsa04060...2 C1 Cytokine-cytokine receptor interaction 3.984338e-04 60.0
## hsa04934 C1 Cushing syndrome 3.153205e-03 40.0
## hsa04612 C2 Antigen processing and presentation 4.189612e-14 87.5
## hsa04650 C2 Natural killer cell mediated cytotoxicity 1.572798e-12 87.5
## hsa05332 C2 Graft-versus-host disease 4.624360e-08 50.0
dim(enrich.kegg) ; nrow(enrich.kegg) # 这里我的kegg，有的组富集不到n=5
visCluster(object = objectst.data,
plot.type = "both",
column_names_rot = 35,
show_row_dend = F,
markGenes = markGenes,
markGenes.side = "left",
annoTerm.data = enrich.go,
go.col = rep( ggsci::pal_d3()(7), each = 5),
annoKegg.data = enrich.kegg,
kegg.col = rep(jjAnno::useMyCol("stallion",n = 7), each = 5)[1:nrow(enrich.kegg)], # 这里我的kegg，有的组富集不到n=5
line.side = "left",
cluster.order = 1:7
#sample.group = rep(c("sample1","sample2","sample3"),each = 2) # 这里我没有分组
)
# 格式化注释面板
visCluster(object = objectst.data,
plot.type = "both",
column_names_rot = 35,
show_row_dend = F,
markGenes = markGenes,
markGenes.side = "left",
annoTerm.data = enrich.go,
go.col = rep(jjAnno::useMyCol("calm",n = 7),each = 5),
by.go = "anno_block",
annoKegg.data = enrich.kegg,
kegg.col = rep(jjAnno::useMyCol("stallion",n = 7),each = 5)[1:nrow(enrich.kegg)], # 这里我的kegg，有的组富集不到n=5
by.kegg = "anno_block",
word_wrap = F,add_new_line = F,
line.side = "left",
cluster.order = 1:7
#sample.group = rep(c("sample1","sample2","sample3"),each = 2)
)
# 添加 GO 和 KEGG 的条形图
pdf('ClusterGVis_GO_KEGG.pdf',height = 16, width = 22,onefile = F)
visCluster(object = objectst.data,
plot.type = "both",
column_names_rot = 35,
show_row_dend = F,
markGenes = markGenes,
markGenes.side = "left",
annoTerm.data = enrich.go,
go.col = rep(ggsci::pal_d3()(7), each = 5),
add.bar = T,
by.go = "anno_block",
annoKegg.data = enrich.kegg,
kegg.col = rep(jjAnno::useMyCol("calm",n = 7),each = 5 )[1:nrow(enrich.kegg)], # # 这里我的kegg，有的组富集不到n=5
add.kegg.bar = T,
by.kegg = "anno_block",
word_wrap = F, add_new_line = F,
line.side = "left",
cluster.order = 1:7
# sample.group = rep(c("sample1","sample2","sample3"),each = 2)
)
dev.off()
GO 和 KEGG
GO 和 KEGG 结果合并在一起
# 将 GO 和 KEGG 结果合并并绘制在一个面板中，使用不同的颜色区分
library(dplyr)
all <- rbind(enrich.go %>% group_by(group) %>% slice_head(n = 5) %>%
mutate(Description = paste("BP: ",Description,sep = "")),
enrich.kegg %>% group_by(group) %>% slice_head(n = 5) %>%
mutate(Description = paste("KEGG: ",Description,sep = ""))) %>%
arrange(group)
dim(all);table(all$group)
# plot
pdf('ClusterGVis_rbind.GO.KEGG.pdf',height = 16,width = 22,onefile = F)
visCluster(object = objectst.data,
plot.type = "both",
column_names_rot = 35,
show_row_dend = F,
markGenes = markGenes,
markGenes.side = "left",
annoTerm.data = all,
go.col = c( rep(ggsci::pal_d3()(7), each = 5),
rep(jjAnno::useMyCol("calm",n = 7),each = 5 )[1:nrow(enrich.kegg)] ),
line.side = "left",
word_wrap = F,add_new_line = F,
by.go = "anno_block",
add.bar = T,
cluster.order = 1:7,
#sample.group = rep(c("sample1","sample2","sample3"),each = 2)
)
dev.off()
#保存
save(enrich.go,enrich.kegg,all,objectst.data, file="ClusterGVis.GO.KEGG.RData")
GO 和 KEGG 合并后的
