#R code for Jardine et al, Shedding light on sporopollenin chemistry, with reference to UV reconstructions
#Assumes all data tables have been saved as .txt files

untreated.ftir <- read.table("Table S1.txt", header = T,
                             row.names = 1)
acetolysis.ftir <- read.table("Table S2.txt", header = T,
                              row.names = 1)
untreated.ftir <- untreated.ftir[,20:1738]
acetolysis.ftir <- acetolysis.ftir[,20:1738] #trim off plunge at end

#Sample info files for plotting etc
sample.info <- read.table("Table S3.txt", header = T, row.names = 1)
untreated.ftir.sample.info <- sample.info[-11,] #No 6A_1 in this dataset

wavenumbers <- as.numeric(unlist(strsplit(colnames(untreated.ftir), "X"))[seq(2, 3438, by = 2)])
#alternative: wavenumbers <- as.numeric(t(as.data.frame(strsplit(colnames(untreated.ftir), "X")))[,2])


#Baseline correction with 1st order polynomial (i.e. straight line)
library(baseline)

#Untreated samples
untreated.ftir.base <- baseline(as.matrix(untreated.ftir),
                                method='modpolyfit', deg=1)
untreated.ftir.cor <- getCorrected(untreated.ftir.base)
colnames(untreated.ftir.cor) <- wavenumbers
untreated.ftir.cor <- as.data.frame(untreated.ftir.cor,
                                  row.names = rownames(untreated.ftir))

#Acetolysed samples
acetolysis.ftir.base <- baseline(as.matrix(acetolysis.ftir),
                                 method='modpolyfit', deg=1)
acetolysis.ftir.cor <- getCorrected(acetolysis.ftir.base)
colnames(acetolysis.ftir.cor) <- wavenumbers
acetolysis.ftir.cor <- as.data.frame(acetolysis.ftir.cor,
                                     row.names = rownames(acetolysis.ftir))


#Standardisation
#Standardisation function
#Each spectrum standardised to mean = 0, variance = 1 (z-scores)
z.stand <- function(x) {
  (x-mean(x))/sd(x)
}

untreated.ftir.z.stand <- t(apply(untreated.ftir.cor, 1, z.stand))
acetolysis.ftir.z.stand <- t(apply(acetolysis.ftir.cor, 1, z.stand))


#Plotting of mean spectrum for each dataset (Figure 1)
#Calculate means for plotting
untreated.ftir.z.stand.mean <- colMeans(untreated.ftir.z.stand)
acetolysis.ftir.z.stand.mean <- colMeans(acetolysis.ftir.z.stand)

#Untreated mean sample plot
par(mar = c(4, 2, 2, 2) + 0.1)
plot(wavenumbers, untreated.ftir.z.stand.mean, type = "l",
     xlim = c(3700, 700), ylim = c(-1.5, 2.5), las = 1, lwd = 1.5,
     yaxt = "n", xaxt = "n", xlab = "", ylab = "")
axis(1, lwd = 0, lwd.ticks = 1.5, tcl = 0.3,
     cex.axis = 0.8, mgp = c(1.5, 0.1, 0))
title(ylab = "Absorbance", line = 0.75)
title(xlab = "Wavenumber (cm-1)", line = 1.5)
title(main = "Z-score standardised, untreated", line = 0.75)

#Untreated mean sample, z-score standardised, 1800 to 800 wavenumbers
par(mar = c(4, 2, 2, 2) + 0.1)
plot(wavenumbers, untreated.ftir.z.stand.mean, type = "n",
     xlim = c(1800, 800), ylim = c(-1.5, 2.5), las = 1,
     yaxt = "n", xaxt = "n", xlab = "", ylab = "")
axis(1, lwd = 0, lwd.ticks = 1.5, tcl = 0.3,
     cex.axis = 0.8, mgp = c(1.5, 0.1, 0),
     xaxp = c(1800, 800, 10))
title(ylab = "Absorbance", line = 0.75)
title(xlab = "Wavenumber (cm-1)", line = 1.5)
title(main = "Z-score standardised, untreated", line = 0.75)
abline(v = 1515, col = "red", lty = 2)
abline(v = 1550, col = "blue", lty = 2)
abline(v = 1650, col = "blue", lty = 2)
lines(wavenumbers, untreated.ftir.z.stand.mean,
      lwd = 1.5)

#Acetolysed mean sample, z-score standardised
par(mar = c(4, 2, 2, 2) + 0.1)
plot(wavenumbers, acetolysis.ftir.z.stand.mean, type = "l",
     xlim = c(3700, 700), ylim = c(-2, 2.5), las = 1, lwd = 1.5,
     yaxt = "n", xaxt = "n", xlab = "", ylab = "")
lines(raman.wavenumbers, acetolysis.raman.z.stand.rep.mean+4,
      col = "grey50", lwd = 1.5)
axis(1, lwd = 0, lwd.ticks = 1.5, tcl = 0.3,
     cex.axis = 0.8, mgp = c(1.5, 0.1, 0))
title(ylab = "Absorbance", line = 0.75)
title(xlab = "Wavenumber (cm-1)", line = 1.5)
title(main = "Z-score standardised, acetolysed", line = 0.75)

#Acetolysed mean sample, z-score standardised, 1800 to 800 wavenumbers
par(mar = c(4, 2, 2, 2) + 0.1)
plot(wavenumbers, acetolysis.ftir.z.stand.mean, type = "n",
     xlim = c(1800, 800), ylim = c(-2, 2.5), las = 1,
     yaxt = "n", xaxt = "n", xlab = "", ylab = "")
axis(1, lwd = 0, lwd.ticks = 1.5, tcl = 0.3,
     cex.axis = 0.8, mgp = c(1.5, 0.1, 0),
     xaxp = c(1800, 800, 10))
title(ylab = "Absorbance", line = 0.75)
title(xlab = "Wavenumber (cm-1)", line = 1.5)
title(main = "Z-score standardised, acetolysed", line = 0.75)
abline(v = 1515, col = "red", lty = 2)
abline(v = 1550, col = "blue", lty = 2)
abline(v = 1650, col = "blue", lty = 2)
lines(wavenumbers, acetolysis.ftir.z.stand.mean,
      lwd = 1.5)


#Extraction and plotting of peak heights across the untreated gradient
#Function to extract peak heights
peaks <- function(x, wavenums, lower, upper) {
  peak.height <- max(x[wavenums>lower & wavenums<upper])
}

#Untreated, z-score standardised
untreated.ftir.z.stand.peaks <- data.frame(arom.1515 = apply(untreated.ftir.z.stand, 1, FUN = peaks,
                                                      wavenums = wavenumbers, lower = 1505,
                                                      upper = 1525),
                                        OH = apply(untreated.ftir.z.stand, 1, FUN = peaks,
                                                   wavenums = wavenumbers, lower = 3190,
                                                   upper = 3550)
)

#Untreated, not standardised
untreated.ftir.peaks <- data.frame(arom.1515 = apply(untreated.ftir.cor, 1, FUN = peaks,
                                                     wavenums = wavenumbers, lower = 1505,
                                                     upper = 1525),
                                   OH = apply(untreated.ftir.cor, 1, FUN = peaks,
                                              wavenums = wavenumbers, lower = 3190,
                                              upper = 3550)
)
untreated.ftir.peaks <- transform(untreated.ftir.peaks, Ph.OH = arom.1515/OH)


#Acetolysis, z-score standardised
acetolysis.ftir.z.stand.peaks <- data.frame(arom.1515 = apply(acetolysis.ftir.z.stand, 1, FUN = peaks,
                                                             wavenums = wavenumbers, lower = 1505,
                                                             upper = 1525),
                                           OH = apply(acetolysis.ftir.z.stand, 1, FUN = peaks,
                                                      wavenums = wavenumbers, lower = 3190,
                                                      upper = 3550)
)

#Acetolysis, not standardised
acetolysis.ftir.peaks <- data.frame(arom.1515 = apply(acetolysis.ftir.cor, 1, FUN = peaks,
                                                        wavenums = wavenumbers, lower = 1505,
                                                        upper = 1525),
                                              OH = apply(acetolysis.ftir.cor, 1, FUN = peaks,
                                                         wavenums = wavenumbers, lower = 3190,
                                                         upper = 3550)
)
acetolysis.ftir.peaks <- transform(acetolysis.ftir.peaks, Ph.OH = arom.1515/OH)


#Boxplots (Figure 2)
sample.names <- c("E10-7B", "E10-5A", "E10-5B", "E10-6A", "E10-6B", "E10-6C")
cols <- c("white", rep("grey80", 2), rep("grey50",3))

par(mfrow = c(3,2), mar = c(0,1,1.5,1), oma = c(3,3,0,0))
boxplot(untreated.ftir.z.stand.peaks$arom.1515 ~ untreated.ftir.sample.info$Plotting.order,
        names = sample.names, las = 1, ylab = "",
        main = "", col = cols, ylim = c(0.1, 0.6),
        yaxt = "n", xaxt = "n")
axis(1, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     labels = F, at = 1:6,
     mgp = c(1.5, -0.1, 0))
axis(2, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

boxplot(acetolysis.ftir.z.stand.peaks$arom.1515 ~ sample.info$Plotting.order,
        names = sample.names, las = 1, ylab = "",
        main = "", col = cols, ylim = c(-0.6, 0.6),
        yaxt = "n", xaxt = "n")
axis(1, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     labels = F, at = 1:6,
     mgp = c(1.5, -0.1, 0))
axis(2, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

boxplot(untreated.ftir.z.stand.peaks$OH ~ untreated.ftir.sample.info$Plotting.order,
        names = sample.names, las = 1, ylab = "",
        main = "", col = cols, ylim = c(0.5, 1.1),
        yaxt = "n", xaxt = "n")
axis(1, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     labels = F, at = 1:6,
     mgp = c(1.5, -0.1, 0))
axis(2, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

boxplot(acetolysis.ftir.z.stand.peaks$OH ~ sample.info$Plotting.order,
        names = sample.names, las = 1, ylab = "",
        main = "", col = cols, ylim = c(0.8, 2),
        yaxt = "n", xaxt = "n")
axis(1, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     labels = F, at = 1:6,
     mgp = c(1.5, -0.1, 0))
axis(2, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

boxplot(untreated.ftir.peaks$Ph.OH ~ untreated.ftir.sample.info$Plotting.order,
        names = sample.names, las = 1, ylab = "",
        main = "", col = cols, ylim = c(0.65, 0.85),
        yaxt = "n", xaxt = "n")
axis(1, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     labels = F, at = 1:6,
     mgp = c(1.5, -0.1, 0))
axis(2, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

boxplot(acetolysis.ftir.peaks$Ph.OH ~ sample.info$Plotting.order,
        names = sample.names, las = 1, ylab = "",
        main = "", col = cols, ylim = c(0.2, 0.9),
        yaxt = "n", xaxt = "n")
axis(1, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     labels = F, at = 1:6,
     mgp = c(1.5, -0.1, 0))
axis(2, lwd = 0, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)
mtext("Sample", side = 1, outer = T, line = 1.5, cex = 0.7)


#Bar charts of shading level medians (Figure 3)
par(mfrow = c(3,2), mar = c(1,1,1.5,1), oma = c(3,3,0,0))
barplot(aggregate(untreated.ftir.z.stand.peaks$arom.1515,
                  by = list(untreated.ftir.sample.info$Shade.level),
                  FUN = median)[c(2,3,1),2],
        col = cols[c(1,2,4)], ylim = c(-0.1, 0.5), yaxt = "n")
axis(2, lwd = 1, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)
abline(h = 0, lty = 2)

barplot(aggregate(acetolysis.ftir.z.stand.peaks$arom.1515,
                  by = list(sample.info$Shade.level),
                  FUN = median)[c(2,3,1),2],
        col = cols[c(1,2,4)], ylim = c(-0.1, 0.5), yaxt = "n")
axis(2, lwd = 1, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)
abline(h = 0, lty = 2)

barplot(aggregate(untreated.ftir.z.stand.peaks$OH,
                  by = list(untreated.ftir.sample.info$Shade.level),
                  FUN = median)[c(2,3,1),2],
        col = cols[c(1,2,4)], ylim = c(0, 2), yaxt = "n")
axis(2, lwd = 1, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

barplot(aggregate(acetolysis.ftir.z.stand.peaks$OH,
                  by = list(sample.info$Shade.level),
                  FUN = median)[c(2,3,1),2],
        col = cols[c(1,2,4)], ylim = c(0, 2), yaxt = "n")
axis(2, lwd = 1, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

barplot(aggregate(untreated.ftir.peaks$Ph.OH,
                  by = list(untreated.ftir.sample.info$Shade.level),
                  FUN = median)[c(2,3,1),2],
        col = cols[c(1,2,4)], ylim = c(0, 0.8), yaxt = "n")
axis(2, lwd = 1, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)

barplot(aggregate(acetolysis.ftir.peaks$Ph.OH,
                  by = list(sample.info$Shade.level),
                  FUN = median)[c(2,3,1),2],
        col = cols[c(1,2,4)], ylim = c(0, 0.8), yaxt = "n")
axis(2, lwd = 1, lwd.ticks = 1, tcl = 0.3,
     cex.axis = 0.5, mgp = c(1.5, 0.2, 0),
     las = 1)
mtext("Shading level", side = 1, outer = T, line = 1.5, cex = 0.7)


#Kruskal-Wallace tests
#Z-score standardised
kruskal.test(untreated.ftir.z.stand.peaks$arom.1515 ~
               untreated.ftir.sample.info$Shade.level)
#Kruskal-Wallis chi-squared = 6.6175, df = 2, p-value = 0.03656

kruskal.test(acetolysis.ftir.z.stand.peaks$arom.1515 ~
               sample.info$Shade.level)
#Kruskal-Wallis chi-squared = 9.4563, df = 2, p-value = 0.008843

kruskal.test(untreated.ftir.z.stand.peaks$OH ~
               untreated.ftir.sample.info$Shade.level)
#Kruskal-Wallis chi-squared = 1.4879, df = 2, p-value = 0.4752

kruskal.test(acetolysis.ftir.z.stand.peaks$OH ~
               sample.info$Shade.level)
#Kruskal-Wallis chi-squared = 2.5161, df = 2, p-value = 0.2842

#Aromatic/OH
kruskal.test(untreated.ftir.peaks$Ph.OH ~ untreated.ftir.sample.info$Shade.level)
#Kruskal-Wallis chi-squared = 10.228, df = 2, p-value = 0.006011

kruskal.test(acetolysis.ftir.peaks$Ph.OH ~ sample.info$Shade.level)
#Kruskal-Wallis chi-squared = 6.4422, df = 2, p-value = 0.03991