24-有参转录组实战10-差异基因KEGG富集分析

##########下面做KEGG########
emapper <- read.delim("out.emapper.annotations")
emapper[emapper=="-"] <- NA#change "-" to "NA"
emapper <- emapper[-(49584:49586),]#remove the final 3 rows
DE <- read.delim("DE_genes_filter.txt")

pathway2gene <- dplyr::select(emapper,Pathway=KEGG_Pathway, GID=query)%>%
  separate_rows(Pathway, sep = ',', convert = F)%>%
  filter(str_detect(Pathway, 'ko'))%>%
  mutate(Pathway= str_remove(Pathway, 'ko'))
pathway2gene <- pathway2gene%>%dplyr::mutate(Pathway=paste0("map",Pathway))#前面加"map"
library(magrittr)
get_path2name <- function(){
  keggpathid2name.df <- clusterProfiler:::kegg_list("pathway")
  keggpathid2name.df[,1] %>% gsub("path:map", "", .)
  colnames(keggpathid2name.df) <- c("path_id","path_name")
  return(keggpathid2name.df)
}
pathway2name <- get_path2name()

#查看下变量pathway2gene和pathway2name

ekp <- enricher(gene = DE$GID,
                TERM2GENE = pathway2gene,
                TERM2NAME = pathway2name,
                pvalueCutoff = 0.05,
                qvalueCutoff = 0.05)
write.table(ekp, file = "Ptri_KEGG_test",sep = '\t',quote = F)
#这里富集出来了，通路很少，可以自己调整P值和Q值的参数，或者在差异基因的筛选条件上放宽一点。Dsecription是对应https://rest.kegg.jp/list/pathway/#另外GeneRatio做成百分比

ekp2 <- read.delim("Ptri_KEGG_test")
library(ggplot2)
pp <- ggplot(ekp2, aes(GeneRatio, Description))
pp + geom_point() +
  geom_point(aes(size = Count)) +
  geom_point(aes(size = Count, color = -1 * log10(qvalue))) +
  scale_colour_gradient(low = "green", high = "red") +
  labs(color = expression(-log[10](Qvalue)), size = "Gene Number", x = "Rich Factor", y = "Pathway Name", title = "KEGG Pathway Enrichment") +
  theme_bw()
#图片自己慢慢修，更多图表自行百度，可参考教程https://zhuanlan.zhihu.com/p/377356510