|
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 }
|
