Heat experiment for TE mutants

Experiment setup

The seeds of mutant lines and control lines were received from Dr. Etienne Bucher. The seeds were stratified at 4C for 24h and germinated from 1st of January on Cornell Mix soil. The plants were grown in BTI’s walk-in growth chamber #03, with the 16h light period (lights on from 7:00 AM to 11:00 PM), 22C throughout the growth period and 60% relative humidity.

On 20th of January (20 days after germination), trays K, N, R, S, T and V were moved at 9:00 AM into growth chamber #35 for heat treatment. In this growth chamber, the temperature gradually increased from 22C to 40C between 9:00 and 11:00 AM, and kept at 40C for 6 h (between 11:00 AM and 5:00 PM). At 5 PM the trays were moved into Growth chamber #03 where the did recover from the heat stress under the same conditions as they did germinate in.

The plant germination and heat stress treatment was performed by Magdalena Julkowska (PI, BTI)

The images were recorded using Raspberry Pi cameras, and processed using PlantCV pipeline by Liang Yu (PostDoc, BTI). The resulting data was processed in R - as described by the data analysis pipeline below by both Liang Yu and Magdalena Julkowska.

Data preparation

First of all - we need to load the libraries that we will use in this analysis:

library(tidyr)
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(corrplot)

## corrplot 0.84 loaded

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──

## ✓ ggplot2 3.3.1     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(ggpubr)

## Warning: package 'ggpubr' was built under R version 4.0.2

library(directlabels)

## Warning: package 'directlabels' was built under R version 4.0.2

library(ggplot2)
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

Then - let’s load all of the files and calculate the time points for each of the images:

# loading the data
sample.list <- list.files(pattern=".clean.csv")
sample.list

##  [1] "RaspiK_CamA.clean.csv" "RaspiK_CamB.clean.csv" "RaspiL_CamA.clean.csv"
##  [4] "RaspiL_CamB.clean.csv" "RaspiM_CamA.clean.csv" "RaspiM_CamB.clean.csv"
##  [7] "RaspiN_CamA.clean.csv" "RaspiN_CamB.clean.csv" "RaspiO_CamB.clean.csv"
## [10] "RaspiP_CamA.clean.csv" "RaspiP_CamB.clean.csv" "RaspiQ_CamA.clean.csv"
## [13] "RaspiQ_CamB.clean.csv" "RaspiS_CamA.clean.csv" "RaspiT_CamA.clean.csv"
## [16] "RaspiT_CamB.clean.csv" "RaspiU_CamA.clean.csv" "RaspiU_CamB.clean.csv"
## [19] "RaspiV_CamA.clean.csv" "RaspiV_CamB.clean.csv"

i=1
temp <- read.csv(sample.list[i])
raspi_combined <- temp
raspi_combined

for (i in 2:length(sample.list)){
  temp <- read.csv(sample.list[i])
  raspi_combined <- rbind(raspi_combined, temp)
}
raspi_combined

#Create Meta phenotype sheet by combining multiple files
raspi_clean <- raspi_combined[,c(2,3,16,17,18,19,20,21,22,23,24,25,26,27)]

raspi_clean

Map <- read.csv("/Users/magdalena/Downloads/Magda-Etienne/Heat_EXP_data/heat_exp_map.csv", header = T)

raspi_clean$Treatment = "treatment"
raspi_clean$Genotype = "genotype"


#replacing each genotype based on Slots
for (i in 1: nrow(Map)){
  if(nrow(raspi_clean[raspi_clean$RasPi.ID == Map[i,1] & raspi_clean$Cam.ID == Map[i,2] & raspi_clean$ID == Map[i,4],]) == 0) next
  raspi_clean[raspi_clean$RasPi.ID == Map[i,1] & raspi_clean$Cam.ID == Map[i,2] & raspi_clean$ID == Map[i,4],]$Treatment = Map[i,3]
  raspi_clean[raspi_clean$RasPi.ID == Map[i,1] & raspi_clean$Cam.ID == Map[i,2] & raspi_clean$ID == Map[i,4],]$Genotype = Map[i,6]
}

#Re-oder dataframe based on Time and genotype
raspi_clean <- raspi_clean[
  order(raspi_clean[,15], raspi_clean[,16], raspi_clean[,14]),
]

select <- read.csv("/Users/magdalena/Downloads/Magda-Etienne/Heat_EXP_data/heat_exp_select.csv", header = T)

raspi_clean$select = "NA"

for (i in 1:nrow(select)){
  raspi_clean[raspi_clean$RasPi.ID == select[i,1] & raspi_clean$Cam.ID == select[i,2] & raspi_clean$ID == select[i,3],]$select = "TRUE"
}

raspi_select <- raspi_clean[raspi_clean$select == "TRUE",]
raspi_select$ID <- as.character(raspi_select$ID)
raspi_select

Then - let’s have a look at the dataset created:

head(raspi_select)

raspi_select$Genotype <- gsub("Line_", "hcLine_", raspi_select$Genotype)
raspi_select$Genotype <- gsub("hcLine_52", "control_drugs", raspi_select$Genotype)
raspi_select$Genotype <- gsub("hcLine_50", "control_heat", raspi_select$Genotype)
head(raspi_select)

Now - let’s plot the data based on each rig:

rig_plot <- ggplot(data = raspi_select, aes(x= all.min, y=area, color = ID)) 
rig_plot <- rig_plot + geom_line(alpha = 0.1,size = 0.8, aes(group= ID)) 
rig_plot <- rig_plot + facet_wrap(~ RasPi.ID, ncol=5) 
rig_plot <- rig_plot + xlab("Total Minutes After Start") + rremove("legend")
rig_plot <- rig_plot +  ylab("Leaf area") + xlim(0,21870) + ylim(0,50000)
rig_plot

## Warning: Removed 3518 row(s) containing missing values (geom_path).

As you can see from the above - there is a large variation between the rigs in terms of growth.

Rigs K, N, R, S, T, V underwent heat stress - and for some of them you can see a drop in their Leaf area which reflects the hyponastic position of the leaves after the plants were moved back from the heat stress treatment. Let’s look at one rig for better insight:

RaspiS <- subset(raspi_select, raspi_select$RasPi.ID == "raspiS")
Srig_plot <- ggplot(data = RaspiS, aes(x= all.min, y=area, color = ID)) 
Srig_plot <- Srig_plot + geom_line(alpha = 0.1,size = 0.8, aes(group= ID)) 
Srig_plot <- Srig_plot + xlab("Total Minutes After Start") + rremove("legend")
Srig_plot <- Srig_plot +  ylab("Leaf area") + xlim(0,21870) + ylim(0,50000)
Srig_plot

## Warning: Removed 1989 row(s) containing missing values (geom_path).

ggplotly(Srig_plot)

## Warning: `group_by_()` is deprecated as of dplyr 0.7.0.
## Please use `group_by()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

Calculate the average and std for each certain plant at each time point for analysis

Control plants:

#Extract all control plants
control <- raspi_select[raspi_select$Treatment == "Control",]

control <- control %>% group_by(day,RasPi.ID,Cam.ID,ID) %>% mutate(Rep.number = n())
control <- control %>% group_by(day,RasPi.ID,Cam.ID,ID) %>% mutate(area.min = min(area))
control <- control %>% group_by(day,RasPi.ID,Cam.ID,ID) %>% mutate(area.max = max(area))
control$delta <- control$area.max-control$area.min
control_plot <- control[,c(2,3,4,6,7,8,9,10,14,15,16,19,20,21)]

#add a new index to represent each plant
control_plot$Index <- paste0(control_plot$RasPi.ID,"_",control_plot$Cam.ID,"_",control_plot$ID)
head(control_plot)

Plot the plant growth over-time for each genotype based on average and each rep

#Each genotype will be split
##Both average and exact rep values will be plotted

ggplot(data = control_plot, aes(x= all.min, y=area, color = ID)) + 
  geom_line(alpha = 0.3,size = 0.8, aes(group= Index)) + facet_wrap(~ Genotype, ncol=5) +
  geom_point(alpha = 0.6,size = 0.4, aes(group= Index)) + facet_wrap(~ Genotype, ncol=5) +
  xlab("Total Minutes After Start") +  ylab("Leaf area") +
  stat_summary(fun.y=mean, size=1, geom="line", linetype = "dashed", color = 'red') +
  xlim(0,21870) + ylim(0,50000) +
  geom_dl(aes(label = Index), method = list(dl.trans(x = x - 1), "last.points", cex = 0.4))

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 14990 rows containing non-finite values (stat_summary).

## Warning: Removed 14990 row(s) containing missing values (geom_path).

## Warning: Removed 14990 rows containing missing values (geom_point).

## Warning: Removed 14990 rows containing missing values (geom_dl).

Calculate the average and std for each genotype at each time point for analysis

#Extract all heat treated plants
Treatment <- raspi_select[raspi_select$Treatment == "Heat",]

Treatment <- Treatment %>% group_by(day,RasPi.ID,Cam.ID,ID) %>% mutate(Rep.number = n())
Treatment <- Treatment %>% group_by(day,RasPi.ID,Cam.ID,ID) %>% mutate(area.min = min(area))
Treatment <- Treatment %>% group_by(day,RasPi.ID,Cam.ID,ID) %>% mutate(area.max = max(area))
Treatment$delta <- Treatment$area.max-Treatment$area.min
Treatment_plot <- Treatment[,c(2,3,4,6,7,8,9,10,14,15,16,19,20,21)]

#add a new index to represent each plant
Treatment_plot$Index <- paste0(Treatment_plot$RasPi.ID,"_",Treatment_plot$Cam.ID,"_",Treatment_plot$ID)
head(Treatment_plot)

Plot the plant growth over-time for each genotype based on average and each rep

#Each genotype will be split
##Both average and exact rep values will be plotted
ggplot(data = Treatment_plot, aes(x= all.min, y=area, color = ID)) + 
  geom_line(alpha = 0.3,size = 0.8, aes(group= Index)) + facet_wrap(~ Genotype, ncol=5) +
  geom_point(alpha = 0.6,size = 0.4, aes(group= Index)) + facet_wrap(~ Genotype, ncol=5) +
  xlab("Total Minutes After Start") +  ylab("Leaf area") +
  stat_summary(fun.y=mean, size=1, geom="line", linetype = "dashed", color = 'red') +
  xlim(0,21870) + ylim(0,50000) +
  geom_dl(aes(label = Index), method = list(dl.trans(x = x - 1), "last.points", cex = 0.4))

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 14804 rows containing non-finite values (stat_summary).

## Warning: Removed 14803 row(s) containing missing values (geom_path).

## Warning: Removed 14804 rows containing missing values (geom_point).

## Warning: Removed 14804 rows containing missing values (geom_dl).

Combine control and treatment data sets:

Raspi_final <- rbind(Treatment_plot,control_plot)

Let’s collapse the data set for stats summary:

#Calculate the mean, standaard deviation for each genotype under certain condition per time stamp
Raspi_final <- Raspi_final %>% group_by(Treatment,all.min,Genotype) %>% mutate(Genotype.Rep.number = n())
Raspi_final <- Raspi_final %>% group_by(Treatment,all.min,Genotype) %>% mutate(Genotype.area.mean = ave(area))
Raspi_final <- Raspi_final %>% group_by(Treatment,all.min,Genotype) %>% mutate(Genotype.area.std = sd(area))

#Simplify the stats dataframe 
Raspi_stats <- unique(Raspi_final[,c(5,6,7,8,9,10,11,15,16,17,18)])

head(Raspi_stats)

Collapse the data set for stats summary

#Calculate the mean, standaard deviation for each genotype under certain condition per time stamp
Raspi_final <- Raspi_final %>% group_by(Treatment,all.min,Genotype) %>% mutate(Genotype.Rep.number = n())
Raspi_final <- Raspi_final %>% group_by(Treatment,all.min,Genotype) %>% mutate(Genotype.area.mean = ave(area))
Raspi_final <- Raspi_final %>% group_by(Treatment,all.min,Genotype) %>% mutate(Genotype.area.std = sd(area))

#Simplify the stats dataframe 
Raspi_stats <- unique(Raspi_final[,c(5,6,7,8,9,10,11,15,16,17,18)])

head(Raspi_stats)

Calculate growth rate

Calculate the growth rate over time based on individual rep under certain genotype

#Build the Table for multiple linear regression to be written
Growth_rate_rep <- data.frame("Month" = integer(0), 
                             "Date" = integer(0), "Treatment" = character(0),"Plant.ID" = character(0), "Genotype"= character(0), 
                             "intercept" = integer(0), "delta" = integer(0), "R.square" = integer(0))

#Build time schedule
sub_raspi_rep <- unique(Raspi_stats[,c(1,2,6,8,7)])


#Loop the calculation of linear fit of Day_time growth
for (x in 1:nrow(sub_raspi_rep)){
  #Name initial four columns
  Growth_rate_rep[x,1:5] <- sub_raspi_rep[x, 1:5]
  
  #Calculate the intercept
  Growth_rate_rep[x,6] <- as.numeric(lm(area ~ all.min, Raspi_final[Raspi_final$month== as.character(sub_raspi_rep[x,1]) & 
                                        Raspi_final$day == as.character(sub_raspi_rep[x,2]) & Raspi_final$Treatment == as.character(sub_raspi_rep[x,3]) &                                                Raspi_final$Genotype == as.character(sub_raspi_rep[x,5]) &
                                        Raspi_final$Index == as.character(sub_raspi_rep[x,4]),])$coefficient[[1]])
  
  #Calculate the delta
  Growth_rate_rep[x,7] <- as.numeric(lm(area ~ all.min, Raspi_final[Raspi_final$month== as.character(sub_raspi_rep[x,1]) & 
                                        Raspi_final$day == as.character(sub_raspi_rep[x,2]) & Raspi_final$Treatment == as.character(sub_raspi_rep[x,3]) &                                                Raspi_final$Genotype == as.character(sub_raspi_rep[x,5]) &
                                        Raspi_final$Index == as.character(sub_raspi_rep[x,4]),])$coefficient[[2]])
  
  #Calculate the R.square
  Growth_rate_rep[x,8] <- summary(lm(area ~ all.min, Raspi_final[Raspi_final$month== as.character(sub_raspi_rep[x,1]) & 
                                        Raspi_final$day == as.character(sub_raspi_rep[x,2]) & Raspi_final$Treatment == as.character(sub_raspi_rep[x,3]) &                                                Raspi_final$Genotype == as.character(sub_raspi_rep[x,5]) &
                                        Raspi_final$Index == as.character(sub_raspi_rep[x,4]),]))$r.squared
}

Clean and write the final table

Growth_rate_rep <- na.omit(Growth_rate_rep)
head(Growth_rate_rep)

write.csv(Growth_rate_rep, "Growth_rate_all_reps.csv", row.names = FALSE)

Plot distribution of Rsquare

ggplot(data = Growth_rate_rep, aes(x = R.square, fill = Treatment))+ 
  geom_histogram() + facet_wrap(~ Genotype, ncol=5) +
  ylab("Count") + xlab("R square of linear fit")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Make comparison for each LINE between treatment and control

ggplot(data = Growth_rate_rep, aes(x = Date, y =delta, color = Treatment))+ 
   
    #plot each rep
    geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("dodgerblue3", "tomato1")) +
    #plot average lines
    stat_summary(fun.y=mean,  size=0.7, geom="line", linetype = "dashed") +
    #Add unit to the Y-axis
    ylab("Growth rate") + xlab("Days after experiment") +
    stat_compare_means(aes(x = Date), label = "p.signif", method = "t.test", hide.ns = T) +
    
    #Delimit the time schedule
    xlim(12,27) + ylim(0, 20)

## Warning: `fun.y` is deprecated. Use `fun` instead.

## Warning: Removed 1205 rows containing non-finite values (stat_summary).

## Warning: Removed 1205 rows containing non-finite values (stat_compare_means).

## Warning: Removed 967 row(s) containing missing values (geom_path).

    #facet based on genotype
    facet_wrap(~ Genotype, ncol=5)

## <ggproto object: Class FacetWrap, Facet, gg>
##     compute_layout: function
##     draw_back: function
##     draw_front: function
##     draw_labels: function
##     draw_panels: function
##     finish_data: function
##     init_scales: function
##     map_data: function
##     params: list
##     setup_data: function
##     setup_params: function
##     shrink: TRUE
##     train_scales: function
##     vars: function
##     super:  <ggproto object: Class FacetWrap, Facet, gg>

However - we want to split the plants per treatment only AFTER 20th of January - because that’s when the treatment took place

Control <- subset(Growth_rate_rep, Growth_rate_rep$Date < 20)
Treatment <- subset(Growth_rate_rep, Growth_rate_rep$Date > 20)
Treatmen <- as.data.frame(Treatment)
head(Treatmen)

Control <- subset(Control, Control$delta > 0)
Treatment2 <- subset(Treatmen, Treatmen$delta > 0)

Control_growth <- ggplot(data = Control, aes(x = Date, y =delta, color = Genotype))
Control_growth <- Control_growth + geom_line(alpha = 0.3, aes(group= Plant.ID))
Control_growth <- Control_growth + stat_summary(fun.y=mean,  size=0.7, geom="line", linetype = "dashed")

## Warning: `fun.y` is deprecated. Use `fun` instead.

Control_growth <- Control_growth + ylab("Growth rate") + xlab("Days after germination")
Control_growth <- Control_growth + facet_wrap(~ Genotype, ncol=5) + rremove("legend")
Control_growth

Treatment2$DAS <- Treatment2$Date - 20

Heat_growth <- ggplot(data = Treatment2, aes(x = DAS, y =delta, color = Genotype))
Heat_growth <- Heat_growth + geom_line(alpha = 0.3, aes(group= Plant.ID))
Heat_growth <- Heat_growth + stat_summary(fun.y=mean,  size=0.7, geom="line", linetype = "dashed")

## Warning: `fun.y` is deprecated. Use `fun` instead.

Heat_growth <- Heat_growth + ylab("Growth rate") + xlab("Days after stress")
Heat_growth <- Heat_growth + facet_wrap(Genotype ~ Treatment) + rremove("legend")
Heat_growth

We do want to compare all of the lines to control_heat (Control for drugs) or control_drugs (Control for heat)

In this case we will have to subset into a series of subsets:

mutants <- c("hcLine_02", "hcLine_04", "hcLine_07", "hcLine_09", "hcLine_18", "hcLine_31", "hcLine_33A", "hcLine_33D", "hcLine_34", "hcLine_45A", "hcLine_45B", "hcLine_45C", "hcLine_45D")

getwd()

## [1] "/Users/magdalena/Downloads/Magda-Etienne/Heat_EXP_data"

for(i in 1:length(mutants)){
  want_now <- c("control_heat", mutants[i])
  want_now2 <- c("control_drugs", mutants[i])
  temp_50 <- subset(Control, Control$Genotype %in% want_now)
  temp_52 <- subset(Control, Control$Genotype %in% want_now2)
  
  Control_l50 <- ggplot(data = temp_50, aes(x = Date, y =delta, color = Genotype))
  Control_l50 <- Control_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Control_l50 <- Control_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Control_l50 <- Control_l50 + ylab("Growth rate") + xlab("Days after germination") + theme(legend.position = c(0.8, 0.8))
  Control_l50 <- Control_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_uncorrected_Control_", mutants[i],"_compared_to_control_heat.pdf", sep =""))
  plot(Control_l50)
  dev.off()
  
  Control_l52 <- ggplot(data = temp_52, aes(x = Date, y =delta, color = Genotype))
  Control_l52 <- Control_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Control_l52 <- Control_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Control_l52 <- Control_l52 + ylab("Growth rate") + xlab("Days after germination")+ theme(legend.position = c(0.8, 0.8))
  Control_l52 <- Control_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_uncorrected_Control_", mutants[i],"_compared_to_control_drugs.pdf", sep =""))
  plot(Control_l52)
  dev.off()
  
  temp_50 <- subset(Treatment2, Treatment2$Genotype %in% want_now)
  temp_52 <- subset(Treatment2, Treatment2$Genotype %in% want_now2)
  temp_50$DAS <- as.numeric(as.character(temp_50$DAS))
  temp_50 <- subset(temp_50, temp_50$DAS < 7)
  temp_52 <- subset(temp_52, temp_52$DAS < 7)
  
  Treatment_l50 <- ggplot(data = temp_50, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l50 <- Treatment_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l50 <- Treatment_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l50 <- Treatment_l50 + facet_wrap(~ Treatment, ncol=2)
  Treatment_l50 <- Treatment_l50 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l50 <- Treatment_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_uncorrected_Treatment_", mutants[i],"_compared_to_control_heat.pdf", sep ="")) 
  plot(Treatment_l50)
  dev.off()
  
  Treatment_l52 <- ggplot(data = temp_52, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l52 <- Treatment_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l52 <- Treatment_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l52 <- Treatment_l52 + facet_wrap(~ Treatment, ncol=2)
  Treatment_l52 <- Treatment_l52 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l52 <- Treatment_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_uncorrected_Treatment_", mutants[i],"_compared_to_control_drugs.pdf", sep =""))
  plot(Treatment_l52)
  dev.off()
}

OK - so it seems we do need to correct for individual trays that have been under different light conditions. Not exactly the perfect solution - but so be it

Normalizing plant size per rig

First - let’s calculate the average plant size PER RIG at 20th of January at 8:00 AM

one_day <- subset(raspi_select, raspi_select$day == 21)
one_day <- subset(one_day, one_day$hour == 8)
#one_day <- subset(one_day, one_day$min == 0)
unique(one_day$RasPi.ID)

##  [1] "raspiL" "raspiM" "raspiO" "raspiP" "raspiQ" "raspiU" "raspiK" "raspiN"
##  [9] "raspiS" "raspiT" "raspiV"

library(doBy)

## 
## Attaching package: 'doBy'

## The following object is masked from 'package:dplyr':
## 
##     order_by

head(one_day)

summed_up <- summaryBy(area ~ RasPi.ID, data = one_day)
summed_up

Now - let’s fuse our Summary into the Raspi_select file to normalize the area:

all_normalized <- merge(raspi_select, summed_up, by="RasPi.ID", all = T)
head(all_normalized)

dim(raspi_select)

## [1] 124921     17

dim(all_normalized)

## [1] 124921     18

all_normalized <- subset(all_normalized, all_normalized$day < 28)

all_normalized$area.norm <- all_normalized$area / all_normalized$area.mean
all_normalized$plantID <- paste(all_normalized$RasPi.ID, all_normalized$Cam.ID, all_normalized$ID, sep="_")

rig_plot <- ggplot(data = all_normalized, aes(x= all.min, y=area.norm, color = Treatment)) 
rig_plot <- rig_plot + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
rig_plot <- rig_plot + facet_wrap(~ Genotype, ncol=5) 
rig_plot <- rig_plot + xlab("Total Minutes After Imaging Start") + rremove("legend")
rig_plot <- rig_plot +  ylab("Normalized leaf area") 
rig_plot

Let’s inspect all of the genotypes for possible outliers and remove them from the overall dataset:

library(plotly)
Line02 <- subset(all_normalized, all_normalized$Genotype == "hcLine_02")
graph_Line02 <- ggplot(data = Line02, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line02 <- graph_Line02 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line02 <- graph_Line02 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line02 <- graph_Line02 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line02)

Line04 <- subset(all_normalized, all_normalized$Genotype == "hcLine_04")
graph_Line04 <- ggplot(data = Line04, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line04 <- graph_Line04 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line04 <- graph_Line04 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line04 <- graph_Line04 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line04)

Line07 <- subset(all_normalized, all_normalized$Genotype == "hcLine_07")
graph_Line07 <- ggplot(data = Line07, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line07 <- graph_Line07 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line07 <- graph_Line07 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line07 <- graph_Line07 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line07)

Line09 <- subset(all_normalized, all_normalized$Genotype == "hcLine_09")
graph_Line09 <- ggplot(data = Line09, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line09 <- graph_Line09 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line09 <- graph_Line09 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line09 <- graph_Line09 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line09)

Line18 <- subset(all_normalized, all_normalized$Genotype == "hcLine_18")
graph_Line18 <- ggplot(data = Line18, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line18 <- graph_Line18 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line18 <- graph_Line18 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line18 <- graph_Line18 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line18)

Line31 <- subset(all_normalized, all_normalized$Genotype == "hcLine_31")
graph_Line31 <- ggplot(data = Line31, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line31 <- graph_Line31 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line31 <- graph_Line31 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line31 <- graph_Line31 +  ylab("Normalized leaf area") + ylim(0,5)
ggplotly(graph_Line31)

Line33A <- subset(all_normalized, all_normalized$Genotype == "hcLine_33A")
graph_Line33A <- ggplot(data = Line33A, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line33A <- graph_Line33A + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line33A <- graph_Line33A + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line33A <- graph_Line33A +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line33A)

Line33D <- subset(all_normalized, all_normalized$Genotype == "hcLine_33D")
graph_Line33D <- ggplot(data = Line33D, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line33D <- graph_Line33D + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line33D <- graph_Line33D + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line33D <- graph_Line33D +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line33D)

Line34 <- subset(all_normalized, all_normalized$Genotype == "hcLine_34")
graph_Line34 <- ggplot(data = Line34, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line34 <- graph_Line34 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line34 <- graph_Line34 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line34 <- graph_Line34 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line34)

Line45A <- subset(all_normalized, all_normalized$Genotype == "hcLine_45A")
graph_Line45A <- ggplot(data = Line45A, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line45A <- graph_Line45A + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line45A <- graph_Line45A + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line45A <- graph_Line45A +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line45A)

Line45B <- subset(all_normalized, all_normalized$Genotype == "hcLine_45B")
graph_Line45B <- ggplot(data = Line45B, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line45B <- graph_Line45B + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line45B <- graph_Line45B + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line45B <- graph_Line45B +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line45B)

Line45C <- subset(all_normalized, all_normalized$Genotype == "hcLine_45C")
graph_Line45C <- ggplot(data = Line45C, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line45C <- graph_Line45C + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line45C <- graph_Line45C + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line45C <- graph_Line45C +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line45C)

Line45D <- subset(all_normalized, all_normalized$Genotype == "hcLine_45D")
graph_Line45D <- ggplot(data = Line45D, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line45D <- graph_Line45D + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line45D <- graph_Line45D + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line45D <- graph_Line45D +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line45D)

Line50 <- subset(all_normalized, all_normalized$Genotype == "control_heat")
graph_Line50 <- ggplot(data = Line50, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line50 <- graph_Line50 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line50 <- graph_Line50 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line50 <- graph_Line50 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line50)

Line52 <- subset(all_normalized, all_normalized$Genotype == "control_drugs")
graph_Line52 <- ggplot(data = Line52, aes(x= all.min, y=area.norm, color = Treatment)) 
graph_Line52 <- graph_Line52 + geom_line(alpha = 0.1,size = 0.8, aes(group= plantID)) 
graph_Line52 <- graph_Line52 + xlab("Total Minutes After Imaging Start") + rremove("legend")
graph_Line52 <- graph_Line52 +  ylab("Normalized leaf area") + ylim(0,12.5)
ggplotly(graph_Line52)

Let’s remove all of the lines that showed very limited growth or abnormal growth

bye <- c("raspiK_cameraA_13", "raspiQ_cameraA_1", "raspiS_cameraA_10", "raspiO_cameraB_9", 
          "raspiO_cameraB_12", "raspiT_camera1_8", "raspiU_cameraA_11", "raspiV_cameraA_10", 
          "raspiN_cameraA_0", "raspiP_cameraB_0", "raspiL_cameraA_05", "raspiT_camera2_11",
          "raspiV_cameraB_13", "raspiK_cameraA_9", "raspiV_cameraB_13", "raspiS_cameraA_4",
          "raspiL_cameraB_12", "raspiK_cameraA_11", "raspiK_cameraA_12", "raspiK_cameraB_10",
          "raspiU_cameraA_6")

curated_norm <- subset(all_normalized, !(all_normalized$plantID %in% bye))

Now - let’s split it again for “Control” and “Treatment” portions of the graphs:

Control_norm <- subset(curated_norm, curated_norm$day < 20)
Treatment_norm <- subset(curated_norm, curated_norm$day > 20)

Then - let’s have a look at how individual genotype size compares with Line 50 or Line 52

mutants <- c("hcLine_02", "hcLine_04", "hcLine_07", "hcLine_09", "hcLine_18", "hcLine_31", "hcLine_33A", "hcLine_33D", "hcLine_34", "hcLine_45A", "hcLine_45B", "hcLine_45C", "hcLine_45D")

getwd()

## [1] "/Users/magdalena/Downloads/Magda-Etienne/Heat_EXP_data"

for(i in 1:length(mutants)){
  #i=1
  want_now <- c("control_heat", mutants[i])
  want_now2 <- c("control_drugs", mutants[i])
  temp_50 <- subset(Control_norm, Control_norm$Genotype %in% want_now)
  temp_52 <- subset(Control_norm, Control_norm$Genotype %in% want_now2)
  
  Control_l50 <- ggplot(data = temp_50, aes(x = all.min, y =area.norm, color = Genotype))
  Control_l50 <- Control_l50 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Control_l50 <- Control_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Control_l50 <- Control_l50 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.2, 0.8))
  Control_l50 <- Control_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Rosette_normalized_Control_", mutants[i],"_compared_to_control_heat.pdf", sep =""))
  plot(Control_l50)
  dev.off()
  
  Control_l52 <- ggplot(data = temp_52, aes(x = all.min, y =area.norm, color = Genotype))
  Control_l52 <- Control_l52 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Control_l52 <- Control_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Control_l52 <- Control_l52 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.2, 0.8))
  Control_l52 <- Control_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Rosette_normalized_Control_", mutants[i],"_compared_to_control_drugs.pdf", sep =""))
  plot(Control_l52)
  dev.off()
  
  temp_50 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now)
  temp_50 <- subset(temp_50, temp_50$Genotype %in% want_now)
  temp_52 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now2)
  
  Treatment_l50 <- ggplot(data = temp_50, aes(x = all.min, y =area.norm, color = Genotype))
  Treatment_l50 <- Treatment_l50 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l50 <- Treatment_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l50 <- Treatment_l50 + facet_wrap(~ Treatment, ncol=2)
  Treatment_l50 <- Treatment_l50 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.2, 0.8))
  Treatment_l50 <- Treatment_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Rosette_normalized_Treatment_", mutants[i],"_compared_to_control_heat.pdf", sep =""))
  plot(Treatment_l50)
  dev.off()
  
  Treatment_l52 <- ggplot(data = temp_52, aes(x = all.min, y =area.norm, color = Genotype))
  Treatment_l52 <- Treatment_l52 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l52 <- Treatment_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l52 <- Treatment_l52 + facet_wrap(~ Treatment, ncol=2)
  Treatment_l52 <- Treatment_l52 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.2, 0.8))
  Treatment_l52 <- Treatment_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Rosette_normalized_Treatment_", mutants[i],"_compared_to_control_drugs.pdf", sep =""))
  plot(Treatment_l52)
  dev.off()
}

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

## Warning: Computation failed in `stat_compare_means()`:
## '*.npc coord for x axis should be either a numeric value in [0-1] or a character strings including one of right, left, center, centre, middle

Calculating normalized growth rate:

now we have normalized rosette area - let’s calculate the growth rate per day:

colnames(curated_norm)

##  [1] "RasPi.ID"  "in_bounds" "area"      "Cam.ID"    "timestamp" "ID"       
##  [7] "month"     "day"       "hour"      "min"       "month.min" "day.min"  
## [13] "hour.min"  "all.min"   "Treatment" "Genotype"  "select"    "area.mean"
## [19] "area.norm" "plantID"

Calculate the growth rate over time based on individual rep under certain genotype

#Build the Table for multiple linear regression to be written
Growth_rate_rep <- data.frame("Month" = integer(0), 
                             "Date" = integer(0), "Treatment" = character(0),"Plant.ID" = character(0), "Genotype"= character(0), 
                             "intercept" = integer(0), "delta" = integer(0), "R.square" = integer(0))

#Build time schedule
sub_raspi_rep <- unique(curated_norm[,c(7, 8, 15, 20, 16)])
sub_raspi_rep

#Loop the calculation of linear fit of Day_time growth
for (x in 1:nrow(sub_raspi_rep)){
  i=1
  #Name initial four columns
  Growth_rate_rep[x,1:5] <- sub_raspi_rep[x, 1:5]
  
  #Calculate the intercept
  Growth_rate_rep[x,6] <- as.numeric(lm(area ~ all.min, curated_norm[curated_norm$month== as.character(sub_raspi_rep[x,1]) & 
                                        curated_norm$day == as.character(sub_raspi_rep[x,2]) & curated_norm$Treatment == as.character(sub_raspi_rep[x,3]) &                                                
                                          curated_norm$Genotype == as.character(sub_raspi_rep[x,5]) &
                                        curated_norm$plantID == as.character(sub_raspi_rep[x,4]),])$coefficient[[1]])
  
  #Calculate the delta
  Growth_rate_rep[x,7] <- as.numeric(lm(area ~ all.min, curated_norm[curated_norm$month== as.character(sub_raspi_rep[x,1]) & 
                                        curated_norm$day == as.character(sub_raspi_rep[x,2]) & curated_norm$Treatment == as.character(sub_raspi_rep[x,3]) &                                                
                                          curated_norm$Genotype == as.character(sub_raspi_rep[x,5]) &
                                        curated_norm$plantID == as.character(sub_raspi_rep[x,4]),])$coefficient[[2]])
  
  #Calculate the R.square
  Growth_rate_rep[x,8] <- summary(lm(area ~ all.min, curated_norm[curated_norm$month== as.character(sub_raspi_rep[x,1]) & 
                                        curated_norm$day == as.character(sub_raspi_rep[x,2]) & curated_norm$Treatment == as.character(sub_raspi_rep[x,3]) &                                                
                                          curated_norm$Genotype == as.character(sub_raspi_rep[x,5]) &
                                        curated_norm$plantID == as.character(sub_raspi_rep[x,4]),]))$r.squared
}
Growth_rate_rep

dim(Growth_rate_rep)

## [1] 3344    8

pos_growth_rate <- subset(Growth_rate_rep, Growth_rate_rep$delta > 0)
dim(pos_growth_rate)

## [1] 3097    8

#Clean the final table
Growth_rate_rep <- na.omit(pos_growth_rate)
head(Growth_rate_rep)

write.csv(Growth_rate_rep, "Positive_Growth_rate_normalized_all_reps.csv", row.names = FALSE)

Plot distribution of R-square

ggplot(data = Growth_rate_rep, aes(x = R.square, fill = Treatment))+ 
  geom_histogram() + facet_wrap(~ Genotype, ncol=5) +
  ylab("Count") + xlab("R spuare of linear fit")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

OK - looks great - now let’s do visualization of these growth rates as previously:

norm_gr_all <- ggplot(data = Growth_rate_rep, aes(x = Date, y = delta, color = Treatment))
norm_gr_all <- norm_gr_all + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("dodgerblue3", "tomato1")) 
norm_gr_all <- norm_gr_all + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed") 
norm_gr_all <- norm_gr_all + ylab("normalized Growth Rate") + xlab("Days after germination") 
norm_gr_all <- norm_gr_all + stat_compare_means(aes(group = Treatment), label = "p.signif", method = "aov", hide.ns = T)
norm_gr_all <- norm_gr_all + facet_wrap(~ Genotype, ncol=5)
norm_gr_all

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

## Warning: Computation failed in `stat_compare_means()`:
## Problem with `mutate()` input `p`.
## x contrasts can be applied only to factors with 2 or more levels
## ℹ Input `p` is `purrr::map(...)`.

Control <- subset(Growth_rate_rep, Growth_rate_rep$Date < 20)
Treatment <- subset(Growth_rate_rep, Growth_rate_rep$Date > 20)
Treatmen <- as.data.frame(Treatment)
head(Treatmen)

Control <- subset(Control, Control$delta > 0)
Treatment2 <- subset(Treatmen, Treatmen$delta > 0)

Control_growth <- ggplot(data = Control, aes(x = Date, y =delta, color = Genotype))
Control_growth <- Control_growth + geom_line(alpha = 0.3, aes(group= Plant.ID))
Control_growth <- Control_growth + stat_summary(fun.y=mean,  size=0.7, geom="line", linetype = "dashed")

## Warning: `fun.y` is deprecated. Use `fun` instead.

Control_growth <- Control_growth + ylab("Growth rate") + xlab("Days after germination")
# Control_growth <- Control_growth + stat_compare_means(aes(x = Date), label = "p.signif", method = "t.test", hide.ns = T)
Control_growth <- Control_growth + facet_wrap(~ Genotype, ncol=5) + rremove("legend")
Control_growth

Treatment2$DAS <- Treatment2$Date - 20

Heat_growth <- ggplot(data = Treatment2, aes(x = DAS, y =delta, color = Genotype))
Heat_growth <- Heat_growth + geom_line(alpha = 0.3, aes(group= Plant.ID))
Heat_growth <- Heat_growth + stat_summary(fun.y=mean,  size=0.7, geom="line", linetype = "dashed")

## Warning: `fun.y` is deprecated. Use `fun` instead.

Heat_growth <- Heat_growth + ylab("Growth rate") + xlab("Days after stress")
Heat_growth <- Heat_growth + facet_wrap(Treatment ~ Genotype, ncol=5) + rremove("legend")
Heat_growth

We do want to compare all of the lines to control_heat (Control for drugs) or control_drugs (Control for heat)

In this case we will have to subset into a series of subsets:

mutants <- c("hcLine_02", "hcLine_04", "hcLine_07", "hcLine_09", "hcLine_18", "hcLine_31", "hcLine_33A", "hcLine_33D", "hcLine_34", "hcLine_45A", "hcLine_45B", "hcLine_45C", "hcLine_45D")

getwd()

## [1] "/Users/magdalena/Downloads/Magda-Etienne/Heat_EXP_data"

for(i in 1:length(mutants)){
  want_now <- c("control_heat", mutants[i])
  want_now2 <- c("control_drugs", mutants[i])
  temp_50 <- subset(Control, Control$Genotype %in% want_now)
  temp_52 <- subset(Control, Control$Genotype %in% want_now2)
  
  Control_l50 <- ggplot(data = temp_50, aes(x = Date, y =delta, color = Genotype))
  Control_l50 <- Control_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Control_l50 <- Control_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Control_l50 <- Control_l50 + ylab("Growth rate") + xlab("Days after germination") + theme(legend.position = c(0.8, 0.8))
  Control_l50 <- Control_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_corrected_Control_", mutants[i],"_compared_to_control_heat.pdf", sep =""))
  plot(Control_l50)
  dev.off()
  
  Control_l52 <- ggplot(data = temp_52, aes(x = Date, y =delta, color = Genotype))
  Control_l52 <- Control_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Control_l52 <- Control_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Control_l52 <- Control_l52 + ylab("Growth rate") + xlab("Days after germination")+ theme(legend.position = c(0.8, 0.8))
  Control_l52 <- Control_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_corrected_Control_", mutants[i],"_compared_to_control_drugs.pdf", sep =""))
  plot(Control_l52)
  dev.off()
  
  temp_50 <- subset(Treatment2, Treatment2$Genotype %in% want_now)
  temp_52 <- subset(Treatment2, Treatment2$Genotype %in% want_now2)
  temp_50$DAS <- as.numeric(as.character(temp_50$DAS))
  temp_50 <- subset(temp_50, temp_50$DAS < 7)
  temp_52 <- subset(temp_52, temp_52$DAS < 7)
  
  Treatment_l50 <- ggplot(data = temp_50, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l50 <- Treatment_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l50 <- Treatment_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l50 <- Treatment_l50 + facet_wrap(~ Treatment, ncol=2)
  Treatment_l50 <- Treatment_l50 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l50 <- Treatment_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_corrected_Treatment_", mutants[i],"_compared_to_control_heat.pdf", sep ="")) 
  plot(Treatment_l50)
  dev.off()
  
  Treatment_l52 <- ggplot(data = temp_52, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l52 <- Treatment_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l52 <- Treatment_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l52 <- Treatment_l52 + facet_wrap(~ Treatment, ncol=2)
  Treatment_l52 <- Treatment_l52 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l52 <- Treatment_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  pdf(paste("Growth_corrected_Treatment_", mutants[i],"_compared_to_control_drugs.pdf", sep =""))
  plot(Treatment_l52)
  dev.off()
}

Final figures

After talking with Dr Etienne, we wanted to have a look at both Controls and compare how their growth is:

temp <- subset(Treatment_norm, Treatment_norm$Genotype == "control_heat" | Treatment_norm$Genotype == "control_drugs")
head(Treatment_norm)

Control_graph <- ggplot(data = temp, aes(x = all.min, y =area.norm, color = Genotype))
Control_graph <- Control_graph + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Control_graph <- Control_graph + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Control_graph <- Control_graph + facet_wrap(~ Treatment, ncol=2)
Control_graph <- Control_graph + ylab("Normalized rosette size") + xlab("Minutes of imaging") 
Control_graph <- Control_graph + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme_classic() + theme(legend.position = c(0.8, 0.8))
Control_graph

Another graph that the reviewer # 1, Etienne and David wanted to make is the overall comparison of variance between Controls and HC lines.

For this purpose - we need to calculate variance for each time point, HC line and condition.

variance <- summaryBy(area.norm ~ Genotype + all.min + Treatment, data = Treatment_norm, FUN = c(mean, var))

variance$HCvsWT <- "n.a."
for(i in 1:nrow(variance)){
  if(variance$Genotype[i] %in% mutants){
    variance$HCvsWT[i] <- "hc_line"
  }
  else{
    variance$HCvsWT[i] <- "control"
  }
}

head(variance)

Summary_graph <- ggplot(data = variance, aes(x = all.min, y =area.norm.var, color = HCvsWT))
Summary_graph <- Summary_graph + geom_line(alpha = 0.4, aes(group= Genotype)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Summary_graph <- Summary_graph + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Summary_graph <- Summary_graph + facet_wrap(~ Treatment, ncol=2)
Summary_graph <- Summary_graph + ylab("Variance in rosette size") + xlab("Minutes of imaging") 
Summary_graph <- Summary_graph + stat_compare_means(aes(group = HCvsWT), label = "p.signif", method = "aov", hide.ns = T) + theme_classic() + theme(legend.position = c(0.8, 0.8)) + theme(axis.text.x = element_text(angle = 90)) +labs(colour="Genotype")

Summary_graph

OK - this doesn’t even look close to what we were envisioning. So let me get to the most interesting lines:

want_now <- c("control_heat", "hcLine_04")
want_now2 <- c("control_drugs", "hcLine_04")

temp_50 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now)  
temp_52 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now2)
  
Line04_l50 <- ggplot(data = temp_50, aes(x = all.min, y =area.norm, color = Genotype))
Line04_l50 <- Line04_l50 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Line04_l50 <- Line04_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Line04_l50 <- Line04_l50 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
Line04_l50 <- Line04_l50 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.8, 0.8)) + ylim(0, 12.5)
Line04_l50 <- Line04_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme(axis.text.x = element_text(angle = 90))

Line04_l52 <- ggplot(data = temp_52, aes(x = all.min, y =area.norm, color = Genotype))
Line04_l52 <- Line04_l52 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Line04_l52 <- Line04_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Line04_l52 <- Line04_l52 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
Line04_l52 <- Line04_l52 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.8, 0.8))+ ylim(0, 12.5)
Line04_l52 <- Line04_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme(axis.text.x = element_text(angle = 90))

want_now <- c("control_heat", "hcLine_31")
want_now2 <- c("control_drugs", "hcLine_31")

temp_50 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now)  
temp_52 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now2)
  
Line31_l50 <- ggplot(data = temp_50, aes(x = all.min, y =area.norm, color = Genotype))
Line31_l50 <- Line31_l50 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Line31_l50 <- Line31_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Line31_l50 <- Line31_l50 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
Line31_l50 <- Line31_l50 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.8, 0.8))+ ylim(0, 12.5)
Line31_l50 <- Line31_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme(axis.text.x = element_text(angle = 90))

Line31_l52 <- ggplot(data = temp_52, aes(x = all.min, y =area.norm, color = Genotype))
Line31_l52 <- Line31_l52 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Line31_l52 <- Line31_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Line31_l52 <- Line31_l52 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
Line31_l52 <- Line31_l52 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.8, 0.8))+ ylim(0, 12.5)
Line31_l52 <- Line31_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme(axis.text.x = element_text(angle = 90))

want_now <- c("control_heat", "hcLine_45C")
want_now2 <- c("control_drugs", "hcLine_45C")

temp_50 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now)  
temp_52 <- subset(Treatment_norm, Treatment_norm$Genotype %in% want_now2)
  
Line45C_l50 <- ggplot(data = temp_50, aes(x = all.min, y =area.norm, color = Genotype))
Line45C_l50 <- Line45C_l50 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Line45C_l50 <- Line45C_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Line45C_l50 <- Line45C_l50 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
Line45C_l50 <- Line45C_l50 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.8, 0.8))+ ylim(0, 12.5)
Line45C_l50 <- Line45C_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme(axis.text.x = element_text(angle = 90))

Line45C_l52 <- ggplot(data = temp_52, aes(x = all.min, y =area.norm, color = Genotype))
Line45C_l52 <- Line45C_l52 + geom_line(alpha = 0.3, aes(group= plantID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
Line45C_l52 <- Line45C_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
Line45C_l52 <- Line45C_l52 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
Line45C_l52 <- Line45C_l52 + ylab("Normalized rosette size") + xlab("Minutes of imaging") + theme(legend.position = c(0.8, 0.8))+ ylim(0, 12.5)
Line45C_l52 <- Line45C_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T) + theme(axis.text.x = element_text(angle = 90))

OK - now let’s combine all of the graphs into one figure:

library(cowplot)

## 
## ********************************************************

## Note: As of version 1.0.0, cowplot does not change the

##   default ggplot2 theme anymore. To recover the previous

##   behavior, execute:
##   theme_set(theme_cowplot())

## ********************************************************

## 
## Attaching package: 'cowplot'

## The following object is masked from 'package:ggpubr':
## 
##     get_legend

mid_row <- plot_grid(Line04_l50, Line31_l50, Line45C_l50, labels = c('B', 'C', 'D'), label_size = 12, ncol = 3)
bottom_row <- plot_grid(Line04_l52, Line31_l52, Line45C_l52, labels = c('E', 'F', 'G'), label_size = 12, ncol = 3)

mid_row

bottom_row

pdf("Final_Figure_Phenotyping.pdf", width = 10, height = 17)
plot_grid(Summary_graph, mid_row, bottom_row, ncol = 1, labels=c("A", "", ""))
dev.off()

## quartz_off_screen 
##                 2

OK - since the reviewer is not convinced by our increase in rosette size over time, we should also add graphs representing the growth rates for all of the above mentioned lines:

mutants <- c("hcLine_04", "hcLine_31", "hcLine_45C")

# Line 04
i=1
  want_now <- c("control_heat", mutants[i])
  want_now2 <- c("control_drugs", mutants[i])
  temp_50 <- subset(Treatment2, Treatment2$Genotype %in% want_now)
  temp_52 <- subset(Treatment2, Treatment2$Genotype %in% want_now2)
  temp_50$DAS <- as.numeric(as.character(temp_50$DAS))
  temp_50 <- subset(temp_50, temp_50$DAS < 7)
  temp_52 <- subset(temp_52, temp_52$DAS < 7)
  
  Treatment_l50 <- ggplot(data = temp_50, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l50 <- Treatment_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l50 <- Treatment_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l50 <- Treatment_l50 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
  Treatment_l50 <- Treatment_l50 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l50 <- Treatment_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  
  Treatment_l52 <- ggplot(data = temp_52, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l52 <- Treatment_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l52 <- Treatment_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l52 <- Treatment_l52 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
  Treatment_l52 <- Treatment_l52 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l52 <- Treatment_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  
top_row <- plot_grid(Treatment_l50, Treatment_l52, labels = c('A', 'B'), label_size = 12, ncol = 2)

# Line 31
i=2
  want_now <- c("control_heat", mutants[i])
  want_now2 <- c("control_drugs", mutants[i])
  temp_50 <- subset(Treatment2, Treatment2$Genotype %in% want_now)
  temp_52 <- subset(Treatment2, Treatment2$Genotype %in% want_now2)
  temp_50$DAS <- as.numeric(as.character(temp_50$DAS))
  temp_50 <- subset(temp_50, temp_50$DAS < 7)
  temp_52 <- subset(temp_52, temp_52$DAS < 7)
  
  Treatment_l50 <- ggplot(data = temp_50, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l50 <- Treatment_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l50 <- Treatment_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l50 <- Treatment_l50 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
  Treatment_l50 <- Treatment_l50 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l50 <- Treatment_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  
  Treatment_l52 <- ggplot(data = temp_52, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l52 <- Treatment_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l52 <- Treatment_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l52 <- Treatment_l52 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
  Treatment_l52 <- Treatment_l52 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l52 <- Treatment_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  
mid_row <- plot_grid(Treatment_l50, Treatment_l52, labels = c('C', 'D'), label_size = 12, ncol = 2)
  
# Line 45C
i=3
  want_now <- c("control_heat", mutants[i])
  want_now2 <- c("control_drugs", mutants[i])
  temp_50 <- subset(Treatment2, Treatment2$Genotype %in% want_now)
  temp_52 <- subset(Treatment2, Treatment2$Genotype %in% want_now2)
  temp_50$DAS <- as.numeric(as.character(temp_50$DAS))
  temp_50 <- subset(temp_50, temp_50$DAS < 7)
  temp_52 <- subset(temp_52, temp_52$DAS < 7)
  
  Treatment_l50 <- ggplot(data = temp_50, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l50 <- Treatment_l50 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l50 <- Treatment_l50 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l50 <- Treatment_l50 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
  Treatment_l50 <- Treatment_l50 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l50 <- Treatment_l50 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  
  Treatment_l52 <- ggplot(data = temp_52, aes(x = DAS, y =delta, color = Genotype))
  Treatment_l52 <- Treatment_l52 + geom_line(alpha = 0.3, aes(group= Plant.ID)) + scale_colour_manual(values = c("tomato1", "dodgerblue3"))
  Treatment_l52 <- Treatment_l52 + stat_summary(fun=mean,  size=0.7, geom="line", linetype = "dashed")
  Treatment_l52 <- Treatment_l52 + facet_wrap(~ Treatment, ncol=2) + theme_classic() 
  Treatment_l52 <- Treatment_l52 + ylab("Growth rate") + xlab("Days after stress") + theme(legend.position = c(0.8, 0.8))
  Treatment_l52 <- Treatment_l52 + stat_compare_means(aes(group = Genotype), label = "p.signif", method = "aov", hide.ns = T)
  
bottom_row <- plot_grid(Treatment_l50, Treatment_l52, labels = c('E', 'F'), label_size = 12, ncol = 2)

Finally - let’s combine all of the figures into one pdf:

pdf("Supplemental_Figure_growth_Phenotyping.pdf", width = 13, height = 17)
plot_grid(top_row, mid_row, bottom_row, ncol = 1, labels=c("", "", ""))
dev.off()

## quartz_off_screen 
##                 2