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

Uploaded

7:2efa1a284546

#!/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 )