aq-pca.in

#!/usr/bin/env Rscript
options(error = quote(dump.frames("pca-debug", TRUE)))

#source("@prefix@/share/@PACKAGE@/biom.R")

pkgTest <- function(x)
{
        if (!require(x,character.only = TRUE))
        {
                install.packages(x,dep=TRUE)
                if(!require(x,character.only = TRUE)) stop("Package not found")
        }
}

pkgTest('getopt')

#Grab arguments
#Arguments required:
#-i input OTU table (tabular format ONLY, JSON libraries much too slow in R)
#-m mapping file
#-e mapping.extra file
#-t headers.txt file
#-o output dir
spec = matrix(c('input', 'i', 1, "character",
'mapping','m',1,"character",
'mapping.extra','e',1,"character",
'headers','t',1,'character',
'output','o',1,"character",
'help','h',2,"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)
}

if ( is.null(opt$input) ) {
  print("Input OTU table required.")
  q(status=1)
}
if ( is.null(opt$mapping) ) {
  print("Mapping file required.")
  q(status=1)
}
if ( is.null(opt$mapping.extra) ) {
  print("mapping.txt file required.")
  q(status=1)
}
if ( is.null(opt$headers) ) {
  print("headers.txt file required.")
  q(status=1)
}
if ( is.null(opt$output) ) {
  opt$output <- getwd()
}

outDir <- opt$output
otuTable <- opt$input
mappingFile <- opt$mapping
mappingExtra <- opt$mapping.extra
headersFile <- opt$headers

dir.create(outDir)

rawtable <- read.table(otuTable, skip = 1,
comment.char = "", header = TRUE, row.names = 1, sep = "\t")
otutable <- t(rawtable[1:(ncol(rawtable) - 1)])
otutable <- otutable[order(as.integer(sub("X","", rownames(otutable)))),]

print("Reading mapping");
mapping <- read.table(mappingFile, header = TRUE, comment.char = "", row.names = "X.SampleID", sep = "\t")
print("Reading mapping extra");
extra <- read.table(mappingExtra, header = TRUE,
    comment.char = "", row.names = "X.SampleID", sep = "\t")
print("Reading header");
interest <- read.table(headersFile)[1,] == "TRUE";
print("Computing PCA");

for(i in 0 : length(interest)) {
	if (interest[i] && length(levels(factor(as.matrix(mapping[, i])))) == 1) {
		interest[i] <- FALSE
		print(paste("Ignoring", colnames(mapping)[i], "because all values are identical."))
	}
}
newmapping <- t(mapping[,interest]);

# For non-numeric sample names, comment out the following line
colnames(newmapping) <- paste("X", rownames(mapping), sep = "");
rownames(newmapping) <- colnames(mapping)[interest];

d <- rbind(t(otutable), newmapping);
p <- prcomp(t(d), scale = TRUE);

pdf(paste(outDir, "/pca-biplot.pdf",sep=""));

print("Making Scree plot");
plot(p);

print("Making MDS plot");
plot(p$x, asp = p$sdev[2]/p$sdev[1], col = as.matrix(extra[ ,"Colour"]));
text(x = p$x[, "PC1"], y = p$x[, "PC2"], labels = extra[ ,"Description"], col = as.matrix(extra[ ,"Colour"]));

print("Making biplot");
contrib <- apply(p$rotation, 1, function(x) {sqrt(sum(x[1:2]^2));});

# What variables (rows) do you want to display as arrows?
#useful <- contrib > (mean(contrib) + 2*sd(contrib));
useful <- colnames(mapping)[interest];
if (length(useful) > 1) {
	biplot(p$x[, c("PC1", "PC2")], p$rotation[useful, c("PC1", "PC2")], xlabs = extra[, "Description"], asp = p$sdev[2]/p$sdev[1]);
}