rm(list = ls())
options(stringsAsFactors = F)

library(clusterProfiler)
library(org.Bt.eg.db)


DEG <- read.csv('output/magenta_muscle_genes.csv')
head(DEG)
DEG = DEG[ ,3]
gene_all <- read.table('ids.txt')
head(gene_all)
gene_all = gene_all[ ,1]

keytypes(org.Bt.eg.db)

# bitr in clusterProfiler
allID <- bitr(gene_all, fromType = "ENSEMBL", 
              toType = c( "ENTREZID" ), 
              OrgDb = org.Bt.eg.db )
head(allID)
degID <- bitr(DEG, fromType = "SYMBOL", 
              toType = c( "ENTREZID" ), 
              OrgDb = org.Bt.eg.db )
head(degID)

class(degID)
# KEGG analysis----
enrich <- enrichKEGG(gene =degID[,2],
                     organism='bta',
                     universe=allID[,2],
                     pvalueCutoff=1,
                     qvalueCutoff=1)


GeneRatio <- as.numeric(lapply(strsplit(enrich$GeneRatio,split="/"),function(x) 
  as.numeric(x[1])/as.numeric(x[2])))
head(GeneRatio)

BgRatio <- as.numeric(lapply(strsplit(enrich$BgRatio,split="/"),function(x) 
  as.numeric(x[1])/as.numeric(x[2])  ))
head(BgRatio)

enrich_factor <- GeneRatio/BgRatio

out <- data.frame(enrich$ID,
                  enrich$Description,
                  enrich$GeneRatio,
                  enrich$BgRatio,
                  round(enrich_factor,2),
                  enrich$pvalue,
                  enrich$qvalue,
                  enrich$geneID)

colnames(out) <- c("ID","Description","GeneRatio","BgRatio","enrich_factor","pvalue","qvalue","geneID")
write.table(out,"magenta-muscle_enrich_KEGG.xls",row.names = F,sep="\t",quote = F)

out_sig0.05 <- out[out$qvalue<0.01,]

# barplot
bar <- barplot(enrich,showCategory=10,title="KEGG Pathway",
               colorBy="p.adjust")
bar

# save
pdf(file = "magenta-muscle_kegg_bar_plot.pdf",width = 6,height = 8)
print(bar)
dev.off()

# dotplot
dot <- dotplot(enrich,x="geneRatio",showCategory=10,font.size=12,title="KEGG Pathway")
dot

# save
pdf(file = "magenta-muscle_dot_plot.pdf",width = 7,height = 7)
print(dot)
dev.off()


# GO 
enrich <- enrichGO(gene =degID[,2],OrgDb='org.Bt.eg.db',
                   ont="BP",universe=allID[,2],pvalueCutoff=1,qvalueCutoff=1)


GeneRatio <- as.numeric(lapply(strsplit(enrich$GeneRatio,split="/"),function(x) 
  as.numeric(x[1])/as.numeric(x[2])))

BgRatio <- as.numeric(lapply(strsplit(enrich$BgRatio,split="/"),function(x) 
  as.numeric(x[1])/as.numeric(x[2])))

enrich_factor <- GeneRatio/BgRatio

out <- data.frame(enrich$ID,
                  enrich$Description,
                  enrich$GeneRatio,
                  enrich$BgRatio,
                  round(enrich_factor,2),
                  enrich$pvalue,
                  enrich$qvalue,
                  enrich$geneID)

colnames(out) <- c("ID","Description","GeneRatio","BgRatio","enrich_factor","pvalue","qvalue","geneID")
write.table(out,"magenta-muscle_enrich_GO_BP.xls",row.names = F,sep="\t",quote = F)

out_sig0.05 <- out[out$qvalue<0.01,]


# barplot
bar <- barplot(enrich,showCategory=10,title="Biological Pathway",colorBy="p.adjust")
bar

# save
pdf(file = "magenta-muscle_BP_bar_plot.pdf",width = 8,height = 7)
print(bar)
dev.off()

# dotplot
dot <- dotplot(enrich,x="geneRatio",showCategory=10,font.size=12,title="Biological Pathway")
dot

# save
pdf(file = "magenta-muscle_BP_dot_plot.pdf",width = 6,height = 6)
print(dot)
dev.off()