# SCRIPT FOR PERFORMING DATA ANALYSIS AND CREATING FIGURES # # Author: Matthew T Wayland # # # CONTENTS # # 1. SETUP # 2. ACCURACY OF HOOK AREA ESTIMATION # 3. GENERATE AND PLOT MERISTOGRAMS # 4. INTERSPECIFIC COMPARISON OF MERISTOGRAMS USING PCA AND CLUSTERING # 5. CREATE FIGURES FOR SUPPLEMENTARY MATERIALS # 6. GET MORPHOMETRIC DATA TO EXPAND DESCRIPTION OF HOOK PATTERNS # 1. SETUP # close all graphics devices graphics.off() # turn off stringsAsFactors options(stringsAsFactors = FALSE) # Load required packages (install if not present) if (!require("RCurl")) { install.packages("RCurl", dependencies = TRUE) library(RCurl) } if (!require("ape")) { install.packages("ape", dependencies = TRUE) library(ape) } if (!require("dendextend")) { install.packages("dendextend", dependencies = TRUE) library(dendextend) } # set working directory - location where figures will be saved if (.Platform$OS.type == "windows") { setwd("C:/Users/mw283/Dropbox/Acanthocephala/meristogram/figures") }else{ setwd("~/Dropbox/Acanthocephala/meristogram/figures") } # source meristogram code eval(expr = parse(text = getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/CLI/meristogram.R", ssl.verifypeer=FALSE) )) #-------------------------------------------------------------------------------------------------------------------------------- # 2. ACCURACY OF HOOK AREA ESTIMATION # comparison of measured and estimated hook area # read in hook area data (N.B. a complete row of hook data is not available for all specimens) hks <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_spp_estimated_and_measured_hk_area.csv")) # split data.frame of hook area data by taxon hksA <- hks[hks$species=="A",] hksB <- hks[hks$species=="B",] hksI <- hks[hks$species=="I",] # report some summary statistics cat(paste("Number of specimens of E. gadi sp. A: ", length(unique(hksA$specimen)), "\n", sep="")) cat(paste("Number of female E. gadi sp. A: ", length(unique(hksA$specimen[hksA$sex=="F"])), "\n", sep="")) cat(paste("Number of male E. gadi sp. A: ", length(unique(hksA$specimen[hksA$sex=="M"])), "\n", sep="")) cat(paste("Total number of hooks for E. gadi sp. A: ", length(hksA[,1]), "\n", sep="")) cat(paste("Number of specimens of E. gadi sp. B: ", length(unique(hksB$specimen)), "\n", sep="")) cat(paste("Number of female E. gadi sp. B: ", length(unique(hksB$specimen[hksB$sex=="F"])), "\n", sep="")) cat(paste("Number of male E. gadi sp. B: ", length(unique(hksB$specimen[hksB$sex=="M"])), "\n", sep="")) cat(paste("Total number of hooks for E. gadi sp. B: ", length(hksB[,1]), "\n", sep="")) cat(paste("Number of specimens of E. gadi sp. I: ", length(unique(hksI$specimen)), "\n", sep="")) cat(paste("Number of female E. gadi sp. I: ", length(unique(hksI$specimen[hksI$sex=="F"])), "\n", sep="")) cat(paste("Number of male E. gadi sp. I: ", length(unique(hksI$specimen[hksI$sex=="M"])), "\n", sep="")) cat(paste("Total number of hooks for E. gadi sp. I: ", length(hksI[,1]), "\n", sep="")) # analyse and plot hook area data for E. gadi sp. A cor.test(hksA$measured_area, hksA$estimated_area, method='pearson') png("E_gadi_sp_A_hk_area_estimated_measured.png", res=1000, units="mm", width=120, height=120) par(mar=c(4, 4.5, 1, 1) + 0.1) plot(hksA$measured_area, hksA$estimated_area, cex=0.7, pch=1, xlab = expression(paste('Measured Area', ' ', mu,m^2)), ylab = expression(paste('Estimated Area', ' ', mu,m^2))) xy=0:1000 lines(xy,xy, col="blue", lty=2) hksAlm <- lm(hksA$estimated_area ~ hksA$measured_area) summary(hksAlm) abline(hksAlm, col="red", lty=1) dev.off() # analyse and plot hook area data for E. gadi sp. B cor.test(hksB$measured_area, hksB$estimated_area, method='pearson') png("E_gadi_sp_B_hk_area_estimated_measured.png", res=1000, units="mm", width=120, height=120) par(mar=c(4, 4.5, 1, 1) + 0.1) plot(hksB$measured_area, hksB$estimated_area, cex=0.7, pch=1, xlab = expression(paste('Measured Area', ' ', mu,m^2)), ylab = expression(paste('Estimated Area', ' ', mu,m^2))) xy=0:1000 lines(xy,xy, col="blue", lty=2) hksBlm <- lm(hksB$estimated_area ~ hksB$measured_area) summary(hksBlm) abline(hksBlm, col="red", lty=1) dev.off() # anlayse and plot hook area data for E. gadi sp. I cor.test(hksI$measured_area, hksI$estimated_area, method='pearson') png("E_gadi_sp_I_hk_area_estimated_measured.png", res=1000, units="mm", width=120, height=120) par(mar=c(4, 4.5, 1, 1) + 0.1) plot(hksI$measured_area, hksI$estimated_area, cex=0.7, pch=1, xlab = expression(paste('Measured Area', ' ', mu,m^2)), ylab = expression(paste('Estimated Area', ' ', mu,m^2))) xy=0:1000 lines(xy,xy, col="blue", lty=2) hksIlm <- lm(hksI$estimated_area ~ hksI$measured_area) summary(hksIlm) abline(hksIlm, col="red", lty=1) dev.off() #-------------------------------------------------------------------------------------------------------------------------------- # 3. GENERATE AND PLOT MERISTOGRAMS # read in hook data gadi_A_female_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_A_female.csv")) gadi_A_male_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_A_male.csv")) gadi_B_female_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_B_female.csv")) gadi_B_male_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_B_male.csv")) gadi_I_female_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_I_female.csv")) gadi_I_male_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/E_gadi_I_male.csv")) bothniensis_1_female_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_bothniensis_1_female.csv")) bothniensis_1_male_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_bothniensis_1_male.csv")) bothniensis_2_female_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_bothniensis_2_female.csv")) brayi_female_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_brayi_female.csv")) brayi_male_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_brayi_male.csv")) salmonis_female_dorsal_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_salmonis_female_dorsal.csv")) salmonis_female_ventral_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_salmonis_female_ventral.csv")) salmonis_male_dorsal_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_salmonis_male_dorsal.csv")) salmonis_male_ventral_raw <- read.csv(text=getURL("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_salmonis_male_ventral.csv")) truttae_female_Drummore_raw <- read.csv("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_truttae_Drummore_female.csv") truttae_male_Drummore_raw <- read.csv("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_truttae_Drummore_male.csv") truttae_female_Carron_raw <- read.csv("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_truttae_Carron_female.csv") truttae_male_Carron_raw <- read.csv("https://raw.githubusercontent.com/WaylandM/meristogram/master/data/Echinorhynchus_truttae_Carron_male.csv") # concatenate female and male datasets for each taxon gadi_A_raw <- rbind(gadi_A_female_raw, gadi_A_male_raw) gadi_B_raw <- rbind(gadi_B_female_raw, gadi_B_male_raw) gadi_I_raw <- rbind(gadi_I_female_raw, gadi_I_male_raw) bothniensis_1_raw <- rbind(bothniensis_1_female_raw, bothniensis_1_male_raw) bothniensis_2_raw<- bothniensis_2_female_raw brayi_raw <- rbind(brayi_female_raw, brayi_male_raw) salmonis_dorsal_raw <- rbind(salmonis_female_dorsal_raw, salmonis_male_dorsal_raw) salmonis_ventral_raw <- rbind(salmonis_female_ventral_raw, salmonis_male_ventral_raw) salmonis_raw <- rbind(salmonis_dorsal_raw, salmonis_ventral_raw) truttae_Drummore_raw <- rbind(truttae_female_Drummore_raw, truttae_male_Drummore_raw) truttae_Carron_raw <- rbind(truttae_female_Carron_raw, truttae_male_Carron_raw) # plot meristograms for each taxon using moving average interval of 17% (minimum moving average interval applicable to all taxa) # brayi (generate meristogram with and without linear interpolation for demonstration purposes) png("brayi_no_lerp.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(brayi_raw, 17, lerp = F), cex=0.8) dev.off() png("brayi_lerp.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(brayi_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # salmonis png("salmonis_dorsal.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(salmonis_dorsal_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() png("salmonis_ventral.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(salmonis_ventral_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # gadi A png("gadi_A.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(gadi_A_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # gadi B png("gadi_B.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(gadi_B_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # gadi I png("gadi_I.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(gadi_I_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # bothniensis png("bothniensis_1.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(bothniensis_1_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # 'bothniensis' png("bothniensis_2.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(bothniensis_2_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # truttae Drummore png("truttae_Drummore.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(truttae_Drummore_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() # truttae Carron png("truttae_Carron.png", res=1000, units="mm", width=180, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plotMeristogram(createMeristogram(truttae_Carron_raw, 17, lerp = T, rm.lerp.na = T), cex=0.8) dev.off() #-------------------------------------------------------------------------------------------------------------------------------- # 4. INTERSPECIFIC COMPARISON OF MERISTOGRAMS USING PCA AND CLUSTERING # E. salmonis has data for dorsal and ventral hook rows pooled to yield a single data-set for this taxon. # First plot phylogenetic relationsips of taxa in study, as inferred by Wayland et al. (2015) phylogramLabels <- expression(italic('E.')*' '*italic("'bothniensis'")*'', italic('E.')*' '*italic('bothniensis')*'', italic('E.')*' '*italic('gadi')*' sp. I', italic('E.')*' '*italic('truttae')*'', italic('E.')*' '*italic('brayi')*'', italic('E.')*' '*italic('salmonis')*'') etree <- read.tree(text = "(((((bothniensis,'bothniensis'),gadi1),truttae),brayi),salmonis);") etree <- ladderize(etree, right=F) etree$tip.label = phylogramLabels png("molecular_phylogeny.png", res=1000, units="mm", width=150, height=120) par(mar=c(1, 1, 1, 1) + 0.1) plot(etree, edge.color=c("blue"), edge.width=3, cex=1.5, cex.main=1.7) dev.off() # Principal Component Analysis (PCA) groups <- c("gadi_A", "gadi_B", "gadi_I", "bothniensis_1", "bothniensis_2", "brayi", "salmonis", "truttae_drummore", "truttae_carron") # create meristogram for each taxon with option rm.lerp.na set to false, so that all meristogram data.frames # have same number of rows (helpful for concatentaing data-sets) gadi_A_meristo <- createMeristogram(gadi_A_raw, 17, lerp=T, rm.lerp.na = F) gadi_B_meristo <- createMeristogram(gadi_B_raw, 17, lerp=T, rm.lerp.na = F) gadi_I_meristo <- createMeristogram(gadi_I_raw, 17, lerp=T, rm.lerp.na = F) bothniensis_1_meristo <- createMeristogram(bothniensis_1_raw, 17, lerp=T, rm.lerp.na = F) bothniensis_2_meristo <- createMeristogram(bothniensis_2_raw, 17, lerp=T, rm.lerp.na = F) brayi_meristo <- createMeristogram(brayi_raw, 17, lerp=T, rm.lerp.na = F) salmonis_dorsal_meristo <- createMeristogram(salmonis_dorsal_raw, 17, lerp=T, rm.lerp.na = F) salmonis_ventral_meristo <- createMeristogram(salmonis_ventral_raw, 17, lerp=T, rm.lerp.na = F) truttae_Drummore_meristo <- createMeristogram(truttae_Drummore_raw, 17, lerp=T, rm.lerp.na=F) truttae_Carron_meristo <- createMeristogram(truttae_Carron_raw, 17, lerp=T, rm.lerp.na=F) salmonis_meristo <- createMeristogram(salmonis_raw, 17, lerp=T, rm.lerp.na=F) # concatenate meristograms in to a single data.frame meristoL <- as.data.frame(cbind(gadi_A_meristo$length, gadi_B_meristo$length, gadi_I_meristo$length, bothniensis_1_meristo$length, bothniensis_2_meristo$length, brayi_meristo$length, salmonis_meristo$length, truttae_Drummore_meristo$length, truttae_Carron_meristo$length)) meristoB <- as.data.frame(cbind(gadi_A_meristo$base, gadi_B_meristo$base, gadi_I_meristo$base, bothniensis_1_meristo$base, bothniensis_2_meristo$base, brayi_meristo$base, salmonis_meristo$base, truttae_Drummore_meristo$base, truttae_Carron_meristo$base)) meristoA <- as.data.frame(cbind(gadi_A_meristo$area, gadi_B_meristo$area, gadi_I_meristo$area, bothniensis_1_meristo$area, bothniensis_2_meristo$area, brayi_meristo$area, salmonis_meristo$area, truttae_Drummore_meristo$area, truttae_Carron_meristo$area)) meristoR <- as.data.frame(cbind(gadi_A_meristo$ratio, gadi_B_meristo$ratio, gadi_I_meristo$ratio, bothniensis_1_meristo$ratio, bothniensis_2_meristo$ratio, brayi_meristo$ratio, salmonis_meristo$ratio, truttae_Drummore_meristo$ratio, truttae_Carron_meristo$ratio)) names(meristoL) <- groups names(meristoB) <- groups names(meristoA) <- groups names(meristoR) <- groups tMeristoL <- t(meristoL) tMeristoB <- t(meristoB) tMeristoA <- t(meristoA) tMeristoR <- t(meristoR) tMeristo <- as.data.frame(cbind(tMeristoL, tMeristoB, tMeristoA, tMeristoR)) colHeaders <- c(unlist(lapply(0:100, function(x){paste("L", x, sep="")})), unlist(lapply(0:100, function(x){paste("B", x, sep="")})), unlist(lapply(0:100, function(x){paste("A", x, sep="")})), unlist(lapply(0:100, function(x){paste("R", x, sep="")}))) names(tMeristo) <- colHeaders tMeristo <- tMeristo[,!apply(tMeristo, 2, function(x){any(is.na(x))})] # PCA # taxon labels of PCA plot pcaTaxonLabels <- expression(italic('E.')*' '*italic('gadi')*' sp. A', italic('E.')*' '*italic('gadi')*' sp. B', italic('E.')*' '*italic('gadi')*' sp. I', italic('E.')*' '*italic('bothniensis')*'', italic('E.')*' '*italic("'bothniensis'")*'', italic('E.')*' '*italic('brayi')*'', italic('E.')*' '*italic('salmonis')*'', italic('E.')*' '*italic('truttae')*' (Drummore)', italic('E.')*' '*italic('truttae')*' (Carron)') # analysis meristoPCA <- prcomp(tMeristo, scale=F, na.action=na.omit) # show variance associate with each PC summary(meristoPCA)$importance[1:3,1:2] # PCA scores plot png("meristogram_pca_scores.png", res=1000, units="mm", width=120, height=120) par(mar=c(4, 4, 1, 1) + 0.1) plot(meristoPCA$x[,1], meristoPCA$x[,2], cex=1.5, col="red", pch=16, xlab="PC1", ylab="PC2") text(meristoPCA$x[,1], meristoPCA$x[,2], labels=pcaTaxonLabels, pos=c(2, 4, 4,4,2 ,4, 4,2,2), col="blue", cex=0.7) dev.off() # PCA loadings plot png("meristogram_pca_loadings.png", res=1000, units="mm", width=120, height=120) par(mar=c(4, 4, 1, 1) + 0.1) loadingsColours <- rep(c(rep("#1B9E77", 18), rep("#D95F02", 18), rep("#7570B3", 18), rep("#E7298A", 17)), 4) plot(meristoPCA$rotation[,1], meristoPCA$rotation[,2], col=loadingsColours, pch=c(rep("L", 71), rep("B", 71), rep("A", 71), rep("R", 71)), xlab="PC1", ylab="PC2", cex=0.8) legend("bottomright", c("15-32%", "33-50%", "51-68%", "69-85%"), title.col="black", col=c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"), pch=16, text.col=c("#1B9E77", "#D95F02", "#7570B3", "#E7298A"), title="Hook Position") dev.off() # create dendrogram labels dendrogramLabels <- expression(italic('E.')*' '*italic('salmonis')*'', italic('E.')*' '*italic('brayi')*'', italic('E.')*' '*italic('truttae')*' (Carron)', italic('E.')*' '*italic('truttae')*' (Drummore)', italic('E.')*' '*italic('gadi')*' sp. I', italic('E.')*' '*italic('bothniensis')*'', italic('E.')*' '*italic('gadi')*' sp. A', italic('E.')*' '*italic("'bothniensis'")*'', italic('E.')*' '*italic('gadi')*' sp. B') # hierarchical clustering of PC1 and PC2 using Euclidean distance measure dendrogram <- as.dendrogram(hclust(dist(meristoPCA$x[,1:2], method='euclidean'), method='average'), hang=0.02) dendrogram <- reorder(dendrogram, c(9:1)) labels(dendrogram) <- dendrogramLabels png("cluster_dendrogram_Euclidean.png", res=1000, units="mm", width=150, height=120) par(mar=c(4, 1, 1.5, 4) + 0.1) plot(dendrogram, ylab='', sub='', xlab='Distance', main="Euclidean", horiz=T) dev.off() # hierarchical clustering of PC1 and PC2 using Manhattan distance measure dendrogram <- as.dendrogram(hclust(dist(meristoPCA$x[,1:2], method='manhattan'), method='average'), hang=0.02) dendrogram <- reorder(dendrogram, c(9:1)) labels(dendrogram) <- dendrogramLabels png("cluster_dendrogram_Manhattan.png", res=1000, units="mm", width=150, height=120) par(mar=c(4, 1, 1.5, 4) + 0.1) plot(dendrogram, ylab='', sub='', xlab='Distance', main="Manhattan", horiz=T) dev.off() # hierarchical clustering of PC1 and PC2 using maximum distance measure dendrogram <- as.dendrogram(hclust(dist(meristoPCA$x[,1:2], method='maximum'), method='average'), hang=0.02) dendrogram <- reorder(dendrogram, c(9:1)) labels(dendrogram) <- dendrogramLabels png("cluster_dendrogram_Maximum.png", res=1000, units="mm", width=150, height=120) par(mar=c(4, 1, 1.5, 4) + 0.1) plot(dendrogram, ylab='', sub='', xlab='Distance', main="Maximum", horiz=T) dev.off() # hierarchical clustering of PC1 and PC2 using Minkowski distance measure dendrogram <- as.dendrogram(hclust(dist(meristoPCA$x[,1:2], method='minkowski'), method='average'), hang=0.02) dendrogram <- reorder(dendrogram, c(9:1)) labels(dendrogram) <- dendrogramLabels png("cluster_dendrogram_Minkowski.png", res=1000, units="mm", width=150, height=120) par(mar=c(4, 1, 1.5, 4) + 0.1) plot(dendrogram, ylab='', sub='', xlab='Distance', main="Minkowski", horiz=T) dev.off() #-------------------------------------------------------------------------------------------------------------------------------- # 5. CREATE FIGURES FOR SUPPLEMENTARY MATERIALS # set working directory for supplementary figures setwd("suppl_figs") # For most taxa meristograms for the supplementary material can be generated using this function supplMeristograms <- function(acanthocephalan) { # create pdf, calculate movAvg seg sizes, etc f <- get(paste(acanthocephalan, "female_raw", sep='_')) m <- get(paste(acanthocephalan, "male_raw", sep='_')) movAvgMin <- max(c(minMovAvgInterval(minHksPerRow(f)), minMovAvgInterval(minHksPerRow(m)))) pdf(paste(acanthocephalan, "suppl_meristograms.pdf", sep='_'), paper="a4", height=11, width=8) par(mar=c(3,4,4,2)) par(mgp=c(2,1,0)) par(mfrow=c(3,2)) plotMeristogram(createMeristogram(f, movAvgMin), main=paste("female ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(m, movAvgMin), main=paste("male ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(f, round(movAvgMin*1.5)), main=paste("female ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(m, round(movAvgMin*1.5)), main=paste("male ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(f, movAvgMin*2), main=paste("female ", movAvgMin*2, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(m, movAvgMin*2), main=paste("male ", movAvgMin*2, "%", sep=''), cex=0.8) dev.off() } # Use supplMeristograms function to generate supplementary meristograms for the following taxa: brayi, gadi A, gadi B, gadi I, bothniensis and truttae lapply(c("brayi", "gadi_A", "gadi_B", "gadi_I", "bothniensis_1"), supplMeristograms) # 'bothniensis' (cannot use supplMeristograms function, because we do not have data for males) movAvgMin <- minMovAvgInterval(minHksPerRow(bothniensis_2_raw)) pdf("bothniensis_2_suppl_meristograms.pdf", paper="a4", height=11, width=8) par(mar=c(3,4,4,2)) par(mgp=c(2,1,0)) par(mfrow=c(3,1)) plotMeristogram(createMeristogram(bothniensis_2_female_raw, movAvgMin), main=paste("female ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(bothniensis_2_female_raw, round(movAvgMin*1.5)), main=paste("female ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(bothniensis_2_female_raw, movAvgMin*2), main=paste("female ", movAvgMin*2, "%", sep=''), cex=0.8) dev.off() # salmonis (cannot use supplMeristograms function, because we have separate datasets for dorsal and ventral rows of hooks) movAvgMin <- max(c(minMovAvgInterval(minHksPerRow(salmonis_female_dorsal_raw)), minMovAvgInterval(minHksPerRow(salmonis_female_ventral_raw)), minMovAvgInterval(minHksPerRow(salmonis_male_dorsal_raw)), minMovAvgInterval(minHksPerRow(salmonis_male_ventral_raw)))) pdf("salmonis_suppl_meristograms.pdf", paper="a4", height=11, width=8, onefile=T) par(mar=c(3,4,4,2)) par(mgp=c(2,1,0)) par(mfrow=c(3,2)) #females plotMeristogram(createMeristogram(salmonis_female_dorsal_raw, movAvgMin), main=paste("female dorsal ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_female_ventral_raw, movAvgMin), main=paste("female ventral ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_female_dorsal_raw, round(movAvgMin*1.5)), main=paste("female dorsal ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_female_ventral_raw, round(movAvgMin*1.5)), main=paste("female ventral ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_female_dorsal_raw, movAvgMin*2), main=paste("female dorsal ", movAvgMin*2, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_female_ventral_raw, movAvgMin*2), main=paste("female ventral ", movAvgMin*2, "%", sep=''), cex=0.8) #males plotMeristogram(createMeristogram(salmonis_male_dorsal_raw, movAvgMin), main=paste("male dorsal ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_male_ventral_raw, movAvgMin), main=paste("male ventral ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_male_dorsal_raw, round(movAvgMin*1.5)), main=paste("male dorsal ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_male_ventral_raw, round(movAvgMin*1.5)), main=paste("male ventral ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_male_dorsal_raw, movAvgMin*2), main=paste("male dorsal ", movAvgMin*2, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(salmonis_male_ventral_raw, movAvgMin*2), main=paste("male ventral ", movAvgMin*2, "%", sep=''), cex=0.8) dev.off() # truttae (cannot use supplMeristograms function, because we have separate datasets for Drummore and River Carron catchment) movAvgMin <- max(c(minMovAvgInterval(minHksPerRow(truttae_female_Drummore_raw)), minMovAvgInterval(minHksPerRow(truttae_female_Carron_raw)), minMovAvgInterval(minHksPerRow(truttae_male_Drummore_raw)), minMovAvgInterval(minHksPerRow(truttae_male_Carron_raw)))) pdf("truttae_suppl_meristograms.pdf", paper="a4", height=11, width=8, onefile=T) par(mar=c(3,4,4,2)) par(mgp=c(2,1,0)) par(mfrow=c(3,2)) #females plotMeristogram(createMeristogram(truttae_female_Drummore_raw, movAvgMin), main=paste("Drummore female ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_female_Carron_raw, movAvgMin), main=paste("Carron female ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_female_Drummore_raw, round(movAvgMin*1.5)), main=paste("Drummore female ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_female_Carron_raw, round(movAvgMin*1.5)), main=paste("Carron female ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_female_Drummore_raw, movAvgMin*2), main=paste("Drummore female ", movAvgMin*2, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_female_Carron_raw, movAvgMin*2), main=paste("Carron female ", movAvgMin*2, "%", sep=''), cex=0.8) #males plotMeristogram(createMeristogram(truttae_male_Drummore_raw, movAvgMin), main=paste("Drummore male ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_male_Carron_raw, movAvgMin), main=paste("Carron male ", movAvgMin, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_male_Drummore_raw, round(movAvgMin*1.5)), main=paste("Drummore male ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_male_Carron_raw, round(movAvgMin*1.5)), main=paste("Carron male ", round(movAvgMin*1.5), "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_male_Drummore_raw, movAvgMin*2), main=paste("Drummore male ", movAvgMin*2, "%", sep=''), cex=0.8) plotMeristogram(createMeristogram(truttae_male_Carron_raw, movAvgMin*2), main=paste("Carron male ", movAvgMin*2, "%", sep=''), cex=0.8) dev.off() #-------------------------------------------------------------------------------------------------------------------------------- # 6. GET MORPHOMETRIC DATA TO EXPAND DESCRIPTION OF HOOK PATTERNS extractKeyHookMorphometrics <- function(hookMeasurements){ area <- hookMeasurements$length * hookMeasurements$base/2 ratio <- hookMeasurements$base * 100/hookMeasurements$length hookMeasurements <- cbind(hookMeasurements, area, ratio) names(hookMeasurements)[5:6] <- c("area", "ratio") specimens <- unique(hookMeasurements$specimen) keyHookMorphometrics <- as.data.frame(cbind(specimens, matrix(nrow=length(specimens), ncol=16))) names(keyHookMorphometrics) <- c("specimen", "firstHkL", "firstHkB", "firstHkA", "firstHkR", "lastHkL", "lastHkB", "lastHkA", "lastHkR", "maxL", "maxB", "maxA", "maxR", "posMaxL", "posMaxB", "posMaxA", "posMaxR") i=1 for (s in specimens) { sIndices <- which(hookMeasurements$specimen==s) nPlusOne <- max(hookMeasurements[sIndices,]$hook)+1 # number of hooks per row + corrective factor of 1 firstHkIndx <- which.min(hookMeasurements[sIndices,]$hook) keyHookMorphometrics$firstHkL[i] <- hookMeasurements$length[sIndices][firstHkIndx] keyHookMorphometrics$firstHkB[i] <- hookMeasurements$base[sIndices][firstHkIndx] keyHookMorphometrics$firstHkA[i] <- hookMeasurements$area[sIndices][firstHkIndx] keyHookMorphometrics$firstHkR[i] <- hookMeasurements$ratio[sIndices][firstHkIndx] lastHkIndx <- which.max(hookMeasurements[sIndices,]$hook) keyHookMorphometrics$lastHkL[i] <- hookMeasurements$length[sIndices][lastHkIndx] keyHookMorphometrics$lastHkB[i] <- hookMeasurements$base[sIndices][lastHkIndx] keyHookMorphometrics$lastHkA[i] <- hookMeasurements$area[sIndices][lastHkIndx] keyHookMorphometrics$lastHkR[i] <- hookMeasurements$ratio[sIndices][lastHkIndx] indxMaxL <- which.max(hookMeasurements[sIndices,]$length) indxMaxB <- which.max(hookMeasurements[sIndices,]$base) indxMaxA <- which.max(hookMeasurements[sIndices,]$area) indxMaxR <- which.max(hookMeasurements[sIndices,]$ratio) keyHookMorphometrics$maxL[i] <- hookMeasurements$length[sIndices][indxMaxL] keyHookMorphometrics$maxB[i] <- hookMeasurements$base[sIndices][indxMaxB] keyHookMorphometrics$maxA[i] <- hookMeasurements$area[sIndices][indxMaxA] keyHookMorphometrics$maxR[i] <- hookMeasurements$ratio[sIndices][indxMaxR] keyHookMorphometrics$posMaxL[i] <- (hookMeasurements$hook[sIndices][indxMaxL]/nPlusOne)*100 keyHookMorphometrics$posMaxB[i] <- (hookMeasurements$hook[sIndices][indxMaxB]/nPlusOne)*100 keyHookMorphometrics$posMaxA[i] <- (hookMeasurements$hook[sIndices][indxMaxA]/nPlusOne)*100 keyHookMorphometrics$posMaxR[i] <- (hookMeasurements$hook[sIndices][indxMaxR]/nPlusOne)*100 i=i+1 } keyHookMorphometrics$firstHkL <- as.numeric(keyHookMorphometrics$firstHkL) keyHookMorphometrics$firstHkB <- as.numeric(keyHookMorphometrics$firstHkB) keyHookMorphometrics$firstHkA <- as.numeric(keyHookMorphometrics$firstHkA) keyHookMorphometrics$firstHkR <- as.numeric(keyHookMorphometrics$firstHkR) keyHookMorphometrics$lastHkL <- as.numeric(keyHookMorphometrics$lastHkL) keyHookMorphometrics$lastHkB <- as.numeric(keyHookMorphometrics$lastHkB) keyHookMorphometrics$lastHkA <- as.numeric(keyHookMorphometrics$lastHkA) keyHookMorphometrics$lastHkR <- as.numeric(keyHookMorphometrics$lastHkR) keyHookMorphometrics$maxL <- as.numeric(keyHookMorphometrics$maxL) keyHookMorphometrics$maxB <- as.numeric(keyHookMorphometrics$maxB) keyHookMorphometrics$maxA <- as.numeric(keyHookMorphometrics$maxA) keyHookMorphometrics$maxR <- as.numeric(keyHookMorphometrics$maxR) keyHookMorphometrics$posMaxL <- as.numeric(keyHookMorphometrics$posMaxL) keyHookMorphometrics$posMaxB <- as.numeric(keyHookMorphometrics$posMaxB) keyHookMorphometrics$posMaxA <- as.numeric(keyHookMorphometrics$posMaxA) keyHookMorphometrics$posMaxR <- as.numeric(keyHookMorphometrics$posMaxR) return(keyHookMorphometrics) } # E. brayi brayi_female_key_morphometrics <- extractKeyHookMorphometrics(brayi_female_raw) brayi_male_key_morphometrics <- extractKeyHookMorphometrics(brayi_male_raw) summary(brayi_female_key_morphometrics) summary(brayi_male_key_morphometrics) # E. salmonis salmonis_female_dorsal_key_morphometrics <- extractKeyHookMorphometrics(salmonis_female_dorsal_raw) salmonis_female_ventral_key_morphometrics <- extractKeyHookMorphometrics(salmonis_female_ventral_raw) salmonis_male_dorsal_key_morphometrics <- extractKeyHookMorphometrics(salmonis_male_dorsal_raw) salmonis_male_ventral_key_morphometrics <- extractKeyHookMorphometrics(salmonis_male_ventral_raw) summary(salmonis_female_dorsal_key_morphometrics) summary(salmonis_female_ventral_key_morphometrics) summary(salmonis_male_dorsal_key_morphometrics) summary(salmonis_male_ventral_key_morphometrics) # E. truttae truttae_female_Carron_key_morphometrics <- extractKeyHookMorphometrics(truttae_female_Carron_raw) truttae_male_Carron_key_morphometrics <- extractKeyHookMorphometrics(truttae_male_Carron_raw) truttae_female_Drummore_key_morphometrics <- extractKeyHookMorphometrics(truttae_female_Drummore_raw) truttae_male_Drummore_key_morphometrics <- extractKeyHookMorphometrics(truttae_male_Drummore_raw) summary(truttae_female_Carron_key_morphometrics) summary(truttae_male_Carron_key_morphometrics) summary(truttae_female_Drummore_key_morphometrics) summary(truttae_male_Drummore_key_morphometrics) # E. bothniensis bothniensis_1_female_key_morphometrics <- extractKeyHookMorphometrics(bothniensis_1_female_raw) bothniensis_1_male_key_morphometrics <- extractKeyHookMorphometrics(bothniensis_1_male_raw) summary(bothniensis_1_female_key_morphometrics) summary(bothniensis_1_male_key_morphometrics) # E. 'bothniensis' bothniensis_2_female_key_morphometrics <- extractKeyHookMorphometrics(bothniensis_2_female_raw) summary(bothniensis_2_female_key_morphometrics) # E. gadi sp. A gadi_A_female_key_morphometrics <- extractKeyHookMorphometrics(gadi_A_female_raw) gadi_A_male_key_morphometrics <- extractKeyHookMorphometrics(gadi_A_male_raw) summary(gadi_A_female_key_morphometrics) summary(gadi_A_male_key_morphometrics) # E. gadi sp. B gadi_B_female_key_morphometrics <- extractKeyHookMorphometrics(gadi_B_female_raw) gadi_B_male_key_morphometrics <- extractKeyHookMorphometrics(gadi_B_male_raw) summary(gadi_B_female_key_morphometrics) summary(gadi_B_male_key_morphometrics) # E. gadi sp. I gadi_I_female_key_morphometrics <- extractKeyHookMorphometrics(gadi_I_female_raw) gadi_I_male_key_morphometrics <- extractKeyHookMorphometrics(gadi_I_male_raw) summary(gadi_I_female_key_morphometrics) summary(gadi_I_male_key_morphometrics)