diff cluster.tools/cluster.2.centroid.R @ 0:0decf3fd54bc draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cluster.tools/cluster.2.centroid.R	Thu Feb 28 01:45:39 2013 -0500
@@ -0,0 +1,137 @@
+#!/usr/bin/env Rscript
+argspec <- c("tab.2.cdt.R converts a data matrix to cdt format
+
+        Usage: 
+
+        Optional:
+
+                \n\n")
+args <- commandArgs(TRUE)
+if ( length( args ) == 1 && args =="--help") { 
+  write(argspec, stderr())
+  q();
+}
+
+lib.load.quiet <- function( package ) {
+   package <- as.character(substitute(package))
+   suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
+}
+
+lib.load.quiet( getopt )
+lib.load.quiet( amap )
+
+if ( any( c( 'flashClust', 'fastcluster' ) %in% installed.packages() ) ) {
+  if ( 'flashClust' %in% installed.packages() ) {
+    lib.load.quiet( flashClust )
+  } else {
+    if ( 'fastcluster' %in% installed.packages() ) {
+      lib.load.quiet( fastcluster )
+    }
+  }
+}
+
+spec <- matrix( c( "dataset",             "d", 1, "character",
+                   "gen.new.dgram",       "g", 2, "character",
+                   "output.fname",        "o", 2, "character"
+                   ),
+                nc=4,
+                byrow=TRUE
+               )
+
+
+opt <- getopt( spec=spec )
+if ( is.null( opt$output.report.dir ) ) {
+  opt$output.report.dir <- "report"
+
+  if (! file.exists( opt$output.report.dir ) ) {
+    dir.create( opt$output.report.dir )
+  } else {
+      if ( ! file.info( 'report' )$isdir ) {
+        opt$output.report.dir <- 'heatmap.report'
+        dir.create( opt$output.report.dir )
+      }
+    }
+}
+if ( is.null( opt$output.fname ) ) { opt$output.fname <- file.path( opt$output.report.dir, paste( "data.RData", sep="." ) ) }
+if ( is.null( opt$gen.new.dgram ) ) {
+  opt$gen.new.dgram <- FALSE
+} else {
+  if ( ! opt$gen.new.dgram %in% c( "no", "yes" ) ) {
+    stop( "invalid input to gen.new.dgram param", opt$gen.new.dgram, "\n" )
+  }
+  ##  set to TRUE/FALSE
+  opt$gen.new.dgram <- ( opt$gen.new.dgram == "yes" )
+}
+ 
+
+load( opt$dataset )  ## should load the cl, treecl.res (or partcl.res) and data
+
+if ( ! exists( 'data' ) ) stop( "No data object in the rdata file provided for", opt$output.format, "format!!\n" )
+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" )
+
+## the rest of this is for the remaining output formats
+##  pre-set the cluster results for rows & cols to NULL
+direction <- NULL
+if ( exists( 'treecl.res' ) ) {
+  cl.res <- treecl.res
+  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
+} else {
+  if ( exists( 'partcl.res' ) ) {
+    cl.res <- partcl.res
+  }
+  else {
+    stop( 'could not find a valid cluster result to use for primary direction\n' )
+  }
+}
+
+if ( all( names( cl ) %in% rownames( data ) ) ) {
+  direction <- "rows"
+} else if ( all( names( cl ) %in% colnames( data ) ) ) {
+  direction <- "cols"
+  data <- t( data )
+} else {
+  stop( "Specified cluster result does not come from this data set\n" )
+}
+
+
+centroids <- NULL
+cl <- sort( cl )
+if ( inherits( cl.res, "kmeans" ) ) {
+  ## already comes pre-calculated for us!!
+  centroids <- cl.res$centers
+} else {
+  data <- data[ names( cl ), ]
+  cl.list <- unique( cl )
+  cl.list <- lapply( cl.list, function(i) cl[ cl %in% i ] )
+  centroids <- sapply( cl.list,
+                       function(x) {
+                         return( apply( data[ names(x), , drop=F ], 2, mean, na.rm=T ) )
+                       }
+                    )
+  centroids <- t( centroids )  ## get them back to the same number of columns that data has now
+}
+
+data <- centroids
+rownames( data ) <- sapply( 1:max( cl ), function(i) sprintf( "cluster-%02d", i ) )
+
+if ( opt$gen.new.dgram ) {
+  distance <- 'euclidean'
+  if ( inherits( cl.res, 'hclust' ) ) {
+    distance <- cl.res$dist.method
+  }
+  amap.distance <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
+                      "pearson", "abspearson", "correlation", "abscorrelation", "spearman", "kendall" )
+  names( amap.distance ) <- c( "euclidean", "maximum", "manhattan", "canberra", "binary",
+                               "cosine", "abscosine", "pearson", "abspearson", "spearman", "kendall" )
+
+  if ( ! distance %in% names( amap.distance ) ) stop("unsupported distance.")
+  dist.mat <- Dist( data, method=as.character( amap.distance[ distance ] ) )
+  treecl.res <- hclust( dist.mat )
+  cl <- cutree( treecl.res, nrow(data) )
+}
+
+if ( direction == "cols" ) {
+  data <- t( data )
+}
+
+save( file=opt$output.fname, treecl.res, cl, data )