第 30 课-KEGG 分析

###############KEGG分析#################################
rm(list=ls())   #清空环境变量
options(stringsAsFactors = F)
load("TCGA-STADDEG_symbol.Rdata")

#1. 添加ENTREZID列
library(clusterProfiler)
library(org.Hs.eg.db)
library(tidyverse)
genelist <- bitr(DEG$Gene, fromType="SYMBOL",
                 toType="ENTREZID", OrgDb='org.Hs.eg.db')
DEG <- inner_join(DEG,genelist,by=c("Gene"="SYMBOL"))
gene_up = DEG[DEG$change == 'UP','ENTREZID'] 
gene_down = DEG[DEG$change == 'DOWN','ENTREZID'] 
gene_diff = c(gene_up,gene_down)

#2. KEGG分析
kk <- enrichKEGG(gene         =  gene_diff,
                 organism     = 'hsa',
                 pvalueCutoff = 0.1,
                 qvalueCutoff =0.1)
write.table(kk,file="KEGG.txt",sep="\t",quote=F,row.names = F)
save(kk, file = "kegg_TCGA-STAD.Rdata")

#3. 可视化
load("kegg_TCGA-STAD.Rdata")
##3.1 柱状图
barplot(kk, showCategory = 20,color = "pvalue")
##3.2 气泡图
dotplot(kk, showCategory = 20)
##3.3 网络图
List = DEG$log2FoldChange
names(List)= DEG$ENTREZID
List = sort(List,decreasing = T)
cnetplot(kk, categorySize="pvalue", foldChange = List,colorEdge = TRUE)
##3.4 网络图2
cnetplot(kk, foldChange = List, circular = TRUE, colorEdge = TRUE)

#4. 绘制通路图
library("pathview")
keggxls=read.table("KEGG.txt",sep="\t",header=T)
pvalue_up = DEG[DEG$change == 'UP','pvalue'] 
pvalue_down = DEG[DEG$change == 'DOWN','pvalue'] 
pvalue_diff = c(pvalue_up,pvalue_down)
names(pvalue_diff)=gene_diff

##4.1 单个样本通路图可视化
hsa04080 <- pathview(gene.data  = -log10(pvalue_diff),
                     pathway.id = "hsa04080",
                     species    = "hsa",
                     out.suffix = "pathview", 
                     limit      = list(gene= 10, cpd= 1))

##4.2 所有通路图可视化
for(i in keggxls$ID){
  pv.out <- pathview(gene.data = -log10(pvalue_diff), 
                     pathway.id = i, 
                     species = "hsa", 
                     out.suffix = "pathview",
                     limit=list(gene=10, cpd=10),
                     bins = list(gene = 10, cpd= 10))
}