clusterProfiler——GO和KEGG分析之R代码

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)) ###信号通路个数
length(unique(hsa_kegg$KEGGPATHID2EXTID$to))  ###所有信号通路中的基因数
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)  ###信号通路个数  
### 一个信号通路一个名字，部分信号通路无基因？？？

###
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三个基因个数为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数 unique的基因数为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,后面为对应差异基因倍数。

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数，y为信号通路，颜色为p值
dotplot(kk, showCategory=12)  ##x为geneRatio数
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")  ###改变基因与基因之间的线条