diff maaslin-4450aa4ecc84/src/lib/scriptBiplotTSV.R @ 1:a87d5a5f2776

Uploaded the version running on the prod server

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/maaslin-4450aa4ecc84/src/lib/scriptBiplotTSV.R	Sun Feb 08 23:08:38 2015 -0500
@@ -0,0 +1,515 @@
+#!/usr/bin/env Rscript
+
+library(vegan)
+library(optparse)
+
+funcGetCentroidForMetadatum <- function(
+### Given a binary metadatum, calculate the centroid of the samples associated with the metadata value of 1
+# 1. Get all samples that have the metadata value of 1
+# 2. Get the x and y coordinates of the selected samples
+# 3. Get the median value for the x and ys
+# 4. Return those coordinates as the centroid's X and Y value
+vfMetadata,
+### Logical or integer (0,1) vector, TRUE or 1 values indicate correspoinding samples in the
+### mSamplePoints which will be used to define the centroid
+mSamplePoints
+### Coordinates (columns;n=2) of samples (rows) corresponding to the vfMetadata
+){
+  # Check the lengths which should be equal
+  if(length(vfMetadata)!=nrow(mSamplePoints))
+  {
+    print(paste("funcGetCentroidForMetadata::Error: Should have received metadata and samples of the same length, received metadata length ",length(vfMetadata)," and sample ",nrow(mSamplePoints)," length.",sep=""))
+    return( FALSE )
+  }
+
+  # Get all the samples that have the metadata value of 1
+  viMetadataSamples = which(as.integer(vfMetadata)==1)
+
+  # Get the x and y coordinates for the selected samples
+  mSelectedPoints = mSamplePoints[viMetadataSamples,]
+
+  # Get the median value for the x and the ys
+  if(!is.null(nrow(mSelectedPoints)))
+  {
+    return( list(x=median(mSelectedPoints[,1],na.rm = TRUE),y=median(mSelectedPoints[,2],na.rm = TRUE)) )
+  } else {
+    return( list(x=mSelectedPoints[1],y=mSelectedPoints[2]) )
+  }
+}
+
+funcGetMaximumForMetadatum <- function(
+### Given a continuous metadata
+### 1. Use the x and ys from mSamplePoints for coordinates and the metadata value as a height (z)
+### 2. Use lowess to smooth the landscape
+### 3. Take the maximum of the landscape
+### 4. Return the coordiantes for the maximum as the centroid
+vdMetadata,
+### Continuous (numeric or integer) metadata
+mSamplePoints
+### Coordinates (columns;n=2) of samples (rows) corresponding to the vfMetadata
+){
+  # Work with data frame
+  if(class(mSamplePoints)=="matrix")
+  {
+    mSamplePoints = data.frame(mSamplePoints)
+  }
+  # Check the lengths of the dataframes and the metadata
+  if(length(vdMetadata)!=nrow(mSamplePoints))
+  {
+    print(paste("funcGetMaximumForMetadatum::Error: Should have received metadata and samples of the same length, received metadata length ",length(vdMetadata)," and sample ",nrow(mSamplePoints)," length.",sep=""))
+    return( FALSE )
+  }
+
+  # Add the metadata value to the points
+  mSamplePoints[3] = vdMetadata
+  names(mSamplePoints) = c("x","y","z") 
+
+  # Create lowess to smooth the surface
+  # And calculate the fitted heights
+  # x = sample coordinate 1
+  # y = sample coordinate 2
+  # z = metadata value
+  loessSamples = loess(z~x*y, data=mSamplePoints, degree = 1, normalize = FALSE, na.action=na.omit)
+
+  # Naively get the max
+  vdCoordinates = loessSamples$x[which(loessSamples$y==max(loessSamples$y)),]
+  return(list(lsmod = loessSamples, x=vdCoordinates[1],y=vdCoordinates[2]))
+}
+
+funcMakeShapes <- function(
+### Takes care of defining shapes for the plot
+dfInput,
+### Data frame of metadata measurements
+sShapeBy,
+### The metadata to shape by
+sShapes,
+### List of custom metadata (per level if factor).
+### Should correspond to the number of levels in shapeBy; the format is level:shape,level:shape for example HighLuminosity:14,LowLuminosity:2,HighPH:10,LowPH:18 
+cDefaultShape
+### Shape to default to if custom shapes are not used
+){
+  lShapes = list()
+  vsShapeValues = c()
+  vsShapeShapes = c()
+  vsShapes = c()
+  sMetadataId = sShapeBy
+
+  # Set default shape, color, and color ranges 
+  if(!is.null(cDefaultShape))
+  {
+    # Default shape should be an int for the int pch options
+    if(!is.na(as.integer(cDefaultShape)))
+    {
+      cDefaultShape = as.integer(cDefaultShape)
+    }
+  } else {
+    cDefaultShape = 16
+  }
+
+  # Make shapes
+  vsShapes = rep(cDefaultShape,nrow(dfInput))
+
+  if(!is.null(sMetadataId))
+  {
+    if(is.null(sShapes))
+    {
+      vsShapeValues = unique(dfInput[[sMetadataId]])
+      vsShapeShapes = 1:length(vsShapeValues)
+    } else {
+      # Put the markers in the order of the values)
+      vsShapeBy = unlist(strsplit(sShapes,","))
+      for(sShapeBy in vsShapeBy)
+      {
+        vsShapeByPieces = unlist(strsplit(sShapeBy,":"))
+        lShapes[vsShapeByPieces[1]] = as.integer(vsShapeByPieces[2])
+      }
+      vsShapeValues = names(lShapes)
+   }
+
+    # Shapes in the correct order
+    if(!is.null(sShapes))
+    {
+      vsShapeShapes = unlist(lapply(vsShapeValues,function(x) lShapes[[x]]))
+    }
+    vsShapeValues = paste(vsShapeValues)
+
+    # Make the list of shapes
+    for(iShape in 1:length(vsShapeValues))
+    {
+      vsShapes[which(paste(dfInput[[sMetadataId]])==vsShapeValues[iShape])]=vsShapeShapes[iShape]
+    }
+
+    # If they are all numeric characters, make numeric
+    viIntNas = which(is.na(as.integer(vsShapes)))
+    viNas = which(is.na(vsShapes))
+    if(length(setdiff(viIntNas,viNas))==0)
+    {
+      vsShapes = as.integer(vsShapes)
+    } else {
+      print("funcMakeShapes::Error: Please supply numbers 1-25 for shape in the -y,--shapeBy option")
+      vsShapeValues = c()
+      vsShapeShapes = c()
+    }
+  }
+  return(list(PlotShapes=vsShapes,Values=vsShapeValues,Shapes=vsShapeShapes,ID=sMetadataId,DefaultShape=cDefaultShape))
+}
+
+### Global defaults
+c_sDefaultColorBy = NULL
+c_sDefaultColorRange = "orange,cyan"
+c_sDefaultTextColor = "black"
+c_sDefaultArrowColor = "cyan"
+c_sDefaultArrowTextColor = "Blue"
+c_sDefaultNAColor = "grey"
+c_sDefaultShapeBy = NULL
+c_sDefaultShapes = NULL
+c_sDefaultMarker = "16"
+c_sDefaultRotateByMetadata = NULL
+c_sDefaultResizeArrow = 1
+c_sDefaultTitle = "Custom Biplot of Bugs and Samples - Metadata Plotted with Centroids"
+c_sDefaultOutputFile = NULL
+
+### Create command line argument parser
+pArgs <- OptionParser( usage = "%prog last_metadata input.tsv" )
+
+# Selecting features to plot
+pArgs <- add_option( pArgs, c("-b", "--bugs"), type="character", action="store", default=NULL, dest="sBugs", metavar="BugsToPlot", help="Comma delimited list of data to plot as text. Bug|1,Bug|2")
+pArgs <- add_option( pArgs, c("-m", "--metadata"), type="character", action="store", default=NULL, dest="sMetadata", metavar="MetadataToPlot", help="Comma delimited list of metadata to plot as arrows. metadata1,metadata2,metadata3")
+
+# Colors
+pArgs <- add_option( pArgs, c("-c", "--colorBy"), type="character", action="store", default=c_sDefaultColorBy, dest="sColorBy", metavar="MetadataToColorBy", help="The id of the metadatum to use to make the marker colors. Expected to be a continuous metadata.")
+pArgs <- add_option( pArgs, c("-r", "--colorRange"), type="character", action="store", default=c_sDefaultColorRange, dest="sColorRange", metavar="ColorRange", help=paste("Colors used to color the samples; a gradient will be formed between the color.Default=", c_sDefaultColorRange))
+pArgs <- add_option( pArgs, c("-t", "--textColor"), type="character", action="store", default=c_sDefaultTextColor, dest="sTextColor", metavar="TextColor", help=paste("The color bug features will be plotted with as text. Default =", c_sDefaultTextColor))
+pArgs <- add_option( pArgs, c("-a", "--arrowColor"), type="character", action="store", default=c_sDefaultArrowColor, dest="sArrowColor", metavar="ArrowColor", help=paste("The color metadata features will be plotted with as an arrow and text. Default", c_sDefaultArrowColor))
+pArgs <- add_option( pArgs, c("-w", "--arrowTextColor"), type="character", action="store", default=c_sDefaultArrowTextColor, dest="sArrowTextColor", metavar="ArrowTextColor", help=paste("The color for the metadata text ploted by the head of the metadata arrow. Default", c_sDefaultArrowTextColor))
+pArgs <- add_option(pArgs, c("-n","--plotNAColor"), type="character", action="store", default=c_sDefaultNAColor, dest="sPlotNAColor", metavar="PlotNAColor", help=paste("Plot NA values as this color. Example -n", c_sDefaultNAColor))
+
+# Shapes
+pArgs <- add_option( pArgs, c("-y", "--shapeby"), type="character", action="store", default=c_sDefaultShapeBy, dest="sShapeBy", metavar="MetadataToShapeBy", help="The metadata to use to make marker shapes. Expected to be a discrete metadatum. An example would be -y Environment")
+pArgs <- add_option( pArgs, c("-s", "--shapes"), type="character", action="store", default=c_sDefaultShapes, dest="sShapes", metavar="ShapesForPlotting", help="This is to be used to specify the shapes to use for plotting. Can use numbers recognized by R as shapes (see pch). Should correspond to the number of levels in shapeBy; the format is level:shape,level:shape for example HighLuminosity:14,LowLuminosity:2,HighPH:10,LowPH:18 . Need to specify -y/--shapeBy for this option to work.")
+pArgs <- add_option( pArgs, c("-d", "--defaultMarker"), type="character", action="store", default=c_sDefaultMarker, dest="sDefaultMarker", metavar="DefaultColorMarker", help="Default shape for markers which are not otherwise indicated in --shapes, can be used for unspecified values or NA. Must not be a shape in --shapes.")
+
+# Plot manipulations
+pArgs <- add_option( pArgs, c("-e","--rotateByMetadata"), type="character", action="store", default=c_sDefaultRotateByMetadata, dest="sRotateByMetadata", metavar="RotateByMetadata", help="Rotate the ordination by a metadata. Give both the metadata and value to weight it by. The larger the weight, the more the ordination is influenced by the metadata. If the metadata is continuous, use the metadata id; if the metadata is discrete, the ordination will be by one of the levels so use the metadata ID and level seperated by a '_'. Discrete example -e Environment_HighLumninosity,100 ; Continuous example -e Environment,100 .")
+pArgs <- add_option( pArgs, c("-z","--resizeArrow"), type="numeric", action="store", default=c_sDefaultResizeArrow, dest="dResizeArrow", metavar="ArrowScaleFactor", help="A constant to multiple the length of the arrow to expand or shorten all arrows together. This will not change the angle of the arrow nor the relative length of arrows to each other.")
+
+# Misc
+pArgs <- add_option( pArgs, c("-i", "--title"), type="character", action="store", default=c_sDefaultTitle, dest="sTitle", metavar="Title", help="This is the title text to add to the plot.")
+pArgs <- add_option( pArgs, c("-o", "--outputfile"), type="character", action="store", default=c_sDefaultOutputFile, dest="sOutputFileName", metavar="OutputFile", help="This is the name for the output pdf file. If an output file is not given, an output file name is made based on the input file name.")
+
+funcDoBiplot <- function(
+### Perform biplot. Samples are markers, bugs are text, and metadata are text with arrows. Markers and bugs are dtermined usiing NMDS and Bray-Curtis dissimilarity. Metadata are placed on the ordination in one of two ways: 1. Factor data - for each level take the ordination points for the samples that have that level and plot the metadata text at the average orindation point. 2. For continuous data - make a landscape (x and y form ordination of the points) and z (height) as the metadata value. Use a lowess line to get the fitted values for z and take the max of the landscape. Plot the metadata text at that smoothed max.
+sBugs,
+### Comma delimited list of data to plot as text. Bug|1,Bug|2
+sMetadata,
+### Comma delimited list of metadata to plot as arrows. metadata1,metadata2,metadata3.
+sColorBy = c_sDefaultColorBy,
+### The id of the metadatum to use to make the marker colors. Expected to be a continuous metadata.
+sColorRange = c_sDefaultColorRange,
+### Colors used to color the samples; a gradient will be formed between the color. Example orange,cyan
+sTextColor = c_sDefaultTextColor,
+### The color bug features will be plotted with as text. Example black
+sArrowColor = c_sDefaultArrowColor,
+### The color metadata features will be plotted with as an arrow and text. Example cyan
+sArrowTextColor = c_sDefaultArrowTextColor,
+### The color for the metadata text ploted by the head of the metadata arrow. Example Blue
+sPlotNAColor = c_sDefaultNAColor,
+### Plot NA values as this color. Example grey
+sShapeBy = c_sDefaultShapeBy,
+### The metadata to use to make marker shapes. Expected to be a discrete metadatum.
+sShapes = c_sDefaultShapes,
+### This is to be used to specify the shapes to use for plotting. Can use numbers recognized by R as shapes (see pch). Should correspond to the number of levels in shapeBy; the format is level:shape,level:shape for example HighLuminosity:14,LowLuminosity:2,HighPH:10,LowPH:18 .  Works with sShapesBy.
+sDefaultMarker = c_sDefaultMarker,
+### The default marker shape to use if shapes are not otherwise indicated.
+sRotateByMetadata = c_sDefaultRotateByMetadata,
+### Metadata and value to rotate by. example Environment_HighLumninosity,100
+dResizeArrow = c_sDefaultResizeArrow,
+### Scale factor to resize tthe metadata arrows
+sTitle = c_sDefaultTitle,
+### The title for the figure.
+sInputFileName,
+### File to input (tsv file: tab separated, row = sample file)
+sLastMetadata,
+### Last metadata that seperates data and metadata
+sOutputFileName = c_sDefaultOutputFile
+### The file name to save the figure.
+){
+  print("IN Biplot")
+  # Define the colors
+  vsColorRange = c("blue","orange")
+  cDefaultColor = "black"
+  if(!is.null(sColorRange))
+  {
+    vsColorRange = unlist(strsplit(sColorRange,","))
+  }
+
+  # List of bugs to plot
+  # If there is a list it needs to be more than one.
+  vsBugsToPlot = c()
+  if(!is.null(sBugs))
+  {
+    vsBugsToPlot = unlist(strsplit(sBugs,","))
+  }
+
+  print("vsBugsToPlot")
+  print(vsBugsToPlot)
+  # Metadata to plot
+  vsMetadata = c()
+  if(!is.null(sMetadata))
+  {
+    vsMetadata = unlist(strsplit(sMetadata,","))
+  }
+
+  print("vsMetadata")
+  print(vsMetadata)
+  ### Load table
+  if(class(sInputFileName)=="character")
+  {
+    dfInput = read.table(sInputFileName, sep = "\t", header=TRUE)
+    names(dfInput) = unlist(lapply(names(dfInput),function(x) gsub(".","|",x,fixed=TRUE)))
+    row.names(dfInput) = dfInput[,1]
+    dfInput = dfInput[-1]
+  } else {dfInput = sInputFileName}
+
+  ### Get positions of all metadata or all data
+  iLastMetadata = which(names(dfInput)==sLastMetadata)
+  viMetadata = 1:iLastMetadata
+  viData = (iLastMetadata+1):ncol(dfInput)
+
+  ### Dummy the metadata if discontinuous
+  ### Leave the continous metadata alone but include
+  listMetadata = list()
+  vsRowNames = c()
+  viContinuousMetadata = c()
+  for(i in viMetadata)
+  {
+    print( names( dfInput )[i] )
+    vCurMetadata = unlist(dfInput[i])
+    if( ( is.numeric(vCurMetadata)||is.integer(vCurMetadata) )  && ( length( unique( vCurMetadata ) ) >= c_iNonFactorLevelThreshold ) )
+    {
+      vCurMetadata[which(is.na(vCurMetadata))] = mean(vCurMetadata,na.rm=TRUE)
+      listMetadata[[length(listMetadata)+1]] = vCurMetadata
+      vsRowNames = c(vsRowNames,names(dfInput)[i])
+      viContinuousMetadata = c(viContinuousMetadata,length(listMetadata))
+    } else {
+      vCurMetadata = as.factor(vCurMetadata)
+      vsLevels = levels(vCurMetadata)
+      for(sLevel in vsLevels)
+      { 
+        vNewMetadata = rep(0,length(vCurMetadata))
+        vNewMetadata[which(vCurMetadata == sLevel)] = 1
+        listMetadata[[length(listMetadata)+1]] = vNewMetadata
+        vsRowNames = c(vsRowNames,paste(names(dfInput)[i],sLevel,sep="_"))
+      }
+    }
+  }
+
+  # Convert to data frame
+  dfDummyMetadata = as.data.frame(sapply(listMetadata,rbind))
+  names(dfDummyMetadata) = vsRowNames
+  iNumberMetadata = ncol(dfDummyMetadata)
+
+  # Data to use in ordination in NMDS
+  # All cleaned bug data
+  dfData = dfInput[viData]
+
+  # If rotating the ordination by a metadata
+  # 1. Add in the metadata as a bug
+  # 2. Multiply the bug by the weight
+  # 3. Push this through the NMDS
+  if(!is.null(sRotateByMetadata))
+  {
+    vsRotateMetadata = unlist(strsplit(sRotateByMetadata,","))
+    sMetadata = vsRotateMetadata[1]
+    dWeight = as.numeric(vsRotateMetadata[2])
+    sOrdinationMetadata = dfDummyMetadata[sMetadata]*dWeight
+    dfData[sMetadata] = sOrdinationMetadata
+  }
+
+  # Run NMDS on bug data (Default B-C)
+  # Will have species and points because working off of raw data
+  mNMDSData = metaMDS(dfData,k=2)
+
+  ## Make shapes
+  # Defines the shapes and the metadata they are based on
+  # Metadata to use as shapes
+  lShapeInfo = funcMakeShapes(dfInput=dfInput, sShapeBy=sShapeBy, sShapes=sShapes, cDefaultShape=sDefaultMarker)
+
+  sMetadataShape = lShapeInfo[["ID"]]
+  vsShapeValues = lShapeInfo[["Values"]]
+  vsShapeShapes = lShapeInfo[["Shapes"]]
+  vsShapes = lShapeInfo[["PlotShapes"]]
+  cDefaultShape = lShapeInfo[["DefaultShape"]]
+
+  # Colors
+  vsColors = rep(cDefaultColor,nrow(dfInput))
+  vsColorValues = c()
+  vsColorRBG = c()
+  if(!is.null(sColorBy))
+  {
+    vsColorValues = paste(sort(unique(unlist(dfInput[[sColorBy]])),na.last=TRUE))
+    iLengthColorValues = length(vsColorValues)
+
+    vsColorRBG = lapply(1:iLengthColorValues/iLengthColorValues,colorRamp(vsColorRange))
+    vsColorRBG = unlist(lapply(vsColorRBG, function(x) rgb(x[1]/255,x[2]/255,x[3]/255)))
+
+    for(iColor in 1:length(vsColorRBG))
+    {
+      vsColors[which(paste(dfInput[[sColorBy]])==vsColorValues[iColor])]=vsColorRBG[iColor]
+    }
+
+    #If NAs are seperately given color, then color here
+    if(!is.null(sPlotNAColor))
+    {
+      vsColors[which(is.na(dfInput[[sColorBy]]))] = sPlotNAColor
+      vsColorRBG[which(vsColorValues=="NA")] = sPlotNAColor
+    }
+  }
+
+  print("names(dfDummyMetadata)")
+  print(names(dfDummyMetadata))
+
+  # Reduce the bugs down to the ones in the list to be plotted
+  viBugsToPlot = which(row.names(mNMDSData$species) %in% vsBugsToPlot)
+  viMetadataDummy = which(names(dfDummyMetadata) %in% vsMetadata)
+
+  print("viBugsToPlot")
+  print(viBugsToPlot)
+  print("viMetadataDummy")
+  print(names(dfDummyMetadata)[viMetadataDummy])
+
+  # Build the matrix of metadata coordinates
+  mMetadataCoordinates = matrix(rep(NA, iNumberMetadata*2),nrow=iNumberMetadata)
+  for( i in 1:iNumberMetadata )
+  {
+    lxReturn = NA
+    if( i %in% viContinuousMetadata )
+    {
+      lxReturn = funcGetMaximumForMetadatum(dfDummyMetadata[[i]],mNMDSData$points)
+    } else {
+      lxReturn = funcGetCentroidForMetadatum(dfDummyMetadata[[i]],mNMDSData$points)
+    }
+    mMetadataCoordinates[i,] = c(lxReturn$x,lxReturn$y)
+  }
+  row.names(mMetadataCoordinates) = vsRowNames
+
+  # Plot the biplot with the centroid constructed metadata coordinates
+  if(length(viMetadataDummy)==0)
+  {
+    viMetadataDummy = 1:nrow(mMetadataCoordinates)
+  }
+
+  # Plot samples
+  # Make output name
+  if(is.null(sOutputFileName))
+  {
+    viPeriods = which(sInputFileName==".")
+    if(length(viPeriods)>0)
+    {
+      sOutputFileName = paste(OutputFileName[1:viPeriods[length(viPeriods)]],"pdf",sep=".")
+    } else {
+      sOutputFileName = paste(sInputFileName,"pdf",sep=".")
+    }
+  }
+
+  pdf(sOutputFileName, useDingbats=FALSE)
+  plot(mNMDSData$points, xlab=paste("NMDS1","Stress=",mNMDSData$stress), ylab="NMDS2", pch=vsShapes, col=vsColors)
+  title(sTitle,line=3)
+
+  # Plot Bugs
+  mPlotBugs = mNMDSData$species[viBugsToPlot,]
+  if(length(viBugsToPlot)==1)
+  {
+    text(x=mPlotBugs[1],y=mPlotBugs[2],labels=row.names(mNMDSData$species)[viBugsToPlot],col=sTextColor)
+  } else if(length(viBugsToPlot)>1){
+    text(x=mPlotBugs[,1],y=mPlotBugs[,2],labels=row.names(mNMDSData$species)[viBugsToPlot],col=sTextColor)
+  }
+
+  # Add alternative axes
+  axis(3, col=sArrowColor)
+  axis(4, col=sArrowColor)
+  box(col = "black")
+
+  # Plot Metadata
+  if(length(viMetadataDummy)>0)
+  {
+    for(i in viMetadataDummy)
+    {
+      curCoordinates = mMetadataCoordinates[i,]
+      curCoordinates = curCoordinates * dResizeArrow
+      # Plot Arrow
+      arrows(0,0, curCoordinates[1] * 0.8, curCoordinates[2] * 0.8, col=sArrowColor, length=0.1 )
+    }
+    # Plot text
+    if(length(viMetadataDummy)==1)
+    {
+      text(x=mMetadataCoordinates[viMetadataDummy,][1]*dResizeArrow*0.8, y=mMetadataCoordinates[viMetadataDummy,][2]*dResizeArrow*0.8, labels=row.names(mMetadataCoordinates)[viMetadataDummy],col=sArrowTextColor)
+    } else {
+      text(x=mMetadataCoordinates[viMetadataDummy,1]*dResizeArrow*0.8, y=mMetadataCoordinates[viMetadataDummy,2]*dResizeArrow*0.8, labels=row.names(mMetadataCoordinates)[viMetadataDummy],col=sArrowTextColor)
+    }
+  }
+
+  # Create Legend
+  # The text default is the colorMetadata_level (one per level) plus the ShapeMetadata_level (one per level)
+  # The color default is already determined colors plus grey for shapes.
+  sLegendText = c(paste(vsColorValues,sColorBy,sep="_"),paste(sMetadataShape,vsShapeValues,sep="_"))
+  sLegendColors = c(vsColorRBG,rep(cDefaultColor,length(vsShapeValues)))
+
+  # If the color values are numeric
+  # Too many values may be given in the legend (given they may be a continuous range of values)
+  # To reduce this they are summarized instead, given the colors and values for the extreme ends.
+  if( !sum( is.na( as.numeric( vsColorValues[ which( !is.na( vsColorValues ) ) ] ) ) ) )
+  {
+    vdNumericColors = as.numeric( vsColorValues )
+    vdNumericColors = vdNumericColors[ which( !is.na( vdNumericColors ) ) ]
+    vdSortedNumericColors = sort( vdNumericColors )
+    sLegendText = c( paste( sColorBy, vdSortedNumericColors[ 1 ], sep="_" ), 
+                     paste( sColorBy, vdSortedNumericColors[ length(vdSortedNumericColors) ], sep="_" ),
+                     paste( sMetadataShape, vsShapeValues, sep="_" ) )
+    sLegendColors = c(vsColorRBG[ which( vdNumericColors == vdSortedNumericColors[ 1 ] )[ 1 ] ],
+                      vsColorRBG[ which( vdNumericColors == vdSortedNumericColors[ length( vdSortedNumericColors ) ] )[ 1 ] ],
+                      rep(cDefaultColor,length(vsShapeValues)))
+  }
+  sLegendShapes = c( rep( cDefaultShape, length( sLegendText ) - length( vsShapeShapes ) ), vsShapeShapes )
+
+  # If any legend text was constructed then make the legend.
+  if( length( sLegendText ) >0 )
+  {
+    legend( "topright", legend = sLegendText, pch = sLegendShapes, col = sLegendColors )
+  }
+
+  # Original biplot call if you want to check the custom plotting of the script
+  # There will be one difference where the biplot call scales an axis, this one does not. In relation to the axes, the points, text and arrows should still match.
+  # Axes to the top and right are for the arrow, others are for markers and bug names.
+  #biplot(mNMDSData$points,mMetadataCoordinates[viMetadataDummy,],xlabs=vsShapes,xlab=paste("MDS1","Stress=",mNMDSData$stress),main="Biplot function Bugs and Sampes - Metadata Plotted with Centroids")
+  dev.off()
+}
+
+# This is the equivalent of __name__ == "__main__" in Python.
+# That is, if it's true we're being called as a command line script;
+# if it's false, we're being sourced or otherwise included, such as for
+# library or inlinedocs.
+if( identical( environment( ), globalenv( ) ) &&
+	!length( grep( "^source\\(", sys.calls( ) ) ) )
+{
+  lsArgs <- parse_args( pArgs, positional_arguments=TRUE )
+
+  funcDoBiplot(
+    sBugs = lsArgs$options$sBugs,
+    sMetadata = lsArgs$options$sMetadata,
+    sColorBy = lsArgs$options$sColorBy,
+    sColorRange = lsArgs$options$sColorRange,
+    sTextColor = lsArgs$options$sTextColor,
+    sArrowColor = lsArgs$options$sArrowColor,
+    sArrowTextColor = lsArgs$options$sArrowTextColor,
+    sPlotNAColor = lsArgs$options$sPlotNAColor,
+    sShapeBy = lsArgs$options$sShapeBy,
+    sShapes = lsArgs$options$sShapes,
+    sDefaultMarker = lsArgs$options$sDefaultMarker,
+    sRotateByMetadata = lsArgs$options$sRotateByMetadata,
+    dResizeArrow = lsArgs$options$dResizeArrow,
+    sTitle = lsArgs$options$sTitle,
+    sInputFileName = lsArgs$args[2],
+    sLastMetadata = lsArgs$args[1],
+    sOutputFileName = lsArgs$options$sOutputFileName)
+}