diff cluster.tools/rdata.2.out.R @ 0:0decf3fd54bc draft

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/cluster.tools/rdata.2.out.R	Thu Feb 28 01:45:39 2013 -0500
@@ -0,0 +1,224 @@
+#!/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();
+}
+
+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",
+                   "dataset2",            "D", 2, "character",
+                   "output.format",       "f", 2, "character",
+                   "output.report.dir",   "p", 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$output.format ) ) { opt$output.format <- "cdt" }
+
+
+load( opt$dataset )  ## should load the cl, treecl.res (or partcl.res) and data
+
+if ( opt$output.format %in% c( "cls-only", "newick" ) ) {
+  if ( opt$output.format == "cls-only" ) {
+
+    cl <- cbind( names( cl ), as.numeric( cl ) )
+    colnames( cl ) <- c( "ID", "Class" )
+
+    opt$output.fname <- gsub( "cls-only$", "tab", opt$output.fname )
+    write.table( cl, opt$output.fname, sep="\t", quote=FALSE, row.names=FALSE, col.names=FALSE )
+  } else {
+    ##if ( opt$output.format == "newick" ) {
+
+    if ( ! exists( "treecl.res" ) ) stop( "no HAC result found in results file proved - necessary to generate a Newick formated file.\n" )
+    write( hc2Newick( treecl.res ), opt$output.fname )
+  }
+} else {
+  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
+  hr <- hr.cl <- hc <- hc.cl <- NULL
+  if ( exists( '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
+    if ( all( names( cl ) %in% rownames( data ) ) ) {
+      hr <- treecl.res
+      hr.cl <- cl
+    } else if ( all( names( cl ) %in% colnames( data ) ) ) {
+      hc <- treecl.res
+      hc.cl <- cl
+    } else {
+      stop( "Specified cluster result does not come from this data set\n" )
+    }
+
+  } else {
+    if ( exists( 'partcl.res' ) ) {
+      if ( all( names( cl ) %in% rownames( data ) ) ) {
+        hr <- NA
+        hr.cl <- cl
+        orig.data <- data
+        data <- data[ names( cl ), ]  ## partcl.res should now be sorted in order of cluster
+      } else if ( all( names( cl ) %in% colnames( data ) ) ) {
+        hc <- NA
+        hc.cl <- cl
+        orig.data <- data
+        data <- data[ , names( cl ) ]  ## partcl.res should now be sorted in order of cluster
+      } else {
+        stop( "Specified cluster result does not come from this data set\n" )
+      }
+    }
+    else {
+      stop( 'could not find a valid cluster result to use for primary direction\n' )
+    }
+  }
+  
+
+  if ( ! is.null( opt$dataset2 ) ) {
+    
+    ## prep for loading new cluster result
+    if ( ! exists( 'orig.data' ) ) orig.data <- data
+    if ( exists( "treecl.res" ) ) {
+      rm( treecl.res )
+    } else if ( exists( "partcl.res" ) ) {
+      rm( partcl.res )
+    } else stop( "no primary clustering found when generating the 2nd\n" )
+    rm( cl, data )
+    
+  
+    load( opt$dataset2 ) ## this should bring in the cl obj for the 2nd direction
+
+    ## check the data 1st
+    if ( length( orig.data ) != length( data ) ) stop( "incompatible cluster results in 2nd results file - matrices are diff lengths\n" )
+    if ( nrow( orig.data ) != nrow( data ) ) stop( "incompatible cluster results in 2nd results file - matrices have diff dimensions\n" )
+    if ( sum( orig.data == data ) != length( orig.data ) )  stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
+    ## looks like data is the same, so drop a copy and start chugging
+    rm( orig.data ); gc()
+    
+    if ( exists( '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
+      
+      if ( is.null( hr ) ) {
+        if ( all( rownames( cl ) %in% rownames( data ) ) ) {
+          hr <- treecl.res
+          hr.cl <- cl
+        } else {
+          stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
+        }
+      } else if ( is.null( hc ) ) {
+        if ( all( rownames( cl ) %in% colnames( data ) ) ) {
+          hc <- treecl.res
+          hc.cl <- cl
+        } else {
+          stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
+        }
+      } else {
+        stop( "should never get here\n" )
+      }
+    } else if ( exists( 'partcl.res' ) ) {
+      if ( is.null( hr ) ) {
+        if ( all( names( cl ) %in% rownames( data ) ) ) {
+          hr <- NA
+          hr.cl <- cl
+          data <- data[ names( cl ), ]  ## partcl.res should now be sorted in order of cluster
+        } else {
+          stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (rows in this case)\n" )
+        }
+      } else if ( is.null( hc ) ) {
+        if ( all( names( cl ) %in% colnames( data ) ) ) {
+          hc <- NA
+          hc.cl <- cl
+          data <- data[ , names( cl ) ]  ## partcl.res should now be sorted in order of cluster
+        } else {
+          stop( "results file for 2nd direction doesn't contain cluster for 2ndary direction (genes in this case)\n" )
+        }
+      } else {
+        stop( "should never get here\n" )
+      }      
+    }
+  }
+  
+  ## Now, re-set hc & nr to NULL if they were set to NA
+  ## we used NA to signify that they were set by kmeans/pam, but now, we need to reset them
+  ## for the following lines (that generate the dendrograms (if there was an hclust result)
+  if ( ( !is.null( hr ) ) && is.na( hr ) ) hr <- NULL
+  if ( ( !is.null( hc ) ) && is.na( hc ) ) hc <- NULL
+  
+  if ( ! exists( 'data' ) ) stop( "No data object in the rdata file provided!!\n" )
+  
+  if ( is.null( hc ) ) hc <- list( order=1:ncol( data ) )
+  if ( is.null( hr ) ) hr <- list( order=1:nrow( data ) )
+
+  if ( opt$output.format == "tabular" ) {
+    write.table( data[ hr$order, hc$order ], opt$output.fname, quote=FALSE, sep="\t", col.names=NA )
+  } else if ( opt$output.format == "cdt" ) {
+    if ( !file.exists( opt$output.report.dir ) ){
+      dir.create(opt$output.report.dir, recursive=T)
+    }
+    
+    treeview.fname.stem <- file.path( opt$output.report.dir, "cluster.heatmap")
+    fnames <- character()
+    if ( inherits( hr, "hclust" ) ) {
+      fname <- paste( treeview.fname.stem, ".gtr", sep="" )
+      ## we manually specify a 'stub' distance b/c o/w it'll try using the attr(hr,"method")
+      ##  and the r2gtr fn's get grumpy if the distance was anything starting with a 'p'
+      r2gtr( hr, file=fname, distance='stub' )  
+      fnames <- c( fnames, fname )
+    }
+    if ( inherits( hc, "hclust" ) ) {
+      fname <- paste( treeview.fname.stem, ".atr", sep="" )
+      r2atr( hc, file=fname, distance='stub' )
+      fnames <- c( fnames, fname )
+    }
+
+    fname <- paste( treeview.fname.stem, ".cdt", sep="" )
+    r2cdt( hr, hc, data, file=fname )
+    fnames <- c( fnames, fname )
+
+    ## jtv file now
+    jtv.str <- '<DocumentConfig><UrlExtractor/><ArrayUrlExtractor/><Views><View type="Dendrogram" dock="1"><ColorExtractor contrast="2.0"><ColorSet zero="#FFFFFF" down="#0000FF"/></ColorExtractor><ArrayDrawer/><GlobalXMap current="Fill"><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></GlobalXMap><GlobalYMap current="Fill"><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></GlobalYMap><ZoomXMap><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></ZoomXMap><ZoomYMap><FixedMap type="Fixed"/><FillMap type="Fill"/><NullMap type="Null"/></ZoomYMap><TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView><TextView face="Monospaced" size="14"><GeneSummary/></TextView></TextView><ArrayNameView face="Monospaced" size="14"><ArraySummary included="0"/></ArrayNameView><AtrSummary/><GtrSummary/></View></Views></DocumentConfig>'
+    fname <- paste( treeview.fname.stem, ".jtv", sep="" )
+    cat( jtv.str, file=fname )
+    fnames <- c( fnames, fname )
+
+    cmd <- paste( "tar -zcf",
+                  opt$output.fname,
+                  paste( "--directory=", opt$output.report.dir, sep="" ),
+                  paste( basename( fnames ), collapse=" " ) )
+    system( cmd )
+  }
+}
+
+
+
+