skbio_ordination_biplot.R

get_ord_dataframe <- function(lines, linestart){
  header_line <- strsplit(lines[linestart], split="\t")[[1]]
  num_rows = as.numeric(header_line[2])
  num_cols = as.numeric(header_line[3])
  if (num_rows == 0){return(data.frame())}
  
  data <- strsplit(lines[(linestart+1):(linestart+num_rows)], split="\t")
  names <- unlist(lapply(data, function(x) x[1]))
  coords <- lapply(data, function(x) strsplit(x[2:(2+num_cols-1)], split="\t"))
  coords <- data.frame(matrix(unlist(coords), nrow=length(coords), byrow=T), stringsAsFactors=FALSE)
  rownames(coords) <- names
  return(coords)
}

shitty_biplot <- function(sites_coords, species_coords, num_arrows=5){
  plot(sites_coords$X1, sites_coords$X2)
  for (i in c(1:num_arrows)){
    arrows(0, 0, as.numeric(species_coords[i, 1]), as.numeric(species_coords[i, 2]))
  }
}

lines <- readLines("results/tutorial/data/ordination.txt")
# format is as follows Eigvals, Proportion explained, Species, Site, Biplot, Site constraints

eigvals <- as.numeric(strsplit(lines[2], split="\t")[[1]])
num_axes <- length(eigvals)
prop_explained <- as.numeric(strsplit(lines[5], split="\t")[[1]])

species_linestart = min(grep("^Species", lines))
species_coords <- get_ord_dataframe(lines, species_linestart)
norm_vec <- function(x) sqrt(sum(as.numeric(x)^2))
species_coords$magnitude <- apply(species_coords, 1, norm_vec)
species_coords <- species_coords[order(-species_coords$magnitude), ]

sites_linestart = min(grep("^Site", lines))
sites_coords <- get_ord_dataframe(lines, sites_linestart)

biplot_linestart = min(grep("^Biplot", lines))
biplot_coords <- get_ord_dataframe(lines, biplot_linestart)

constraints_linestart = min(grep("^Site constraints", lines))
constraints_coords <- get_ord_dataframe(lines, constraints_linestart)

shitty_biplot(sites_coords, species_coords)