##########下面做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
#加个尾图
#魔法少女小圆~
