diff cluster.tools/dichotomize.sample.clusters.R @ 7:2efa1a284546 draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cluster.tools/dichotomize.sample.clusters.R	Mon Mar 04 04:11:28 2013 -0500
@@ -0,0 +1,121 @@
+#!/usr/bin/env Rscript
+argspec <- c("tab.2.cdt.R converts a data matrix to cdt format
+
+        Usage: 
+                tab.2.cdt.R -d <data.file> 
+        Optional:
+                            -o <output_file>
+                \n\n")
+args <- commandArgs(TRUE)
+if ( length( args ) == 1 && args =="--help") { 
+  write(argspec, stderr())
+  q();
+}
+
+## some helper fn's
+write.2.tab <- function( mat,
+                         fname ) {
+  mat <- rbind( colnames( mat ), mat )
+  mat <- cbind( c( "ID", rownames( mat )[-1] ),
+                      mat )
+  write.table( mat, fname, sep="\t", row.names=FALSE, col.names=FALSE, quote=FALSE )
+}
+
+lib.load.quiet <- function( package ) {
+   package <- as.character(substitute(package))
+   suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
+}
+
+lib.load.quiet( getopt )
+lib.load.quiet( ctc )
+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",
+                   "num.k",               "k", 2, "character",
+                   "output.fname",        "o", 2, "character"
+                   ),
+                nc=4,
+                byrow=TRUE
+               )
+
+
+opt <- getopt( spec=spec )
+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$num.k ) ) {
+  opt$num.k <- -1
+} else {
+  num.k <- as.integer( eval( parse( text=paste( "c(",gsub( "-", ":", gsub( ", |,", ",", opt$num.k ) ), ")" ) ) ) )
+  num.k <- num.k[ ! is.na( num.k ) ]
+  if ( length( opt$num.k ) == 0 ) stop( 'invalid input for k_range specified:', opt$num.k, "\n" )
+
+  num.k <- num.k[ ! num.k %in% 1 ]  # strip out a k==1 since that doesn't make any sense
+  opt$num.k <- num.k; rm( num.k )
+}
+
+load( opt$dataset )
+## if this is a clustering result w/cluster assignments ('raw' CCPLUS does not)
+if ( exists( 'cl' ) ) {
+  k <- max( as.numeric( cl ) )
+  cl <- matrix( cl, nc=1, dimnames=list( names(cl), k ) )
+  if ( (length(opt$num.k)==1) && (opt$num.k == -1 ) ) opt$num.k <- k
+  
+  ## if this is a one-off to produce a phenotype for the number of clusters that the user originally proposed
+  if ( !opt$num.k[1] %in% c( -1, k ) ) {
+
+    if ( exists( 'partcl.res' ) || exists( 'select.result' ) ) {
+    
+      if ( exists( 'partcl.res' ) ) {
+        warning( 'The k_range value(s) specified are:',
+                 opt$num.k,
+                 "however k_range vals can not specify alternate k values for partition clusters. Using the K value that corresponds to this result instead\n" )
+      } else {
+        warning( 'The k_range value(s) specified are:',
+                 opt$num.k,
+                 "however k_range vals can not specify alternate k values for specific cluster results from CCPLUS (i.e. those from the Select K or Extract tools). To get alternate K values, re-run the dichotomizer on the 'raw' CCPLUS results. Using the K value that corresponds to this result instead\n" )
+      }
+      
+      opt$num.k <- k
+      cl <- matrix( cl, nr=1, dimnames=list( k, names(cl) ) )
+    } else {
+      ## handle if this is a hclust result
+
+      opt$num.k <- opt$num.k[ opt$num.k < length( cl ) ]
+      cl.samps <- rownames( cl )
+      cl <- sapply( opt$num.k, function(i) cutree( treecl.res, i )[ cl.samps ] )
+      colnames( cl ) <- opt$num.k
+    }
+  }
+} else if ( exists( 'results' ) ) {
+  ## handle if this is a ccplus-raw result
+  opt$num.k <- opt$num.k[ opt$num.k <= length( results ) ]
+  cl <- sapply( results[ opt$num.k ], '[[', 'consensusClass' )
+  colnames( cl ) <- opt$num.k
+}
+
+pheno.mat <- lapply( 1:ncol(cl),
+                     function(i) {
+                       x <- cl[,i]
+                       cls <- ks <- sort( unique(x) )
+                       cls <- sapply( cls, function(y) as.numeric( x %in% y ) )
+                       colnames(cls) <- paste( "CLeq", ks, sep="" )
+                       rownames(cls) <- names(x)
+                       return(cls)
+                     }
+                    )
+names( pheno.mat ) <- opt$num.k
+
+final.mat <- matrix( NA, nc=0, nrow=nrow(cl), dimnames=list( names(cl), NULL ) )
+for ( i in names( pheno.mat ) ) {
+  colnames( pheno.mat[[i]] ) <- paste( "Keq", i, "_", colnames( pheno.mat[[i]] ), sep="" )
+  final.mat <- cbind( final.mat, pheno.mat[[i]] )
+}
+
+write.2.tab( final.mat, opt$output.fname )