/
enrichMap.R
111 lines (97 loc) · 3.07 KB
/
enrichMap.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
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
##' enrichment map
##'
##' enrichment map
##' @title enrichMap
##' @param x gseaResult or enrichResult object
##' @param n maximum number of category to shown
##' @param fixed if set to FALSE, will invoke tkplot
##' @param vertex.label.font font size of vertex label
##' @param ... additional parameter
##' @return figure
##' @importFrom reshape2 melt
##' @importFrom igraph delete.edges
##' @importFrom igraph get.edgelist
##' @importFrom igraph add_vertices
##' @importFrom igraph graph.empty
##' @export
##' @author Guangchuang Yu
enrichMap <- function(x, n = 50, fixed=TRUE, vertex.label.font=1, ...) {
if (is(x, "gseaResult")) {
geneSets <- x@geneSets
}
if (is(x, "enrichResult")) {
geneSets <- geneInCategory(x)
}
y <- as.data.frame(x)
if (nrow(y) < n) {
n <- nrow(y)
}
y <- y[1:n,]
if (n == 0) {
stop("no enriched term found...")
} else if (n == 1) {
g <- graph.empty(0, directed=FALSE)
g <- add_vertices(g, nv = 1)
V(g)$name <- y$Description
V(g)$color <- "red"
} else {
pvalue <- y$pvalue
id <- y[,1]
geneSets <- geneSets[id]
n <- nrow(y) #
w <- matrix(NA, nrow=n, ncol=n)
colnames(w) <- rownames(w) <- y$Description
for (i in 1:n) {
for (j in i:n) {
w[i,j] = overlap_ratio(geneSets[id[i]], geneSets[id[j]])
}
}
wd <- melt(w)
wd <- wd[wd[,1] != wd[,2],]
wd <- wd[!is.na(wd[,3]),]
g <- graph.data.frame(wd[,-3], directed=F)
E(g)$width=sqrt(wd[,3]*20)
g <- delete.edges(g, E(g)[wd[,3] < 0.2])
idx <- unlist(sapply(V(g)$name, function(x) which(x == y$Description)))
cols <- color_scale("red", "#E5C494")
V(g)$color <- cols[sapply(pvalue, getIdx, min(pvalue), max(pvalue))]
## seq_gradient_pal("red", "grey")(pvalue[idx])
## data can be exported to view in Cytoscape or other tools
##
##
## Edata <- as.data.frame(get.edgelist(g))
## Edata$edgewidth <- E(g)$width
## Vdata <- data.frame(pathway=V(g)$name, color=V(g)$color)
## map_data <- list(edge_data=Edata, vertex_data=Vdata)
if (is(x, "gseaResult")) {
cnt <- y$setSize / 10
}
if (is(x, "enrichResult")) {
cnt <- y$Count
}
names(cnt) <- y$Description
cnt2 <- cnt[V(g)$name]
V(g)$size <- log(cnt2, base=10) * 10 ## cnt2/sum(cnt2) * 100
}
netplot(g, vertex.label.font=vertex.label.font, vertex.label.color="black", fixed=fixed, ...)
## invisible(map_data)
invisible(g)
}
overlap_ratio <- function(x, y) {
x <- unlist(x)
y <- unlist(y)
length(intersect(x, y))/length(unique(c(x,y)))
}
##' @importFrom grDevices colorRampPalette
color_scale <- function(c1="grey", c2="red") {
pal <- colorRampPalette(c(c1, c2))
colors <- pal(100)
return(colors)
}
getIdx <- function(v, MIN, MAX) {
if ( MIN == MAX ) {
return(100)
}
intervals <- seq(MIN, MAX, length.out=100)
max(which(intervals <= v))
}