|
0
|
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 "survival.mode", "M", 2, "character",
|
|
3
|
49 "survival.title", "T", 2, "character"
|
|
0
|
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 opt$out.dir <- 'report'
|
|
|
65 if (! file.exists( opt$out.dir ) ) {
|
|
|
66 dir.create( opt$out.dir )
|
|
|
67 } else {
|
|
|
68 if ( ! file.info( 'report' )$isdir ) {
|
|
|
69 opt$out.dir <- 'heatmap.report'
|
|
|
70 dir.create( opt$out.dir )
|
|
|
71 }
|
|
|
72 }
|
|
|
73
|
|
|
74 if ( opt$image.format == "pdf" ) opt$output.report.html <- file.path( opt$out.dir ,"heatmap.pdf" )
|
|
|
75 if ( opt$image.format == "png" ) opt$output.report.html <- file.path( opt$out.dir ,"index.html" )
|
|
|
76 }
|
|
|
77 if ( is.null( opt$plot.kms ) ) {
|
|
|
78 opt$plot.kms <- FALSE
|
|
|
79 } else {
|
|
|
80 if ( ! opt$plot.kms %in% c( "no", "yes" ) ) {
|
|
|
81 stop( "invalid input to plot.kms param", opt$plot.kms, "\n" )
|
|
|
82 }
|
|
|
83 ## set to TRUE/FALSE
|
|
|
84 opt$plot.kms <- ( opt$plot.kms == "yes" )
|
|
|
85 if ( opt$plot.kms ) {
|
|
|
86 opt$cluster.fname <- opt$dataset
|
|
|
87 if ( is.null( opt$survival.fname ) || ( !file.exists( opt$survival.fname ) ) ) stop( 'must provide a valid file w/clinical data\n' )
|
|
|
88 if ( is.null( opt$survival.script ) || ( !file.exists( opt$survival.script ) ) ) stop( 'must provide a valid path to the gen.survival.curves.R file\n' )
|
|
|
89 if ( is.null(opt$mode ) ) {
|
|
|
90 opt$mode <- "all"
|
|
|
91 } else {
|
|
|
92 if ( ! opt$mode %in% c( 'all', 'one', 'both' ) ) {
|
|
|
93 stop( "invalid mode specified,' -m", opt$mode, "'. must be either {all, one, both}\n" )
|
|
|
94 }
|
|
|
95 }
|
|
|
96 if ( is.null( opt$title ) ) {
|
|
|
97 opt$title <- opt$cluster.fname
|
|
|
98 opt$title <- strsplit( opt$title, "\\/" )[[1]]
|
|
|
99 opt$title <- opt$title[ length( opt$title ) ]
|
|
|
100 }
|
|
|
101 }
|
|
|
102 }
|
|
|
103 if ( is.null( opt$output.treeview ) ) {
|
|
|
104 opt$output.treeview <- FALSE
|
|
|
105 } else {
|
|
|
106 if ( ! opt$output.treeview %in% c( "no", "yes" ) ) {
|
|
|
107 stop( "invalid input to output.treeview param", opt$output.treeview, "\n" )
|
|
|
108 }
|
|
|
109 ## set to TRUE/FALSE
|
|
|
110 opt$output.treeview <- ( opt$output.treeview == "yes" )
|
|
|
111 }
|
|
|
112 if ( is.null( opt$reverse.rows ) ) {
|
|
|
113 opt$reverse.rows <- TRUE
|
|
|
114 } else {
|
|
|
115 if ( ! opt$reverse.rows %in% c( "no", "yes" ) ) {
|
|
|
116 stop( "invalid input to reverse.rows param", opt$reverse.rows, "\n" )
|
|
|
117 }
|
|
|
118
|
|
|
119 ## set to TRUE/FALSE
|
|
|
120 opt$reverse.rows <- ( opt$reverse.rows == "yes" )
|
|
|
121 }
|
|
|
122
|
|
|
123 if ( is.null( opt$second.dir ) ) { opt$second.dir <- "no" }
|
|
|
124 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" )
|
|
|
125 if ( ( opt$image.format == "png" ) || opt$output.treeview ) {
|
|
|
126 if ( !file.exists( opt$output.report.dir ) ){
|
|
|
127 dir.create(opt$output.report.dir, recursive=T)
|
|
|
128 }
|
|
|
129 }
|
|
|
130
|
|
|
131
|
|
|
132
|
|
|
133 load( opt$dataset ) ## should load the cl, treecl.res (or partcl.res) and data
|
|
|
134 ## pre-set the cluster results for rows & cols to NULL
|
|
|
135 hr <- hr.cl <- hc <- hc.cl <- row.ddr <- col.ddr <- NULL
|
|
|
136 if ( exists( 'treecl.res' ) ) {
|
|
|
137
|
|
|
138 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
|
|
|
139 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
140 hr <- treecl.res
|
|
|
141 hr.cl <- cl
|
|
|
142 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
143 hc <- treecl.res
|
|
|
144 hc.cl <- cl
|
|
|
145 } else {
|
|
|
146 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
147 }
|
|
|
148
|
|
|
149 } else {
|
|
|
150 if ( exists( 'partcl.res' ) ) {
|
|
|
151 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
152 hr <- NA
|
|
|
153 hr.cl <- cl
|
|
|
154 orig.data <- data
|
|
|
155 data <- data[ names( cl ), ] ## partcl.res should now be sorted in order of cluster
|
|
|
156 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
157 hc <- NA
|
|
|
158 hc.cl <- cl
|
|
|
159 orig.data <- data
|
|
|
160 data <- data[ , names( cl ) ] ## partcl.res should now be sorted in order of cluster
|
|
|
161 } else {
|
|
|
162 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
163 }
|
|
|
164 }
|
|
|
165 else {
|
|
|
166 stop( 'could not find a valid cluster result to use for primary direction\n' )
|
|
|
167 }
|
|
|
168 }
|
|
|
169
|
|
|
170
|
|
|
171 if ( opt$second.dir %in% c( "yes", "prev" ) ) {
|
|
|
172
|
|
|
173 if ( opt$second.dir == "yes" ) {
|
|
|
174 if ( is.null( hr ) ) {
|
|
|
175 hr <- hclust( dist( data ) )
|
|
|
176 } else if ( is.null( hc ) ) {
|
|
|
177 hc <- hclust( dist( t( data ) ) )
|
|
|
178 }
|
|
|
179 } else { ## opt$second.dir == "prev"
|
|
|
180
|
|
|
181 ## prep for loading new cluster result
|
|
|
182 if ( ! exists( 'orig.data' ) ) orig.data <- data
|
|
|
183 if ( exists( "treecl.res" ) ) {
|
|
|
184 rm( treecl.res )
|
|
|
185 } else if ( exists( "partcl.res" ) ) {
|
|
|
186 rm( partcl.res )
|
|
|
187 } else stop( "no primary clustering found when generating the 2nd\n" )
|
|
|
188 rm( cl, data )
|
|
|
189
|
|
|
190
|
|
|
191 load( opt$dataset2 ) ## this should bring in the cl obj for the 2nd direction
|
|
|
192
|
|
|
193 ## check the data 1st
|
|
|
194 if ( length( orig.data ) != length( data ) ) stop( "incompatible cluster results in 2nd results file - matrices are diff lengths\n" )
|
|
|
195 if ( nrow( orig.data ) != nrow( data ) ) stop( "incompatible cluster results in 2nd results file - matrices have diff dimensions\n" )
|
|
|
196 if ( any( is.na( orig.data ) ) ) {
|
|
|
197 nas <- which( is.na( orig.data ) )
|
|
|
198 num.nas <- length( nas )
|
|
|
199 ## 1st, chk the NAs
|
|
|
200 if ( sum( which( is.na( data ) ) != nas ) == num.nas ) stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
201 if ( ( sum( orig.data == data, na.rm=T )+num.nas ) != length( orig.data ) ) {
|
|
|
202 stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
203 }
|
|
|
204
|
|
|
205 } else {
|
|
|
206 if ( sum( orig.data == data ) != length( orig.data ) ) stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
207 }
|
|
|
208 ## looks like data is the same, so drop a copy and start chugging
|
|
|
209 rm( orig.data ); gc()
|
|
|
210
|
|
|
211 if ( exists( 'treecl.res' ) ) {
|
|
|
212 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
|
|
|
213
|
|
|
214 if ( is.null( hr ) ) {
|
|
|
215 if ( all( rownames( cl ) %in% rownames( data ) ) ) {
|
|
|
216 hr <- treecl.res
|
|
|
217 hr.cl <- cl
|
|
|
218 } else {
|
|
|
219 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
|
|
|
220 }
|
|
|
221 } else if ( is.null( hc ) ) {
|
|
|
222 if ( all( rownames( cl ) %in% colnames( data ) ) ) {
|
|
|
223 hc <- treecl.res
|
|
|
224 hc.cl <- cl
|
|
|
225 } else {
|
|
|
226 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
|
|
|
227 }
|
|
|
228 } else {
|
|
|
229 stop( "should never get here\n" )
|
|
|
230 }
|
|
|
231 } else if ( exists( 'partcl.res' ) ) {
|
|
|
232 if ( is.null( hr ) ) {
|
|
|
233 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
234 hr <- NA
|
|
|
235 hr.cl <- cl
|
|
|
236 data <- data[ names( cl ), ] ## partcl.res should now be sorted in order of cluster
|
|
|
237 } else {
|
|
|
238 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
|
|
|
239 }
|
|
|
240 } else if ( is.null( hc ) ) {
|
|
|
241 if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
242 hc <- NA
|
|
|
243 hc.cl <- cl
|
|
|
244 data <- data[ , names( cl ) ] ## partcl.res should now be sorted in order of cluster
|
|
|
245 } else {
|
|
|
246 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
|
|
|
247 }
|
|
|
248 } else {
|
|
|
249 stop( "should never get here\n" )
|
|
|
250 }
|
|
|
251 }
|
|
|
252 }
|
|
|
253 }
|
|
|
254
|
|
|
255 ## Now, re-set hc & nr to NULL if they were set to NA
|
|
|
256 ## we used NA to signify that they were set by kmeans/pam, but now, we need to reset them
|
|
|
257 ## for the following lines (that generate the dendrograms (if there was an hclust result)
|
|
|
258 if ( ( !is.null( hr ) ) && is.na( hr ) ) hr <- NULL
|
|
|
259 if ( ( !is.null( hc ) ) && is.na( hc ) ) hc <- NULL
|
|
|
260
|
|
|
261 if ( ! is.null( hr ) ) {
|
|
|
262 row.ddr <- as.dendrogram( hr )
|
|
|
263 ## need this to make sure that the heatmap is oriented the same way as it is in TreeView
|
|
|
264 if ( opt$reverse.rows ) row.ddr <- rev( row.ddr )
|
|
|
265 }
|
|
|
266
|
|
|
267 if ( ! is.null( hc ) ) {
|
|
|
268 col.ddr <- as.dendrogram( hc )
|
|
|
269 }
|
|
|
270
|
|
|
271
|
|
|
272 hmcols<-colorRampPalette(c("blue","white","red"))(256)
|
|
|
273
|
|
|
274 if ( ( ! is.null( row.ddr ) ) && ( ! is.null( col.ddr ) ) ) {
|
|
|
275 dend.param <- "both"
|
|
|
276 } else {
|
|
|
277 dend.param <- "none"
|
|
|
278 if ( ! is.null( row.ddr ) ) dend.param <- "row"
|
|
|
279 if ( ! is.null( col.ddr ) ) dend.param <- "column"
|
|
|
280 }
|
|
|
281
|
|
|
282
|
|
|
283 param.list <- list( x=data,
|
|
|
284 Rowv=row.ddr,
|
|
|
285 Colv=col.ddr,
|
|
|
286 dendrogram=dend.param,
|
|
|
287 trace="none",
|
|
|
288 col=hmcols,
|
|
|
289 symbreaks=TRUE,
|
|
|
290 scale="none",
|
|
|
291 labRow="",
|
|
|
292 labCol="",
|
|
|
293 na.color='grey' ) #,
|
|
|
294 ##key=FALSE )
|
|
|
295
|
|
|
296 if ( ! is.null( hr.cl ) ) {
|
|
|
297 hrcols <- rainbow( max( as.numeric( hr.cl ) ) )
|
|
|
298 names( hrcols ) <- sort( unique( as.numeric( hr.cl ) ) )
|
|
|
299 rowColLabs <- hrcols[ as.character( as.numeric( hr.cl ) ) ]
|
|
|
300 param.list <- c( param.list, list( RowSideColors=rowColLabs ) )
|
|
|
301 }
|
|
|
302 if ( ! is.null( hc.cl ) ) {
|
|
|
303 hccols <- rainbow( max( as.numeric( hc.cl ) ) )
|
|
|
304 names( hccols ) <- sort( unique( as.numeric( hc.cl ) ) )
|
|
|
305 colColLabs <- hccols[ as.character( as.numeric( hc.cl ) ) ]
|
|
|
306 param.list <- c( param.list, list( ColSideColors=colColLabs ) )
|
|
|
307 }
|
|
|
308
|
|
|
309
|
|
|
310 if ( opt$image.format == 'png' ) {
|
|
|
311 png.fname <- file.path( opt$output.report.dir, "cluster.heatmap.png")
|
|
|
312 plot.dev <- png( png.fname,
|
|
|
313 width=8.5,
|
|
|
314 height=11,
|
|
|
315 units='in',
|
|
|
316 res=72 )
|
|
|
317 } else {
|
|
|
318 pdf.fname <- opt$output.report.html
|
|
|
319 pdf( opt$output.report.html,
|
|
|
320 paper="letter" )
|
|
|
321 }
|
|
|
322
|
|
|
323 do.call( "heatmap.2", param.list )
|
|
|
324
|
|
|
325 dev.off() ## close the previous device
|
|
|
326
|
|
|
327 if ( opt$plot.kms ) {
|
|
|
328
|
|
|
329 cmd.string <- opt$survival.script
|
|
|
330
|
|
|
331 ## get the consensusClass file that's associated with the k.select
|
|
|
332 cmd.string <- paste( cmd.string, "-C", opt$dataset )
|
|
|
333 cmd.string <- paste( cmd.string, "-S", opt$survival.fname )
|
|
|
334 cmd.string <- paste( cmd.string, "-M", opt$survival.mode )
|
|
3
|
335 if ( ! is.null( opt$suvival.title ) )
|
|
|
336 cmd.string <- paste( cmd.string, "-T", opt$suvival.title )
|
|
0
|
337
|
|
|
338 ## only call kms if we're the image is png
|
|
|
339 if ( opt$image.format=="png" ) {
|
|
|
340 png.fname <- file.path( opt$output.report.dir, "kaplan.meier.survival.png")
|
|
|
341 cmd.string <- paste( cmd.string, "-I", "png" )
|
|
|
342 cmd.string <- paste( cmd.string, "-O", png.fname )
|
|
|
343 system( cmd.string )
|
|
|
344 }
|
|
|
345 }
|
|
|
346
|
|
|
347
|
|
|
348 if ( opt$image.format == 'png' ) {
|
|
|
349 pngs = list.files(path=opt$output.report.dir, patt="png")
|
|
|
350 html.out <- paste( "<html>",
|
|
|
351 paste( paste( paste( "<div><img src=\'", pngs, sep="" ), "\'/></div>", sep="" ), collapse=""),
|
|
|
352 "</html>" )
|
|
|
353 cat( html.out, file=opt$output.report.html )
|
|
|
354 }
|
|
|
355
|
|
|
356
|
|
|
357 if ( opt$output.treeview ) {
|
|
|
358 treeview.fname.stem <- file.path( opt$output.report.dir, "cluster.heatmap")
|
|
|
359 fnames <- character()
|
|
|
360 if ( ! is.null( hr ) ) {
|
|
|
361 fname <- paste( treeview.fname.stem, ".gtr", sep="" )
|
|
|
362 r2gtr( hr, file=fname )
|
|
|
363 fnames <- c( fnames, fname )
|
|
|
364 } else {
|
|
|
365 hr <- list( order=1:nrow( data ) )
|
|
|
366 }
|
|
|
367 if ( ! is.null( hc ) ) {
|
|
|
368 fname <- paste( treeview.fname.stem, ".atr", sep="" )
|
|
|
369 r2atr( hc, file=fname )
|
|
|
370 fnames <- c( fnames, fname )
|
|
|
371 } else {
|
|
|
372 hc <- list( order=1:ncol( data ) )
|
|
|
373 }
|
|
|
374
|
|
|
375
|
|
|
376 fname <- paste( treeview.fname.stem, ".cdt", sep="" )
|
|
|
377 r2cdt( hr, hc, data, file=fname )
|
|
|
378 fnames <- c( fnames, fname )
|
|
|
379
|
|
|
380 ## jtv file now
|
|
|
381 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>'
|
|
|
382 fname <- paste( treeview.fname.stem, ".jtv", sep="" )
|
|
|
383 cat( jtv.str, file=fname )
|
|
|
384 fnames <- c( fnames, fname )
|
|
|
385
|
|
|
386 cmd <- paste( "tar -zcf", opt$output.fname, paste( "--directory=", opt$output.report.dir, sep="" ), paste( basename( fnames ), collapse=" " ) )
|
|
|
387 system( cmd )
|
|
|
388 }
|
|
|
389
|
|
|
390
|
