PlantAI-dev/codeMiaoshu/enrichment/GOendKEGG.md

#==============================================

五十四、基因功能注释

#============================================== library(dplyr) library(tidyr) library(DOSE) library(GO.db) library(org.Hs.eg.db) library(GSEABase) library(clusterProfiler)

dta <- read.csv("res0.05.csv",header = T,row.names = 1,stringsAsFactors = F)

------------------------------------------------------------------------

x <- rownames(dta) length(x) keytypes(org.Hs.eg.db) eg <- bitr(x, fromType="ENSEMBL", toType="ENTREZID", OrgDb="org.Hs.eg.db") head(eg) ids <- bitr(x, fromType="ENSEMBL", toType=c( "SYMBOL"), OrgDb="org.Hs.eg.db") head(ids) go <- bitr(x,fromType = "ENSEMBL",toType = c("SYMBOL","GO","ONTOLOGY"),OrgDb = "org.Hs.eg.db") head(go) #====================================

GO功能富集分析，输入Gene ID 为GI号

#==================================== gene <- eg$ENTREZID gene.df <- bitr(gene, fromType = "ENTREZID", toType = c("ENSEMBL", "SYMBOL"), OrgDb = org.Hs.eg.db) head(gene.df) ggo <- groupGO(gene = gene, OrgDb = org.Hs.eg.db, ont = "MF", level = 3, readable = TRUE) head(ggo) View(as.data.frame(ggo))

ego <- enrichGO(gene = gene, OrgDb = org.Hs.eg.db, ont = "CC", pAdjustMethod = "BH", readable = TRUE) head(ego) View(as.data.frame(ego))

GO功能富集可视化

barplot(ggo, drop=TRUE, showCategory=12)

----fig.height=5, fig.width=8-------------------------------------------

barplot(ego, showCategory=15)

dotplot(ego)

goplot(ego)

#====================================

KEGG富集分析

#==================================== search_kegg_organism('ece', by='kegg_code') ecoli <- search_kegg_organism('Escherichia coli', by='scientific_name') dim(ecoli) head(ecoli)

kk <- enrichKEGG(gene = gene, organism = 'hsa', pvalueCutoff = 0.05) head(kk)

----eval = FALSE--------------------------------------------------------

mkk <- enrichMKEGG(gene = gene,organism = 'hsa')

browseKEGG(kk, 'hsa04110') barplot(kk) dotplot(kk)

#====================================

#单独绘制GO条目图

#==================================== library(ggplot2) go <- read.csv("go.csv",header = T) head(go) go_sort <- go[order(go$Ontology,-go$Percentage),] m <- go_sort[go_sort$Ontology=="Molecular function",][1:10,] c <- go_sort[go_sort$Ontology=="Cellular component",][1:10,] b <- go_sort[go_sort$Ontology=="Biological process",][1:10,] slimgo <- rbind(b,c,m)

#首先需要将Trem转换为因子 slimgo$Term=factor(slimgo$Term,levels=slimgo$Term) head(slimgo) colnames(slimgo) p <- ggplot(data = slimgo, mapping = aes(x=Term,y=Percentage,fill=Ontology)) p p+geom_bar(stat="identity") p+geom_bar(stat="identity")+coord_flip() p+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(slimgo$Term))) p+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(slimgo$Term)))+guides(fill=FALSE) p+geom_bar(stat="identity")+coord_flip()+scale_x_discrete(limits=rev(levels(slimgo$Term)))+guides(fill=FALSE)+theme_bw()

#====================================

单独绘制kegg条目图

#==================================== library(ggplot2)

pathway <- read.csv("kegg.csv",header=T) head(pathway) colnames(pathway)

p <- ggplot(data=pathway,mapping = aes(x=richFactor,y=Pathway)) p p + geom_point() p + geom_point(aes(size=AvsB)) p + geom_point(aes(size=AvsB,color=Qvalue)) p + geom_point(aes(size=AvsB,color=Qvalue)) + scale_colour_gradient(low="green",high="red") p + geom_point(aes(size=AvsB,color=Qvalue)) + scale_colour_gradient(low="green",high="red")+labs(title="Top20 of pathway enrichment",x="Rich factor",y="Pathway Name",color="-log10(Qvalue)",size="Gene Numbers") p + geom_point(aes(size=AvsB,color=Qvalue)) + scale_colour_gradient(low="green",high="red")+labs(title="Top20 of pathway enrichment",x="Rich factor",y="Pathway Name",color="-log10(Qvalue)",size="Gene Numbers")+theme_bw()