|
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 = NULL
|
|
|
164 c_sDefaultShapeBy = NULL
|
|
|
165 c_sDefaultShapes = NULL
|
|
|
166 c_sDefaultMarker = "16"
|
|
|
167 c_fDefaultPlotArrows = TRUE
|
|
|
168 c_sDefaultRotateByMetadata = NULL
|
|
|
169 c_sDefaultResizeArrow = 1
|
|
|
170 c_sDefaultTitle = "Custom Biplot of Bugs and Samples - Metadata Plotted with Centroids"
|
|
|
171 c_sDefaultOutputFile = NULL
|
|
|
172
|
|
|
173 ### Create command line argument parser
|
|
|
174 pArgs <- OptionParser( usage = "%prog last_metadata input.tsv" )
|
|
|
175
|
|
|
176 # Selecting features to plot
|
|
|
177 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. To not plot metadata pass a blank or empty space. Bug|1,Bug|2")
|
|
|
178 # The default needs to stay null for metadata or code needs to be changed in the plotting for a check to see if the metadata was default. Search for "#!# metadata default dependent"
|
|
|
179 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. To not plot metadata pass a blank or empty space. metadata1,metadata2,metadata3")
|
|
|
180
|
|
|
181 # Colors
|
|
|
182 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.")
|
|
|
183 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))
|
|
|
184 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))
|
|
|
185 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))
|
|
|
186 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))
|
|
|
187 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))
|
|
|
188
|
|
|
189 # Shapes
|
|
|
190 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")
|
|
|
191 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.")
|
|
|
192 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.")
|
|
|
193
|
|
|
194 # Plot manipulations
|
|
|
195 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 .")
|
|
|
196 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.")
|
|
|
197 pArgs <- add_option( pArgs, c("-A", "--noArrows"), type="logical", action="store_true", default=FALSE, dest="fNoPlotMetadataArrows", metavar="DoNotPlotArrows", help="Adding this flag allows one to plot metadata labels without the arrows.")
|
|
|
198
|
|
|
199 # Misc
|
|
|
200 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.")
|
|
|
201 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.")
|
|
|
202
|
|
|
203 funcDoBiplot <- function(
|
|
|
204 ### 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.
|
|
|
205 sBugs,
|
|
|
206 ### Comma delimited list of data to plot as text. Bug|1,Bug|2
|
|
|
207 sMetadata,
|
|
|
208 ### Comma delimited list of metadata to plot as arrows. metadata1,metadata2,metadata3.
|
|
|
209 sColorBy = c_sDefaultColorBy,
|
|
|
210 ### The id of the metadatum to use to make the marker colors. Expected to be a continuous metadata.
|
|
|
211 sColorRange = c_sDefaultColorRange,
|
|
|
212 ### Colors used to color the samples; a gradient will be formed between the color. Example orange,cyan
|
|
|
213 sTextColor = c_sDefaultTextColor,
|
|
|
214 ### The color bug features will be plotted with as text. Example black
|
|
|
215 sArrowColor = c_sDefaultArrowColor,
|
|
|
216 ### The color metadata features will be plotted with as an arrow and text. Example cyan
|
|
|
217 sArrowTextColor = c_sDefaultArrowTextColor,
|
|
|
218 ### The color for the metadata text ploted by the head of the metadat arrow. Example Blue
|
|
|
219 sPlotNAColor = c_sDefaultNAColor,
|
|
|
220 ### Plot NA values as this color. Example grey
|
|
|
221 sShapeBy = c_sDefaultShapeBy,
|
|
|
222 ### The metadata to use to make marker shapes. Expected to be a discrete metadatum.
|
|
|
223 sShapes = c_sDefaultShapes,
|
|
|
224 ### 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.
|
|
|
225 sDefaultMarker = c_sDefaultMarker,
|
|
|
226 ### The default marker shape to use if shapes are not otherwise indicated.
|
|
|
227 sRotateByMetadata = c_sDefaultRotateByMetadata,
|
|
|
228 ### Metadata and value to rotate by. example Environment_HighLumninosity,100
|
|
|
229 dResizeArrow = c_sDefaultResizeArrow,
|
|
|
230 ### Scale factor to resize tthe metadata arrows
|
|
|
231 fPlotArrow = c_fDefaultPlotArrows,
|
|
|
232 ### A flag which can be used to turn off arrow plotting.
|
|
|
233 sTitle = c_sDefaultTitle,
|
|
|
234 ### The title for the figure.
|
|
|
235 sInputFileName,
|
|
|
236 ### File to input (tsv file: tab seperated, row = sample file)
|
|
|
237 sLastMetadata,
|
|
|
238 ### Last metadata that seperates data and metadata
|
|
|
239 sOutputFileName = c_sDefaultOutputFile
|
|
|
240 ### The file name to save the figure.
|
|
|
241 ){
|
|
|
242 # Define the colors
|
|
|
243 vsColorRange = c("blue","orange")
|
|
|
244 if(!is.null(sColorRange))
|
|
|
245 {
|
|
|
246 vsColorRange = unlist(strsplit(sColorRange,","))
|
|
|
247 }
|
|
|
248 cDefaultColor = "grey"
|
|
|
249 if(!is.null(vsColorRange) && length(vsColorRange)>0)
|
|
|
250 {
|
|
|
251 cDefaultColor = vsColorRange[1]
|
|
|
252 }
|
|
|
253
|
|
|
254 # List of bugs to plot
|
|
|
255 # If there is a list it needs to be more than one.
|
|
|
256 vsBugsToPlot = c()
|
|
|
257 if(!is.null(sBugs))
|
|
|
258 {
|
|
|
259 vsBugsToPlot = unlist(strsplit(sBugs,","))
|
|
|
260 }
|
|
|
261 # Metadata to plot
|
|
|
262 vsMetadata = c()
|
|
|
263 if(!is.null(sMetadata))
|
|
|
264 {
|
|
|
265 vsMetadata = unlist(strsplit(sMetadata,","))
|
|
|
266 }
|
|
|
267
|
|
|
268 ### Load table
|
|
|
269 dfInput = sInputFileName
|
|
|
270 if(class(sInputFileName)=="character")
|
|
|
271 {
|
|
|
272 dfInput = read.table(sInputFileName, sep = "\t", header=TRUE)
|
|
|
273 names(dfInput) = unlist(lapply(names(dfInput),function(x) gsub(".","|",x,fixed=TRUE)))
|
|
|
274 row.names(dfInput) = dfInput[,1]
|
|
|
275 dfInput = dfInput[-1]
|
|
|
276 }
|
|
|
277
|
|
|
278 iLastMetadata = which(names(dfInput)==sLastMetadata)
|
|
|
279 viMetadata = 1:iLastMetadata
|
|
|
280 viData = (iLastMetadata+1):ncol(dfInput)
|
|
|
281
|
|
|
282 ### Dummy the metadata if discontinuous
|
|
|
283 ### Leave the continous metadata alone but include
|
|
|
284 listMetadata = list()
|
|
|
285 vsRowNames = c()
|
|
|
286 viContinuousMetadata = c()
|
|
|
287 for(i in viMetadata)
|
|
|
288 {
|
|
|
289 vCurMetadata = unlist(dfInput[i])
|
|
|
290 if(is.numeric(vCurMetadata)||is.integer(vCurMetadata))
|
|
|
291 {
|
|
|
292 vCurMetadata[which(is.na(vCurMetadata))] = mean(vCurMetadata,na.rm=TRUE)
|
|
|
293 listMetadata[[length(listMetadata)+1]] = vCurMetadata
|
|
|
294 vsRowNames = c(vsRowNames,names(dfInput)[i])
|
|
|
295 viContinuousMetadata = c(viContinuousMetadata,length(listMetadata))
|
|
|
296 } else {
|
|
|
297 vCurMetadata = as.factor(vCurMetadata)
|
|
|
298 vsLevels = levels(vCurMetadata)
|
|
|
299 for(sLevel in vsLevels)
|
|
|
300 {
|
|
|
301 vNewMetadata = rep(0,length(vCurMetadata))
|
|
|
302 vNewMetadata[which(vCurMetadata == sLevel)] = 1
|
|
|
303 listMetadata[[length(listMetadata)+1]] = vNewMetadata
|
|
|
304 vsRowNames = c(vsRowNames,paste(names(dfInput)[i],sLevel,sep="_"))
|
|
|
305 }
|
|
|
306 }
|
|
|
307 }
|
|
|
308
|
|
|
309 # Convert to data frame
|
|
|
310 dfDummyMetadata = as.data.frame(sapply(listMetadata,rbind))
|
|
|
311 names(dfDummyMetadata) = vsRowNames
|
|
|
312 iNumberMetadata = ncol(dfDummyMetadata)
|
|
|
313
|
|
|
314 # Data to use in ordination in NMDS
|
|
|
315 dfData = dfInput[viData]
|
|
|
316
|
|
|
317 # If rotating the ordination by a metadata
|
|
|
318 # 1. Add in the metadata as a bug
|
|
|
319 # 2. Multiply the bug by the weight
|
|
|
320 # 3. Push this through the NMDS
|
|
|
321 if(!is.null(sRotateByMetadata))
|
|
|
322 {
|
|
|
323 vsRotateMetadata = unlist(strsplit(sRotateByMetadata,","))
|
|
|
324 sMetadata = vsRotateMetadata[1]
|
|
|
325 dWeight = as.numeric(vsRotateMetadata[2])
|
|
|
326 sOrdinationMetadata = dfDummyMetadata[sMetadata]*dWeight
|
|
|
327 dfData[sMetadata] = sOrdinationMetadata
|
|
|
328 }
|
|
|
329
|
|
|
330 # Run NMDS on bug data (Default B-C)
|
|
|
331 # Will have species and points because working off of raw data
|
|
|
332 mNMDSData = metaMDS(dfData,k=2)
|
|
|
333
|
|
|
334 ## Make shapes
|
|
|
335 # Defines thes shapes and the metadata they are based on
|
|
|
336 # Metadata to use as shapes
|
|
|
337 lShapeInfo = funcMakeShapes(dfInput=dfInput, sShapeBy=sShapeBy, sShapes=sShapes, cDefaultShape=sDefaultMarker)
|
|
|
338
|
|
|
339 sMetadataShape = lShapeInfo[["ID"]]
|
|
|
340 vsShapeValues = lShapeInfo[["Values"]]
|
|
|
341 vsShapeShapes = lShapeInfo[["Shapes"]]
|
|
|
342 vsShapes = lShapeInfo[["PlotShapes"]]
|
|
|
343 cDefaultShape = lShapeInfo[["DefaultShape"]]
|
|
|
344
|
|
|
345 # Colors
|
|
|
346 vsColors = rep(cDefaultColor,nrow(dfInput))
|
|
|
347 vsColorValues = c()
|
|
|
348 vsColorRBG = c()
|
|
|
349 if(!is.null(sColorBy))
|
|
|
350 {
|
|
|
351 vsColorValues = paste(sort(unique(unlist(dfInput[[sColorBy]])),na.last=TRUE))
|
|
|
352 iLengthColorValues = length(vsColorValues)
|
|
|
353
|
|
|
354 vsColorRBG = lapply(1:iLengthColorValues/iLengthColorValues,colorRamp(vsColorRange))
|
|
|
355 vsColorRBG = unlist(lapply(vsColorRBG, function(x) rgb(x[1]/255,x[2]/255,x[3]/255)))
|
|
|
356
|
|
|
357 for(iColor in 1:length(vsColorRBG))
|
|
|
358 {
|
|
|
359 vsColors[which(paste(dfInput[[sColorBy]])==vsColorValues[iColor])]=vsColorRBG[iColor]
|
|
|
360 }
|
|
|
361
|
|
|
362 #If NAs are seperately given color, then color here
|
|
|
363 if(!is.null(sPlotNAColor))
|
|
|
364 {
|
|
|
365 vsColors[which(is.na(dfInput[[sColorBy]]))] = sPlotNAColor
|
|
|
366 vsColorRBG[which(vsColorValues=="NA")] = sPlotNAColor
|
|
|
367 }
|
|
|
368 }
|
|
|
369
|
|
|
370 # Reduce the bugs down to the ones in the list to be plotted
|
|
|
371 viBugsToPlot = which(row.names(mNMDSData$species) %in% vsBugsToPlot)
|
|
|
372 viMetadataDummy = which(names(dfDummyMetadata) %in% vsMetadata)
|
|
|
373
|
|
|
374 # Build the matrix of metadata coordinates
|
|
|
375 mMetadataCoordinates = matrix(rep(NA, iNumberMetadata*2),nrow=iNumberMetadata)
|
|
|
376
|
|
|
377 for( i in 1:iNumberMetadata )
|
|
|
378 {
|
|
|
379 lxReturn = NA
|
|
|
380 if( i %in% viContinuousMetadata )
|
|
|
381 {
|
|
|
382 lxReturn = funcGetMaximumForMetadatum(dfDummyMetadata[[i]],mNMDSData$points)
|
|
|
383 } else {
|
|
|
384 lxReturn = funcGetCentroidForMetadatum(dfDummyMetadata[[i]],mNMDSData$points)
|
|
|
385 }
|
|
|
386 mMetadataCoordinates[i,] = c(lxReturn$x,lxReturn$y)
|
|
|
387 }
|
|
|
388 row.names(mMetadataCoordinates) = vsRowNames
|
|
|
389
|
|
|
390 #!# metadata default dependent
|
|
|
391 # Plot the biplot with the centroid constructed metadata coordinates
|
|
|
392 if( ( length( viMetadataDummy ) == 0 ) && ( is.null( sMetadata ) ) )
|
|
|
393 {
|
|
|
394 viMetadataDummy = 1:nrow(mMetadataCoordinates)
|
|
|
395 }
|
|
|
396
|
|
|
397 # Plot samples
|
|
|
398 # Make output name
|
|
|
399 if(is.null(sOutputFileName))
|
|
|
400 {
|
|
|
401 viPeriods = which(sInputFileName==".")
|
|
|
402 if(length(viPeriods)>0)
|
|
|
403 {
|
|
|
404 sOutputFileName = paste(OutputFileName[1:viPeriods[length(viPeriods)]],"pdf",sep=".")
|
|
|
405 } else {
|
|
|
406 sOutputFileName = paste(sInputFileName,"pdf",sep=".")
|
|
|
407 }
|
|
|
408 }
|
|
|
409
|
|
|
410 pdf(sOutputFileName,useDingbats=FALSE)
|
|
|
411 plot(mNMDSData$points, xlab=paste("NMDS1","Stress=",mNMDSData$stress), ylab="NMDS2", pch=vsShapes, col=vsColors)
|
|
|
412 title(sTitle,line=3)
|
|
|
413 # Plot Bugs
|
|
|
414 mPlotBugs = mNMDSData$species[viBugsToPlot,]
|
|
|
415 if(length(viBugsToPlot)==1)
|
|
|
416 {
|
|
|
417 text(x=mPlotBugs[1],y=mPlotBugs[2],labels=row.names(mNMDSData$species)[viBugsToPlot],col=sTextColor)
|
|
|
418 } else if(length(viBugsToPlot)>1){
|
|
|
419 text(x=mPlotBugs[,1],y=mPlotBugs[,2],labels=row.names(mNMDSData$species)[viBugsToPlot],col=sTextColor)
|
|
|
420 }
|
|
|
421
|
|
|
422 # Add alternative axes
|
|
|
423 axis(3, col=sArrowColor)
|
|
|
424 axis(4, col=sArrowColor)
|
|
|
425 box(col = "black")
|
|
|
426
|
|
|
427 # Plot Metadata
|
|
|
428 if(length(viMetadataDummy)>0)
|
|
|
429 {
|
|
|
430 if(fPlotArrow)
|
|
|
431 {
|
|
|
432 # Plot arrows
|
|
|
433 for(i in viMetadataDummy)
|
|
|
434 {
|
|
|
435 curCoordinates = mMetadataCoordinates[i,]
|
|
|
436 curCoordinates = curCoordinates * dResizeArrow
|
|
|
437 # Plot Arrow
|
|
|
438 arrows(0,0, curCoordinates[1] * 0.8, curCoordinates[2] * 0.8, col=sArrowColor, length=0.1 )
|
|
|
439 }
|
|
|
440 }
|
|
|
441 # Plot text
|
|
|
442 if(length(viMetadataDummy)==1)
|
|
|
443 {
|
|
|
444 text(x=mMetadataCoordinates[viMetadataDummy,][1]*dResizeArrow*0.8, y=mMetadataCoordinates[viMetadataDummy,][2]*dResizeArrow*0.8, labels=row.names(mMetadataCoordinates)[viMetadataDummy],col=sArrowTextColor)
|
|
|
445 } else {
|
|
|
446 text(x=mMetadataCoordinates[viMetadataDummy,1]*dResizeArrow*0.8, y=mMetadataCoordinates[viMetadataDummy,2]*dResizeArrow*0.8, labels=row.names(mMetadataCoordinates)[viMetadataDummy],col=sArrowTextColor)
|
|
|
447 }
|
|
|
448 }
|
|
|
449
|
|
|
450 sLegendText = c(paste(vsColorValues,sColorBy,sep="_"),paste(vsShapeValues,sMetadataShape,sep="_"))
|
|
|
451 sLegendShapes = c(rep(cDefaultShape,length(vsColorValues)),vsShapeShapes)
|
|
|
452 sLegendColors = c(vsColorRBG,rep(cDefaultColor,length(vsShapeValues)))
|
|
|
453 if(length(sLegendText)>0)
|
|
|
454 {
|
|
|
455 legend("topright",legend=sLegendText,pch=sLegendShapes,col=sLegendColors)
|
|
|
456 }
|
|
|
457
|
|
|
458 # Original biplot call if you want to check the custom ploting of the script
|
|
|
459 # 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.
|
|
|
460 # Axes to the top and right are for the arrow, otherse are for markers and bug names.
|
|
|
461 #biplot(mNMDSData$points,mMetadataCoordinates[viMetadataDummy,],xlabs=vsShapes,xlab=paste("MDS1","Stress=",mNMDSData$stress),main="Biplot function Bugs and Sampes - Metadata Plotted with Centroids")
|
|
|
462 dev.off()
|
|
|
463 }
|
|
|
464
|
|
|
465 # This is the equivalent of __name__ == "__main__" in Python.
|
|
|
466 # That is, if it's true we're being called as a command line script;
|
|
|
467 # if it's false, we're being sourced or otherwise included, such as for
|
|
|
468 # library or inlinedocs.
|
|
|
469 if( identical( environment( ), globalenv( ) ) &&
|
|
|
470 !length( grep( "^source\\(", sys.calls( ) ) ) )
|
|
|
471 {
|
|
|
472 lsArgs <- parse_args( pArgs, positional_arguments=TRUE )
|
|
|
473
|
|
|
474 print("lsArgs")
|
|
|
475 print(lsArgs)
|
|
|
476
|
|
|
477 funcDoBiplot(
|
|
|
478 sBugs = lsArgs$options$sBugs,
|
|
|
479 sMetadata = lsArgs$options$sMetadata,
|
|
|
480 sColorBy = lsArgs$options$sColorBy,
|
|
|
481 sColorRange = lsArgs$options$sColorRange,
|
|
|
482 sTextColor = lsArgs$options$sTextColor,
|
|
|
483 sArrowColor = lsArgs$options$sArrowColor,
|
|
|
484 sArrowTextColor = lsArgs$options$sArrowTextColor,
|
|
|
485 sPlotNAColor = lsArgs$options$sPlotNAColor,
|
|
|
486 sShapeBy = lsArgs$options$sShapeBy,
|
|
|
487 sShapes = lsArgs$options$sShapes,
|
|
|
488 sDefaultMarker = lsArgs$options$sDefaultMarker,
|
|
|
489 sRotateByMetadata = lsArgs$options$sRotateByMetadata,
|
|
|
490 dResizeArrow = lsArgs$options$dResizeArrow,
|
|
|
491 fPlotArrow = !lsArgs$options$fNoPlotMetadataArrows,
|
|
|
492 sTitle = lsArgs$options$sTitle,
|
|
|
493 sInputFileName = lsArgs$args[2],
|
|
|
494 sLastMetadata = lsArgs$args[1],
|
|
|
495 sOutputFileName = lsArgs$options$sOutputFileName)
|
|
|
496 }
|
