-
Notifications
You must be signed in to change notification settings - Fork 0
/
enrichment.R
89 lines (72 loc) · 3.32 KB
/
enrichment.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
enrichment <- function(geneList_, fileName, outputPath = F){
# # library(UniProt.ws)
# up <- UniProt.ws(taxId=9606)
# select(up, geneList.uniprot, "ENTREZ_GENE")
convert2Symbol <- function(x, map){
l<-unlist(strsplit(x, split="/"))
map <- select(org.Hs.eg.db, l, "SYMBOL", "ENTREZID")
return(paste(map$SYMBOL, collapse = "/"))
}
map <- select(org.Hs.eg.db, geneList_, "ENTREZID", "SYMBOL")
geneList <- na.omit(map$ENTREZID)
ego <- enrichGO(gene = geneList,
OrgDb = org.Hs.eg.db,
ont = "BP",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05)
tmp <- data.frame(ego)
genes <- sapply(tmp$geneID,convert2Symbol)
tmp$geneSymbols <- genes
write.csv(tmp , file = paste0("./functional_analysis/",fileName,"_GO_BP.csv"), row.names=FALSE)
ego <- enrichGO(gene = geneList,
OrgDb = org.Hs.eg.db,
ont = "MF",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05)
tmp <- data.frame(ego)
genes <- sapply(tmp$geneID,convert2Symbol)
tmp$geneSymbols <- genes
write.csv(tmp , file = paste0("./functional_analysis/",fileName,"_GO_MF.csv"), row.names=FALSE)
ego <- enrichGO(gene = geneList,
OrgDb = org.Hs.eg.db,
ont = "CC",
pAdjustMethod = "fdr",
pvalueCutoff = 0.05)
tmp <- data.frame(ego)
genes <- sapply(tmp$geneID,convert2Symbol)
tmp$geneSymbols <- genes
write.csv(tmp , file = paste0("./functional_analysis/",fileName,"_GO_CC.csv"), row.names=FALSE)
# -------------------------------------------------------------------------
ekegg <- enrichKEGG(gene = geneList,
#keyType = 'uniprot',
pAdjustMethod = "fdr",
pvalueCutoff = 0.05)
tmp <- data.frame(ekegg)
genes <- sapply(tmp$geneID,convert2Symbol)
tmp$geneSymbols <- genes
write.csv(tmp , file = paste0("./functional_analysis/",fileName,"_Kegg.csv") , row.names=FALSE)
# PathView ----------------------------------------------------------------
# Not needed at this stage
pathViewOut <- function(keggID, geneEntrez){
gene.data <- rep(1, length(geneEntrez))
names(gene.data) <- geneEntrez
msg <- capture.output(pathview::pathview(gene.data = geneEntrez,
pathway.id = keggID,
limit = list(gene = 2, cpd = 1), bins = list(gene = 20, cpd= 10),
low = list(gene = "green", cpd = "blue"),
species = "hsa"), type = "message")
msg <- grep("image file", msg, value = T)
filename <- sapply(strsplit(msg, " "), function(x) x[length(x)])
img <- png::readPNG(filename)
png::writePNG(img, target = paste0("./functional_analysis/",fileName,"_", filename))
invisible(file.remove(filename))
}
# if(outputPath & nrow(tmp)>0){
# N = 5
# if(N>nrow(tmp))N=nrow(tmp)
# for (i in 1:nrow(tmp)) {
# print(tmp$ID[i])
# pathViewOut(tmp$ID[i], geneEntrez = geneList)
# }
# }
}