This notebook summarizes Magda’s analysis of data for the paper

Library load

library(ggplot2)
library(ggpubr)
library(reshape2)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
## 
## Attaching package: 'tidyr'
## The following object is masked from 'package:reshape2':
## 
##     smiths
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(readxl)
library(tinytex)
library(ggthemes)
library(wesanderson)
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats   1.0.0     ✔ readr     2.1.5
## ✔ lubridate 1.9.3     ✔ stringr   1.5.1
## ✔ purrr     1.0.2     ✔ tibble    3.2.1
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ plotly::filter() masks dplyr::filter(), stats::filter()
## ✖ dplyr::lag()     masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(ggsci)
library(ggsignif)
library(doBy)
## 
## Attaching package: 'doBy'
## 
## The following object is masked from 'package:dplyr':
## 
##     order_by
require(styler)
## Loading required package: styler
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(cowplot)
## 
## Attaching package: 'cowplot'
## 
## The following object is masked from 'package:lubridate':
## 
##     stamp
## 
## The following object is masked from 'package:ggthemes':
## 
##     theme_map
## 
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
library(ggforce)
library(see)
## 
## Attaching package: 'see'
## 
## The following objects are masked from 'package:ggsci':
## 
##     scale_color_material, scale_colour_material, scale_fill_material

Load RasPi imaging data

img_final_Ex1 <- read.csv("final_Ex1.csv")
img_final_Ex2 <- read.csv("final_Ex2.csv")
img_final_Ex3 <- read.csv("final_Ex3.csv")
img_final_Ex4 <- read.csv("final_Ex4.csv")
img_final_Ex5 <- read.csv("final_Ex5.csv")
img_final_Ex6 <- read.csv("final_Ex6.csv")
img_final_Ex7 <- read.csv("final_Ex7.csv")
img_final_Ex8 <- read.csv("final_Ex8.csv")

analyse <- rbind(img_final_Ex1, img_final_Ex2, img_final_Ex3, img_final_Ex4, img_final_Ex5, img_final_Ex6, img_final_Ex7, img_final_Ex8)
analyse$sid <- paste(analyse$Ex, analyse$sid, sep = "_") # Create a new column 'sid' by pasting together 'Ex' and 'sid' columns with '_' separator
row.names(analyse) <- NULL
if ("X" %in% colnames(analyse)) {
  analyse <- analyse[, -which(colnames(analyse) == "X")]
}
analyse
la_gen <- ggplot(data = analyse, aes(x = days, y = Area.sum)) +

  # Add a line geometry with transparency (alpha), size, and grouping by 'sid'
  geom_line(alpha = 0.3, size = 0.8, aes(group = sid)) +

  # Create facets (small multiples) based on the 'Treatment' variable with 3 columns
  facet_wrap(~Treatment, ncol = 3) +

  # Add a y-axis label
  ylab("Average Leaf area across each Accession") +

  # Apply a black-and-white theme
  theme_bw()
## 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.
# Plot the 'la_gen' ggplot object
plot(la_gen)

Smooth spline modeling of plant size throughout the experiment

# Create a vector 'days' containing values from 0 to 14
days <- 0:14

# Subset the 'analyse' data frame to select rows where 'sid' matches the first unique 'sid'
temp <- subset(analyse, analyse$sid == unique(analyse$sid)[1])

# Convert the 'days' column in 'temp' from a factor to numeric
temp$days <- as.numeric(as.character(temp$days))

# Display the 'days' column from the 'temp' data frame
temp$days
## [1]  0  4  2  6  8 10 12 14
temp
# Fit a cubic smoothing spline to the 'days' and 'Area.sum' data in the 'temp' data frame
plot.spl <- with(temp, smooth.spline(days, Area.sum, df = 4))
# Create a scatterplot of 'Area.sum' against 'days' using data from 'temp'
plot(Area.sum ~ days, data = temp)

# Add a blue line to the plot based on the cubic smoothing spline ('plot.spl')
lines(plot.spl, col = "blue")

# Add a red line to the plot by predicting values using the cubic smoothing spline ('plot.spl')
lines(predict(plot.spl, days), col = "red")

# Define a vector 'names' containing column names
names <- c(text = "ID", "Accession", "Treatment", "day", "Area.SUM")

# Create an empty data frame 'spline_data' to store the data
spline_data <- data.frame()

# Loop through each element 'k' in the 'names' vector
for (k in names) {
  # Create a new column in 'spline_data' with the name specified by 'k' and initialize it as character type
  spline_data[[k]] <- as.character()
}

spline_data
pred_temp <- predict(plot.spl, days)

# Display the predictions stored in 'pred_temp'
pred_temp
## $x
##  [1]  0  1  2  3  4  5  6  7  8  9 10 11 12 13 14
## 
## $y
##  [1]  127997.5  177043.5  225338.8  272391.9  318754.1  366282.7  421017.1
##  [8]  489098.1  572890.5  673335.1  789453.1  918413.3 1051893.0 1182245.2
## [15] 1310016.8
# Update the first 15 rows of the 'spline_data' data frame, setting the 4th column to 'pred_temp$x'
spline_data[1:15, 4] <- pred_temp$x

# Update the first 15 rows of the 'spline_data' data frame, setting the 5th column to 'pred_temp$y'
spline_data[1:15, 5] <- pred_temp$y

# Update the first 15 rows of the 'spline_data' data frame, setting the 1st column to 'temp$sid[1]'
spline_data[1:15, 1] <- temp$sid[1]

# Update the first 15 rows of the 'spline_data' data frame, setting the 2nd column to 'temp$Accession[1]'
spline_data[1:15, 2] <- temp$Accession[1]

# Update the first 15 rows of the 'spline_data' data frame, setting the 3rd column to 'temp$Treatment[1]'
spline_data[1:15, 3] <- temp$Treatment[1]

spline_data
# Create a copy of 'spline_data' and store it in 'final_spline'
final_spline <- spline_data

# Get unique 'sid' values from the 'analyse' data frame and calculate their count
all_plants <- unique(analyse$sid)
rownames(all_plants) <- NULL
length(all_plants)
## [1] 586
for (i in 2:586) {
  # Subset 'temp' data frame for the current 'sid'
  temp <- subset(analyse, analyse$sid == all_plants[i])

  # Check if 'temp' has more than 4 rows
  if (dim(temp)[1] > 4) {
    # Fit a cubic smoothing spline to 'days' and 'Area.sum' in 'temp'
    plot.spl <- with(temp, smooth.spline(days, Area.sum, df = 4))

    # Generate predictions using the smoothing spline
    pred_temp <- predict(plot.spl, days)

    # Update specific columns in 'spline_data' with 'pred_temp' and 'temp' values
    spline_data[1:15, 4] <- pred_temp$x
    spline_data[1:15, 5] <- pred_temp$y
    spline_data[1:15, 1] <- temp$sid[1]
    spline_data[1:15, 2] <- temp$Accession[1]
    spline_data[1:15, 3] <- temp$Treatment[1]

    # Append 'spline_data' to 'final_spline'
    final_spline <- rbind(final_spline, spline_data)
  } else {
    # Set specific columns in 'spline_data' when 'temp' has <= 4 rows
    spline_data[1:15, 4] <- days
    spline_data[1:15, 5] <- 0
    spline_data[1:15, 1] <- temp$sid[1]
    spline_data[1:15, 2] <- temp$Accession[1]
    spline_data[1:15, 3] <- temp$Treatment[1]
  }
}

final_spline
# Convert the 'day' column in 'final_spline' to numeric type
final_spline$day <- as.numeric(as.character(final_spline$day))

# Convert the 'Area.SUM' column in 'final_spline' to numeric type
final_spline$Area.SUM <- as.numeric(as.character(final_spline$Area.SUM))

# Filter 'final_spline' to keep rows where 'Area.SUM' is greater than 1
final_spline <- final_spline[final_spline$Area.SUM > 1, ]

# Calculate the count of unique 'ID' values in the filtered 'final_spline'
length(unique(final_spline$ID))
## [1] 582

plotting increase in digital biomass over the course of the experiment:

# Convert 'Area.SUM' column to numeric type
final_spline$Area.SUM <- as.numeric(as.character(final_spline$Area.SUM))

# Convert 'day' column to a factor with numeric levels
final_spline$day <- as.factor(as.numeric(as.character(final_spline$day)))

# Create a line graph (Area_lgraph) using ggplot2
Area_lgraph <- ggplot(data = final_spline, aes(x = day, y = Area.SUM, group = ID, color = Treatment)) +
  geom_line(alpha = 0.1) +
  stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = Treatment), alpha = 0.3) + 
  stat_summary(fun = mean, aes(group = Treatment), size = 0.7, geom = "line", linetype = "dashed") +
  stat_compare_means(aes(group = Treatment), label = "p.signif", method = "t.test", hide.ns = TRUE) +
  labs(x = "Days After Stress", y = "Shoot Size (7 x SV pixels)") +
  scale_color_d3("category10") +
  theme(legend.position = "bottom")

Area_lgraph

# Convert the 'day' column in 'final_spline' to numeric
final_spline$day <- as.numeric(as.character(final_spline$day))

# Subset 'final_spline' to select rows where 'ID' matches the first element of 'all_plants'
temp <- subset(final_spline, final_spline$ID == all_plants[1])
temp
# Create a scatterplot of 'Area.SUM' against 'day' in the 'temp' data frame
plot(temp$Area.SUM ~ temp$day) +
# Linear regression line is added to the scatterplot using the lm() function
abline(lm(temp$Area.SUM ~ temp$day))

## integer(0)
# Fit a linear regression model to predict 'Area.SUM' based on 'day' in the 'temp' data frame
model <- lm(temp$Area.SUM ~ temp$day)

# Extract the growth rate (GR) coefficient from the linear regression model
GR <- model$coefficients[2]

# Calculate the coefficient of determination (R-squared) for the model
R2 <- summary(model)$r.squared
# Define vector 'names' containing column names for the 'growth_data' data frame
names <- c(text = "ID", "Accession", "Treatment", "GR", "R2")

# Create an empty data frame 'growth_data' to store growth-related data
growth_data <- data.frame()

# Loop through each element 'k' in the 'names' vector
for (k in names) {
  # Create a new column in 'growth_data' with the name specified by 'k' and set as character type
  growth_data[[k]] <- as.character()
}

# Assign the 'ID' value from the 'temp' data frame to the first row, first column of 'growth_data'
growth_data[1, 1] <- temp$ID[1]

# Assign the 'Accession' value from the 'temp' data frame to the first row, second column of 'growth_data'
growth_data[1, 2] <- temp$Accession[1]

# Assign the 'Treatment' value from the 'temp' data frame to the first row, third column of 'growth_data'
growth_data[1, 3] <- temp$Treatment[1]

# Assign the growth rate (GR) to the first row, fourth column of 'growth_data'
growth_data[1, 4] <- GR

# Assign the coefficient of determination (R-squared, R2) to the first row, fifth column of 'growth_data'
growth_data[1, 5] <- R2

growth_data
# Create vector 'all_plants' containing unique 'ID' values from 'final_spline'
all_plants <- unique(final_spline$ID)

# Loop through each unique 'ID' in 'all_plants'
for (i in 1:length(all_plants)) {
  # Subset 'final_spline' to select data for the current 'ID'
  temp <- subset(final_spline, final_spline$ID == all_plants[i])

  # Fit a linear regression model to predict 'Area.SUM' based on 'day'
  model <- lm(temp$Area.SUM ~ temp$day)

  # Extract the growth rate (GR) coefficient from the linear regression model
  GR <- model$coefficients[2]

  # Calculate the coefficient of determination (R-squared, R2) for the model
  R2 <- summary(model)$r.squared

  # Populate 'growth_data' with relevant information for the current 'ID'
  growth_data[i, 1] <- temp$ID[1]
  growth_data[i, 2] <- temp$Accession[1]
  growth_data[i, 3] <- temp$Treatment[1]
  growth_data[i, 4] <- GR
  growth_data[i, 5] <- R2
}

growth_data
# Convert the 'GR' and 'R2' columns in 'growth_data' to numeric
growth_data$GR <- as.numeric(as.character(growth_data$GR))
growth_data$R2 <- as.numeric(as.character(growth_data$R2))

growth_data_clean <- subset(growth_data, growth_data$R2 > 0.7)
growth_data_clean <- subset(growth_data_clean, growth_data_clean$GR > 0)
growth_data_clean
decode <- read.csv("Supplemental-Table1-Germplasm-info - Table S1.csv")
decode <- decode[,c(1,3)]
growth_data_clean
colnames(decode)[1] <- "Accession"
decode$Accession <- gsub("T", "TDP", decode$Accession)
decode$Accession %in% growth_data_clean$Accession
##   [1] FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE
##  [13]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE
##  [25] FALSE FALSE FALSE  TRUE FALSE  TRUE FALSE  TRUE FALSE FALSE  TRUE  TRUE
##  [37]  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [49]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE
##  [61] FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
##  [73]  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [85]  TRUE  TRUE FALSE  TRUE  TRUE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
##  [97]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE
## [109] FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [121]  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## [133] FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [145] FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE FALSE FALSE
## [157]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE
## [169]  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE
## [181] FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE  TRUE
## [193] FALSE  TRUE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [205]  TRUE  TRUE  TRUE FALSE  TRUE FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## [217]  TRUE  TRUE  TRUE FALSE  TRUE  TRUE FALSE FALSE FALSE  TRUE  TRUE  TRUE
## [229]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [241]  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE
## [253]  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
## [277]  TRUE FALSE FALSE  TRUE  TRUE  TRUE FALSE  TRUE  TRUE  TRUE FALSE  TRUE
## [289] FALSE FALSE  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE
## [301] FALSE  TRUE  TRUE FALSE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE
## [313]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE
## [325]  TRUE  TRUE FALSE  TRUE  TRUE FALSE  TRUE FALSE  TRUE  TRUE FALSE  TRUE
## [337]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE  TRUE
## [349]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE
growth_data2 <- merge(decode, growth_data_clean, by="Accession", all = T)
growth_data3 <- na.omit(growth_data2)
growth_data3$Type <- gsub("wild", "Wild", growth_data3$Type)
want <- c("Wild", "Cultivated", "Weedy", "Breeding Line")
growth_data4 <- subset(growth_data3, growth_data3$Type %in% want)
growth_data4 <- subset(growth_data4, growth_data4$GR > 0)
growth_data4$Type <- factor(growth_data4$Type, levels = c("Breeding Line", "Cultivated", "Weedy", "Wild"))

growth_data4$Treatment <- factor(growth_data4$Treatment, levels = c("Control", "Drought"))

GR_overall <- ggerrorplot(growth_data4, x = "Treatment", y = "GR", color = "Treatment", fill = "Treatment", facet.by = "Type", ncol = 4,
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Growth Rate (pix / day)", add.params = list(color = "darkgray")) + scale_color_d3("category10") 
GR_overall <- GR_overall + stat_compare_means(method = "aov", label.y = 180000)
GR_overall <- GR_overall + rremove("legend")
GR_overall

growth_data_mean <- growth_data4 %>% 
  group_by(Treatment, Type) %>% 
  summarize(class_mean = mean(GR))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
growth_data_mean
32188.10/72436.56
## [1] 0.4443626
11016.08/29075.07
## [1] 0.378884
growth_data_clean2 <- growth_data4[, c(1:2, 4:5)]
# Compute a summary of 'GR' values grouped by 'Treatment' and 'Accession' using the 'summaryBy' function
growth_sum <- summaryBy(GR ~ Treatment + Accession + Type, data = growth_data_clean2)
growth_sum
# Reshape the 'growth_sum' data frame using the 'dcast' function to create a new data frame 'growth_cast'
# where 'Accession' becomes the rows, and 'Treatment' values become columns
growth_cast <- dcast(growth_sum, Accession + Type ~ Treatment)
## Using GR.mean as value column: use value.var to override.
growth_cast
# Convert 'Control' and 'Drought' columns to numeric type
growth_cast <- na.omit(growth_cast)
# Calculate the ratio of 'GR.Drought' to 'GR.Control' and store it in the 'GR.STI' column
growth_cast$STI <- growth_cast$Drought / growth_cast$Control
growth_cast
growth_cast <- subset(growth_cast, growth_cast$STI < 2 )
growth_cast$Type <- factor(growth_cast$Type, levels = c("Breeding Line", "Cultivated", "Weedy", "Wild"))


GR_STI <- ggerrorplot(growth_cast, x = "Type", y = "STI", 
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Relative Growth Rate (Fraction of Control)", add.params = list(color = "darkgray"))
GR_STI <- GR_STI + stat_compare_means(method = "aov")
GR_STI <- GR_STI + rremove("legend")
GR_STI

STI_mean <- growth_cast %>% 
  group_by(Type) %>% 
  summarize(class_mean = mean(STI), class_max = max(STI), class_min = min(STI))

STI_mean

FW data:

FW1 <- read_xlsx("FWDWWC_1.xlsx")
FW2 <- read_xlsx("FWDWWC_2.xlsx")
FW3 <- read_xlsx("FWDWWC_3.xlsx")
FW4 <- read_xlsx("FWDWWC_4.xlsx")
FW5 <- read_xlsx("FWDWWC_5.xlsx")
FW6 <- read_xlsx("FWDWWC_6.xlsx")
FW7 <- read_xlsx("FWDWWC_7.xlsx")
FW8 <- read_xlsx("FWDWWC_8.xlsx")

FW <- rbind(FW1, FW2, FW3, FW4, FW5, FW6, FW7, FW8)
FW
last_day <- subset(final_spline, final_spline$day == 14)
colnames(last_day)[2] <- "Accession"
colnames(last_day)[3] <- "Treatment"
colnames(last_day)[4] <- "Day"
FW$ID <- paste(FW$Ex, "_", FW$Accession, "_", substr(FW$Treatment, 1, 1), sep="")

unique(FW$Accession)%in% unique(last_day$Accession)
##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [166] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [196] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [226] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [256] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [271] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
growth_data4$Accession %in% last_day$Accession
##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [166] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [196] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [226] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [256] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [271] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [286] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [301] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [316] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [331] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [346] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [361] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [376] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [391] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [406] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [421] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [436] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
last_day$ID <- gsub("C_1", "C", last_day$ID)
last_day$ID <- gsub("D_1", "D", last_day$ID)

growth_and_LD <- merge(last_day, growth_data4, by = c("ID", "Accession", "Treatment"), all = T)
growth_and_LD
bye <- c("TDP_359_c", "TDP_22_c")

growth_and_LD2 <- subset(growth_and_LD, !(growth_and_LD$Accession %in% bye))
growth_and_LD2
FW_GR_LD <- merge(growth_and_LD2, FW, by = c("ID", "Accession", "Treatment"), all=T)
FW_GR_LD
FW_GR_LD2 <- subset(FW_GR_LD, !(FW_GR_LD$Accession %in% bye))
FW_GR_LD2
FW_GR_LD3 <- subset(FW_GR_LD2, FW_GR_LD2$Type %in% want)
FW_GR_LD3
unique(FW_GR_LD3$Type)
## [1] Cultivated    Wild          Weedy         Breeding Line
## Levels: Breeding Line Cultivated Weedy Wild
FW_GR_LD3$FW <- as.numeric(as.character(FW_GR_LD3$FW))
FW_GR_LD3$DW <- as.numeric(as.character(FW_GR_LD3$DW))
FW_GR_LD3$WC <- as.numeric(as.character(FW_GR_LD3$WC))
FW_GR_LD4 <- na.omit(FW_GR_LD3)
FW_overall <- ggerrorplot(FW_GR_LD4, x = "Treatment", y = "FW", color = "Treatment", fill = "Treatment", facet.by = "Type", ncol = 4,
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Fresh Weight (g)", add.params = list(color = "darkgray")) + scale_color_d3("category10") 
FW_overall <- FW_overall + stat_compare_means(method = "aov", label.y = 30)
FW_overall <- FW_overall + rremove("legend")
FW_overall

FW_data_mean <- FW_GR_LD4 %>% 
  group_by(Treatment, Type) %>% 
  summarize(class_mean = mean(FW), max = max(FW), min = min(FW))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
FW_data_mean
# Rel decrease BL
6.812500 / 15.100000
## [1] 0.4511589
# Rel decrease Wild
2.885937 / 6.610526
## [1] 0.4365669
DW_overall <- ggerrorplot(FW_GR_LD4, x = "Treatment", y = "DW", color = "Treatment", fill = "Treatment", facet.by = "Type", ncol = 4,
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Dry Weight (g)", add.params = list(color = "darkgray")) + scale_color_d3("category10") 
DW_overall <- DW_overall + stat_compare_means(method = "aov", label.y = 7.5)
DW_overall <- DW_overall + rremove("legend")
DW_overall

DW_data_mean <- FW_GR_LD4 %>% 
  group_by(Treatment, Type) %>% 
  summarize(class_mean = mean(DW), max = max(DW), min = min(DW))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
DW_data_mean
# Rel decrease for BL
1.6200000 / 3.4027273
## [1] 0.4760887
# Rel decrease for Wild
0.7142969 / 1.5312030
## [1] 0.4664939
WC_overall <- ggerrorplot(FW_GR_LD4, x = "Treatment", y = "WC", color = "Treatment", fill = "Treatment", facet.by = "Type", ncol = 4,
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Water Content (g)", add.params = list(color = "darkgray")) + scale_color_d3("category10") 
WC_overall <- WC_overall + stat_compare_means(method = "aov", label.y = 30)
WC_overall <- WC_overall + rremove("legend")
WC_overall

WC_data_mean <- FW_GR_LD4 %>% 
  group_by(Treatment, Type) %>% 
  summarize(class_mean = mean(WC), max = max(WC), min = min(WC))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
WC_data_mean
# BL
5.192500 / 11.697273
## [1] 0.4439069
# Cultivated
3.789704 / 9.037162
## [1] 0.4193467
# Wild
2.227745 / 4.815000
## [1] 0.4626677
FW_vs_Area <- ggscatter(FW_GR_LD4, x = "FW", y = "Area.SUM", xlab = "Shoot Fresh Weight (g)", ylab = "Projected shoot area (pixels)", color = "Treatment", rug = TRUE) +
  stat_cor() +
  scale_color_d3() +
  theme_minimal() # Set custom colors
FW_vs_Area

DW_vs_Area <- ggscatter(FW_GR_LD4, x = "DW", y = "Area.SUM", xlab = "Shoot Dry weight (g)", ylab = "Projected shoot area (pixels)", color = "Treatment", rug = TRUE) +
  stat_cor() +
  scale_color_d3() + theme_minimal()
DW_vs_Area

WC_vs_Area <- ggscatter(FW_GR_LD4, x = "Area.SUM", y = "WC", color = "Treatment", rug = TRUE) +
  stat_cor() +
  scale_color_d3()
WC_vs_Area

PhotoSynQ Data:

1st experiment:

PSQ_all <- list.files(pattern = "PQ")
PSQ1 <- read.csv(PSQ_all[1])
PSQ1 <- PSQ1 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ1 <- PSQ1[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
unique(PSQ1$time)
##   [1] "12/16/2022 15:22" "12/16/2022 15:21" "12/16/2022 15:17"
##   [4] "12/16/2022 15:16" "12/16/2022 15:15" "12/16/2022 15:14"
##   [7] "12/16/2022 15:12" "12/16/2022 15:11" "12/16/2022 15:10"
##  [10] "12/16/2022 15:09" "12/16/2022 15:08" "12/16/2022 15:07"
##  [13] "12/16/2022 15:06" "12/16/2022 15:05" "12/16/2022 15:03"
##  [16] "12/16/2022 15:02" "12/16/2022 15:01" "12/16/2022 15:00"
##  [19] "12/16/2022 14:59" "12/16/2022 14:57" "12/16/2022 14:56"
##  [22] "12/16/2022 14:55" "12/16/2022 14:54" "12/16/2022 14:53"
##  [25] "12/16/2022 14:52" "12/16/2022 14:51" "12/16/2022 14:50"
##  [28] "12/16/2022 14:49" "12/16/2022 14:48" "12/16/2022 14:47"
##  [31] "12/16/2022 14:46" "12/16/2022 14:45" "12/16/2022 14:44"
##  [34] "12/16/2022 14:43" "12/16/2022 14:42" "12/16/2022 14:41"
##  [37] "12/16/2022 14:40" "12/16/2022 14:39" "12/16/2022 14:38"
##  [40] "12/16/2022 14:37" "12/16/2022 14:36" "12/16/2022 14:34"
##  [43] "12/16/2022 14:33" "12/16/2022 14:32" "12/16/2022 14:31"
##  [46] "12/16/2022 14:30" "12/16/2022 14:29" "12/16/2022 14:27"
##  [49] "12/16/2022 14:26" "12/16/2022 14:25" "12/16/2022 14:24"
##  [52] "12/16/2022 14:23" "12/16/2022 14:22" "12/16/2022 14:21"
##  [55] "12/16/2022 14:20" "12/9/2022 12:26"  "12/9/2022 12:25" 
##  [58] "12/9/2022 12:22"  "12/9/2022 12:21"  "12/9/2022 12:20" 
##  [61] "12/9/2022 12:19"  "12/9/2022 12:18"  "12/9/2022 12:17" 
##  [64] "12/9/2022 12:16"  "12/9/2022 12:15"  "12/9/2022 12:14" 
##  [67] "12/9/2022 12:13"  "12/9/2022 12:12"  "12/9/2022 12:11" 
##  [70] "12/9/2022 12:10"  "12/9/2022 12:08"  "12/9/2022 12:07" 
##  [73] "12/9/2022 12:06"  "12/9/2022 12:05"  "12/9/2022 12:04" 
##  [76] "12/9/2022 12:03"  "12/9/2022 12:01"  "12/9/2022 12:00" 
##  [79] "12/9/2022 11:59"  "12/9/2022 11:58"  "12/9/2022 11:57" 
##  [82] "12/9/2022 11:53"  "12/9/2022 11:52"  "12/9/2022 11:51" 
##  [85] "12/9/2022 11:50"  "12/9/2022 11:49"  "12/9/2022 11:48" 
##  [88] "12/9/2022 11:47"  "12/9/2022 11:46"  "12/9/2022 11:44" 
##  [91] "12/9/2022 11:43"  "12/9/2022 11:42"  "12/9/2022 11:41" 
##  [94] "12/9/2022 11:40"  "12/9/2022 11:39"  "12/9/2022 11:35" 
##  [97] "12/9/2022 11:34"  "12/9/2022 11:33"  "12/9/2022 11:32" 
## [100] "12/9/2022 11:31"  "12/9/2022 11:30"  "12/9/2022 11:29" 
## [103] "12/9/2022 11:28"  "12/9/2022 11:27"  "12/9/2022 11:26" 
## [106] "12/9/2022 11:21"  "12/9/2022 11:20"  "12/9/2022 11:19" 
## [109] "12/9/2022 11:18"  "12/9/2022 11:17"
colnames(PSQ1) <- gsub("\\.", "_", colnames(PSQ1))
colnames(PSQ1)[1] <- ("Pot.No")
colnames(PSQ1)
##  [1] "Pot.No"                "Treatment"             "sid"                  
##  [4] "Accession"             "AmbientHumidity"       "FmPrime"              
##  [7] "FoPrime"               "Fs"                    "FvP"                  
## [10] "LeafTemperature"       "Leafthickness"         "LightIntensity"       
## [13] "NPQt"                  "Phi2"                  "PhiNO"                
## [16] "PhiNPQ"                "PS1ActiveCenters"      "PS1OpenCenters"       
## [19] "PS1OverReducedCenters" "PS1OxidizedCenters"    "SPAD"                 
## [22] "time"
PSQ1$time <- gsub(".*12/9/2022.*", "day6", PSQ1$time)
PSQ1$time <- gsub(".*12/16/2022.*", "day13", PSQ1$time)
PSQ1 <- subset(PSQ1, PSQ1$Leafthickness < 2)
mPSQ1 <- melt(PSQ1, id = c("Pot.No", "time"))
PSQ1_cast <- dcast(mPSQ1, Pot.No ~ time + variable)
PSQ1_cast$exp <- 1

2nd experiment:

PSQ2 <- read.csv(PSQ_all[2])
PSQ2 <- PSQ2 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ2 <- PSQ2[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
colnames(PSQ2) <- gsub("\\.", "_", colnames(PSQ2))
colnames(PSQ2)[1] <- "Pot.No"
PSQ2$time <- gsub(".*12/29/2022.*", "day6", PSQ2$time)
PSQ2$time <- gsub(".*1/5/2023.*", "day13", PSQ2$time)
PSQ2 <- subset(PSQ2, PSQ2$Leafthickness < 2)
mPSQ2 <- melt(PSQ2, id = c("Pot.No", "time"))
PSQ2_cast <- dcast(mPSQ2, Pot.No ~ time + variable)
PSQ2_cast$exp <- 2

3rd experiment:

PSQ3 <- read.csv(PSQ_all[3])
PSQ3 <- PSQ3 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ3 <- PSQ3[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
PSQ3$time <- gsub(".*1/17/2023.*", "day6", PSQ3$time)
PSQ3$time <- gsub(".*1/24/2023.*", "day13", PSQ3$time)
colnames(PSQ3) <- gsub("\\.", "_", colnames(PSQ3))
colnames(PSQ3)[1] <- "Pot.No"
mPSQ3 <- melt(PSQ3, id = c("Pot.No", "time"))
PSQ3_cast <- dcast(mPSQ3, Pot.No ~ time + variable)
PSQ3_cast$exp <- 3

4th experiment:

PSQ4 <- read.csv(PSQ_all[4])
PSQ4 <- PSQ4 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ4 <- PSQ4[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
colnames(PSQ4) <- gsub("\\.", "_", colnames(PSQ4))
colnames(PSQ4)[1] <- "Pot.No"
PSQ4$time <- gsub(".*4/21/2023.*", "day6", PSQ4$time)
PSQ4$time <- gsub(".*4/28/2023.*", "day13", PSQ4$time)
PSQ4 <- subset(PSQ4, PSQ4$Leafthickness < 2)
mPSQ4 <- melt(PSQ4, id = c("Pot.No", "time"))
PSQ4_cast <- dcast(mPSQ4, Pot.No ~ time + variable)
PSQ4_cast$exp <- 4

5th experiment:

PSQ5 <- read.csv(PSQ_all[5])
PSQ5 <- PSQ5 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ5 <- PSQ5[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
colnames(PSQ5) <- gsub("\\.", "_", colnames(PSQ5))
colnames(PSQ5)[1] <- "Pot.No"
PSQ5$time <- gsub(".*5/5/2023.*", "day6", PSQ5$time)
PSQ5$time <- gsub(".*5/12/2023.*", "day13", PSQ5$time)
PSQ5 <- subset(PSQ5, PSQ5$Leafthickness < 2)
mPSQ5 <- melt(PSQ5, id = c("Pot.No", "time"))
PSQ5_cast <- dcast(mPSQ5, Pot.No ~ time + variable)
PSQ5_cast$exp <- 5

6th experiment:

PSQ6 <- read.csv(PSQ_all[6])
PSQ6 <- PSQ6 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ6 <- PSQ6[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
colnames(PSQ6) <- gsub("\\.", "_", colnames(PSQ6))
colnames(PSQ6)[1] <- "Pot.No"
PSQ6$time <- gsub(".*6/2/2023.*", "day6", PSQ6$time)
PSQ6$time <- gsub(".*6/9/2023.*", "day13", PSQ6$time)
PSQ6 <- subset(PSQ6, PSQ6$Leafthickness < 2)
mPSQ6 <- melt(PSQ6, id = c("Pot.No", "time"))
PSQ6_cast <- dcast(mPSQ6, Pot.No ~ time + variable)
PSQ6_cast$exp <- 6
PSQ7 <- read.csv(PSQ_all[7])
PSQ7 <- PSQ7 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ7 <- PSQ7[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
colnames(PSQ7) <- gsub("\\.", "_", colnames(PSQ7))
colnames(PSQ7)[1] <- "Pot.No"
PSQ7$time <- gsub(".*6/23/2023.*", "day6", PSQ7$time)
PSQ7$time <- gsub(".*6/30/2023.*", "day13", PSQ7$time)
PSQ7 <- subset(PSQ7, PSQ7$Leafthickness < 2)
mPSQ7 <- melt(PSQ7, id = c("Pot.No", "time"))
PSQ7_cast <- dcast(mPSQ7, Pot.No ~ time + variable)
PSQ7_cast$exp <- 7
PSQ8 <- read.csv(PSQ_all[8])
PSQ8 <- PSQ8 %>%
  rename(
    AmbientHumidity = "Ambient.Humidity",
    FvP = "FvP_over_FmP",
    LeafTemperature = "Leaf.Temperature",
    Leafthickness = "leaf_thickness",
    LightIntensity = "Light.Intensity..PAR.",
    PS1ActiveCenters = "PS1.Active.Centers",
    PS1OpenCenters = "PS1.Open.Centers",
    PS1OverReducedCenters = "PS1.Over.Reduced.Centers",
    PS1OxidizedCenters = "PS1.Oxidized.Centers"
  )
PSQ8 <- PSQ8[, c(4, 5, 7, 8, 11, 19:22, 28, 32, 35, 36, 42:48, 52, 54)]
colnames(PSQ8) <- gsub("\\.", "_", colnames(PSQ8))
colnames(PSQ8)[1] <- "Pot.No"
PSQ8$time <- gsub(".*7/14/2023.*", "day6", PSQ8$time)
PSQ8$time <- gsub(".*7/21/2023.*", "day13", PSQ8$time)
PSQ8 <- subset(PSQ8, PSQ8$Leafthickness < 2)
mPSQ8 <- melt(PSQ8, id = c("Pot.No", "time"))
PSQ8_cast <- dcast(mPSQ8, Pot.No ~ time + variable)
PSQ8_cast$exp <- 8
PSQ <- rbind(PSQ1_cast, PSQ2_cast)
PSQ <- rbind(PSQ, PSQ3_cast)
PSQ <- rbind(PSQ, PSQ4_cast)
PSQ <- rbind(PSQ, PSQ5_cast)
PSQ <- rbind(PSQ, PSQ6_cast)
PSQ <- rbind(PSQ, PSQ7_cast)
PSQ <- rbind(PSQ, PSQ8_cast)
PSQ2 <- PSQ[,c(1,42,22:41,6:21)]
colnames(PSQ2)[3] <- "Treatment"
colnames(PSQ2)[4] <- "ID"
colnames(PSQ2)[5] <- "Accession"
PSQ2$ID <- paste(PSQ2$exp, PSQ2$Accession, PSQ2$Treatment, sep="_")
PSQ2$ID <- gsub("Control", "C", PSQ2$ID)
PSQ2$ID <- gsub("Drought", "D", PSQ2$ID)
PSQ2
mPSQ2 <- melt(PSQ2, id.vars = c("Pot.No", "exp", "Treatment", "ID", "Accession"))
unique(mPSQ2$variable)
##  [1] day6_AmbientHumidity        day6_FmPrime               
##  [3] day6_FoPrime                day6_Fs                    
##  [5] day6_FvP                    day6_LeafTemperature       
##  [7] day6_Leafthickness          day6_LightIntensity        
##  [9] day6_NPQt                   day6_Phi2                  
## [11] day6_PhiNO                  day6_PhiNPQ                
## [13] day6_PS1ActiveCenters       day6_PS1OpenCenters        
## [15] day6_PS1OverReducedCenters  day6_PS1OxidizedCenters    
## [17] day6_SPAD                   day13_FmPrime              
## [19] day13_FoPrime               day13_Fs                   
## [21] day13_FvP                   day13_LeafTemperature      
## [23] day13_Leafthickness         day13_LightIntensity       
## [25] day13_NPQt                  day13_Phi2                 
## [27] day13_PhiNO                 day13_PhiNPQ               
## [29] day13_PS1ActiveCenters      day13_PS1OpenCenters       
## [31] day13_PS1OverReducedCenters day13_PS1OxidizedCenters   
## [33] day13_SPAD                 
## 33 Levels: day6_AmbientHumidity day6_FmPrime day6_FoPrime day6_Fs ... day13_SPAD
TOI <- c("day6_FvP", "day6_LeafTemperature",  "day6_Leafthickness", "day6_NPQt", "day6_Phi2", "day6_PhiNO", "day6_PhiNPQ", 
         "day6_SPAD", "day13_FvP", "day13_LeafTemperature", "day13_Leafthickness", "day13_NPQt", "day13_Phi2", "day13_PhiNO", 
         "day13_PhiNPQ", "day13_SPAD" )
mPSQ3 <- subset(mPSQ2, mPSQ2$variable %in% TOI)

mPSQ3$day <- mPSQ3$variable
mPSQ3$day <- gsub("_FvP", "", mPSQ3$day)
mPSQ3$day <- gsub("_LeafTemperature", "", mPSQ3$day)
mPSQ3$day <- gsub("_Leafthickness", "", mPSQ3$day)
mPSQ3$day <- gsub("_NPQt", "", mPSQ3$day)
mPSQ3$day <- gsub("_PhiNO", "", mPSQ3$day)
mPSQ3$day <- gsub("_Phi2", "", mPSQ3$day)
mPSQ3$day <- gsub("_PhiNPQ", "", mPSQ3$day)
mPSQ3$day <- gsub("_SPAD", "", mPSQ3$day)

mPSQ3$variable <- gsub("day6_", "", mPSQ3$variable)
mPSQ3$variable <- gsub("day13_", "", mPSQ3$variable)
mPSQ3$value <- as.numeric(as.character(mPSQ3$value))
library(reshape2)
cPSQ <- dcast(mPSQ3, Pot.No + exp + Treatment + ID + Accession + day ~ variable)
PSQ <- na.omit(cPSQ)
PSQ <- subset(PSQ, PSQ$Accession != "#N/A")
PSQ$day <- factor(PSQ$day, levels = c("day6", "day13"))
PSQ_decoded <- merge(PSQ, decode, by = "Accession")

PSQ_decoded$Treatment <- factor(PSQ_decoded$Treatment, levels = c("Control", "Drought"))
PSQ_decoded$Type <- factor(PSQ_decoded$Type, levels = c("Breeding Line", "Cultivated", "Weedy", "Wild"))
PSQ_decoded
PSQ_decoded <- subset(PSQ_decoded, PSQ_decoded$Type %in% want)
PSQ_decoded <- subset(PSQ_decoded, PSQ_decoded$FvP >= 0)
PSQ_decoded
FvFm_plot <- ggerrorplot(PSQ_decoded,
  y = "FvP", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "Fv'/Fm'"
)
FvFm_plot <- FvFm_plot + rremove("legend")
FvFm_plot <- FvFm_plot + stat_compare_means(method = "aov", label.y = 0.8)
FvFm_plot <- FvFm_plot + scale_color_d3("category10")
FvFm_plot

FvFm_data_mean <- PSQ_decoded %>% 
  group_by(Treatment, Type, day) %>% 
  summarize(class_mean = mean(FvP))
## `summarise()` has grouped output by 'Treatment', 'Type'. You can override using
## the `.groups` argument.
FvFm_data_mean
Temp_plot <- ggerrorplot(PSQ_decoded, y = "LeafTemperature", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "Leaf Temperature (C)")
Temp_plot <- Temp_plot + rremove("legend")
Temp_plot <- Temp_plot + stat_compare_means(method = "aov", label.y = 31)
Temp_plot <- Temp_plot + scale_color_d3("category10")
Temp_plot

Temp_data_mean <- PSQ_decoded %>% 
  group_by(Treatment, Type, day) %>% 
  summarize(class_mean = mean(LeafTemperature))
## `summarise()` has grouped output by 'Treatment', 'Type'. You can override using
## the `.groups` argument.
Temp_data_mean
# BL d6
29.35545 / 28.92913
## [1] 1.014737
# BL d13
29.22150 / 28.54300
## [1] 1.023771
# Cult d6
28.89705 / 28.14674
## [1] 1.026657
# Cult d13
28.58835 / 27.71385
## [1] 1.031555
# Wild d6
28.35085 / 27.62402
## [1] 1.026312
# Wild d13
28.30453 / 27.23772
## [1] 1.039167
Thick_plot <- ggerrorplot(PSQ_decoded, y = "Leafthickness", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "Leaf Thickness (mm)")
Thick_plot <- Thick_plot + rremove("legend")
Thick_plot <- Thick_plot + stat_compare_means(method = "aov", label.y = 1.7)
Thick_plot <- Thick_plot + scale_color_d3("category10")
Thick_plot

NPQt_plot <- ggerrorplot(PSQ_decoded, y = "NPQt", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "Non-Photochemical Quenching (a.u.)")
NPQt_plot <- NPQt_plot + rremove("legend")
NPQt_plot <- NPQt_plot + stat_compare_means(method = "aov", label.y = 22)
NPQt_plot <- NPQt_plot + scale_color_d3("category10")
NPQt_plot

Phi2_plot <- ggerrorplot(PSQ_decoded, y = "Phi2", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "Phi2 (PSII efficiency, a.u.)")
Phi2_plot <- Phi2_plot + rremove("legend")
Phi2_plot <- Phi2_plot + stat_compare_means(method = "aov", label.y = 0.7)
Phi2_plot <- Phi2_plot + scale_color_d3("category10")
Phi2_plot

PhiNO_plot <- ggerrorplot(PSQ_decoded, y = "PhiNO", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "PhiNO (lost:incoming light, a.u.)")
PhiNO_plot <- PhiNO_plot + rremove("legend")
PhiNO_plot <- PhiNO_plot + stat_compare_means(method = "aov", label.y = 0.25)
PhiNO_plot <- PhiNO_plot + scale_color_d3("category10")
PhiNO_plot

PhiNPQ_plot <- ggerrorplot(PSQ_decoded, y = "PhiNPQ", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "PhiNPQ (NPQ:incoming light, a.u.)")
PhiNPQ_plot <- PhiNPQ_plot + rremove("legend")
PhiNPQ_plot <- PhiNPQ_plot + stat_compare_means(method = "aov", label.y = 0.75)
PhiNPQ_plot <- PhiNPQ_plot + scale_color_d3("category10")
PhiNPQ_plot

SPAD_plot <- ggerrorplot(PSQ_decoded, y = "SPAD", x = "Treatment", fill = "Treatment", color = "Treatment",
  desc_stat = "mean_sd", add = "jitter", facet.by = c("day", "Type"),
  add.params = list(color = "darkgrey"),
  xlab = "", ylab = "Relative Chrolophyll Content (SPAD)")
SPAD_plot <- SPAD_plot + rremove("legend")
SPAD_plot <- SPAD_plot + stat_compare_means(method = "aov", label.y = 75)
SPAD_plot <- SPAD_plot + scale_color_d3("category10")
SPAD_plot

Evapotranspiration data

EVT_all <- list.files(pattern = "EVT")
melted_data <- list()

# Loop over ET1 to ET8
for (i in 1:8) {
  ET <- read_xlsx(EVT_all[i])
  ET <- ET %>% select(-1)
  ET_melted <- melt(ET, id.vars = c("PotTag", "Accession", "Treatment"))
  colnames(ET_melted)[4] <- "day"
  colnames(ET_melted)[5] <- "ET"
  ET_melted$day <- gsub("Day", "", ET_melted$day)
  ET_melted <- ET_melted %>%
    mutate(exp = i)
  melted_data[[i]] <- ET_melted
}

merged_data <- bind_rows(melted_data)
unique_days <- merged_data %>%
  group_by(exp) %>%
  distinct(day, .keep_all = TRUE) %>%
  select(exp, day)

merged_data$ET <- as.numeric(as.character(merged_data$ET))
merged_data$ID <- paste(merged_data$exp, merged_data$Accession, substr(merged_data$Treatment, 1, 1), sep = "_")
colnames(merged_data)[5] <- "EVT"
ET_clean <- subset(merged_data, merged_data$EVT > 0)
ET_clean <- subset(ET_clean, ET_clean$EVT < 200)
ET_clean$day <- as.numeric(ET_clean$day)
ET_clean$day <- as.factor(as.numeric(as.character(ET_clean$day)))

ET_lgraph <- ggplot(data = ET_clean, aes(x = day, y = EVT, group = ID, color = Treatment)) +
  geom_line(alpha = 0.1) +
  stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = Treatment), alpha = 0.3) +
  stat_summary(fun = mean, aes(group = Treatment), size = 0.7, geom = "line", linetype = "dashed") +
  stat_compare_means(aes(group = Treatment), label = "p.signif", method = "t.test", hide.ns = TRUE) +
  ylab("Evapotranspiration (g H2O / day)") +
  xlab("Days After Stress") +
  scale_color_d3()

ET_lgraph

Smooth spline modeling of EVT:

days <- 1:14
temp <- subset(ET_clean, ET_clean$ID == unique(ET_clean$ID)[1])
temp$day <- as.numeric(as.character(temp$day))
# Fit a cubic smoothing spline ('smooth.spline') to the 'day' and 'EVT' data in 'temp'
plot.spl <- with(temp, smooth.spline(day, EVT, df = 5))
# Create a scatter plot of 'EVT' against 'day' using data from 'temp'
plot(EVT ~ day, data = temp)
# Add a blue line representing the cubic smoothing spline to the scatter plot
lines(plot.spl, col = "blue")
# Add a red line representing the predicted values from the cubic smoothing spline for the 'days' vector
lines(predict(plot.spl, days), col = "red")

# Create a vector 'names' specifying column names for a new data frame 'EVT_spline_data'
names <- c(text = "ID", "Accession", "Treatment", "day", "EVT")

# Create an empty data frame 'EVT_spline_data' with columns set as character vectors
EVT_spline_data <- data.frame()
for (k in names) EVT_spline_data[[k]] <- as.character()
pred_temp <- predict(plot.spl, days)
EVT_spline_data[1:14, 4] <- pred_temp$x
EVT_spline_data[1:14, 5] <- pred_temp$y
EVT_spline_data[1:14, 1] <- temp$ID[1]
EVT_spline_data[1:14, 2] <- temp$Accession[1]
EVT_spline_data[1:14, 3] <- temp$Treatment[1]
EVT_spline_data
#Create a 'EVT_spline_data' data frame as 'final_EVT_spline'
final_EVT_spline <- EVT_spline_data
# Extract unique 'ID' values from 'ET_clean' and calculate the length
all_plants <- unique(ET_clean$ID)
length(all_plants)
## [1] 609
# Loop through unique 'ID' values starting from the second one (index 2) to the 609th
for (i in 2:609) {
  # Create a subset 'temp' of 'ET_clean' for the current 'ID'
  temp <- subset(ET_clean, ET_clean$ID == all_plants[i])

  # Check if the number of rows in 'temp' is greater than 4
  if (dim(temp)[1] > 4) {
    # Fit a cubic smoothing spline to the 'day' and 'EVT' data in 'temp'
    plot.spl <- with(temp, smooth.spline(day, EVT, df = 5))

    # Generate predictions for 'EVT' using the smoothing spline and 'days'
    pred_temp <- predict(plot.spl, days)

    # Update the first 14 rows of 'EVT_spline_data' with predictions and other information
    EVT_spline_data[1:14, 4] <- pred_temp$x
    EVT_spline_data[1:14, 5] <- pred_temp$y
    EVT_spline_data[1:14, 1] <- temp$ID[1]
    EVT_spline_data[1:14, 2] <- temp$Accession[1]
    EVT_spline_data[1:14, 3] <- temp$Treatment[1]

    # Append the updated 'EVT_spline_data' to 'final_EVT_spline'
    final_EVT_spline <- rbind(final_EVT_spline, EVT_spline_data)
  } else {
    # For 'temp' with fewer than 5 rows, set specific values in 'spline_data'
    spline_data[1:14, 4] <- days
    spline_data[1:14, 5] <- 0
    spline_data[1:14, 1] <- temp$ID[1]
    spline_data[1:14, 2] <- temp$Accession[1]
    spline_data[1:14, 3] <- temp$Treatment[1]
  }
}

final_EVT_spline
# Create a subset 'clean_EVT_spline' from 'final_EVT_spline' where 'EVT' is greater than 0
clean_EVT_spline <- subset(final_EVT_spline, final_EVT_spline$EVT > 0)

clean_EVT_spline
# Convert the 'day' column in 'clean_EVT_spline' to a factor with numeric values
clean_EVT_spline$day <- as.factor(as.numeric(clean_EVT_spline$day))

# Convert the 'EVT' column in 'clean_EVT_spline' to numeric type
clean_EVT_spline$EVT <- as.numeric(as.character(clean_EVT_spline$EVT))

# Exclude plants from 'clean_EVT_spline' that were previously excluded from 'growth_data_clean'
clean_EVT_spline <- subset(clean_EVT_spline, clean_EVT_spline$ID %in% growth_data_clean$ID)

# Create a line graph 'ET_lgraph_spline' for 'EVT' data with grouping by 'ID' and coloring by 'Treatment'
ET_lgraph_spline <- ggplot(data = clean_EVT_spline, aes(x = day, y = EVT, group = ID, color = Treatment)) +
  geom_line(alpha = 0.1) +
  stat_summary(fun.data = mean_se, geom = "ribbon", linetype = 0, aes(group = Treatment), alpha = 0.2) +
  stat_summary(fun = mean, aes(group = Treatment), size = 0.7, geom = "line", linetype = "dashed") +
  stat_compare_means(aes(group = Treatment), label = "p.signif", method = "t.test", hide.ns = TRUE) +
  ylab("Evapotranspiration (g H2O / day)") +
  xlab("Days After Stress") +
  scale_color_d3()+
  guides(color=FALSE)+
  theme(
    legend.position = "none",
    panel.background = element_rect(fill = "transparent"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank(),
    panel.border = element_rect(color = "black", fill = NA),
    axis.line = element_line(color = "black"),
    axis.text.x = element_text(size = 14),
    axis.text.y = element_text(size = 14)
  )
## Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
## of ggplot2 3.3.4.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
ET_lgraph_spline

# Calculate the summary statistics for 'EVT' grouped by 'ID,' 'Treatment,' and 'Accession' in 'clean_EVT_spline'
sum_ET <- summaryBy(EVT ~ ID + Treatment + Accession, data = clean_EVT_spline, FUN = function(x) c(sum = sum(x)))
# Reset row names to avoid index-based row names
rownames(sum_ET) <- NULL
sum_ET
# Convert the 'day' column in 'clean_EVT_spline' to numeric type
clean_EVT_spline$day <- as.numeric(as.character(clean_EVT_spline$day))

# Create a subset 'sub_EVT' of 'clean_EVT_spline' for days between 4 and 11
sub_EVT <- subset(clean_EVT_spline, clean_EVT_spline$day > 3)
sub_EVT <- subset(sub_EVT, sub_EVT$day < 12)

# Calculate the median 'EVT' values grouped by 'ID,' 'Treatment,' and 'Accession' in 'sub_EVT'
median_ET <- summaryBy(EVT ~ ID + Treatment + Accession, data = sub_EVT, FUN = function(x) c(med = median(x)))

# Reset row names to avoid index-based row names
rownames(median_ET) <- NULL

median_ET
# Merge the summary statistics data frames 'sum_ET' and 'median_ET' based on common columns: "ID," "Treatment," and "Accession"
EVT_final <- merge(sum_ET, median_ET, by = c("ID", "Treatment", "Accession"))

EVT_final
EVT_data2 <- merge(decode, EVT_final, by="Accession", all = T)
EVT_data3 <- na.omit(EVT_data2)
EVT_data3$Type <- gsub("wild", "Wild", EVT_data3$Type)
want <- c("Wild", "Cultivated", "Weedy", "Breeding Line")
EVT_data4 <- subset(EVT_data3, EVT_data3$Type %in% want)


EVT_data4$Treatment <- factor(EVT_data4$Treatment, levels = c("Control", "Drought"))

EVT_med_overall <- ggerrorplot(EVT_data4, x = "Treatment", y = "EVT.med", color = "Treatment", fill = "Treatment", facet.by = "Type", ncol = 4,
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Median evapotranspiration (g H2O / day)", add.params = list(color = "darkgray")) + scale_color_d3("category10") 
EVT_med_overall <- EVT_med_overall + stat_compare_means(method = "aov", label.y = 120)
EVT_med_overall <- EVT_med_overall + rremove("legend")
EVT_med_overall

EVT.med_data_mean <- EVT_data4 %>% 
  group_by(Treatment, Type) %>% 
  summarize(class_mean = mean(EVT.med), max = max(EVT.med), min = min(EVT.med))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
EVT.med_data_mean
# Rel decrease for Cultivated
27.02711 / 63.60887
## [1] 0.4248953
# Breedling Line:
21.99963 / 57.81133
## [1] 0.3805418
# Rel decrease for Wild
25.06244 / 50.43050
## [1] 0.4969699
EVT_sum_overall <- ggerrorplot(EVT_data4, x = "Treatment", y = "EVT.sum", color = "Treatment", fill = "Treatment", facet.by = "Type", ncol = 4,
                          desc_stat = "mean_sd", add = "jitter", 
                          xlab="", ylab= "Cumulative evapotranspiration (g H2O / experiment)", add.params = list(color = "darkgray")) + scale_color_d3("category10") 
EVT_sum_overall <- EVT_sum_overall + stat_compare_means(method = "aov", label.y = 1800)
EVT_sum_overall <- EVT_sum_overall + rremove("legend")
EVT_sum_overall

EVT.sum_data_mean <- EVT_data4 %>% 
  group_by(Treatment, Type) %>% 
  summarize(class_mean = mean(EVT.sum), max = max(EVT.sum), min = min(EVT.sum))
## `summarise()` has grouped output by 'Treatment'. You can override using the
## `.groups` argument.
EVT.sum_data_mean
# Rel decrease for Cultivated
375.9083 / 868.3132
## [1] 0.4329179
# Breedling Line:
312.8012 / 763.4955
## [1] 0.4096962
# Rel decrease for Wild
354.1226 / 699.6058
## [1] 0.5061745

Data Merge

Let’s bind this data with other data that we have

unique(PSQ_decoded$ID)
##   [1] "2_TDP_101_D" "2_TDP_101_C" "2_TDP_102_C" "2_TDP_102_D" "2_TDP_103_C"
##   [6] "2_TDP_103_D" "2_TDP_104_D" "2_TDP_104_C" "2_TDP_105_C" "2_TDP_105_D"
##  [11] "2_TDP_108_C" "2_TDP_108_D" "2_TDP_109_D" "2_TDP_109_C" "2_TDP_110_D"
##  [16] "2_TDP_110_C" "3_TDP_112_C" "7_TDP_113_C" "7_TDP_113_D" "7_TDP_114_D"
##  [21] "7_TDP_114_C" "7_TDP_115_C" "7_TDP_115_D" "7_TDP_120_C" "7_TDP_127_C"
##  [26] "7_TDP_127_D" "7_TDP_129_D" "7_TDP_129_C" "3_TDP_130_D" "3_TDP_130_C"
##  [31] "8_TDP_132_C" "8_TDP_132_D" "3_TDP_133_C" "3_TDP_134_D" "3_TDP_136_D"
##  [36] "3_TDP_136_C" "3_TDP_137_D" "3_TDP_137_C" "3_TDP_138_C" "3_TDP_138_D"
##  [41] "3_TDP_143_C" "3_TDP_143_D" "3_TDP_144_C" "4_TDP_147_C" "4_TDP_147_D"
##  [46] "3_TDP_149_D" "3_TDP_150_C" "3_TDP_150_D" "7_TDP_151_D" "7_TDP_151_C"
##  [51] "4_TDP_153_C" "4_TDP_153_D" "3_TDP_154_C" "3_TDP_154_D" "3_TDP_156_C"
##  [56] "3_TDP_156_D" "3_TDP_158_D" "3_TDP_158_C" "3_TDP_159_C" "3_TDP_159_D"
##  [61] "3_TDP_161_D" "3_TDP_161_C" "3_TDP_162_C" "3_TDP_162_D" "4_TDP_164_C"
##  [66] "4_TDP_164_D" "3_TDP_165_C" "3_TDP_165_D" "3_TDP_169_D" "3_TDP_169_C"
##  [71] "3_TDP_170_C" "3_TDP_170_D" "4_TDP_171_C" "4_TDP_171_D" "4_TDP_172_D"
##  [76] "4_TDP_172_C" "4_TDP_173_D" "4_TDP_173_C" "4_TDP_174_D" "4_TDP_174_C"
##  [81] "4_TDP_176_D" "4_TDP_176_C" "4_TDP_177_D" "4_TDP_177_C" "4_TDP_180_C"
##  [86] "4_TDP_180_D" "4_TDP_181_C" "4_TDP_181_D" "4_TDP_183_C" "4_TDP_183_D"
##  [91] "4_TDP_184_C" "4_TDP_184_D" "4_TDP_185_D" "4_TDP_185_C" "4_TDP_186_D"
##  [96] "4_TDP_186_C" "4_TDP_187_D" "4_TDP_187_C" "4_TDP_188_C" "4_TDP_188_D"
## [101] "4_TDP_189_C" "4_TDP_189_D" "4_TDP_190_D" "4_TDP_190_C" "4_TDP_191_C"
## [106] "4_TDP_191_D" "4_TDP_194_C" "4_TDP_194_D" "4_TDP_195_C" "4_TDP_195_D"
## [111] "4_TDP_197_D" "4_TDP_197_C" "4_TDP_198_D" "4_TDP_198_C" "5_TDP_199_C"
## [116] "5_TDP_199_D" "1_TDP_2_C"   "1_TDP_2_D"   "5_TDP_200_C" "5_TDP_200_D"
## [121] "5_TDP_201_D" "5_TDP_201_C" "5_TDP_202_D" "5_TDP_202_C" "4_TDP_203_D"
## [126] "4_TDP_203_C" "5_TDP_204_D" "5_TDP_204_C" "8_TDP_205_D" "8_TDP_205_C"
## [131] "4_TDP_206_D" "4_TDP_206_C" "4_TDP_207_D" "4_TDP_207_C" "5_TDP_209_C"
## [136] "5_TDP_209_D" "5_TDP_211_C" "5_TDP_211_D" "5_TDP_214_D" "5_TDP_214_C"
## [141] "5_TDP_215_D" "5_TDP_215_C" "4_TDP_216_C" "4_TDP_216_D" "5_TDP_218_C"
## [146] "5_TDP_218_D" "5_TDP_219_C" "5_TDP_219_D" "1_TDP_22_D"  "1_TDP_22_C" 
## [151] "5_TDP_220_D" "5_TDP_220_C" "5_TDP_221_D" "5_TDP_221_C" "5_TDP_222_C"
## [156] "5_TDP_222_D" "5_TDP_223_C" "5_TDP_223_D" "4_TDP_224_D" "5_TDP_226_D"
## [161] "5_TDP_226_C" "5_TDP_227_C" "5_TDP_227_D" "5_TDP_228_D" "5_TDP_228_C"
## [166] "5_TDP_229_D" "5_TDP_229_C" "5_TDP_230_C" "5_TDP_230_D" "5_TDP_231_D"
## [171] "5_TDP_231_C" "7_TDP_233_D" "5_TDP_234_D" "5_TDP_234_C" "5_TDP_236_C"
## [176] "5_TDP_236_D" "5_TDP_237_C" "5_TDP_237_D" "6_TDP_238_C" "6_TDP_238_D"
## [181] "6_TDP_239_D" "6_TDP_239_C" "6_TDP_240_C" "6_TDP_240_D" "5_TDP_241_C"
## [186] "5_TDP_241_D" "6_TDP_242_C" "6_TDP_242_D" "6_TDP_243_D" "6_TDP_243_C"
## [191] "6_TDP_244_C" "6_TDP_244_D" "6_TDP_245_D" "6_TDP_245_C" "6_TDP_248_D"
## [196] "6_TDP_248_C" "1_TDP_25_D"  "1_TDP_25_C"  "5_TDP_252_C" "5_TDP_252_D"
## [201] "6_TDP_258_C" "6_TDP_258_D" "6_TDP_259_D" "6_TDP_259_C" "6_TDP_262_C"
## [206] "6_TDP_262_D" "6_TDP_263_D" "6_TDP_263_C" "6_TDP_264_D" "6_TDP_264_C"
## [211] "6_TDP_265_C" "6_TDP_265_D" "6_TDP_266_D" "6_TDP_266_C" "6_TDP_267_C"
## [216] "6_TDP_267_D" "6_TDP_269_C" "6_TDP_269_D" "1_TDP_27_D"  "1_TDP_27_C" 
## [221] "5_TDP_271_C" "5_TDP_271_D" "6_TDP_272_C" "6_TDP_272_D" "6_TDP_273_C"
## [226] "6_TDP_273_D" "6_TDP_274_D" "6_TDP_274_C" "6_TDP_275_D" "6_TDP_275_C"
## [231] "6_TDP_276_C" "6_TDP_276_D" "6_TDP_279_D" "6_TDP_279_C" "6_TDP_284_D"
## [236] "6_TDP_284_C" "6_TDP_285_C" "6_TDP_285_D" "6_TDP_288_C" "6_TDP_288_D"
## [241] "6_TDP_293_C" "6_TDP_293_D" "6_TDP_294_C" "6_TDP_294_D" "6_TDP_295_C"
## [246] "6_TDP_295_D" "6_TDP_296_C" "6_TDP_296_D" "6_TDP_299_D" "6_TDP_299_C"
## [251] "7_TDP_302_C" "7_TDP_302_D" "7_TDP_303_C" "7_TDP_303_D" "7_TDP_307_C"
## [256] "7_TDP_307_D" "7_TDP_308_C" "7_TDP_308_D" "7_TDP_311_C" "7_TDP_315_C"
## [261] "1_TDP_33_D"  "1_TDP_33_C"  "1_TDP_34_C"  "1_TDP_34_D"  "8_TDP_354_C"
## [266] "8_TDP_354_D" "8_TDP_357_C" "8_TDP_357_D" "8_TDP_358_C" "8_TDP_358_D"
## [271] "1_TDP_36_D"  "1_TDP_36_C"  "8_TDP_360_C" "8_TDP_360_D" "8_TDP_361_D"
## [276] "8_TDP_361_C" "8_TDP_362_C" "8_TDP_362_D" "8_TDP_363_C" "8_TDP_363_D"
## [281] "8_TDP_364_C" "8_TDP_364_D" "7_TDP_365_D" "7_TDP_365_C" "8_TDP_367_D"
## [286] "8_TDP_367_C" "8_TDP_368_C" "8_TDP_368_D" "8_TDP_369_C" "8_TDP_369_D"
## [291] "1_TDP_37_D"  "1_TDP_37_C"  "8_TDP_370_C" "8_TDP_370_D" "8_TDP_372_C"
## [296] "8_TDP_372_D" "8_TDP_373_C" "8_TDP_373_D" "8_TDP_374_C" "8_TDP_374_D"
## [301] "8_TDP_375_D" "8_TDP_375_C" "8_TDP_376_C" "8_TDP_376_D" "8_TDP_377_C"
## [306] "8_TDP_377_D" "8_TDP_379_C" "8_TDP_379_D" "1_TDP_38_D"  "1_TDP_38_C" 
## [311] "8_TDP_380_C" "8_TDP_380_D" "8_TDP_381_D" "8_TDP_381_C" "8_TDP_382_C"
## [316] "8_TDP_382_D" "8_TDP_383_C" "8_TDP_383_D" "8_TDP_384_D" "8_TDP_384_C"
## [321] "8_TDP_385_C" "8_TDP_385_D" "8_TDP_386_C" "8_TDP_386_D" "8_TDP_388_D"
## [326] "8_TDP_388_C" "1_TDP_39_D"  "1_TDP_39_C"  "8_TDP_392_C" "8_TDP_392_D"
## [331] "8_TDP_394_D" "8_TDP_394_C" "8_TDP_395_C" "8_TDP_395_D" "8_TDP_396_C"
## [336] "8_TDP_396_D" "8_TDP_398_D" "8_TDP_398_C" "8_TDP_400_C" "8_TDP_400_D"
## [341] "5_TDP_402_D" "5_TDP_402_C" "5_TDP_403_C" "5_TDP_403_D" "5_TDP_404_D"
## [346] "5_TDP_404_C" "8_TDP_406_D" "8_TDP_406_C" "1_TDP_41_D"  "1_TDP_41_C" 
## [351] "8_TDP_416_C" "8_TDP_419_D" "8_TDP_419_C" "1_TDP_42_C"  "1_TDP_42_D" 
## [356] "8_TDP_423_D" "8_TDP_423_C" "1_TDP_43_C"  "1_TDP_43_D"  "8_TDP_431_D"
## [361] "8_TDP_431_C" "8_TDP_432_D" "8_TDP_432_C" "8_TDP_433_C" "8_TDP_433_D"
## [366] "8_TDP_434_D" "8_TDP_434_C" "8_TDP_435_D" "8_TDP_435_C" "8_TDP_436_C"
## [371] "8_TDP_436_D" "8_TDP_437_C" "8_TDP_437_D" "8_TDP_439_D" "8_TDP_439_C"
## [376] "1_TDP_44_C"  "1_TDP_44_D"  "8_TDP_440_C" "8_TDP_440_D" "8_TDP_441_C"
## [381] "8_TDP_441_D" "8_TDP_442_D" "8_TDP_442_C" "8_TDP_443_C" "8_TDP_443_D"
## [386] "8_TDP_445_C" "8_TDP_445_D" "8_TDP_446_C" "8_TDP_446_D" "8_TDP_447_C"
## [391] "8_TDP_447_D" "8_TDP_448_D" "8_TDP_448_C" "8_TDP_450_C" "8_TDP_450_D"
## [396] "8_TDP_451_D" "8_TDP_451_C" "8_TDP_452_D" "8_TDP_452_C" "8_TDP_453_C"
## [401] "8_TDP_453_D" "1_TDP_49_C"  "1_TDP_49_D"  "1_TDP_54_C"  "1_TDP_54_D" 
## [406] "1_TDP_55_C"  "1_TDP_55_D"  "1_TDP_56_D"  "1_TDP_56_C"  "1_TDP_57_D" 
## [411] "1_TDP_57_C"  "1_TDP_58_C"  "1_TDP_58_D"  "1_TDP_60_C"  "1_TDP_60_D" 
## [416] "1_TDP_61_C"  "1_TDP_61_D"  "1_TDP_62_D"  "1_TDP_62_C"  "1_TDP_63_C" 
## [421] "1_TDP_63_D"  "4_TDP_64_D"  "4_TDP_64_C"  "1_TDP_65_D"  "1_TDP_65_C" 
## [426] "1_TDP_68_C"  "1_TDP_68_D"  "1_TDP_69_D"  "1_TDP_69_C"  "1_TDP_70_D" 
## [431] "1_TDP_70_C"  "1_TDP_71_C"  "1_TDP_71_D"  "1_TDP_74_C"  "1_TDP_74_D" 
## [436] "1_TDP_78_D"  "1_TDP_78_C"  "1_TDP_79_C"  "1_TDP_79_D"  "1_TDP_80_D" 
## [441] "1_TDP_80_C"  "4_TDP_81_D"  "4_TDP_81_C"  "4_TDP_82_C"  "4_TDP_82_D" 
## [446] "2_TDP_83_C"  "2_TDP_83_D"  "2_TDP_84_D"  "2_TDP_84_C"  "4_TDP_85_C" 
## [451] "4_TDP_85_D"  "4_TDP_86_D"  "4_TDP_86_C"  "4_TDP_87_C"  "4_TDP_87_D" 
## [456] "2_TDP_88_D"  "2_TDP_88_C"  "7_TDP_89_D"  "7_TDP_91_C"  "2_TDP_92_C" 
## [461] "2_TDP_92_D"  "2_TDP_94_D"  "2_TDP_94_C"  "2_TDP_95_C"  "2_TDP_95_D" 
## [466] "4_TDP_96_C"  "4_TDP_96_D"  "2_TDP_97_D"  "2_TDP_97_C"  "2_TDP_98_D" 
## [471] "2_TDP_98_C"  "2_TDP_99_D"  "2_TDP_99_C"
colnames(FW_GR_LD4) %in% colnames(PSQ_decoded)
##  [1]  TRUE  TRUE  TRUE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## [13] FALSE FALSE
PSQ_decoded
mPSQ_decoded <- melt(PSQ_decoded, id.vars = c("Accession", "Pot.No", "exp", "Treatment", "ID", "day", "Type"))
mPSQ_decoded
cPSQ <- dcast(mPSQ_decoded, Pot.No + exp + Treatment + ID + Accession + Type ~ day + variable)
cPSQ
FW_PSQ <- merge(FW_GR_LD4, cPSQ, by=c("Accession", "ID", "Treatment", "Type"), all = T)
FW_PSQ2 <- na.omit(FW_PSQ) 
FW_PSQ2 <- FW_PSQ2[,c(1:7,11:13,17:32)]
FW_PSQ2
FW_PSQ_EVT <- merge(FW_PSQ2, EVT_data4, by=c("Accession", "ID", "Treatment", "Type"), all = T)
FW_PSQ_EVT

Correlation plots:

require(ggcorrplot)
## Loading required package: ggcorrplot
Control <- subset(FW_PSQ_EVT, FW_PSQ_EVT$Treatment == "Control")
dim(Control)
## [1] 198  28
Control <- Control[,c(6:28)]
Control <- na.omit(Control)
corr2 <- cor(Control)
head(corr2[, 1:6])
##           Area.SUM        GR        FW        DW        WC  day6_FvP
## Area.SUM 1.0000000 0.9757442 0.9559241 0.9020268 0.9593795 0.2035648
## GR       0.9757442 1.0000000 0.9134261 0.8517202 0.9195874 0.2098869
## FW       0.9559241 0.9134261 1.0000000 0.9659177 0.9973662 0.1826586
## DW       0.9020268 0.8517202 0.9659177 1.0000000 0.9445990 0.1342429
## WC       0.9593795 0.9195874 0.9973662 0.9445990 1.0000000 0.1939995
## day6_FvP 0.2035648 0.2098869 0.1826586 0.1342429 0.1939995 1.0000000
p.mat2 <- cor_pmat(Control)
head(p.mat2[, 1:5])
##               Area.SUM            GR            FW           DW            WC
## Area.SUM  0.000000e+00 4.066315e-110  3.377315e-89 1.010642e-61  4.813581e-92
## GR       4.066315e-110  0.000000e+00  6.411255e-66 6.799566e-48  1.947197e-68
## FW        3.377315e-89  6.411255e-66  0.000000e+00 3.535503e-98 8.336118e-189
## DW        1.010642e-61  6.799566e-48  3.535503e-98 0.000000e+00  2.956683e-81
## WC        4.813581e-92  1.947197e-68 8.336118e-189 2.956683e-81  0.000000e+00
## day6_FvP  8.525128e-03  6.646237e-03  1.849979e-02 8.464589e-02  1.226380e-02
p.mat2 <- cor_pmat(Control)
head(p.mat2[, 1:5])
##               Area.SUM            GR            FW           DW            WC
## Area.SUM  0.000000e+00 4.066315e-110  3.377315e-89 1.010642e-61  4.813581e-92
## GR       4.066315e-110  0.000000e+00  6.411255e-66 6.799566e-48  1.947197e-68
## FW        3.377315e-89  6.411255e-66  0.000000e+00 3.535503e-98 8.336118e-189
## DW        1.010642e-61  6.799566e-48  3.535503e-98 0.000000e+00  2.956683e-81
## WC        4.813581e-92  1.947197e-68 8.336118e-189 2.956683e-81  0.000000e+00
## day6_FvP  8.525128e-03  6.646237e-03  1.849979e-02 8.464589e-02  1.226380e-02
Correlation_Control <- ggcorrplot(
  corr2, # Correlation matrix to plot
  p.mat = p.mat2, # Optional p-values matrix for significance indicators
  type = "upper", # Type of plot (upper triangle)
  outline.col = "white", # Color of the outline around the plot
  colors = c("#6D9EC1", "white", "#E46726"), # Color scale for correlations
)

Correlation_Control

Drought <- subset(FW_PSQ_EVT, FW_PSQ_EVT$Treatment == "Drought")
dim(Drought)
## [1] 183  28
Drought <- Drought[,c(6:28)]
Drought <- na.omit(Drought)
corr2 <- cor(Drought)
head(corr2[, 1:6])
##           Area.SUM        GR        FW        DW        WC  day6_FvP
## Area.SUM 1.0000000 0.9005607 0.9438190 0.9041407 0.9479051 0.2675481
## GR       0.9005607 1.0000000 0.8275179 0.7800056 0.8352448 0.2774171
## FW       0.9438190 0.8275179 1.0000000 0.9782628 0.9977163 0.2523177
## DW       0.9041407 0.7800056 0.9782628 1.0000000 0.9620220 0.2303787
## WC       0.9479051 0.8352448 0.9977163 0.9620220 1.0000000 0.2571010
## day6_FvP 0.2675481 0.2774171 0.2523177 0.2303787 0.2571010 1.0000000
p.mat2 <- cor_pmat(Drought)
head(p.mat2[, 1:5])
##              Area.SUM           GR            FW            DW            WC
## Area.SUM 0.000000e+00 3.172332e-57  1.817263e-75  2.208075e-58  6.585031e-78
## GR       3.172332e-57 0.000000e+00  3.413147e-40  5.844139e-33  1.401674e-41
## FW       1.817263e-75 3.413147e-40  0.000000e+00 1.906058e-106 5.460401e-181
## DW       2.208075e-58 5.844139e-33 1.906058e-106  0.000000e+00  3.560500e-88
## WC       6.585031e-78 1.401674e-41 5.460401e-181  3.560500e-88  0.000000e+00
## day6_FvP 7.639003e-04 4.745793e-04  1.538557e-03  3.928206e-03  1.240347e-03
p.mat2 <- cor_pmat(Drought)
head(p.mat2[, 1:5])
##              Area.SUM           GR            FW            DW            WC
## Area.SUM 0.000000e+00 3.172332e-57  1.817263e-75  2.208075e-58  6.585031e-78
## GR       3.172332e-57 0.000000e+00  3.413147e-40  5.844139e-33  1.401674e-41
## FW       1.817263e-75 3.413147e-40  0.000000e+00 1.906058e-106 5.460401e-181
## DW       2.208075e-58 5.844139e-33 1.906058e-106  0.000000e+00  3.560500e-88
## WC       6.585031e-78 1.401674e-41 5.460401e-181  3.560500e-88  0.000000e+00
## day6_FvP 7.639003e-04 4.745793e-04  1.538557e-03  3.928206e-03  1.240347e-03
Correlation_Drought <- ggcorrplot(
  corr2, # Correlation matrix to plot
  p.mat = p.mat2, # Optional p-values matrix for significance indicators
  type = "lower", # Type of plot (upper triangle)
  outline.col = "white", # Color of the outline around the plot
  colors = c("#6D9EC1", "white", "#E46726"), # Color scale for correlations
)

Correlation_Drought

library(cowplot)

pdf("Figure.3.Natural_diversity_Drought_Responses.pdf", width = 20, height = 15 )
plot_grid(Area_lgraph, GR_overall, GR_STI,
          ET_lgraph_spline, EVT_med_overall,
          FW_overall, Temp_plot, labels = "AUTO")
dev.off()
## quartz_off_screen 
##                 2
pdf("SFigure.Correlation_btw_traits_Control_and_Drought.pdf", width = 15, height = 7 )
plot_grid(Correlation_Control, Correlation_Drought, labels = "AUTO")
dev.off()
## quartz_off_screen 
##                 2
DW_vs_Area <- DW_vs_Area + rremove("legend")
FW_vs_Area <- FW_vs_Area + rremove("legend") 
pdf("SFigure.Correlation_FW_DW_vs_Area.pdf", width = 12, height = 7) 
plot_grid(FW_vs_Area, DW_vs_Area , labels = "AUTO")
dev.off()
## quartz_off_screen 
##                 2