使用的是Y叔的包
主要是構(gòu)建geneList數(shù)據(jù)結(jié)構(gòu)
不使用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)
#需要將差異倍數(shù)logFC按從高到底排序脾歇,同時將gene name轉(zhuǎn)化為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)