R语言Pathway可视化

使用的是Y叔的包
主要是构建geneList数据结构
不使用Y叔的包参照这个

#pathway analysis function
go_analysis<-function(gene_symbols){
  options(connectionObserver = NULL)
  suppressMessages(library(clusterProfiler))
  suppressMessages(library(enrichplot))
  suppressMessages(library(DOSE))
  suppressMessages(library(org.Hs.eg.db))
  gene_symbols<-gene_symbols[!duplicated(gene_symbols)]
  gene_list<-AnnotationDbi::select(org.Hs.eg.db, keys=as.character(gene_symbols), columns=c("SYMBOL","ENTREZID"), keytype="SYMBOL")
  gene_id<-as.character(gene_list$ENTREZID)
  ego <- enrichGO(gene_id, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)
  return(ego)
}  

#栗子一
library(clusterProfiler)
library(enrichplot)
#需要将差异倍数logFC按从高到底排序，同时将gene name转化为NCBI的ID
data<-read.csv("~/Desktop/DEGs.csv",header = T)
geneList<-data$logFC
names(geneList)<-data$geneID
de<-as.character(data$geneID)
ego <- enrichGO(de, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)
goplot(ego)
barplot(ego, showCategory=20)
dotplot(ego, showCategory=30)
ego2 <- simplify(ego)
cnetplot(ego2, foldChange=geneList)
cnetplot(ego2, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
heatplot(ego2, foldChange=geneList)
upsetplot(ego)
emapplot(ego2)
kk <- gseKEGG(geneList, nPerm=1000)
ridgeplot(kk)
#栗子二 有点繁琐
library(tidyverse)
library(org.Mm.eg.db)
library(clusterProfiler)
df<-read.csv("～/diff_expr_result.csv")
df_dig<-df %>% filter(logFC>1&padj<0.01) 

deg<-bitr(df_dig$X,fromType = "SYMBOL",toType = "ENTREZID",OrgDb = org.Mm.eg.db)%>%
  left_join(df_dig,by=c("SYMBOL"="X"))%>%
  distinct(ENTREZID,.keep_all = TRUE)

genelist<-deg$logFC
names(genelist)<-deg$ENTREZID
genelist<-sort(genelist,decreasing = T)
ego<-enrichGO(
  gene=names(genelist),#should use this 
  OrgDb = org.Mm.eg.db,
  readable = T,
  ont = "BP",#MF,CC
  pvalueCutoff = 0.05,
  qvalueCutoff = 0.05)

cnetplot(ego,
         #showCategory = 5,
         foldChange = genelist,
         circular=TRUE,
         colorEdge=TRUE)

我的一个例子

rt<-human_dif
filename<-"human_dif"

human_fasting_fed_dif<-human_dif[order(human_dif$log2FoldChange,decreasing = T),]
new_names<-unlist(lapply(row.names(rt), FUN = function(x) {return(strsplit(x, split = ".", fixed=T)[[1]][1])}))
row.names(rt)<-new_names
gene_list<-select(org.Hs.eg.db, keys=as.character(new_names), columns=c("SYMBOL","ENTREZID"), keytype="ENSEMBL") 
gene_list<-gene_list[!duplicated(gene_list$ENSEMBL),]
row.names(gene_list)<-as.character(gene_list$ENSEMBL)
a<-intersect(row.names(rt),row.names(gene_list))
data<-cbind(rt[a,],gene_list[a,])
data<-na.omit(data)
geneList<-data[,8]
names(geneList)<-as.character(data$ENTREZID)
de<-as.character(data$ENTREZID)
ego <- enrichGO(de, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)

goplot<-goplot(ego)
barplot<-barplot(ego, showCategory=20)
dotplot<-dotplot(ego, showCategory=30)
ego2 <- simplify(ego)
cnetplot<-cnetplot(ego2, foldChange=geneList)
cnetplot2<-cnetplot(ego2, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
heatplot<-heatplot(ego2, foldChange=geneList)
upsetplot<-upsetplot(ego)
emapplot<-emapplot(ego2)
ggsave(plot=goplot,paste0(filename,"_goplot.pdf"),device = "pdf")
ggsave(plot=barplot,paste0(filename,"_barplot.pdf"),device = "pdf")
ggsave(plot=dotplot,paste0(filename,"_dotplot.pdf"),device = "pdf")
ggsave(plot=cnetplot,paste0(filename,"_cnetplot.pdf"),device = "pdf")
ggsave(plot=cnetplot2,paste0(filename,"_cnetplot2.pdf"),device = "pdf")
ggsave(plot=heatplot,paste0(filename,"_heatplot.pdf"),device = "pdf")
ggsave(plot=upsetplot,paste0(filename,"_upsetplot.pdf"),device = "pdf")
ggsave(plot=upsetplot,paste0(filename,"_emapplot.pdf"),device = "pdf")
kk <- gseKEGG(geneList, nPerm=1000)
ridgeplot<-ridgeplot(kk)
ggsave(plot=ridgeplot,paste0(filename,"_ridgeplot"),device = "pdf")

Group之间比较
参照clusterProfiler-book第十一章

#use ENTREZID to generate the pathway_data
pathway_data<-data.frame(c(as.character(entrezid1),as.character(entrezid2)))
colnames(pathway_data)<-"ENTREZID"
pathway_data$"Group"<-c(rep("Group1",nrow(human_fasting)),
                              rep("Group2",nrow(mouse_fasting)),
                              rep("Group1",nrow(human_refed)),
                              rep("Group2",nrow(mouse_refed)))

pathway_data$"Treatment"<-c(rep("Treatment1",nrow(human_fasting)),
                            rep("Treatment1",nrow(mouse_fasting)),
                            rep("Treatment2",nrow(human_refed)),
                            rep("Treatment2",nrow(mouse_refed)))
kegg_pathway <- compareCluster(ENTREZID~Species+Treatment, data=pathway_data,
                                       fun='enrichKEGG')

all_dot_plot<-dotplot(kegg_pathway, x=~Group,showCategory=20,color = "p.adjust") + ggplot2::facet_grid(~Treatment)

ggsave(filename = "/out_dir/Combind_KEGG_plot.pdf",height = 10,width = 10,units = "in")

例子3

bp_analysis<-function(x,FC,p){
  suppressMessages(library(clusterProfiler))
  suppressMessages(library(enrichplot))
  suppressMessages(library(org.Hs.eg.db))
  suppressMessages(library(tidyverse))
  new_out_dir<-paste0(out_dir,x)
  dir.create(new_out_dir,recursive = T)
  degs<-get(x)
  degs_sig<-degs %>% filter(pvalue<p & abs(log2FoldChange)>FC &gene_type=="protein_coding")
  degs_sig<-degs_sig[!duplicated(degs_sig$gene_name),]
  degs_sig_pick<-degs_sig[,c("gene_name","log2FoldChange")]
  colnames(degs_sig_pick)<-c("SYMBOL","log2FoldChange")
  row.names(degs_sig)<-as.character(degs_sig$gene_name)
  gene_list<-AnnotationDbi::select(org.Hs.eg.db, keys=as.character(row.names(degs_sig)), columns=c("SYMBOL","ENTREZID"), keytype="SYMBOL")
  rt<-degs_sig_pick %>% left_join(gene_list,by="SYMBOL")
  rt<-rt[order(rt$log2FoldChange,decreasing = T),]
  geneList<-rt$log2FoldChange
  names(geneList)<-as.character(rt$ENTREZID)
  de<-as.character(rt$ENTREZID)
  ego <- enrichGO(de, OrgDb = "org.Hs.eg.db", ont="BP", readable=TRUE)
  ego2 <- simplify(ego)
  cnetplot<-cnetplot(ego2, foldChange=geneList, circular = TRUE, colorEdge = TRUE)
  return(cnetplot)
}
h6<-read.csv(paste0(data_dir,"h6_v_v5_GenePass_DESeq2.csv"))
h6_bp<-bp_analysis(x="h6",FC=1.5,p=0.05)