|
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
|
|
8
|
15 ## some helper fn's
|
|
|
16 write.2.tab <- function( mat,
|
|
|
17 fname ) {
|
|
|
18 mat <- rbind( colnames( mat ), mat )
|
|
|
19 mat <- cbind( c( "ID", rownames( mat )[-1] ),
|
|
|
20 mat )
|
|
|
21 write.table( mat, fname, sep="\t", row.names=FALSE, col.names=FALSE, quote=FALSE )
|
|
|
22 }
|
|
|
23
|
|
0
|
24 lib.load.quiet <- function( package ) {
|
|
|
25 package <- as.character(substitute(package))
|
|
|
26 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
27 }
|
|
|
28
|
|
|
29 lib.load.quiet( getopt )
|
|
|
30 lib.load.quiet( ctc )
|
|
|
31 if ( any( c( 'flashClust', 'fastcluster' ) %in% installed.packages() ) ) {
|
|
|
32 if ( 'flashClust' %in% installed.packages() ) {
|
|
|
33 lib.load.quiet( flashClust )
|
|
|
34 } else {
|
|
|
35 if ( 'fastcluster' %in% installed.packages() ) {
|
|
|
36 lib.load.quiet( fastcluster )
|
|
|
37 }
|
|
|
38 }
|
|
|
39 }
|
|
|
40
|
|
|
41 spec <- matrix( c( "dataset", "d", 1, "character",
|
|
|
42 "dataset2", "D", 2, "character",
|
|
|
43 "output.format", "f", 2, "character",
|
|
|
44 "output.report.dir", "p", 2, "character",
|
|
|
45 "output.fname", "o", 2, "character"
|
|
|
46 ),
|
|
|
47 nc=4,
|
|
|
48 byrow=TRUE
|
|
|
49 )
|
|
|
50
|
|
|
51
|
|
|
52 opt <- getopt( spec=spec )
|
|
|
53 if ( is.null( opt$output.report.dir ) ) { opt$output.report.dir <- "report" }
|
|
|
54 if ( is.null( opt$output.fname ) ) { opt$output.fname <- file.path( opt$output.report.dir, paste( "data", opt$output.format, sep="." ) ) }
|
|
|
55 if ( is.null( opt$output.format ) ) { opt$output.format <- "cdt" }
|
|
|
56
|
|
|
57
|
|
|
58 load( opt$dataset ) ## should load the cl, treecl.res (or partcl.res) and data
|
|
|
59
|
|
|
60 if ( opt$output.format %in% c( "cls-only", "newick" ) ) {
|
|
|
61 if ( opt$output.format == "cls-only" ) {
|
|
|
62
|
|
8
|
63 cl <- matrix( as.numeric( cl ), nc=1, dimnames=list( names(cl), "Class" ) )
|
|
|
64 opt$output.fname <- gsub( "cls-only$", "tab", opt$output.fname )
|
|
0
|
65
|
|
8
|
66 write.2.tab( cl, opt$output.fname )
|
|
0
|
67 } else {
|
|
|
68 ##if ( opt$output.format == "newick" ) {
|
|
|
69
|
|
|
70 if ( ! exists( "treecl.res" ) ) stop( "no HAC result found in results file proved - necessary to generate a Newick formated file.\n" )
|
|
|
71 write( hc2Newick( treecl.res ), opt$output.fname )
|
|
|
72 }
|
|
|
73 } else {
|
|
|
74 if ( ! exists( 'data' ) ) stop( "No data object in the rdata file provided for", opt$output.format, "format!!\n" )
|
|
|
75 if ( inherits( data, "dist" ) ) stop( "data provided is a distance matrix - not a data matrix. Can't generate TreeView or Tab-delimited files w/distance matrices!\n" )
|
|
|
76
|
|
|
77 ## the rest of this is for the remaining output formats
|
|
|
78 ## pre-set the cluster results for rows & cols to NULL
|
|
|
79 hr <- hr.cl <- hc <- hc.cl <- NULL
|
|
|
80 if ( exists( 'treecl.res' ) ) {
|
|
|
81
|
|
|
82 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
|
|
|
83 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
84 hr <- treecl.res
|
|
|
85 hr.cl <- cl
|
|
|
86 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
87 hc <- treecl.res
|
|
|
88 hc.cl <- cl
|
|
|
89 } else {
|
|
|
90 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
91 }
|
|
|
92
|
|
|
93 } else {
|
|
|
94 if ( exists( 'partcl.res' ) ) {
|
|
|
95 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
96 hr <- NA
|
|
|
97 hr.cl <- cl
|
|
|
98 orig.data <- data
|
|
|
99 data <- data[ names( cl ), ] ## partcl.res should now be sorted in order of cluster
|
|
|
100 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
101 hc <- NA
|
|
|
102 hc.cl <- cl
|
|
|
103 orig.data <- data
|
|
|
104 data <- data[ , names( cl ) ] ## partcl.res should now be sorted in order of cluster
|
|
|
105 } else {
|
|
|
106 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
107 }
|
|
|
108 }
|
|
|
109 else {
|
|
|
110 stop( 'could not find a valid cluster result to use for primary direction\n' )
|
|
|
111 }
|
|
|
112 }
|
|
|
113
|
|
|
114
|
|
|
115 if ( ! is.null( opt$dataset2 ) ) {
|
|
|
116
|
|
|
117 ## prep for loading new cluster result
|
|
|
118 if ( ! exists( 'orig.data' ) ) orig.data <- data
|
|
|
119 if ( exists( "treecl.res" ) ) {
|
|
|
120 rm( treecl.res )
|
|
|
121 } else if ( exists( "partcl.res" ) ) {
|
|
|
122 rm( partcl.res )
|
|
|
123 } else stop( "no primary clustering found when generating the 2nd\n" )
|
|
|
124 rm( cl, data )
|
|
|
125
|
|
|
126
|
|
|
127 load( opt$dataset2 ) ## this should bring in the cl obj for the 2nd direction
|
|
|
128
|
|
|
129 ## check the data 1st
|
|
|
130 if ( length( orig.data ) != length( data ) ) stop( "incompatible cluster results in 2nd results file - matrices are diff lengths\n" )
|
|
|
131 if ( nrow( orig.data ) != nrow( data ) ) stop( "incompatible cluster results in 2nd results file - matrices have diff dimensions\n" )
|
|
|
132 if ( sum( orig.data == data ) != length( orig.data ) ) stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
|
|
|
133 ## looks like data is the same, so drop a copy and start chugging
|
|
|
134 rm( orig.data ); gc()
|
|
|
135
|
|
|
136 if ( exists( 'treecl.res' ) ) {
|
|
|
137 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
|
|
|
138
|
|
|
139 if ( is.null( hr ) ) {
|
|
|
140 if ( all( rownames( cl ) %in% rownames( data ) ) ) {
|
|
|
141 hr <- treecl.res
|
|
|
142 hr.cl <- cl
|
|
|
143 } else {
|
|
|
144 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
|
|
|
145 }
|
|
|
146 } else if ( is.null( hc ) ) {
|
|
|
147 if ( all( rownames( cl ) %in% colnames( data ) ) ) {
|
|
|
148 hc <- treecl.res
|
|
|
149 hc.cl <- cl
|
|
|
150 } else {
|
|
|
151 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
|
|
|
152 }
|
|
|
153 } else {
|
|
|
154 stop( "should never get here\n" )
|
|
|
155 }
|
|
|
156 } else if ( exists( 'partcl.res' ) ) {
|
|
|
157 if ( is.null( hr ) ) {
|
|
|
158 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
159 hr <- NA
|
|
|
160 hr.cl <- cl
|
|
|
161 data <- data[ names( cl ), ] ## partcl.res should now be sorted in order of cluster
|
|
|
162 } else {
|
|
|
163 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
|
|
|
164 }
|
|
|
165 } else if ( is.null( hc ) ) {
|
|
|
166 if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
167 hc <- NA
|
|
|
168 hc.cl <- cl
|
|
|
169 data <- data[ , names( cl ) ] ## partcl.res should now be sorted in order of cluster
|
|
|
170 } else {
|
|
|
171 stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
|
|
|
172 }
|
|
|
173 } else {
|
|
|
174 stop( "should never get here\n" )
|
|
|
175 }
|
|
|
176 }
|
|
|
177 }
|
|
|
178
|
|
|
179 ## Now, re-set hc & nr to NULL if they were set to NA
|
|
|
180 ## we used NA to signify that they were set by kmeans/pam, but now, we need to reset them
|
|
|
181 ## for the following lines (that generate the dendrograms (if there was an hclust result)
|
|
|
182 if ( ( !is.null( hr ) ) && is.na( hr ) ) hr <- NULL
|
|
|
183 if ( ( !is.null( hc ) ) && is.na( hc ) ) hc <- NULL
|
|
|
184
|
|
|
185 if ( ! exists( 'data' ) ) stop( "No data object in the rdata file provided!!\n" )
|
|
|
186
|
|
|
187 if ( is.null( hc ) ) hc <- list( order=1:ncol( data ) )
|
|
|
188 if ( is.null( hr ) ) hr <- list( order=1:nrow( data ) )
|
|
|
189
|
|
|
190 if ( opt$output.format == "tabular" ) {
|
|
|
191 write.table( data[ hr$order, hc$order ], opt$output.fname, quote=FALSE, sep="\t", col.names=NA )
|
|
|
192 } else if ( opt$output.format == "cdt" ) {
|
|
|
193 if ( !file.exists( opt$output.report.dir ) ){
|
|
|
194 dir.create(opt$output.report.dir, recursive=T)
|
|
|
195 }
|
|
|
196
|
|
|
197 treeview.fname.stem <- file.path( opt$output.report.dir, "cluster.heatmap")
|
|
|
198 fnames <- character()
|
|
|
199 if ( inherits( hr, "hclust" ) ) {
|
|
|
200 fname <- paste( treeview.fname.stem, ".gtr", sep="" )
|
|
|
201 ## we manually specify a 'stub' distance b/c o/w it'll try using the attr(hr,"method")
|
|
|
202 ## and the r2gtr fn's get grumpy if the distance was anything starting with a 'p'
|
|
|
203 r2gtr( hr, file=fname, distance='stub' )
|
|
|
204 fnames <- c( fnames, fname )
|
|
|
205 }
|
|
|
206 if ( inherits( hc, "hclust" ) ) {
|
|
|
207 fname <- paste( treeview.fname.stem, ".atr", sep="" )
|
|
|
208 r2atr( hc, file=fname, distance='stub' )
|
|
|
209 fnames <- c( fnames, fname )
|
|
|
210 }
|
|
|
211
|
|
|
212 fname <- paste( treeview.fname.stem, ".cdt", sep="" )
|
|
|
213 r2cdt( hr, hc, data, file=fname )
|
|
|
214 fnames <- c( fnames, fname )
|
|
|
215
|
|
|
216 ## jtv file now
|
|
|
217 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>'
|
|
|
218 fname <- paste( treeview.fname.stem, ".jtv", sep="" )
|
|
|
219 cat( jtv.str, file=fname )
|
|
|
220 fnames <- c( fnames, fname )
|
|
|
221
|
|
|
222 cmd <- paste( "tar -zcf",
|
|
|
223 opt$output.fname,
|
|
|
224 paste( "--directory=", opt$output.report.dir, sep="" ),
|
|
|
225 paste( basename( fnames ), collapse=" " ) )
|
|
|
226 system( cmd )
|
|
|
227 }
|
|
|
228 }
|
|
|
229
|
|
|
230
|
|
|
231
|
|
|
232
|
