20220824_RSA_M248M058LA1511_ABA_Salt_2days_only

This is a RSA calculation for 3 tomatoes experienced ABA and salt treatments on plate. The seeds were germinated on 1/4 ms plate for 4 complete days, and then transferred to treatment plates at d5, and being scanned for 5 consecutive days. The root tracing was done by REU student YALMARIE NUMAN-VAZQUEZ in summer 2022.

getwd()

## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20220809_RSA_M248M058LA1511_ABA_100mM_Salt"

list.files(pattern = ".csv")

## [1] "2022Nov_3acc_ABA_growth_factors.csv" "d5-ABA2.csv"                        
## [3] "d9-ABA2.csv"

d5<-read.csv("d5-ABA2.csv")
d9<-read.csv("d9-ABA2.csv")
head(d5)

head(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)
library(cowplot)
d5<-d5%>%mutate(day=5,d5 )
head(d5)

d9<-d9%>%mutate(day=9,d9 )
head(d9)

all_data<-rbind(d5,d9)
all_data

unique(all_data$day)

## [1] 5 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, 14: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 am not sure what is the purpose of the 2 following temporary files? 
temporary <- subset(Lateral_root, Lateral_root$parent == MR_data2$root[1])
temporary

temporary <- subset(temporary, temporary$day == MR_data2$day[3])
temporary

dim(temporary)

## [1]  8 22

total_LRL <- sum(temporary$length)
angle_LR.avg <- mean(temporary$insertion_angle)
angle_LR.sd <- sd(temporary$insertion_angle)
LR_number <- dim(temporary)[1]
total_LRL

## [1] 4.751957

angle_LR.avg

## [1] 65.94405

angle_LR.sd

## [1] 34.4859

LR_number

## [1] 8

MR_data2$LRL <- 0
MR_data2$LRno <- 0
MR_data2$LRangle.avg <- 0
MR_data2$LRangle.sd <- 0
MR_data2

MR_data2$LRL[1] <- total_LRL
MR_data2$LRno[1] <- LR_number
MR_data2$LRangle.avg[1] <- angle_LR.avg
MR_data2$LRangle.sd[1] <- angle_LR.sd
MR_data2

MR_data_noChild <- subset(MR_data, MR_data$n_child == 0)
MR_data_noChild

length(MR_data2$root)

## [1] 167

for(i in 2:167){
  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)
  angle_LR.avg <- mean(temporary$insertion_angle)
  angle_LR.sd <- sd(temporary$insertion_angle)
  LR_number <- dim(temporary)[1]
  MR_data2$LRL[i] <- total_LRL
  MR_data2$LRno[i] <- LR_number
  MR_data2$LRangle.avg[i] <- angle_LR.avg
  MR_data2$LRangle.sd[i] <- angle_LR.sd
    
}

MR_data2$check <- MR_data2$n_child - MR_data2$LRno
MR_data2

unique(MR_data2$check)

## [1] 0

colnames(MR_data2)

##  [1] "image"                 "root_name"             "root"                 
##  [4] "length"                "root_order"            "insertion_position"   
##  [7] "insertion_angle"       "n_child"               "child_density"        
## [10] "insertion_first_child" "last_child"            "insertion_last_child" 
## [13] "day"                   "LRL"                   "LRno"                 
## [16] "LRangle.avg"           "LRangle.sd"            "check"

MR_data_Child2 <- MR_data2[,1:17]
MR_data_Child2

MR_data_noChild

MR_data_noChild$LRL <- 0
MR_data_noChild$LRno <- 0
MR_data_noChild$LRangle.avg <- 0
MR_data_noChild$LRangle.sd <- 0
MR_all <-rbind(MR_data_Child2, MR_data_noChild)
MR_all

dim(MR_data_Child2)

## [1] 167  17

dim(MR_data_noChild)

## [1] 176  17

dim(MR_all)

## [1] 343  17

MR_all$root_name[50]

## [1] " pl1-1aba-s-M058"

text <- strsplit(x = MR_all$root_name[50], split = "-")
text

## [[1]]
## [1] " pl1" "1aba" "s"    "M058"

plate <- text[[1]][1]
plate <- gsub(" ", "", plate)
plate

## [1] "pl1"

cond <- text[[1]][2:3]
cond

## [1] "1aba" "s"

genotype <- text[[1]][4]
genotype

## [1] "M058"

dim(MR_all)

## [1] 343  17

#I guess here I am having issue with number of components in root naming , 3 vs 4???  
for(i in 1:343){
  text <- strsplit(x = MR_all$root_name[i], split = "-")
  plate <- text[[1]][1]
  cond <- paste(text[[1]][2], text[[1]][3], sep="_")
  genotype <- text[[1]][4]
   
  MR_all$genotype[i] <- genotype
  MR_all$condition[i] <- cond
  MR_all$plate[i] <- plate
}

MR_all

It’s a bit of a mess - let’s look at unique condition - so we can fix that:

unique(MR_all$condition)

##  [1] "10aba_M248"   "10aba_M058"   "10aba_s"      "1aba_s"       "1aba_M248"   
##  [6] "1aba_M058"    "s_M058"       "s_M248"       "c_M248"       "c_M058"      
## [11] "10aba_LA1511" "1aba_LA1511"  "s_LA1511"     "c_LA1511"

MR_all$condition1 <- "EMPTY"
MR_all$condition2 <- "EMPTY"


for(i in 1:343){
  temp <- strsplit(MR_all$condition[i], split="_")
  MR_all$condition1[i] <- temp[[1]][1]
  MR_all$condition2[i] <- temp[[1]][2]
}

unique(MR_all$condition1)

## [1] "10aba" "1aba"  "s"     "c"

unique(MR_all$condition2)

## [1] "M248"   "M058"   "s"      "LA1511"

genotypes <- c("M248", "M058", "LA1511")
conditions <- c("10aba", "1aba", "s", "c")

MR_all$cond_new2 <- "EMPTY"

for(i in 1:343){
  
  if(MR_all$condition2[i] %in% conditions){
    MR_all$cond_new2[i] <- MR_all$condition2[i]
  }
  if(MR_all$condition2[i] %in% genotypes){
    MR_all$genotype[i] <- MR_all$condition2[i]
  }
  
}
i

## [1] 343

MR_all$condition2[i]

## [1] "s"

MR_all

MR_all$cond <- paste(MR_all$condition1, MR_all$cond_new2, sep="_")
unique(MR_all$cond)

## [1] "10aba_EMPTY" "10aba_s"     "1aba_s"      "1aba_EMPTY"  "s_EMPTY"    
## [6] "c_EMPTY"

MR_all$cond <- gsub("10aba_EMPTY", "10aba_c", MR_all$cond)
MR_all$cond <- gsub("1aba_EMPTY", "1aba_c", MR_all$cond)
MR_all$cond <- gsub("s_EMPTY", "noaba_s", MR_all$cond)
MR_all$cond <- gsub("c_EMPTY", "noaba_c", MR_all$cond)
MR_all$cond <- factor(MR_all$cond, levels = c("noaba_c", "1aba_c", "10aba_c",  "noaba_s", "1aba_s", "10aba_s"))
unique(MR_all$cond)

## [1] 10aba_c 10aba_s 1aba_s  1aba_c  noaba_s noaba_c
## Levels: noaba_c 1aba_c 10aba_c noaba_s 1aba_s 10aba_s

MR_all$TRS <- MR_all$length + MR_all$LRL

for(i in 1:nrow(MR_all)){
  if(MR_all$LRno[i] == 0){
    MR_all$aLRL[i] <- 0
  } else {
  MR_all$aLRL[i] <- MR_all$LRL[i]/ MR_all$LRno[i]  
  }
}


MR_all$MRpLRL <- MR_all$length / MR_all$LRL
MR_all

length(unique(MR_all$root_name))

## [1] 169

unique(MR_all$day)

## [1] 9 5

MR_all$day <- as.numeric(as.character(MR_all$day))
unique(MR_all$day)

## [1] 9 5

MR_d9 <- subset(MR_all, MR_all$day == 5)
MR_d5 <- subset(MR_all, MR_all$day == 9)

suspish <- subset(MR_d9, !(MR_d9$root_name %in% MR_d5$root_name))
suspish

MR_clean <- subset(MR_all, (MR_all$root_name %in% MR_d5$root_name))

MR_d9 <- subset(MR_clean, MR_clean$day == 5)
MR_d5 <- subset(MR_clean, MR_clean$day == 9)

suspish <- subset(MR_d9, !(MR_d9$root_name %in% MR_d5$root_name))
suspish

suspish <- subset(MR_d9, !(MR_d5$root_name %in% MR_d9$root_name))
suspish

library(ggplot2)
library(ggpubr)

## 
## Attaching package: 'ggpubr'

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

MR_all$cond<- factor(MR_all$cond, levels=c("noaba_c", "1aba_c", "10aba_c",  "noaba_s", "1aba_s", "10aba_s"))

unique(MR_all$cond)

## [1] 10aba_c 10aba_s 1aba_s  1aba_c  noaba_s noaba_c
## Levels: noaba_c 1aba_c 10aba_c noaba_s 1aba_s 10aba_s

histogram_TRS_ABA <- ggdensity(MR_all, x = "TRS",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "cond", fill = "cond")
                           

histogram_TRS_ABA

histogram_LRno_ABA <- ggdensity(MR_all, x = "LRno",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "cond", fill = "cond")
                           #palette = c("#00AFBB", "#E7B800","#FC4E07", "#FF0000"))

histogram_LRno_ABA

pdf("histogram.TRS.ABA.pdf")
plot(histogram_TRS_ABA)
# if plotting multiple graphs - this command is extremely important 
dev.off()

## png 
##   2

pdf("histogram.LRno.ABA.pdf")
plot(histogram_LRno_ABA)
# if plotting multiple graphs - this command is extremely important 
dev.off()

## png 
##   2

library(cowplot)

pdf("Figure_MAIN_1.pdf", height = 15, width = 12)
plot_grid(histogram_TRS_ABA, histogram_LRno_ABA, ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()

## png 
##   2

unique(MR_all$genotype)

## [1] "M248"   "M058"   "LA1511"

TRS_lgraph <- ggplot(data=MR_all, aes(x= genotype, y=TRS, color = cond)) 
TRS_lgraph <- TRS_lgraph + geom_boxplot()
TRS_lgraph <- TRS_lgraph + facet_grid(day ~ cond, scales = "free") + scale_color_manual(values=c("blue","turquoise3", "cyan","red", "maroon3", "dark orange"))
TRS_lgraph <- TRS_lgraph + ylab("Total root size (cm)") + xlab("Genotype") + theme(legend.position='none')
TRS_lgraph

unique(MR_all$day)

## [1] 9 5

TRS_data <- subset(MR_all, MR_all$day == 9)
TRS_data

TRS_data$cond <- factor(TRS_data$cond, levels = c("noaba_c", "1aba_c", "10aba_c",  "noaba_s", "1aba_s", "10aba_s"))
TRS_data$genotype <- factor(TRS_data$genotype, levels = c("LA1511", "M058", "M248"))

TRS_graph_d9 <- ggerrorplot(TRS_data, y = "TRS", x = "cond", fill="cond", color="cond", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "genotype", palette = "jco",
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Total Root Size (cm)") 
TRS_graph_d9 <- TRS_graph_d9 + rremove("legend") + stat_compare_means(method="t.test", ref.group = "noaba_c", 
                                                              label = "p.signif", hide.ns = T)
TRS_graph_d9 <- TRS_graph_d9 + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
TRS_graph_d9

library(ggsci)
library(ggbeeswarm)
library(gapminder)
library(RColorBrewer)
library(ggridges)
better_TRS_graph <- ggplot(data=MR_all, 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 ~ cond, 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 = "LA1511", hide.ns = TRUE) 
better_TRS_graph

### Increase in MRL / LRL / LRno through time - first look at the TRS graphs

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

my_graph <- ggplot(data=MR_all, aes(x= day, y=length, group = root_name, color = cond)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ genotype) 
my_graph <- my_graph + ylab("Main Root Length (cm)") + xlab("time (day)") 
my_graph

ggplotly(my_graph)

my_graph <- ggerrorplot(MR_all, y = "length", x = "cond", fill="cond", color="cond", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "genotype", palette = "jco",
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Main Root Length (cm)") 
my_graph <- my_graph + rremove("legend") + stat_compare_means(method="t.test", ref.group = "noaba_c", 
                                                              label = "p.signif", hide.ns = T)
my_graph <- my_graph + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
my_graph

### Curating the data

geno_1511 <- subset(MR_all, MR_all$genotype == "LA1511")

my_graph <- ggplot(data=geno_1511, 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(~ cond) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

my_graph <- ggplot(data=geno_1511, 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(~ cond) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Total Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_1511, 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(~ cond) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_1511, 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(~ cond) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

geno_248 <- subset(MR_all, MR_all$genotype == "M248")
geno_248

my_graph <- ggplot(data=geno_248, 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(~ cond) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

my_graph <- ggplot(data=geno_248, 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(~ cond) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Total Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_248, 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(~ cond) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_248, 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(~ cond) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

geno_058 <- subset(MR_all, MR_all$genotype == "M058")

my_graph <- ggplot(data=geno_058, 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(~ cond) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_058, 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(~ cond) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Total Root Length") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_058, 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(~ cond) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

my_graph <- ggplot(data=geno_058, 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(~ cond) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

Increase in MRL / LRL / LRno through time - final graphs

MR_all2 <- rbind(geno_1511, geno_248)
MR_all2 <- rbind(MR_all2, geno_058)

MR_time_graph <- ggplot(data=MR_all2, aes(x= day, y=length, group = root_name, color = cond)) 
MR_time_graph <- MR_time_graph + geom_line(alpha = 0.2) 
MR_time_graph <- MR_time_graph + facet_grid(~ genotype) 
MR_time_graph <- MR_time_graph + ylab("Main root length (cm)") + xlab("Time (days after germination)") #+ ggtitle("Main Root Length")
MR_time_graph <- MR_time_graph + stat_summary(fun=mean, aes(group= cond),  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

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

LRno_time_graph <- ggplot(data=MR_all2, aes(x= day, y=LRno, group = root_name, color = cond)) 
LRno_time_graph <- LRno_time_graph + geom_line(alpha = 0.2) 
LRno_time_graph <- LRno_time_graph + facet_grid(~ genotype) 
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= cond),  size=0.7, geom="line", linetype = "dashed")
LRno_time_graph

aLRL_time_graph <- ggplot(data=MR_all2, aes(x= day, y=aLRL, group = root_name, color = cond)) 
aLRL_time_graph <- aLRL_time_graph + geom_line(alpha = 0.2) 
aLRL_time_graph <- aLRL_time_graph + facet_grid(~ genotype) 
aLRL_time_graph <- aLRL_time_graph + ylab("Average lateral root 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= cond),  size=0.7, geom="line", linetype = "dashed")

aLRL_time_graph

TRS_time_graph <- ggplot(data=MR_all2, aes(x= day, y=TRS, group = root_name, color = cond)) 
TRS_time_graph <- TRS_time_graph + geom_line(alpha = 0.2) 
TRS_time_graph <- TRS_time_graph + facet_grid(~ genotype) 
TRS_time_graph <- TRS_time_graph + ylab("Total root length (cm)") + xlab("Time (days after germination)") #+ ggtitle("Total Root Length")
TRS_time_graph <- TRS_time_graph + stat_summary(fun.y=mean, aes(group= cond),  size=0.7, geom="line", linetype = "dashed")
TRS_time_graph

growth factor calculations - setup

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_all2, MR_all2$root_name == unique(MR_all2$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")
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  
##      3.39487       0.06053

summary(MR_model)

## 
## Call:
## lm(formula = temp_MR2$length ~ temp_MR2$day)
## 
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)   3.39487        NaN     NaN      NaN
## temp_MR2$day  0.06053        NaN     NaN      NaN
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 1 and 0 DF,  p-value: NA

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  
##        -1.25          0.25

summary(LRno_model)

## 
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
## 
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)     -1.25        NaN     NaN      NaN
## LR_temp2$day     0.25        NaN     NaN      NaN
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 1 and 0 DF,  p-value: NA

LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase

## [1] 0.25

 #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  
##     -0.38996       0.07799

summary(aLRL_model)

## 
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
## 
## Residuals:
## ALL 2 residuals are 0: no residual degrees of freedom!
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)  -0.38996        NaN     NaN      NaN
## LR_temp2$day  0.07799        NaN     NaN      NaN
## 
## Residual standard error: NaN on 0 degrees of freedom
## Multiple R-squared:      1,  Adjusted R-squared:    NaN 
## F-statistic:   NaN on 1 and 0 DF,  p-value: NA

aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth

## [1] 0.07799167

temp2

names <- c(text="root_name", "genotype", "cond.final", "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] <- as.character(temp2$genotype[1])
growth_factors[1,3] <- as.character(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_all2$root_name))

## [1] 169

growth factor calculations - loop

# 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:169){

  temp1 <- subset(MR_all2, MR_all2$root_name == unique(MR_all2$root_name)[e])
  temp2 <- temp1[order(temp1$day),]
  temp2

############ MR calculations
  temp_MR <- temp2[,c("day", "length")]
  if(dim(temp_MR)[1] > 1){
    MR_model <- lm(temp_MR$length~ temp_MR$day)
    MR_growth_rate <- MR_model$coefficients[[2]]} else{MR_growth_rate <- 0}
    
  ############ LRno calculations
  LR_temp <- temp2[,c("day", "LRno", "aLRL")]
  LR_temp$aLRL <- as.numeric(as.character(LR_temp$aLRL))
  #LR_temp2 <- na.omit(LR_temp)
  dim(LR_temp)
  
  ####################### safety precaution to calculate LR growth rate only for the plants that have LR at #least for two days: 
  if(dim(LR_temp)[1] > 1){
  LRno_model <- lm(LR_temp$LRno ~ LR_temp$day)
  LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))

  ############ aLRL calculations
  aLRL_model <- lm(LR_temp$aLRL ~ LR_temp$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] <- as.character(temp2$genotype[1])
  growth_factors[e,3] <- as.character(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 <- subset(growth_factors, growth_factors$MR.delta > 0)
growth_factors

growth_factors$cond.final <- factor(growth_factors$cond.final, levels = c("noaba_c", "1aba_c", "10aba_c",  "noaba_s", "1aba_s", "10aba_s"))
growth_factors$genotype <- factor(growth_factors$genotype, levels = c("LA1511", "M058", "M248"))


write.csv(growth_factors, "2022Nov_3acc_ABA_growth_factors.csv", row.names = FALSE)

growth factor calculations - graphs

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 = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Main root growth rate (cm / day)") 
MR.delta_p_geno <- MR.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "noaba_c", 
                                                              label = "p.signif", hide.ns = T)
MR.delta_p_geno

LRno.delta_p_geno <- ggerrorplot(growth_factors, y = "LRno.delta", x = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Lateral root number increase rate (# LR / day)") 
LRno.delta_p_geno <- LRno.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "noaba_c", 
                                                              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 = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Average lateral root length growth rate (cm / day)") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "noaba_c", 
                                                              label = "p.signif", hide.ns = T)
aLRL.delta_p_geno <- aLRL.delta_p_geno #+ ggtitle("average Lateral Root Growth")
aLRL.delta_p_geno

growth factor calculations - Salt Stress Tolerance Indexes calculations. in theory, I have 3 STI, noaba_s/noaba_c, 1aba_s/1aba_c, and 10aba_s/10aba_c.

library(doBy)

## 
## Attaching package: 'doBy'

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

unique(growth_factors$cond.final)

## [1] 10aba_s 10aba_c 1aba_s  1aba_c  noaba_s noaba_c
## Levels: noaba_c 1aba_c 10aba_c noaba_s 1aba_s 10aba_s

avg_growth <- summaryBy(data = growth_factors, MR.delta + aLRL.delta + LRno.delta ~ genotype + cond.final + aba)
avg_growth

control_treatments <- c("10aba_c", "1aba_c", "noaba_c")
avg_growth2 <- subset(avg_growth, avg_growth$cond.final %in% control_treatments)
colnames(avg_growth2) <- gsub(".mean", ".control", colnames(avg_growth2))
avg_growth2$aba <- avg_growth2$cond.final
avg_growth2$aba <- gsub("_c", "", avg_growth2$aba)

avg_growth2 <- avg_growth2[,c(1,3:6)]
avg_growth2

growth_factors2 <- growth_factors
growth_factors2$aba <- growth_factors2$cond.final
growth_factors2$aba <- gsub("_c", "", growth_factors2$aba)
growth_factors2$aba <- gsub("_s", "", growth_factors2$aba)
unique(growth_factors2$aba)

## [1] "10aba" "1aba"  "noaba"

unique(growth_factors2$cond.final)

## [1] 10aba_s 10aba_c 1aba_s  1aba_c  noaba_s noaba_c
## Levels: noaba_c 1aba_c 10aba_c noaba_s 1aba_s 10aba_s

growth_factors3 <- merge(growth_factors2, avg_growth2, by = c("genotype", "aba"))
growth_factors3

growth_factors3$MR.STI <- growth_factors3$MR.delta / growth_factors3$MR.delta.control
growth_factors3$aLRL.STI <- growth_factors3$aLRL.delta / growth_factors3$aLRL.delta.control
growth_factors3$LRno.STI <- growth_factors3$LRno.delta / growth_factors3$LRno.delta.control
head(growth_factors3)

growth_factors3$aba <- factor(growth_factors3$aba, levels=c("noaba", "1aba", "10aba"))

growth factor calculations - Salt Stress Tolerance Indexes graphs….

MR.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = aba, y = MR.STI, colour = aba)) 
MR.STI_plot <- MR.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
MR.STI_plot <- MR.STI_plot + facet_grid(~ genotype)
MR.STI_plot <- MR.STI_plot + stat_summary(fun=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("Salt tolerance index based on main root growth (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 = "noaba_c", hide.ns = T) 
MR.STI_plot <- MR.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
MR.STI_plot

## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `if (ref.group == ".all.") ...`:
## ! missing value where TRUE/FALSE needed

aLRL.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = aba, y = aLRL.STI, colour = aba)) 
aLRL.STI_plot <- aLRL.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
aLRL.STI_plot <- aLRL.STI_plot + facet_grid(~ genotype)
aLRL.STI_plot <- aLRL.STI_plot + stat_summary(fun=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("Salt tolerance index based on average lateral root growth (Fraction of control)")
aLRL.STI_plot <- aLRL.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "noaba_c", hide.ns = T) 
aLRL.STI_plot <- aLRL.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
aLRL.STI_plot

## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `if (ref.group == ".all.") ...`:
## ! missing value where TRUE/FALSE needed

LRno.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = aba, y = LRno.STI, colour = aba)) 
LRno.STI_plot <- LRno.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
LRno.STI_plot <- LRno.STI_plot + facet_grid(~ genotype)
LRno.STI_plot <- LRno.STI_plot + stat_summary(fun=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("Salt tolerance index based on average increase in lateral root number (Fraction of control)")
LRno.STI_plot <- LRno.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "noaba_c", hide.ns = T) 
LRno.STI_plot <- LRno.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
LRno.STI_plot

## Warning: Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Computation failed in `stat_compare_means()`
## Caused by error in `if (ref.group == ".all.") ...`:
## ! missing value where TRUE/FALSE needed

# to continue: I included traits such as insertion position and angle in the data to analyze. Do mean and median and SD and SE for that traits as well………………….

MR_all2_D9 <- subset(MR_all2, MR_all2$day == 9)
MR_all2_D9

MR_all2_D9$condition<- factor(MR_all2_D9$condition, levels=c("noaba_c", "1aba_c", "10aba_c","noaba_s","1aba_s","10aba_s"))

LRangle.sd_graph <- ggplot(data=MR_all2_D9, aes(x= cond, y=LRangle.sd, color = cond))
LRangle.sd_graph <- LRangle.sd_graph + geom_beeswarm(alpha=0.6, priority = "density")
LRangle.sd_graph <- LRangle.sd_graph + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
LRangle.sd_graph <- LRangle.sd_graph + facet_grid(~ genotype, scales = "free") + scale_color_manual(values=c("blue","turquoise3", "cyan","red", "maroon3", "dark orange"))
LRangle.sd_graph <- LRangle.sd_graph + ylab("Standard deviation for lateral root angle") + xlab("") + theme(legend.position='none')
LRangle.sd_graph <- LRangle.sd_graph + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
LRangle.sd_graph <- LRangle.sd_graph + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "noaba_c", hide.ns = TRUE) 
LRangle.sd_graph

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

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

## Warning: Removed 3 rows containing missing values (`position_beeswarm()`).

## Warning: Removed 5 rows containing missing values (`position_beeswarm()`).

LRangle.avg_graph <- ggplot(data=MR_all2_D9, aes(x= cond, y=LRangle.avg, color = cond))
LRangle.avg_graph <- LRangle.avg_graph + geom_beeswarm(alpha=0.6, priority = "density")
LRangle.avg_graph <- LRangle.avg_graph + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
LRangle.avg_graph <- LRangle.avg_graph + facet_grid(~ genotype, scales = "free") + scale_color_manual(values=c("blue","turquoise3", "cyan","red", "maroon3", "dark orange"))

LRangle.avg_graph <- LRangle.avg_graph + ylab("Average of lateral root angle per plant") + xlab("") + theme(legend.position='none')
LRangle.avg_graph <- LRangle.avg_graph + theme(axis.text.x = element_text(angle=90, hjust=0.9, vjust=0.5))
LRangle.avg_graph <- LRangle.avg_graph + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "noaba_c", hide.ns = TRUE) 
LRangle.avg_graph