Same as RSA for duf160-5 as Maryam did before.Seeds were sterilized with 50% bleach for 10 min, rinsed 5 times with MQ sterile water, and incubated at 4°C for 24 h. ½ ms media + sucrose (including vitamins) was made for both control and salt with 75 and 125 mM NaCl. Salt added after PH adjustment.The seeds were germinated for 4 complete days at control plate, and then at d5, they were transferred to +/- salt plates, and scanned for 5 times, every other day. and then the samples were harvested after 2 weeks on salt plate for root and shoot FW. The Arabidopsis growth chamber (#34 at BTI) had 21 °C and 16 h light/8 h dark cycle, with 60% humidity.

getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/202201_DUF160-4_RSA_traced_by_HS"
list.files(pattern = ".csv")
## [1] "202310_duf-4_growth_factors.csv" "Maryam_DUF_roots.csv"           
## [3] "Maryam_DUF_roots_D11.csv"        "Maryam_DUF_roots_D13.csv"       
## [5] "Maryam_DUF_roots_D5.csv"         "Maryam_DUF_roots_D7.csv"        
## [7] "Maryam_DUF_roots_D9.csv"
my_data_d5 <-  read.csv("Maryam_DUF_roots_D5.csv" )
my_data_d7 <- read.csv("Maryam_DUF_roots_D7.csv")
my_data_d9 <- read.csv("Maryam_DUF_roots_D9.csv" )
my_data_d11 <- read.csv("Maryam_DUF_roots_D11.csv" )
my_data_d13 <- read.csv("Maryam_DUF_roots_D13.csv" )

head(my_data_d5)
my_data_d5$day <- 5
my_data_d7$day <- 7
my_data_d9$day <- 9
my_data_d11$day <- 11
my_data_d13$day <- 13

head(my_data_d13)
all_data <- rbind(my_data_d5, my_data_d7)
all_data <- rbind(all_data, my_data_d9)
all_data <- rbind(all_data, my_data_d11)
all_data <- rbind(all_data, my_data_d13)
all_data
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"
Main_root[,c(3:5,10)]
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_2 <- Main_root[,c(1:3,4,16,17,19:22)]
MR_data_2
Lateral_root <- subset(all_data, all_data$root_order == 1)
Lateral_root
unique(Lateral_root$root_order)
## [1] 1
MR_data_2$root[1]
## [1] " e53862cc-ed8f-4f24-9807-a50fc85f64b9"
temporary <- subset(Lateral_root, Lateral_root$parent == MR_data_2$root[1])
temporary
MR_data_3 <- subset(MR_data_2, MR_data_2$n_child >0)
temporary <- subset(Lateral_root, Lateral_root$parent == MR_data_3$root[1])
temporary <- subset(temporary, temporary$day == MR_data_3$day[1])
temporary
dim(temporary)
## [1]  1 22
dim(temporary)[1]
## [1] 1
dim(temporary)[2]
## [1] 22
total_LRL <- sum(temporary$length)
LR_number <- dim(temporary)[1]
total_LRL
## [1] 0.04473915
LR_number
## [1] 1
MR_data_3$LRL <- 0
MR_data_3$LRno <- 0
MR_data_3
MR_data_3$LRL[1] <- total_LRL
MR_data_3$LRno[1] <- LR_number
MR_data_3
MR_data_noChild <- subset(MR_data_2, MR_data_2$n_child == 0)
MR_data_noChild

lets loop

dim(MR_data_3)
## [1] 236  12
for(i in 1:236){
  temporary <- subset(Lateral_root, Lateral_root$parent == MR_data_3$root[i])
  temporary <- subset(temporary, temporary$day == MR_data_3$day[i])
  total_LRL <- sum(temporary$length)
  LR_number <- dim(temporary)[1]
  MR_data_3$LRL[i] <- total_LRL
  MR_data_3$LRno[i] <- LR_number
}

let have a look at data now

MR_data_3$check <- MR_data_3$n_child - MR_data_3$LRno
MR_data_3
unique(MR_data_3$check)
## [1] 0
#MR_all <- rbind(MR_data_3,MR_data_noChild)
head(MR_data_3)
#let's remove check collumn from Child:

MR_data_Child2 <- MR_data_3[,1:12]
MR_data_Child2
MR_data_noChild
MR_data_noChild$LRL <- 0
MR_data_noChild$LRno <- 0


MR_all <- rbind(MR_data_Child2, MR_data_noChild)
unique(MR_all$day)
## [1]  7  9 11 13  5
?strsplit()
## starting httpd help server ... done
MR_all$root_name[1]
## [1] " pl5-0-col2"
text <- strsplit(x = MR_all$root_name[1], split = "-")
text
## [[1]]
## [1] " pl5" "0"    "col2"
plate <- text[[1]][1]
plate
## [1] " pl5"
plate <- gsub(" ", "", plate)
plate
## [1] "pl5"
cond <- text[[1]][2]
cond
## [1] "0"
genotype <- text[[1]][3]
genotype
## [1] "col2"
#OK - now let's loop it:

dim(MR_all)
## [1] 613  12
for(i in 1:613){
  text <- strsplit(x = MR_all$root_name[i], split = "-")
  plate <- text[[1]][1]
  cond <- text[[1]][2]
  genotype <- text[[1]][3]
   
  MR_all$genotype[i] <- genotype
  MR_all$condition[i] <- cond
  MR_all$plate[i] <- plate
}

MR_all
unique(MR_all$genotype)
## [1] "col2" "duf2" "duf1" "col1" NA
MR_all.nona <- na.omit(MR_all)
unique(MR_all.nona$genotype)
## [1] "col2" "duf2" "duf1" "col1"
MR_all.nona$genotype <- gsub("col1", "Col-0", MR_all.nona$genotype)
MR_all.nona$genotype <- gsub("col2", "Col-0", MR_all.nona$genotype)
MR_all.nona$genotype <- gsub("duf1", "duf-4", MR_all.nona$genotype)
MR_all.nona$genotype <- gsub("duf2", "duf-4", MR_all.nona$genotype)
unique(MR_all.nona$genotype)
## [1] "Col-0" "duf-4"
unique(MR_all.nona$condition)
## [1] "0"   "75"  "125"
dim(MR_all.nona)
## [1] 597  15
unique(MR_all$day)
## [1]  7  9 11 13  5
MR_all.nona$TRS <- MR_all.nona$length + MR_all.nona$LRL
MR_all.nona$aLRL <- MR_all.nona$LRL/ MR_all.nona$LRno
MR_all.nona$MRpLRL <- MR_all.nona$length / MR_all.nona$LRL
MR_all.nona
library(ggplot2)
library(ggpubr)

MR_all.nona$condition<- factor(MR_all.nona$condition, levels=c("0","75", "125"))


histogram_TRS <- ggdensity(MR_all.nona, x = "TRS",
                           add = "mean", rug = TRUE,facet.by = "day",
                           color = "condition", fill = "condition",
                           palette = c("#00AFBB", "#E7B800", "#FF0000"))
histogram_TRS

histogram_LRL <- ggdensity(MR_all.nona, x = "LRL",
                           add = "mean", rug = TRUE,facet.by = "day",
                           color = "condition", fill = "condition",
                           palette = c("#00AFBB", "#E7B800", "#FF0000"))
histogram_LRL

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

library(ggsci)
library(ggbeeswarm)
library(gapminder)
library(RColorBrewer)
library(ggridges)
better_TRS_graph <- ggplot(data=MR_all.nona, aes(x= genotype, y=TRS, color = condition))
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) + 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 = "Col-0") 
better_TRS_graph
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 5, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 6, 3
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 5, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 6, 3
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 5, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 3, 6
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 5, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 6, 3
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! argument "x" is missing, with no default
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 5, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 6, 3

#install.packages("plotly")
library(ggplot2)
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
head(MR_all)
my_graph <- ggplot(data=MR_all.nona, aes(x= day, y=length, group = root_name, color = genotype))
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ condition) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph

ggplotly(my_graph)

####now I want to inspect each particular genotype to see whether all replicates behaving well….basically curating the data!I start with col first.

geno_c <- subset(MR_all.nona, MR_all.nona$genotype == "Col-0")

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

ggplotly(my_graph)
bye <- c("pl3-0-col1", "pl9-75-col2", "pl1-125-col2", "pl1-125-col1")
geno_c <- subset(geno_c, !(geno_c$root_name %in% bye))
ggplotly(my_graph)

now calculate TRS and LRno as well

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

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

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

####NOW do the same thing for the duf genotype

geno_duf <- subset(MR_all.nona, MR_all.nona$genotype == "duf-4")
my_graph <- ggplot(data=geno_duf, 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_duf, 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_duf, 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_duf, 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

####lets combine all curated data from both genotypes together in one graph

MR_both_geno <- rbind(geno_c, geno_duf)

MR_time_graph <- ggplot(data=MR_both_geno, aes(x= day, y=length, group = root_name, color = genotype)) 
MR_time_graph <- MR_time_graph + geom_line(alpha = 0.2) 
MR_time_graph <- MR_time_graph + facet_grid(~ condition)
MR_time_graph <- MR_time_graph + scale_color_manual(values= c("blue", "red"))
MR_time_graph <- MR_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("Main Root Length")
MR_time_graph <- MR_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## 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_both_geno, aes(x= day, y=LRL, group = root_name, color = genotype)) 
LRL_time_graph <- LRL_time_graph + geom_line(alpha = 0.2) 
LRL_time_graph <- LRL_time_graph + facet_grid(~ condition) 
LRL_time_graph <- LRL_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("Lateral Root Length")
LRL_time_graph <- LRL_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
LRL_time_graph

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

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

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

### Growth calculations

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

temp1 <- subset(MR_both_geno, MR_both_geno$root_name == unique(MR_both_geno$root_name)[1])
temp1
temp2 <- temp1[order(temp1$day),]
temp2
temp2$condition[1]
## [1] 0
## Levels: 0 75 125
# 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)

so in this case - we should remove day 8 and day 9 MR, but we won;t have time to look at individual pictures, and thus let’s make it into a logical removal loop:

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
plot(temp_MR$length~ temp_MR$day)
# let's add the regression line to this graph
abline(lm(temp_MR$length~ temp_MR$day))

MR_model <- lm(temp_MR$length~ temp_MR$day)
MR_model
## 
## Call:
## lm(formula = temp_MR$length ~ temp_MR$day)
## 
## Coefficients:
## (Intercept)  temp_MR$day  
##      -4.287        1.044
summary(MR_model)
## 
## Call:
## lm(formula = temp_MR$length ~ temp_MR$day)
## 
## Residuals:
##       1       2       3       4       5 
##  0.2234 -0.2690 -0.2390  0.3914 -0.1067 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.28733    0.50510  -8.488 0.003434 ** 
## temp_MR$day  1.04402    0.05354  19.499 0.000295 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3386 on 3 degrees of freedom
## Multiple R-squared:  0.9922, Adjusted R-squared:  0.9896 
## F-statistic: 380.2 on 1 and 3 DF,  p-value: 0.0002946
MR_model$coefficients[[2]]
## [1] 1.044018
MR_growth_rate <- MR_model$coefficients[[2]]

Cool - now we have to do this for LRno and aLRL, which are much easier, because we just need to remove all the ‘n.a.’

# Let's remove all NA for LRno and aLRL in temp2
temp2
LR_temp <- temp2[,c("day", "LRno", "aLRL")]
LR_temp2 <- na.omit(LR_temp)
LR_temp2
# 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  
##       -33.75          4.55
summary(LRno_model)
## 
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
## 
## Residuals:
##    1    2    3    4 
##  2.9 -3.2 -2.3  2.6 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -33.7500     8.9771  -3.760   0.0640 .
## LR_temp2$day   4.5500     0.8761   5.194   0.0351 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.918 on 2 degrees of freedom
## Multiple R-squared:  0.931,  Adjusted R-squared:  0.8965 
## F-statistic: 26.97 on 1 and 2 DF,  p-value: 0.03513
LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase
## [1] 4.55

Now let’s move on to aLRL

# 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.8933        0.1311
summary(aLRL_model)
## 
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
## 
## Residuals:
##         1         2         3         4 
##  0.019980 -0.007045 -0.045850  0.032915 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  -0.893253   0.097675  -9.145  0.01175 * 
## LR_temp2$day  0.131145   0.009532  13.758  0.00524 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.04263 on 2 degrees of freedom
## Multiple R-squared:  0.9895, Adjusted R-squared:  0.9843 
## F-statistic: 189.3 on 1 and 2 DF,  p-value: 0.005241
aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth
## [1] 0.1311447

OK - so now let’s make a new empty table to save all of the data from these calculations

temp2$condition[1]
## [1] 0
## Levels: 0 75 125
temp2$condition[1][1]
## [1] 0
## Levels: 0 75 125
names <- c(text="root_name", "genotype", "condition", "MR.delta", "LRno.delta", "aLRL.delta")
growth_factors <- data.frame()
growth_factors
for (k in names) growth_factors[[k]] <- as.character()

growth_factors[1,1] <- temp2$root_name[1]
growth_factors[1,2] <- temp2$genotype[1]
growth_factors[1,3] <- temp2$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

Ok - before we loop it - lets see how many plants we have in our dataset:

length(unique(MR_both_geno$root_name))
## [1] 120

Now - lets loop this loop

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

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

Now - let’s subset for all plants that actually have growth:

growth_factors <- subset(growth_factors, growth_factors$MR.delta > 0)
growth_factors
unique(growth_factors$condition)
## [1] "1" "3" "2"
growth_factors$condition <- gsub("1", "0", growth_factors$condition)
growth_factors$condition <- gsub("2", "75", growth_factors$condition)
growth_factors$condition <- gsub("3", "125", growth_factors$condition)
growth_factors$condition <- as.factor(growth_factors$condition)
growth_factors
write.csv(growth_factors, "202310_duf-4_growth_factors.csv", row.names = FALSE)

Visualizing the growth factors:

growth_factors2 <- subset(growth_factors, growth_factors$MR.delta >= 0)

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

growth_factors$condition <- factor(growth_factors$condition, levels = c("0", "75", "125"))

MR.delta_p_geno <- ggerrorplot(growth_factors2, y = "MR.delta", x = "genotype", fill="genotype", color="genotype", 
                        facet.by = c("condition"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Growth Rate (cm / day)") 
MR.delta_p_geno <- MR.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Col-0", 
                                                              label = "p.signif")
MR.delta_p_geno <- MR.delta_p_geno + ggtitle("Main Root Growth")
MR.delta_p_geno
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 53, 45, 39, 40, 47, 50, 51, 52, 48, 55, 49, 54, 41, 43, 44, 46, 42, 38, 15, 106, 98, 92, 99, 104, 109, 107, 108, 103, 102, 100, 105, 96, 94, 97, 93, 83, 101, 95, 91
## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `mutate()`:
## ℹ In argument: `p = purrr::map(...)`.
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## ℹ With name: x.1.
## Caused by error in `t.test.default()`:
## ! not enough 'x' observations
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 27, 37, 26, 28, 24, 29, 34, 30, 32, 35, 25, 31, 36, 19, 20, 33, 23, 22, 18, 21, 73, 82, 85, 89, 88, 86, 84, 80, 90, 87, 78, 79, 81, 71, 76, 75, 77, 72, 74
## Warning: Removed 39 rows containing missing values (`geom_segment()`).
## Warning: Removed 31 rows containing missing values (`geom_segment()`).
## Warning: Removed 39 rows containing missing values (`geom_segment()`).

growth_factors2 <- subset(growth_factors, growth_factors$LRno.delta >= 0)

LRno.delta_p_geno <- ggerrorplot(growth_factors2, y = "LRno.delta", x = "genotype", fill="genotype", color="genotype", 
                        facet.by = c("condition"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="LR no increase (#LR / day)") 
LRno.delta_p_geno <- LRno.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Col-0", 
                                                              label = "p.signif")
LRno.delta_p_geno <- LRno.delta_p_geno + ggtitle("Lateral Root Increase")
LRno.delta_p_geno

growth_factors2 <- subset(growth_factors, growth_factors$aLRL.delta >= 0)

aLRL.delta_p_geno2 <- ggerrorplot(growth_factors2, y = "aLRL.delta", x = "genotype", fill="genotype", color="genotype", 
                        facet.by = c("condition"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Lateral Root Growth Rate (cm / day)") 
aLRL.delta_p_geno2 <- aLRL.delta_p_geno2 + rremove("legend") + stat_compare_means(method="t.test", ref.group = "Col-0", 
                                                              label = "p.signif")
aLRL.delta_p_geno2

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

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

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

library(doBy)
avg_growth <- summaryBy(data = growth_factors2, MR.delta + aLRL.delta + LRno.delta ~ genotype + condition)
avg_growth

Then - we subset the data into only control condition, and merge the average control with the growth_factors in all conditions:

avg_growth_C <- subset(avg_growth, avg_growth$condition == "0")
avg_growth_C
avg_growth_C <- avg_growth_C[,c(1,3:5)]
avg_growth_C
colnames(avg_growth_C) <- gsub(".mean", ".avg.Control", colnames(avg_growth_C))
avg_growth_C

now let’s merge it into the growth_factors:

STI_growth_factors <- merge(growth_factors2, avg_growth_C, id="genotype")
STI_growth_factors

Now - let’s calculate Salt Tolerance Indexes (STI) by dividing individual growth rates by their avg.Control:

STI_growth_factors$MR.STI <- STI_growth_factors$MR.delta / STI_growth_factors$MR.delta.avg.Control
STI_growth_factors$aLRL.STI <- STI_growth_factors$aLRL.delta / STI_growth_factors$aLRL.delta.avg.Control
STI_growth_factors$LRno.STI <- STI_growth_factors$LRno.delta / STI_growth_factors$LRno.delta.avg.Control
head(STI_growth_factors)
library(reshape2)
STI <- subset(STI_growth_factors, STI_growth_factors$cond == "125")
STI125 <- STI[,c(1,10:12)]
STI125
mSTI125 <- melt(STI125)
## Using genotype as id variables
colnames(mSTI125)[3] <- "fraction.of.control"
mSTI125
unique(mSTI125$variable)
## [1] MR.STI   aLRL.STI LRno.STI
## Levels: MR.STI aLRL.STI LRno.STI

Cool - now let’s visualize this thing for MR:

STI125_plot <- ggplot(data = mSTI125, aes(x = genotype, y = fraction.of.control, colour = genotype)) 
STI125_plot <- STI125_plot + geom_beeswarm(alpha=0.6, priority = "density")
STI125_plot <- STI125_plot + facet_wrap(~ variable, ncol = 3)
STI125_plot <- STI125_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
STI125_plot <- STI125_plot + theme(legend.position = "none")
STI125_plot <- STI125_plot + xlab("")
STI125_plot <- STI125_plot + ylab("Fraction of control") + ggtitle("Salt Tolerance Index for roots grown at 125 mM NaCl")
STI125_plot <- STI125_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "Col-0") 
STI125_plot <- STI125_plot + scale_color_manual(values=c("royalblue", "coral3"))
STI125_plot

STI <- subset(STI_growth_factors, STI_growth_factors$cond == "75")
STI75 <- STI[,c(1,10:12)]
STI75
mSTI75 <- melt(STI75)
## Using genotype as id variables
colnames(mSTI75)[3] <- "fraction.of.control"
mSTI75
unique(mSTI75$variable)
## [1] MR.STI   aLRL.STI LRno.STI
## Levels: MR.STI aLRL.STI LRno.STI
STI75_plot <- ggplot(data = mSTI75, aes(x = genotype, y = fraction.of.control, colour = genotype)) 
STI75_plot <- STI75_plot + geom_beeswarm(alpha=0.6, priority = "density")
STI75_plot <- STI75_plot + facet_wrap(~ variable, ncol = 3)
STI75_plot <- STI75_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
STI75_plot <- STI75_plot + theme(legend.position = "none")
STI75_plot <- STI75_plot + xlab("")
STI75_plot <- STI75_plot + ylab("Fraction of control") + ggtitle("Salt Tolerance Index for roots grown at 75 mM NaCl")
STI75_plot <- STI75_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "Col-0") 
STI75_plot <- STI75_plot + scale_color_manual(values=c("royalblue", "coral3"))
STI75_plot

library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
pdf("202310_duf-4_RSA_Growth_factors_STI_Graphs.pdf", width = 13, height = 13)
plot_grid(MR.delta_p_geno , aLRL.delta_p_geno2, LRno.delta_p_geno, STI75_plot, STI125_plot, labels = c("AUTO"), ncol = 2)
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 53, 45, 39, 40, 47, 50, 51, 52, 48, 55, 49, 54, 41, 43, 44, 46, 42, 38, 15, 106, 98, 92, 99, 104, 109, 107, 108, 103, 102, 100, 105, 96, 94, 97, 93, 83, 101, 95, 91
## Warning: Computation failed in `stat_compare_means()`
## Caused by error in `mutate()`:
## ℹ In argument: `p = purrr::map(...)`.
## Caused by error in `purrr::map()`:
## ℹ In index: 1.
## ℹ With name: x.1.
## Caused by error in `t.test.default()`:
## ! not enough 'x' observations
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 27, 37, 26, 28, 24, 29, 34, 30, 32, 35, 25, 31, 36, 19, 20, 33, 23, 22, 18, 21, 73, 82, 85, 89, 88, 86, 84, 80, 90, 87, 78, 79, 81, 71, 76, 75, 77, 72, 74
## Warning: Removed 39 rows containing missing values (`geom_segment()`).
## Warning: Removed 31 rows containing missing values (`geom_segment()`).
## Warning: Removed 39 rows containing missing values (`geom_segment()`).
dev.off()
## png 
##   2