Mercurial > repos > ethevenot > heatmap

comparison heatmap_script.R @ 0:ad06aeed02c9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/workflow4metabolomics/heatmap.git commit 7e599d006e53fefb7e1b923ba8894b4fb19f9cfa-dirty
author ethevenot
date Tue, 02 Aug 2016 06:26:41 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:ad06aeed02c9
1 ## Etienne Thevenot
2 ## CEA, MetaboHUB
3 ## W4M Core Development Team
4 ## etienne.thevenot@cea.fr
5 ## 2015-05-30
6
7
8 heatmapF <- function(proMN,
9 obsDF,
10 feaDF,
11 disC, ## dissimilarity
12 cutSamN, ## number of sample clusters
13 cutVarN, ## number of variable clusters
14 fig.pdfC,
15 corMetC, ## correlation method
16 aggMetC, ## agglomeration method
17 colC, ## color scale
18 scaL,
19 cexN) {
20
21 ncaN <- 14 ## Sample and variable name truncature for display
22
23 if(aggMetC == "ward") {
24
25 rvsLs <- R.Version()
26 aggMetC <- paste0(aggMetC,
27 ifelse(as.numeric(rvsLs[["major"]]) > 3 ||
28 as.numeric(rvsLs[["major"]]) == 3 && as.numeric(rvsLs[["minor"]]) > 0.3,
29 ".D",
30 ""))
31
32 }
33
34 if(disC %in% c("euclidean", "maximum", "manhattan", "canberra", "binary", "minkowski")) {
35
36 obsHcl <- hclust(dist(proMN, method = disC),
37 method = aggMetC)
38
39 feaHcl <- hclust(dist(t(proMN), method = disC),
40 method = aggMetC)
41
42 } else if(disC == "1-cor") {
43
44 obsHcl <- hclust(as.dist(1-cor(t(proMN),
45 method = corMetC,
46 use = "pairwise.complete.obs")),
47 method = aggMetC)
48
49 feaHcl <- hclust(as.dist(1-cor(proMN,
50 method = corMetC,
51 use = "pairwise.complete.obs")),
52 method = aggMetC)
53
54 } else if(disC == "1-abs(cor)") {
55
56 obsHcl <- hclust(as.dist(1-abs(cor(t(proMN),
57 method = corMetC,
58 use = "pairwise.complete.obs"))),
59 method = aggMetC)
60
61 feaHcl <- hclust(as.dist(1-abs(cor(proMN,
62 method = corMetC,
63 use = "pairwise.complete.obs"))),
64 method = aggMetC)
65
66 }
67
68 heaMN <- proMN <- proMN[obsHcl[["order"]], feaHcl[["order"]]]
69
70 if(scaL)
71 heaMN <- scale(heaMN)
72
73 heaMN <- heaMN[, rev(1:ncol(heaMN)), drop = FALSE]
74
75 switch(colC,
76 blueOrangeRed = {
77 imaPalVn <- colorRampPalette(c("blue", "orange", "red"),
78 space = "rgb")(5)[1:5]
79 },
80 redBlackGreen = {
81 imaPalVn <- colorRampPalette(c("red", "black", "green"),
82 space = "rgb")(5)[1:5]
83 })
84
85
86 ## figure
87 ##-------
88
89 pdf(fig.pdfC,
90 width = 14,
91 height = 14)
92
93 layout(matrix(1:4, nrow = 2),
94 widths = c(1, 9), heights = c(1, 9))
95
96 ## Color scale
97
98 scaN <- length(imaPalVn)
99
100 par(mar = c(0.6, 0.6, 0.6, 4.1))
101
102 ylimVn <- c(0, scaN)
103 ybottomVn <- 0:(scaN - 1)
104 ytopVn <- 1:scaN
105
106 plot(x = 0,
107 y = 0,
108 bty = "n",
109 font.axis = 2,
110 font.lab = 2,
111 type = "n",
112 xlim = c(0, 1),
113 ylim = ylimVn,
114 xlab = "",
115 ylab = "",
116 xaxs = "i",
117 yaxs = "i",
118 xaxt = "n",
119 yaxt = "n")
120
121 rect(xleft = 0.8,
122 ybottom = ybottomVn,
123 xright = 1,
124 ytop = ytopVn,
125 col = imaPalVn,
126 border = NA)
127
128 prtVn <- pretty(range(heaMN, na.rm = TRUE))
129 axis(at = scaN / diff(range(prtVn)) * (prtVn - min(prtVn)),
130 font = 2,
131 font.axis = 2,
132 labels = prtVn,
133 las = 1,
134 lwd = 2,
135 lwd.ticks = 2,
136 side = 4,
137 xpd = TRUE)
138
139 arrows(par("usr")[2],
140 par("usr")[4],
141 par("usr")[2],
142 par("usr")[3],
143 code = 0,
144 lwd = 2,
145 xpd = TRUE)
146
147 ## Feature dendrogram
148
149 par(mar = c(7.1, 0.6, 0, 0.1),
150 lwd = 2)
151
152 plot(rev(as.dendrogram(feaHcl)), horiz = TRUE,
153 leaflab = "none",
154 main = "", xaxs = "i", yaxs = "i",
155 xaxt = "n", yaxt = "n", xlab = "", ylab = "")
156
157 revFeaHcl <- list(merge = cbind(feaHcl[["merge"]][, 2], feaHcl[["merge"]][, 1]),
158 height = feaHcl[["height"]],
159 order = rev(feaHcl[["order"]]),
160 labels = feaHcl[["labels"]])
161
162 if(cutVarN > 1) {
163 cluFeaVn <- cutree(revFeaHcl, k = cutVarN)[revFeaHcl[["order"]]]
164 cutFeaVn <- which(abs(diff(cluFeaVn)) > 0)
165 cutFeaTxtVn <- c(cutFeaVn[1] / 2, cutFeaVn + diff(c(cutFeaVn, length(cluFeaVn))) / 2) + 0.5
166 cutFeaLinVn <- cutFeaVn + 0.5
167 text(par("usr")[1] + 0.2 * diff(par("usr")[1:2]),
168 cutFeaTxtVn,
169 labels = unique(cluFeaVn),
170 cex = 2,
171 font = 2,
172 las = 2)
173 }
174
175 ## Observation dendrogram
176
177 par(mar = c(0.1, 0, 0.6, 7.1),
178 lwd = 2)
179
180 plot(as.dendrogram(obsHcl), leaflab = "none",
181 main = "", xaxs = "i", yaxs = "i",
182 yaxt = "n", xlab = "", ylab = "")
183
184 if(cutSamN > 1) {
185 cluObsVn <- cutree(obsHcl, k = cutSamN)[obsHcl[["order"]]]
186 cutObsVn <- which(abs(diff(cluObsVn)) > 0)
187 cutObsTxtVn <- c(cutObsVn[1] / 2, cutObsVn + diff(c(cutObsVn, length(cluObsVn))) / 2) + 0.5
188 cutObsLinVn <- cutObsVn + 0.5
189 text(cutObsTxtVn,
190 0.8 * par("usr")[4],
191 labels = unique(cluObsVn),
192 cex = 2,
193 font = 2)
194 }
195
196 ## Heatmap
197
198 par(mar = c(7.1, 0, 0, 7.1))
199
200 image(x = 1:nrow(heaMN),
201 y = 1:ncol(heaMN),
202 z = round(heaMN),
203 col = imaPalVn,
204 font.axis = 2,
205 font.lab = 2,
206 xaxt = "n",
207 yaxt = "n",
208 xlab = "",
209 ylab = "")
210
211 obsOrdVc <- obsHcl[["labels"]][obsHcl[["order"]]]
212 obsOrdLenVn <- sapply(obsOrdVc, nchar)
213 obsOrdVc <- substr(obsOrdVc, 1, ncaN)
214 obsOrdVc <- paste0(obsOrdVc, ifelse(obsOrdLenVn > ncaN, ".", ""), " ")
215
216 mtext(obsOrdVc,
217 at = 1:nrow(heaMN),
218 cex = cexN,
219 las = 2,
220 side = 1)
221
222 feaOrdVc <- feaHcl[["labels"]][feaHcl[["order"]]]
223 feaOrdLenVn <- sapply(feaOrdVc, nchar)
224 feaOrdVc <- substr(feaOrdVc, 1, ncaN)
225 feaOrdVc <- paste0(" ", feaOrdVc, ifelse(feaOrdLenVn > ncaN, ".", ""))
226
227 mtext(feaOrdVc,
228 at = ncol(heaMN):1,
229 cex = cexN,
230 las = 2,
231 side = 4)
232
233 if(cutVarN > 1)
234 abline(h = cutFeaLinVn)
235 if(cutSamN > 1)
236 abline(v = cutObsLinVn)
237
238 box()
239
240 dev.off()
241
242 ## Returning
243 ##----------
244
245 if(cutSamN > 1)
246 obsDF[, "heat_clust"] <- cutree(obsHcl, k = cutSamN)
247 obsDF <- obsDF[obsHcl[["order"]], , drop = FALSE]
248
249 if(cutVarN > 1)
250 feaDF[, "heat_clust"] <- cutree(feaHcl, k = cutVarN)
251 feaDF <- feaDF[feaHcl[["order"]], , drop = FALSE]
252
253 return(invisible(list(proMN = proMN,
254 obsDF = obsDF,
255 feaDF = feaDF)))
256
257
258 } ## end of heatmapF