Ipomopsis temperature experiments
Carrie Wu, John Powers
Feb 26 2020
#setwd("/media/C/GCMSsolution/Data/heather Ipo/")
#data18 <- read.shimadzu("./Temperature expt CWu 2018/CWu2018.txt")
#data19 <- read.shimadzu("./Temperature expt CWu 2019/CWu2019..txt")
#ipo.data <- rbind(data18, data19)
#save(ipo.data, file="CWu_2018_2019.Rdata")
load("CWu_2018_2019.Rdata")
ipo.all <- dcast(ipo.data, Filename~Name, sum, value.var="Area")
rownames(ipo.all) <- ipo.all[,1]
ipo.all[,1] <- NULL
ipo.cut <- ipo.all[,colSums(ipo.all)>4e6]
write(rownames(ipo.all), "CWu_samples.txt")
ipos <- read.csv("CWu_samples.csv", colClasses=list(RunDate="Date"))Overview - no metadata
Inventory
with(ipos, table(Year, Type)) Type
Year A Air Bake bb Blank F G GNA H I J NatArea OTC SC Schiedea T
2018 35 27 3 13 7 3 3 0 2 3 3 48 0 0 0 34
2019 29 21 0 0 10 0 0 19 0 0 0 0 53 27 18 27with(ipos, table(Type, DNLeaf)) DNLeaf
Type day leaf night
A 35 12 8 9
Air 31 9 0 8
Bake 3 0 0 0
bb 13 0 0 0
Blank 17 0 0 0
F 3 0 0 0
G 3 0 0 0
GNA 7 3 9 0
H 2 0 0 0
I 3 0 0 0
J 3 0 0 0
NatArea 48 0 0 0
OTC 10 15 19 9
SC 0 10 9 8
Schiedea 18 0 0 0
T 34 10 7 10NMDS
set.seed(1)
nmds.ipo <- metaMDS(decostand(ipo.cut, "hellinger"), dist="bray", autotransform = FALSE, try=5, trymax=5) Run 0 stress 0.1441741
Run 1 stress 0.1600284
Run 2 stress 0.1949266
Run 3 stress 0.1801057
Run 4 stress 0.2101388
Run 5 stress 0.2009994
*** No convergence -- monoMDS stopping criteria:
5: stress ratio > sratmaxpar(bg="grey40")
ordiplot(nmds.ipo, type = "n")
ipos$rs <- rowSums(ipo.all)
ipos$nameBlank <- ipos$Type %in% c("Blank","Bake")
ipos$nameAir <- ipos$Type == "Air"
thisyear <- ipos$Year == 2019
#with(ipos, points(nmds.ipo, display="sites", col=viridis(200)[round(200*sqrt(rs/max(rs)))], pch=as.integer(isblank)+1))
with(ipos, points(nmds.ipo, display="sites", col=rainbow(nlevels(ipos$Type))[ipos$Type], pch=as.integer(nameAir*2+nameBlank)+1))
#with(ipos, points(nmds.ipo, display="sites", col=rainbow(nlevels(ipos$Type))[ipos$Type], pch=ipos$Year-2017))
legend("topright", levels(ipos$Type), fill=rainbow(nlevels(ipos$Type)), cex=0.8)
Kmeans clustering
k <- 8
set.seed(1)
km <- kmeans(decostand(ipo.cut, "log"), k, nstart=3)
ipos$Cluster <- km$cluster
ipos$kBlank <- kblank <- km$cluster %in% c(4)
ipos$Mixup <- ipos$nameBlank != ipos$kBlank
with(ipos, table(kBlank, nameBlank)) nameBlank
kBlank FALSE TRUE
FALSE 354 0
TRUE 11 20ordiplot(nmds.ipo, type = "n")
points(nmds.ipo, display="sites", col=ifelse(kblank, "black", rainbow(k)[km$cluster]), pch=with(ipos,as.integer(nameAir*2+nameBlank))+1)
text(nmds.ipo, display="species",col="grey80",labels=ifelse(colSums(ipo.cut)>4e8, colnames(ipo.cut), ""))
text(nmds.ipo, display="sites", col=with(ipos, as.integer(thisyear)*4+as.integer(nameBlank)*2+as.integer(nameAir)+1), labels=as.character(km$cluster))
CAP - blanks and years
ipo.cap <- capscale(ipo.cut ~ as.factor(ipos$kBlank) * as.factor(thisyear), distance="bray", metaMDSdist = F)
plot(ipo.cap, type="n")
points(ipo.cap, display="sites", col=ifelse(ipos$kBlank, "black", rainbow(k)[ipos$Cluster]), pch=as.integer(ipos$nameAir*2+ipos$nameBlank)+1) 
#View(ipo.cap$CCA$v)CAP - DNLeaf
dnl <- ipos$DNLeaf!="" & ipos$Type !="Air"
ipo.cap.dnl <- capscale(ipo.cut[dnl,] ~ ipos[dnl,"DNLeaf"], distance="bray", metaMDSdist = F)
plot(ipo.cap.dnl, type="n")
points(ipo.cap.dnl, display="sites", pch=as.integer(ipos[dnl,"DNLeaf"])-1, col=as.integer(factor(ipos[dnl,"Type"])))
legend("bottomleft", levels(factor(ipos[dnl,"Type"])), fill=1:nlevels(ipos$Type), cex=0.8)
legend("left", levels(ipos[dnl,"DNLeaf"]), pch=1:nlevels(ipos[dnl,"DNLeaf"])-1, cex=0.8)
#View(ipo.cap.dnl$CCA$v)Read metadata
files <- read_tsv("CWu_samples_text2.csv")
perc <- read_tsv("Percivals2018_2019.csv")
otc <- read_tsv("OTC2018_2019.csv")
perc2 <- left_join(perc,files, by=c("Date","DN","Tissue","Sample2"), suffix=c("",".f"))
otc2 <- left_join(otc,files, by=c("Date","Site","OI","DN","Tissue","Sample2"), suffix=c("",".f"))
meta <- bind_rows(otc=otc2, perc=perc2, .id="expt") %>%
mutate_if(is.character, as.factor) %>%
mutate(Equil=Pump-Bag, Duration=End-Pump, Total=End-Bag)
meta %>% write.csv("CWu_meta.csv")
with(meta, table(paste(Temp, Tissue, OI), DN, Expt)) , , Expt = otc
DN
day night
cool air NA 0 0
cool flower NA 0 0
cool leaf NA 0 0
NA air air 15 3
NA flower in 67 9
NA flower out 49 8
NA leaf in 21 0
NA leaf out 18 0
warm air NA 0 0
warm flower NA 0 0
warm leaf NA 0 0
, , Expt = perc
DN
day night
cool air NA 8 7
cool flower NA 29 27
cool leaf NA 7 0
NA air air 0 0
NA flower in 0 0
NA flower out 0 0
NA leaf in 0 0
NA leaf out 0 0
warm air NA 8 6
warm flower NA 26 25
warm leaf NA 7 0Missing files
#perc2 %>% filter(is.na(Filename)) %>% select("n") %>% write.csv("perc_nomatch.csv")
#otc2 %>% filter(is.na(Filename)) %>% select("n") %>% write.csv("otc_nomatch.csv")
meta %>% filter(is.na(Filename) & nomatch!="NR") %>% select(1:24) %>% kable| expt | n | Year | Date | DN | Site | PlantType | Tissue | OI | Sample | Sample2 | nomatch | Bag | Pump | End | NumFlrs | NumBuds | NumLeaf | PumpID | SoilVWC | SoilDepth | Comments | Collector | LeafDryWtGrams |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| otc | 184 | 2018 | 2018-07-25 | day | GNA | agg | flower | in | KM10-B1 | KM10B1 | 0 | 10:13:00 | 10:43:00 | 10:58:00 | 1 | 0 | NA | DRC5 | 3.7 | 1.171 | NA | Rachel | NA |
| otc | 37 | 2019 | 2019-07-31 | day | GNA | agg | flower | in | OTC-W1 | W1 | 0 | 11:53:00 | 12:26:00 | 12:41:00 | 1 | 0 | 0 | CW2 | 8.6 | 1.258 | NA | Carrie | NA |
| otc | 83 | 2019 | 2019-08-02 | day | SC | ten | leaf | in | OTC-A2-leaf | A2 | 0 | 12:50:00 | 13:20:00 | 13:35:00 | NA | NA | 1 | . | NA | NA | NA | Carrie | 0.01030 |
| otc | 84 | 2019 | 2019-08-02 | day | SC | ten | leaf | in | OTC-C3-leaf | C3 | 0 | 13:08:00 | 13:41:00 | 13:56:00 | NA | NA | 1 | CW2 | NA | NA | NA | Carrie | 0.00829 |
| otc | 65 | 2019 | 2019-08-02 | day | SC | ten | leaf | in | OTC-E1-leaf | E1 | 0 | 09:42:00 | 10:12:00 | 10:27:00 | NA | NA | 1 | CW2 | 6.1 | 1.191 | NA | Carrie | 0.00849 |
| otc | 5 | 2019 | 2019-08-02 | night | SC | ten | flower | out | SC9 | 9 | 0 | 20:21:00 | 20:51:00 | 21:06:00 | 1 | 0 | 0 | CW2 | NA | NA | NA | Diane | NA |
| perc | 47 | 2019 | 2019-07-25 | night | NA | agg | flower | NA | A79 | A79 | nm | 21:02:00 | 21:35:00 | 21:50:00 | 1 | 0 | NA | CW2 | NA | NA | NA | Carrie | NA |
Missing metadata
#files %>% filter(!(Filename %in% meta$Filename)) %>% select("n") %>% write.csv("CWu_samples_nomatch.csv")
files %>% filter(!(Filename %in% meta$Filename) & Expt!="") %>% select(c("Filename", "Date","Site","OI","DN","Tissue","Sample2")) %>% kable| Filename | Date | Site | OI | DN | Tissue | Sample2 |
|---|---|---|---|---|---|---|
| AirlastJuly25_08312018.qgd | 2018-07-25 | GNA | air | day | air | Airlast |
| NatAreaXUP01July25_08312018.qgd | 2018-07-25 | GNA | in | day | flower | XUP01 |
| GNA5redo_leaf_190805_8112019_17.qgd | 2019-08-05 | GNA | out | day | leaf | 5redo |
| A188DredoAug16renamed_09022018.qgd | 2018-08-16 | NA | NA | day | flower | A188redo |
| T196DredoAug16_09022018.qgd | 2018-08-16 | NA | NA | day | flower | T196redo |
| A38redo_day_190725_872019_03.qgd | 2019-07-25 | NA | NA | day | flower | A38redo |
| A91_leaf_190725_882019_02.qgd | 2019-07-25 | NA | NA | day | leaf | A91 |
Times
ggplot(meta, aes(x=paste(Expt, DN,Year), y=Bag, color=DN)) +
geom_boxplot() + geom_point()+
theme_minimal() + theme(axis.text.x=element_text(angle=90))+
ggtitle("Sampling time of day")
ggplot(meta, aes(x=paste(Expt, DN,Year), y=Total/60, color=DN)) +
geom_boxplot() + geom_point()+
theme_minimal() + theme(axis.text.x=element_text(angle=90)) +
ggtitle("Total sampling duration")+
ylab("Equilibration plus pumping time (min)")
OTC Soil moisture
ggplot(meta %>% filter(Expt=="otc", OI !="air"), aes(x=paste(Year, Site, OI), y=SoilVWC, color=OI)) +
geom_boxplot()+#geom_violin(draw_quantiles=c(0.25, 0.5, 0.75)) +
geom_point()+
theme_minimal() + theme(axis.text.x=element_text(angle=90)) +
scale_color_brewer(palette = "Set1")+
ggtitle("Effects of OTCs on Soil Moisture")+
ylab("Soil Volumetric Water Content (%)") +
xlab("Year, in or out of the OTC")
# Requirements: .csv files named with the site name or serial number....
read.format.hobo <- function(dataFile, skiprows,...){
df4 <- read.csv(dataFile, header = TRUE, skip = skiprows, sep = ",")
nametemp <- sub(".*temp.*", "TempC", colnames(df4), ignore.case = TRUE)
colnames(df4) <- sub(".*date.*", "DateTime", nametemp, ignore.case = TRUE)
df3<- df4[,c("DateTime", "TempC")] nameSep<-unlist(strsplit(dataFile, "[_.]")) # unlist/split filename.
df1<- cbind(SiteName = paste(nameSep[1]), df2) # paste SiteName into 1st col.
df<- na.omit(df1) }
setwd("C:/Temperature/RawCSVfiles")
Files <- Sys.glob("*.csv") # list all .csv files in working directory.
out <- lapply(Files, FUN = read.format.hobo, skiprows = 1) # apply custom function to each file.
output <- ldply(out, data.frame) # changes 'out' list to dataframe.Filtering criteria
sv <- meta %>% drop_na(Filename)
ipo.gc <- ipo.all[sv$Filename,]
library(bouquet)
ipo.data.cut <- ipo.data[ipo.data$Filename %in% sv$Filename,]
ipo.data.cut$Filename <- sv$Filename[match(ipo.data.cut$Filename, sv$Filename)]
ipo.data.cut <- ipo.data.cut[ipo.data.cut$Name != "", ]
ipo.data.cut$Name <- droplevels(ipo.data.cut$Name)
longdata <- load_longdata(ipo.data.cut, sample="Filename", RT="Ret.Time", name="Name", area="Area", match = "SI", maxmatch=100)
metadata <- load_metadata(sv, date="Date", sample="Filename", group="PlantType", type="Tissue")
levels(metadata$type) <- c("ambient","floral","leaf")
vol.all <- make_sampletable(longdata)
chems <- make_chemtable(longdata, metadata)
### Flower filtering
chemsf <- chems %>%
filter_RT(2, 17) %>%
filter_match(0.8) %>%
filter_freq(0.1, group = FALSE) %>%
filter_contaminant(cont.list = c("Caprolactam")) %>%
filter_area(min_maximum = 4e6) %>%
filter_ambient_ratio(vol.all, metadata, ratio = 4) %>%
filter_ambient_ttest(vol.all, metadata, alpha = 0.05, adjust = "fdr") %>%
combine_filters()
chemsf$filter_final <- with(chemsf, filter_RT == "OK" & filter_ambient_ratio == "OK" & filters_passed > 5 & filter_contaminant == "OK")
### Leaf filtering
meta.leaf <- metadata
levels(meta.leaf$type) <- c("ambient","flower","floral") #trick bouquet to think the leaf samples are floral TODO: hacky!
chems.leaf <- chems %>%
filter_RT(2, 17) %>%
filter_match(0.8) %>%
filter_freq(0.1, group = FALSE) %>%
filter_contaminant(cont.list = c("Caprolactam")) %>%
filter_area(min_maximum = 4e6) %>%
filter_ambient_ratio(vol.all, meta.leaf, ratio = 4) %>%
filter_ambient_ttest(vol.all, meta.leaf, alpha = 0.05, adjust = "fdr") %>%
combine_filters()
chems.leaf$filter_final <- with(chems.leaf, filter_RT == "OK" & filter_ambient_ratio == "OK" & filters_passed > 5 & filter_contaminant == "OK")
chemsf <- chemsf %>% mutate(filter_ambient_leaf_ratio = chems.leaf$filter_ambient_ratio,
filter_ambient_leaf_ttest = chems.leaf$filter_ambient_ttest)
#chemsf$leaf_not_flower <- with(chemsf, filter_ambient_leaf_ratio=="OK" & filter_ambient_leaf_ttest=="OK" & !(filter_ambient_ratio=="OK" & filter_ambient_ttest=="OK"))
### Manual filtering
lapply(chemsf[grepl("filter", names(chemsf))], table)
write.csv(chemsf, "chemsf_ipo.csv")
#annotated this in a google sheet
chemsf2 <- read.csv("Ipo volatile compounds - chemsf_ipo.csv.csv")
chemsf3 <- chemsf %>% filter_contaminant(cont.list=as.character(chemsf2[chemsf2$verdict %in% c("ambi","bis14","cont","rare"),"name"]))
keep <- as.character(chemsf2[chemsf2$verdict %in% c("fvoc","keep","merge","keep-leaf"),"name"])
#keep_leaf <- as.character(chemsf2[chemsf2$verdict =="keep_leaf","name"]) # > 4 | chemsf3$name %in% keep_leaf)
chemsf3$filter_final <- with(chemsf3, filters_passed > 4 & filter_contaminant == "OK" & name %in% keep)#keeps with passed<=4 still excluded
chemsf3$filter_final_leaf <- chemsf3$filter_contaminant == "OK" & chems.leaf$filter_final & chemsf2$shortname!=""#four big alkanes
#shorten names
colnames(vol.all) <- chemsf3$name2 <- ifelse(chemsf2$shortname=="", as.character(chemsf2$name), as.character(chemsf2$shortname))
#TODO merge the duplicate name2s - the following doesn't work
#longdata3 <- longdata %>% inner_join(chemsf3 %>% select(c("name", "name2")), by = "name") %>%
#mutate(name = factor(name2))
#vol.all3 <- make_sampletable(longdata3)
vol <- prune_sampletable(vol.all, chemsf3, metadata)
files_exclude <- c("T196DAug16_09012018.qgd", "A188DAug16_09012018.qgd", "OTC_Z3_190805_8102019_24.qgd", "T139_day_190725_862019_06.qgd")#no filtered volatiles
svf <- metadata[metadata$type == "floral" & !(metadata$sample %in% files_exclude),]
svf$Year <- factor(svf$Year)
svf <- svf %>% mutate_if(is.factor, droplevels)
svf$Treat <- factor(paste(svf$OI,svf$Temp))
levels(svf$Treat) <- c("warm or in","cool or out", "warm or in","cool or out")
vol <- vol[!(rownames(vol) %in% files_exclude) ,]
vol <- vol / 0.75 / 1 #30 min of equilibration plus 15 minpumping, one flower
#TODO fix svf$Total data entry errors
vol.leaf <- prune_sampletable(vol.all, chemsf3 %>% mutate(filter_final=filter_final_leaf), meta.leaf)
files_exclude_leaf <- c("CWH1LeafJuly27_09022018.qgd","NatArea69leafrenamedJuly25_08312018.qgd", "T164_leaf_190731_882019_18.qgd","OTC_Z2_leaf_190805_8102019_23.qgd")#fist two lack leaf mass, secodn two have no filtered volatiles
svl <- metadata[metadata$type == "leaf" & !(metadata$sample %in% files_exclude_leaf),]
svl <- svl %>% mutate_if(is.factor, droplevels)
svl$Treat <- factor(paste(svl$OI,svl$Temp))
levels(svl$Treat) <- c("warm or in","cool or out", "warm or in","cool or out")
vol.leaf <- vol.leaf[!(rownames(vol.leaf) %in% files_exclude_leaf) ,]
svl$LeafDryWtGrams[c(44,3)] <- svl$LeafDryWtGrams[c(44,3)]/10
vol.leaf <- vol.leaf / 0.75 / svl$LeafDryWtGrams #area per hour per dry gram
#save.image("ipo_workspace.Rdata")load("ipo_workspace.Rdata")Flower heatmap
#cairo_pdf("ipo_heatmap.pdf", width=12, height=12)
ph <- pheatmap(as.matrix(t(vol))^(1/3),
cluster_cols=T, show_colnames=F,
clustering_method="mcquitty", clustering_distance_rows="correlation",
clustering_distance_cols=vegdist(vol, method = "bray"),
clustering_callback = function(hc, ...){dendsort(hc, type="average")},
scale="none", color=inferno(512),
annotation_col = data.frame(svf %>% select("Year", "PlantType", "expt", "Treat","DN"), row.names=rownames(vol)),
fontsize = 10, border_color = NA, legend=F, annotation_legend=T, cutree_rows=6
)
#dev.off()NMDS of all filtered samples
set.seed(1)
nmds.vol <- metaMDS(sqrt(vol), dist="bray", autotransform = FALSE) Run 0 stress 0.2231919
Run 1 stress 0.2247714
Run 2 stress 0.2240661
Run 3 stress 0.2255607
Run 4 stress 0.2330845
Run 5 stress 0.2262488
Run 6 stress 0.2244393
Run 7 stress 0.2270576
Run 8 stress 0.2285393
Run 9 stress 0.2298174
Run 10 stress 0.2322534
Run 11 stress 0.2274148
Run 12 stress 0.2330405
Run 13 stress 0.2274875
Run 14 stress 0.2269826
Run 15 stress 0.2279473
Run 16 stress 0.2332582
Run 17 stress 0.2310793
Run 18 stress 0.2291831
Run 19 stress 0.2383751
Run 20 stress 0.2262468
*** No convergence -- monoMDS stopping criteria:
1: no. of iterations >= maxit
19: stress ratio > sratmaxTEY <- with(svf, factor(paste(PlantType,expt,Year)))
TEYcols <- brewer.pal(8, "Paired")[c(1:4,8,5:6)]
ordiplot(nmds.vol, type = "n")
text(nmds.vol, display="species", cex=0.7, col="grey20")
points(nmds.vol, display="sites", col=TEYcols[TEY], pch=c(13,19)[svf$DN])
#points(nmds.vol, display="sites", col=viridis(100)[sqrt(rowSums(vol))/100], pch=19)
legend("bottomleft", levels(TEY), fill=TEYcols, cex=0.9)
OTCs vs. Ambient
CAP analysis of each OTC site and year
#inventory of samples
with(svf %>% filter(expt=="otc"), table(paste(PlantType, DN, OI),Year, Site)) , , Site = bb
Year
2018 2019
agg day in 13 0
agg day out 13 0
ten day in 0 0
ten day out 0 0
ten night in 0 0
ten night out 0 0
, , Site = GNA
Year
2018 2019
agg day in 30 12
agg day out 16 10
ten day in 0 0
ten day out 0 0
ten night in 0 0
ten night out 0 0
, , Site = SC
Year
2018 2019
agg day in 0 0
agg day out 0 0
ten day in 0 11
ten day out 0 10
ten night in 0 9
ten night out 0 8plot.otc <- function(svf.otc, plot.DN=F, plot.chems=T) {
vol.otc <- vol[svf.otc$sample,]
cap.otc <- capscale(sqrt(vol.otc) ~ OI, distance="bray", data=svf.otc)
if(plot.DN) cap.otc <- capscale(sqrt(vol.otc) ~ DN+OI, distance="bray", data=svf.otc)
summ.cap <- summary(cap.otc)
prop.expl <- summ.cap$cont$importance
ax.labs <- paste0(colnames(prop.expl)[1:2], " (",round(100*prop.expl[2,1:2]),"% explained)")
anova.cap <- anova(cap.otc, by="term")
plot(cap.otc, type="n", xlim=range(cap.otc$CCA$wa[,1]), xlab=ax.labs[1], ylab=ax.labs[2])
with(svf.otc, title(paste("OTC warming:","I.",unique(PlantType), "at", unique(Site),"site,",unique(Year),
do.call(paste0, list(unique(DN), collapse="/")))))
legend("topleft", title=paste("OTC warming","P =",anova.cap["OI","Pr(>F)"]),
levels(svf.otc$OI), fill=2:3)
if(plot.DN) legend("topright", title=paste("Time","P =",anova.cap["DN","Pr(>F)"]),
levels(svf.otc$DN), pch=c(1,19))
points(cap.otc, display="sites", col=c(2:3)[svf.otc$OI], pch=c(1,19)[svf.otc$DN])
text(cap.otc, display="cn")
if(plot.chems) {
if(plot.DN) text(cap.otc, display="species", cex=0.8, srt=45,
col=ifelse(sqrt(cap.otc$CCA$v[,1]^2 + cap.otc$CCA$v[,2]^2)>0.25, alpha("purple",1),alpha("purple",0)))
else text(cap.otc, display="species", cex=0.8, srt=45,
col=ifelse(cap.otc$CCA$v[,1]>0.1, alpha("purple",1),alpha("purple",0)))
}
return(cap.otc)
}
cap.otcbb18 <- plot.otc(svf %>% filter(expt=="otc",Year=="2018",Site=="bb",PlantType=="agg"))
cap.otcGNA18 <- plot.otc(svf %>% filter(expt=="otc",Year=="2018",Site=="GNA",PlantType=="agg"))
cap.otcGNA19 <- plot.otc(svf %>% filter(expt=="otc",Year=="2019",Site=="GNA",PlantType=="agg"))
cap.otcSC19 <- plot.otc(svf %>% filter(expt=="otc",Year=="2019",Site=="SC",PlantType=="ten"),plot.DN=T)
#cap.otc.all <- list(bb18=cap.otcbb18, GNA18=cap.otcGNA18, GNA19=cap.otcGNA19, SC19=cap.otcSC19)CAP of all I. aggregata OTC datasets
#Combined I. aggregata OTC CAP
svf.otc <- svf %>% filter(expt=="otc", PlantType=="agg")
vol.otc <- vol[svf.otc$sample,]
cap.otc.comb <- capscale(sqrt(vol.otc) ~ Year+Site+OI, distance="bray", data=svf.otc)
summ.cap <- summary(cap.otc.comb)
prop.expl <- summ.cap$cont$importance
ax.labs <- paste0(colnames(prop.expl)[1:2], " (",round(100*prop.expl[2,1:2]),"% explained)")
anova.cap <- anova(cap.otc, by="term")
plot(cap.otc.comb, type="n",xlim=range(cap.otc.comb$CCA$wa[,1]), xlab=ax.labs[1], ylab=ax.labs[2])
title("OTC warming: I. agg at GNA or bb sites, 2018/2019")
legend("topleft", title=paste("OTC warming","P =",anova.cap["OI","Pr(>F)"]), levels(svf.otc$OI), fill=2:3)
legend("topright", title=paste("Year","P =",anova.cap["Year","Pr(>F)"]), levels(svf.otc$Year), pch=c(15,19))
legend("bottomleft", title=paste("Site","P =",anova.cap["Site","Pr(>F)"]), levels(svf.otc$Site)[1:2])
text(cap.otc.comb, display="species", cex=0.7, srt=45,
col=ifelse(sqrt(cap.otc.comb$CCA$v[,1]^2 + cap.otc.comb$CCA$v[,2]^2)>0.25, alpha("purple",0.7),alpha("purple",0)))
points(cap.otc.comb, display="sites", col=c(2:3)[svf.otc$OI], pch=c(15,19)[svf.otc$Year])
text(cap.otc.comb, display="cn")
cap.otc.comb Call: capscale(formula = sqrt(vol.otc) ~ Year + Site + OI, data =
svf.otc, distance = "bray")
Inertia Proportion Rank
Total 22.38982 1.00000
Constrained 1.86625 0.08335 3
Unconstrained 25.32007 1.13087 41
Imaginary -4.79649 -0.21423 52
Inertia is squared Bray distance
Species scores projected from 'sqrt' 'vol.otc'
Eigenvalues for constrained axes:
CAP1 CAP2 CAP3
1.2424 0.3193 0.3045
Eigenvalues for unconstrained axes:
MDS1 MDS2 MDS3 MDS4 MDS5 MDS6 MDS7 MDS8
4.464 3.288 2.415 1.927 1.471 1.276 1.202 0.944
(Showing 8 of 41 unconstrained eigenvalues)anova(cap.otc.comb, by="term") Permutation test for capscale under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(vol.otc) ~ Year + Site + OI, data = svf.otc, distance = "bray")
Df SumOfSqs F Pr(>F)
Year 1 1.1512 4.0920 0.001 ***
Site 1 0.3539 1.2578 0.235
OI 1 0.3612 1.2838 0.188
Residual 90 25.3201
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Percivals warm vs. cool
CAP analysis of each Percival dataset
#inventory of samples
with(svf %>% filter(expt=="perc"), table(paste(PlantType, DN, Temp),Year)) Year
2018 2019
agg day cool 9 5
agg day warm 8 6
agg night cool 9 4
agg night warm 7 5
ten day cool 7 5
ten day warm 8 4
ten night cool 9 5
ten night warm 8 5plot.perc <- function(svf.perc, plot.chems=T) {
vol.perc <- vol[svf.perc$sample,]
cap.perc <- capscale(sqrt(vol.perc) ~ DN+Temp, distance="bray", data=svf.perc)
summ.cap <- summary(cap.perc)
prop.expl <- summ.cap$cont$importance[2,c("CAP1","CAP2")]
ax.labs <- paste0("CAP", 1:2, " (",round(100*prop.expl),"% explained)")
anova.cap <- anova(cap.perc, by="term")
plot(cap.perc, type="n", xlim=range(cap.perc$CCA$wa[,1]), xlab=ax.labs[1], ylab=ax.labs[2])
with(svf.perc, title(paste0("Percivals: I. ", unique(PlantType)," in ", unique(Year))))
#, " (",round(100*summ.cap$constr.chi / summ.cap$tot.chi),"% explained)")
legend("topleft", title=paste("Temperature","P =",anova.cap["Temp","Pr(>F)"]),
levels(svf.perc$Temp), fill=c(4,2))
legend("topright", title=paste("Time","P =",anova.cap["DN","Pr(>F)"]),
levels(svf.perc$DN), pch=c(1,19))
if(plot.chems) text(cap.perc, display="species", cex=0.7, srt=45,
col=ifelse(sqrt(cap.perc$CCA$v[,1]^2 + cap.perc$CCA$v[,2]^2)>0.25, alpha("purple",0.7),alpha("purple",0)))
points(cap.perc, display="sites", col=c(4,2)[svf.perc$Temp], pch=c(1,19)[svf.perc$DN])
text(cap.perc, display="cn")
return(cap.perc)
}
cap.perc.agg18 <- plot.perc(svf %>% filter(expt=="perc",Year=="2018",PlantType=="agg"))
cap.perc.agg19 <- plot.perc(svf %>% filter(expt=="perc",Year=="2019",PlantType=="agg"))
cap.perc.ten18 <- plot.perc(svf %>% filter(expt=="perc",Year=="2018",PlantType=="ten"))
cap.perc.ten19 <- plot.perc(svf %>% filter(expt=="perc",Year=="2019",PlantType=="ten"))
#cap.perc.all <- list(agg18=cap.perc.agg18, agg19=cap.perc.agg19, ten18=cap.perc.ten18, ten19=cap.perc.ten19)
#lapply(cap.perc.all, anova, by="term")CAP of Percival dataset combined for each species
plot.perc.comb <- function(svf.perc, plot.chems=F) {
vol.perc <- vol[svf.perc$sample,]
cap.perc.comb <- capscale(sqrt(vol.perc) ~ Year+DN+Temp, distance="bray", data=svf.perc)
summ.cap <- summary(cap.perc.comb)
prop.expl <- summ.cap$cont$importance[2,c("CAP1","CAP2")]
ax.labs <- paste0("CAP", 1:2, " (",round(100*prop.expl),"% explained)")
anova.cap <- anova(cap.perc.comb, by="term")
plot(cap.perc.comb, type="n", xlim=range(cap.perc.comb$CCA$wa[,1]), xlab=ax.labs[1], ylab=ax.labs[2])
title(paste("Percivals: I.",unique(svf.perc$PlantType), "2018/2019"))
legend("bottomleft", title=paste("Temperature","P =",anova.cap["Temp","Pr(>F)"]), levels(svf.perc$Temp), fill=c(4,2))
legend("bottomright", title=paste("Time","P =",anova.cap["DN","Pr(>F)"]), levels(svf.perc$DN), fil=c(0,1))
legend("topright", title=paste("Year","P =",anova.cap["Year","Pr(>F)"]), levels(svf.perc$Year), pch=c(19,15))
points(cap.perc.comb, display="sites", col=c(4,2)[svf.perc$Temp], pch=c(1,0,19,15)[factor(paste(svf.perc$DN, svf.perc$Year))])
text(cap.perc.comb, display="cn")
text(cap.perc.comb, display="species", cex=0.7, srt=45, col=ifelse(sqrt(cap.perc.comb$CCA$v[,1]^2 + cap.perc.comb$CCA$v[,2]^2)>0.2, alpha("purple",1),alpha("purple",0)))
return(cap.perc.comb)
}
svf.perc.agg <- plot.perc.comb(svf %>% filter(expt=="perc", PlantType=="agg"))
svf.perc.ten <- plot.perc.comb(svf %>% filter(expt=="perc", PlantType=="ten"))
#svf.perc.both <- plot.perc.comb(svf %>% filter(expt=="perc"))Barplots
vols.plot <- chemsf3 %>% filter(filter_final, freq.floral>0.21) %>% pull("name2") %>% as.character
vols.plot[vols.plot=="b-caryophyllene"] <- "b-caryophyllene.1"#c("b-ocimene","a-pinene","a-farnesene","b-caryophyllene.1","petasitene","b-myrcene","longifolene","tridec-1-ene ","benzophenone")
vol.long <- gather(cbind(sample=rownames(vol[,vols.plot]), total=rowSums(vol), vol[,vols.plot]), key="chem",value="area", 2:(length(vols.plot)+2), factor_key=T)
vol.plot <- left_join(vol.long, svf)
specDNcols <- c( "#943390","#FF5056", "#A894FF","#FC8C97")
#Treatcols <- c("#DA9E18","#6B744F")
DNcols <- c("grey90","grey20")
ggplot(vol.plot, aes(x=paste(PlantType,expt),y=area, fill=DN, color=paste(PlantType,Treat))) +
geom_boxplot(size=1, position=position_dodge2(preserve="single")) +
facet_wrap(vars(chem), scales="free") +
scale_y_sqrt("Peak area per flower per hour") +
scale_color_manual("Species, Treatment", values=specDNcols) +
scale_fill_manual("Time of day",values=DNcols) +
theme_minimal() + theme(axis.text.x=element_text(angle=20, vjust=0.5))
Leaf volatiles
Inventory
Dropped 2 samples from 2018 without leaf mass. Volatiles units are peak area per hour per gram dry leaf.
with(svl, table(PlantType, Treat, expt)) , , expt = otc
Treat
PlantType warm or in cool or out
agg 10 8
ten 9 9
, , expt = perc
Treat
PlantType warm or in cool or out
agg 4 4
ten 2 3Heatmap
#cairo_pdf("ipo_heatmap_leaf.pdf", width=12, height=12)
ph_leaf <- pheatmap(as.matrix(t(vol.leaf))^(1/3),
cluster_cols=T, show_colnames=F,
clustering_method="mcquitty", clustering_distance_rows="correlation",
clustering_distance_cols=vegdist(vol.leaf, method = "bray"),
clustering_callback = function(hc, ...){dendsort(hc, type="average")},
scale="none", color=inferno(512),
annotation_col = data.frame(svl %>% select("PlantType", "expt", "Treat"), row.names=rownames(vol.leaf)),
fontsize = 10, border_color = NA, legend=F, annotation_legend=T, cutree_rows=6
)
#dev.off()CAP - OTCs
svl.otc <- svl %>% filter(expt=="otc")
cap.leaf <- capscale(sqrt(vol.leaf[svl$expt=="otc",]) ~ PlantType+OI, distance="bray", data=svl.otc)
summ.cap <- summary(cap.leaf)
prop.expl <- summ.cap$cont$importance[2,c("CAP1","CAP2")]
ax.labs <- paste0("CAP", 1:2, " (",round(100*prop.expl),"% explained)")
(anova.cap <- anova(cap.leaf, by="term")) Permutation test for capscale under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(vol.leaf[svl$expt == "otc", ]) ~ PlantType + OI, data = svl.otc, distance = "bray")
Df SumOfSqs F Pr(>F)
PlantType 1 0.7636 5.0178 0.001 ***
OI 1 0.1273 0.8363 0.575
Residual 33 5.0222
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(cap.leaf, type="n", xlim=range(cap.leaf$CCA$wa[,1]), xlab=ax.labs[1], ylab=ax.labs[2])
title(paste("OTCs: I. agg and I. ten leaves 2019"))
legend("bottomleft", title=paste("Temperature","P =",anova.cap["OI","Pr(>F)"]), levels(svl.otc$OI), fill=c(3,2))
legend("bottomright", title=paste("Species","P =",anova.cap["PlantType","Pr(>F)"]), levels(svl.otc$PlantType), pch=c(19,15))
points(cap.leaf, display="sites", col=c(3,2)[svl.otc$OI], pch=c(19,15)[svl.otc$PlantType])
text(cap.leaf, display="cn")
text(cap.leaf, display="species", cex=0.7, col=ifelse(sqrt(cap.leaf$CCA$v[,1]^2 + cap.leaf$CCA$v[,2]^2)>0.2, alpha("purple",1),alpha("purple",0)))
CAP - Percivals
svl.perc <- svl %>% filter(expt=="perc")
cap.leaf.perc <- capscale(sqrt(vol.leaf[svl$expt=="perc",]) ~ PlantType+Temp, distance="bray", data=svl.perc)
summ.cap <- summary(cap.leaf.perc)
prop.expl <- summ.cap$cont$importance[2,c("CAP1","CAP2")]
ax.labs <- paste0("CAP", 1:2, " (",round(100*prop.expl),"% explained)")
(anova.cap <- anova(cap.leaf.perc, by="term")) Permutation test for capscale under reduced model
Terms added sequentially (first to last)
Permutation: free
Number of permutations: 999
Model: capscale(formula = sqrt(vol.leaf[svl$expt == "perc", ]) ~ PlantType + Temp, data = svl.perc, distance = "bray")
Df SumOfSqs F Pr(>F)
PlantType 1 0.09968 0.3418 0.925
Temp 1 0.20294 0.6958 0.575
Residual 10 2.91672