view cluster.tools/determine.IPL.threshold.R @ 0:0decf3fd54bc draft

Uploaded

#!/usr/bin/env Rscript

##usage, options and doc goes here
argspec <- c("determine.IPL.threshold.R takes an IPL result, and determines a statistically sound threshold to use

        Usage: 
                determine.IPL.threshold.R -d <IPL_data_file>
        Optional:
                -o output.rdata  ## rdata output file (contains variables used for calculation, for those who want to review them
                -f filter type # must be either modulated, active, or inactive
                -p percent of samples passing (must be value on [0,1]
\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)

spec <- matrix( c( "data.fname",          "d", 1, "character",
                   "output.rdata",        "o", 2, "character",
                   "filter.type",         "f", 2, "character",
                   "perc.pass",           "p", 2, "numeric",
                   "selection.criteria",  "s", 2, "character",
                   "output.report.dir",   "r", 2, "character",
                   "output.report.html",  "h", 2, "character"
                  ),
                nc=4,
                byrow=T
               )
opt <- getopt( spec=spec )

## default params for non-required params
if ( is.null( opt$filter.type ) ) { opt$filter.type <- 'modulated' }
if ( is.null( opt$perc.pass ) ) { opt$perc.pass <- 1/3 }
if ( is.null( opt$selection.criteria ) ) { opt$selection.criteria <- 'max_diffs' }
if ( is.null( opt$output.report.dir ) ) { opt$output.report.dir <- "report" }
if ( is.null( opt$output.report.html ) ) { opt$output.report.html <- "report/index.html" }
if ( is.null( opt$output.rdata ) ) { opt$output.rdata <- "output.rdata" }
if ( opt$perc.pass < 0  ) {
  stop( "please specify a positive number for the percentage of samples that pass the filter (if applicable)" )
}

if (!file.exists(opt$output.report.dir)){
    dir.create(opt$output.report.dir)
}

                               
data <- as.matrix( read.delim( opt$data.fname, row.names=1, check.names=FALSE ) )
genes <- rownames( data )
genes <- genes[ !grepl( "abstract|family|complex", genes ) ]
data <- data[ genes, ]

nulls.mat <- grepl( "na_", colnames( data ) )
reals <- ! nulls.mat
nulls.mat <- data[ , nulls.mat ]
reals.mat <- data[, reals ]
if ( ncol( nulls.mat ) == 0 ) stop( "no nulls were in the file provided!\n" )
if ( ncol( reals.mat ) == 0 ) stop( "no reals were in the file provided!\n" )


if ( opt$filter.type == 'modulated' ) {
  reals.mat <- abs( reals.mat )
  nulls.mat <- abs( nulls.mat )
} else {
  if ( opt$filter.type == "inactive" ) {
    reals.mat <- -reals.mat
    nulls.mat <- -nulls.mat
  }
} 


##  we only look at the larger 50% of the possible IPL values
##  as possible  thresholds to use (since the lower 50% are almost
##  always uninformative)
thresholds <- unique( quantile( reals.mat,
                                seq( 0.5, 1, by=0.001 ) ) )
thresholds <- c( thresholds,
                 quantile( nulls.mat,
                           seq( 0.5, 1, by=0.001 ) ) )
thresholds <- unique( sort( thresholds ) )


get.num.filtered.feats <- function( mat,
                                    threshold,
                                    perc.samples.passing=1/3 ) {
  feat.vect <- apply( mat,
                      1,
                      function(x) {
                        tmp <- sum( x > threshold )
                        if ( perc.samples.passing >=1 ) {
                          return( tmp >= perc.samples.passing )
                        } else {
                          return( tmp > floor( perc.samples.passing * length(x) ) )
                        }
                      }
                     )
  return( sum( feat.vect ) )
}


real.feats <- null.feats <- length( genes )
chisq.pvals <- binom.pvals <- numeric()

for ( i in 1:length( thresholds ) ) {

  nul.feats.this.thresh <- get.num.filtered.feats( mat=nulls.mat, threshold=thresholds[i], perc.samples.passing=opt$perc.pass )
  ## limit the maximum threshold to one where there are at least 75 valid points
  ##  because if there are fewer nulls than that, it heavily skews the probability
  if ( nul.feats.this.thresh < 50 ) break
  
  null.feats[ i ] <- nul.feats.this.thresh
  real.feats[ i ] <- get.num.filtered.feats( mat=reals.mat, threshold=thresholds[i], perc.samples.passing=opt$perc.pass )  
  
  ## only calculate if there are more real features than nulls, otherwise, give a p-value of 1
  if ( null.feats[i] < real.feats[i] ) {
    p <- null.feats[i]/nrow( nulls.mat )
    sd <- ( nrow( nulls.mat ) * p * (1-p ) )^0.5

    ## binomial test
    p <- -pnorm( q=real.feats[i],
                 mean=null.feats[i],
                 sd=sd,
                 log.p=TRUE,
                 lower.tail=FALSE )

    ##chisq test
    chi <- ( real.feats[i] - null.feats[i] )^2
    chi <- chi/(null.feats[i])^2
    chi <- -pchisq( chi, 1, log.p=TRUE, lower=FALSE )
  } else {
    p <- chi <- 0 ## 0 == -log(1)
  }

  binom.pvals <- c( binom.pvals, p )
  chisq.pvals <- c( chisq.pvals, chi )

  if ( length( chisq.pvals ) != i ) {
    stop( "lengths differ\n" )
  }
}



##names( binom.pvals ) <- names( chisq.pvals ) <- thresholds
diffs <- real.feats - null.feats
if ( opt$selection.criteria == "max_diffs" ) {
  max.diff <- max( diffs )
  opt.thresh <- which( diffs %in% max.diff )
} else if ( opt$selection.criteria == "binomial" ) {
  max.bin <- max( binom.pvals )
  opt.thresh <- which( binom.pvals %in% max.bin )
} else if ( opt$selection.criteria == "chisq" ) {
  max.chi <- max( chisq.pvals )
  opt.thresh <- which( chisq.pvals %in% max.chi )
}

opt.thresh <- mean( c( thresholds[ opt.thresh ], thresholds[ (opt.thresh-1) ] ) )
opt.thresh <- signif( opt.thresh, 4 )

  
##corrected.binom.pvals <- binom.pvals + log( length(thresholds) )
##binom.pvals <- binom.pvals - log( length(thresholds) )
##corrected.chisq.pvals <- chisq.pvals + log( length(thresholds) )
##chisq.pvals <- chisq.pvals - log( length(thresholds) )


eval.thresh <- thresholds[ 1:length( real.feats ) ]
##plot.new(); screens <- split.screen( c( 4,1 ) )
##postscript( "threshold.comparison.ps", paper='letter', horizontal=F )
##png.fname <- file.path( opt$output.report.dir, "IPL.threshold.determination.png")
##plot.dev <- png( png.fname,
##                 width=11,
##                 height=8.5,
##                 units='in',
##                 res=72 )
##par( mar=rep(0,4) )
##screens <- split.screen( c( 4,1 ) )

png.fname <- file.path( opt$output.report.dir, "01.num.feats.IPL.threshold.determination.png")
plot.dev <- png( png.fname,
                 width=11,
                 height=8.5,
                 units='in',
                 res=72 )
par( mar=c(2.25,3,1.5,0.5) )
plot( eval.thresh, null.feats, type='l', lwd=2, col='blue', cex.axis=0.75 )
lines( eval.thresh, real.feats, type='l', lwd=2, col='black', cex.axis=0.75 )
abline( v=opt.thresh )
legend( "topright", c( "Real", "Null" ), lwd=2, col=c('black', 'blue' ) )
mtext( "Number of Genes Passing Threshold", font=2 )
mtext( "IPL Threshold", 1, font=2, line=1.5 )
mtext( "Number of Genes", 2, font=2, line=1.8 )
dev.off()


png.fname <- file.path( opt$output.report.dir, "02.diffs.IPL.threshold.determination.png")
plot.dev <- png( png.fname,
                 width=11,
                 height=8.5,
                 units='in',
                 res=72 )
##screen( screen()+1 )
par( mar=c(2.25,3,1.5,0.5) )
plot( eval.thresh, diffs, type='l', lwd=2, col='black', cex.axis=0.75 )
abline( v=opt.thresh )
mtext( "Difference between number of Real & Null genes passing Threshold", font=2 )
mtext( "IPL Threshold", 1, font=2, line=1.5 )
mtext( "Number of Genes", 2, font=2, line=1.8 )
dev.off()



png.fname <- file.path( opt$output.report.dir, "03.chisq.IPL.threshold.determination.png")
plot.dev <- png( png.fname,
                 width=11,
                 height=8.5,
                 units='in',
                 res=72 )
##screen( screen()+1 )
par( mar=c(2.25,3,1.5,0.5) )
plot( eval.thresh, chisq.pvals, type='l', lwd=2, col='red', cex.axis=0.75 )
abline( v=opt.thresh )
mtext( "Chi-sq p-values", font=2 )
mtext( "IPL Threshold", 1, font=2, line=1.5 )
mtext( "-Log p-value", 2, font=2, line=1.8 )
dev.off()



png.fname <- file.path( opt$output.report.dir, "04.binom.IPL.threshold.determination.png")
plot.dev <- png( png.fname,
                 width=11,
                 height=8.5,
                 units='in',
                 res=72 )
##screen( screen()+1 )
par( mar=c(2.25,3,1.5,0.5) )
plot( eval.thresh, binom.pvals, type='l', lwd=2, col='green', cex.axis=0.75 )
abline( v=opt.thresh )
mtext( "Binomial p-values", font=2 )
mtext( "IPL Threshold", 1, font=2, line=1.5 )
mtext( "-Log p-value", 2, font=2, line=1.8 )
dev.off()

##close.screen( all=T ); dev.off()

report_str = paste( "The threshold to use for consensus clustering filtering is ", opt.thresh, "\n", sep="" )

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

filter.type <- opt$filter.type
perc.pass <- opt$perc.pass
save( file=opt$output.rdata, thresholds, diffs, binom.pvals, chisq.pvals, real.feats, null.feats, data, filter.type, perc.pass, opt.thresh )