|
2
|
1 #!/usr/bin/env Rscript
|
|
|
2 argspec <- c("tab.2.cdt.R converts a data matrix to cdt format
|
|
|
3
|
|
|
4 Usage:
|
|
|
5 tab.2.cdt.R -d <data.file>
|
|
|
6 Optional:
|
|
|
7 -o <output_file>
|
|
|
8 \n\n")
|
|
|
9 args <- commandArgs(TRUE)
|
|
|
10 if ( length( args ) == 1 && args =="--help") {
|
|
|
11 write(argspec, stderr())
|
|
|
12 q();
|
|
|
13 }
|
|
|
14
|
|
|
15
|
|
|
16 lib.load.quiet <- function( package ) {
|
|
|
17 package <- as.character(substitute(package))
|
|
|
18 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
19 }
|
|
|
20 lib.load.quiet(getopt)
|
|
|
21 lib.load.quiet( gplots )
|
|
|
22 lib.load.quiet( ctc )
|
|
|
23 if ( any( c( 'flashClust', 'fastcluster' ) %in% installed.packages() ) ) {
|
|
|
24 if ( 'flashClust' %in% installed.packages() ) {
|
|
|
25 lib.load.quiet( flashClust )
|
|
|
26 } else {
|
|
|
27 if ( 'fastcluster' %in% installed.packages() ) {
|
|
|
28 lib.load.quiet( fastcluster )
|
|
|
29 }
|
|
|
30 }
|
|
|
31 }
|
|
|
32
|
|
|
33
|
|
|
34 spec <- matrix( c( "dataset", "d", 1, "character",
|
|
|
35 "second.dir", "s", 2, "character",
|
|
|
36 "dataset2", "D", 2, "character",
|
|
|
37 "reverse.rows", "r", 2, "character",
|
|
|
38 "image.format", "i", 2, "character",
|
|
|
39 "plot.kms", "k", 2, "character",
|
|
|
40 "output.fname", "o", 2, "character",
|
|
|
41 "output.report.html", "h", 2, "character",
|
|
|
42 "output.report.dir", "p", 2, "character",
|
|
|
43 "output.treeview", "t", 2, "character",
|
|
|
44 "survival.script", "z", 2, "character",
|
|
|
45 "cluster.fname", "C", 2, "character",
|
|
|
46 "survival.fname", "S", 2, "character",
|
|
|
47 "survival.image", "I", 2, "character",
|
|
|
48 "mode", "M", 2, "character",
|
|
|
49 "title", "T", 2, "character"
|
|
|
50 ),
|
|
|
51 nc=4,
|
|
|
52 byrow=TRUE
|
|
|
53 )
|
|
|
54
|
|
|
55
|
|
|
56 opt <- getopt( spec=spec )
|
|
|
57 if ( is.null( opt$image.format ) ){
|
|
|
58 opt$image.format <- "png"
|
|
|
59 } else {
|
|
|
60 if ( ! opt$image.format %in% c( "pdf", "png" ) ) stop( 'invalid image format specified\n' )
|
|
|
61 }
|
|
|
62 if ( is.null( opt$output.report.dir ) ) { opt$output.report.dir <- "report" }
|
|
|
63 if ( is.null( opt$output.report.html ) ) {
|
|
|
64 if ( opt$image.format == "pdf" ) opt$output.report.html <- "report/heatmap.pdf"
|
|
|
65 if ( opt$image.format == "png" ) opt$output.report.html <- "report/index.html"
|
|
|
66 }
|
|
|
67 if ( is.null( opt$plot.kms ) ) {
|
|
|
68 opt$plot.kms <- TRUE
|
|
|
69 } else {
|
|
|
70 if ( ! opt$plot.kms %in% c( "no", "yes" ) ) {
|
|
|
71 stop( "invalid input to plot.kms param", opt$plot.kms, "\n" )
|
|
|
72 }
|
|
|
73 ## set to TRUE/FALSE
|
|
|
74 opt$plot.kms <- ( opt$plot.kms == "yes" )
|
|
|
75 if ( opt$plot.kms ) {
|
|
|
76 opt$cluster.fname <- opt$dataset
|
|
|
77 if ( is.null( opt$survival.script ) || ( !file.exists( opt$survival.script ) ) ) stop( 'must provide a valid file w/clinical data\n' )
|
|
|
78 if ( is.null(opt$mode ) ) {
|
|
|
79 opt$mode <- "all"
|
|
|
80 } else {
|
|
|
81 if ( ! opt$mode %in% c( 'all', 'one', 'both' ) ) {
|
|
|
82 stop( "invalid mode specified,' -m", opt$mode, "'. must be either {all, one, both}\n" )
|
|
|
83 }
|
|
|
84 }
|
|
|
85 if ( is.null( opt$title ) ) {
|
|
|
86 opt$title <- opt$cluster.fname
|
|
|
87 opt$title <- strsplit( opt$title, "\\/" )[[1]]
|
|
|
88 opt$title <- opt$title[ length( opt$title ) ]
|
|
|
89 }
|
|
|
90 }
|
|
|
91 }
|
|
|
92 if ( is.null( opt$output.treeview ) ) {
|
|
|
93 opt$output.treeview <- FALSE
|
|
|
94 } else {
|
|
|
95 if ( ! opt$output.treeview %in% c( "no", "yes" ) ) {
|
|
|
96 stop( "invalid input to output.treeview param", opt$output.treeview, "\n" )
|
|
|
97 }
|
|
|
98 ## set to TRUE/FALSE
|
|
|
99 opt$output.treeview <- ( opt$output.treeview == "yes" )
|
|
|
100 }
|
|
|
101 if ( is.null( opt$reverse.rows ) ) {
|
|
|
102 opt$reverse.rows <- TRUE
|
|
|
103 } else {
|
|
|
104 if ( ! opt$reverse.rows %in% c( "no", "yes" ) ) {
|
|
|
105 stop( "invalid input to reverse.rows param", opt$reverse.rows, "\n" )
|
|
|
106 }
|
|
|
107
|
|
|
108 ## set to TRUE/FALSE
|
|
|
109 opt$reverse.rows <- ( opt$reverse.rows == "yes" )
|
|
|
110 }
|
|
|
111
|
|
|
112 if ( is.null( opt$second.dir ) ) { opt$second.dir <- "no" }
|
|
|
113 if ( is.null( opt$second.dir ) && is.null( opt$dataset2 ) ) stop( "must specify an rdata file to load if a previous result is to be used to cluster the 2nd direction\n" )
|
|
|
114 if ( ( opt$image.format == "png" ) || opt$output.treeview ) {
|
|
|
115 if ( !file.exists( opt$output.report.dir ) ){
|
|
|
116 dir.create(opt$output.report.dir, recursive=T)
|
|
|
117 }
|
|
|
118 }
|
|
|
119
|
|
|
120
|
|
|
121
|
|
|
122 load( opt$dataset ) ## should load the cl, treecl.res (or partcl.res) and data
|
|
|
123 ## pre-set the cluster results for rows & cols to NULL
|
|
|
124 hr <- hr.cl <- hc <- hc.cl <- row.ddr <- col.ddr <- NULL
|
|
|
125 if ( exists( 'treecl.res' ) ) {
|
|
|
126
|
|
|
127 if ( is.null( treecl.res$dist.method ) ) treecl.res$dist.method <- 'euclidean' # just set it to some stub so that the ctc fn's don't complain
|
|
|
128 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
129 hr <- treecl.res
|
|
|
130 hr.cl <- cl
|
|
|
131 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
132 hc <- treecl.res
|
|
|
133 hc.cl <- cl
|
|
|
134 } else {
|
|
|
135 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
136 }
|
|
|
137
|
|
|
138 } else {
|
|
|
139 if ( exists( 'partcl.res' ) ) {
|
|
|
140 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
141 hr <- NA
|
|
|
142 hr.cl <- cl
|
|
|
143 orig.data <- data
|
|
|
144 data <- data[ names( cl ), ] ## partcl.res should now be sorted in order of cluster
|
|
|
145 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
146 hc <- NA
|
|
|
147 hc.cl <- cl
|
|
|
148 orig.data <- data
|
|
|
149 data <- data[ , names( cl ) ] ## partcl.res should now be sorted in order of cluster
|
|
|
150 } else {
|
|
|
151 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
152 }
|
|
|
153 }
|
|
|
154 else {
|
|
|
155 stop( 'could not find a valid cluster result to use for primary direction\n' )
|
|
|
156 }
|
|
|
157 }
|
|
|
158
|
|
|
159
|
|
|
160 if ( opt$second.dir %in% c( "yes", "prev" ) ) {
|
|
|
161
|
|
|
162 if ( opt$second.dir == "yes" ) {
|
|
|
163 if ( is.null( hr ) ) {
|
|
|
164 hr <- hclust( dist( data ) )
|
|
|
165 } else if ( is.null( hc ) ) {
|
|
|
166 hc <- hclust( dist( t( data ) ) )
|
|
|
167 }
|
|
|
168 } else { ## opt$second.dir == "prev"
|
|
|
169
|
|
|
170 ## prep for loading new cluster result
|
|
|
171 if ( ! exists( 'orig.data' ) ) orig.data <- data
|
|
|
172 if ( exists( "treecl.res" ) ) {
|
|
|
173 rm( treecl.res )
|
|
|
174 } else if ( exists( "partcl.res" ) ) {
|
|
|
175 rm( partcl.res )
|
|
|
176 } else stop( "no primary clustering found when generating the 2nd\n" )
|
|
|
177 rm( cl, data )
|
|
|
178
|
|
|
179
|
|
|
180 load( opt$dataset2 ) ## this should bring in the cl obj for the 2nd direction
|
|
|
181
|
|
|
182 ## check the data 1st
|
|
|
183 if ( length( orig.data ) != length( data ) ) stop( "incompatible cluster results in 2nd results file - matrices are diff lengths\n" )
|
|
|
184 if ( nrow( orig.data ) != nrow( data ) ) stop( "incompatible cluster results in 2nd results file - matrices have diff dimensions\n" )
|
|
|
185 if ( any( is.na( orig.data ) ) ) {
|
|
|
186 nas <- which( is.na( orig.data ) )
|
|
|
187 num.nas <- length( nas )
|
|
|
188 ## 1st, chk the NAs
|
|
|
189 if ( sum( which( is.na( data ) ) != nas ) == num.nas ) stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
190 if ( ( sum( orig.data == data, na.rm=T )+num.nas ) != length( orig.data ) ) {
|
|
|
191 stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
192 }
|
|
|
193
|
|
|
194 } else {
|
|
|
195 if ( sum( orig.data == data ) != length( orig.data ) ) stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
196 }
|
|
|
197 ## looks like data is the same, so drop a copy and start chugging
|
|
|
198 rm( orig.data ); gc()
|
|
|
199
|
|
|
200 if ( exists( 'treecl.res' ) ) {
|
|
|
201 if ( is.null( treecl.res$dist.method ) ) treecl.res$dist.method <- 'euclidean' # just set it to some stub so that the ctc fn's don't complain
|
|
|
202
|
|
|
203 if ( is.null( hr ) ) {
|
|
|
204 if ( all( rownames( cl ) %in% rownames( data ) ) ) {
|
|
|
205 hr <- treecl.res
|
|
|
206 hr.cl <- cl
|
|
|
207 } else {
|
|
|
208 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
|
|
|
209 }
|
|
|
210 } else if ( is.null( hc ) ) {
|
|
|
211 if ( all( rownames( cl ) %in% colnames( data ) ) ) {
|
|
|
212 hc <- treecl.res
|
|
|
213 hc.cl <- cl
|
|
|
214 } else {
|
|
|
215 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
|
|
|
216 }
|
|
|
217 } else {
|
|
|
218 stop( "should never get here\n" )
|
|
|
219 }
|
|
|
220 } else if ( exists( 'partcl.res' ) ) {
|
|
|
221 if ( is.null( hr ) ) {
|
|
|
222 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
223 hr <- NA
|
|
|
224 hr.cl <- cl
|
|
|
225 data <- data[ names( cl ), ] ## partcl.res should now be sorted in order of cluster
|
|
|
226 } else {
|
|
|
227 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
|
|
|
228 }
|
|
|
229 } else if ( is.null( hc ) ) {
|
|
|
230 if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
231 hc <- NA
|
|
|
232 hc.cl <- cl
|
|
|
233 data <- data[ , names( cl ) ] ## partcl.res should now be sorted in order of cluster
|
|
|
234 } else {
|
|
|
235 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
|
|
|
236 }
|
|
|
237 } else {
|
|
|
238 stop( "should never get here\n" )
|
|
|
239 }
|
|
|
240 }
|
|
|
241 }
|
|
|
242 }
|
|
|
243
|
|
|
244 ## Now, re-set hc & nr to NULL if they were set to NA
|
|
|
245 ## we used NA to signify that they were set by kmeans/pam, but now, we need to reset them
|
|
|
246 ## for the following lines (that generate the dendrograms (if there was an hclust result)
|
|
|
247 if ( ( !is.null( hr ) ) && is.na( hr ) ) hr <- NULL
|
|
|
248 if ( ( !is.null( hc ) ) && is.na( hc ) ) hc <- NULL
|
|
|
249
|
|
|
250 if ( ! is.null( hr ) ) {
|
|
|
251 row.ddr <- as.dendrogram( hr )
|
|
|
252 ## need this to make sure that the heatmap is oriented the same way as it is in TreeView
|
|
|
253 if ( opt$reverse.rows ) row.ddr <- rev( row.ddr )
|
|
|
254 }
|
|
|
255
|
|
|
256 if ( ! is.null( hc ) ) {
|
|
|
257 col.ddr <- as.dendrogram( hc )
|
|
|
258 }
|
|
|
259
|
|
|
260
|
|
|
261 hmcols<-colorRampPalette(c("blue","white","red"))(256)
|
|
|
262
|
|
|
263 if ( ( ! is.null( row.ddr ) ) && ( ! is.null( col.ddr ) ) ) {
|
|
|
264 dend.param <- "both"
|
|
|
265 } else {
|
|
|
266 dend.param <- "none"
|
|
|
267 if ( ! is.null( row.ddr ) ) dend.param <- "row"
|
|
|
268 if ( ! is.null( col.ddr ) ) dend.param <- "column"
|
|
|
269 }
|
|
|
270
|
|
|
271
|
|
|
272 param.list <- list( x=data,
|
|
|
273 Rowv=row.ddr,
|
|
|
274 Colv=col.ddr,
|
|
|
275 dendrogram=dend.param,
|
|
|
276 trace="none",
|
|
|
277 col=hmcols,
|
|
|
278 scale="none",
|
|
|
279 labRow="",
|
|
|
280 labCol="" ) #,
|
|
|
281 ##key=FALSE )
|
|
|
282
|
|
|
283 if ( ! is.null( hr.cl ) ) {
|
|
|
284 hrcols <- rainbow( max( as.numeric( hr.cl ) ) )
|
|
|
285 names( hrcols ) <- sort( unique( as.numeric( hr.cl ) ) )
|
|
|
286 rowColLabs <- hrcols[ as.character( as.numeric( hr.cl ) ) ]
|
|
|
287 param.list <- c( param.list, list( RowSideColors=rowColLabs ) )
|
|
|
288 }
|
|
|
289 if ( ! is.null( hc.cl ) ) {
|
|
|
290 hccols <- rainbow( max( as.numeric( hc.cl ) ) )
|
|
|
291 names( hccols ) <- sort( unique( as.numeric( hc.cl ) ) )
|
|
|
292 colColLabs <- hccols[ as.character( as.numeric( hc.cl ) ) ]
|
|
|
293 param.list <- c( param.list, list( ColSideColors=colColLabs ) )
|
|
|
294 }
|
|
|
295
|
|
|
296
|
|
|
297 if ( opt$image.format == 'png' ) {
|
|
|
298 png.fname <- file.path( opt$output.report.dir, "cluster.heatmap.png")
|
|
|
299 plot.dev <- png( png.fname,
|
|
|
300 width=8.5,
|
|
|
301 height=11,
|
|
|
302 units='in',
|
|
|
303 res=72 )
|
|
|
304 } else {
|
|
|
305 pdf.fname <- opt$output.report.html
|
|
|
306 pdf( opt$output.report.html,
|
|
|
307 paper="letter" )
|
|
|
308 }
|
|
|
309
|
|
|
310 do.call( "heatmap.2", param.list )
|
|
|
311 if ( opt$plot.kms ) {
|
|
|
312
|
|
|
313 if ( opt$image.format == 'png' ) {
|
|
|
314 dev.off() ## close the previous device
|
|
|
315 png.fname <- file.path( opt$output.report.dir, "kaplan.meier.survival.png")
|
|
|
316 plot.dev <- png( png.fname,
|
|
|
317 width=8.5,
|
|
|
318 height=11,
|
|
|
319 units='in',
|
|
|
320 res=72 )
|
|
|
321 }
|
|
|
322
|
|
|
323 bak <- opt$image.format; opt$image.format <- 'none'
|
|
|
324 source( opt$survival.script )
|
|
|
325 opt$image.format <- bak
|
|
|
326 }
|
|
|
327
|
|
|
328 dev.off()
|
|
|
329
|
|
|
330 if ( opt$image.format == 'png' ) {
|
|
|
331 pngs = list.files(path=opt$output.report.dir, patt="png")
|
|
|
332 html.out <- paste( "<html>",
|
|
|
333 paste( paste( paste( "<div><img src=\'", pngs, sep="" ), "\'/></div>", sep="" ), collapse=""),
|
|
|
334 "</html>" )
|
|
|
335 cat( html.out, file=opt$output.report.html )
|
|
|
336 }
|
|
|
337
|
|
|
338
|
|
|
339 if ( opt$output.treeview ) {
|
|
|
340 treeview.fname.stem <- file.path( opt$output.report.dir, "cluster.heatmap")
|
|
|
341 fnames <- character()
|
|
|
342 if ( ! is.null( hr ) ) {
|
|
|
343 fname <- paste( treeview.fname.stem, ".gtr", sep="" )
|
|
|
344 r2gtr( hr, file=fname )
|
|
|
345 fnames <- c( fnames, fname )
|
|
|
346 } else {
|
|
|
347 hr <- list( order=1:nrow( data ) )
|
|
|
348 }
|
|
|
349 if ( ! is.null( hc ) ) {
|
|
|
350 fname <- paste( treeview.fname.stem, ".atr", sep="" )
|
|
|
351 r2atr( hc, file=fname )
|
|
|
352 fnames <- c( fnames, fname )
|
|
|
353 } else {
|
|
|
354 hc <- list( order=1:ncol( data ) )
|
|
|
355 }
|
|
|
356
|
|
|
357
|
|
|
358 fname <- paste( treeview.fname.stem, ".cdt", sep="" )
|
|
|
359 r2cdt( hr, hc, data, file=fname )
|
|
|
360 fnames <- c( fnames, fname )
|
|
|
361
|
|
|
362 ## jtv file now
|
|
|
363 jtv.str <- '<DocumentConfig><UrlExtractor/><ArrayUrlExtractor/><Views><View type="Dendrogram" dock="1"><ColorExtractor contrast="2.0"><ColorSet zero="#FFFFFF" down="#0000FF"/></ColorExtractor><ArrayDrawer/><GlobalXMap current="Fill"><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></GlobalXMap><GlobalYMap current="Fill"><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></GlobalYMap><ZoomXMap><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></ZoomXMap><ZoomYMap><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></ZoomYMap><TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView></TextView><ArrayNameView face="Monospaced" size="14"><ArraySummary included="0"/></ArrayNameView><AtrSummary/><GtrSummary/></View></Views></DocumentConfig>'
|
|
|
364 fname <- paste( treeview.fname.stem, ".jtv", sep="" )
|
|
|
365 cat( jtv.str, file=fname )
|
|
|
366 fnames <- c( fnames, fname )
|
|
|
367
|
|
|
368 cmd <- paste( "tar -zcf", opt$output.fname, paste( "--directory=", opt$output.report.dir, sep="" ), paste( basename( fnames ), collapse=" " ) )
|
|
|
369 system( cmd )
|
|
|
370 }
|
|
|
371
|
|
|
372
|
