Cowpea miniCore screen - overall
This notebook combines all the data collected for cowpea experiments - starting with side-view images, FW/DW, PhotoSynQ data, and evapotranspiration. It is meant to perform final examination, as well as relationship between the recorded traits.
Side-view data
Here, I combine the data for sideview imaging from all of the experiments, and calculate the overall growth rate, as well as determine which plants will be excluded based on the growth rate / missing datapoints.
Combining the experiments
SV_data <- list.files(pattern="Clean_data.csv")
SV_temp <- read.csv(SV_data[1])
SV_temp$exp <- 1
SV_all <- SV_temp
SV_temp <- read.csv(SV_data[2])
SV_temp$exp <- 2
SV_all <- rbind(SV_all, SV_temp)
SV_temp <- read.csv(SV_data[3])
SV_temp$exp <- 3
SV_all <- rbind(SV_all, SV_temp)
SV_temp <- read.csv(SV_data[4])
SV_temp$exp <- 4
SV_all <- rbind(SV_all, SV_temp)
SV_temp <- read.csv(SV_data[5])
SV_temp$exp <- 5
SV_all <- rbind(SV_all, SV_temp)
SV_temp <- read.csv(SV_data[6])
SV_temp$exp <- 6
SV_all <- rbind(SV_all, SV_temp)
SV_all$Plant.ID <- paste(SV_all$exp, SV_all$pot.no, SV_all$Genotype, SV_all$Treatment, sep="_")
SV_allBecause we have a non-overlapping days between experiment 1 and all the other experiment - in order to have the same amount of data for each day - we need to model each plant using the smooth spline.
days <- 0:14
temp <- subset(SV_all, SV_all$Plant.ID == unique(SV_all$Plant.ID)[1])
temp$day <- as.numeric(as.character(temp$day))
temp$day## [1] 0 2 7 5 11 14 9temp$area.sum## [1] 160323 300031 301616 286674 563936 862429 443940plot.spl <- with(temp, smooth.spline(day, area.sum, df = 4))
plot(area.sum ~ day, data = temp)
lines(plot.spl, col = "blue")
lines(predict(plot.spl, days), col="red")
Now - let’s create a dataframe to keep all of the sline-fitted data:
names <- c(text="ID", "genotype", "treatment", "day", "Area.SUM")
spline_data <- data.frame()
for (k in names) spline_data[[k]] <- as.character()
spline_dataAnd add the predicted values into the dataframe:
pred_temp <- predict(plot.spl, days)
pred_temp## $x
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
##
## $y
## [1] 186517.2 220133.7 247923.8 265975.2 278126.1 290152.3 307698.1 335882.1
## [9] 378420.6 433946.5 500194.6 576392.1 661674.5 753318.3 848135.0# make empty dataframe
spline_data[1:15,4] <- pred_temp$x
spline_data[1:15,5] <- pred_temp$y
spline_data[1:15,1] <- temp$Plant.ID[1]
spline_data[1:15,2] <- temp$Genotype[1]
spline_data[1:15,3] <- temp$Treatment[1]
spline_datafinal_spline <- spline_data
all_plants <- unique(SV_all$Plant.ID)
length(all_plants)## [1] 830Now - let’s loop it:
for(i in 2:830){
temp <- subset(SV_all, SV_all$Plant.ID == all_plants[i])
if(dim(temp)[1]>4){
plot.spl <- with(temp, smooth.spline(day, area.sum, df = 4))
pred_temp <- predict(plot.spl, days)
spline_data[1:15,4] <- pred_temp$x
spline_data[1:15,5] <- pred_temp$y
spline_data[1:15,1] <- temp$Plant.ID[1]
spline_data[1:15,2] <- temp$Genotype[1]
spline_data[1:15,3] <- temp$Treatment[1]
final_spline <- rbind(final_spline, spline_data)
}
else{
spline_data[1:15,4] <- days
spline_data[1:15,5] <- 0
spline_data[1:15,1] <- temp$Plant.ID[1]
spline_data[1:15,2] <- temp$Genotype[1]
spline_data[1:15,3] <- temp$Treatment[1]
}
}
final_splinefinal_spline$day <- as.numeric(as.character(final_spline$day))
final_spline$Area.SUM <- as.numeric(as.character(final_spline$Area.SUM))
final_spline <- subset(final_spline, final_spline$Area.SUM > 1)
length(unique(final_spline$ID))## [1] 827Growth plots:
library(ggplot2)
library(ggpubr)
library("ggsci")
final_spline$Area.SUM <- as.numeric(as.character(final_spline$Area.SUM))
final_spline$day <- as.factor(as.numeric(as.character(final_spline$day)))
Area_lgraph <- ggplot(data=final_spline, aes(x= day, y=Area.SUM, group = ID, color = treatment))
Area_lgraph <- Area_lgraph + geom_line(alpha = 0.1)
Area_lgraph <- Area_lgraph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= treatment), alpha=0.3)
Area_lgraph <- Area_lgraph + stat_summary(fun=mean, aes(group= treatment), size=0.7, geom="line", linetype = "dashed")
Area_lgraph <- Area_lgraph + stat_compare_means(aes(group = treatment), label = "p.signif", method = "t.test", hide.ns = T)
Area_lgraph <- Area_lgraph + ylab("Shoot Size (7 x SV pixels)") + xlab("Days After Stress") + scale_color_jco()
Area_lgraph
Calculating running growth rate
Now we got some spline fitted data - we can try and calculate growth rate that is moving throughout the experiment - to assess how the plant growth dynamics is changing.
First - let’s isolate one plant:
final_spline$day <- as.numeric(as.character(final_spline$day))
temp <- subset(final_spline, final_spline$ID == all_plants[1])
tempThe total number of timepoints for which we will be able to extract growth rate are number of days - 2 > since we will need at least 3 timepoints for each growth rate calculation:
days## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14growth_timeline <- days[1:13]
growth_timeline## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12Then - let’s subset only first three timepoints:
temp_now <- subset(temp, temp$day < 3)
temp_nowNow - let’s fit the linear growth rate to observed changes in shoot size:
plot(temp_now$Area.SUM ~ temp_now$day) + abline(lm(temp_now$Area.SUM ~ temp_now$day))
## integer(0)model_now <- lm(temp_now$Area.SUM ~ temp_now$day)
summary(model_now)##
## Call:
## lm(formula = temp_now$Area.SUM ~ temp_now$day)
##
## Residuals:
## 1 2 3
## -971.1 1942.1 -971.1
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 187488 2171 86.34 0.00737 **
## temp_now$day 30703 1682 18.25 0.03484 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2379 on 1 degrees of freedom
## Multiple R-squared: 0.997, Adjusted R-squared: 0.994
## F-statistic: 333.2 on 1 and 1 DF, p-value: 0.03484GR <- model_now$coefficients[2]Now - let’s save it into an empty dataframe
names <- c(text="ID", "genotype", "treatment", "day", "GR")
growth_data <- data.frame()
for (k in names) growth_data[[k]] <- as.character()
growth_datagrowth_data[1,1] <- temp_now$ID[1]
growth_data[1,2] <- temp_now$genotype[1]
growth_data[1,3] <- temp_now$treatment[1]
growth_data[1,4] <- min(temp_now$day)
growth_data[1,5] <- GR
growth_dataNow - let’s circle through all the possible timepoints for that plant:
for(t in 1:length(growth_timeline)){
min = t-1
max = min + 3
temp_now <- subset(temp, temp$day < max)
temp_now <- subset(temp_now, temp_now$day >= min)
temp_now
model_now <- lm(temp_now$Area.SUM ~ temp_now$day)
growth_data[t,1] <- temp_now$ID[1]
growth_data[t,2] <- temp_now$genotype[1]
growth_data[t,3] <- temp_now$treatment[1]
growth_data[t,4] <- min(temp_now$day)
growth_data[t,5] <- model_now$coefficients[2]
}
growth_dataOK - now let’s loop it for ALL the plants within the experiment:
counter = 1
for(r in 1:length(all_plants)){
temp <- subset(final_spline, final_spline$ID == all_plants[r])
if(dim(temp)[1]>0){
for(t in 1:length(growth_timeline)){
min = (t-1)
max = min + 3
temp_now <- subset(temp, temp$day < max)
temp_now <- subset(temp_now, temp_now$day >= min)
model_now <- lm(temp_now$Area.SUM ~ temp_now$day)
growth_data[counter,1] <- temp_now$ID[1]
growth_data[counter,2] <- temp_now$genotype[1]
growth_data[counter,3] <- temp_now$treatment[1]
growth_data[counter,4] <- min(temp_now$day)
growth_data[counter,5] <- model_now$coefficients[2]
counter <- counter + 1
}
}
}
growth_datalength(unique(growth_data$ID))## [1] 827Awesome! Now - let’s see how this looks like:
growth_data$GR <- as.numeric(as.character(growth_data$GR))
growth_data$day <- as.factor(as.numeric(as.character(growth_data$day)))
Daily_GR_lgraph <- ggplot(data=growth_data, aes(x= day, y=GR, group = ID, color = treatment))
Daily_GR_lgraph <- Daily_GR_lgraph + geom_line(alpha = 0.1)
Daily_GR_lgraph <- Daily_GR_lgraph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= treatment), alpha=0.3)
Daily_GR_lgraph <- Daily_GR_lgraph + stat_summary(fun=mean, aes(group= treatment), size=0.7, geom="line", linetype = "dashed")
Daily_GR_lgraph <- Daily_GR_lgraph + stat_compare_means(aes(group = treatment), label = "p.signif", method = "t.test",
hide.ns = T)
Daily_GR_lgraph <- Daily_GR_lgraph + ylab("Growth Rate (pix / day)") + xlab("Days After Stress") + scale_color_jco()
Daily_GR_lgraph## Warning: Removed 1 rows containing non-finite values (stat_summary).
## Removed 1 rows containing non-finite values (stat_summary).## Warning: Removed 1 rows containing non-finite values (stat_compare_means).## Warning: Removed 1 row(s) containing missing values (geom_path).
OK - the daily growth rate doesn’t seem to change much throughout the
observed period - and it leads us to have WAY too many -ve growth rates
- so let’s keep it simple and calculate the growth rate for each line
just over the entire experimental period based on all of the time
points.
Calculating overall growth rate
final_spline$day <- as.numeric(as.character(final_spline$day))
temp <- subset(final_spline, final_spline$ID == all_plants[1])
tempplot(temp$Area.SUM ~ temp$day) + abline(lm(temp$Area.SUM ~ temp$day))
## integer(0)model <- lm(temp$Area.SUM ~ temp$day)
GR <- model$coefficients[2]
R2 <- summary(model)$r.squarednames <- c(text="ID", "genotype", "treatment", "GR", "R2")
growth_data <- data.frame()
for (k in names) growth_data[[k]] <- as.character()
growth_datagrowth_data[1,1] <- temp$ID[1]
growth_data[1,2] <- temp$genotype[1]
growth_data[1,3] <- temp$treatment[1]
growth_data[1,4] <- GR
growth_data[1,5] <- R2
growth_dataCool - now let’s loop it for all the plants within the experiment:
all_plants <- unique(final_spline$ID)
for(i in 1:length(all_plants)){
temp <- subset(final_spline, final_spline$ID == all_plants[i])
model <- lm(temp$Area.SUM ~ temp$day)
GR <- model$coefficients[2]
R2 <- summary(model)$r.squared
growth_data[i,1] <- temp$ID[1]
growth_data[i,2] <- temp$genotype[1]
growth_data[i,3] <- temp$treatment[1]
growth_data[i,4] <- GR
growth_data[i,5] <- R2
}
growth_data Let’s see what is the distribution of the accuracy - to determine a sane R2 cutoff:
growth_data$GR <- as.numeric(as.character(growth_data$GR))
growth_data$R2 <- as.numeric(as.character(growth_data$R2))
gghistogram(growth_data, x = "R2", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(growth_data, x = "GR", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
Seems like most of the lines have pretty good fit. Let’s do a cutoff for everyone below 0.7
growth_data_clean <- subset(growth_data, growth_data$R2 > 0.7)
gghistogram(growth_data_clean, x = "GR", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
OK - now let’s look at how this overall growth rate looks like when we split it per treatment:
GR_overall <- ggerrorplot(growth_data_clean, y = "GR", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="Growth Rate (pix / day)")
GR_overall <- GR_overall + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
GR_overall <- GR_overall + scale_color_jco()
GR_overall
Now - let’s also calculate the Stress Index for the overall growth rate for each accession:
library(doBy)
growth_data_clean2 <- growth_data_clean[,2:4]
growth_sum <- summaryBy(GR ~ treatment + genotype, data = growth_data_clean2)
growth_sumlibrary(reshape2)
growth_cast <- dcast(growth_sum, genotype ~ treatment)## Using GR.mean as value column: use value.var to override.colnames(growth_cast)[2] <- "GR.Control"
colnames(growth_cast)[3] <- "GR.Drought"
growth_cast$GR.Control <- as.numeric(as.character(growth_cast$GR.Control))
growth_cast$GR.Drought <- as.numeric(as.character(growth_cast$GR.Drought))
growth_cast$GR.STI <- growth_cast$GR.Drought / growth_cast$GR.Control
growth_castgghistogram(growth_cast, x = "GR.STI", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.## Warning: Removed 49 rows containing non-finite values (stat_bin).
Awesome - before we move on - let’s combine the last day of SV imaging and see the correlations with FW, DW and WC:
Fresh weight data
FW_data <- list.files(pattern="FWDW")
FW_data## [1] "FWDW_exp1.csv" "FWDW_exp2.csv" "FWDW_exp3.csv" "FWDW_exp4.csv"
## [5] "FWDW_exp5.csv" "FWDW_exp6.csv"FW_temp <- read.csv(FW_data[1])
FW_temp$exp <- 1
FW_all <- FW_temp
FW_temp <- read.csv(FW_data[2])
FW_temp$exp <- 2
FW_all <- rbind(FW_all, FW_temp)
FW_temp <- read.csv(FW_data[3])
FW_temp$exp <- 3
FW_all <- rbind(FW_all, FW_temp)
FW_temp <- read.csv(FW_data[4])
FW_temp$exp <- 4
FW_all <- rbind(FW_all, FW_temp)
FW_temp <- read.csv(FW_data[5])
FW_temp$exp <- 5
FW_all <- rbind(FW_all, FW_temp)
FW_temp <- read.csv(FW_data[6])
FW_temp$exp <- 6
FW_all <- rbind(FW_all, FW_temp)
FW_all$Treatment <- gsub(" ", "", FW_all$Treatment)
FW_all$ID <- paste(FW_all$exp, FW_all$Pot.number, FW_all$Genotype, FW_all$Treatment, sep="_")
FW_allfinal_splinelast_day <- subset(final_spline, final_spline$day == 14)
colnames(FW_all)[2] <- "genotype"
colnames(FW_all)[3] <- "treatment"
FW_last_day <- merge(last_day, FW_all, id=c("ID", "genotype", "treatment"))
FW_last_dayAlso - let’s take advantage of all the data curation we did with the growth rate - and let’s consider only the plants that had good growth fit:
clean_FW <- subset(FW_last_day, FW_last_day$ID %in% growth_data_clean$ID)
clean_FWclean_FW$FW <- as.numeric(as.character(clean_FW$FW))
clean_FW$DW <- as.numeric(as.character(clean_FW$DW))
clean_FW$WC <- as.numeric(as.character(clean_FW$WC))
FW_Area <- ggscatter(clean_FW, x = "Area.SUM", y = "FW", color = "treatment", rug = TRUE) + stat_cor() + scale_color_jco()
FW_Area
DW_Area <- ggscatter(clean_FW, x = "Area.SUM", y = "DW", color = "treatment", rug = TRUE) + stat_cor() + scale_color_jco()
DW_Area
WC_Area <- ggscatter(clean_FW, x = "Area.SUM", y = "WC", color = "treatment", rug = TRUE) + stat_cor() + scale_color_jco()
WC_Area
PhotosynQ Data
PSQ_all <- list.files(pattern = "photosynq")
PSQ_all## [1] "photosynq_exp1.csv" "photosynq_exp2.csv" "photosynq_exp3.csv"
## [4] "photosynq_exp4.csv" "photosynq_exp5.csv" "photosynq_exp6.csv"PSQ experiment 1
PSQ1 <- read.csv(PSQ_all[1])
colnames(PSQ1)## [1] "ID" "Series"
## [3] "Repeat" "Pot.Number"
## [5] "Treatment" "air_temp_kinetics"
## [7] "Ambient.Humidity" "Ambient.Pressure"
## [9] "Ambient.Temperature" "data_raw_PAM"
## [11] "ECS_averaged_trace" "ECS_tau"
## [13] "ECSt.mAU" "fitinput"
## [15] "FmPrime" "FoPrime"
## [17] "Fs" "FvP_over_FmP"
## [19] "gH." "humidity2_K"
## [21] "humidity_K" "kP700"
## [23] "leaf.angle" "Leaf.Temperature"
## [25] "Leaf.Temperature.Differenial" "Leaf.Temperature.Differential"
## [27] "LEAF_temp" "leaf_thickness"
## [29] "LEF" "LEFd_trace"
## [31] "Light.Intensity..PAR." "NPQt"
## [33] "outdata" "P700_DIRK_ampl"
## [35] "P700_DIRK_averaged_trace" "P700_fitinput"
## [37] "P700_outdata" "Phi2"
## [39] "PhiNO" "PhiNPQ"
## [41] "PS1.Active.Centers" "PS1.Open.Centers"
## [43] "PS1.Over.Reduced.Centers" "PS1.Oxidized.Centers"
## [45] "PSI_data_absorbance" "pump"
## [47] "qL" "SPAD"
## [49] "test_data_raw_PAM" "time"
## [51] "Time.of.Day" "tP700"
## [53] "v_initial_P700" "vH."
## [55] "User" "Device.ID"
## [57] "Latitude" "Longitude"
## [59] "Issues"PSQ1 <- PSQ1[,c(4,50, 15:16, 18, 24, 28, 32, 40:42, 48)]
unique(PSQ1$time)## [1] "3/16/2022 0:00" "3/9/2022 0:00"PSQ1$time<- gsub("3/9/2022 0:00", "day6", PSQ1$time)
PSQ1$time<- gsub("3/16/2022 0:00", "day13", PSQ1$time)
gghistogram(PSQ1, x = "FvP_over_FmP", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ1, x = "Leaf.Temperature", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ1, x = "leaf_thickness", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
PSQ1 <- subset(PSQ1, PSQ1$leaf_thickness < 2)
PSQ1mPSQ1 <- melt(PSQ1, id=c("Pot.Number", "time"))
mPSQ1PSQ1_cast <- dcast(mPSQ1, Pot.Number ~ time + variable)
PSQ1_cast$exp <- 1
PSQ1_castPSQ experiment 2
PSQ2 <- read.csv(PSQ_all[2])
colnames(PSQ2)## [1] "Datum.ID" "Repeat"
## [3] "Pot.Number" "Treatment"
## [5] "time" "leaf.angle"
## [7] "test_data_raw_PAM" "pump"
## [9] "ECS_averaged_trace" "fitinput"
## [11] "outdata" "ECSt.mAU"
## [13] "ECS_tau" "gH."
## [15] "vH." "P700_DIRK_averaged_trace"
## [17] "P700_fitinput" "P700_outdata"
## [19] "P700_DIRK_ampl" "tP700"
## [21] "kP700" "v_initial_P700"
## [23] "LEFd_trace" "data_raw_PAM"
## [25] "Fs" "FoPrime"
## [27] "Phi2" "PhiNPQ"
## [29] "qL" "NPQt"
## [31] "PhiNO" "FvP_over_FmP"
## [33] "FmPrime" "PSI_data_absorbance"
## [35] "PS1.Active.Centers" "PS1.Open.Centers"
## [37] "PS1.Over.Reduced.Centers" "PS1.Oxidized.Centers"
## [39] "humidity_K" "humidity2_K"
## [41] "air_temp_kinetics" "leaf_thickness"
## [43] "LEAF_temp" "Light.Intensity..PAR."
## [45] "Ambient.Temperature" "Ambient.Humidity"
## [47] "Ambient.Pressure" "Leaf.Temperature"
## [49] "Leaf.Temperature.Differential" "LEF"
## [51] "Leaf.Temperature.Differenial" "SPAD"
## [53] "User" "Device.ID"
## [55] "Status" "Notes"
## [57] "Latitude" "Longitude"PSQ2 <- PSQ2[,c(3, 5, 33, 26, 32, 48, 42, 30, 28, 35:36, 52)]
unique(PSQ2$time)## [1] "4/20/2022" "4/13/2022"PSQ2$time<- gsub("4/13/2022", "day6", PSQ2$time)
PSQ2$time<- gsub("4/20/2022", "day13", PSQ2$time)
gghistogram(PSQ2, x = "FvP_over_FmP", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ2, x = "Leaf.Temperature", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ2, x = "leaf_thickness", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
mPSQ2 <- melt(PSQ2, id=c("Pot.Number", "time"))
mPSQ2PSQ2_cast <- dcast(mPSQ2, Pot.Number ~ time + variable)
PSQ2_cast$exp <- 2
PSQ2_castPSQ3 <- read.csv(PSQ_all[3])
colnames(PSQ3)## [1] "Datum.ID" "Repeat"
## [3] "Pot.Number" "Treatment"
## [5] "time" "leaf.angle"
## [7] "test_data_raw_PAM" "pump"
## [9] "ECS_averaged_trace" "fitinput"
## [11] "outdata" "ECSt.mAU"
## [13] "ECS_tau" "gH."
## [15] "vH." "P700_DIRK_averaged_trace"
## [17] "P700_fitinput" "P700_outdata"
## [19] "P700_DIRK_ampl" "tP700"
## [21] "kP700" "v_initial_P700"
## [23] "LEFd_trace" "data_raw_PAM"
## [25] "Fs" "FoPrime"
## [27] "Phi2" "PhiNPQ"
## [29] "qL" "NPQt"
## [31] "PhiNO" "FvP_over_FmP"
## [33] "FmPrime" "PSI_data_absorbance"
## [35] "PS1.Active.Centers" "PS1.Open.Centers"
## [37] "PS1.Over.Reduced.Centers" "PS1.Oxidized.Centers"
## [39] "humidity_K" "humidity2_K"
## [41] "air_temp_kinetics" "leaf_thickness"
## [43] "LEAF_temp" "Light.Intensity..PAR."
## [45] "Ambient.Temperature" "Ambient.Humidity"
## [47] "Ambient.Pressure" "Leaf.Temperature"
## [49] "Leaf.Temperature.Differential" "LEF"
## [51] "Leaf.Temperature.Differenial" "SPAD"
## [53] "User" "Device.ID"
## [55] "Status" "Notes"
## [57] "Latitude" "Longitude"PSQ3 <- PSQ3[,c(3, 5, 33, 26, 32, 48, 42, 30, 28, 35:36, 52)]
unique(PSQ3$time)## [1] "5/4/2022" "4/27/2022"PSQ3$time<- gsub("4/27/2022", "day6", PSQ3$time)
PSQ3$time<- gsub("5/4/2022", "day13", PSQ3$time)
gghistogram(PSQ3, x = "FvP_over_FmP", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ3, x = "Leaf.Temperature", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ3, x = "leaf_thickness", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
mPSQ3 <- melt(PSQ3, id=c("Pot.Number", "time"))
mPSQ3PSQ3_cast <- dcast(mPSQ3, Pot.Number ~ time + variable)
PSQ3_cast$exp <- 3
PSQ3_castPSQ4 <- read.csv(PSQ_all[4])
colnames(PSQ4)## [1] "ID" "Series"
## [3] "Repeat" "Pot.Number"
## [5] "Treatment" "air_temp_kinetics"
## [7] "Ambient.Humidity" "Ambient.Pressure"
## [9] "Ambient.Temperature" "data_raw_PAM"
## [11] "ECS_averaged_trace" "ECS_tau"
## [13] "ECSt.mAU" "fitinput"
## [15] "FmPrime" "FoPrime"
## [17] "Fs" "FvP_over_FmP"
## [19] "gH." "humidity2_K"
## [21] "humidity_K" "kP700"
## [23] "leaf.angle" "Leaf.Temperature"
## [25] "Leaf.Temperature.Differenial" "Leaf.Temperature.Differential"
## [27] "LEAF_temp" "leaf_thickness"
## [29] "LEF" "LEFd_trace"
## [31] "Light.Intensity..PAR." "NPQt"
## [33] "outdata" "P700_DIRK_ampl"
## [35] "P700_DIRK_averaged_trace" "P700_fitinput"
## [37] "P700_outdata" "Phi2"
## [39] "PhiNO" "PhiNPQ"
## [41] "PS1.Active.Centers" "PS1.Open.Centers"
## [43] "PS1.Over.Reduced.Centers" "PS1.Oxidized.Centers"
## [45] "PSI_data_absorbance" "pump"
## [47] "qL" "SPAD"
## [49] "test_data_raw_PAM" "time"
## [51] "Time.of.Day" "tP700"
## [53] "v_initial_P700" "vH."
## [55] "User" "Device.ID"
## [57] "Latitude" "Longitude"
## [59] "Issues"PSQ4 <- PSQ4[,c(4,50, 15:16, 18, 24, 28, 32, 40:42, 48)]
unique(PSQ4$time)## [1] "5/18/2022 15:05" "5/18/2022 15:06" "5/18/2022 15:07" "5/18/2022 15:08"
## [5] "5/18/2022 15:09" "5/18/2022 15:10" "5/18/2022 15:11" "5/18/2022 15:12"
## [9] "5/18/2022 15:13" "5/18/2022 15:14" "5/18/2022 15:15" "5/18/2022 15:16"
## [13] "5/18/2022 15:17" "5/18/2022 15:18" "5/18/2022 15:19" "5/18/2022 15:20"
## [17] "5/18/2022 15:21" "5/18/2022 15:22" "5/18/2022 15:24" "5/18/2022 15:25"
## [21] "5/18/2022 15:26" "5/18/2022 15:27" "5/18/2022 15:28" "5/18/2022 15:29"
## [25] "5/18/2022 15:30" "5/18/2022 15:31" "5/18/2022 15:32" "5/18/2022 15:33"
## [29] "5/18/2022 15:34" "5/18/2022 15:35" "5/18/2022 15:36" "5/18/2022 15:37"
## [33] "5/18/2022 15:38" "5/18/2022 15:39" "5/18/2022 15:40" "5/18/2022 15:41"
## [37] "5/18/2022 15:42" "5/18/2022 15:43" "5/18/2022 15:44" "5/18/2022 15:45"
## [41] "5/18/2022 15:46" "5/18/2022 15:47" "5/18/2022 15:48" "5/18/2022 15:49"
## [45] "5/18/2022 15:50" "5/18/2022 15:51" "5/18/2022 15:52" "5/18/2022 15:53"
## [49] "5/18/2022 15:54" "5/18/2022 15:55" "5/18/2022 15:56" "5/18/2022 15:57"
## [53] "5/18/2022 15:58" "5/18/2022 15:59" "5/18/2022 16:00" "5/18/2022 16:01"
## [57] "5/18/2022 16:02" "5/18/2022 16:03" "5/18/2022 16:04" "5/18/2022 16:05"
## [61] "5/18/2022 16:06" "5/18/2022 16:07" "5/18/2022 16:08" "5/18/2022 16:09"
## [65] "5/18/2022 16:10" "5/18/2022 16:11" "5/18/2022 16:12" "5/18/2022 16:13"
## [69] "5/18/2022 16:14" "5/18/2022 16:15" "5/18/2022 14:56" "5/18/2022 14:55"
## [73] "5/18/2022 14:54" "5/18/2022 14:53" "5/18/2022 14:52" "5/18/2022 14:47"
## [77] "5/18/2022 14:46" "5/18/2022 14:45" "5/18/2022 14:44" "5/18/2022 14:43"
## [81] "5/18/2022 14:42" "5/18/2022 14:41" "5/18/2022 14:40" "5/11/2022 15:20"
## [85] "5/11/2022 15:19" "5/11/2022 15:18" "5/11/2022 15:17" "5/11/2022 15:16"
## [89] "5/11/2022 15:15" "5/11/2022 15:14" "5/11/2022 15:13" "5/11/2022 15:11"
## [93] "5/11/2022 15:10" "5/11/2022 15:09" "5/11/2022 15:08" "5/11/2022 15:07"
## [97] "5/11/2022 15:06" "5/11/2022 15:05" "5/11/2022 15:04" "5/11/2022 15:03"
## [101] "5/11/2022 15:02" "5/11/2022 15:01" "5/11/2022 15:00" "5/11/2022 14:59"
## [105] "5/11/2022 14:58" "5/11/2022 14:57" "5/11/2022 14:56" "5/11/2022 14:55"
## [109] "5/11/2022 14:54" "5/11/2022 14:53" "5/11/2022 14:52" "5/11/2022 14:51"
## [113] "5/11/2022 14:50" "5/11/2022 14:49" "5/11/2022 14:48" "5/11/2022 14:47"
## [117] "5/11/2022 14:46" "5/11/2022 14:45" "5/11/2022 14:44" "5/11/2022 14:43"
## [121] "5/11/2022 14:42" "5/11/2022 14:41" "5/11/2022 14:40" "5/11/2022 14:39"
## [125] "5/11/2022 14:38" "5/11/2022 14:37" "5/11/2022 14:36" "5/11/2022 14:35"
## [129] "5/11/2022 14:34" "5/11/2022 14:33" "5/11/2022 14:32" "5/11/2022 14:31"
## [133] "5/11/2022 14:30" "5/11/2022 14:29" "5/11/2022 14:28" "5/11/2022 14:27"
## [137] "5/11/2022 14:26" "5/11/2022 14:25" "5/11/2022 14:23" "5/11/2022 14:22"
## [141] "5/11/2022 14:21" "5/11/2022 14:20" "5/11/2022 14:19" "5/11/2022 14:14"
## [145] "5/11/2022 14:13" "5/11/2022 14:12" "5/11/2022 14:11" "5/11/2022 14:10"
## [149] "5/11/2022 14:09" "5/11/2022 14:08" "5/11/2022 14:07" "5/11/2022 14:06"
## [153] "5/11/2022 14:05" "5/11/2022 14:04" "5/11/2022 14:03" "5/11/2022 14:02"
## [157] "5/11/2022 14:01" "5/11/2022 14:00" "5/11/2022 13:59"for(i in 1:nrow(PSQ4)){
PSQ4$time[i] <- strsplit(PSQ4$time[i], " ")[[1]][1]
}
unique(PSQ4$time)## [1] "5/18/2022" "5/11/2022"PSQ4$time<- gsub("5/11/2022", "day6", PSQ4$time)
PSQ4$time<- gsub("5/18/2022", "day13", PSQ4$time)
gghistogram(PSQ4, x = "FvP_over_FmP", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ4, x = "Leaf.Temperature", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ4, x = "leaf_thickness", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
mPSQ4 <- melt(PSQ4, id=c("Pot.Number", "time"))
mPSQ4PSQ4_cast <- dcast(mPSQ4, Pot.Number ~ time + variable)
PSQ4_cast$exp <- 4
PSQ4_castPSQ5 <- read.csv(PSQ_all[5])
colnames(PSQ5)## [1] "ID" "Series"
## [3] "Repeat" "Pot.number"
## [5] "Treatment" "air_temp_kinetics"
## [7] "Ambient.Humidity" "Ambient.Pressure"
## [9] "Ambient.Temperature" "data_raw_PAM"
## [11] "ECS_averaged_trace" "ECS_tau"
## [13] "ECSt.mAU" "fitinput"
## [15] "FmPrime" "FoPrime"
## [17] "Fs" "FvP_over_FmP"
## [19] "gH." "humidity2_K"
## [21] "humidity_K" "kP700"
## [23] "leaf.angle" "Leaf.Temperature"
## [25] "Leaf.Temperature.Differenial" "Leaf.Temperature.Differential"
## [27] "LEAF_temp" "leaf_thickness"
## [29] "LEF" "LEFd_trace"
## [31] "Light.Intensity..PAR." "NPQt"
## [33] "outdata" "P700_DIRK_ampl"
## [35] "P700_DIRK_averaged_trace" "P700_fitinput"
## [37] "P700_outdata" "Phi2"
## [39] "PhiNO" "PhiNPQ"
## [41] "PS1.Active.Centers" "PS1.Open.Centers"
## [43] "PS1.Over.Reduced.Centers" "PS1.Oxidized.Centers"
## [45] "PSI_data_absorbance" "pump"
## [47] "qL" "SPAD"
## [49] "test_data_raw_PAM" "time"
## [51] "Time.of.Day" "tP700"
## [53] "v_initial_P700" "vH."
## [55] "User" "Device.ID"
## [57] "Latitude" "Longitude"
## [59] "Issues"PSQ5 <- PSQ5[,c(4,50, 15:16, 18, 24, 28, 32, 40:42, 48)]
unique(PSQ5$time)## [1] "6/22/2022 15:39" "6/22/2022 15:38" "6/22/2022 15:37" "6/22/2022 15:36"
## [5] "6/22/2022 15:35" "6/22/2022 15:34" "6/22/2022 15:33" "6/22/2022 15:32"
## [9] "6/22/2022 15:31" "6/22/2022 15:30" "6/22/2022 15:29" "6/22/2022 15:28"
## [13] "6/22/2022 15:27" "6/22/2022 15:26" "6/22/2022 15:25" "6/22/2022 15:24"
## [17] "6/22/2022 15:23" "6/22/2022 15:22" "6/22/2022 15:21" "6/22/2022 15:20"
## [21] "6/22/2022 15:19" "6/22/2022 15:18" "6/22/2022 15:17" "6/22/2022 15:16"
## [25] "6/22/2022 15:15" "6/22/2022 15:14" "6/22/2022 15:13" "6/22/2022 15:12"
## [29] "6/22/2022 15:11" "6/22/2022 15:10" "6/22/2022 15:09" "6/22/2022 15:08"
## [33] "6/22/2022 15:07" "6/22/2022 15:06" "6/22/2022 15:05" "6/22/2022 15:04"
## [37] "6/22/2022 15:02" "6/22/2022 15:01" "6/22/2022 15:00" "6/22/2022 14:59"
## [41] "6/22/2022 14:58" "6/22/2022 14:57" "6/22/2022 14:56" "6/22/2022 14:55"
## [45] "6/22/2022 14:54" "6/22/2022 14:53" "6/22/2022 14:52" "6/22/2022 14:51"
## [49] "6/22/2022 14:50" "6/22/2022 14:49" "6/22/2022 14:48" "6/22/2022 14:47"
## [53] "6/22/2022 14:46" "6/22/2022 14:44" "6/22/2022 14:41" "6/22/2022 14:40"
## [57] "6/22/2022 14:39" "6/22/2022 14:37" "6/22/2022 14:36" "6/22/2022 14:34"
## [61] "6/22/2022 14:33" "6/22/2022 14:32" "6/22/2022 14:31" "6/22/2022 14:30"
## [65] "6/22/2022 14:29" "6/22/2022 14:24" "6/22/2022 14:23" "6/22/2022 14:22"
## [69] "6/22/2022 14:21" "6/22/2022 14:20" "6/22/2022 14:19" "6/22/2022 14:18"
## [73] "6/22/2022 14:17" "6/22/2022 14:16" "6/22/2022 14:13" "6/22/2022 14:12"
## [77] "6/22/2022 14:11" "6/22/2022 14:10" "6/22/2022 14:09" "6/22/2022 14:08"
## [81] "6/15/2022 15:02" "6/15/2022 15:01" "6/15/2022 15:00" "6/15/2022 14:59"
## [85] "6/15/2022 14:58" "6/15/2022 14:57" "6/15/2022 14:56" "6/15/2022 14:55"
## [89] "6/15/2022 14:54" "6/15/2022 14:53" "6/15/2022 14:52" "6/15/2022 14:51"
## [93] "6/15/2022 14:50" "6/15/2022 14:49" "6/15/2022 14:48" "6/15/2022 14:47"
## [97] "6/15/2022 14:46" "6/15/2022 14:45" "6/15/2022 14:44" "6/15/2022 14:43"
## [101] "6/15/2022 14:42" "6/15/2022 14:41" "6/15/2022 14:40" "6/15/2022 14:39"
## [105] "6/15/2022 14:38" "6/15/2022 14:37" "6/15/2022 14:36" "6/15/2022 14:35"
## [109] "6/15/2022 14:33" "6/15/2022 14:32" "6/15/2022 14:31" "6/15/2022 14:30"
## [113] "6/15/2022 14:29" "6/15/2022 14:28" "6/15/2022 14:27" "6/15/2022 14:26"
## [117] "6/15/2022 14:25" "6/15/2022 14:24" "6/15/2022 14:23" "6/15/2022 14:22"
## [121] "6/15/2022 14:21" "6/15/2022 14:20" "6/15/2022 14:19" "6/15/2022 14:18"
## [125] "6/15/2022 14:16" "6/15/2022 14:15" "6/15/2022 14:14" "6/15/2022 14:13"
## [129] "6/15/2022 14:12" "6/15/2022 14:11" "6/15/2022 14:10" "6/15/2022 14:09"
## [133] "6/15/2022 14:08" "6/15/2022 14:07" "6/15/2022 14:06" "6/15/2022 14:05"
## [137] "6/15/2022 14:04" "6/15/2022 14:03" "6/15/2022 14:02" "6/15/2022 14:01"
## [141] "6/15/2022 14:00" "6/15/2022 13:59" "6/15/2022 13:58" "6/15/2022 13:57"
## [145] "6/15/2022 13:55" "6/15/2022 13:54" "6/15/2022 13:53" "6/15/2022 13:52"
## [149] "6/15/2022 13:51" "6/15/2022 13:50" "6/15/2022 13:49" "6/15/2022 13:48"
## [153] "6/15/2022 13:47" "6/15/2022 13:46" "6/15/2022 13:45" "6/15/2022 13:44"
## [157] "6/15/2022 13:43" "6/15/2022 13:42" "6/15/2022 13:41" "6/15/2022 13:40"for(i in 1:nrow(PSQ5)){
PSQ5$time[i] <- strsplit(PSQ5$time[i], " ")[[1]][1]
}
unique(PSQ5$time)## [1] "6/22/2022" "6/15/2022"PSQ5$time<- gsub("6/15/2022", "day6", PSQ5$time)
PSQ5$time<- gsub("6/22/2022", "day13", PSQ5$time)
gghistogram(PSQ5, x = "FvP_over_FmP", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ5, x = "Leaf.Temperature", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ5, x = "leaf_thickness", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
mPSQ5 <- melt(PSQ5, id=c("Pot.number", "time"))
mPSQ5PSQ5_cast <- dcast(mPSQ5, Pot.number ~ time + variable)
PSQ5_cast$exp <- 5
PSQ5_castPSQ_all[6]## [1] "photosynq_exp6.csv"PSQ6 <- read.csv(PSQ_all[6])
colnames(PSQ6)## [1] "ID" "Series"
## [3] "Repeat" "Pot.number"
## [5] "Treatment" "air_temp_kinetics"
## [7] "Ambient.Humidity" "Ambient.Pressure"
## [9] "Ambient.Temperature" "data_raw_PAM"
## [11] "ECS_averaged_trace" "ECS_tau"
## [13] "ECSt.mAU" "fitinput"
## [15] "FmPrime" "FoPrime"
## [17] "Fs" "FvP_over_FmP"
## [19] "gH." "humidity2_K"
## [21] "humidity_K" "kP700"
## [23] "leaf.angle" "Leaf.Temperature"
## [25] "Leaf.Temperature.Differenial" "Leaf.Temperature.Differential"
## [27] "LEAF_temp" "leaf_thickness"
## [29] "LEF" "LEFd_trace"
## [31] "Light.Intensity..PAR." "NPQt"
## [33] "outdata" "P700_DIRK_ampl"
## [35] "P700_DIRK_averaged_trace" "P700_fitinput"
## [37] "P700_outdata" "Phi2"
## [39] "PhiNO" "PhiNPQ"
## [41] "PS1.Active.Centers" "PS1.Open.Centers"
## [43] "PS1.Over.Reduced.Centers" "PS1.Oxidized.Centers"
## [45] "PSI_data_absorbance" "pump"
## [47] "qL" "SPAD"
## [49] "test_data_raw_PAM" "time"
## [51] "Time.of.Day" "tP700"
## [53] "v_initial_P700" "vH."
## [55] "User" "Device.ID"
## [57] "Latitude" "Longitude"
## [59] "Issues"PSQ6 <- PSQ6[,c(4,50, 15:16, 18, 24, 28, 32, 40:42, 48)]
unique(PSQ6$time)## [1] "07/20/2022 3:23 PM" "07/20/2022 3:22 PM" "07/20/2022 3:21 PM"
## [4] "07/20/2022 3:20 PM" "07/20/2022 3:19 PM" "07/20/2022 3:18 PM"
## [7] "07/20/2022 3:17 PM" "07/20/2022 3:16 PM" "07/20/2022 1:19 PM"
## [10] "07/20/2022 1:20 PM" "07/20/2022 1:21 PM" "07/20/2022 1:22 PM"
## [13] "07/20/2022 1:23 PM" "07/20/2022 1:24 PM" "07/20/2022 1:25 PM"
## [16] "07/20/2022 1:26 PM" "07/20/2022 1:27 PM" "07/20/2022 1:28 PM"
## [19] "07/20/2022 3:03 PM" "07/20/2022 3:02 PM" "07/20/2022 3:01 PM"
## [22] "07/20/2022 3:00 PM" "07/20/2022 2:59 PM" "07/20/2022 2:57 PM"
## [25] "07/20/2022 2:55 PM" "07/20/2022 2:54 PM" "07/20/2022 2:53 PM"
## [28] "07/20/2022 2:52 PM" "07/20/2022 2:51 PM" "07/20/2022 2:50 PM"
## [31] "07/20/2022 2:49 PM" "07/20/2022 2:48 PM" "07/20/2022 2:47 PM"
## [34] "07/20/2022 2:46 PM" "07/20/2022 2:45 PM" "07/20/2022 2:44 PM"
## [37] "07/20/2022 2:43 PM" "07/20/2022 2:42 PM" "07/20/2022 2:41 PM"
## [40] "07/20/2022 2:40 PM" "07/20/2022 2:39 PM" "07/20/2022 2:38 PM"
## [43] "07/20/2022 2:37 PM" "07/20/2022 2:35 PM" "07/20/2022 2:34 PM"
## [46] "07/20/2022 2:33 PM" "07/20/2022 2:32 PM" "07/20/2022 2:31 PM"
## [49] "07/20/2022 2:30 PM" "07/20/2022 2:29 PM" "07/20/2022 2:28 PM"
## [52] "07/20/2022 2:27 PM" "07/20/2022 2:26 PM" "07/20/2022 2:25 PM"
## [55] "07/20/2022 2:24 PM" "07/20/2022 2:23 PM" "07/20/2022 2:22 PM"
## [58] "07/20/2022 2:21 PM" "07/20/2022 2:20 PM" "07/20/2022 2:19 PM"
## [61] "07/20/2022 2:18 PM" "07/20/2022 2:17 PM" "07/20/2022 2:16 PM"
## [64] "07/20/2022 2:15 PM" "07/20/2022 2:14 PM" "07/20/2022 2:13 PM"
## [67] "07/20/2022 2:12 PM" "07/20/2022 2:11 PM" "07/20/2022 2:10 PM"
## [70] "07/20/2022 1:41 PM" "07/20/2022 1:40 PM" "07/20/2022 1:38 PM"
## [73] "07/20/2022 1:37 PM" "07/20/2022 1:36 PM" "07/20/2022 1:35 PM"
## [76] "07/20/2022 1:34 PM" "07/20/2022 1:33 PM" "07/20/2022 1:32 PM"
## [79] "07/20/2022 1:31 PM" "07/20/2022 1:30 PM" "07/20/2022 1:29 PM"
## [82] "07/20/2022 1:18 PM" "07/20/2022 1:17 PM" "07/20/2022 1:16 PM"
## [85] "07/20/2022 1:13 PM" "07/20/2022 1:12 PM" "07/20/2022 1:10 PM"
## [88] "07/20/2022 1:08 PM" "07/20/2022 1:07 PM" "07/20/2022 1:06 PM"
## [91] "07/13/2022 2:03 PM" "07/13/2022 2:02 PM" "07/13/2022 2:01 PM"
## [94] "07/13/2022 2:00 PM" "07/13/2022 1:59 PM" "07/13/2022 1:58 PM"
## [97] "07/13/2022 1:57 PM" "07/13/2022 1:56 PM" "07/13/2022 1:53 PM"
## [100] "07/13/2022 1:52 PM" "07/13/2022 1:51 PM" "07/13/2022 1:50 PM"
## [103] "07/13/2022 1:49 PM" "07/13/2022 1:48 PM" "07/13/2022 1:47 PM"
## [106] "07/13/2022 1:46 PM" "07/13/2022 1:45 PM" "07/13/2022 1:44 PM"
## [109] "07/13/2022 1:43 PM" "07/13/2022 1:42 PM" "07/13/2022 1:41 PM"
## [112] "07/13/2022 1:40 PM" "07/13/2022 1:39 PM" "07/13/2022 1:37 PM"
## [115] "07/13/2022 1:36 PM" "07/13/2022 1:35 PM" "07/13/2022 1:34 PM"
## [118] "07/13/2022 1:33 PM" "07/13/2022 1:32 PM" "07/13/2022 1:31 PM"
## [121] "07/13/2022 1:30 PM" "07/13/2022 1:29 PM" "07/13/2022 1:28 PM"
## [124] "07/13/2022 1:27 PM" "07/13/2022 1:26 PM" "07/13/2022 1:25 PM"
## [127] "07/13/2022 1:24 PM" "07/13/2022 1:21 PM" "07/13/2022 1:20 PM"
## [130] "07/13/2022 1:19 PM" "07/13/2022 1:10 PM" "07/13/2022 1:09 PM"
## [133] "07/13/2022 1:08 PM" "07/13/2022 1:06 PM" "07/13/2022 1:05 PM"
## [136] "07/13/2022 1:04 PM" "07/13/2022 1:02 PM" "07/13/2022 1:01 PM"
## [139] "07/13/2022 1:00 PM" "07/13/2022 12:59 PM" "07/13/2022 12:58 PM"
## [142] "07/13/2022 12:57 PM" "07/13/2022 12:55 PM" "07/13/2022 12:54 PM"
## [145] "07/13/2022 12:52 PM" "07/13/2022 12:51 PM" "07/13/2022 12:50 PM"
## [148] "07/13/2022 12:49 PM" "07/13/2022 12:48 PM" "07/13/2022 12:46 PM"
## [151] "07/13/2022 12:45 PM" "07/13/2022 12:43 PM" "07/13/2022 12:42 PM"
## [154] "07/13/2022 12:40 PM" "07/13/2022 12:39 PM" "07/13/2022 12:38 PM"
## [157] "07/13/2022 12:36 PM" "07/13/2022 12:35 PM" "07/13/2022 12:34 PM"
## [160] "07/13/2022 12:33 PM" "07/13/2022 12:32 PM" "07/13/2022 12:31 PM"
## [163] "07/13/2022 12:30 PM" "07/13/2022 12:29 PM" "07/13/2022 12:28 PM"
## [166] "07/13/2022 12:27 PM" "07/13/2022 12:26 PM" "07/13/2022 12:25 PM"for(i in 1:nrow(PSQ6)){
PSQ6$time[i] <- strsplit(PSQ6$time[i], " ")[[1]][1]
}
unique(PSQ6$time)## [1] "07/20/2022" "07/13/2022"PSQ6PSQ6$time<- gsub(unique(PSQ6$time)[1], "day13", PSQ6$time)
PSQ6$time<- gsub(unique(PSQ6$time)[2], "day6", PSQ6$time)
PSQ6gghistogram(PSQ6, x = "FvP_over_FmP", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ6, x = "Leaf.Temperature", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
gghistogram(PSQ6, x = "leaf_thickness", fill = "lightgray",
add = "mean", rug = TRUE)## Warning: Using `bins = 30` by default. Pick better value with the argument
## `bins`.## Warning: geom_vline(): Ignoring `mapping` because `xintercept` was provided.## Warning: geom_vline(): Ignoring `data` because `xintercept` was provided.
mPSQ6 <- melt(PSQ6, id=c("Pot.number", "time"))
mPSQ6PSQ6_cast <- dcast(mPSQ6, Pot.number ~ time + variable, fun=mean)
PSQ6_cast$exp <- 6
PSQ6_castOK - let’s merge the data together.
colnames(PSQ1_cast) %in% colnames(PSQ2_cast)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUEPSQ <- rbind(PSQ1_cast, PSQ2_cast)
colnames(PSQ1_cast) %in% colnames(PSQ3_cast)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUEPSQ <- rbind(PSQ, PSQ3_cast)
colnames(PSQ1_cast) %in% colnames(PSQ4_cast)## [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUEPSQ <- rbind(PSQ, PSQ4_cast)
colnames(PSQ1_cast) %in% colnames(PSQ5_cast)## [1] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEcolnames(PSQ5_cast)[1] <- "Pot.Number"
PSQ <- rbind(PSQ, PSQ5_cast)
colnames(PSQ1_cast) %in% colnames(PSQ6_cast)## [1] FALSE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [13] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUEcolnames(PSQ6_cast)## [1] "Pot.number" "day13_FmPrime"
## [3] "day13_FoPrime" "day13_FvP_over_FmP"
## [5] "day13_Leaf.Temperature" "day13_leaf_thickness"
## [7] "day13_NPQt" "day13_PhiNPQ"
## [9] "day13_PS1.Active.Centers" "day13_PS1.Open.Centers"
## [11] "day13_SPAD" "day6_FmPrime"
## [13] "day6_FoPrime" "day6_FvP_over_FmP"
## [15] "day6_Leaf.Temperature" "day6_leaf_thickness"
## [17] "day6_NPQt" "day6_PhiNPQ"
## [19] "day6_PS1.Active.Centers" "day6_PS1.Open.Centers"
## [21] "day6_SPAD" "exp"colnames(PSQ5_cast)## [1] "Pot.Number" "day13_FmPrime"
## [3] "day13_FoPrime" "day13_FvP_over_FmP"
## [5] "day13_Leaf.Temperature" "day13_leaf_thickness"
## [7] "day13_NPQt" "day13_PhiNPQ"
## [9] "day13_PS1.Active.Centers" "day13_PS1.Open.Centers"
## [11] "day13_SPAD" "day6_FmPrime"
## [13] "day6_FoPrime" "day6_FvP_over_FmP"
## [15] "day6_Leaf.Temperature" "day6_leaf_thickness"
## [17] "day6_NPQt" "day6_PhiNPQ"
## [19] "day6_PS1.Active.Centers" "day6_PS1.Open.Centers"
## [21] "day6_SPAD" "exp"colnames(PSQ6_cast) <- colnames(PSQ5_cast)
PSQ <- rbind(PSQ, PSQ6_cast)
PSQEvapotranspiration
Now it’s time for evapo-transpiration. It will take a little time to process - so let’s look at each file individually, and fuse them together at the end. I dont think we need to get all time data - but rather sum and mean evapotranspiration per plant.
EVT_all <- list.files(pattern = "EVT")
EVT_all## [1] "EVT_exp1.csv" "EVT_exp2.csv"
## [3] "EVT_exp3.csv" "EVT_exp4.csv"
## [5] "EVT_exp5.csv" "EVT_exp6.csv"
## [7] "Fig.Cowpea2022.Growth.EVT.PSQ.pdf"ET1 <- read.csv(EVT_all[1])
ET1ET1m <- melt(ET1, id.vars = c("Pot.number", "Genotype", "Treatment"))
colnames(ET1m)[2] <- "genotype"
colnames(ET1m)[3] <- "treatment"
colnames(ET1m)[4] <- "day"
colnames(ET1m)[5] <- "ET"
ET1m$day <- gsub("Day.", "", ET1m$day)
ET1m$treatment <- gsub(" ", "", ET1m$treatment)
ET1m$exp <- 1
unique(ET1m$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15"ET2 <- read.csv(EVT_all[2])
ET2ET2m <- melt(ET2, id.vars = c("Pot.number", "Genotype", "Treatment"))
colnames(ET2m)[2] <- "genotype"
colnames(ET2m)[3] <- "treatment"
colnames(ET2m)[4] <- "day"
colnames(ET2m)[5] <- "ET"
ET2m$day <- gsub("Day", "", ET2m$day)
ET2m$treatment <- gsub(" ", "", ET2m$treatment)
ET2m$exp <- 2
ET2munique(ET2m$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"ET3 <- read.csv(EVT_all[3])
ET3ET3m <- melt(ET3, id.vars = c("Pot.number", "Genotype", "Treatment"))
colnames(ET3m)[2] <- "genotype"
colnames(ET3m)[3] <- "treatment"
colnames(ET3m)[4] <- "day"
colnames(ET3m)[5] <- "ET"
ET3m$day <- gsub("Day", "", ET3m$day)
ET3m$treatment <- gsub(" ", "", ET3m$treatment)
ET3m$exp <- 3
ET3munique(ET3m$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"ET4 <- read.csv(EVT_all[4])
ET4ET4m <- melt(ET4, id.vars = c("Pot.number", "Genotype", "Treatment"))
colnames(ET4m)[2] <- "genotype"
colnames(ET4m)[3] <- "treatment"
colnames(ET4m)[4] <- "day"
colnames(ET4m)[5] <- "ET"
ET4m$day <- gsub("X", "", ET4m$day)
ET4m$treatment <- gsub(" ", "", ET4m$treatment)
ET4m$exp <- 4
ET4munique(ET4m$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"ET5 <- read.csv(EVT_all[5])
ET5ET5m <- melt(ET5, id.vars = c("Pot.number", "Genotype", "Treatment"))
colnames(ET5m)[2] <- "genotype"
colnames(ET5m)[3] <- "treatment"
colnames(ET5m)[4] <- "day"
colnames(ET5m)[5] <- "ET"
ET5m$day <- gsub("X", "", ET5m$day)
ET5m$treatment <- gsub(" ", "", ET5m$treatment)
ET5m$exp <- 5
ET5munique(ET5m$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14"ET6 <- read.csv(EVT_all[6])
ET6ET6m <- melt(ET1, id.vars = c("Pot.number", "Genotype", "Treatment"))
colnames(ET6m)[2] <- "genotype"
colnames(ET6m)[3] <- "treatment"
colnames(ET6m)[4] <- "day"
colnames(ET6m)[5] <- "ET"
ET6m$day <- gsub("Day.", "", ET6m$day)
ET6m$treatment <- gsub(" ", "", ET6m$treatment)
ET6m$exp <- 6
ET6munique(ET6m$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15"OK - now let’s combine all the data into one data frame
ET <- rbind(ET1m, ET2m)
ET <- rbind(ET, ET3m)
ET <- rbind(ET, ET4m)
ET <- rbind(ET, ET5m)
ET <- rbind(ET, ET6m)
unique(ET$exp)## [1] 1 2 3 4 5 6unique(ET$day)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12" "13" "14" "15"unique(ET$treatment)## [1] "Control" "Drought"Now - let’s explore the evapotranspiration throughout time and when does the difference between the conditions starts to become significant:
ET$ET <- as.numeric(as.character(ET$ET))
ET$ID <- paste(ET$exp, ET$Pot.number, ET$genotype, ET$treatment, sep="_")
colnames(ET)[5] <- "EVT"
ET_clean <- subset(ET, ET$EVT > 0)
ET_clean <- subset(ET_clean, ET_clean$EVT < 200)
ET_clean$day <- as.numeric(ET_clean$day)
ET_clean <- subset(ET_clean, ET_clean$day < 15)
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))
ET_lgraph <- ET_lgraph + geom_line(alpha = 0.1)
ET_lgraph <- ET_lgraph + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= treatment), alpha=0.3)
ET_lgraph <- ET_lgraph + stat_summary(fun=mean, aes(group= treatment), size=0.7, geom="line", linetype = "dashed")
ET_lgraph <- ET_lgraph + stat_compare_means(aes(group = treatment), label = "p.signif", method = "t.test", hide.ns = T)
ET_lgraph <- ET_lgraph + ylab("Evapotranspiration (g H2O / day)") + xlab("Days After Stress") + scale_color_jco()
ET_lgraph
OK - since the data looks noisy again - let’s apply smooth splines here
too to clean it up just a little bit.
Smooth splines for Evapotranspiration
days <- 1:14
temp <- subset(ET_clean, ET_clean$ID == unique(ET_clean$ID)[1])
temp$day <- as.numeric(as.character(temp$day))
temp$day## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14temp$EVT## [1] 9.35 50.73 58.30 49.85 83.23 57.07 82.14 74.79 80.94 86.51
## [11] 94.90 99.90 97.60 116.39plot.spl <- with(temp, smooth.spline(day, EVT, df = 5))
plot(EVT ~ day, data = temp)
lines(plot.spl, col = "blue")
lines(predict(plot.spl, days), col="red")
Now - let’s create a dataframe to keep all of the spline-fitted
data:
names <- c(text="ID", "genotype", "treatment", "day", "EVT")
EVT_spline_data <- data.frame()
for (k in names) EVT_spline_data[[k]] <- as.character()
EVT_spline_dataAnd add the predicted values into the dataframe:
pred_temp <- predict(plot.spl, days)
pred_temp## $x
## [1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14
##
## $y
## [1] 20.78686 38.63288 52.19954 61.33375 67.77721 71.41976 74.88252
## [8] 78.15368 82.08269 86.96162 92.44815 98.16196 104.56285 112.29653# make empty dataframe
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$genotype[1]
EVT_spline_data[1:14,3] <- temp$treatment[1]
EVT_spline_datafinal_EVT_spline <- EVT_spline_data
all_plants <- unique(ET_clean$ID)
length(all_plants)## [1] 850Now - let’s loop it:
for(i in 2:850){
temp <- subset(ET_clean, ET_clean$ID == all_plants[i])
if(dim(temp)[1]>4){
plot.spl <- with(temp, smooth.spline(day, EVT, df = 5))
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$genotype[1]
EVT_spline_data[1:14,3] <- temp$treatment[1]
final_EVT_spline <- rbind(final_EVT_spline, EVT_spline_data)
}
else{
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$genotype[1]
spline_data[1:14,3] <- temp$treatment[1]
}
}
final_EVT_splineclean_EVT_spline <- subset(final_EVT_spline, final_EVT_spline$EVT > 0)
clean_EVT_splineNow - let’s replot the graph and see if it looks better:
clean_EVT_spline$day <- as.factor(as.numeric(clean_EVT_spline$day))
clean_EVT_spline$EVT <- as.numeric(as.character(clean_EVT_spline$EVT))
# let's clean up all of the odd plants that we also excluded from SV
clean_EVT_spline <- subset(clean_EVT_spline, clean_EVT_spline$ID %in% growth_data_clean$ID)
ET_lgraph_spline <- ggplot(data=clean_EVT_spline, aes(x= day, y=EVT, group = ID, color = treatment))
ET_lgraph_spline <- ET_lgraph_spline + geom_line(alpha = 0.1)
ET_lgraph_spline <- ET_lgraph_spline + stat_summary(fun.data = mean_se, geom="ribbon", linetype=0, aes(group= treatment), alpha=0.3)
ET_lgraph_spline <- ET_lgraph_spline + stat_summary(fun=mean, aes(group= treatment), size=0.7, geom="line", linetype = "dashed")
ET_lgraph_spline <- ET_lgraph_spline + stat_compare_means(aes(group = treatment), label = "p.signif", method = "t.test", hide.ns = T)
ET_lgraph_spline <- ET_lgraph_spline + ylab("Evapotranspiration (g H2O / day)") + xlab("Days After Stress") + scale_color_jco()
ET_lgraph_spline
Now - from this data let’s calculate some summary statistics. From the spline graph - it seems like mot pots are ar steady evapo-transpiration rate between 3 and 11th days after stress. Let’s calculate the median transpiration per pot, as well as cumulative transpiration per entire observation period.
sum_ET <- summaryBy(EVT ~ ID + treatment + genotype, data = clean_EVT_spline, FUN = function(x) c(sum=sum(x)))
sum_ETclean_EVT_spline$day <- as.numeric(as.character(clean_EVT_spline$day))
sub_EVT <- subset(clean_EVT_spline, clean_EVT_spline$day > 3)
sub_EVT <- subset(sub_EVT, sub_EVT$day < 12)
median_ET <- summaryBy(EVT ~ ID + treatment + genotype, data = sub_EVT, FUN = function(x) c(med=median(x)))
median_ETEVT_final <- merge(sum_ET, median_ET, by = c("ID", "treatment", "genotype"))
EVT_finalFinally - let’s plot this final data and see if we can see clear differences between control and drought stress:
EVT_med_overall <- ggerrorplot(EVT_final, y = "EVT.med", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="Median evapotranspiration (g H2O / day)")
EVT_med_overall <- EVT_med_overall + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
EVT_med_overall <- EVT_med_overall + scale_color_jco()
EVT_med_overall
EVT_sum_overall <- ggerrorplot(EVT_final, y = "EVT.sum", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="Cummulative evapotranspiration (g H2O / experiment)")
EVT_sum_overall <- EVT_sum_overall + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
EVT_sum_overall <- EVT_sum_overall + scale_color_jco()
EVT_sum_overall
Fusing all the data together
First - let’s wrangle the side-view data into the format that it can be merged with single time point data
final_splineSV_cast <- dcast(final_spline, ID ~ day)## Using Area.SUM as value column: use value.var to override.colnames(SV_cast) <- paste("Area.SUM_day", colnames(SV_cast), sep="")
colnames(SV_cast)[1] <- "ID"
SV_castLooks pretty good - now - let’s fuse it together with curated growth rate and FW, DW and WC data into one place:
curated_data <- merge(growth_data_clean, clean_FW, by=c("ID", "genotype", "treatment"))
curated_data <- merge(curated_data, SV_cast, by=c("ID"))
curated_dataLet’s merge to this PhotoSynQ data:
colnames(PSQ)[1] <- "Pot.number"
curated_data <- merge(curated_data, PSQ, by=c("Pot.number", "exp"))
curated_dataLast but not least - let’ also merge the evapo-transpiration data. For every day data - we need to do a similar trick as for the SV data:
clean_EVT_splineEVT_cast <- dcast(clean_EVT_spline, ID + treatment + genotype ~ day)## Using EVT as value column: use value.var to override.EVT_castcolnames(EVT_cast) <- paste("EVT_day", colnames(EVT_cast), sep="")
colnames(EVT_cast)[1] <- "ID"
colnames(EVT_cast)[2] <- "treatment"
colnames(EVT_cast)[3] <- "genotype"
EVT_castcomplete_EVT <- merge(EVT_cast, EVT_final, by=c("ID", "treatment", "genotype"))
complete_EVTcurated_data <- merge(curated_data, complete_EVT, by = c("ID", "treatment", "genotype"))
colnames(curated_data)## [1] "ID" "treatment"
## [3] "genotype" "Pot.number"
## [5] "exp" "GR"
## [7] "R2" "day"
## [9] "Area.SUM" "FW"
## [11] "DW" "WC"
## [13] "Area.SUM_day0" "Area.SUM_day1"
## [15] "Area.SUM_day2" "Area.SUM_day3"
## [17] "Area.SUM_day4" "Area.SUM_day5"
## [19] "Area.SUM_day6" "Area.SUM_day7"
## [21] "Area.SUM_day8" "Area.SUM_day9"
## [23] "Area.SUM_day10" "Area.SUM_day11"
## [25] "Area.SUM_day12" "Area.SUM_day13"
## [27] "Area.SUM_day14" "day13_FmPrime"
## [29] "day13_FoPrime" "day13_FvP_over_FmP"
## [31] "day13_Leaf.Temperature" "day13_leaf_thickness"
## [33] "day13_NPQt" "day13_PhiNPQ"
## [35] "day13_PS1.Active.Centers" "day13_PS1.Open.Centers"
## [37] "day13_SPAD" "day6_FmPrime"
## [39] "day6_FoPrime" "day6_FvP_over_FmP"
## [41] "day6_Leaf.Temperature" "day6_leaf_thickness"
## [43] "day6_NPQt" "day6_PhiNPQ"
## [45] "day6_PS1.Active.Centers" "day6_PS1.Open.Centers"
## [47] "day6_SPAD" "EVT_day1"
## [49] "EVT_day2" "EVT_day3"
## [51] "EVT_day4" "EVT_day5"
## [53] "EVT_day6" "EVT_day7"
## [55] "EVT_day8" "EVT_day9"
## [57] "EVT_day10" "EVT_day11"
## [59] "EVT_day12" "EVT_day13"
## [61] "EVT_day14" "EVT.sum"
## [63] "EVT.med"curated_data <- curated_data[,c(1:6, 10:63)]
colnames(curated_data) <- gsub(".x", "", colnames(curated_data))
curated_datawrite.csv(curated_data, "Cowpea_curated_data_all.csv", row.names = FALSE)Before we can go on - I noticed that some columns containing trait values are not numeric. So let’s rectify that:
sapply(curated_data, class)## ID treatment genotype
## "character" "character" "character"
## Pot.number p GR
## "integer" "numeric" "numeric"
## FW DW WC
## "numeric" "numeric" "numeric"
## Area.SUM_day0 Area.SUM_day1 Area.SUM_day2
## "numeric" "numeric" "numeric"
## Area.SUM_day3 Area.SUM_day4 Area.SUM_day5
## "numeric" "numeric" "numeric"
## Area.SUM_day6 Area.SUM_day7 Area.SUM_day8
## "numeric" "numeric" "numeric"
## Area.SUM_day9 Area.SUM_day10 Area.SUM_day11
## "numeric" "numeric" "numeric"
## Area.SUM_day12 Area.SUM_day13 Area.SUM_day14
## "numeric" "numeric" "numeric"
## day13_FmPrime day13_FoPrime day13_FvP_over_FmP
## "character" "character" "character"
## day13_Leaf.Temperature day13_leaf_thickness day13_NPQt
## "character" "character" "character"
## day13_PhiNPQ day13_PS1.Active.Centers day13_PS1.Open.Centers
## "character" "character" "character"
## day13_SPAD day6_FmPrime day6_FoPrime
## "character" "character" "character"
## day6_FvP_over_FmP day6_Leaf.Temperature day6_leaf_thickness
## "character" "character" "character"
## day6_NPQt day6_PhiNPQ day6_PS1.Active.Centers
## "character" "character" "character"
## day6_PS1.Open.Centers day6_SPAD EVT_day1
## "character" "character" "numeric"
## EVT_day2 EVT_day3 EVT_day4
## "numeric" "numeric" "numeric"
## EVT_day5 EVT_day6 EVT_day7
## "numeric" "numeric" "numeric"
## EVT_day8 EVT_day9 EVT_day10
## "numeric" "numeric" "numeric"
## EVT_day11 EVT_day12 EVT_day13
## "numeric" "numeric" "numeric"
## EVT_day14 EVT.sum EVT.med
## "numeric" "numeric" "numeric"cols.num <- c("GR", "FW", "DW", "WC", "Area.SUM_day0", "Area.SUM_day1", "Area.SUM_day2", "Area.SUM_day3", "Area.SUM_day4", "Area.SUM_day5", "Area.SUM_day6", "Area.SUM_day7", "Area.SUM_day8", "Area.SUM_day9", "Area.SUM_day10", "Area.SUM_day11", "Area.SUM_day12", "Area.SUM_day13", "Area.SUM_day14", "day13_FmPrime", "day13_FoPrime", "day13_FvP_over_FmP", "day13_Leaf.Temperature", "day13_leaf_thickness", "day13_NPQt", "day13_PhiNPQ", "day13_PS1.Active.Centers", "day13_PS1.Open.Centers", "day13_SPAD", "day6_FmPrime", "day6_FoPrime", "day6_FvP_over_FmP", "day6_Leaf.Temperature", "day6_leaf_thickness", "day6_NPQt", "day6_PhiNPQ", "day6_PS1.Active.Centers", "day6_PS1.Open.Centers", "day6_SPAD", "EVT_day1", "EVT_day2", "EVT_day3", "EVT_day4", "EVT_day5", "EVT_day6", "EVT_day7", "EVT_day8", "EVT_day9", "EVT_day10", "EVT_day11", "EVT_day12", "EVT_day13", "EVT_day14", "EVT.sum", "EVT.med")
curated_data[cols.num] <- sapply(curated_data[cols.num],as.numeric)
sapply(curated_data, class)## ID treatment genotype
## "character" "character" "character"
## Pot.number p GR
## "integer" "numeric" "numeric"
## FW DW WC
## "numeric" "numeric" "numeric"
## Area.SUM_day0 Area.SUM_day1 Area.SUM_day2
## "numeric" "numeric" "numeric"
## Area.SUM_day3 Area.SUM_day4 Area.SUM_day5
## "numeric" "numeric" "numeric"
## Area.SUM_day6 Area.SUM_day7 Area.SUM_day8
## "numeric" "numeric" "numeric"
## Area.SUM_day9 Area.SUM_day10 Area.SUM_day11
## "numeric" "numeric" "numeric"
## Area.SUM_day12 Area.SUM_day13 Area.SUM_day14
## "numeric" "numeric" "numeric"
## day13_FmPrime day13_FoPrime day13_FvP_over_FmP
## "numeric" "numeric" "numeric"
## day13_Leaf.Temperature day13_leaf_thickness day13_NPQt
## "numeric" "numeric" "numeric"
## day13_PhiNPQ day13_PS1.Active.Centers day13_PS1.Open.Centers
## "numeric" "numeric" "numeric"
## day13_SPAD day6_FmPrime day6_FoPrime
## "numeric" "numeric" "numeric"
## day6_FvP_over_FmP day6_Leaf.Temperature day6_leaf_thickness
## "numeric" "numeric" "numeric"
## day6_NPQt day6_PhiNPQ day6_PS1.Active.Centers
## "numeric" "numeric" "numeric"
## day6_PS1.Open.Centers day6_SPAD EVT_day1
## "numeric" "numeric" "numeric"
## EVT_day2 EVT_day3 EVT_day4
## "numeric" "numeric" "numeric"
## EVT_day5 EVT_day6 EVT_day7
## "numeric" "numeric" "numeric"
## EVT_day8 EVT_day9 EVT_day10
## "numeric" "numeric" "numeric"
## EVT_day11 EVT_day12 EVT_day13
## "numeric" "numeric" "numeric"
## EVT_day14 EVT.sum EVT.med
## "numeric" "numeric" "numeric"corrections to curated data that became apparent from the PSQ graphs
curated_data <- subset(curated_data, curated_data$day6_PS1.Active.Centers < 30)
curated_data <- subset(curated_data, curated_data$day6_PS1.Open.Centers > -100)
curated_data <- subset(curated_data, curated_data$day13_PS1.Open.Centers < 100)Because we didn’t have decoded values for all PhotosynQ data before - now after decoding we should provide these graphs and how they change +/- Drought and 6 vs 13 days after stress:
colnames(curated_data)## [1] "ID" "treatment"
## [3] "genotype" "Pot.number"
## [5] "p" "GR"
## [7] "FW" "DW"
## [9] "WC" "Area.SUM_day0"
## [11] "Area.SUM_day1" "Area.SUM_day2"
## [13] "Area.SUM_day3" "Area.SUM_day4"
## [15] "Area.SUM_day5" "Area.SUM_day6"
## [17] "Area.SUM_day7" "Area.SUM_day8"
## [19] "Area.SUM_day9" "Area.SUM_day10"
## [21] "Area.SUM_day11" "Area.SUM_day12"
## [23] "Area.SUM_day13" "Area.SUM_day14"
## [25] "day13_FmPrime" "day13_FoPrime"
## [27] "day13_FvP_over_FmP" "day13_Leaf.Temperature"
## [29] "day13_leaf_thickness" "day13_NPQt"
## [31] "day13_PhiNPQ" "day13_PS1.Active.Centers"
## [33] "day13_PS1.Open.Centers" "day13_SPAD"
## [35] "day6_FmPrime" "day6_FoPrime"
## [37] "day6_FvP_over_FmP" "day6_Leaf.Temperature"
## [39] "day6_leaf_thickness" "day6_NPQt"
## [41] "day6_PhiNPQ" "day6_PS1.Active.Centers"
## [43] "day6_PS1.Open.Centers" "day6_SPAD"
## [45] "EVT_day1" "EVT_day2"
## [47] "EVT_day3" "EVT_day4"
## [49] "EVT_day5" "EVT_day6"
## [51] "EVT_day7" "EVT_day8"
## [53] "EVT_day9" "EVT_day10"
## [55] "EVT_day11" "EVT_day12"
## [57] "EVT_day13" "EVT_day14"
## [59] "EVT.sum" "EVT.med"PSQ_decoded <- melt(curated_data[,c(1:3, 25:44)], id.vars = c("ID", "treatment", "genotype"))
PSQ_decoded$variable <- as.character(PSQ_decoded$variable)
for(i in 1:nrow(PSQ_decoded)){
PSQ_decoded$day[i] <- strsplit(PSQ_decoded$variable[i], "_")[[1]][1]
PSQ_decoded$trait[i] <- strsplit(PSQ_decoded$variable[i], "_")[[1]][2]
}
PSQ_decoded <- PSQ_decoded[,c(1:3,6:7,5)]
PSQ_decodedcPSQ <- dcast(PSQ_decoded, genotype + treatment + ID + day ~ trait)
cPSQcPSQ$day <- factor(cPSQ$day, levels = c("day6", "day13"))Fm_plot <- ggerrorplot(cPSQ, y = "FmPrime", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="Fm'")
Fm_plot <- Fm_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
Fm_plot <- Fm_plot + scale_color_jco()
Fm_plot
Fo_plot <- ggerrorplot(cPSQ, y = "FoPrime", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="Fo'")
Fo_plot <- Fo_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
Fo_plot <- Fo_plot + scale_color_jco()
Fo_plot
QY_plot <- ggerrorplot(cPSQ, y = "FvP", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="Quantum Yield (Fv' / Fm')")
QY_plot <- QY_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
QY_plot <- QY_plot + scale_color_jco()
QY_plot
leaf_thick_plot <- ggerrorplot(cPSQ, y = "leaf", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="leaf thicknes")
leaf_thick_plot <- leaf_thick_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
leaf_thick_plot <- leaf_thick_plot + scale_color_jco()
leaf_thick_plot
SPAD_plot <- ggerrorplot(cPSQ, y = "SPAD", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="SPAD")
SPAD_plot <- SPAD_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
SPAD_plot <- SPAD_plot + scale_color_jco()
SPAD_plot
PS1.Open.Centers_plot <- ggerrorplot(cPSQ, y = "PS1.Open.Centers", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="PS1.Open.Centers")
PS1.Open.Centers_plot <- PS1.Open.Centers_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
PS1.Open.Centers_plot <- PS1.Open.Centers_plot + scale_color_jco()
PS1.Open.Centers_plot
PS1.Active.Centers_plot <- ggerrorplot(cPSQ, y = "PS1.Active.Centers", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="PS1.Active.Centers")
PS1.Active.Centers_plot <- PS1.Active.Centers_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
PS1.Active.Centers_plot <- PS1.Active.Centers_plot + scale_color_jco()
PS1.Active.Centers_plot
NPQt_plot <- ggerrorplot(cPSQ, y = "NPQt", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="NPQt")
NPQt_plot <- NPQt_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
NPQt_plot <- NPQt_plot + scale_color_jco()
NPQt_plot
PhiNPQ_plot <- ggerrorplot(cPSQ, y = "PhiNPQ", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="PhiNPQ")
PhiNPQ_plot <- PhiNPQ_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
PhiNPQ_plot <- PhiNPQ_plot + scale_color_jco()
PhiNPQ_plot
leaf_temp_plot <- ggerrorplot(cPSQ, y = "Leaf.Temperature", x = "treatment", fill="treatment", color="treatment",
desc_stat = "mean_sd", add = "jitter", facet.by = "day",
add.params = list(color = "darkgray"),
xlab="Treatment", ylab="leaf temperature (C)")
leaf_temp_plot <- leaf_temp_plot + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Control",
label = "p.signif", hide.ns = T)
leaf_temp_plot <- leaf_temp_plot + scale_color_jco()
leaf_temp_plot
Extracting genotype specific values
Before moving on - let’s do a summary stats on all scored traits - and combine all of the replicates / remove the ones that were supposed to be removed due to problems.
library(dplyr)##
## Attaching package: 'dplyr'## The following object is masked from 'package:doBy':
##
## order_by## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unioncurated_data %>% group_by(genotype, treatment) %>% tally()From the above it seems that we only have replicates for the core accession lines that were included within each experiment. So we can go ahead and calculate the mean value per accession + treatment.
curated_sum <- summaryBy(GR + FW + DW + WC + Area.SUM_day0 + Area.SUM_day1 + Area.SUM_day2 + Area.SUM_day3 + Area.SUM_day4 + Area.SUM_day5 + Area.SUM_day6 + Area.SUM_day7 + Area.SUM_day8 + Area.SUM_day9 + Area.SUM_day10 + Area.SUM_day11 + Area.SUM_day12 + Area.SUM_day13 + Area.SUM_day14 + day6_FmPrime + day6_FoPrime + day6_FvP_over_FmP + day6_Leaf.Temperature + day6_leaf_thickness + day6_NPQt + day6_PhiNPQ + day6_PS1.Active.Centers + day6_PS1.Open.Centers + day6_SPAD + day13_FmPrime + day13_FoPrime + day13_FvP_over_FmP + day13_Leaf.Temperature + day13_leaf_thickness + day13_NPQt + day13_PhiNPQ + day13_PS1.Active.Centers + day13_PS1.Open.Centers + day13_SPAD + EVT_day1 + EVT_day2 + EVT_day3 + EVT_day4 + EVT_day5 + EVT_day6 + EVT_day7 + EVT_day8 + EVT_day9 + EVT_day10 + EVT_day11 + EVT_day12 + EVT_day13 + EVT_day14 + EVT.sum + EVT.med ~ treatment + genotype, curated_data)
curated_sumNow - let’s split the data for each treatment:
# just add this one trait - that might be interesting to look at:
curated_sum$EVT.p.FW.mean <- curated_sum$EVT.sum.mean / curated_sum$FW.mean
# prep data
Control_sum <- subset(curated_sum, curated_sum$treatment == "Control")
Drought_sum <- subset(curated_sum, curated_sum$treatment == "Drought")
dim(Control_sum)## [1] 321 58dim(Drought_sum)## [1] 274 58colnames(Control_sum)<- gsub(".mean", ".Cntrl", colnames(Control_sum))
colnames(Drought_sum)<- gsub(".mean", ".Drght", colnames(Drought_sum))
Control_sum <- Control_sum[,c(2:58)]
Drought_sum <- Drought_sum[,c(2:58)]
all_sum <- merge(Control_sum, Drought_sum, by = "genotype", all=TRUE)
head(all_sum)dim(all_sum)## [1] 322 113this data isn now ready for GWAS
write.csv(all_sum, "Cowpea2022.GWASready.csv", row.names = FALSE)Correlations between the traits:
Now we have everything combined - let’s have a look at the correlations between the traits. But let’s examine the correlations for Control and Drought separately
#install.packages("ggcorrplot")
library(ggcorrplot)Drght <- Drought_sum[,2:57]
Drght <- na.omit(Drght)
corr2 <- cor(Drght)
head(corr2[, 1:6])## GR.Drght FW.Drght DW.Drght WC.Drght
## GR.Drght 1.00000000 0.4722305 0.1697122 0.5096016
## FW.Drght 0.47223054 1.0000000 0.7818121 0.9916218
## DW.Drght 0.16971224 0.7818121 1.0000000 0.6947193
## WC.Drght 0.50960163 0.9916218 0.6947193 1.0000000
## Area.SUM_day0.Drght -0.04228666 0.5458052 0.5860622 0.5082203
## Area.SUM_day1.Drght 0.04572535 0.6214716 0.6239648 0.5876561
## Area.SUM_day0.Drght Area.SUM_day1.Drght
## GR.Drght -0.04228666 0.04572535
## FW.Drght 0.54580520 0.62147162
## DW.Drght 0.58606220 0.62396478
## WC.Drght 0.50822034 0.58765613
## Area.SUM_day0.Drght 1.00000000 0.97710932
## Area.SUM_day1.Drght 0.97710932 1.00000000p.mat2 <- cor_pmat(Drght)
head(p.mat2[,1:5])## GR.Drght FW.Drght DW.Drght WC.Drght
## GR.Drght 0.000000e+00 1.264371e-16 4.849614e-03 1.631416e-19
## FW.Drght 1.264371e-16 0.000000e+00 9.735001e-58 8.484803e-244
## DW.Drght 4.849614e-03 9.735001e-58 0.000000e+00 8.263522e-41
## WC.Drght 1.631416e-19 8.484803e-244 8.263522e-41 0.000000e+00
## Area.SUM_day0.Drght 4.857524e-01 1.131204e-22 1.148297e-26 2.117259e-19
## Area.SUM_day1.Drght 4.509553e-01 1.141334e-30 5.703505e-31 7.768780e-27
## Area.SUM_day0.Drght
## GR.Drght 4.857524e-01
## FW.Drght 1.131204e-22
## DW.Drght 1.148297e-26
## WC.Drght 2.117259e-19
## Area.SUM_day0.Drght 0.000000e+00
## Area.SUM_day1.Drght 7.378780e-185Correlation_plot_Drought <- ggcorrplot(corr2, p.mat = p.mat2, type = "upper", outline.col = "white",
colors = c("#6D9EC1", "white", "#E46726"),
ggtheme = ggplot2::theme_gray, insig = "blank")
Correlation_plot_Drought
colnames(p.mat2) %in% colnames(corr2)## [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 TRUErownames(p.mat2) %in% rownames(corr2)## [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 TRUEpdf("Fig.Cowpea2022.Correlation_plot_Drought_all_traits.pdf", width = 13, height = 13)
plot(Correlation_plot_Drought)
dev.off()## quartz_off_screen
## 2Cntrl <- Control_sum[,2:57]
corr <- cor(Cntrl)
head(corr[, 1:6])## GR.Cntrl FW.Cntrl DW.Cntrl WC.Cntrl Area.SUM_day0.Cntrl
## GR.Cntrl 1.0000000 0.7077127 0.4728750 0.7269195 0.1716607
## FW.Cntrl 0.7077127 1.0000000 0.8313913 0.9933683 0.5415490
## DW.Cntrl 0.4728750 0.8313913 1.0000000 0.7619872 0.6838490
## WC.Cntrl 0.7269195 0.9933683 0.7619872 1.0000000 0.4896220
## Area.SUM_day0.Cntrl 0.1716607 0.5415490 0.6838490 0.4896220 1.0000000
## Area.SUM_day1.Cntrl 0.2872299 0.6431997 0.7494108 0.5945194 0.9729493
## Area.SUM_day1.Cntrl
## GR.Cntrl 0.2872299
## FW.Cntrl 0.6431997
## DW.Cntrl 0.7494108
## WC.Cntrl 0.5945194
## Area.SUM_day0.Cntrl 0.9729493
## Area.SUM_day1.Cntrl 1.0000000p.mat <- cor_pmat(Cntrl)
head(p.mat[, 1:4])## GR.Cntrl FW.Cntrl DW.Cntrl WC.Cntrl
## GR.Cntrl 0.000000e+00 4.627625e-50 2.741959e-19 5.251592e-54
## FW.Cntrl 4.627625e-50 0.000000e+00 2.139113e-83 9.681814e-302
## DW.Cntrl 2.741959e-19 2.139113e-83 0.000000e+00 3.732208e-62
## WC.Cntrl 5.251592e-54 9.681814e-302 3.732208e-62 0.000000e+00
## Area.SUM_day0.Cntrl 2.025285e-03 7.383867e-26 1.385966e-45 9.333658e-21
## Area.SUM_day1.Cntrl 1.633558e-07 7.143254e-39 4.463349e-59 4.610805e-32Correlation_plot_Control <- ggcorrplot(corr, type = "lower", outline.col = "white", colors = c("#6D9EC1", "white", "#E46726"),
p.mat = p.mat, ggtheme = ggplot2::theme_gray, insig = "blank")
Correlation_plot_Control
pdf("Fig.Cowpea2022.Correlation_plot_Control_all_traits.pdf", width = 13, height = 13)
plot(Correlation_plot_Control)
dev.off()## quartz_off_screen
## 2Extracting figure files
finally - let’s extract other most important figures from this analysis so far:
library(cowplot)##
## Attaching package: 'cowplot'## The following object is masked from 'package:ggpubr':
##
## get_legendtime_plot <- plot_grid(Area_lgraph, ET_lgraph_spline, ncol = 2, labels = c("A", "B"))
time_plot
PSQ_plot <- plot_grid(QY_plot, NPQt_plot, leaf_temp_plot, SPAD_plot, ncol=4, labels = c("C", "D", "E", "F"))
PSQ_plot
final_plot <- plot_grid(time_plot, PSQ_plot, ncol = 1, labels="")
final_plot
pdf("Fig.Cowpea2022.Growth.EVT.PSQ.pdf", width = 15, height = 10)
plot(final_plot)
dev.off()## quartz_off_screen
## 2pdf("Fig.Cowpea2022.Area.vs.weight.pdf", width = 22.5, height = 7.5)
plot_grid(FW_Area, DW_Area, WC_Area, ncol = 3)
dev.off()## quartz_off_screen
## 2