|
KEGG相关包-clusterProfiler，Pathview的学习
在分析出差异基因后，需对差异基因进行GO和KEGG富集分析，可用clusterProfiler，Pathview包完成，相关学习代码如下。。。
能举一反三即可！！！
相关代码如下：
setwd("c:/Users/Administrator/Desktop/kegg)
rm(list=ls())
source("https:///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") ###改变基因与基因之间的线条
|
