|
0
|
1 #!/usr/bin/env Rscript
|
|
|
2 argspec <- c("tab.2.cdt.R converts a data matrix to cdt format
|
|
|
3
|
|
|
4 Usage:
|
|
|
5
|
|
|
6 Optional:
|
|
|
7
|
|
|
8 \n\n")
|
|
|
9 args <- commandArgs(TRUE)
|
|
|
10 if ( length( args ) == 1 && args =="--help") {
|
|
|
11 write(argspec, stderr())
|
|
|
12 q();
|
|
|
13 }
|
|
|
14
|
|
|
15 lib.load.quiet <- function( package ) {
|
|
|
16 package <- as.character(substitute(package))
|
|
|
17 suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
|
|
|
18 }
|
|
|
19
|
|
|
20 lib.load.quiet( getopt )
|
|
|
21 lib.load.quiet( amap )
|
|
|
22
|
|
|
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 spec <- matrix( c( "dataset", "d", 1, "character",
|
|
|
34 "gen.new.dgram", "g", 2, "character",
|
|
|
35 "output.fname", "o", 2, "character"
|
|
|
36 ),
|
|
|
37 nc=4,
|
|
|
38 byrow=TRUE
|
|
|
39 )
|
|
|
40
|
|
|
41
|
|
|
42 opt <- getopt( spec=spec )
|
|
|
43 if ( is.null( opt$output.report.dir ) ) {
|
|
|
44 opt$output.report.dir <- "report"
|
|
|
45
|
|
|
46 if (! file.exists( opt$output.report.dir ) ) {
|
|
|
47 dir.create( opt$output.report.dir )
|
|
|
48 } else {
|
|
|
49 if ( ! file.info( 'report' )$isdir ) {
|
|
|
50 opt$output.report.dir <- 'heatmap.report'
|
|
|
51 dir.create( opt$output.report.dir )
|
|
|
52 }
|
|
|
53 }
|
|
|
54 }
|
|
|
55 if ( is.null( opt$output.fname ) ) { opt$output.fname <- file.path( opt$output.report.dir, paste( "data.RData", sep="." ) ) }
|
|
|
56 if ( is.null( opt$gen.new.dgram ) ) {
|
|
|
57 opt$gen.new.dgram <- FALSE
|
|
|
58 } else {
|
|
|
59 if ( ! opt$gen.new.dgram %in% c( "no", "yes" ) ) {
|
|
|
60 stop( "invalid input to gen.new.dgram param", opt$gen.new.dgram, "\n" )
|
|
|
61 }
|
|
|
62 ## set to TRUE/FALSE
|
|
|
63 opt$gen.new.dgram <- ( opt$gen.new.dgram == "yes" )
|
|
|
64 }
|
|
|
65
|
|
|
66
|
|
|
67 load( opt$dataset ) ## should load the cl, treecl.res (or partcl.res) and data
|
|
|
68
|
|
|
69 if ( ! exists( 'data' ) ) stop( "No data object in the rdata file provided for", opt$output.format, "format!!\n" )
|
|
|
70 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" )
|
|
|
71
|
|
|
72 ## the rest of this is for the remaining output formats
|
|
|
73 ## pre-set the cluster results for rows & cols to NULL
|
|
|
74 direction <- NULL
|
|
|
75 if ( exists( 'treecl.res' ) ) {
|
|
|
76 cl.res <- treecl.res
|
|
|
77 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
|
|
|
78 } else {
|
|
|
79 if ( exists( 'partcl.res' ) ) {
|
|
|
80 cl.res <- partcl.res
|
|
|
81 }
|
|
|
82 else {
|
|
|
83 stop( 'could not find a valid cluster result to use for primary direction\n' )
|
|
|
84 }
|
|
|
85 }
|
|
|
86
|
|
|
87 if ( all( names( cl ) %in% rownames( data ) ) ) {
|
|
|
88 direction <- "rows"
|
|
|
89 } else if ( all( names( cl ) %in% colnames( data ) ) ) {
|
|
|
90 direction <- "cols"
|
|
|
91 data <- t( data )
|
|
|
92 } else {
|
|
|
93 stop( "Specified cluster result does not come from this data set\n" )
|
|
|
94 }
|
|
|
95
|
|
|
96
|
|
|
97 centroids <- NULL
|
|
|
98 cl <- sort( cl )
|
|
|
99 if ( inherits( cl.res, "kmeans" ) ) {
|
|
|
100 ## already comes pre-calculated for us!!
|
|
|
101 centroids <- cl.res$centers
|
|
|
102 } else {
|
|
|
103 data <- data[ names( cl ), ]
|
|
|
104 cl.list <- unique( cl )
|
|
|
105 cl.list <- lapply( cl.list, function(i) cl[ cl %in% i ] )
|
|
|
106 centroids <- sapply( cl.list,
|
|
|
107 function(x) {
|
|
|
108 return( apply( data[ names(x), , drop=F ], 2, mean, na.rm=T ) )
|
|
|
109 }
|
|
|
110 )
|
|
|
111 centroids <- t( centroids ) ## get them back to the same number of columns that data has now
|
|
|
112 }
|
|
|
113
|
|
|
114 data <- centroids
|
|
|
115 rownames( data ) <- sapply( 1:max( cl ), function(i) sprintf( "cluster-%02d", i ) )
|
|
|
116
|
|
|
117 if ( opt$gen.new.dgram ) {
|
|
|
118 distance <- 'euclidean'
|
|
|
119 if ( inherits( cl.res, 'hclust' ) ) {
|
|
|
120 distance <- cl.res$dist.method
|
|
|
121 }
|
|
|
122 amap.distance <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
|
|
|
123 "pearson", "abspearson", "correlation", "abscorrelation", "spearman", "kendall" )
|
|
|
124 names( amap.distance ) <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
|
|
|
125 "cosine", "abscosine", "pearson", "abspearson", "spearman", "kendall" )
|
|
|
126
|
|
|
127 if ( ! distance %in% names( amap.distance ) ) stop("unsupported distance.")
|
|
|
128 dist.mat <- Dist( data, method=as.character( amap.distance[ distance ] ) )
|
|
|
129 treecl.res <- hclust( dist.mat )
|
|
|
130 cl <- cutree( treecl.res, nrow(data) )
|
|
|
131 }
|
|
|
132
|
|
|
133 if ( direction == "cols" ) {
|
|
|
134 data <- t( data )
|
|
|
135 }
|
|
|
136
|
|
|
137 save( file=opt$output.fname, treecl.res, cl, data )
|
