aq-pcoa.in

#!/usr/bin/env Rscript
options(error = quote(dump.frames("pcoa-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
#-d distance method
#-p plot ellipsoids
spec = matrix(c('input', 'i', 1, "character",
'distance','d',2,'character',
'plot_ellipsoids','p',2,'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$distance) ) {
	opt$distance <- "bray"
}
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
dmethod <- opt$distance
ellipsoidConf <- opt$plot_ellipsoids
#Make sure we have a valid distance method
dlist = c('manhattan','euclidean','canberra','bray','kulczynski','jaccard','gower','altGower','morisita','horn','mountford','raup','binomial','chao','cao')

if (! dmethod %in% dlist ) {
	print(paste("Invalid distance method:", dmethod))
	q(status=1)
} 

dir.create(outDir)

pkgTest("vegan")
pkgTest("ape")
if (!is.null(ellipsoidConf) ) {
	pkgTest("car")
}
#Read in and format otu table
print("Reading OTU table")

rawtable <- read.table(otuTable, skip = 1,
comment.char = "", header = TRUE, row.names = 1, sep = "\t")
otutable <- t(rawtable[1:(ncol(rawtable) - 1)])
#Sort numerically the samples
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")
interest <- read.table(headersFile)[1,] == "TRUE";

print("Computing PCoA");

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 <- vegdist(otutable, method=dmethod);
pdf(paste(outDir, "/pcoa-",dmethod,"-biplot.pdf",sep=""),useDingbats=F);

print("Making MDS Plot");

p <- pcoa(d);
for (x in 1 : ncol(mapping)) {
  name <- colnames(mapping)[x]
	metadata <- factor(as.matrix(mapping[,x]))
  fac.len <- length(levels(metadata))
  if (fac.len < 2 || fac.len >= nrow(mapping)) {
    print(paste("Ignoring", name, "because all values are identical/different."))
    next
  }
	metadata_colours <- metadata
	levels(metadata_colours) <- rainbow(fac.len);
	metadata_colours<-unlist(lapply(metadata_colours, function(x) { substr(x, 0, 7); }));
	plot(p$vectors[,1:2], main = paste("PCoA ordination:", name, "\nMethod: '", dmethod, "'", sep=""), col = metadata_colours);
	text(x = p$vectors[,"Axis.1"], y = p$vectors[,"Axis.2"], adj = c(-0.2,0.5), labels = extra[, "Description"], col = metadata_colours);
	#If selected, plot the ellipses around the a priori clusters
	if (!is.null(ellipsoidConf)) {
		ellipse <- dataEllipse(p$vectors[,1],p$vectors[,2],groups=metadata,levels=as.numeric(ellipsoidConf),add=TRUE,plot.points=FALSE,grid=FALSE,center.pch=FALSE,col=rainbow(fac.len))
	}
	legend("bottomleft", legend = unique(metadata), pch = 21, col = unique(metadata_colours));
}

useful <- colnames(mapping)[interest];
if ( length(useful) > 1 ) {
	print("Making biplot");
	biplot(p, apply(log1p(t(newmapping[useful,])), 2, scale, center=TRUE, scale=TRUE), xlabs = extra[, "Description"]);
}


sink(paste(outDir, "/pcoa-",dmethod,".txt",sep=""))
print("Eigenvalues:")
print(as.matrix(p$values$Eigenvalues))
print("Relative Eigenvalues:")
print(as.matrix(p$values$Relative_eig))
sink()