view cluster.tools/gen.survival.curves.R @ 0:0decf3fd54bc draft

Uploaded

#!/usr/bin/env Rscript
## 
## Calculates the log-rank test for a given clustering, in the output format from ConsensusClusterPlus
##
## Input (required):
##    - consensus cluster file (consensusClass.csv file)
##    - survival data
## Input (optional):
## Output: a KM plot, with the most significant p-value.  Output to stdout can be captured by re-direction
##
## Uses: survival library
## Date: August 21, 2012
## Author: Peter Waltman
##

##usage, options and doc goes here
argspec <- c("gen.survival.curves.R takes a clustering from ConsensusClusterPlus and clinical survival data
and generates a KM-plot, along with the log-rank p-values

        Usage: 
                gen.survival.curves.R -c <cluster.file> -s <clinical.file> 
        Options:
                -o <output file> (postscript)
                -m <mode>        (all, one, both)
                                  \"all\" - perform all-vs-all log-rank test
                                  \"one\" - perform one-vs-others log-rank test (returns min)
                                  \"both\" - perform both \"all\" and \"one\" tests
                -t <title>
                -p <pval.only>  ( only return the p-value for this given mode - no plotting at all (to screen or postscript))
                -v <verbose>
              ")

lib.load.quiet <- function( package ) {
   package <- as.character(substitute(package))
   suppressPackageStartupMessages( do.call( "library", list( package=package ) ) )
}
lib.load.quiet(getopt)
lib.load.quiet( survival )


args <- commandArgs(TRUE)
if ( length( args ) == 1 && args =="--help") { 
  write(argspec, stderr())
  q();
}

spec <- matrix( c( "cluster.fname",  "C", 1, "character",
                   "survival.fname", "S", 1, "character",  
                   "mode",           "M", 2, "character",
                   "title",          "T", 2, "character",
                   "myplots.rda",    "R", 2, "character",
                   "image.format",   "I", 2, "character",
                   "output.fname",   "O", 2, "character",
                   "pval.only",      "P", 0, "logical",
                   "verbose",        "V", 0, "logical"
                  ),
                ncol=4,
                byrow=TRUE
               )
opt <- getopt( spec=spec )


#set some reasonable defaults for the options that are needed,
#but were not specified.
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$image.format ) ){
  opt$image.format <- "png"
} else {
  if ( ! opt$image.format %in% c( "pdf", "png", "none" ) ) stop( 'invalid image format specified\n' )
}
if ( is.null(opt$output.fname ) ) { opt$output.fname <- paste( opt$mode, "survival.curve", opt$image.format, sep="." ) }
if ( is.null(opt$cluster.header ) ) { opt$cluster.header = FALSE }
if ( is.null(opt$pval.only ) ) { opt$pval.only = FALSE }
if ( is.null(opt$verbose ) ) { opt$verbose = FALSE }

##print some progress messages to stderr, if requested.
if ( opt$verbose ) { write("writing...",stderr()); }

load( opt$cluster.fname )
cluster.data <- cbind( names( cl ), as.numeric( cl ) )
colnames( cluster.data ) <- c( "id", "group_num" )
rownames( cluster.data ) <- names( cl )

survival.data <- read.delim( opt$survival.fname, as.is=TRUE, row.names=1 )
survival.data <- cbind( rownames( survival.data ), survival.data ) ## add in the ids, so we can merge on them
if ( length( colnames( survival.data ) ) == 3 ) {
  ## we have to left-shift the current colanmes to drop the 1st one
  ##  b/c cbind will add one for the column we just  added
  colnames( survival.data ) <- c( "id", colnames( survival.data )[-1] )  
}
if ( length( colnames( survival.data ) ) == 2 ) {
  ## added just in case there's a change to cbind as R is prone to doing
  colnames( survival.data ) <- c( "id", colnames( survival.data ) )  
}
survival.data$id <- as.character( survival.data$id )


## Now, filter so we only contain the same samples
n.clust.data.samps <- nrow( cluster.data )
n.surv.data.samps <- nrow( survival.data )
if ( n.clust.data.samps > n.surv.data.samps ) {
  ovp.samples <- rownames( cluster.data )
  ovp.samples <- ovp.samples[ ovp.samples %in% survival.data$id ]
} else {
  ovp.samples <- survival.data$id
  ovp.samples <- ovp.samples[ ovp.samples %in% rownames( cluster.data ) ]
}

cluster.data <- cluster.data[ ovp.samples, , drop=FALSE]
survival.data <- survival.data[ ovp.samples, ]
survival.data <- merge( survival.data, cluster.data )


calc.all.pval <- function( survival.data ) {
  survdiff( Surv( time, status )~group_num, data=survival.data )
  surv.res <- survdiff( Surv( time, status )~group_num, data=survival.data )
  pval <- surv.res$chisq
  df <- length( surv.res$n ) - 1
  pval <- pchisq( pval, df=df, lower=F )
  return( pval )
}

calc.one.v.others.pval <- function( survival.data ) {
  grps <- sort( unique( as.numeric( survival.data$group_num ) ) )

  retval <- numeric()
  for ( g in grps ) {
    one.v.all.survival.data <- survival.data
    tmp <- as.numeric( one.v.all.survival.data$group_num )
    tmp[ ! tmp %in% g ] <- -1
    tmp[ tmp %in% g ] <- 1
    tmp[ tmp %in% -1 ] <- 2
    one.v.all.survival.data$group_num <- tmp
    surv.res <- survdiff( Surv( time, status )~group_num, data=one.v.all.survival.data )
    pval <- surv.res$chisq
    df <- length( surv.res$n ) - 1
    pval <- pchisq( pval, df=df, lower=F )
    retval <- c( retval, pval )
  }
  names( retval ) <- grps
  return( retval )
}


if ( opt$mode == "all" ) {

  pval <- calc.all.pval( survival.data )
  log.rank <- paste( "Log Rank p-value:", sprintf( "%1.2e",pval ) )
  opt$title <- paste( opt$title, log.rank, sep="\n" )
} else {
  if ( opt$mode == "one" ) {

    pvals <- calc.one.v.others.pval( survival.data )
    min.p <- min( pvals, na.rm=T )
    if ( length( min.p ) == 0 ) {
      stop( 'no valid p-value returned from the one-v-others test\n' )
    }
    cluster.num <- names( pvals )[ pvals == min.p ]
    pval <- pvals[ cluster.num ]
    log.rank <- paste( "Log Rank p-value for cluster", cluster.num,"is:", sprintf( "%1.2e",pval ) )
    opt$title <- paste( opt$title, log.rank, sep="\n" )
  } else {
    if ( opt$mode== "both" ) {
      ##  add the all-v-all p-value
      bak <- pval <- calc.all.pval( survival.data )
      log.rank <- paste( "Log Rank p-value:", sprintf( "%1.2e",pval ) )
      opt$title <- paste( opt$title, log.rank, sep="\n" )

      ## now add the one-v-all p-value
      pvals <- calc.one.v.others.pval( survival.data )
      min.p <- min( pvals, na.rm=T )
      if ( length( min.p ) == 0 ) {
        stop( 'no valid p-value returned from the one-v-others test\n' )
      }
      cluster.num <- names( pvals )[ pvals == min.p ]
      pval <- pvals[ cluster.num ]
      log.rank <- paste( "Log Rank p-value for cluster", cluster.num,"is:", sprintf( "%1.2e",pval ) )
      opt$title <- paste( opt$title, log.rank, sep="\n" )

      if ( opt$pval.only ) {
        pval <- min( c( bak, pval ), na.rm=T )
      }
    }
    else {
      stop( "invalid mode specified, mode = ", opt$mode, "\n" )
    }
  }
}

if ( opt$pval.only ) {
  cat( paste(pval, "\n", sep="" ), file=stdout() )
}


if ( ! opt$pval.only ) {
  ngrps <- length( unique( survival.data$group_num ) )
  col.map <- rainbow( ngrps )



  ##postscript( opt$output.fname, horizontal=T, paper='letter' )
  if ( opt$image.format == 'png' ) {
    plot.dev <- png( opt$output.fname,
                     width=11,
                     height=8.5,
                     units='in',
                     res=72 )
  } else if ( opt$image.format == 'pdf' ) {
    pdf( opt$output.fname,
         paper="letter" )
  } else if ( opt$image.format == 'none' ) {
    ## do nothing - this allows other scripts to call this and hopefully plot into them
    ## NOPE, this doesn't work.  see what I do with the myplots.rda file
  }
  
  plot( survfit( Surv( time, status )~group_num, data=survival.data ),
        main = opt$title,
        ##lty = 1:ngrps,
        lty=1,
        col=col.map,      
        ylab = "Probability",
        xlab = "Survival Time in Days",
       )


  ## set the legend.labels if they're still not set yet
  if( ! exists( "legend.labels" ) ) {
    grp.counts <- table( as.factor( survival.data[, "group_num" ] ) )
    legend.labels <- paste( "Cluster", 1:ngrps, paste( "(n=", as.integer(grp.counts), ")", sep="" ) )
  }

  legend( "topright",
          lty=1,
          col=col.map,
          bty = "n",
          legend=legend.labels
         )
  
  if( opt$image.format != "none" ) dev.off()
}