Packages used

library(mdatools)
library(vegan)
library(ggplot2)
library(plotly)
library(gapminder)
library(readr)
library(knitr)
library(rsconnect)

Import and organize data

phloem            <- read_csv("phloem_avg_no_outliers.csv")
phloem$Familytype <- as.factor(phloem$Familytype)
phloem            <- subset(phloem, Familytype=="Full")

Second-derivative transformation

expl <- as.numeric(phloem$BV)
spec <- phloem[, 6:262]

derv <- prep.savgol(as.matrix(spec), width=15, porder=2, dorder=2)
fam  <- as.factor(phloem$Family)

NMDS + output

dist      <- vegdist(derv, method = 'euclidian')
peak.NMDS <- metaMDS(dist, distance = "euclidian", k = 2, try = 1000)
## Run 0 stress 0.02448102 
## Run 1 stress 0.02448102 
## ... New best solution
## ... Procrustes: rmse 6.376946e-05  max resid 0.0005055442 
## ... Similar to previous best
## Run 2 stress 0.02448102 
## ... Procrustes: rmse 1.653629e-05  max resid 0.0001032412 
## ... Similar to previous best
## Run 3 stress 0.02448101 
## ... New best solution
## ... Procrustes: rmse 6.12285e-05  max resid 0.0005046605 
## ... Similar to previous best
## Run 4 stress 0.02448098 
## ... New best solution
## ... Procrustes: rmse 4.934414e-05  max resid 0.0003204511 
## ... Similar to previous best
## Run 5 stress 0.02448097 
## ... New best solution
## ... Procrustes: rmse 1.043448e-05  max resid 7.176914e-05 
## ... Similar to previous best
## Run 6 stress 0.06186769 
## Run 7 stress 0.0879983 
## Run 8 stress 0.06186778 
## Run 9 stress 0.02448097 
## ... New best solution
## ... Procrustes: rmse 5.151479e-06  max resid 4.375675e-05 
## ... Similar to previous best
## Run 10 stress 0.06186767 
## Run 11 stress 0.02448099 
## ... Procrustes: rmse 1.531291e-05  max resid 9.777323e-05 
## ... Similar to previous best
## Run 12 stress 0.06650732 
## Run 13 stress 0.09353139 
## Run 14 stress 0.09500943 
## Run 15 stress 0.07906823 
## Run 16 stress 0.02448101 
## ... Procrustes: rmse 2.212642e-05  max resid 0.000140086 
## ... Similar to previous best
## Run 17 stress 0.06186769 
## Run 18 stress 0.06650726 
## Run 19 stress 0.02448102 
## ... Procrustes: rmse 4.183481e-05  max resid 0.0002685304 
## ... Similar to previous best
## Run 20 stress 0.02448103 
## ... Procrustes: rmse 3.978623e-05  max resid 0.0002786159 
## ... Similar to previous best
## *** Best solution repeated 5 times
stressplot(peak.NMDS)

peak.NMDS$stress
## [1] 0.02448097
adonis2(dist ~ fam,  permutations = 9999)
## Significance indicates that the families are significantly grouped by distance

scores <- as.data.frame(scores(peak.NMDS))
all.scores   <- cbind.data.frame(expl, fam, scores)
colnames(all.scores) <- c("IBV", 'family', 'NMDS1', 'NMDS2')

ibv    <- as.data.frame(unique(all.scores$IBV))

centroids <- aggregate(cbind(NMDS1, NMDS2) ~ family, all.scores, mean)
f         <- function(z)sd(z)/sqrt(length(z))
se        <- aggregate(cbind(se.x=NMDS1, se.y=NMDS2) ~ family, all.scores, f)
centroids <- merge(centroids, se, by = "family")
centroids <- cbind(centroids, ibv)
colnames(centroids) <- c("family", "NMDS1", "NMDS2", "se.1", "se.2", "IBV")
p <- ggplot(centroids, aes(NMDS1, NMDS2, label = family)) + geom_point(aes(col = IBV)) + theme_bw() +
  geom_point(data = centroids, aes(col=IBV)) +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=centroids, aes(ymin=NMDS2-se.2, ymax=NMDS2+se.2, col = IBV)) +
  geom_errorbarh(data=centroids, aes(xmin=NMDS1-se.1, xmax=NMDS1+se.1, col = IBV)) +
  geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label = family), alpha=.5) +
  labs(title = "Full sibs, phloem") + scale_color_viridis_c() 
## Warning in geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label =
## family), : Ignoring unknown aesthetics: label
ggplotly(p) ## interactive plot

Extract average lesion length per family

Add data to NMDS by family. Because these are just averages, we can’t run a PERMANOVA, so this is mostly for visualiation purposes.

lesion <- read_csv("M053_Gclav_lesion_data_tidy_YrAvg.csv")


avg.lesion <- aggregate(lesion_mm_avgforyear ~ family, lesion, mean)
f         <- function(z)sd(z)/sqrt(length(z))
se        <- aggregate(lesion_mm_avgforyear ~ family, lesion, f)

lesion.df <- merge(avg.lesion, se, by = 'family')

colnames(lesion.df) <- c("family", "mean", "se")


merged.df <- merge(centroids, lesion.df, by = 'family')

colnames(merged.df) <- c("family", "NMDS1", "NMDS2", "se.1", "se.2", "IBV", 
                         "lesion.avg", 'se.avg')

Is there a relation between spectra, family, and lesion? Note: cannot run a PERMANOVA on lesion lengths since they are averaged by family

lesion.p <- ggplot(merged.df, aes(NMDS1, NMDS2, label = family)) + geom_point(aes(col = lesion.avg)) + theme_bw() +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=merged.df, aes(ymin=NMDS2-se.2, ymax=NMDS2+se.2, col = lesion.avg), alpha = 0.8) +
  geom_errorbarh(data=merged.df, aes(xmin=NMDS1-se.1, xmax=NMDS1+se.1, col = lesion.avg), alpha = 0.8) + scale_color_viridis_c()

ggplotly(lesion.p)

To better visualize the relationship between breeding value, lesion, and spectra (i.e. NDMS distances),

merged.df$phenotype[merged.df$IBV < 0.5] <- 'Susceptible'
merged.df$phenotype[merged.df$IBV > 0.5] <- 'Resistant'


pheno.p <- ggplot(merged.df, aes(NMDS1, NMDS2, label = family)) + 
  geom_point(aes(col = IBV, size = lesion.avg)) + theme_bw() +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=merged.df, aes(ymin=NMDS2-se.2, ymax=NMDS2+se.2, col = IBV)) +
  geom_errorbarh(data=merged.df, aes(xmin=NMDS1-se.1, xmax=NMDS1+se.1, col = IBV)) + scale_color_viridis_b() +
  labs(title = "Phloem, Full sibs", subtitle = "size indicates average lesion size") + 
  geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label = family), alpha=.5)
## Warning in geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label =
## family), : Ignoring unknown aesthetics: label
plot(pheno.p)

ggplotly(pheno.p)

This is just for visualization purposes to see if there are any patterns we can parse out –

I combined all the time data, did some data analysis to see if there were any associations between proportion of time spent and family/breeding value/susceptibility.

Then I paired the average time proportions for each family in the NMDS. We bypass a lot of assumptions here so I wouldn’t use this as an official analysis (especially since I’m not sure if you wanted to combine the two experiments M055 and M059) but just as a visualization or a guiding tool.

m055 <- read_csv("M055_pine_olfactometer_data_FINAL.csv")
m059 <- read_csv("M059_pine_olfactometer_data_FINAL.csv")

m055_sub <- m055[, c(4, 7, 10, 11)]
m059_sub <- m059[, c(7, 12, 15, 16)]

time_df <- rbind(m055_sub, m059_sub)
colnames(time_df) <- c("time", "family", "phenotype", "IBV")

By family

avg.time <- aggregate(time ~ family*IBV, time_df, mean)
f         <- function(z)sd(z)/sqrt(length(z))
se        <- aggregate(time ~ family, time_df, f)

time.df <- merge(avg.time, se, by = 'family')

colnames(time.df) <- c("family", "IBV", "mean", "se")
time.df$family <- as.factor(time.df$family)

time.df$phenotype[time.df$IBV < 0.5] <- 'Susceptible'
time.df$phenotype[time.df$IBV > 0.5] <- 'Resistant'


ggplot(time.df, aes(x= family, y = mean )) + geom_bar(stat='identity') +
  geom_errorbar(aes(ymin=mean - se, ymax = mean + se), size = .3, width = .3)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

By phenotype

avg.time.p <- aggregate(time ~ phenotype, time_df, mean)
f         <- function(z)sd(z)/sqrt(length(z))
se        <- aggregate(time ~ phenotype, time_df, f)

time.df.p <- merge(avg.time.p, se, by = 'phenotype')

colnames(time.df.p) <- c("phenotype", "mean", "se")


ggplot(time.df.p, aes(x= phenotype, y = mean )) + geom_bar(stat='identity') +
  geom_errorbar(aes(ymin=mean - se, ymax = mean + se), size = .3, width = .3)

Add average beetle time for average family centroids

time.w.centroid <- merge(avg.time[,c(1,3)], centroids, by ='family')



time.plot <- ggplot(time.w.centroid, aes(NMDS1, NMDS2, label = family)) + 
  geom_point(aes(col = IBV, size = time), show.legend = T) + theme_bw() +
  geom_text(nudge_x = 0.005, nudge_y = 0.005) +
  geom_errorbar(data=time.w.centroid, aes(ymin=NMDS2-se.2, 
                                          ymax=NMDS2+se.2, col = IBV)) +
  geom_errorbarh(data=time.w.centroid, aes(xmin=NMDS1-se.1, 
                                           xmax=NMDS1+se.1, col = IBV)) + scale_color_viridis_b()  + 
  geom_point(data = all.scores, aes(NMDS1, NMDS2,col = IBV, label = family))
## Warning in geom_point(data = all.scores, aes(NMDS1, NMDS2, col = IBV, label =
## family)): Ignoring unknown aesthetics: label
plot(time.plot)

ggplotly(time.plot)