view cluster.tools/cluster.2.centroid.R @ 2:b442996b66ae draft

Uploaded

#!/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 ( is.null( opt$output.fname ) ) { opt$output.fname <- file.path( opt$output.report.dir, paste( "data", opt$output.format, 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), ], 2, mean, na.rm=T ) )
                       }
                    )
  centroids <- t( centroids )  ## get them back to the same number of columns that data has now
}

data <- centroids
colnames( 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
  }
  dmat <- Dist( data, distance )
  treecl.res <- hclust( dmat )
  cl <- cutree( treecl.res )
}

if ( direction == "cols" ) {
  data <- t( data )
}

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