| 
0
 | 
     1 #!/usr/bin/env Rscript
 | 
| 
 | 
     2 
 | 
| 
 | 
     3 library(vegan)
 | 
| 
 | 
     4 library(optparse)
 | 
| 
 | 
     5 
 | 
| 
 | 
     6 funcGetCentroidForMetadatum <- function(
 | 
| 
 | 
     7 ### Given a binary metadatum, calculate the centroid of the samples associated with the metadata value of 1
 | 
| 
 | 
     8 # 1. Get all samples that have the metadata value of 1
 | 
| 
 | 
     9 # 2. Get the x and y coordinates of the selected samples
 | 
| 
 | 
    10 # 3. Get the median value for the x and ys
 | 
| 
 | 
    11 # 4. Return those coordinates as the centroid's X and Y value
 | 
| 
 | 
    12 vfMetadata,
 | 
| 
 | 
    13 ### Logical or integer (0,1) vector, TRUE or 1 values indicate correspoinding samples in the
 | 
| 
 | 
    14 ### mSamplePoints which will be used to define the centroid
 | 
| 
 | 
    15 mSamplePoints
 | 
| 
 | 
    16 ### Coordinates (columns;n=2) of samples (rows) corresponding to the vfMetadata
 | 
| 
 | 
    17 ){
 | 
| 
 | 
    18   # Check the lengths which should be equal
 | 
| 
 | 
    19   if(length(vfMetadata)!=nrow(mSamplePoints))
 | 
| 
 | 
    20   {
 | 
| 
 | 
    21     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=""))
 | 
| 
 | 
    22     return( FALSE )
 | 
| 
 | 
    23   }
 | 
| 
 | 
    24 
 | 
| 
 | 
    25   # Get all the samples that have the metadata value of 1
 | 
| 
 | 
    26   viMetadataSamples = which(as.integer(vfMetadata)==1)
 | 
| 
 | 
    27 
 | 
| 
 | 
    28   # Get the x and y coordinates for the selected samples
 | 
| 
 | 
    29   mSelectedPoints = mSamplePoints[viMetadataSamples,]
 | 
| 
 | 
    30 
 | 
| 
 | 
    31   # Get the median value for the x and the ys
 | 
| 
 | 
    32   if(!is.null(nrow(mSelectedPoints)))
 | 
| 
 | 
    33   {
 | 
| 
 | 
    34     return( list(x=median(mSelectedPoints[,1],na.rm = TRUE),y=median(mSelectedPoints[,2],na.rm = TRUE)) )
 | 
| 
 | 
    35   } else {
 | 
| 
 | 
    36     return( list(x=mSelectedPoints[1],y=mSelectedPoints[2]) )
 | 
| 
 | 
    37   }
 | 
| 
 | 
    38 }
 | 
| 
 | 
    39 
 | 
| 
 | 
    40 funcGetMaximumForMetadatum <- function(
 | 
| 
 | 
    41 ### Given a continuous metadata
 | 
| 
 | 
    42 ### 1. Use the x and ys from mSamplePoints for coordinates and the metadata value as a height (z)
 | 
| 
 | 
    43 ### 2. Use lowess to smooth the landscape
 | 
| 
 | 
    44 ### 3. Take the maximum of the landscape
 | 
| 
 | 
    45 ### 4. Return the coordiantes for the maximum as the centroid
 | 
| 
 | 
    46 vdMetadata,
 | 
| 
 | 
    47 ### Continuous (numeric or integer) metadata
 | 
| 
 | 
    48 mSamplePoints
 | 
| 
 | 
    49 ### Coordinates (columns;n=2) of samples (rows) corresponding to the vfMetadata
 | 
| 
 | 
    50 ){
 | 
| 
 | 
    51   # Work with data frame
 | 
| 
 | 
    52   if(class(mSamplePoints)=="matrix")
 | 
| 
 | 
    53   {
 | 
| 
 | 
    54     mSamplePoints = data.frame(mSamplePoints)
 | 
| 
 | 
    55   }
 | 
| 
 | 
    56   # Check the lengths of the dataframes and the metadata
 | 
| 
 | 
    57   if(length(vdMetadata)!=nrow(mSamplePoints))
 | 
| 
 | 
    58   {
 | 
| 
 | 
    59     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=""))
 | 
| 
 | 
    60     return( FALSE )
 | 
| 
 | 
    61   }
 | 
| 
 | 
    62 
 | 
| 
 | 
    63   # Add the metadata value to the points
 | 
| 
 | 
    64   mSamplePoints[3] = vdMetadata
 | 
| 
 | 
    65   names(mSamplePoints) = c("x","y","z") 
 | 
| 
 | 
    66 
 | 
| 
 | 
    67   # Create lowess to smooth the surface
 | 
| 
 | 
    68   # And calculate the fitted heights
 | 
| 
 | 
    69   # x = sample coordinate 1
 | 
| 
 | 
    70   # y = sample coordinate 2
 | 
| 
 | 
    71   # z = metadata value
 | 
| 
 | 
    72   loessSamples = loess(z~x*y, data=mSamplePoints, degree = 1, normalize = FALSE, na.action=na.omit)
 | 
| 
 | 
    73 
 | 
| 
 | 
    74   # Naively get the max
 | 
| 
 | 
    75   vdCoordinates = loessSamples$x[which(loessSamples$y==max(loessSamples$y)),]
 | 
| 
 | 
    76   return(list(lsmod = loessSamples, x=vdCoordinates[1],y=vdCoordinates[2]))
 | 
| 
 | 
    77 }
 | 
| 
 | 
    78 
 | 
| 
 | 
    79 funcMakeShapes <- function(
 | 
| 
 | 
    80 ### Takes care of defining shapes for the plot
 | 
| 
 | 
    81 dfInput,
 | 
| 
 | 
    82 ### Data frame of metadata measurements
 | 
| 
 | 
    83 sShapeBy,
 | 
| 
 | 
    84 ### The metadata to shape by
 | 
| 
 | 
    85 sShapes,
 | 
| 
 | 
    86 ### List of custom metadata (per level if factor).
 | 
| 
 | 
    87 ### 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 
 | 
| 
 | 
    88 cDefaultShape
 | 
| 
 | 
    89 ### Shape to default to if custom shapes are not used
 | 
| 
 | 
    90 ){
 | 
| 
 | 
    91   lShapes = list()
 | 
| 
 | 
    92   vsShapeValues = c()
 | 
| 
 | 
    93   vsShapeShapes = c()
 | 
| 
 | 
    94   vsShapes = c()
 | 
| 
 | 
    95   sMetadataId = sShapeBy
 | 
| 
 | 
    96 
 | 
| 
 | 
    97   # Set default shape, color, and color ranges 
 | 
| 
 | 
    98   if(!is.null(cDefaultShape))
 | 
| 
 | 
    99   {
 | 
| 
 | 
   100     # Default shape should be an int for the int pch options
 | 
| 
 | 
   101     if(!is.na(as.integer(cDefaultShape)))
 | 
| 
 | 
   102     {
 | 
| 
 | 
   103       cDefaultShape = as.integer(cDefaultShape)
 | 
| 
 | 
   104     }
 | 
| 
 | 
   105   } else {
 | 
| 
 | 
   106     cDefaultShape = 16
 | 
| 
 | 
   107   }
 | 
| 
 | 
   108 
 | 
| 
 | 
   109   # Make shapes
 | 
| 
 | 
   110   vsShapes = rep(cDefaultShape,nrow(dfInput))
 | 
| 
 | 
   111 
 | 
| 
 | 
   112   if(!is.null(sMetadataId))
 | 
| 
 | 
   113   {
 | 
| 
 | 
   114     if(is.null(sShapes))
 | 
| 
 | 
   115     {
 | 
| 
 | 
   116       vsShapeValues = unique(dfInput[[sMetadataId]])
 | 
| 
 | 
   117       vsShapeShapes = 1:length(vsShapeValues)
 | 
| 
 | 
   118     } else {
 | 
| 
 | 
   119       # Put the markers in the order of the values)
 | 
| 
 | 
   120       vsShapeBy = unlist(strsplit(sShapes,","))
 | 
| 
 | 
   121       for(sShapeBy in vsShapeBy)
 | 
| 
 | 
   122       {
 | 
| 
 | 
   123         vsShapeByPieces = unlist(strsplit(sShapeBy,":"))
 | 
| 
 | 
   124         lShapes[vsShapeByPieces[1]] = as.integer(vsShapeByPieces[2])
 | 
| 
 | 
   125       }
 | 
| 
 | 
   126       vsShapeValues = names(lShapes)
 | 
| 
 | 
   127    }
 | 
| 
 | 
   128 
 | 
| 
 | 
   129     # Shapes in the correct order
 | 
| 
 | 
   130     if(!is.null(sShapes))
 | 
| 
 | 
   131     {
 | 
| 
 | 
   132       vsShapeShapes = unlist(lapply(vsShapeValues,function(x) lShapes[[x]]))
 | 
| 
 | 
   133     }
 | 
| 
 | 
   134     vsShapeValues = paste(vsShapeValues)
 | 
| 
 | 
   135 
 | 
| 
 | 
   136     # Make the list of shapes
 | 
| 
 | 
   137     for(iShape in 1:length(vsShapeValues))
 | 
| 
 | 
   138     {
 | 
| 
 | 
   139       vsShapes[which(paste(dfInput[[sMetadataId]])==vsShapeValues[iShape])]=vsShapeShapes[iShape]
 | 
| 
 | 
   140     }
 | 
| 
 | 
   141 
 | 
| 
 | 
   142     # If they are all numeric characters, make numeric
 | 
| 
 | 
   143     viIntNas = which(is.na(as.integer(vsShapes)))
 | 
| 
 | 
   144     viNas = which(is.na(vsShapes))
 | 
| 
 | 
   145     if(length(setdiff(viIntNas,viNas))==0)
 | 
| 
 | 
   146     {
 | 
| 
 | 
   147       vsShapes = as.integer(vsShapes)
 | 
| 
 | 
   148     } else {
 | 
| 
 | 
   149       print("funcMakeShapes::Error: Please supply numbers 1-25 for shape in the -y,--shapeBy option")
 | 
| 
 | 
   150       vsShapeValues = c()
 | 
| 
 | 
   151       vsShapeShapes = c()
 | 
| 
 | 
   152     }
 | 
| 
 | 
   153   }
 | 
| 
 | 
   154   return(list(PlotShapes=vsShapes,Values=vsShapeValues,Shapes=vsShapeShapes,ID=sMetadataId,DefaultShape=cDefaultShape))
 | 
| 
 | 
   155 }
 | 
| 
 | 
   156 
 | 
| 
 | 
   157 ### Global defaults
 | 
| 
 | 
   158 c_sDefaultColorBy = NULL
 | 
| 
 | 
   159 c_sDefaultColorRange = "orange,cyan"
 | 
| 
 | 
   160 c_sDefaultTextColor = "black"
 | 
| 
 | 
   161 c_sDefaultArrowColor = "cyan"
 | 
| 
 | 
   162 c_sDefaultArrowTextColor = "Blue"
 | 
| 
 | 
   163 c_sDefaultNAColor = "grey"
 | 
| 
 | 
   164 c_sDefaultShapeBy = NULL
 | 
| 
 | 
   165 c_sDefaultShapes = NULL
 | 
| 
 | 
   166 c_sDefaultMarker = "16"
 | 
| 
 | 
   167 c_sDefaultRotateByMetadata = NULL
 | 
| 
 | 
   168 c_sDefaultResizeArrow = 1
 | 
| 
 | 
   169 c_sDefaultTitle = "Custom Biplot of Bugs and Samples - Metadata Plotted with Centroids"
 | 
| 
 | 
   170 c_sDefaultOutputFile = NULL
 | 
| 
 | 
   171 
 | 
| 
 | 
   172 ### Create command line argument parser
 | 
| 
 | 
   173 pArgs <- OptionParser( usage = "%prog last_metadata input.tsv" )
 | 
| 
 | 
   174 
 | 
| 
 | 
   175 # Selecting features to plot
 | 
| 
 | 
   176 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")
 | 
| 
 | 
   177 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")
 | 
| 
 | 
   178 
 | 
| 
 | 
   179 # Colors
 | 
| 
 | 
   180 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.")
 | 
| 
 | 
   181 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))
 | 
| 
 | 
   182 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))
 | 
| 
 | 
   183 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))
 | 
| 
 | 
   184 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))
 | 
| 
 | 
   185 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))
 | 
| 
 | 
   186 
 | 
| 
 | 
   187 # Shapes
 | 
| 
 | 
   188 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")
 | 
| 
 | 
   189 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.")
 | 
| 
 | 
   190 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.")
 | 
| 
 | 
   191 
 | 
| 
 | 
   192 # Plot manipulations
 | 
| 
 | 
   193 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 .")
 | 
| 
 | 
   194 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.")
 | 
| 
 | 
   195 
 | 
| 
 | 
   196 # Misc
 | 
| 
 | 
   197 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.")
 | 
| 
 | 
   198 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.")
 | 
| 
 | 
   199 
 | 
| 
 | 
   200 funcDoBiplot <- function(
 | 
| 
 | 
   201 ### 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.
 | 
| 
 | 
   202 sBugs,
 | 
| 
 | 
   203 ### Comma delimited list of data to plot as text. Bug|1,Bug|2
 | 
| 
 | 
   204 sMetadata,
 | 
| 
 | 
   205 ### Comma delimited list of metadata to plot as arrows. metadata1,metadata2,metadata3.
 | 
| 
 | 
   206 sColorBy = c_sDefaultColorBy,
 | 
| 
 | 
   207 ### The id of the metadatum to use to make the marker colors. Expected to be a continuous metadata.
 | 
| 
 | 
   208 sColorRange = c_sDefaultColorRange,
 | 
| 
 | 
   209 ### Colors used to color the samples; a gradient will be formed between the color. Example orange,cyan
 | 
| 
 | 
   210 sTextColor = c_sDefaultTextColor,
 | 
| 
 | 
   211 ### The color bug features will be plotted with as text. Example black
 | 
| 
 | 
   212 sArrowColor = c_sDefaultArrowColor,
 | 
| 
 | 
   213 ### The color metadata features will be plotted with as an arrow and text. Example cyan
 | 
| 
 | 
   214 sArrowTextColor = c_sDefaultArrowTextColor,
 | 
| 
 | 
   215 ### The color for the metadata text ploted by the head of the metadata arrow. Example Blue
 | 
| 
 | 
   216 sPlotNAColor = c_sDefaultNAColor,
 | 
| 
 | 
   217 ### Plot NA values as this color. Example grey
 | 
| 
 | 
   218 sShapeBy = c_sDefaultShapeBy,
 | 
| 
 | 
   219 ### The metadata to use to make marker shapes. Expected to be a discrete metadatum.
 | 
| 
 | 
   220 sShapes = c_sDefaultShapes,
 | 
| 
 | 
   221 ### 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.
 | 
| 
 | 
   222 sDefaultMarker = c_sDefaultMarker,
 | 
| 
 | 
   223 ### The default marker shape to use if shapes are not otherwise indicated.
 | 
| 
 | 
   224 sRotateByMetadata = c_sDefaultRotateByMetadata,
 | 
| 
 | 
   225 ### Metadata and value to rotate by. example Environment_HighLumninosity,100
 | 
| 
 | 
   226 dResizeArrow = c_sDefaultResizeArrow,
 | 
| 
 | 
   227 ### Scale factor to resize tthe metadata arrows
 | 
| 
 | 
   228 sTitle = c_sDefaultTitle,
 | 
| 
 | 
   229 ### The title for the figure.
 | 
| 
 | 
   230 sInputFileName,
 | 
| 
 | 
   231 ### File to input (tsv file: tab separated, row = sample file)
 | 
| 
 | 
   232 sLastMetadata,
 | 
| 
 | 
   233 ### Last metadata that seperates data and metadata
 | 
| 
 | 
   234 sOutputFileName = c_sDefaultOutputFile
 | 
| 
 | 
   235 ### The file name to save the figure.
 | 
| 
 | 
   236 ){
 | 
| 
 | 
   237   print("IN Biplot")
 | 
| 
 | 
   238   # Define the colors
 | 
| 
 | 
   239   vsColorRange = c("blue","orange")
 | 
| 
 | 
   240   cDefaultColor = "black"
 | 
| 
 | 
   241   if(!is.null(sColorRange))
 | 
| 
 | 
   242   {
 | 
| 
 | 
   243     vsColorRange = unlist(strsplit(sColorRange,","))
 | 
| 
 | 
   244   }
 | 
| 
 | 
   245 
 | 
| 
 | 
   246   # List of bugs to plot
 | 
| 
 | 
   247   # If there is a list it needs to be more than one.
 | 
| 
 | 
   248   vsBugsToPlot = c()
 | 
| 
 | 
   249   if(!is.null(sBugs))
 | 
| 
 | 
   250   {
 | 
| 
 | 
   251     vsBugsToPlot = unlist(strsplit(sBugs,","))
 | 
| 
 | 
   252   }
 | 
| 
 | 
   253 
 | 
| 
 | 
   254   print("vsBugsToPlot")
 | 
| 
 | 
   255   print(vsBugsToPlot)
 | 
| 
 | 
   256   # Metadata to plot
 | 
| 
 | 
   257   vsMetadata = c()
 | 
| 
 | 
   258   if(!is.null(sMetadata))
 | 
| 
 | 
   259   {
 | 
| 
 | 
   260     vsMetadata = unlist(strsplit(sMetadata,","))
 | 
| 
 | 
   261   }
 | 
| 
 | 
   262 
 | 
| 
 | 
   263   print("vsMetadata")
 | 
| 
 | 
   264   print(vsMetadata)
 | 
| 
 | 
   265   ### Load table
 | 
| 
 | 
   266   if(class(sInputFileName)=="character")
 | 
| 
 | 
   267   {
 | 
| 
 | 
   268     dfInput = read.table(sInputFileName, sep = "\t", header=TRUE)
 | 
| 
 | 
   269     names(dfInput) = unlist(lapply(names(dfInput),function(x) gsub(".","|",x,fixed=TRUE)))
 | 
| 
 | 
   270     row.names(dfInput) = dfInput[,1]
 | 
| 
 | 
   271     dfInput = dfInput[-1]
 | 
| 
 | 
   272   } else {dfInput = sInputFileName}
 | 
| 
 | 
   273 
 | 
| 
 | 
   274   ### Get positions of all metadata or all data
 | 
| 
 | 
   275   iLastMetadata = which(names(dfInput)==sLastMetadata)
 | 
| 
 | 
   276   viMetadata = 1:iLastMetadata
 | 
| 
 | 
   277   viData = (iLastMetadata+1):ncol(dfInput)
 | 
| 
 | 
   278 
 | 
| 
 | 
   279   ### Dummy the metadata if discontinuous
 | 
| 
 | 
   280   ### Leave the continous metadata alone but include
 | 
| 
 | 
   281   listMetadata = list()
 | 
| 
 | 
   282   vsRowNames = c()
 | 
| 
 | 
   283   viContinuousMetadata = c()
 | 
| 
 | 
   284   for(i in viMetadata)
 | 
| 
 | 
   285   {
 | 
| 
 | 
   286     print( names( dfInput )[i] )
 | 
| 
 | 
   287     vCurMetadata = unlist(dfInput[i])
 | 
| 
 | 
   288     if( ( is.numeric(vCurMetadata)||is.integer(vCurMetadata) )  && ( length( unique( vCurMetadata ) ) >= c_iNonFactorLevelThreshold ) )
 | 
| 
 | 
   289     {
 | 
| 
 | 
   290       vCurMetadata[which(is.na(vCurMetadata))] = mean(vCurMetadata,na.rm=TRUE)
 | 
| 
 | 
   291       listMetadata[[length(listMetadata)+1]] = vCurMetadata
 | 
| 
 | 
   292       vsRowNames = c(vsRowNames,names(dfInput)[i])
 | 
| 
 | 
   293       viContinuousMetadata = c(viContinuousMetadata,length(listMetadata))
 | 
| 
 | 
   294     } else {
 | 
| 
 | 
   295       vCurMetadata = as.factor(vCurMetadata)
 | 
| 
 | 
   296       vsLevels = levels(vCurMetadata)
 | 
| 
 | 
   297       for(sLevel in vsLevels)
 | 
| 
 | 
   298       { 
 | 
| 
 | 
   299         vNewMetadata = rep(0,length(vCurMetadata))
 | 
| 
 | 
   300         vNewMetadata[which(vCurMetadata == sLevel)] = 1
 | 
| 
 | 
   301         listMetadata[[length(listMetadata)+1]] = vNewMetadata
 | 
| 
 | 
   302         vsRowNames = c(vsRowNames,paste(names(dfInput)[i],sLevel,sep="_"))
 | 
| 
 | 
   303       }
 | 
| 
 | 
   304     }
 | 
| 
 | 
   305   }
 | 
| 
 | 
   306 
 | 
| 
 | 
   307   # Convert to data frame
 | 
| 
 | 
   308   dfDummyMetadata = as.data.frame(sapply(listMetadata,rbind))
 | 
| 
 | 
   309   names(dfDummyMetadata) = vsRowNames
 | 
| 
 | 
   310   iNumberMetadata = ncol(dfDummyMetadata)
 | 
| 
 | 
   311 
 | 
| 
 | 
   312   # Data to use in ordination in NMDS
 | 
| 
 | 
   313   # All cleaned bug data
 | 
| 
 | 
   314   dfData = dfInput[viData]
 | 
| 
 | 
   315 
 | 
| 
 | 
   316   # If rotating the ordination by a metadata
 | 
| 
 | 
   317   # 1. Add in the metadata as a bug
 | 
| 
 | 
   318   # 2. Multiply the bug by the weight
 | 
| 
 | 
   319   # 3. Push this through the NMDS
 | 
| 
 | 
   320   if(!is.null(sRotateByMetadata))
 | 
| 
 | 
   321   {
 | 
| 
 | 
   322     vsRotateMetadata = unlist(strsplit(sRotateByMetadata,","))
 | 
| 
 | 
   323     sMetadata = vsRotateMetadata[1]
 | 
| 
 | 
   324     dWeight = as.numeric(vsRotateMetadata[2])
 | 
| 
 | 
   325     sOrdinationMetadata = dfDummyMetadata[sMetadata]*dWeight
 | 
| 
 | 
   326     dfData[sMetadata] = sOrdinationMetadata
 | 
| 
 | 
   327   }
 | 
| 
 | 
   328 
 | 
| 
 | 
   329   # Run NMDS on bug data (Default B-C)
 | 
| 
 | 
   330   # Will have species and points because working off of raw data
 | 
| 
 | 
   331   mNMDSData = metaMDS(dfData,k=2)
 | 
| 
 | 
   332 
 | 
| 
 | 
   333   ## Make shapes
 | 
| 
 | 
   334   # Defines the shapes and the metadata they are based on
 | 
| 
 | 
   335   # Metadata to use as shapes
 | 
| 
 | 
   336   lShapeInfo = funcMakeShapes(dfInput=dfInput, sShapeBy=sShapeBy, sShapes=sShapes, cDefaultShape=sDefaultMarker)
 | 
| 
 | 
   337 
 | 
| 
 | 
   338   sMetadataShape = lShapeInfo[["ID"]]
 | 
| 
 | 
   339   vsShapeValues = lShapeInfo[["Values"]]
 | 
| 
 | 
   340   vsShapeShapes = lShapeInfo[["Shapes"]]
 | 
| 
 | 
   341   vsShapes = lShapeInfo[["PlotShapes"]]
 | 
| 
 | 
   342   cDefaultShape = lShapeInfo[["DefaultShape"]]
 | 
| 
 | 
   343 
 | 
| 
 | 
   344   # Colors
 | 
| 
 | 
   345   vsColors = rep(cDefaultColor,nrow(dfInput))
 | 
| 
 | 
   346   vsColorValues = c()
 | 
| 
 | 
   347   vsColorRBG = c()
 | 
| 
 | 
   348   if(!is.null(sColorBy))
 | 
| 
 | 
   349   {
 | 
| 
 | 
   350     vsColorValues = paste(sort(unique(unlist(dfInput[[sColorBy]])),na.last=TRUE))
 | 
| 
 | 
   351     iLengthColorValues = length(vsColorValues)
 | 
| 
 | 
   352 
 | 
| 
 | 
   353     vsColorRBG = lapply(1:iLengthColorValues/iLengthColorValues,colorRamp(vsColorRange))
 | 
| 
 | 
   354     vsColorRBG = unlist(lapply(vsColorRBG, function(x) rgb(x[1]/255,x[2]/255,x[3]/255)))
 | 
| 
 | 
   355 
 | 
| 
 | 
   356     for(iColor in 1:length(vsColorRBG))
 | 
| 
 | 
   357     {
 | 
| 
 | 
   358       vsColors[which(paste(dfInput[[sColorBy]])==vsColorValues[iColor])]=vsColorRBG[iColor]
 | 
| 
 | 
   359     }
 | 
| 
 | 
   360 
 | 
| 
 | 
   361     #If NAs are seperately given color, then color here
 | 
| 
 | 
   362     if(!is.null(sPlotNAColor))
 | 
| 
 | 
   363     {
 | 
| 
 | 
   364       vsColors[which(is.na(dfInput[[sColorBy]]))] = sPlotNAColor
 | 
| 
 | 
   365       vsColorRBG[which(vsColorValues=="NA")] = sPlotNAColor
 | 
| 
 | 
   366     }
 | 
| 
 | 
   367   }
 | 
| 
 | 
   368 
 | 
| 
 | 
   369   print("names(dfDummyMetadata)")
 | 
| 
 | 
   370   print(names(dfDummyMetadata))
 | 
| 
 | 
   371 
 | 
| 
 | 
   372   # Reduce the bugs down to the ones in the list to be plotted
 | 
| 
 | 
   373   viBugsToPlot = which(row.names(mNMDSData$species) %in% vsBugsToPlot)
 | 
| 
 | 
   374   viMetadataDummy = which(names(dfDummyMetadata) %in% vsMetadata)
 | 
| 
 | 
   375 
 | 
| 
 | 
   376   print("viBugsToPlot")
 | 
| 
 | 
   377   print(viBugsToPlot)
 | 
| 
 | 
   378   print("viMetadataDummy")
 | 
| 
 | 
   379   print(names(dfDummyMetadata)[viMetadataDummy])
 | 
| 
 | 
   380 
 | 
| 
 | 
   381   # Build the matrix of metadata coordinates
 | 
| 
 | 
   382   mMetadataCoordinates = matrix(rep(NA, iNumberMetadata*2),nrow=iNumberMetadata)
 | 
| 
 | 
   383   for( i in 1:iNumberMetadata )
 | 
| 
 | 
   384   {
 | 
| 
 | 
   385     lxReturn = NA
 | 
| 
 | 
   386     if( i %in% viContinuousMetadata )
 | 
| 
 | 
   387     {
 | 
| 
 | 
   388       lxReturn = funcGetMaximumForMetadatum(dfDummyMetadata[[i]],mNMDSData$points)
 | 
| 
 | 
   389     } else {
 | 
| 
 | 
   390       lxReturn = funcGetCentroidForMetadatum(dfDummyMetadata[[i]],mNMDSData$points)
 | 
| 
 | 
   391     }
 | 
| 
 | 
   392     mMetadataCoordinates[i,] = c(lxReturn$x,lxReturn$y)
 | 
| 
 | 
   393   }
 | 
| 
 | 
   394   row.names(mMetadataCoordinates) = vsRowNames
 | 
| 
 | 
   395 
 | 
| 
 | 
   396   # Plot the biplot with the centroid constructed metadata coordinates
 | 
| 
 | 
   397   if(length(viMetadataDummy)==0)
 | 
| 
 | 
   398   {
 | 
| 
 | 
   399     viMetadataDummy = 1:nrow(mMetadataCoordinates)
 | 
| 
 | 
   400   }
 | 
| 
 | 
   401 
 | 
| 
 | 
   402   # Plot samples
 | 
| 
 | 
   403   # Make output name
 | 
| 
 | 
   404   if(is.null(sOutputFileName))
 | 
| 
 | 
   405   {
 | 
| 
 | 
   406     viPeriods = which(sInputFileName==".")
 | 
| 
 | 
   407     if(length(viPeriods)>0)
 | 
| 
 | 
   408     {
 | 
| 
 | 
   409       sOutputFileName = paste(OutputFileName[1:viPeriods[length(viPeriods)]],"pdf",sep=".")
 | 
| 
 | 
   410     } else {
 | 
| 
 | 
   411       sOutputFileName = paste(sInputFileName,"pdf",sep=".")
 | 
| 
 | 
   412     }
 | 
| 
 | 
   413   }
 | 
| 
 | 
   414 
 | 
| 
 | 
   415   pdf(sOutputFileName, useDingbats=FALSE)
 | 
| 
 | 
   416   plot(mNMDSData$points, xlab=paste("NMDS1","Stress=",mNMDSData$stress), ylab="NMDS2", pch=vsShapes, col=vsColors)
 | 
| 
 | 
   417   title(sTitle,line=3)
 | 
| 
 | 
   418 
 | 
| 
 | 
   419   # Plot Bugs
 | 
| 
 | 
   420   mPlotBugs = mNMDSData$species[viBugsToPlot,]
 | 
| 
 | 
   421   if(length(viBugsToPlot)==1)
 | 
| 
 | 
   422   {
 | 
| 
 | 
   423     text(x=mPlotBugs[1],y=mPlotBugs[2],labels=row.names(mNMDSData$species)[viBugsToPlot],col=sTextColor)
 | 
| 
 | 
   424   } else if(length(viBugsToPlot)>1){
 | 
| 
 | 
   425     text(x=mPlotBugs[,1],y=mPlotBugs[,2],labels=row.names(mNMDSData$species)[viBugsToPlot],col=sTextColor)
 | 
| 
 | 
   426   }
 | 
| 
 | 
   427 
 | 
| 
 | 
   428   # Add alternative axes
 | 
| 
 | 
   429   axis(3, col=sArrowColor)
 | 
| 
 | 
   430   axis(4, col=sArrowColor)
 | 
| 
 | 
   431   box(col = "black")
 | 
| 
 | 
   432 
 | 
| 
 | 
   433   # Plot Metadata
 | 
| 
 | 
   434   if(length(viMetadataDummy)>0)
 | 
| 
 | 
   435   {
 | 
| 
 | 
   436     for(i in viMetadataDummy)
 | 
| 
 | 
   437     {
 | 
| 
 | 
   438       curCoordinates = mMetadataCoordinates[i,]
 | 
| 
 | 
   439       curCoordinates = curCoordinates * dResizeArrow
 | 
| 
 | 
   440       # Plot Arrow
 | 
| 
 | 
   441       arrows(0,0, curCoordinates[1] * 0.8, curCoordinates[2] * 0.8, col=sArrowColor, length=0.1 )
 | 
| 
 | 
   442     }
 | 
| 
 | 
   443     # Plot text
 | 
| 
 | 
   444     if(length(viMetadataDummy)==1)
 | 
| 
 | 
   445     {
 | 
| 
 | 
   446       text(x=mMetadataCoordinates[viMetadataDummy,][1]*dResizeArrow*0.8, y=mMetadataCoordinates[viMetadataDummy,][2]*dResizeArrow*0.8, labels=row.names(mMetadataCoordinates)[viMetadataDummy],col=sArrowTextColor)
 | 
| 
 | 
   447     } else {
 | 
| 
 | 
   448       text(x=mMetadataCoordinates[viMetadataDummy,1]*dResizeArrow*0.8, y=mMetadataCoordinates[viMetadataDummy,2]*dResizeArrow*0.8, labels=row.names(mMetadataCoordinates)[viMetadataDummy],col=sArrowTextColor)
 | 
| 
 | 
   449     }
 | 
| 
 | 
   450   }
 | 
| 
 | 
   451 
 | 
| 
 | 
   452   # Create Legend
 | 
| 
 | 
   453   # The text default is the colorMetadata_level (one per level) plus the ShapeMetadata_level (one per level)
 | 
| 
 | 
   454   # The color default is already determined colors plus grey for shapes.
 | 
| 
 | 
   455   sLegendText = c(paste(vsColorValues,sColorBy,sep="_"),paste(sMetadataShape,vsShapeValues,sep="_"))
 | 
| 
 | 
   456   sLegendColors = c(vsColorRBG,rep(cDefaultColor,length(vsShapeValues)))
 | 
| 
 | 
   457 
 | 
| 
 | 
   458   # If the color values are numeric
 | 
| 
 | 
   459   # Too many values may be given in the legend (given they may be a continuous range of values)
 | 
| 
 | 
   460   # To reduce this they are summarized instead, given the colors and values for the extreme ends.
 | 
| 
 | 
   461   if( !sum( is.na( as.numeric( vsColorValues[ which( !is.na( vsColorValues ) ) ] ) ) ) )
 | 
| 
 | 
   462   {
 | 
| 
 | 
   463     vdNumericColors = as.numeric( vsColorValues )
 | 
| 
 | 
   464     vdNumericColors = vdNumericColors[ which( !is.na( vdNumericColors ) ) ]
 | 
| 
 | 
   465     vdSortedNumericColors = sort( vdNumericColors )
 | 
| 
 | 
   466     sLegendText = c( paste( sColorBy, vdSortedNumericColors[ 1 ], sep="_" ), 
 | 
| 
 | 
   467                      paste( sColorBy, vdSortedNumericColors[ length(vdSortedNumericColors) ], sep="_" ),
 | 
| 
 | 
   468                      paste( sMetadataShape, vsShapeValues, sep="_" ) )
 | 
| 
 | 
   469     sLegendColors = c(vsColorRBG[ which( vdNumericColors == vdSortedNumericColors[ 1 ] )[ 1 ] ],
 | 
| 
 | 
   470                       vsColorRBG[ which( vdNumericColors == vdSortedNumericColors[ length( vdSortedNumericColors ) ] )[ 1 ] ],
 | 
| 
 | 
   471                       rep(cDefaultColor,length(vsShapeValues)))
 | 
| 
 | 
   472   }
 | 
| 
 | 
   473   sLegendShapes = c( rep( cDefaultShape, length( sLegendText ) - length( vsShapeShapes ) ), vsShapeShapes )
 | 
| 
 | 
   474 
 | 
| 
 | 
   475   # If any legend text was constructed then make the legend.
 | 
| 
 | 
   476   if( length( sLegendText ) >0 )
 | 
| 
 | 
   477   {
 | 
| 
 | 
   478     legend( "topright", legend = sLegendText, pch = sLegendShapes, col = sLegendColors )
 | 
| 
 | 
   479   }
 | 
| 
 | 
   480 
 | 
| 
 | 
   481   # Original biplot call if you want to check the custom plotting of the script
 | 
| 
 | 
   482   # 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.
 | 
| 
 | 
   483   # Axes to the top and right are for the arrow, others are for markers and bug names.
 | 
| 
 | 
   484   #biplot(mNMDSData$points,mMetadataCoordinates[viMetadataDummy,],xlabs=vsShapes,xlab=paste("MDS1","Stress=",mNMDSData$stress),main="Biplot function Bugs and Sampes - Metadata Plotted with Centroids")
 | 
| 
 | 
   485   dev.off()
 | 
| 
 | 
   486 }
 | 
| 
 | 
   487 
 | 
| 
 | 
   488 # This is the equivalent of __name__ == "__main__" in Python.
 | 
| 
 | 
   489 # That is, if it's true we're being called as a command line script;
 | 
| 
 | 
   490 # if it's false, we're being sourced or otherwise included, such as for
 | 
| 
 | 
   491 # library or inlinedocs.
 | 
| 
 | 
   492 if( identical( environment( ), globalenv( ) ) &&
 | 
| 
 | 
   493 	!length( grep( "^source\\(", sys.calls( ) ) ) )
 | 
| 
 | 
   494 {
 | 
| 
 | 
   495   lsArgs <- parse_args( pArgs, positional_arguments=TRUE )
 | 
| 
 | 
   496 
 | 
| 
 | 
   497   funcDoBiplot(
 | 
| 
 | 
   498     sBugs = lsArgs$options$sBugs,
 | 
| 
 | 
   499     sMetadata = lsArgs$options$sMetadata,
 | 
| 
 | 
   500     sColorBy = lsArgs$options$sColorBy,
 | 
| 
 | 
   501     sColorRange = lsArgs$options$sColorRange,
 | 
| 
 | 
   502     sTextColor = lsArgs$options$sTextColor,
 | 
| 
 | 
   503     sArrowColor = lsArgs$options$sArrowColor,
 | 
| 
 | 
   504     sArrowTextColor = lsArgs$options$sArrowTextColor,
 | 
| 
 | 
   505     sPlotNAColor = lsArgs$options$sPlotNAColor,
 | 
| 
 | 
   506     sShapeBy = lsArgs$options$sShapeBy,
 | 
| 
 | 
   507     sShapes = lsArgs$options$sShapes,
 | 
| 
 | 
   508     sDefaultMarker = lsArgs$options$sDefaultMarker,
 | 
| 
 | 
   509     sRotateByMetadata = lsArgs$options$sRotateByMetadata,
 | 
| 
 | 
   510     dResizeArrow = lsArgs$options$dResizeArrow,
 | 
| 
 | 
   511     sTitle = lsArgs$options$sTitle,
 | 
| 
 | 
   512     sInputFileName = lsArgs$args[2],
 | 
| 
 | 
   513     sLastMetadata = lsArgs$args[1],
 | 
| 
 | 
   514     sOutputFileName = lsArgs$options$sOutputFileName)
 | 
| 
 | 
   515 }
 |