基因功能富集分析-R语言

##安装bioconductor上的包；
source(“http://bioconductor.org/biocLite.R”)
biocLite（“clusterprofiler”）
biocLite("org.Hs.eg.db")#人基因名称等信息包；
##加载clusterprofiler包到当前工作路径；
library（clusterprofiler）#基因富集分析用；
library（org.Hs.eg.db）
#读入需要分析的数据，包含一列基因名称的列表；
a <- read.table(file.choose(),header = F,colClasses = "character")
#a <- read.table(file.choose(),header = F,colClasses = c("V1"= "character")),只设置第一列值为字符型；
###选取基因列的所有行
b <- a[,1]
###利用bitr函数将基因名称转换为ENTREZID号；物种是人org.Hs.eg.db；
eg = bitr(b,fromType = "SYMBOL",toType = "ENTREZID",OrgDb = "org.Hs.eg.db")
#可能会有部分基因对应不到ENTREZID，0.4% of input gene IDs are fail to map...
###转换后的基因名称保存为文档；
write.table(eg,file = "test_id.txt")
gene <- eg[,2]
###进行GO和KEGG分析；
library(clusterProfiler)
library(org.Hs.eg.db)
a <- read.table(file.choose(),header = F,colClasses = c("V1"= "character"))
b <- a[,1]
eg <- bitr(b,fromType = "SYMBOL",toType = "ENTREZID",OrgDb = "org.Hs.eg.db")
gene <- eg[,2]
ego_CC <- enrichGO(gene = gene,
                   OrgDb=org.Hs.eg.db,
                   ont = "CC",
                   pAdjustMethod = "BH",
                   minGSSize = 1,
                   pvalueCutoff = 0.01,
                   qvalueCutoff = 0.01,
                   readable = TRUE)
write.csv(as.data.frame(ego_CC),row.names = F, file = "ego_CC.csv")
barplot(ego_CC,drop = TRUE,title = "enrichment_CC",showCategory = 12)

ego_BP <- enrichGO(gene = gene,
                   OrgDb=org.Hs.eg.db,
                   ont = "BP",
                   pAdjustMethod = "BH",
                   minGSSize = 1,
                   pvalueCutoff = 0.01,
                   qvalueCutoff = 0.01,
                   readable = TRUE)
write.csv(as.data.frame(ego_BP),row.names = F, file = "ego_BP.csv")
barplot(ego_BP,drop = TRUE,title = "enrichment_BP",showCategory = 12)

ego_MF <- enrichGO(gene = gene,
                   OrgDb=org.Hs.eg.db,
                   ont = "MF",
                   pAdjustMethod = "BH",
                   minGSSize = 1,
                   pvalueCutoff = 0.01,
                   qvalueCutoff = 0.01,
                   readable = TRUE)
write.csv(as.data.frame(ego_MF),row.names = F, file = "ego_MF.csv")
barplot(ego_MF,drop = TRUE,title = "enrichment_MF",showCategory = 12)

kk <- enrichKEGG(gene = gene,
                 organism ="hsa",
                 pvalueCutoff = 0.01,
                 qvalueCutoff = 0.01,
                 minGSSize = 1,
                 #readable = TRUE ,
                 use_internal_data = FALSE)
write.csv(as.data.frame(kk),row.names = F, file = "kk.csv")
barplot(kk,drop = TRUE,title = "enrichment_kegg",showCategory = 12)

###DisGeNET4 is an integrative and comprehensive resources of gene-disease associations from several public data sources and the literature. It contains gene-disease associations and snp-gene-disease associations.
###The enrichment analysis of disease-gene associations is supported by the enrichDGN function and analysis of snp-gene-disease associations is supported by the enrichDGNv function.
dgn <- enrichDGN(gene = gene,
                 pAdjustMethod = "BH",
                 pvalueCutoff = 0.05,
                 qvalueCutoff = 0.05,
                 readable = TRUE)
head(dgn)

write.csv(as.data.frame(dgn),row.names = F, file = "dgn.csv")
barplot(dgn,drop = TRUE,title = "enrichment_disease",showCategory = 12)

###条行图，按p值从小到大排列；
barplot(ego_CC,showCategory = 24,title = "EnrichmentGO_CC")
###点状图，按富集数从大到小进行排列；
dotplot(ego_CC,title = "EnrichenmentGo_CC")