Experimental description:

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

Genotype definition code: M248=a, M058=b, LA1511=c.

getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/8-2-21 RSA M248 058 LA1511_analyzed"
setwd("C:/Users/Julkowska Lab/Desktop/R codes by Maryam/8-2-21 RSA M248 058 LA1511_analyzed")
list.files(pattern = ".csv")
## [1] "20211215_M248M058LA1511_ICPMS_analysis_on_Plate_Ranalysis.csv"
## [2] "2021Feb_ER_growth_factors.csv"                                
## [3] "2021Nov_M248M058LA1511_growth_factors.csv"                    
## [4] "M248M058LA1511_d5.csv"                                        
## [5] "M248M058LA1511_d6.csv"                                        
## [6] "M248M058LA1511_d7.csv"                                        
## [7] "M248M058LA1511_d8.csv"                                        
## [8] "M248M058LA1511_d9.csv"
my_data_d5 <-  read.csv("M248M058LA1511_d5.csv")
my_data_d6 <-  read.csv("M248M058LA1511_d6.csv")
my_data_d7 <-  read.csv("M248M058LA1511_d7.csv")
my_data_d8 <-  read.csv("M248M058LA1511_d8.csv")
my_data_d9 <-  read.csv("M248M058LA1511_d9.csv")
head(my_data_d5)
head(my_data_d6)
head(my_data_d7)
head(my_data_d8)
head(my_data_d9)

Adding the day after germination

my_data_d5$day <- 5
my_data_d6$day <- 6
my_data_d7$day <- 7
my_data_d8$day <- 8
my_data_d9$day <- 9
head(my_data_d8)
head(my_data_d5)

fuse days together

all_data <- rbind(my_data_d5, my_data_d6)
all_data <- rbind(all_data, my_data_d7)
all_data <- rbind(all_data, my_data_d8)
all_data <- rbind(all_data, my_data_d9)
all_data
unique(all_data$day)
## [1] 5 6 7 8 9
Main_root <- subset(all_data, all_data$root_order == 0)
Main_root

now we want to navigate the columns in the fused file

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:4,10,16,17,19:22)]
MR_data_2
Lateral_root <- subset(all_data, all_data$root_order != 0)
Lateral_root
unique(Lateral_root$root_order)
## [1] 1 2
MR_data_3 <- subset(MR_data_2, MR_data_2$n_child >0)

Because you also have secondary LR in there - we should remove them from calculations of LR

temporary <- subset(Lateral_root, Lateral_root$parent == MR_data_3$root[1])
temporary

######after runing above code, the number of rows decreased from 1700 to 44, what does that mean, means I removed lots of secondary LR???? in my RSA…..11/24/21 In the above - we also have more LR because the same parent was traced over multiple days - and thus we should also match it with the MR_data_3$day collumn:

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

now this looks much better - let;s move on to calculating LRL and LRno

dim(temporary)
## [1]  5 22
dim(temporary)[1]
## [1] 5
dim(temporary)[2]
## [1] 22
total_LRL <- sum(temporary$length)
LR_number <- dim(temporary)[1]
total_LRL
## [1] 1.358813
LR_number
## [1] 5
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
length(MR_data_3$root)
## [1] 210
for(i in 2:210){
  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
}

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

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


MR_all <- rbind(MR_data_Child2, MR_data_noChild)
MR_all
?strsplit()
## starting httpd help server ... done
MR_all$root_name[1]
## [1] " pl1_c_a"
text <- strsplit(x = MR_all$root_name[1], split = "_")
text
## [[1]]
## [1] " pl1" "c"    "a"
plate <- text[[1]][1]
plate <- gsub(" ", "", plate)
plate
## [1] "pl1"
cond <- text[[1]][2]
cond
## [1] "c"
genotype <- text[[1]][3]
genotype
## [1] "a"

####what does this mean? “plate <- gsub(” “,”“, plate) plate”11-24-21

I continue with looping

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

MR_all
MR_all$TRS <- MR_all$length + MR_all$LRL
MR_all$aLRL <- MR_all$LRL/ MR_all$LRno
MR_all$MRpLRL <- MR_all$length / MR_all$LRL
MR_all

Let’s check few things:

length(unique(MR_all$root_name))
## [1] 60
unique(MR_all$day)
## [1] 6 7 8 9 5

Let’s make sure that our day is saved as character…..I

#####what does 19 mean?

MR_all$day <- as.numeric(as.character(MR_all$day))
unique(MR_all$day)
## [1] 6 7 8 9 5
# let;s check if all of the days are complete
# we have so much unique plants that we traced:
length(unique(MR_all$root_name))
## [1] 60
# and so much unique images:
dim(MR_all)
## [1] 300  19
# thus dividing images / plants should give us number of days:
300/60
## [1] 5
#Correct!!! ok - let's move on!

Data visualization

start visualizing

#install.packages("ggplot")
#install.packages("ggpubr")
library(ggplot2)
library(ggpubr)
histogram_TRS <- ggdensity(MR_all, x = "TRS",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition",
                           palette = c("#00AFBB", "#E7B800"))

histogram_TRS

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

histogram_LRno

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

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

from here I followed Magda’s “R analysis of salt induced changes in Root Architecture of tomato plants”in Rpubs…..I think I need to graph the total root length over the time for different genotypes……

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

ggplotly(my_graph)

Curating the data

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

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

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

ggplotly(my_graph)

Let’s additionally check the TRS and LRno and so on:

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

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

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

Cool - looks good - so let’s move to the next genotype:

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

ggplotly(my_graph)

Now - lets do the other plots to see whether we can find anything else wrong:

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

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

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

And finally let’s look at genotype C:

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

ggplotly(my_graph)

Now - let’s have a look at all the graphs:

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

### Visualizing the time series:

ok - all look well to me. So let’s combine all the curated datasets into one and let;s replot them into one big graph all together.

Pls note that at this stage it is too few samples to do stats on it

MR <- rbind(geno_a, geno_b)
MR <- rbind(MR, geno_c)

MR_time_graph <- ggplot(data=MR, aes(x= day, y=length, group = root_name, color = genotype)) 
MR_time_graph <- MR_time_graph + geom_line(alpha = 0.2) 
MR_time_graph <- MR_time_graph + facet_grid(~ condition) 
MR_time_graph <- MR_time_graph + ylab("Main root length (cm)") + xlab("Time (days after germination)") #+ ggtitle("Main Root Length")
MR_time_graph <- MR_time_graph + stat_summary(fun.y=mean, aes(group= genotype),  size=0.7, geom="line", linetype = "dashed")
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
MR_time_graph

LRL_time_graph <- ggplot(data=MR, 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("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= genotype),  size=0.7, geom="line", linetype = "dashed")
LRL_time_graph

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

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

TRS_time_graph <- ggplot(data=MR, 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

OK - but because we have the line b that is a slow grower anyways, it would be good to have the relative effect of salt, and its best to calculate it using the growth rates of individual organs. Let’s do that then!

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, MR$root_name == unique(MR$root_name)[1])
temp2 <- temp1[order(temp1$day),]
temp2
# For Main Root Growth Rate - we want to remove all the data points that are repeating, because that indicates root hitting the plate edge:
temp_MR <- temp2[,c("day", "length")]
plot(temp_MR$length~ temp_MR$day)

Although I do not have growth inhibition for the M248 M058 LA1511 in all days but I still would follow the following steps to make sure this is true for all plants……“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

OK - the above looks good - so now we have to subset into MR temp where MRdouble == no, and calculate the growth rate:

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

This looks good! now we just have to get the model paramerers out of the linear model (lm) that was used to draw the regression line:

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.177         0.735
summary(MR_model)
## 
## Call:
## lm(formula = temp_MR2$length ~ temp_MR2$day)
## 
## Residuals:
##        1        2        3        4        5 
## -0.04894  0.10583 -0.04897 -0.02377  0.01586 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -1.17659    0.16904   -6.96  0.00608 ** 
## temp_MR2$day  0.73503    0.02367   31.05 7.34e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.07485 on 3 degrees of freedom
## Multiple R-squared:  0.9969, Adjusted R-squared:  0.9959 
## F-statistic: 964.2 on 1 and 3 DF,  p-value: 7.338e-05
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
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  
##          -19             4
summary(LRno_model)
## 
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
## 
## Residuals:
##  1  2  3  4 
##  0 -1  2 -1 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  -19.0000     5.8737  -3.235   0.0837 .
## LR_temp2$day   4.0000     0.7746   5.164   0.0355 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.732 on 2 degrees of freedom
## Multiple R-squared:  0.9302, Adjusted R-squared:  0.8953 
## F-statistic: 26.67 on 1 and 2 DF,  p-value: 0.03551
LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase
## [1] 4

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  
##      -2.7275        0.4906
summary(aLRL_model)
## 
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
## 
## Residuals:
##         1         2         3         4 
##  0.055671 -0.005927 -0.155159  0.105415 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -2.7275     0.4694  -5.810   0.0284 *
## LR_temp2$day   0.4906     0.0619   7.925   0.0156 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1384 on 2 degrees of freedom
## Multiple R-squared:  0.9691, Adjusted R-squared:  0.9537 
## F-statistic: 62.81 on 1 and 2 DF,  p-value: 0.01555
aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth
## [1] 0.4905961

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

names <- c(text="root_name", "genotype", "cond", "MR.delta", "LRno.delta", "aLRL.delta")
growth_factors <- data.frame()

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

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

growth_factors

OK - the above data table looks good - let’s loop it for all the root_names in the MR_OK dataset:

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

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

growth_factors

There are SOME negative values in here - let’s remove them:

growth_factors <- subset(growth_factors, growth_factors$MR.delta > 0)
growth_factors <- subset(growth_factors, growth_factors$LRno.delta > 0)
growth_factors <- subset(growth_factors, growth_factors$aLRL.delta > 0)
dim(growth_factors)
## [1] 54  6

Looks good - now let’s save this file and have a look at how the growth factos compare between stress and genotypes

write.csv(growth_factors, "2021Nov_M248M058LA1511_growth_factors.csv", row.names = FALSE)

Visualizing the growth factors:

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


library(multcompView)
growth_factors$GenoCond <- paste(growth_factors$cond, "_", growth_factors$genotype, sep = "")
growth_factors
Output <- TukeyHSD(aov(MR.delta ~ GenoCond, data = growth_factors))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = MR.delta ~ GenoCond, data = growth_factors)
## 
## $GenoCond
##               diff        lwr          upr     p adj
## c_b-c_a -0.5310910 -1.1696289  0.107446791 0.1541215
## c_c-c_a  0.2749453 -0.3465618  0.896452369 0.7764969
## s_a-c_a -0.6129007 -1.2514385  0.025637092 0.0666480
## s_b-c_a -0.7627821 -1.4804367 -0.045127568 0.0312112
## s_c-c_a -0.3517869 -0.9732940  0.269720146 0.5513705
## c_c-c_b  0.8060363  0.1674985  1.444574143 0.0060504
## s_a-c_b -0.0818097 -0.7369357  0.573316282 0.9990283
## s_b-c_b -0.2316911 -0.9641442  0.500762020 0.9342980
## s_c-c_b  0.1793041 -0.4592337  0.817841920 0.9597402
## s_a-c_c -0.8878460 -1.5263838 -0.249308200 0.0019183
## s_b-c_c -1.0377274 -1.7553820 -0.320072860 0.0011445
## s_c-c_c -0.6267322 -1.2482393 -0.005225146 0.0470154
## s_b-s_a -0.1498814 -0.8823345  0.582571719 0.9900028
## s_c-s_a  0.2611138 -0.3774240  0.899651619 0.8280799
## s_c-s_b  0.4109952 -0.3066594  1.128649749 0.5387172
P4 = Output$GenoCond[,'p adj']
P4
##     c_b-c_a     c_c-c_a     s_a-c_a     s_b-c_a     s_c-c_a     c_c-c_b 
## 0.154121480 0.776496855 0.066648001 0.031211227 0.551370532 0.006050363 
##     s_a-c_b     s_b-c_b     s_c-c_b     s_a-c_c     s_b-c_c     s_c-c_c 
## 0.999028308 0.934298029 0.959740162 0.001918290 0.001144490 0.047015370 
##     s_b-s_a     s_c-s_a     s_c-s_b 
## 0.990002795 0.828079883 0.538717242
stat.test<- multcompLetters(P4)
stat.test
##  c_b  c_c  s_a  s_b  s_c  c_a 
## "ab"  "c" "ab"  "a" "ab" "bc"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
growth_factors$GenoCond <- factor(growth_factors$GenoCond, levels = c("c_c", "c_b", "c_a", "s_c", "s_b", "s_a"))

MR.delta_p_geno <- ggerrorplot(growth_factors, y = "MR.delta", x = "GenoCond", fill="genotype", 
                               color="genotype", 
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="Genotype", ylab="Growth Rate (cm / day)") 
MR.delta_p_geno <- MR.delta_p_geno + rremove("legend")
MR.delta_p_geno <- MR.delta_p_geno + stat_pvalue_manual(test, label = "Tukey", y.position = 3)
MR.delta_p_geno <- MR.delta_p_geno + ggtitle("Main Root Growth")
MR.delta_p_geno

growth_factors$GenoCond <- factor(growth_factors$GenoCond, levels = c("c_c","s_c", "c_b", "s_b", "c_a", "s_a"))
growth_factors$genotype <- factor(growth_factors$genotype, levels = c("la1511", "m058", "m248"))
growth_factors$genotype <- NA
growth_factors$info <- strsplit(growth_factors$root_name, "_")
growth_factors$info[[1]][2]
## [1] "c"
for(i in 1:nrow(growth_factors)){
  growth_factors$genotype[i] <- growth_factors$info[[i]][3]
}
growth_factors
growth_factors$genotype <- factor(growth_factors$genotype, levels = c("c", "b", "a"))

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", fill="cond", color="cond", 
                        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", 
                                                              label = "p.signif", hide.ns = T)
MR.delta_p_geno

Output <- TukeyHSD(aov(LRno.delta ~ GenoCond, data = growth_factors))
P4 = Output$GenoCond[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"


LRno.delta_p_geno <- ggerrorplot(growth_factors, y = "LRno.delta", x = "GenoCond", fill="genotype", 
                                 color="genotype", 
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="Genotype", ylab="LR increase Rate (# LR / day)") 
LRno.delta_p_geno <- LRno.delta_p_geno + rremove("legend") 
LRno.delta_p_geno <- LRno.delta_p_geno + stat_pvalue_manual(test, label = "Tukey", y.position = 11)
LRno.delta_p_geno <- LRno.delta_p_geno + ggtitle("Lateral Root Number Increase")
LRno.delta_p_geno

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

Output <- TukeyHSD(aov(aLRL.delta ~ GenoCond, data = growth_factors))
P4 = Output$GenoCond[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

aLRL.delta_p_geno <- ggerrorplot(growth_factors, y = "aLRL.delta", x = "GenoCond", fill="genotype",
                                   color="genotype", 
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="Genotype", ylab="growth (cm / day)") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + stat_pvalue_manual(test, label = "Tukey", y.position = 1.3)
aLRL.delta_p_geno <- aLRL.delta_p_geno + ggtitle("average Lateral Root Growth")
aLRL.delta_p_geno

aLRL.delta_p_geno <- ggerrorplot(growth_factors, y = "aLRL.delta", x = "cond", fill="cond", color="cond", 
                        facet.by = c("genotype"), ncol=3,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Average lateral root growth") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "c", 
                                                              label = "p.signif", hide.ns = T)
aLRL.delta_p_geno

Calculating the relative growth rates (Salt Tolerance Indexes):

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

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

library(doBy)
avg_growth <- summaryBy(data = growth_factors, MR.delta + aLRL.delta + LRno.delta ~ genotype + cond)
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$cond == "c")
avg_growth_C <- avg_growth_C[,c(1,3:5)]
colnames(avg_growth_C) <- gsub(".mean", ".avg.Control", colnames(avg_growth_C))
avg_growth_C

now let’s merge it into the growth_factors:

STI_growth_factors <- merge(growth_factors, 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)

and because we are only interested in STI under SALT conditions - we subset this dataset for cond == s

STI <- subset(STI_growth_factors, STI_growth_factors$cond == "s")
STI

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

Output <- TukeyHSD(aov(MR.STI ~ genotype, data = STI))
P4 = Output$genotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
STI$genotype <- factor(STI$genotype, levels = c("c", "b", "a"))

MR.STI_plot <- ggplot(data = STI, mapping = aes(x = genotype, y = MR.STI, colour = genotype)) 
MR.STI_plot <- MR.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
MR.STI_plot <- MR.STI_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
MR.STI_plot <- MR.STI_plot + theme(legend.position = "none")
MR.STI_plot <- MR.STI_plot + xlab("")
MR.STI_plot <- MR.STI_plot + ylab("Fraction of control") + ggtitle("Salt Tolerance Index based on Main Root Growth")
MR.STI_plot <- MR.STI_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 1.2)
MR.STI_plot <- MR.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
MR.STI_plot

Then for aLRL:

Output <- TukeyHSD(aov(aLRL.STI ~ genotype, data = STI))
P4 = Output$genotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

aLRL.STI_plot <- ggplot(data = STI, mapping = aes(x = genotype, y = aLRL.STI, colour = genotype)) 
aLRL.STI_plot <- aLRL.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
aLRL.STI_plot <- aLRL.STI_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
aLRL.STI_plot <- aLRL.STI_plot + theme(legend.position = "none")
aLRL.STI_plot <- aLRL.STI_plot + xlab("")
aLRL.STI_plot <- aLRL.STI_plot + ylab("Salt tolerance index based on average lateral root growth (Fraction of control)")  #+ ggtitle("Salt Tolerance Index based on average Lateral Root Growth")
aLRL.STI_plot <- aLRL.STI_plot + stat_pvalue_manual(test, label = "Tukey", y.position = 2.2) 
aLRL.STI_plot <- aLRL.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
aLRL.STI_plot

Output <- TukeyHSD(aov(LRno.STI ~ genotype, data = STI))
P4 = Output$genotype[,'p adj']
stat.test<- multcompLetters(P4)
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"

LRno.STI_plot <- ggplot(data = STI, mapping = aes(x = genotype, y = LRno.STI, colour = genotype)) 
LRno.STI_plot <- LRno.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
LRno.STI_plot <- LRno.STI_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
LRno.STI_plot <- LRno.STI_plot + theme(legend.position = "none")
LRno.STI_plot <- LRno.STI_plot + xlab("")
LRno.STI_plot <- LRno.STI_plot + ylab("Salt tolerance index based on average increase in lateral root number (Fraction of control)")  #+ ggtitle("Salt Tolerance Index based on average Increase in Lateral Root Number")
LRno.STI_plot <- LRno.STI_plot  + stat_pvalue_manual(test, label = "Tukey", y.position = 2.2)  
LRno.STI_plot <- LRno.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
LRno.STI_plot

####I decided to plot Main root length, Lateral root length, and lateral root number, and main root growth together in one pdf. Let’s see

library(cowplot)

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

lets try to create a single pdf for all 4 graphs here

library(cowplot)
plot_grid(MR_time_graph, LRL_time_graph, LRno_time_graph, MR.delta_p_geno, labels = c("AUTO"), ncol = 2)

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

pdf("20211128_M248M058LA1511_RSA_analysis_STI.pdf", width = 13, height = 10)
plot_grid(MR.STI_plot, LRno.STI_plot, aLRL.STI_plot, labels = c("AUTO"), ncol = 2)
dev.off()
## png 
##   2
plot_grid(MR_time_graph, LRL_time_graph, LRno_time_graph, MR.delta_p_geno,LRno.delta_p_geno,aLRL.delta_p_geno, labels = c("AUTO"), ncol = 3)

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