Experimental description:

The seeds were sterilized in 50% bleach for 10 min with gentle shaking, then rinsed 5 times with deionized water, and incubated at 4 °C overnight. Then seeds were germinated in ¼ MS + sucrose (with vitamins) control plate for 4 complete days in growth chamber #35. At d5, the seedlings were transferred to either 100 mM salt or control plates, and started being scanned for 5 total consecutive days (starting from d5 to d9). The growth chamber condition was 26 °C and 20 h light /4 h dark cycle, with 60% humidity (Chamber #35).The root tracing was done using ImageJ SamrtRoot plugin. Shoot and root fresh weight were measured after 10 d in treatment plates.

getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20211222_ER_3_tomato_accessions_RSA_2nd_replicate"
setwd("C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20211222_ER_3_tomato_accessions_RSA_2nd_replicate/")
list.files(pattern = ".csv")
## [1] "2022jan_ER_growth_factors.csv" "ER_d5.csv"                    
## [3] "ER_d6.csv"                     "ER_d7.csv"                    
## [5] "ER_d8.csv"                     "ER_d9.csv"
my_data_d5<- read.csv("ER_d5.csv")
my_data_d6<- read.csv("ER_d6.csv")
my_data_d7<- read.csv("ER_d7.csv")
my_data_d8<- read.csv("ER_d8.csv")
my_data_d9<- read.csv("ER_d9.csv")
head(my_data_d5)
head(my_data_d6)
head(my_data_d9)
library(readr)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
my_data_d5<-my_data_d5%>%mutate(day=5,my_data_d5 )
head(my_data_d5)
my_data_d6<-my_data_d6%>%mutate(day=6,my_data_d6 )
head(my_data_d6)
my_data_d7<-my_data_d7%>%mutate(day=7,my_data_d7 )
head(my_data_d7)
my_data_d8<-my_data_d8%>%mutate(day=8,my_data_d8 )
head(my_data_d8)
my_data_d9<-my_data_d9%>%mutate(day=9,my_data_d9 )
head(my_data_d9)

now I fuse days

all_data<- rbind(my_data_d5,my_data_d6)
all_data<-rbind(all_data, my_data_d7)
all_data<-rbind(all_data, my_data_d8)
all_data<-rbind(all_data, my_data_d9)
all_data
unique(all_data$day)
## [1] 5 6 7 8 9
Main_root <- subset(all_data, all_data$root_order == 0)
Main_root
colnames(Main_root)
##  [1] "image"                 "root_name"             "root"                 
##  [4] "length"                "vector_length"         "surface"              
##  [7] "volume"                "direction"             "diameter"             
## [10] "root_order"            "root_ontology"         "parent_name"          
## [13] "parent"                "insertion_position"    "insertion_angle"      
## [16] "n_child"               "child_density"         "first_child"          
## [19] "insertion_first_child" "last_child"            "insertion_last_child" 
## [22] "day"
MR_data<-Main_root[, c(1:4, 10, 16,17, 19:22)]
MR_data
Lateral_root <- subset(all_data, all_data$root_order != 0)
Lateral_root
unique(Lateral_root$root_order)
## [1] 1 2
colnames(Lateral_root)
##  [1] "image"                 "root_name"             "root"                 
##  [4] "length"                "vector_length"         "surface"              
##  [7] "volume"                "direction"             "diameter"             
## [10] "root_order"            "root_ontology"         "parent_name"          
## [13] "parent"                "insertion_position"    "insertion_angle"      
## [16] "n_child"               "child_density"         "first_child"          
## [19] "insertion_first_child" "last_child"            "insertion_last_child" 
## [22] "day"
MR_data2 <- subset(MR_data, MR_data$n_child >0)
MR_data2

I have started to confuse now….by [1], etc…..

temporary <- subset(Lateral_root, Lateral_root$parent == MR_data2$root[1])
temporary
temporary <- subset(temporary, temporary$day == MR_data2$day[3])
temporary

again what is dim here tries to do????

dim(temporary)
## [1]  1 22
dim(temporary)[1]
## [1] 1
dim(temporary)[2]
## [1] 22
total_LRL <- sum(temporary$length)
LR_number <- dim(temporary)[1]
total_LRL
## [1] 0.1701659
LR_number
## [1] 1
MR_data2$LRL <- 0
MR_data2$LRno <- 0
MR_data2
MR_data2$LRL[1] <- total_LRL
MR_data2$LRno[1] <- LR_number
MR_data2
MR_data_noChild <- subset(MR_data, MR_data$n_child == 0)
MR_data_noChild
length(MR_data2$root)
## [1] 182
for(i in 2:182){
  temporary <- subset(Lateral_root, Lateral_root$parent == MR_data2$root[i])
  temporary <- subset(temporary, temporary$day == MR_data2$day[i])
  total_LRL <- sum(temporary$length)
  LR_number <- dim(temporary)[1]
  MR_data2$LRL[i] <- total_LRL
  MR_data2$LRno[i] <- LR_number
}

MR_data2$check <- MR_data2$n_child - MR_data2$LRno
MR_data2
unique(MR_data2$check)
## [1] 0
#let's remove check column from Child:

MR_data_Child2 <- MR_data2[,1:13]
MR_data_Child2
MR_data_noChild
MR_data_noChild$LRL <- 0
MR_data_noChild$LRno <- 0


MR_all <- rbind(MR_data_Child2, MR_data_noChild)
MR_all
MR_all$root_name[1]
## [1] " pl3_c_a"
text <- strsplit(x = MR_all$root_name[1], split = "_")
text
## [[1]]
## [1] " pl3" "c"    "a"
plate <- text[[1]][1]
plate <- gsub(" ", "", plate)
plate
## [1] "pl3"
cond <- text[[1]][2]
cond
## [1] "c"
genotype <- text[[1]][3]
genotype
## [1] "a"
dim(MR_all)
## [1] 240  13
for(i in 1:240){
  text <- strsplit(x = MR_all$root_name[i], split = "_")
  plate <- text[[1]][1]
  cond <- text[[1]][2]
  genotype <- text[[1]][3]
   
  MR_all$genotype[i] <- genotype
  MR_all$condition[i] <- cond
  MR_all$plate[i] <- plate
}

MR_all
MR_all$TRS <- MR_all$length + MR_all$LRL
MR_all$aLRL <- MR_all$LRL/ MR_all$LRno
MR_all$MRpLRL <- MR_all$length / MR_all$LRL
MR_all
length(unique(MR_all$root_name))
## [1] 49
unique(MR_all$day)
## [1] 5 6 7 8 9
MR_all$day <- as.numeric(as.character(MR_all$day))
unique(MR_all$day)
## [1] 5 6 7 8 9

thus dividing images / plants should give us number of days: but it does not make sense….cannot be 4.8, it should be 5.

240/49
## [1] 4.897959
library(ggplot2)
library(ggpubr)

unique(MR_all$condition)
## [1] "c" "s" "0"
histogram_TRS_ER <- ggdensity(MR_all, x = "TRS",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition")
                           #palette = c("#00AFBB", "#E7B800"))

histogram_TRS_ER
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf

we have a problem - since we have three conditions - rather than two expected ones - lets check what are the roots with “0” as condition (should be “c”):

MR_all2 <- subset(MR_all, MR_all$condition != "0")
MR_funky <- subset(MR_all, MR_all$condition == "0")
MR_funky
# wait - but they also have no genotype - let's check images
temp <- subset(MR_all2, MR_all2$image == "EricsTomatoLines_202102015011.rsml")
temp
temp2 <- subset(MR_all2, MR_all2$image == "EricsTomatoLines_202102017025.rsml")
temp2

conclusion >> the two roots with condition “0” are actually lateral roots that were not classified as such (whatever reason) - and thus we can ignore them and remove them from the dataset - continuing with MR_all2 onwards

MR_all2$TRS <- as.numeric(MR_all2$TRS)
unique(MR_all2$condition)
## [1] "c" "s"
histogram_TRS_ER <- ggdensity(MR_all2, x = "TRS",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition",
                           palette = c("#00AFBB", "#E7B800"))

histogram_TRS_ER

histogram_LRno_ER <- ggdensity(MR_all2, x = "LRno",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition",
                           palette = c("#00AFBB", "#E7B800"))

histogram_LRno_ER

pdf("histogram.TRS.ER.pdf")
plot(histogram_TRS_ER)
# if plotting multiple graphs - this command is extremely important 
dev.off()
## png 
##   2
pdf("histogram.LRno.ER.pdf")
plot(histogram_LRno_ER)
# if plotting multiple graphs - this command is extremely important 
dev.off()
## png 
##   2
library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
pdf("Figure_MAIN_1.pdf", height = 15, width = 12)
plot_grid(histogram_TRS_ER, histogram_LRno_ER, ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## png 
##   2
unique(MR_all2$genotype)
## [1] "a" "b" "c"
TRS_lgraph <- ggplot(data=MR_all2, aes(x= genotype, y=TRS, color = condition)) 
TRS_lgraph <- TRS_lgraph + geom_boxplot()
TRS_lgraph <- TRS_lgraph + facet_grid(day ~ condition, scales = "free") + scale_color_manual(values=c("turquoise3", "maroon3", "dark orange"))
TRS_lgraph <- TRS_lgraph + ylab("Total root size (cm)") + xlab("Genotype") + theme(legend.position='none')
TRS_lgraph

library(ggsci)
library(ggbeeswarm)
library(gapminder)
library(RColorBrewer)
library(ggridges)
better_TRS_graph <- ggplot(data=MR_all2, aes(x= genotype, y=TRS, color = genotype))
better_TRS_graph <- better_TRS_graph + geom_beeswarm(alpha=0.6, priority = "density")
better_TRS_graph <- better_TRS_graph + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
better_TRS_graph <- better_TRS_graph + facet_grid(day ~ condition, scales = "free") + scale_color_manual(values=c("turquoise3", "maroon3", "dark orange"))
better_TRS_graph <- better_TRS_graph + ylab("Total root size (cm)") + xlab("Genotype") + theme(legend.position='none')
better_TRS_graph <- better_TRS_graph + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
better_TRS_graph <- better_TRS_graph + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "a", hide.ns = TRUE) 
better_TRS_graph

#install.packages("plotly")
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
head(MR_all2)
my_graph <- ggplot(data=MR_all2, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

Curating the data

OK - now let’s subset the data into individual genotypes and have a look at which of the plants we want to remove from the dataset:

geno_a <- subset(MR_all2, MR_all2$genotype == "a")

my_graph <- ggplot(data=geno_a, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

There is one pant that didnot grow but then stops to grow (pl6_c_a), so let;s remove them:

geno_a$root_name <- gsub(" ", "", geno_a$root_name)
bye_a <- c("pl6_c_a")
OK_geno_a <- subset(geno_a, !(geno_a$root_name %in% bye_a))
unique(OK_geno_a$root_name)
##  [1] "pl3_c_a" "pl4_c_a" "pl4_s_a" "pl8_c_a" "pl1_c_a" "pl2_s_a" "pl2_c_a"
##  [8] "pl5_s_a" "pl5_c_a" "pl7_s_a" "pl7_c_a" "pl1_s_a" "pl3_s_a" "pl6_s_a"
## [15] "pl8_s_a"
my_graph <- ggplot(data=OK_geno_a, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_a, aes(x= day, y=TRS, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Total Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_a, aes(x= day, y=LRno, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_a, aes(x= day, y=LRL, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

so let’s move to the next genotype:

geno_b <- subset(MR_all2, MR_all2$genotype == "b")
geno_b
my_graph <- ggplot(data=geno_b, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

Most of the above look good except plant “pl2_c_b” and “pl4_c_b” and “pl1_s_b”so let’s remove them:

geno_b$root_name <- gsub(" ", "", geno_b$root_name)
bye_b <- c("pl2_c_b","pl4_c_b","pl1_s_b")
OK_geno_b <- subset(geno_b, !(geno_b$root_name %in% bye_b))
unique(OK_geno_b$root_name)
##  [1] "pl2_s_b" "pl3_s_b" "pl1_c_b" "pl3_c_b" "pl4_s_b" "pl5_c_b" "pl6_s_b"
##  [8] "pl6_c_b" "pl7_c_b" "pl8_s_b" "pl8_c_b" "pl5_s_b" "pl7_s_b"
my_graph <- ggplot(data=OK_geno_b, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_b, aes(x= day, y=TRS, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Total Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_b, aes(x= day, y=LRno, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_b, aes(x= day, y=LRL, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

And finally let’s look at genotype C:

geno_c <- subset(MR_all2, MR_all2$genotype == "c")
my_graph <- ggplot(data=geno_c, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

library(plotly)
ggplotly(my_graph)

There is one pant that didnot grow at all (pl8_s_c), and one plant that grows very odd (pl5_s_c), so let;s remove them:

geno_c$root_name <- gsub(" ", "", geno_c$root_name)
bye_c <- c("pl8_s_c", "pl5_s_c")
OK_geno_c <- subset(geno_c, !(geno_c$root_name %in% bye_c))
unique(OK_geno_c$root_name)
##  [1] "pl1_c_c" "pl2_c_c" "pl3_c_c" "pl4_c_c" "pl5_c_c" "pl6_c_c" "pl7_c_c"
##  [8] "pl8_c_c" "pl3_s_c" "pl4_s_c" "pl7_s_c" "pl1_s_c" "pl6_s_c" "pl2_s_c"
my_graph <- ggplot(data=OK_geno_c, aes(x= day, y=length, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_c, aes(x= day, y=TRS, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Total Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_c, aes(x= day, y=LRno, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=OK_geno_c, aes(x= day, y=LRL, group = root_name, color = genotype)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

MR_OK <- rbind(OK_geno_a, OK_geno_b)
MR_OK <- rbind(MR_OK, OK_geno_c)

MR_time_graph <- ggplot(data=MR_OK, aes(x= day, y=length, group = root_name, color = genotype)) 
MR_time_graph <- MR_time_graph + geom_line(alpha = 0.2) 
MR_time_graph <- MR_time_graph + facet_grid(~ condition) 
MR_time_graph <- MR_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("Main Root Length")
MR_time_graph <- MR_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
MR_time_graph <- MR_time_graph + scale_y_continuous(name = " length (cm)", breaks = c(0:10))
MR_time_graph <- MR_time_graph + scale_color_manual(values=c("blue", "dark orange", "red"))

MR_time_graph

LRL_time_graph <- ggplot(data=MR_OK, aes(x= day, y=LRL, group = root_name, color = genotype)) 
LRL_time_graph <- LRL_time_graph + geom_line(alpha = 0.2) 
LRL_time_graph <- LRL_time_graph + facet_grid(~ condition) 
LRL_time_graph <- LRL_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("Lateral Root Length")
LRL_time_graph <- LRL_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
LRL_time_graph

LRno_time_graph <- ggplot(data=MR_OK, aes(x= day, y=LRno, group = root_name, color = genotype)) 
LRno_time_graph <- LRno_time_graph + geom_line(alpha = 0.2) 
LRno_time_graph <- LRno_time_graph + facet_grid(~ condition) 
LRno_time_graph <- LRno_time_graph + ylab("Lateral Root Number") + xlab("time (days after germination)") + ggtitle("Lateral Root Number")
LRno_time_graph <- LRno_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
LRno_time_graph

aLRL_time_graph <- ggplot(data=MR_OK, aes(x= day, y=aLRL, group = root_name, color = genotype)) 
aLRL_time_graph <- aLRL_time_graph + geom_line(alpha = 0.2) 
aLRL_time_graph <- aLRL_time_graph + facet_grid(~ condition) 
aLRL_time_graph <- aLRL_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("average Lateral Root Length")
aLRL_time_graph <- aLRL_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
aLRL_time_graph
## Warning: Removed 52 rows containing non-finite values (`stat_summary()`).
## Warning: Removed 52 rows containing missing values (`geom_line()`).

TRS_time_graph <- ggplot(data=MR_OK, aes(x= day, y=TRS, group = root_name, color = genotype)) 
TRS_time_graph <- TRS_time_graph + geom_line(alpha = 0.2) 
TRS_time_graph <- TRS_time_graph + facet_grid(~ condition) 
TRS_time_graph <- TRS_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("Total Root Length")
TRS_time_graph <- TRS_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
TRS_time_graph <- TRS_time_graph + scale_color_manual(values=c("blue", "dark orange", "red"))

TRS_time_graph

### Growth calculations

So for calculating growth rate - lets first establish calculations on one plant. Let’s take the first plant that is within the experiment:

temp1 <- subset(MR_OK, MR_OK$root_name == unique(MR_OK$root_name)[1])
temp2 <- temp1[order(temp1$day),]
temp2
# For Main Root Growth Rate - we want to remove all the datapoints that are repeating, because that indicates root hitting the plate edge:
temp_MR <- temp2[,c("day", "length")]
plot(temp_MR$length~ temp_MR$day)

temp_MR$MRdouble <- "no"
  for(i in 2:nrow(temp_MR)){
    # we want the root to be at least 1 mm larger than the previous day - all the other ones will just indicate noise:
   if(temp_MR$length[i] <= temp_MR$length[i-1]+0.09){
    temp_MR$MRdouble[i] <- "yes"   
   } 
    else{
    temp_MR$MRdouble[i] <- "no"   
   } 
  }
temp_MR
temp_MR2 <- subset(temp_MR, temp_MR$MRdouble == "no")
temp_MR2 <- subset(temp_MR2, temp_MR2$day < 8)
temp_MR2
plot(temp_MR2$length~ temp_MR2$day)
# let's add the regression line to this graph
abline(lm(temp_MR2$length~ temp_MR2$day))

MR_model <- lm(temp_MR2$length~ temp_MR2$day)
MR_model
## 
## Call:
## lm(formula = temp_MR2$length ~ temp_MR2$day)
## 
## Coefficients:
##  (Intercept)  temp_MR2$day  
##       -6.317         2.163
summary(MR_model)
## 
## Call:
## lm(formula = temp_MR2$length ~ temp_MR2$day)
## 
## Residuals:
##         1         2         3 
## -0.002641  0.005281 -0.002641 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -6.316808   0.027696  -228.1  0.00279 **
## temp_MR2$day  2.162885   0.004574   472.9  0.00135 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.006468 on 1 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:      1 
## F-statistic: 2.236e+05 on 1 and 1 DF,  p-value: 0.001346
MR_growth_rate <- MR_model$coefficients[[2]]
# Let's remove all NA for LRno and aLRL in temp2
LR_temp <- temp2[,c("day", "LRno", "aLRL")]
LR_temp2 <- na.omit(LR_temp)

# Let's start with LRno
plot(LR_temp2$LRno ~ LR_temp2$day)
abline(lm(LR_temp2$LRno ~ LR_temp2$day))

LRno_model <- lm(LR_temp2$LRno ~ LR_temp2$day)
LRno_model
## 
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
## 
## Coefficients:
##  (Intercept)  LR_temp2$day  
##        -22.4           4.6
summary(LRno_model)
## 
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
## 
## Residuals:
##    1    2    3    4    5 
##  0.4 -2.2  2.2  0.6 -1.0 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   -22.400      4.364  -5.134  0.01432 * 
## LR_temp2$day    4.600      0.611   7.529  0.00486 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.932 on 3 degrees of freedom
## Multiple R-squared:  0.9497, Adjusted R-squared:  0.933 
## F-statistic: 56.68 on 1 and 3 DF,  p-value: 0.004858
LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase
## [1] 4.6
# Let's start with LRno
plot(LR_temp2$aLRL ~ LR_temp2$day)
abline(lm(LR_temp2$aLRL ~ LR_temp2$day))

aLRL_model <- lm(LR_temp2$aLRL ~ LR_temp2$day)
aLRL_model
## 
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
## 
## Coefficients:
##  (Intercept)  LR_temp2$day  
##      -1.4744        0.3326
summary(aLRL_model)
## 
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
## 
## Residuals:
##        1        2        3        4        5 
## -0.01834  0.20510 -0.24279 -0.05636  0.11239 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -1.47440    0.44628  -3.304   0.0456 *
## LR_temp2$day  0.33258    0.06249   5.322   0.0130 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1976 on 3 degrees of freedom
## Multiple R-squared:  0.9042, Adjusted R-squared:  0.8723 
## F-statistic: 28.32 on 1 and 3 DF,  p-value: 0.01296
aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth
## [1] 0.3325816
names <- c(text="root_name", "genotype", "cond", "MR.delta", "LRno.delta", "aLRL.delta")
growth_factors <- data.frame()

for (k in names) growth_factors[[k]] <- as.character()

growth_factors[1,1] <- temp2$root_name[1]
growth_factors[1,2] <- temp2$genotype[1]
growth_factors[1,3] <- temp2$cond[1]
growth_factors[1,4] <- as.numeric(as.character(MR_growth_rate))
growth_factors[1,5] <- as.numeric(as.character(LRno_increase))
growth_factors[1,6] <- as.numeric(as.character(aLRL_growth))

growth_factors
length(unique(MR_OK$root_name))
## [1] 42
# we are starting the loop from 2nd plant (i in 2:...), because we already calculated growth rates for the 1st plant:
for(e in 1:42){
  temp1 <- subset(MR_OK, MR_OK$root_name == unique(MR_OK$root_name)[e])
  temp2 <- temp1[order(temp1$day),]
  temp2

############ MR calculations
  temp_MR <- temp2[,c("day", "length")]
  temp_MR$MRdouble <- "no"
  for(i in 2:nrow(temp_MR)){
    # we want the root to be at least 1 mm larger than the previous day - all the other ones will just indicate noise:
   if(temp_MR$length[i] <= temp_MR$length[i-1]+0.09){
    temp_MR$MRdouble[i] <- "yes"   
   } else{temp_MR$MRdouble[i] <- "no"}}
  temp_MR2 <- subset(temp_MR, temp_MR$MRdouble == "no")
  temp_MR2 <- subset(temp_MR2, temp_MR2$day < 8)
  temp_MR2
  MR_model <- lm(temp_MR2$length~ temp_MR2$day)
  MR_growth_rate <- MR_model$coefficients[[2]]
  MR_growth_rate
  ############ LRno calculations
  LR_temp <- temp2[,c("day", "LRno", "aLRL")]
  LR_temp2 <- na.omit(LR_temp)
  LR_temp2
  dim(LR_temp2)
  
  ####################### safety precaution to calculate LR growth rate only for the plants that have LR at least for two days: 
  if(dim(LR_temp2)[1] > 1){
  LRno_model <- lm(LR_temp2$LRno ~ LR_temp2$day)
  LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))

  ############ aLRL calculations
  aLRL_model <- lm(LR_temp2$aLRL ~ LR_temp2$day)
  aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
  } else{
  ####################### safety precaution continued:
  ####################### so if you only have one day where LR are there - this wont be good enough to calculate LRno or LRL rate
  ####################### and thus:
  LRno_increase <- 0
  aLRL_growth <- 0
  }
  LRno_increase
  aLRL_growth
  ############ adding the information to the table:
  growth_factors[e,1] <- temp2$root_name[1]
  growth_factors[e,2] <- temp2$genotype[1]
  growth_factors[e,3] <- temp2$cond[1]
  growth_factors[e,4] <- as.numeric(as.character(MR_growth_rate))
  growth_factors[e,5] <- as.numeric(as.character(LRno_increase))
  growth_factors[e,6] <- as.numeric(as.character(aLRL_growth))
}

growth_factors
write.csv(growth_factors, "2022jan_ER_growth_factors.csv", row.names = FALSE)
growth_factors$MR.delta <- as.numeric(as.character(growth_factors$MR.delta))
growth_factors$aLRL.delta <- as.numeric(as.character(growth_factors$aLRL.delta))
growth_factors$LRno.delta <- as.numeric(as.character(growth_factors$LRno.delta))

MR.delta_p_geno <- ggerrorplot(growth_factors, y = "MR.delta", x = "genotype", fill="genotype", color="genotype", 
                        facet.by = c("cond"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="Genotype", ylab="Growth Rate (cm / day)") 
MR.delta_p_geno <- MR.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "a", 
                                                              label = "p.signif", hide.ns = T)
MR.delta_p_geno <- MR.delta_p_geno + ggtitle("Main Root Growth")
MR.delta_p_geno <- MR.delta_p_geno + scale_color_manual(values=c("blue", "dark orange", "red"))
MR.delta_p_geno

LRno.delta_p_geno <- ggerrorplot(growth_factors, y = "LRno.delta", x = "genotype", fill="genotype", color="genotype", 
                        facet.by = c("cond"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="Genotype", ylab="LR increase Rate (# LR / day)") 
LRno.delta_p_geno <- LRno.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "a", 
                                                              label = "p.signif", hide.ns = T)
LRno.delta_p_geno <- LRno.delta_p_geno + ggtitle("Lateral Root Number Increase")
LRno.delta_p_geno

aLRL.delta_p_geno <- ggerrorplot(growth_factors, y = "aLRL.delta", x = "genotype", fill="genotype", color="genotype", 
                        facet.by = c("cond"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="Genotype", ylab="Growth Rate (cm / day)") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "a", 
                                                              label = "p.signif", hide.ns = T)
aLRL.delta_p_geno <- aLRL.delta_p_geno + ggtitle("average Lateral Root Growth")
aLRL.delta_p_geno

### Calculating the relative growth rates (Salt Tolerance Indexes):

So in order to calculate relative plant performance at salt, we need to divide all the growth rates by the average value for that specific genotype under control conditions.

Let’s start with calculating the average growth rate per accession per condition then:

library(doBy)
## 
## Attaching package: 'doBy'
## The following object is masked from 'package:dplyr':
## 
##     order_by
avg_growth <- summaryBy(data = growth_factors, MR.delta + aLRL.delta + LRno.delta ~ genotype + cond)
avg_growth
avg_growth_C <- subset(avg_growth, avg_growth$cond == "c")
avg_growth_C <- avg_growth_C[,c(1,3:5)]
colnames(avg_growth_C) <- gsub(".mean", ".avg.Control", colnames(avg_growth_C))
avg_growth_C
STI_growth_factors <- merge(growth_factors, avg_growth_C, id="genotype")
STI_growth_factors
STI_growth_factors$MR.STI <- STI_growth_factors$MR.delta / STI_growth_factors$MR.delta.avg.Control
STI_growth_factors$aLRL.STI <- STI_growth_factors$aLRL.delta / STI_growth_factors$aLRL.delta.avg.Control
STI_growth_factors$LRno.STI <- STI_growth_factors$LRno.delta / STI_growth_factors$LRno.delta.avg.Control
head(STI_growth_factors)
STI <- subset(STI_growth_factors, STI_growth_factors$cond == "s")
MR.STI_plot <- ggplot(data = STI, mapping = aes(x = genotype, y = MR.STI, colour = genotype)) 
MR.STI_plot <- MR.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
MR.STI_plot <- MR.STI_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
MR.STI_plot <- MR.STI_plot + theme(legend.position = "none")
MR.STI_plot <- MR.STI_plot + xlab("")
MR.STI_plot <- MR.STI_plot + ylab("Fraction of control") + ggtitle("Salt Tolerance Index based on Main Root Growth")
MR.STI_plot <- MR.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "a") 
MR.STI_plot <- MR.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
MR.STI_plot

aLRL.STI_plot <- ggplot(data = STI, mapping = aes(x = genotype, y = aLRL.STI, colour = genotype)) 
aLRL.STI_plot <- aLRL.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
aLRL.STI_plot <- aLRL.STI_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
aLRL.STI_plot <- aLRL.STI_plot + theme(legend.position = "none")
aLRL.STI_plot <- aLRL.STI_plot + xlab("")
aLRL.STI_plot <- aLRL.STI_plot + ylab("Fraction of control")  + ggtitle("Salt Tolerance Index based on average Lateral Root Growth")
aLRL.STI_plot <- aLRL.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "a") 
aLRL.STI_plot <- aLRL.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
aLRL.STI_plot

LRno.STI_plot <- ggplot(data = STI, mapping = aes(x = genotype, y = LRno.STI, colour = genotype)) 
LRno.STI_plot <- LRno.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
LRno.STI_plot <- LRno.STI_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
LRno.STI_plot <- LRno.STI_plot + theme(legend.position = "none")
LRno.STI_plot <- LRno.STI_plot + xlab("")
LRno.STI_plot <- LRno.STI_plot + ylab("Fraction of control")  + ggtitle("Salt Tolerance Index based on average Increase in Lateral Root Number")
LRno.STI_plot <- LRno.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "a") 
LRno.STI_plot <- LRno.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
LRno.STI_plot

library(cowplot)

pdf("Figure1_MainRootLength.pdf", height = 6, width = 12)
plot_grid(MR_time_graph,
          align = "hv", labels=c(""), 
          label_size = 24)
dev.off()
## png 
##   2
pdf("Figure2_LateralRootLength.pdf", height = 6, width = 12)
plot_grid(LRL_time_graph,
          align = "hv", labels=c(""), 
          label_size = 24)
dev.off()
## png 
##   2
pdf("Figure3_LateralRootNumber.pdf", height = 6, width = 12)
plot_grid(LRno_time_graph,
          align = "hv", labels=c(""), 
          label_size = 24)
dev.off()
## png 
##   2
pdf("Figure4_MainRootGrowth.pdf", height = 6, width = 12)
plot_grid(MR.delta_p_geno,
          align = "hv", labels=c(""), 
          label_size = 24)
dev.off()
## png 
##   2
library(cowplot)
plot_grid(MR_time_graph, LRL_time_graph, LRno_time_graph, labels = c("AUTO"), ncol = 2)

pdf("20220114_ER_3accessions_RSA_analysis_Fig1.pdf", width = 13, height = 10)
plot_grid(MR_time_graph, LRL_time_graph, LRno_time_graph, labels = c("AUTO"), ncol = 2)
dev.off()
## png 
##   2
library(cowplot)
plot_grid(MR.delta_p_geno, aLRL.delta_p_geno, LRno.delta_p_geno, labels = c("AUTO"), ncol = 2)

pdf("20220114_ER_3accessions_RSA_analysis_Fig2.pdf", width = 13, height = 10)
plot_grid(MR.delta_p_geno, aLRL.delta_p_geno, LRno.delta_p_geno, labels = c("AUTO"), ncol = 2)
dev.off()
## png 
##   2
library(cowplot)
plot_grid(MR.STI_plot, aLRL.STI_plot, LRno.STI_plot, labels = c("AUTO"), ncol = 2)

pdf("20220114_ER_3accessions_RSA_analysis_Fig3.pdf", width = 13, height = 10)
plot_grid(MR.STI_plot, aLRL.STI_plot, LRno.STI_plot, labels = c("AUTO"), ncol = 2)
dev.off()
## png 
##   2
pdf("20230128_ER_3accessions_TRS_Fig4.pdf", width = 13, height = 10)
plot_grid(TRS_time_graph, aLRL.STI_plot, LRno.STI_plot, labels = c("AUTO"), ncol = 2)
dev.off()
## png 
##   2

When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).

The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.