20230504_RSA_M248M058LA1511_ACC_100mM_Salt

This is a RSA analysis for 3 tomatoes experienced ACC 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 April 2023.

getwd()

## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20230504_RSA_M248M058LA1511_ACC_100mM_Salt-Stacey"

list.files(pattern = ".csv")

## [1] "2023May_3acc_ACC_growth_factors.csv" "d5-ACC-nodes.csv"                   
## [3] "d5-ACC.csv"                          "d9-ACC-nodes.csv"                   
## [5] "d9-ACC.csv"

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

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]  1 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.09437667

angle_LR.avg

## [1] 64.34508

angle_LR.sd

## [1] NA

LR_number

## [1] 1

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] 207

for(i in 2:207){
  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] 207  17

dim(MR_data_noChild)

## [1] 129  17

dim(MR_all)

## [1] 336  17

MR_all$root_name[50]

## [1] " pl2-s-5acc-la1511"

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

## [[1]]
## [1] " pl2"   "s"      "5acc"   "la1511"

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

## [1] "pl2"

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

## [1] "s"    "5acc"

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

## [1] "la1511"

dim(MR_all)

## [1] 336  17

for(i in 1:336){
  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] 162

unique(MR_all$day)

## [1] 5 9

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

## [1] 5 9

336/162

## [1] 2.074074

unique(MR_all$root_name)

##   [1] " pl7-s-5acc-la1511"  " pl1-s-5acc-m248"    " pl2-s-5acc-m248"   
##   [4] " pl3-s-5acc-m248"    " pl4-s-5acc-m248"    " pl5-s-5acc-m248"   
##   [7] " pl6-s-5acc-m248"    " pl7-s-5acc-m248"    " pl8-s-5acc-m248"   
##  [10] " pl9-s-5acc-m248"    " pl10-s-5acc-m248"   " pl2-s-5acc-m058"   
##  [13] " pl3-c-5acc-m248"    " pl4-s-noacc-m248"   " pl5-s-noacc-m248"  
##  [16] " pl6-s-noacc-m248"   " pl7-s-noacc-m248"   " pl8-s-noacc-m248"  
##  [19] " pl9-s-noacc-m248"   " pl10-s-noacc-m248"  " pl1-c-noacc-m248"  
##  [22] " pl3-c-noacc-m248"   " pl4-c-noacc-m248"   " pl5-c-noacc-m248"  
##  [25] " pl6-c-noacc-m248"   " pl7-c-noacc-m248"   " pl10-c-noacc-m248" 
##  [28] " pl2-c-noacc-m058"   " pl2-s-1acc-la1511"  " pl2-s-1acc-m248"   
##  [31] " pl3-s-1acc-m248"    " pl4-s-1acc-m248"    " pl5-s-1acc-m248"   
##  [34] " pl6-s-1acc-m248"    " pl7-s-1acc-m248"    " pl8-s-1acc-m248"   
##  [37] " pl9-s-1acc-m248"    " pl10-s-1acc-m248"   " pl6-s-1acc-m058"   
##  [40] " pl1-c-1acc-m248"    " pl3-c-1acc-m248"    " pl4-c-1acc-m248"   
##  [43] " pl5-c-1acc-m248"    " pl6-c-1acc-m248"    " pl7-c-1acc-m248"   
##  [46] " pl9-c-1acc-m248"    " pl10-c-1acc-m248"   " pl9-c-1acc-m058"   
##  [49] " pl1-s-5acc-la1511"  " pl2-s-5acc-la1511"  " pl3-s-5acc-la1511" 
##  [52] " pl4-s-5acc-la1511"  " pl6-s-5acc-la1511"  " pl1-s-5acc-m058"   
##  [55] " pl3-s-5acc-m058"    " pl4-s-5acc-m058"    " pl5-s-5acc-m058"   
##  [58] " pl6-s-5acc-m058"    " pl7-s-5acc-m058"    " pl8-5acc-s-M058"   
##  [61] " pl9-s-5acc-m058"    " pl10-s-5acc-m058"   " pl1-c-5acc-la1511" 
##  [64] " pl2-c-5acc-la1511"  " pl3-c-5acc-la1511"  " pl4-c-5acc-la1511" 
##  [67] " pl5-c-5acc-la1511"  " pl6-c-5acc-la1511"  " pl7-c-5acc-la1511" 
##  [70] " pl1-c-5acc-m248"    " pl2-c-5acc-m248"    " pl4-c-5acc-m248"   
##  [73] " pl5-c-5acc-m248"    " pl6-c-5acc-m248"    " pl7-c-5acc-m248"   
##  [76] " pl8-c-5acc-m248"    " pl9-c-5acc-m248"    " pl10-c-5acc-m248"  
##  [79] " pl1-c-5acc-m058"    " pl2-c-5acc-m058"    " pl3-c-5acc-m058"   
##  [82] " pl4-c-5acc-m058"    " pl5-c-5acc-m058"    " pl6-c-5acc-m058"   
##  [85] " pl7-c-5acc-m058"    " pl8-c-5acc-m058"    " pl9-c-5acc-m058"   
##  [88] " pl10-c-5acc-m058"   " pl1-s-noacc-la1511" " pl2-s-noacc-la1511"
##  [91] " pl3-s-noacc-la1511" " pl4-s-noacc-la1511" " pl5-s-noacc-la1511"
##  [94] " pl6-s-noacc-la1511" " pl7-s-noacc-la1511" " pl1-s-noacc-m248"  
##  [97] " pl2-s-noacc-m248"   " pl3-s-noacc-m248"   " pl1-s-noacc-m058"  
## [100] " pl2-s-noacc-m058"   " pl3-s-noacc-m058"   " pl4-s-noacc-m058"  
## [103] " pl5-s-noacc-m058"   " pl6-s-noacc-m058"   " pl7-s-noacc-m058"  
## [106] " pl8-s-noacc-m058"   " pl9-s-noacc-m058"   " pl10-s-noacc-m058" 
## [109] " pl1-c-noacc-la1511" " pl2-c-noacc-la1511" " pl3-c-noacc-la1511"
## [112] " pl4-c-noacc-la1511" " pl5-c-noacc-la1511" " pl6-c-noacc-la1511"
## [115] " pl7-c-noacc-la1511" " pl2-c-noacc-m248"   " pl8-c-noacc-m248"  
## [118] " pl9-c-noacc-m248"   " pl1-c-noacc-m058"   " pl3-c-noacc-m058"  
## [121] " pl4-c-noacc-m058"   " pl5-c-noacc-m058"   " pl6-c-noacc-m058"  
## [124] " pl7-c-noacc-m058"   " pl8-c-noacc-m058"   " pl9-c-noacc-m058"  
## [127] " pl10-c-noacc-m058"  " pl1-s-1acc-la1511"  " pl3-s-1acc-la1511" 
## [130] " pl4-s-1acc-la1511"  " pl5-s-1acc-la1511"  " pl6-s-1acc-la1511" 
## [133] " pl1-s-1acc-m058"    " pl2-s-1acc-m058"    " pl3-s-1acc-m058"   
## [136] " pl4-s-1acc-m058"    " pl5-s-1acc-m058"    " pl7-s-1acc-m058"   
## [139] " pl8-s-1acc-m058"    " pl9-s-1acc-m058"    " pl10-s-1acc-m058"  
## [142] " pl2-c-1acc-la1511"  " pl3-c-1acc-la1511"  " pl4-c-1acc-la1511" 
## [145] " pl5-c-1acc-la1511"  " pl6-c-1acc-la1511"  " pl7-c-1acc-la1511" 
## [148] " pl2-c-1acc-m248"    " pl8-c-1acc-m248"    " pl1-c-1acc-m058"   
## [151] " pl2-c-1acc-m058"    " pl3-c-1acc-m058"    " pl4-c-1acc-m058"   
## [154] " pl5-c-1acc-m058"    " pl6-c-1acc-m058"    " pl7-c-1acc-m058"   
## [157] " pl8-c-1acc-m058"    " pl5-s-5acc-la1511"  " pl8-s-5acc-m058"   
## [160] " root_1"             " pl1-c-1acc-la1511"  " root_0"

unique(MR_all$genotype)

## [1] "la1511" "m248"   "m058"   "M058"   NA

MR_all$genotype <- gsub( "M058", "m058", MR_all$genotype)
MR_all$root_name <- gsub( "M058", "m058", MR_all$root_name)

MR_all$condition <- gsub( "5acc_s", "s_5acc", MR_all$condition)
MR_all$root_name <- gsub( "5acc_s", "s_5acc", MR_all$root_name)

good_stuff <- c("m248", "m058", "la1511")
funk <- subset(MR_all, !(MR_all$genotype %in% good_stuff))
funk

good_stuff2<- c("s_5acc",  "c_5acc",  "s_noacc", "c_noacc" ,"s_1acc" , "c_1acc")
funk2 <- subset(MR_all, !(MR_all$condition %in% good_stuff2))
funk

unique(MR_all$genotype)

## [1] "la1511" "m248"   "m058"   NA

unique(MR_all$condition)

## [1] "s_5acc"  "c_5acc"  "s_noacc" "c_noacc" "s_1acc"  "c_1acc"  "NA_NA"

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$condition)

## [1] "s_5acc"  "c_5acc"  "s_noacc" "c_noacc" "s_1acc"  "c_1acc"  "NA_NA"

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_noacc","c_1acc", "c_5acc", "s_noacc", "s_1acc", "s_5acc"))


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

histogram_TRS_ACC

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

histogram_LRno_ACC

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

## png 
##   2

pdf("histogram.LRno.ACC.pdf")
plot(histogram_LRno_ACC)
# 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_ACC, histogram_LRno_ACC, ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()

## png 
##   2

unique(MR_all2$genotype)

## [1] "la1511" "m248"   "m058"

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_noacc", 
                                                              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)

## Warning: package 'ggbeeswarm' was built under R version 4.3.3

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.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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_noacc", 
                                                              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

geno_1511 <- subset(MR_all2, MR_all2$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(~ 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_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

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)

bye <- c("pl10-c-1acc-m248", "pl10-c-5acc-m248", "pl9-c-5acc-m248", "pl3-s-5acc-m248")
geno_248 <- subset(geno_248, !(geno_248$root_name %in% bye))

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

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

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

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

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.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

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

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.6398        0.4267

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.6398        NaN     NaN      NaN
## temp_MR2$day   0.4267        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  
##       -22.75          4.75

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)    -22.75        NaN     NaN      NaN
## LR_temp2$day     4.75        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] 4.75

 #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.6221        0.1433

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.6221        NaN     NaN      NaN
## LR_temp2$day   0.1433        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.1432954

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] 160

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

  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_noacc","c_1acc", "c_5acc", "s_noacc", "s_1acc", "s_5acc"))
growth_factors$genotype <- factor(growth_factors$genotype, levels = c("la1511", "m058", "m248"))


write.csv(growth_factors, "2023May_3acc_ACC_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_noacc", 
                                                              label = "p.signif", hide.ns = T)
MR.delta_p_geno

#TO EXTRACT control and salt (CS) data only

control_salt_only <- subset(growth_factors, cond.final %in% c("c_noacc","s_noacc"))


MR.delta_p_geno_CS <- ggerrorplot(control_salt_only, y = "MR.delta", x = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=9,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Main root growth rate (cm / day)") 
MR.delta_p_geno_CS <- MR.delta_p_geno_CS + rremove("legend") + stat_compare_means(method="t.test", ref.group = "c_noacc", 
                                                              label = "p.signif", hide.ns = T)
MR.delta_p_geno_CS

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_noacc", 
                                                              label = "p.signif", hide.ns = T)
LRno.delta_p_geno <- LRno.delta_p_geno #+ ggtitle("Lateral Root Number Increase")
LRno.delta_p_geno

#TO EXTRACT control and salt (CS) data only

LRno.delta_p_geno_CS <- ggerrorplot(control_salt_only, 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_CS <- LRno.delta_p_geno_CS + rremove("legend") + stat_compare_means(method="t.test", ref.group = "c_noacc", 
                                                              label = "p.signif", hide.ns = T)
LRno.delta_p_geno_CS <- LRno.delta_p_geno_CS #+ ggtitle("Lateral Root Number Increase")
LRno.delta_p_geno_CS

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_noacc", 
                                                              label = "p.signif", hide.ns = T)
aLRL.delta_p_geno <- aLRL.delta_p_geno #+ ggtitle("average Lateral Root Growth")
aLRL.delta_p_geno

#TO EXTRACT control and salt (CS) data only

aLRL.delta_p_geno_CS <- ggerrorplot(control_salt_only, 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_CS <- aLRL.delta_p_geno_CS + rremove("legend") + stat_compare_means(method="t.test", ref.group = "c_noacc", 
                                                              label = "p.signif", hide.ns = T)
aLRL.delta_p_geno_CS <- aLRL.delta_p_geno_CS #+ ggtitle("average Lateral Root Growth")
aLRL.delta_p_geno_CS

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_5acc  s_1acc  c_5acc  s_noacc c_noacc c_1acc 
## Levels: c_noacc c_1acc c_5acc s_noacc s_1acc s_5acc

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

control_treatments <- c("c_noacc", "c_1acc", "c_5acc")
avg_growth2 <- subset(avg_growth, avg_growth$cond.final %in% control_treatments)
colnames(avg_growth2) <- gsub(".mean", ".control", colnames(avg_growth2))
avg_growth2$acc <- avg_growth2$cond.final
avg_growth2$acc <- gsub("c_", "", avg_growth2$acc)

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

growth_factors2 <- growth_factors
growth_factors2$acc <- growth_factors2$cond.final
growth_factors2$acc <- gsub("c_", "", growth_factors2$acc)
growth_factors2$acc <- gsub("s_", "", growth_factors2$acc)
unique(growth_factors2$acc)

## [1] "5acc"  "1acc"  "noacc"

unique(growth_factors2$cond.final)

## [1] s_5acc  s_1acc  c_5acc  s_noacc c_noacc c_1acc 
## Levels: c_noacc c_1acc c_5acc s_noacc s_1acc s_5acc

growth_factors3 <- merge(growth_factors2, avg_growth2, by = c("genotype", "acc"))
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$acc <- factor(growth_factors3$acc, levels=c("noacc", "1acc", "5acc"))

### growth factor calculations - Salt Stress Tolerance Indexes graphs....

growth_factors3 <- subset(growth_factors3, growth_factors3$MR.STI > 0.5& MR.STI < 1.25 )

MR.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = acc, y = MR.STI, colour = acc)) 
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 = "noacc") #hide.ns = T) 
MR.STI_plot <- MR.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
MR.STI_plot

growth_factors3 <- subset(growth_factors3, growth_factors3$aLRL.STI > 0.25& aLRL.STI < 1.5 )

aLRL.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = acc, y = aLRL.STI, colour = acc)) 
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 = "noacc")#, hide.ns = T) 
aLRL.STI_plot <- aLRL.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
aLRL.STI_plot

growth_factors3 <- subset(growth_factors3, growth_factors3$LRno.STI > 0.25& LRno.STI < 1.6 )

LRno.STI_plot <- ggplot(data = growth_factors3, mapping = aes(x = acc, y = LRno.STI, colour = acc)) 
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 = "noacc")#, hide.ns = T) 
LRno.STI_plot <- LRno.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
LRno.STI_plot

pdf("Figure_STI.pdf", height = 5, width = 12)
plot_grid(MR.STI_plot, aLRL.STI_plot, LRno.STI_plot, ncol=3,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()

## png 
##   2

#TO EXTRACT control and salt (CS) data only

pdf("Figure_control_salt_only.pdf", height = 5, width = 12)
plot_grid(MR.delta_p_geno_CS,aLRL.delta_p_geno_CS, LRno.delta_p_geno_CS, ncol=3,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()

## png 
##   2