view cluster.tools/heatmap.from.cluster.result.R @ 9:a3c03541fe6f draft default tip

Uploaded

#!/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( gplots )
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",
                   "second.dir",          "s", 2, "character",
                   "dataset2",            "D", 2, "character",
                   "reverse.rows",        "r", 2, "character",
                   "image.format",        "i", 2, "character",
                   "plot.kms",            "k", 2, "character",
                   "output.fname",        "o", 2, "character",
                   "output.report.html",  "h", 2, "character",
                   "output.report.dir",   "p", 2, "character",
                   "output.treeview",     "t", 2, "character",
                   "survival.script",     "z", 2, "character",
                   "cluster.fname",       "C", 2, "character",
                   "survival.fname",      "S", 2, "character",
                   "survival.image",      "I", 2, "character",
                   "survival.mode",       "M", 2, "character",
                   "survival.title",      "T", 2, "character"
                   ),
                nc=4,
                byrow=TRUE
               )


opt <- getopt( spec=spec )
if ( is.null( opt$image.format ) ){
  opt$image.format <- "png"
} else {
  if ( ! opt$image.format %in% c( "pdf", "png" ) ) stop( 'invalid image format specified\n' )
}
if ( is.null( opt$output.report.dir ) ) { opt$output.report.dir <- "report" }
if ( is.null( opt$output.report.html ) ) {
  opt$out.dir <- 'report'
  if (! file.exists( opt$out.dir ) ) {
    dir.create( opt$out.dir )
  } else {
      if ( ! file.info( 'report' )$isdir ) {
        opt$out.dir <- 'heatmap.report'
        dir.create( opt$out.dir )
      }
    }

  if ( opt$image.format == "pdf" ) opt$output.report.html <- file.path( opt$out.dir ,"heatmap.pdf" )
  if ( opt$image.format == "png" ) opt$output.report.html <- file.path( opt$out.dir ,"index.html" )
}
if ( is.null( opt$plot.kms ) ) {
  opt$plot.kms <- FALSE
}  else {
  if ( ! opt$plot.kms %in% c( "no", "yes" ) ) {
    stop( "invalid input to plot.kms param", opt$plot.kms, "\n" )
  }
  ##  set to TRUE/FALSE
  opt$plot.kms <- ( opt$plot.kms == "yes" )
  if ( opt$plot.kms ) {
    opt$cluster.fname <- opt$dataset
    if ( is.null( opt$survival.fname ) || ( !file.exists( opt$survival.fname ) ) ) stop( 'must provide a valid file w/clinical data\n' )
    if ( is.null( opt$survival.script ) || ( !file.exists( opt$survival.script ) ) ) stop( 'must provide a valid path to the gen.survival.curves.R file\n' )
    if ( is.null(opt$mode ) ) {
      opt$mode <- "all"
    } else {
      if ( ! opt$mode %in% c( 'all', 'one', 'both' ) ) {
        stop( "invalid mode specified,' -m", opt$mode, "'.  must be either {all, one, both}\n" )
      }
    }
    if ( is.null( opt$title ) ) {
      opt$title <- opt$cluster.fname
      opt$title <- strsplit( opt$title, "\\/" )[[1]]
      opt$title <- opt$title[ length( opt$title ) ]
    }
  }
}
if ( is.null( opt$output.treeview ) ) {
  opt$output.treeview <- FALSE
} else {
  if ( ! opt$output.treeview %in% c( "no", "yes" ) ) {
    stop( "invalid input to output.treeview param", opt$output.treeview, "\n" )
  }
  ##  set to TRUE/FALSE
  opt$output.treeview <- ( opt$output.treeview == "yes" )
}
if ( is.null( opt$reverse.rows ) ) {
  opt$reverse.rows <- TRUE
}  else {
  if ( ! opt$reverse.rows %in% c( "no", "yes" ) ) {
    stop( "invalid input to reverse.rows param", opt$reverse.rows, "\n" )
  }

  ##  set to TRUE/FALSE
  opt$reverse.rows <- ( opt$reverse.rows == "yes" )
}

if ( is.null( opt$second.dir ) ) { opt$second.dir <- "no" }
if ( is.null( opt$second.dir ) && is.null( opt$dataset2 ) ) stop( "must specify an rdata file to load if a previous result is to be used to cluster the 2nd direction\n" )
if ( ( opt$image.format == "png" ) || opt$output.treeview ) {
  if ( !file.exists( opt$output.report.dir ) ){
    dir.create(opt$output.report.dir, recursive=T)
  }
}



load( opt$dataset )  ## should load the cl, treecl.res (or partcl.res) and data
##  pre-set the cluster results for rows & cols to NULL
hr <- hr.cl <- hc <- hc.cl <- row.ddr <- col.ddr <- 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 ( opt$second.dir %in% c( "yes", "prev" ) ) {

  if ( opt$second.dir == "yes" ) {
    if ( is.null( hr ) ) {
      hr <- hclust( dist( data ) )
    } else if ( is.null( hc ) ) {
      hc <- hclust( dist( t( data ) ) )
    } 
  } else {  ## opt$second.dir == "prev"

    ## 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 ( any( is.na( orig.data ) ) ) {
      nas <- which( is.na( orig.data ) )
      num.nas <- length( nas )
      ## 1st, chk the NAs
      if ( sum( which( is.na( data ) ) != nas ) == num.nas ) stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
      if ( ( sum( orig.data == data, na.rm=T )+num.nas ) != length( orig.data ) ) {
        stop( "incompatible cluster results in 2nd results file - matrices contain diff contents\n" )
      }
      
    } else {
      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 ( ! is.null( hr ) ) {
  row.ddr <- as.dendrogram( hr )
  ##  need this to make sure that the heatmap is oriented the same way as it is in TreeView
  if ( opt$reverse.rows ) row.ddr <- rev( row.ddr )  
}

if ( ! is.null( hc ) ) {
  col.ddr <- as.dendrogram( hc )
}
  

hmcols<-colorRampPalette(c("blue","white","red"))(256)

if ( ( ! is.null( row.ddr ) ) && ( ! is.null( col.ddr ) ) ) {
  dend.param <- "both"
} else {
  dend.param <- "none"
  if ( ! is.null( row.ddr ) ) dend.param <- "row"
  if ( ! is.null( col.ddr ) ) dend.param <- "column"
}


param.list <- list( x=data,
                    Rowv=row.ddr,
                    Colv=col.ddr,
                    dendrogram=dend.param,
                    trace="none",
                    col=hmcols,
                    symbreaks=TRUE,
                    scale="none",
                    labRow="",
                    labCol="",
                    na.color='grey' ) #,
                    ##key=FALSE )

if ( ! is.null( hr.cl ) ) {
  hrcols <- rainbow( max( as.numeric( hr.cl ) ) )
  names( hrcols ) <- sort( unique( as.numeric( hr.cl ) ) )
  rowColLabs <- hrcols[ as.character( as.numeric( hr.cl ) ) ]
  param.list <- c( param.list, list( RowSideColors=rowColLabs ) )
}
if ( ! is.null( hc.cl ) ) {
  hccols <- rainbow( max( as.numeric( hc.cl ) ) )
  names( hccols ) <- sort( unique( as.numeric( hc.cl ) ) )
  colColLabs <- hccols[ as.character( as.numeric( hc.cl ) ) ]
  param.list <- c( param.list, list( ColSideColors=colColLabs ) )
}


if ( opt$image.format == 'png' ) {
  png.fname <- file.path( opt$output.report.dir, "cluster.heatmap.png")
  plot.dev <- png( png.fname,
                   width=8.5,
                   height=11,
                   units='in',
                   res=72 )
} else {
  pdf.fname <- opt$output.report.html
  pdf( opt$output.report.html,
       paper="letter" )
}

do.call( "heatmap.2", param.list )

dev.off() ## close the previous device

if ( opt$plot.kms ) {

  cmd.string <- opt$survival.script

  ## get the consensusClass file that's associated with the k.select
  cmd.string <- paste( cmd.string, "-C", opt$dataset )
  cmd.string <- paste( cmd.string, "-S", opt$survival.fname )
  cmd.string <- paste( cmd.string, "-M", opt$survival.mode )
  if ( ! is.null( opt$suvival.title ) )
    cmd.string <- paste( cmd.string, "-T", opt$suvival.title )

  ##  only call kms if we're the image is png
  if ( opt$image.format=="png" ) {
    png.fname <- file.path( opt$output.report.dir, "kaplan.meier.survival.png")
    cmd.string <- paste( cmd.string, "-I", "png" ) 
    cmd.string <- paste( cmd.string, "-O", png.fname )
    system( cmd.string )
  }
}


if ( opt$image.format == 'png' ) {
  pngs = list.files(path=opt$output.report.dir, patt="png")
  html.out <- paste( "<html>",
                     paste( paste( paste( "<div><img src=\'", pngs, sep="" ), "\'/></div>", sep="" ), collapse=""),
                     "</html>" )
  cat( html.out, file=opt$output.report.html )
}


if ( opt$output.treeview ) {
  treeview.fname.stem <- file.path( opt$output.report.dir, "cluster.heatmap")
  fnames <- character()
  if ( ! is.null( hr ) ) {
    fname <- paste( treeview.fname.stem, ".gtr", sep="" )
    r2gtr( hr, file=fname )
    fnames <- c( fnames, fname )
  } else {
    hr <- list( order=1:nrow( data ) )
  }
  if ( ! is.null( hc ) ) {
    fname <- paste( treeview.fname.stem, ".atr", sep="" )
    r2atr( hc, file=fname )
    fnames <- c( fnames, fname )
  } else {
    hc <- list( order=1:ncol( data ) )
  }


  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 )
}