This is a RSA analysis for 3 tomatoes experienced GA 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 undergraduate student Stacey Na in Dec.2022-January2023.

getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20230201RSA_M248M058LA1511_GA_100mM_Salt_Stacey"
list.files(pattern = ".csv")
## [1] "2023Feb_3acc_GA_growth_factors.csv" "d5.csv"                            
## [3] "d9-updated.csv"                     "d9.csv"
d5<-read.csv("d5.csv")
d9<-read.csv("d9-updated.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
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
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]  4 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] 0.4969276
angle_LR.avg
## [1] 83.18485
angle_LR.sd
## [1] 20.85216
LR_number
## [1] 4
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] 241
for(i in 2:241){
  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] 241  17
dim(MR_data_noChild)
## [1] 132  17
dim(MR_all)
## [1] 373  17
MR_all$root_name[50]
## [1] " pl2-c-higa-M058"
text <- strsplit(x = MR_all$root_name[50], split = "-")
text
## [[1]]
## [1] " pl2" "c"    "higa" "M058"
plate <- text[[1]][1]
plate <- gsub(" ", "", plate)
plate
## [1] "pl2"
cond <- text[[1]][2:3]
cond
## [1] "c"    "higa"
genotype <- text[[1]][4]
genotype
## [1] "M058"
dim(MR_all)
## [1] 373  17
for(i in 1:373){
  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
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] 232
unique(MR_all$day)
## [1] 5 9
MR_all$day <- as.numeric(as.character(MR_all$day))
unique(MR_all$day)
## [1] 5 9
373/232
## [1] 1.607759
MR_all$genotype <- gsub( "la1511", "LA1511", MR_all$genotype)
MR_all$root_name <- gsub( "la1511", "LA1511", MR_all$root_name)


good_stuff <- c("M248", "M058", "LA1511")
funk <- subset(MR_all, !(MR_all$genotype %in% good_stuff))
funk
unique(MR_all$genotype)
## [1] "M248"   "M058"   "LA1511" NA
MR_all <- subset(MR_all, (MR_all$genotype %in% good_stuff))

This is the same scenario I had for ABA…..more than 2 days!!!!! here I have less than 2 days

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
unique(MR_all$cond)
## [1] "c_noga" "s_noga" "c_loga" "s_loga" "c_higa" "s_higa"
MR_all2 <- subset(MR_all, MR_all$condition != "NA_NA")


library(ggplot2)
library(ggpubr)
## 
## Attaching package: 'ggpubr'
## The following object is masked from 'package:cowplot':
## 
##     get_legend
MR_all2$condition<- factor(MR_all2$condition, levels=c("c_noga", "c_loga", "c_higa",  "s_noga", "s_loga", "s_higa"))


histogram_TRS_GA <- ggdensity(MR_all2, x = "TRS",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition")
                           

histogram_TRS_GA

histogram_LRno_GA <- ggdensity(MR_all2, x = "LRno",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition")
                           

histogram_LRno_GA

pdf("histogram.TRS.GA.pdf")
plot(histogram_TRS_GA)
# if plotting multiple graphs - this command is extremely important 
dev.off()
## png 
##   2
pdf("histogram.LRno.GA.pdf")
plot(histogram_LRno_GA)
# if plotting multiple graphs - this command is extremely important 
dev.off()
## png 
##   2
library(cowplot)
pdf("Figure_Histogram.pdf", height = 15, width = 12)
plot_grid(histogram_TRS_GA, histogram_LRno_GA, ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## png 
##   2
unique(MR_all2$genotype)
## [1] "M248"   "M058"   "LA1511"

#check wether this peice is correct to use instead of gsub, if yes…what can I add to this peice to get actual genotype name, instead of TRUE……..

MR_all2$genotype <- gsub( "la1511", "LA1511", MR_all2$genotype)

#MR_all2$genotype<- "la1511" != "LA1511"
unique(MR_all2$genotype)
## [1] "M248"   "M058"   "LA1511"
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("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_all2$day)
## [1] 5 9
TRS_data <- subset(MR_all2, MR_all2$day == 9)
TRS_data
TRS_graph_d9 <- ggerrorplot(TRS_data, y = "TRS", x = "condition", fill="condition", color="condition", 
                        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 = "c_noga", 
                                                              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_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 = "LA1511", hide.ns = TRUE) 
better_TRS_graph

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_all2, aes(x= day, y=length, group = root_name, color = condition)) 
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)") 
ggplotly(my_graph)
my_graph <- ggerrorplot(MR_all2, y = "length", x = "condition", fill="cond", color="condition", 
                        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 = "c_noga", 
                                                              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

MR_all2$root_name <- gsub(" ", "", MR_all2$root_name)
geno_1511 <- subset(MR_all2, MR_all2$genotype == "LA1511")

bye <- c("pl2-s-loga-LA1511", "pl4-s-loga-LA1511", "pl3-s-higa-LA1511", "pl8-s-higa-LA1511")
geno_1511 <- subset(geno_1511, !(geno_1511$root_name %in% bye))
unique(geno_1511$root_name)
##  [1] "pl4-s-higa-LA1511"  "pl10-s-higa-LA1511" "pl5-c-higa-LA1511" 
##  [4] "pl1-s-noga-LA1511"  "pl1-s-higa-LA1511"  "pl2-s-higa-LA1511" 
##  [7] "pl5-s-higa-LA1511"  "pl6-s-higa-LA1511"  "pl7-s-higa-LA1511" 
## [10] "pl9-s-higa-LA1511"  "pl11-s-higa-LA1511" "pl1-c-higa-LA1511" 
## [13] "pl1-c-loga-LA1511"  "pl2-c-loga-LA1511"  "pl3-c-loga-LA1511" 
## [16] "pl4-c-loga-LA1511"  "pl5-c-loga-LA1511"  "pl6-c-loga-LA1511" 
## [19] "pl7-c-loga-LA1511"  "pl8-c-loga-LA1511"  "pl9-c-loga-LA1511" 
## [22] "pl10-c-loga-LA1511" "pl2-c-higa-LA1511"  "pl3-c-higa-LA1511" 
## [25] "pl4-c-higa-LA1511"  "pl6-c-higa-LA1511"  "pl7-c-higa-LA1511" 
## [28] "pl8-c-higa-LA1511"  "pl9-c-higa-LA1511"  "pl10-c-higa-LA1511"
## [31] "pl11-c-higa-LA1511" "pl1-s-loga-LA1511"  "pl3-s-loga-LA1511" 
## [34] "pl5-s-loga-LA1511"  "pl6-s-loga-LA1511"  "pl7-s-loga-LA1511" 
## [37] "pl8-s-loga-LA1511"  "pl9-s-loga-LA1511"  "pl10-s-loga-LA1511"
## [40] "pl11-s-loga-LA1511" "pl1-c-noga-LA1511"  "pl2-c-noga-LA1511" 
## [43] "pl3-c-noga-LA1511"  "pl4-c-noga-LA1511"  "pl5-c-noga-LA1511" 
## [46] "pl6-c-noga-LA1511"  "pl7-c-noga-LA1511"  "pl8-c-noga-LA1511" 
## [49] "pl9-c-noga-LA1511"  "pl10-c-noga-LA1511" "pl2-s-noga-LA1511" 
## [52] "pl3-s-noga-LA1511"  "pl4-s-noga-LA1511"  "pl5-s-noga-LA1511" 
## [55] "pl6-s-noga-LA1511"  "pl7-s-noga-LA1511"  "pl8-s-noga-LA1511" 
## [58] "pl9-s-noga-LA1511"  "pl10-s-noga-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(~ condition) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length (cm)") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

######try to remove outliers here after determining them through ggplotly

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

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

bye <- c("pl5-c-loga-M248", "pl6-c-loga-M248", "pl10-c-higa-M248", "pl8-s-noga-M248", "pl1-s-higa-M248")
geno_248 <- subset(geno_248, !(geno_248$root_name %in% bye))
unique(geno_248$root_name)
##  [1] "pl1-c-noga-M248"  "pl2-c-noga-M248"  "pl3-c-noga-M248"  "pl5-c-noga-M248" 
##  [5] "pl6-c-noga-M248"  "pl7-c-noga-M248"  "pl8-c-noga-M248"  "pl10-c-noga-M248"
##  [9] "pl2-s-noga-M248"  "pl3-s-noga-M248"  "pl4-s-noga-M248"  "pl5-s-noga-M248" 
## [13] "pl6-s-noga-M248"  "pl7-s-noga-M248"  "pl9-s-noga-M248"  "pl1-c-loga-M248" 
## [17] "pl2-c-loga-M248"  "pl3-c-loga-M248"  "pl4-c-loga-M248"  "pl7-c-loga-M248" 
## [21] "pl8-c-loga-M248"  "pl9-c-loga-M248"  "pl10-c-loga-M248" "pl1-s-loga-M248" 
## [25] "pl2-s-loga-M248"  "pl3-s-loga-M248"  "pl4-s-loga-M248"  "pl5-s-loga-M248" 
## [29] "pl6-s-loga-M248"  "pl7-s-loga-M248"  "pl8-s-loga-M248"  "pl9-s-loga-M248" 
## [33] "pl10-s-loga-M248" "pl1-c-higa-M248"  "pl2-c-higa-M248"  "pl3-c-higa-M248" 
## [37] "pl4-c-higa-M248"  "pl5-c-higa-M248"  "pl6-c-higa-M248"  "pl9-c-higa-M248" 
## [41] "pl2-s-higa-M248"  "pl4-s-higa-M248"  "pl6-s-higa-M248"  "pl4-c-noga-M248" 
## [45] "pl9-c-noga-M248"  "pl1-s-noga-M248"  "pl10-s-noga-M248" "pl7-c-higa-M248" 
## [49] "pl8-c-higa-M248"  "pl3-s-higa-M248"  "pl5-s-higa-M248"  "pl7-s-higa-M248" 
## [53] "pl8-s-higa-M248"  "pl9-s-higa-M248"  "pl10-s-higa-M248"
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(~ condition) 
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(~ 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=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(~ condition) 
my_graph <- my_graph + ylab("number") + xlab("time (day)") + ggtitle("Lateral Root number") + theme(legend.position='none')
my_graph

bye <- c("pl3-s-higa-M248", "pl6-s-higa-M248")
geno_248 <- subset(geno_248, !(geno_248$root_name %in% bye))
unique(geno_248$root_name)
##  [1] "pl1-c-noga-M248"  "pl2-c-noga-M248"  "pl3-c-noga-M248"  "pl5-c-noga-M248" 
##  [5] "pl6-c-noga-M248"  "pl7-c-noga-M248"  "pl8-c-noga-M248"  "pl10-c-noga-M248"
##  [9] "pl2-s-noga-M248"  "pl3-s-noga-M248"  "pl4-s-noga-M248"  "pl5-s-noga-M248" 
## [13] "pl6-s-noga-M248"  "pl7-s-noga-M248"  "pl9-s-noga-M248"  "pl1-c-loga-M248" 
## [17] "pl2-c-loga-M248"  "pl3-c-loga-M248"  "pl4-c-loga-M248"  "pl7-c-loga-M248" 
## [21] "pl8-c-loga-M248"  "pl9-c-loga-M248"  "pl10-c-loga-M248" "pl1-s-loga-M248" 
## [25] "pl2-s-loga-M248"  "pl3-s-loga-M248"  "pl4-s-loga-M248"  "pl5-s-loga-M248" 
## [29] "pl6-s-loga-M248"  "pl7-s-loga-M248"  "pl8-s-loga-M248"  "pl9-s-loga-M248" 
## [33] "pl10-s-loga-M248" "pl1-c-higa-M248"  "pl2-c-higa-M248"  "pl3-c-higa-M248" 
## [37] "pl4-c-higa-M248"  "pl5-c-higa-M248"  "pl6-c-higa-M248"  "pl9-c-higa-M248" 
## [41] "pl2-s-higa-M248"  "pl4-s-higa-M248"  "pl4-c-noga-M248"  "pl9-c-noga-M248" 
## [45] "pl1-s-noga-M248"  "pl10-s-noga-M248" "pl7-c-higa-M248"  "pl8-c-higa-M248" 
## [49] "pl5-s-higa-M248"  "pl7-s-higa-M248"  "pl8-s-higa-M248"  "pl9-s-higa-M248" 
## [53] "pl10-s-higa-M248"
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(~ condition) 
my_graph <- my_graph + ylab("length (cm)") + xlab("time (day)") + ggtitle("Lateral Root Length") + theme(legend.position='none')
my_graph

ggplotly(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(~ condition) 
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(~ 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=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(~ condition) 
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(~ condition) 
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_all3 <- rbind(geno_1511, geno_248)
MR_all3 <- rbind(MR_all3, geno_058)

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

aLRL_time_graph <- ggplot(data=MR_all3, aes(x= day, y=aLRL, group = root_name, color = condition)) 
aLRL_time_graph <- aLRL_time_graph + geom_line(alpha = 0.2) 
aLRL_time_graph <- aLRL_time_graph + facet_grid(genotype ~ condition) 
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_all3, aes(x= day, y=TRS, group = root_name, color = condition)) 
TRS_time_graph <- TRS_time_graph + geom_line(alpha = 0.2) 
TRS_time_graph <- TRS_time_graph + facet_grid(genotype ~ condition) 
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_all3, MR_all3$root_name == unique(MR_all3$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  
##       -1.212         1.113
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)    -1.212        NaN     NaN      NaN
## temp_MR2$day    1.113        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  
##          -24             5
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)       -24        NaN     NaN      NaN
## LR_temp2$day        5        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] 5
 #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.6351        0.1583
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.6351        NaN     NaN      NaN
## LR_temp2$day   0.1583        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.1582848
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$condition[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_all3$root_name))
## [1] 172

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:172){

  temp1 <- subset(MR_all3, MR_all3$root_name == unique(MR_all3$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$condition[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("c_noga", "c_loga", "c_higa",  "s_noga", "s_loga", "s_higa"))
growth_factors$genotype <- factor(growth_factors$genotype, levels = c("LA1511", "M058", "M248"))


write.csv(growth_factors, "2023Feb_3acc_GA_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 = "c_noga", 
                                                              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 = "c_noga", 
                                                              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 = "c_noga", 
                                                              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.

library(doBy)
## 
## Attaching package: 'doBy'
## The following object is masked from 'package:dplyr':
## 
##     order_by
unique(growth_factors$cond.final)
## [1] s_higa c_higa s_noga c_loga s_loga c_noga
## Levels: c_noga c_loga c_higa s_noga s_loga s_higa
avg_growth <- summaryBy(data = growth_factors, MR.delta + aLRL.delta + LRno.delta ~ genotype + cond.final + ga)
avg_growth
control_treatments <- c("c_noga", "c_loga", "c_higa")
avg_growth2 <- subset(avg_growth, avg_growth$cond.final %in% control_treatments)
colnames(avg_growth2) <- gsub(".mean", ".control", colnames(avg_growth2))
avg_growth2$ga <- avg_growth2$cond.final
avg_growth2$ga <- gsub("c_", "", avg_growth2$ga)

avg_growth2 <- avg_growth2[,c(1,3:6)]
avg_growth2
growth_factors2 <- growth_factors
growth_factors2$ga <- growth_factors2$cond.final
growth_factors2$ga <- gsub("c_", "", growth_factors2$ga)
growth_factors2$ga <- gsub("s_", "", growth_factors2$ga)
unique(growth_factors2$ga)
## [1] "higa" "noga" "loga"
unique(growth_factors2$cond.final)
## [1] s_higa c_higa s_noga c_loga s_loga c_noga
## Levels: c_noga c_loga c_higa s_noga s_loga s_higa
growth_factors3 <- merge(growth_factors2, avg_growth2, by = c("genotype", "ga"))
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$ga <- factor(growth_factors3$ga, levels=c("noga", "loga", "higa"))
### growth factor calculations - Salt Stress Tolerance Indexes graphs....

MR.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = ga, y = MR.STI, colour = ga)) 
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 = "c_noga", 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 = ga, y = aLRL.STI, colour = ga)) 
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 = "c_noga", 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 = ga, y = LRno.STI, colour = ga)) 
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 = "c_noga", 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