aq-biplot

#!/usr/bin/env Rscript
# Turn QIIME biplot data into something that looks pretty
options(error = quote(dump.frames("biplot-debug", TRUE)))
# Scaling factor for arrow width
arrowScale = 0.1
# Load the mapping, the principle components, the biplot
#   vectors, and the OTUs for the plot
levelname <- tail(commandArgs(trailingOnly = TRUE), 2)[1];
flavour <- tail(commandArgs(trailingOnly = TRUE), 1);
mapping <- read.table("mapping.extra", header = TRUE, 
    comment.char = "", row.names = "X.SampleID", sep = "\t")
pcoa <- read.table(paste("beta_div_pcoa", flavour, "/pcoa_weighted_unifrac_otu_table", flavour, ".txt", sep = ""), 
    comment.char = "", sep = "\t", header = T, row.names = 1)
bipltvec <- read.table(paste("biplot_coords_", levelname, flavour, ".txt", sep = ""), row.names = 1, sep = "\t")
otutable <- read.table(paste("otu_table_summarized_", levelname, flavour, ".txt", sep=""), 
    row.names = 1, header = T, sep = "\t", comment.char = "");
# Drop useless rows (eigenvalues, etc)
pcoa <- pcoa[1:nrow(mapping), ]
svg(paste("biplot_", levelname, flavour, ".svg", sep=""))
# Plot with large points, using the colours determined
#   previouslya
plot(pcoa[, 1], pcoa[, 2], pch = 16, cex = 2, col = as.matrix(mapping["Colour"]), 
    xaxt = "n", yaxt = "n", xlab = "PC1", ylab = "PC2", bty = "n")
# Add a second solid white plotting point, offset slightly
#   to mimic a light source
points(pcoa[, 1] - 0.002, pcoa[, 2] + 0.003, pch = 16, 
    cex = 0.2, col = "white")
for (v in 1:nrow(pcoa)) {
    text(pcoa[v, 1], pcoa[v, 2], mapping[v, "Description"])
}
box(bty = "l", lwd = 3)
# Add arrows, using abundance data to determine linewidth
arrows(rep(0, nrow(bipltvec)), rep(0, nrow(bipltvec)), 
    bipltvec[, 1], bipltvec[, 2], lwd = rowSums(otutable) * arrowScale, 
    len = 0.1, xlab="", ylab="")
text(bipltvec[, 1:2], rownames(otutable))
dev.off() 
warnings()