KEGG相關包-clusterProfiler，Pathview的學習
在分析出差異基因后，需對差異基因進行GO和KEGG富集分析，可用clusterProfiler，Pathview包完成，相關學習代碼如下。。。
能舉一反三即可?。?！
相關代碼如下：
setwd("c:/Users/Administrator/Desktop/kegg)
rm(list=ls())
source("https://bioconductor.org/biocLite.R")
options(BioC_mirror="http://mirrors.ustc.edu.cn/bioc/")
biocLite("DOSE")
biocLite("topGO")
biocLite("clusterProfiler")
biocLite("pathview")
##
library(DOSE)
library(GO.db)
library(org.Hs.eg.db)
library(topGO)
library(GSEABase)
library(clusterProfiler)
hsa_kegg<-download_KEGG("hsa")
str(hsa_kegg)
length(unique(hsa_kegg$KEGGPATHID2EXTID$from)) ###信號通路個數(shù)
length(unique(hsa_kegg$KEGGPATHID2EXTID$to)) ###所有信號通路中的基因數(shù)
length(unique(hsa_kegg$KEGGPATHID2NAME$from)) ###信號通路對應的名字 為啥與1個length不等；
length(unique(hsa_kegg$KEGGPATHID2NAME$to)) ###信號通路對應的名字，部分信號通路無基因？？
x<-data.frame(hsa_kegg$KEGGPATHID2EXTID) ###兩列為pathwayid,基因id extid 擴增標識符
#xxx#一個信號通路含有多個基因，一個基因在多個信號通路
y<-data.frame(hsa_kegg$KEGGPATHID2NAME) ###信號通路個數(shù)
### 一個信號通路一個名字，部分信號通路無基因？？？
###
hsa_go<-org.Hs.egGO ###GO分析
str(hsa_go)
head(hsa_go@L2Rchain)
mapped_Go_genes<-mappedkeys(hsa_go)
length(mapped_Go_genes) #GO（BP,CC，MF三個基因個數(shù)為18622個gene）
hsa_kegg2<-org.Hs.egPATH
str(hsa_kegg2)
mapped_Kegg_genes<-mappedkeys(hsa_kegg2)
length(mapped_Kegg_genes)
length(unique(mapped_Kegg_genes))
###KEGG里面PATH 的Gene數(shù) unique的基因數(shù)為5869個基因
###
data(geneList, package="DOSE")
str(geneList)
# a<-as.data.frame(data(geneList, package="DOSE")) ##錯誤
# b<-data(geneList, package="DOSE")
# ##錯誤
c<-as.data.frame(geneList) ##可以，基因為entrez gene id,后面為對應差異基因倍數(shù)。
gene <- names(geneList)[abs(geneList) > 2]
d<-as.data.frame(gene)
kk <- enrichKEGG(gene = gene,
organism = 'hsa',
pvalueCutoff = 0.05)
head(kk)
length(kk$ID)
barplot(kk, drop=TRUE, showCategory=12) ##x軸為基因counts數(shù)，y為信號通路，顏色為p值
dotplot(kk, showCategory=12) ##x為geneRatio數(shù)
enrichMap(kk)
cnetplot(kk, categorySize="pvalue", foldChange=geneList)
MM <- setReadable(kk, OrgDb = org.Hs.eg.db,keytype = "ENTREZID")
cnetplot(MM, categorySize="pvalue", foldChange=geneList)
##
library(pathview)
gene.with.fc<-geneList[abs(geneList)>2]
e<-as.data.frame(geneList)
hsa04110 <- pathview(gene.data = gene.with.fc,
pathway.id = "hsa04110",
species = "hsa",
out.suffix = "fc",
kegg.native=T)
hsa04110.21layer<-pathview(gene.data = gene.with.fc,
pathway.id = "hsa04110",
species = "hsa",
out.suffix = "fc.21layer",
limit=list(gene=max(abs(gene.with.fc)),cpd=1),
kegg.native=TRUE,same.layer=F) ###將部分基因ID轉換成gene name;
hsa04110.graphviz<-pathview(gene.data = gene.with.fc,
pathway.id = "hsa04110",
species = "hsa",
out.suffix = "fc.graphviz",
limit=list(gene=max(abs(gene.with.fc)),cpd=1),
kegg.native=F,sign.pos="bottomleft") ###改變基因與基因之間的線條
