diff cluster.tools/partition.R @ 2:b442996b66ae draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cluster.tools/partition.R	Wed Feb 27 20:17:04 2013 -0500
@@ -0,0 +1,168 @@
+#!/usr/bin/env Rscript
+
+argspec <- c("partition.R help TBD
+                \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 )
+lib.load.quiet( cluster )
+
+##  we're going to use the amap Dist function, but they misname their correlation
+##  functions, so re-name them correctly
+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" )
+
+spec <- matrix( c( "data.fname",         "d", 1, "character",
+                   "algorithm",        "a", 2, "character",
+                   "distance.metric",    "m", 2, "character", ## must be one supported by R's dist function
+                   "dist.obj",           "D", 2, "logical",
+                   "direction",          "n", 2, "character",
+                   "num.k",              "k", 2, "integer",
+                   "output.name",        "o", 2, "character"
+                   ),
+                nc=4,
+                byrow=TRUE
+               )
+
+opt <- getopt( spec=spec )
+
+if ( is.null( opt$distance.metric ) ) { opt$distance.metric <- "euclidean" }
+if ( is.null( opt$algorithm ) ) { opt$algorithm <- "km" }
+if ( is.null( opt$dist.obj ) ) { opt$dist.obj <- FALSE }
+if ( is.null( opt$direction ) ) { opt$direction <- "cols"  }
+if ( is.null( opt$num.k ) ) { opt$num.k <- 10 }
+if ( is.null( opt$output.name ) ) { opt$output.name <- "partition.result" }
+
+data <- as.matrix( read.delim( opt$data.fname, header=T, row.names=1 , check.names=FALSE ) )
+
+if ( opt$direction == "cols" ) {
+  ## need to transpose b/c both kmeans & pam cluster the rows
+  ## this shouldn't have an effect upon a distance matrix
+  data <- t( data )
+}
+if ( opt$num.k > nrow( data ) ) {
+  err.msg <- paste( "K specified is greater than the number of elements (", opt$direction, ") in data matrix to be clustereed\n", sep="" )
+  stop( err.msg )
+}
+
+mat.2.b.clustered <- data
+if ( opt$dist.obj ) {
+  ## To be updated
+
+  mat.2.b.clustered <- as.dist( data )
+
+  if ( opt$algorithm=="km" ) {
+    ##clusterAlg is kmeans
+    if (verbose) warning()
+  }
+} else {
+  ## this is a data matrix -- we always generate a dist.mat object (b/c we need it
+  ##  in case this result is used with a heatmap
+  
+  ## PAM clustering
+  if ( opt$algorithm != "km" ) {
+
+    if ( ! opt$distance.metric %in% names( amap.distance ) ) stop("unsupported distance.")
+    mat.2.b.clustered <- Dist( data, method=as.character( amap.distance[ opt$distance.metric ] ) )
+    attr( mat.2.b.clustered, "method" ) <- opt$distance.metric
+
+  } else {
+    mat.2.b.clustered <- data
+  }
+}
+
+## now run the clustering
+partcl.res <- cl <- NA
+if ( opt$algorithm=="pam" ) {
+  partcl.res <- pam( x=mat.2.b.clustered,
+                     k=opt$num.k,
+                     metric=( ifelse( inherits(data, "dist" ), 
+                                      attr( data, "method" ),  ## this is ok if data is a dist object (b/c pam will ignore it)
+                                      opt$distance.metric) ) ) ##,
+                     ##cluster.only=TRUE )
+  cl <- partcl.res$clustering
+
+  if ( is.character( partcl.res$medoids ) ) {
+    medoids <- data[ partcl.res$medoids, ]
+  } else {
+    ##partcl.res$medoids is a matrix -- we shouldn't get this (only if mat.2.b.clustered is a data matrix)
+    medoids <- partcl.res$medoids
+  }
+  med.names <- rownames( medoids )
+  med.hc <- hclust( as.dist( as.matrix( mat.2.b.clustered )[ med.names, med.names ] ) )
+  med.cls <- as.numeric( cl[ med.names[ med.hc$order ] ] )
+
+  cl.list <- lapply( med.cls, function(i) names( cl[ cl %in% i ] ) )
+  names( cl.list ) <- med.cls 
+
+  cl.list <- lapply( cl.list,
+                     function( elts ) {
+                       if ( length( elts ) == 1 ) {
+                         retval <- 1
+                         names( retval ) <- elts
+                       } else {
+                         subdist <- as.dist( as.matrix( mat.2.b.clustered )[ elts, elts ] )
+                         sub.hc <- hclust( subdist )
+                         retval <- sub.hc$order
+                         names( retval ) <- sub.hc$labels
+                         retval <- sort( retval )
+                       }
+                       return( retval )
+                     }
+                    )
+  
+  fnl.ord <- as.character( unlist( lapply( cl.list, names ) ) )
+  cl <- cl[ fnl.ord ]
+} else {
+  partcl.res <- kmeans( x=mat.2.b.clustered,
+                        centers=opt$num.k )
+  cl <- partcl.res$cluster
+  centroids <- partcl.res$centers
+  cent.hc <- hclust( Dist( centroids, method=as.character( amap.distance[ opt$distance.metric ] ) ) )
+  cent.cls <- as.numeric( cent.hc$labels[ cent.hc$order ] )
+
+  cl.list <- lapply( cent.cls, function(i) names( cl[ cl %in% i ] ) )
+  names( cl.list ) <- cent.cls 
+
+  cl.list <- lapply( cl.list,
+                     function( elts ) {
+                       if ( length( elts ) == 1 ) {
+                         retval <- 1
+                         names( retval ) <- elts
+                       } else {
+                         if ( all( elts %in% colnames( mat.2.b.clustered ) ) ) {
+                           submat <- mat.2.b.clustered[ elts, elts ]
+                         } else {
+                           submat <- mat.2.b.clustered[ elts, ]
+                         }
+                         subdist <- Dist( submat, method=as.character( amap.distance[ opt$distance.metric ] ) )
+                         sub.hc <- hclust( subdist )
+                         retval <- sub.hc$order
+                         names( retval ) <- sub.hc$labels
+                         retval <- sort( retval )
+                       }
+                       return( retval )
+                     }
+                    )
+  
+  fnl.ord <- as.character( unlist( lapply( cl.list, names ) ) )
+  cl <- cl[ fnl.ord ]  
+}
+
+if ( opt$direction == "cols" ) {
+  ## need to transpose back
+  data <- t( data )
+}
+save( file=opt$output.name, partcl.res, cl, data )
+