Mercurial > repos > peter-waltman > ucsc_cluster_tools
diff cluster.tools/cluster.2.centroid.R @ 2:b442996b66ae draft
Uploaded
author | peter-waltman |
---|---|
date | Wed, 27 Feb 2013 20:17:04 -0500 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/cluster.tools/cluster.2.centroid.R Wed Feb 27 20:17:04 2013 -0500 @@ -0,0 +1,119 @@ +#!/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 )