view cluster.tools/hclust.R @ 9:a3c03541fe6f draft default tip

Uploaded

#!/usr/bin/env Rscript

argspec <- c("hclust.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 )
##  if any of the faster clustering methods are available on this system, load them
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( "data.fname",         "d", 1, "character",
                   "num.k",              "k", 1, "integer",
                   "distance.metric",    "m", 2, "character", 
                   "dist.obj",           "D", 2, "logical",
                   "direction",          "n", 2, "character",
                   "linkage",            "l", 2, "character",
                   "output.name",        "o", 2, "character"
                   ),
                nc=4,
                byrow=TRUE
               )

opt <- getopt( spec=spec )

data <- as.matrix( read.delim( opt$data.fname, header=T, row.names=1 , check.names=FALSE ) )
if ( is.null( opt$distance.metric ) ) { opt$distance.metric <- "euclidean" }
if ( is.null( opt$dist.obj ) ) { opt$dist.obj <- FALSE }
if ( is.null( opt$direction ) ) { opt$direction <- "cols"  }
if ( is.null( opt$linkage ) ) { opt$linkage <- "average" }
if ( is.null( opt$output.name ) ) { opt$output.name <- "hclust.result.rda" }
if ( is.null( opt$num.k ) || ( opt$num.k == -1 )) {
  if ( opt$direction == 'cols' ) {
    opt$num.k <- 5
  } else if ( opt$direction == 'rows' ) {
    opt$num.k <- nrow( data ) / 30  ## we use an estimated average size of gene clusters to be 30
    if ( opt$num.k > 1000 ) {
      opt$num.k <- ( opt$num.k %/% 10 ) * 10
    } else {
      opt$num.k <- ( opt$num.k %/% 5 ) * 5
    }
  }
}

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 )
}

if ( opt$dist.obj ) {
  dist.mat <- as.dist( data )  
} else {
  ##  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" )

  if ( ! opt$distance.metric %in% names( amap.distance ) ) stop("unsupported distance.")
  dist.mat <- Dist( data, method=as.character( amap.distance[ opt$distance.metric ] ) )
  attr( dist.mat, "method" ) <- opt$distance.metric
}

##  now, do the clustering
treecl.res <- hclust( dist.mat, method=opt$linkage )
cutree.res <- cutree( treecl.res, k=opt$num.k )
##cl <- cbind( names( cutree.res ), as.numeric( cutree.res ) )
##colnames( cl ) <- c( "ID", "class" )

if ( opt$direction == "cols" ) {
  ## need to re-transpose the data back to it's original dimensionality
  data <- t( data )
}

cl <- cutree.res
save( file=opt$output.name, treecl.res, cl, data )