clusterProfiler进行富集分析并使用ggplot2进
2022-11-30 本文已影响0人
yingyonghui
library(clusterProfiler)
library(org.Mm.eg.db)
###富集分析
sig.entrezid.id <- bitr(sig.genes, fromType="SYMBOL", toType="ENTREZID", OrgDb="org.Mm.eg.db")[, 2]
GO.res <- enrichGO(gene=sig.entrezid.id, OrgDb="org.Mm.eg.db", ont="ALL", readable=T)
GO.res <- as.data.frame(GO.res)
catch.count <- as.numeric(sapply(GO.res$GeneRatio, function(x){strsplit(x, split='/')[[1]][1]}))
gene.count <- as.numeric(sapply(GO.res$GeneRatio, function(x){strsplit(x, split='/')[[1]][2]}))
GO.res$GeneRatio <- catch.count/gene.count
write.csv(GO.res,paste0('celltype.marker.GO.FC1.P05.csv'))
###可视化
GO.res.dat <- GO.res.dat[order(GO.res.dat$p.adjust, decreasing=F), ]
# to select GO terms based on ONTOLOGY (CC, BP, and MF)
# GO.term <- unlist(by(data=rownames(GO.res.dat), INDICES=GO.res.dat$ONTOLOGY, function(x){x[1:10]}))
# names(GO.term) <- NULL
# GO.res.dat <- GO.res.dat[GO.term, ]
GO.res.dat <- GO.res.dat[1:20, ]
all(GO.res.dat$p.adjust < 0.05)
GO.res.dat <- GO.res.dat[order(GO.res.dat$GeneRatio, decreasing=T), ]
GO.res.dat$Description <- factor(GO.res.dat$Description, levels=rev(c(GO.res.dat$Description)))
pdf(paste0('celltype.marker.GO.FC1.P05.pdf'),height=6,width=9.5)
ggplot(GO.res.dat, aes(GeneRatio, Description)) +
geom_point(aes(size=Count, color=-log10(p.adjust))) +
scale_colour_gradient(low="green",high="red") +
# scale_x_continuous(limits=c(1,4),breaks=c(1,2,3,4)) +
# scale_size_continuous(breaks=c(5,10,15), limits=c(4,15)) +
# facet_grid(ONTOLOGY~., scales="free") +
labs(x='Gene ratio',y='') +
theme(legend.key=element_blank(),
axis.text=element_text(colour='black', size=13),
axis.title=element_text(colour='black', size=13),
panel.background=element_rect(fill=NA, color="black"),
panel.grid=element_line(colour='gray')
)
dev.off()