view roi_pca.Rscript @ 0:61d9bdb6d519 draft

Uploaded

#!/usr/bin/env Rscript

# Compute PCA metric for ROI files.
#
# USAGE: roi_pca.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.x
roi = readROI(opt$infile)

# Compute PCA.
data=roi[,8:length(roi)-1]
pca=princomp(data, na.action=na.exclude)
prj=predict(pca, data)
roi$pca1=prj[,1]
roi$pca2=prj[,2]

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

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