-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathDataset 10. PCA script.R
339 lines (279 loc) · 19.3 KB
/
Dataset 10. PCA script.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
### Load packages
require( geomorph )
require( ape )
#Clear workspace
rm( list = ls() )
#Set working directory to where turtle landmark data (from supplemental folder: Dataset 3. Turtle labyrinth data) is stored
setwd( "INSERT DIRECTORY PATH" )
#Load landmark data and transform to correct format for GPA commands
temp.file <- list.files(pattern = ".csv")
landmark.data <- lapply (temp.file, read.csv, row.names=1)
names(landmark.data) <- gsub(".csv","",temp.file)
landmark.data.temp <- array(as.numeric(unlist(landmark.data)), dim = c(123, 3, 184)) #number of landmarks, number of dimensions (3D coordinates), number of specimens
dimnames(landmark.data.temp)[[3]] <- gsub(".csv","",temp.file)
dimnames(landmark.data.temp)[[1]] <- rownames(landmark.data[[1]])
dimnames(landmark.data.temp)[[2]] <- c("x","y","z")
#Load slider information and colour information for deformation plots
setwd( "INSERT DIRECTORY PATH" )
sliders <- read.csv("Dataset 4. sliders.turtles.csv", row.names=1)
colours <- as.character( read.csv("Dataset 5. landmark_colours.csv", row.names=1)[,1] )
#Load specimen information
setwd( "INSERT DIRECTORY PATH" )
specimen.info <- read.csv( "Dataset 2. Specimen info.csv", header = TRUE )
rownames( specimen.info ) <- specimen.info[ , "Specimen_name" ]
#Load tree
setwd( "INSERT DIRECTORY PATH" )
tree <- read.nexus( "Dataset 7. cal3tree.calibrated.txt" )
alternative.tree <- read.nexus("Dataset 8. mbltree.calibrated.txt")
##GPA analysis and PCA
#Do GPA of labyrinth shape for all taxa available
GPA.data <- landmark.data.temp
#exclude all membranous specimens from main dataset
membranous.specimens <- c(which(grepl("membranous", dimnames( GPA.data )[[ 3 ]] ) == TRUE))
membranous.specimens.check <- dimnames( GPA.data )[[ 3 ]][membranous.specimens]
GPA.data <- GPA.data[,, - membranous.specimens]
tree.names <- as.character( specimen.info[ dimnames( GPA.data )[[ 3 ]] , "Tree_names" ] )
dimnames( GPA.data )[[ 3 ]][ !is.na( tree.names ) ] <- tree.names[ !is.na( tree.names ) ]
GPA.labyrinth.all <- gpagen( GPA.data , curves = sliders , ProcD = F )
#Extract centroid sizes for later scaling
labyrinth.Csize <- GPA.labyrinth.all$Csize
labyrinth.Csize[ labyrinth.Csize > 5000 ] <- labyrinth.Csize[ labyrinth.Csize > 5000 ] / 1000
#DO PCA
PCA.labyrinth <- plotTangentSpace( GPA.labyrinth.all$coords , warpgrids = F )
superposed.lms <- GPA.labyrinth.all
#Extract order of species in PCA data
datarownames <- rownames(PCA.labyrinth$pc.scores)
get.two <- function( X ) { X[1:2] }
speciesPCAorder <- unlist( lapply( lapply( strsplit( datarownames , "_" ) , get.two ) , paste , collapse = "_" ) )
#Make reduced info file with columns relevant for plotting
reduced.info <- specimen.info[, c("Species_ID", "Plotting_group", "Plotting_habitat", "Replicates_endosseous", "Tree_node_number")]
#Delete all entries of specimens based on membranous labyrinths
reduced.info <- reduced.info[ - c(which(grepl("membranous", rownames( reduced.info )) == TRUE)) , ]
#Check if these are in the same order/sequence
speciesPCAorder == reduced.info$Species_ID
#These should now contain the same number of species
length(reduced.info[, "Species_ID"]) == length(speciesPCAorder)
#Extract families from reduced info sheet, which are now ordered according the PCA data
plot.group <- as.character( reduced.info[, "Plotting_group" ] )
names(plot.group) <- rownames(reduced.info)
unique(plot.group)
#Colour sheme from Roger
good.colours <- function( n ) {
c( "firebrick2" , "seagreen3" , "aquamarine2" , "slategray3" , "darkorchid2" , "cadetblue3" , "darkred" , "coral2" , "chartreuse3" , "grey48" , "orange2" , "springgreen2" , "salmon2" , "dodgerblue2" )[ 1:n ]
}
#For colour assignment
specimen.order <- c( "early-diverging" , "Perichelydia" , "Xinjiangchelyidae" , "Macrobaenidae/Sinemydidae" , "Angolachelonia" , "Chelidae" , "Pelomedusoides" , "Trionychia" , "Testudinidae" , "Geoemydidae", "Emysternia", "Chelydridae" , "Kinosternoidea" , "Chelonioidea" )
#Assign colours according to families
group.bg <- good.colours( length( unique( plot.group ) ) )[ match( plot.group , specimen.order ) ]
names( group.bg ) <- plot.group
#Tests if colours indeed match families (works)
#reduced.info[,"Colours"]<-group.bg
#reduced.info[,-c(3:4)]
#Colour sheme for ecological categories
ecology.bg <- c( "dodgerblue3" , "lightblue2" , "tan" )
names( ecology.bg ) <- c( "marine" , "freshwater" , "terrestrial" )
ecology.bg <- ecology.bg[ as.character( reduced.info[ , "Plotting_habitat" ] ) ]
#extract numbers for point labels on plot
numbers.for.plot <- reduced.info[, "Tree_node_number"]
#Delete all entries of specimens based on membranous labyrinths
reduced.info.two <- specimen.info[ - c(which(grepl("membranous", rownames( specimen.info )) == TRUE)) , ]
#Check if these are in the same order/sequence
speciesPCAorder == reduced.info.two$Species_ID
#These should now contain the same number of species
length(reduced.info.two[, "Species_ID"]) == length(speciesPCAorder)
#Compute skull heigth/width ratio
skull_geometry <- reduced.info.two$Skull_height_mm / reduced.info.two$Skull_width_mm
names(skull_geometry) <- names(labyrinth.Csize)
#Object that contains species with incomplete data
no.skull.geometry.data <- names(which( is.na(skull_geometry) ) )
#Computes average values for skull geometry from species with data
average.numeric.values <- sum( skull_geometry[- c( which( is.na(skull_geometry) ) )] ) / length(skull_geometry[- c( which( is.na(skull_geometry) ) )])
#Assigns average values to species without data
skull_geometry.temp <- skull_geometry
skull_geometry.temp[which( is.na(skull_geometry.temp) )] <- average.numeric.values
#Scales all values between 0 and 1
skull.geometry.relative <- (as.numeric( skull_geometry.temp ) - min( as.numeric( skull_geometry.temp ))) / max(as.numeric( skull_geometry.temp ) - min( as.numeric( skull_geometry.temp )))
names(skull.geometry.relative) <- names(labyrinth.Csize)
#Re-assigns NAs for speciess without data
skull.geometry.relative[no.skull.geometry.data] <- NA
##Plot PC space
PCA.results <- PCA.labyrinth
labelled = T #Select if you want taxon names to appear on plot
labelled = F #Select if you do not want taxon names to appear on plot
#Choose which PC axes to show and create axis labels that state the proportion of variance explained
X <- "PC1"
XLAB <- paste( X , " (" , round( PCA.results$pc.summary$importance[ 2 , X ] , 3 ) * 100 , "%" , ")" , sep = "")
Y <- "PC2"
YLAB <- paste( Y , " (" , round( PCA.results$pc.summary$importance[ 2 , Y ] , 3 ) * 100 , "%" , ")" , sep = "")
Z <- "PC3"
ZLAB <- paste( Z , " (" , round( PCA.results$pc.summary$importance[ 2 , Z ] , 3 ) * 100 , "%" , ")" , sep = "")
dev.new( width = 8 , height = 8 )
close.screen( all.screens = T )
split.screen( c( 2 , 2 ) )
screen( 1 ) ; par( mar = c( 4, 4, 1, 1 ) )
#I could use plot( PCA.results$pc.scores[ , c(X,Y) ] , parameters), but wanted to flip the PC1 axis so that my previous descriptive text is still true
plot( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , bty = "n" , pch = 21 , bg = group.bg, lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = YLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Y] ) + c( 0 , 0.07 ) , col = "darkgrey" )
if(labelled) {text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , rownames( PCA.results$pc.scores ) , cex = 0.2 )}
text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , labels = numbers.for.plot , cex = 0.3 )
legend( "topleft" , legend = specimen.order , pch = 21 , pt.bg = good.colours( length( unique( plot.group ) ) ) , cex = 0.5 , pt.cex = 1.5 , bty = "n" )
screen( 2 ) ; par( mar = c( 4, 4, 1, 1 ) )
plot( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , bty = "n" , pch = 21 , bg = group.bg, lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = ZLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Z] ) + c( 0 , 0.07 ) , col = "darkgrey" )
if(labelled) {text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , rownames( PCA.results$pc.scores ) , cex = 0.2 )}
text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , labels = numbers.for.plot , cex = 0.3 )
screen( 3 ) ; par( mar = c( 4, 4, 1, 1 ) )
plot( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , bty = "n" , pch = 21 , bg = ecology.bg, lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = YLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Y] ) + c( 0 , 0.07 ) , col = "darkgrey" )
if(labelled) {text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , rownames( PCA.results$pc.scores ) , cex = 0.2 )}
text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , labels = numbers.for.plot , cex = 0.3 )
legend( "topleft" , legend = c( "marine" , "freshwater" , "terrestrial" ) , pch = 21 , pt.bg = c( "dodgerblue3" , "lightblue2" , "tan" ) , cex = 0.5 , pt.cex = 1.5 , bty = "n" )
screen( 4 ) ; par( mar = c( 4, 4, 1, 1 ) )
plot( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , bty = "n" , pch = 21 , bg = ecology.bg, lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = ZLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Z] ) + c( 0 , 0.07 ) , col = "darkgrey" )
if(labelled) {text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , rownames( PCA.results$pc.scores ) , cex = 0.2 )}
text( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , labels = numbers.for.plot , cex = 0.3 )
dev.new()
plot( PCA.results$pc.scores[ , X ] , skull_geometry , , bty = "n" , pch = 21 , bg = ecology.bg, lwd = 1.2 , cex = 1.9 , , xlab = XLAB , ylab = "skull height-width ratio" , cex.lab = 0.9 , cex.axis = 0.8)
### Write source data file
source.data <- reduced.info.two[, c("Specimen_name", "Species_ID", "Plotting_habitat", "Tree_node_number")]
PC1.values <- -PCA.results$pc.scores[ , X ]
PC2.values <- PCA.results$pc.scores[ , Y ]
PC3.values <- PCA.results$pc.scores[ , Z ]
source.data[,"PC1.values"] <- PC1.values
source.data[,"PC2.values"] <- PC2.values
source.data[,"PC3.values"] <- PC3.values
write.table(source.data, file = "Source_Data_file_1.txt")
###
# Plot individual families against all other points
#Define colours for single families and the legend
i <- specimen.order[1]
single.group.bg <- group.bg
single.group.bg[ which(names(single.group.bg[]) != i ) ] <- "grey"
for.labels <- c(which(single.group.bg[] != "grey"))
legend.single.group.text <- c(i, "other")
legend.single.group.bg <- c(single.group.bg[[i]], "grey")
for.points.x <- -PCA.results$pc.scores[ , X ][which(plot.group == i) ]
for.points.y <- PCA.results$pc.scores[ , Y ][which(plot.group == i) ]
for.points.z <- PCA.results$pc.scores[ , Z ][which(plot.group == i) ]
#Plot
dev.new( width = 8 , height = 4 )
close.screen( all.screens = T )
split.screen( c( 1 , 2 ) )
screen( 1 ) ; par( mar = c( 4, 4, 1, 1 ) )
plot( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Y ] , bty = "n" , pch = 21 , bg = single.group.bg, lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = YLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Y] ) + c( 0 , 0.07 ) , col = "darkgrey" )
points( for.points.x , for.points.y , bty = "n" , pch = 21 , bg = group.bg[which(names(group.bg[]) == i )], lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = YLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Y] ) + c( 0 , 0.07 ) , col = "darkgrey" )
if(labelled) {text( for.points.x , for.points.y , rownames( PCA.results$pc.scores )[for.labels] , cex = 0.2 )}
text( for.points.x , for.points.y , labels = numbers.for.plot[for.labels] , cex = 0.3 )
legend( "topleft" , legend = legend.single.group.text , pch = 21 , pt.bg = legend.single.group.bg , cex = 0.5 , pt.cex = 1.5 , bty = "n" )
screen( 2 ) ; par( mar = c( 4, 4, 1, 1 ) )
plot( -PCA.results$pc.scores[ , X ], PCA.results$pc.scores[ , Z ] , bty = "n" , pch = 21 , bg = single.group.bg, lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = ZLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Z] ) + c( 0 , 0.07 ) , col = "darkgrey" )
points( for.points.x , for.points.z , bty = "n" , pch = 21 , bg = group.bg[which(names(group.bg[]) == i )], lwd = 1.2 , cex = 1.9 , xlab = XLAB , ylab = ZLAB , cex.lab = 0.9 , cex.axis = 0.8 , ylim = range( PCA.results$pc.scores[,Z] ) + c( 0 , 0.07 ) , col = "darkgrey" )
if(labelled) {text( for.points.x , for.points.z , rownames( PCA.results$pc.scores )[for.labels] , cex = 0.2 )}
text( for.points.x , for.points.z , labels = numbers.for.plot[for.labels] , cex = 0.3 )
####Plot PC shapes
##It all runs much faster if you omit the lines that plot lines connecting landmarks together
labelled = FALSE
lines.plot <- rbind( as.matrix(sliders[,c(1,2)]) , as.matrix(sliders[,c(2,3)]) )
open3d()
par3d(windowRect = c(0,0,1400,900))
Sys.sleep(1)
mfrow3d(nr = 3, nc = 4, byrow = TRUE, sharedMouse = TRUE)
#to match reversed PC1 axis I plit the max shape first and label it as min
choose <- "PC1max" #choose <- "PC1max"; reversed
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main="PC1min")} )
aspect3d( "iso" )
next3d()
#to match reversed PC1 axis I plit the min shape second and label it as max
choose <- "PC1min" #choose <- "PC1min"; reversed so this fits with my reversed axes in the plot
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main="PC1max")} )
aspect3d( "iso" )
next3d()
choose <- "PC4min"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC4max"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC2min"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC2max"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
choose <- "PC5min"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC5max"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC3min"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC3max"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC6min"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
next3d()
choose <- "PC6max"
plot3d( PCA.labyrinth$pc.shapes[[choose]] , col = colours , size = 10 , box = "n")
for( i in 1:nrow( lines.plot ) ) { lines3d( PCA.labyrinth$pc.shapes[[choose]][ lines.plot[i,] , ] , lwd = 2 , col = "darkgrey" ) }
points3d( superposed.lms$consensus , col = "grey" , size = 5 , box = "n")
if( labelled) { title3d( main = choose, line=1 , cex = 1) }
bgplot3d( { plot.new() ;title(main=choose)} )
aspect3d( "iso" )
##Screenshots for PC shapes.
rgl.snapshot( "PC shapes 1-6 lateral view turtles.png" )
rgl.snapshot( "PC shapes 1-6 dorsal view turtles.png" )
rgl.snapshot( "PC shapes 1-6 anterior canal view turtles.png" )
rgl.snapshot( "PC shapes 1-6 posterior canal view turtles.png" )
#