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)

## Registered S3 methods overwritten by 'tibble':
##   method     from  
##   format.tbl pillar
##   print.tbl  pillar

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.2     ✓ purrr   0.3.4
## ✓ tibble  3.0.3     ✓ 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)
library(directlabels)
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 <- read.table("/Users/magda/Documents/DataAndAnalysis/collaborations/Magda-Etienne/Heat_EXP_data/sample.list", 
                          sep ="\t", header = F)
sample.list

for (i in 1:nrow(sample.list)){
  sample.df <- read.csv(paste0("/Users/magda/Documents/DataAndAnalysis/collaborations/Magda-Etienne/Heat_EXP_data/",
                               sample.list[i,1]), header = T)
  assign(sample.list[i,1], sample.df)
}

rm(sample.df)

#Create Meta phenotype sheet by combining multiple files
merge_list <- lapply(ls(pattern = "Raspi"), get)
raspi_combined <- bind_rows(merge_list)

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/magda/Documents/DataAndAnalysis/collaborations/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/magda/Documents/DataAndAnalysis/collaborations/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(all.min ~ area, 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(all.min ~ area, 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(all.min ~ area, 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 spuare 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) + 
    #facet based on genotype
    facet_wrap(~ Genotype, ncol=5)

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

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

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

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

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/magda/Documents/DataAndAnalysis/collaborations/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/magda/Documents/DataAndAnalysis/collaborations/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"