view roi_metrics.Rscript @ 12:4cf71b199381 draft default tip

Uploaded

#!/usr/bin/env Rscript

# Compute metrics (all except PCA) for ROI files.
#
# USAGE: roi_metrics.Rscript -i IN.roi -o OUT.roi

# TODO(holtgrew): Bernd, can you give the English names of the metrics below in the comments?

# Author: Bernd Jagla <[e-mail]>
# Author: Manuel Holtgrewe <[e-mail]>

library('getopt')  # parsing args
library('ngsroi')  # ROI I/O
library('boot')    # k3.linear()

# ----------------------------------------------------------------------------
# Command Line Parsing
# ----------------------------------------------------------------------------

# Parse command line options.
spec = matrix(c(
  "help",     "h", 0, "logical",
  "infile",   "i", 1, "character",
  "outfile",  "o", 1, "character"
), byrow=TRUE, ncol=4);
opt = getopt(spec);

# If help was asked for print a friendly message and exit with a non-zero
# error code.
if (!is.null(opt$help) )
{
  cat(getopt(spec, usage=TRUE));
  q(status=1);
}

# Check arguments.
if (is.null(opt$infile))
{
    write("-i/--infile missing.", stderr());
    q(status=1)
}
if (is.null(opt$out))
{
    write("-o/--outfile missing.", stderr());
    q(status=1)
}

# ----------------------------------------------------------------------------
# Compute metrics.
# ----------------------------------------------------------------------------

# Read input files.
roi = readROI(opt$infile)

# Compute some basic statistics.
roi$min = as.numeric(lapply(roi$counts, min))
roi$median =  as.numeric(lapply(roi$counts, median))
roi$mean = as.numeric(lapply(roi$counts, mean))
roi$quantile75 = as.numeric(lapply(roi$counts, quantile, probs=0.75))
roi$quantile95 = as.numeric(lapply(roi$counts, quantile, probs=0.95))

# Compute area under curve metric.
aoc <- function(v){
  m = max(v)
  return (sum(v/m)/length(v))
}
roi$aoc = as.numeric(lapply(roi$counts, aoc))

# Compute xreaXX metric.
xreaXX <- function(v){
  vo=order(v)
  mx = max(v)
  mi = min(v)
  if((mx-mi)==0) return(0)
  vn=(v-mi)/(mx-mi)
  x=c(1:length(v))/length(v)
  vno=vn[vo]
  return(sum(x-vno)/length(v))
}
roi$xreaXX  =  as.numeric(lapply(roi$counts, xreaXX))

# Compute r3linear metric.
roi$r3linear = as.numeric(lapply(roi$counts, k3.linear))

# Compute distMax3p metric.
distMax3p <- function(v)
{
   return(length(v)-which.max(rev(v))+1)
}
roi$distMax3p = as.numeric(lapply(roi$count, distMax3p))

# Compute distMax5p metric.
distMax5p <- function(v)
{
   return (which.max(v))
}
roi$distMax5p = as.numeric(lapply(roi$count, distMax5p))

# Write output.
writeROI(roi, opt$outfile)

# Quit with success code.
q(status=0)