Mercurial > repos > ethevenot > heatmap
comparison heatmap_script.R @ 0:ad06aeed02c9 draft
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 |