This is an analysis of the tomato root system architecture dynamics. The experiment was performed between 27th of March 2019 and 4th of April 2019, where eight tomato lines were screened for their dynamic response to salt stress, focusing on the root growth dynamics and changes in morphology.

The material used in this experiment is:

“Regarding the NILs 157-14 and 157-17, we have purified them into 4 NILs because we detected segregation at marker SSR66 in chromosome 2 (a region involved in plant vigor). Their names and genotypes are (L, lycopersicum alelle; C, cheesmaniae allele):

NIL HKT1;1 HKT1;2 SSR66

14-8-1 CC CC LL

14-8-2 CC CC CC

17-2-4 LL LL LL

17-2-2 LL LL CC

"

Pimpinellifolium accessions selected based on the sodium accumulation in the shoot and the salt tolerance index of biomass. Data here >> https://mmjulkowska.shinyapps.io/pimp_app/

The experiment was performed as described in protocol >> “Studying Root System Architecture changes in tomato (S.lycopersicum and S.pimpinellifolium) in response to stress” on protocols.io at Magda’s profile (https://www.protocols.io/researchers/magdalena-julkowska)

The seeds sterilized at 27/03/2019 The seeds in GC for germination at 26/24C 12/12h day/night at 28/03/2019 The seedlings transplanted at 31/03/2019 Plates scanned and analyzed at 31/03/2019, 01/04/2019, 02/04/2019, 03/04/2019

The RSA was analyzed using SmartRoot on scanned plate pictures at 300 dpi resolution.

The data for this project is stored here:

setwd("~/Dropbox/Data and analysis/Tomato/tomato all 201903/")
my_data <- read.csv("tomato_201903_day0_to_day4.csv")
head(my_data)

As you see above, the SmartRoot provides us with the data on every single root - whether it is a main root or lateral root. First of all, we need to convert this into a datafile which contains only the data we want (length) of individual plants as well as the lateral root number, density, length and average length of lateral roots.

# create a list of root ID only for the main roots:
onlyMR <- subset(my_data, my_data$root_order == 0)
# create a data frame containing only LR:
onlyLR <- subset(my_data, my_data$root_order > 0)
# add the other neccessary collumns in the Main Root file:
onlyMR$total.LR.length <- "0"
onlyMR$average.LR.length <- "0"


# then - let's run the loop to collapse the information into the onlyMR file:
for(i in 1:length(unique(onlyLR$parent))){
  parent <- unique(onlyLR$parent)[i] 
  temp2 <- onlyLR[onlyLR$parent %in% parent,]
  totalLRL <- sum(temp2$length)
  
  number_LR <- dim(temp2)[[1]]
  avgLRL <- totalLRL / number_LR
  
  onlyMR[onlyMR$root %in% parent,]$total.LR.length <- totalLRL
  onlyMR[onlyMR$root %in% parent,]$average.LR.length <- avgLRL
}

head(onlyMR)

Then - let’s make sure that columns that are supposed to be numeric are indeed numeric, add other neccessary collumns and remove the not neccessary collumns:

onlyMR$total.LR.length <- as.numeric(onlyMR$total.LR.length)
onlyMR$Total.Root.Length <- onlyMR$length + onlyMR$total.LR.length
colnames(onlyMR)
##  [1] "image"                 "root_name"            
##  [3] "root"                  "length"               
##  [5] "vector_length"         "surface"              
##  [7] "volume"                "direction"            
##  [9] "diameter"              "root_order"           
## [11] "root_ontology"         "parent_name"          
## [13] "parent"                "insertion_position"   
## [15] "insertion_angle"       "n_child"              
## [17] "child_density"         "first_child"          
## [19] "insertion_first_child" "last_child"           
## [21] "insertion_last_child"  "total.LR.length"      
## [23] "average.LR.length"     "Total.Root.Length"
final_data <- onlyMR[,c(1,2,4:9,16:17,19,21:24)]
colnames(final_data)
##  [1] "image"                 "root_name"            
##  [3] "length"                "vector_length"        
##  [5] "surface"               "volume"               
##  [7] "direction"             "diameter"             
##  [9] "n_child"               "child_density"        
## [11] "insertion_first_child" "insertion_last_child" 
## [13] "total.LR.length"       "average.LR.length"    
## [15] "Total.Root.Length"
colnames(final_data)[3] <- "MR.path.length"
colnames(final_data)[4] <- "MR.vector.length"
colnames(final_data)[5] <- "MR.surface"
colnames(final_data)[6] <- "MR.volume"
colnames(final_data)[7] <- "MR.direction"
colnames(final_data)[8] <- "MR.diameter"
colnames(final_data)[9] <- "LR.number"
colnames(final_data)[10] <- "LR.density"
colnames(final_data)[11] <- "first.LR.position"
colnames(final_data)[12] <- "last.LR.position"
colnames(final_data)[13] <- "LRL"
colnames(final_data)[14] <- "aLRL"

colnames(final_data)
##  [1] "image"             "root_name"         "MR.path.length"   
##  [4] "MR.vector.length"  "MR.surface"        "MR.volume"        
##  [7] "MR.direction"      "MR.diameter"       "LR.number"        
## [10] "LR.density"        "first.LR.position" "last.LR.position" 
## [13] "LRL"               "aLRL"              "Total.Root.Length"
head(final_data)

OK - the data looks like I expected it to look, so now I can save it in a csv file:

write.csv(final_data, "Tomato.RSA.experiment1.csv", row.names = F)

Then - I have to brush up the csv file in excell - to add “time” and de-construct my sample coding magic.

Once I am done, I save the file as “Tomato.RSA.brushed.experiment1.csv” and load it up again for the analysis:

new_data <- read.csv("Tomato.RSA.brushed.experiment1.csv")
head(new_data)

As we would like to examine the dynamic change in most of the RSA, we need now to match the data from individual plates.

Let’s have a look at the data and how it looks for individual genotypes first:

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
NIL <- subset(new_data, new_data$Genotype == "NIL14.8.1")

my_graph <- ggplot(data=NIL, aes(x= Day, y=MR.path.length, group = plate.no)) 
my_graph <- my_graph + geom_line() 
my_graph <- my_graph + scale_colour_manual(values = c("grey")) 
my_graph <- my_graph + facet_grid(~ Treatment) 
my_graph <- my_graph + ylab("length") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
ggplotly(my_graph)

However, we can also add some modelled data to this graph:

# Let's add an average Growth Rate added to this graph

library(doBy)
sum_ski <- summaryBy(MR.path.length ~ Treatment + Day, data = NIL)
sum_ski
# let's split the dataset so it will be easier to modell
Control <- subset(sum_ski, sum_ski$Treatment == "1/4MS")
Control <- na.omit(Control)
Control
Salt <- subset(sum_ski, sum_ski$Treatment == "100 mM NaCl")
Salt <- na.omit(Salt)

MRL_c <- Control$MR.path.length.mean
MRL_c
## [1] 4.506864 6.870674 8.775666 9.455610
time_d <- Control$Day
model_C <- lm(MRL_c~ time_d)
model_C
## 
## Call:
## lm(formula = MRL_c ~ time_d)
## 
## Coefficients:
## (Intercept)       time_d  
##       4.890        1.675
MRL_s <- Salt$MR.path.length.mean
time_d <- Salt$Day
model_S <- lm(MRL_s~ time_d)
model_S
## 
## Call:
## lm(formula = MRL_s ~ time_d)
## 
## Coefficients:
## (Intercept)       time_d  
##       4.611        1.749
timevalues <- unique(time_d)

MRL.C <- predict(model_C,list(Time=timevalues))
MRL.S <- predict(model_S,list(Time=timevalues))

# create a dataframe with 4 rows - because we have 4 unique time points
Open_for_model <-data.frame(Day=numeric(4),plate.no=character(4), Genotype=character(4), Treatment=character(4), MR.path.length=numeric(4))

# add the model values into the data frame:
Control <- Open_for_model
Salt <- Open_for_model

unique(time_d)
## [1] 0 1 2 3
Control$Day <- unique(time_d)
Salt$Day <- unique(time_d)

Control$plate.no <- "model"
Salt$plate.no <- "model"

Control$Genotype <- "modelled"
Salt$Genotype <- "modelled"

Control$Treatment <- "1/4MS"
Salt$Treatment <- "100 mM NaCl"

Control$MR.path.length <- MRL.C
Salt$MR.path.length <- MRL.S

model_total <- rbind(Control, Salt)
head(model_total)
head(NIL)
my_data_relevant <- subset(NIL, select = c("Day", "plate.no", "Genotype", "Treatment", "MR.path.length"))
for_graph <- rbind(my_data_relevant, model_total)

ggplot(data=for_graph, aes(x= Day, y=MR.path.length, group = plate.no, color = Genotype)) + 
  geom_line() + 
  scale_colour_manual(values = c("grey", "blue")) +
  facet_grid(~ Treatment) + 
  ylab("length (cm)") + xlab("time (days)") + theme(legend.position='none')

quartz.save("NIL14.8.1_MR_growth_fitted.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2

Let’s loop it for all the other genotypes too:

unique(new_data$Genotype)[1]
## [1] NIL14.8.1
## 9 Levels: M042.4 M050.11 M218.3 M234.1 NIL14.8.1 NIL14.8.2 ... S.lyco
for(i in 1:9){  

  geno <- unique(new_data$Genotype)[i]
  NIL <- subset(new_data, new_data$Genotype == geno)
  sum_ski <- summaryBy(MR.path.length ~ Treatment + Day, data = NIL)

# let's split the dataset so it will be easier to modell
  Control <- subset(sum_ski, sum_ski$Treatment == "1/4MS")
  Control <- na.omit(Control)
  Control

  Salt <- subset(sum_ski, sum_ski$Treatment == "100 mM NaCl")
  Salt <- na.omit(Salt)

  MRL_c <- Control$MR.path.length.mean
  time_d <- Control$Day
  model_C <- lm(MRL_c~ time_d)

  MRL_s <- Salt$MR.path.length.mean
  time_d <- Salt$Day
  model_S <- lm(MRL_s~ time_d)


  timevalues <- unique(time_d)
  MRL.C <- predict(model_C,list(Time=timevalues))
  MRL.S <- predict(model_S,list(Time=timevalues))

# create a dataframe with 4 rows - because we have 4 unique time points
  Open_for_model <-data.frame(Day=numeric(4),plate.no=character(4), Genotype=character(4), Treatment=character(4), MR.path.length=numeric(4))

# add the model values into the data frame:
  Control <- Open_for_model
  Salt <- Open_for_model

  Control$Day <- unique(time_d)
  Salt$Day <- unique(time_d)

  Control$plate.no <- "model"
  Salt$plate.no <- "model"

  Control$Genotype <- "modelled"
  Salt$Genotype <- "modelled"

  Control$Treatment <- "1/4MS"
  Salt$Treatment <- "100 mM NaCl"

  Control$MR.path.length <- MRL.C
  Salt$MR.path.length <- MRL.S

  model_total <- rbind(Control, Salt)
  my_data_relevant <- subset(NIL, select = c("Day", "plate.no", "Genotype", "Treatment", "MR.path.length"))
  for_graph <- rbind(my_data_relevant, model_total)

  pdf(paste(geno,"_MR_growth_fitted.pdf", sep=""))

  my_plot <- ggplot(data=for_graph, aes(x= Day, y=MR.path.length, group = plate.no, color = Genotype)) 
  my_plot <- my_plot + geom_line() + scale_colour_manual(values = c("grey", "blue")) 
  my_plot <- my_plot + ylim(0,12.5)
  my_plot <- my_plot + facet_grid(~ Treatment) + ggtitle ("Main Root Length") + ylab("length (cm)") + xlab("time (days)")
  my_plot <- my_plot + theme(legend.position='none')
  print(my_plot)
  dev.off()

}
## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

OK - now let’s do the same for the Lateral Root Number:

for(i in 1:9){  

  geno <- unique(new_data$Genotype)[i]
  NIL <- subset(new_data, new_data$Genotype == geno)
  sum_ski <- summaryBy(LR.number ~ Treatment + Day, data = NIL)

# let's split the dataset so it will be easier to modell
  Control <- subset(sum_ski, sum_ski$Treatment == "1/4MS")
  Control <- na.omit(Control)
  Control

  Salt <- subset(sum_ski, sum_ski$Treatment == "100 mM NaCl")
  Salt <- na.omit(Salt)

  MRL_c <- Control$LR.number.mean
  time_d <- Control$Day
  model_C <- lm(MRL_c~ time_d)

  MRL_s <- Salt$LR.number.mean
  time_d <- Salt$Day
  model_S <- lm(MRL_s~ time_d)


  timevalues <- unique(time_d)
  MRL.C <- predict(model_C,list(Time=timevalues))
  MRL.S <- predict(model_S,list(Time=timevalues))

# create a dataframe with 4 rows - because we have 4 unique time points
  Open_for_model <-data.frame(Day=numeric(4),plate.no=character(4), Genotype=character(4), Treatment=character(4), LR.number=numeric(4))

# add the model values into the data frame:
  Control <- Open_for_model
  Salt <- Open_for_model

  Control$Day <- unique(time_d)
  Salt$Day <- unique(time_d)

  Control$plate.no <- "model"
  Salt$plate.no <- "model"

  Control$Genotype <- "modelled"
  Salt$Genotype <- "modelled"

  Control$Treatment <- "1/4MS"
  Salt$Treatment <- "100 mM NaCl"

  Control$LR.number <- MRL.C
  Salt$LR.number <- MRL.S

  model_total <- rbind(Control, Salt)
  my_data_relevant <- subset(NIL, select = c("Day", "plate.no", "Genotype", "Treatment", "LR.number"))
  for_graph <- rbind(my_data_relevant, model_total)

  pdf(paste(geno,"_LR.number_growth_fitted.pdf", sep=""))

  my_plot <- ggplot(data=for_graph, aes(x= Day, y=LR.number, group = plate.no, color = Genotype)) 
  my_plot <- my_plot + geom_line() + scale_colour_manual(values = c("grey", "blue")) 
  my_plot <- my_plot + ylim(0,23)
  my_plot <- my_plot + facet_grid(~ Treatment) + ggtitle ("Lateral Root Number") + ylab("length (cm)") + xlab("time (days)")
  my_plot <- my_plot + theme(legend.position='none')
  print(my_plot)
  dev.off()

}
## Warning: Removed 3 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 2 rows containing missing values (geom_path).

## Warning: Removed 2 rows containing missing values (geom_path).
## Warning: Removed 4 rows containing missing values (geom_path).
## Warning: Removed 1 rows containing missing values (geom_path).

And for average Lateral Root Length:

for(i in 1:9){  

  geno <- unique(new_data$Genotype)[i]
  NIL <- subset(new_data, new_data$Genotype == geno)
  sum_ski <- summaryBy(aLRL ~ Treatment + Day, data = NIL)

# let's split the dataset so it will be easier to modell
  Control <- subset(sum_ski, sum_ski$Treatment == "1/4MS")
  Control <- na.omit(Control)
  Control

  Salt <- subset(sum_ski, sum_ski$Treatment == "100 mM NaCl")
  Salt <- na.omit(Salt)

  MRL_c <- Control$aLRL.mean
  time_d <- Control$Day
  model_C <- lm(MRL_c~ time_d)

  MRL_s <- Salt$aLRL.mean
  time_d <- Salt$Day
  model_S <- lm(MRL_s~ time_d)


  timevalues <- unique(time_d)
  MRL.C <- predict(model_C,list(Time=timevalues))
  MRL.S <- predict(model_S,list(Time=timevalues))

# create a dataframe with 4 rows - because we have 4 unique time points
  Open_for_model <-data.frame(Day=numeric(4),plate.no=character(4), Genotype=character(4), Treatment=character(4), aLRL=numeric(4))

# add the model values into the data frame:
  Control <- Open_for_model
  Salt <- Open_for_model

  Control$Day <- unique(time_d)
  Salt$Day <- unique(time_d)

  Control$plate.no <- "model"
  Salt$plate.no <- "model"

  Control$Genotype <- "modelled"
  Salt$Genotype <- "modelled"

  Control$Treatment <- "1/4MS"
  Salt$Treatment <- "100 mM NaCl"

  Control$aLRL <- MRL.C
  Salt$aLRL <- MRL.S

  model_total <- rbind(Control, Salt)
  my_data_relevant <- subset(NIL, select = c("Day", "plate.no", "Genotype", "Treatment", "aLRL"))
  for_graph <- rbind(my_data_relevant, model_total)

  pdf(paste(geno,"_aLRL_growth_fitted.pdf", sep=""))

  my_plot <- ggplot(data=for_graph, aes(x= Day, y=aLRL, group = plate.no, color = Genotype)) 
  my_plot <- my_plot + geom_line() + scale_colour_manual(values = c("grey", "blue")) 
  my_plot <- my_plot + ylim(0,2)
  my_plot <- my_plot + facet_grid(~ Treatment) + ggtitle ("average Lateral Root Length") + ylab("length (cm)") + xlab("time (days)")
  my_plot <- my_plot + theme(legend.position='none')
  print(my_plot)
  dev.off()

}
## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 1 rows containing missing values (geom_path).

OK - this is all nice representation, but we dont have the Growth Rate numbers for all of these plants. In order to calculate the growth rate per plant, we will run the following loop:

head(new_data)
unique(new_data$plate.no)
##   [1] 144 71  42  13  70  100 141 113 128 55  14  99  86  143 85  12  69 
##  [18] 129 93  27  39  15  59  53  140 87  115 84  54  82  112 127 98  23 
##  [35] 38  16  116 101 41  142 114 52  60  111 139 81  126 24  97  37  51 
##  [52] 61  110 138 7   96  35  19  76  65  28  80  125 25  36  8   124 9  
##  [69] 17  47  121 134 102 26  50  109 137 62  79  78  63  106 92  32  130
##  [86] 95  49  108 135 123 94  34  10  2   20  66  89  44  68  118 103 29 
## [103] 131 5   72  43  56  105 133 91  3   74  119 90  104 4   22  11  46 
## [120] 75  120 31  21  67  45  132 30  73  117 136 18  48  77  64  122 107
## [137] 33  1   88  40  58  57  6   148 215 DN 
## 146 Levels: 1 10 100 101 102 103 104 105 106 107 108 109 11 110 111 ... DN
new_data$PlantID <- paste(new_data$plate.no, new_data$Treatment, new_data$Genotype, sep="_")
length(unique(new_data$PlantID))
## [1] 150
head(new_data)
i=1

the_one <- unique(new_data$PlantID)[i]
uni <- subset(new_data, new_data$PlantID == the_one)

names <- c(text="Plant.ID", "Genotype", "Treatment", "plate.no", "intercept", "delta")
growth_factors <- data.frame()
growth_factors
for (k in names) growth_factors[[k]] <- as.character()

# let's get also all the individual identifiers in the table:
growth_factors[i,1] <- as.character(unique(uni$PlantID))
growth_factors[i,2] <- as.character(unique(uni$Genotype))
growth_factors[i,3] <- as.character(unique(uni$Treatment))
growth_factors[i,4] <- as.character(unique(uni$plate.no))

# let's calculate the model:  
Area_mm2 <- uni$MR.path.length
time_d <- uni$Day
model_C <- lm(Area_mm2~ time_d)
# add model parts into the main table with growth factors
growth_factors[i,5] <- as.numeric(model_C$coefficients[[1]])
growth_factors[i,6] <- as.numeric(model_C$coefficients[[2]])

# calculate predicted Area values for this specific sample:
timevalues <- unique(time_d)
timevalues
## [1] 0 1 2 3
Area.pred <- predict(model_C,list(Time=timevalues))
Area.pred
##         1         2         3         4 
##  5.283955  6.873323  8.462690 10.052058
uni$MRL_pred <- Area.pred
done <- uni

growth_factors2 <- growth_factors
growth_factors2
for(i in 2:length(unique(new_data$PlantID))){
  the_one <- unique(new_data$PlantID)[i]
  uni <- subset(new_data, new_data$PlantID == the_one)

# let's get also all the individual identifiers in the table:
  growth_factors[i,1] <- as.character(unique(uni$PlantID))
  growth_factors[i,2] <- as.character(unique(uni$Genotype))
  growth_factors[i,3] <- as.character(unique(uni$Treatment))
  growth_factors[i,4] <- as.character(unique(uni$plate.no))

# let's calculate the model:  
  Area_mm2 <- uni$MR.path.length
  time_d <- uni$Day
  model_C <- lm(Area_mm2~ time_d)
# add model parts into the main table with growth factors
  growth_factors[i,5] <- as.numeric(model_C$coefficients[[1]])
  growth_factors[i,6] <- as.numeric(model_C$coefficients[[2]])

# calculate predicted Area values for this specific sample:
  timevalues <- unique(time_d)
  timevalues

  Area.pred <- predict(model_C,list(Time=timevalues))
  Area.pred

  uni$MRL_pred <- Area.pred
  done <- rbind(done,uni)
  growth_factors

}
## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading
dim(growth_factors2)
## [1] 1 6
head(growth_factors)
dim(done)
## [1] 567  21
dim(growth_factors)
## [1] 150   6
growth_factors
growth_factors <- na.omit(growth_factors)
write.csv(growth_factors, "Main_Root_GR.csv", row.names = FALSE)

Now we can visualize the differences between the individual lines in a graph:

#Let's have an order in which we would like to plot things first:
unique(growth_factors$Genotype)
## [1] "NIL14.8.1" "M042.4"    "NIL14.8.2" "NIL17.2.4" "M234.1"    "M218.3"   
## [7] "S.lyco"    "NIL17.2.2" "M050.11"
growth_factors$Genotype <- factor(growth_factors$Genotype,levels = c("S.lyco", "NIL14.8.1", "NIL14.8.2", "NIL17.2.2", "NIL17.2.4", "M042.4", "M050.11", "M218.3", "M234.1"))
growth_factors$intercept <- as.numeric(growth_factors$intercept)
growth_factors$delta <- as.numeric(growth_factors$delta)

library(ggpubr)
## Loading required package: magrittr
library(ggbeeswarm)
library(gapminder)
library(RColorBrewer)
library(ggridges)
## 
## Attaching package: 'ggridges'
## The following object is masked from 'package:ggplot2':
## 
##     scale_discrete_manual
my_box_plot <- ggplot(data = growth_factors, mapping = aes(x = Genotype, y = intercept, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + facet_wrap(~ Treatment)
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("GR intercept")

my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("MR_GR_all_Intercept.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2
my_box_plot <- ggplot(data = growth_factors, mapping = aes(x = Genotype, y = delta, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + facet_wrap(~ Treatment)
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("GR delta (cm / day)")
my_box_plot <- my_box_plot + ggtitle("Main Root GR")
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("MR_GR_all_Delta.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2

OK - let’s do exactly the same for the LR number:

i=1

the_one <- unique(new_data$PlantID)[i]
uni <- subset(new_data, new_data$PlantID == the_one)

names <- c(text="Plant.ID", "Genotype", "Treatment", "plate.no", "intercept", "delta")
growth_factors <- data.frame()
growth_factors
for (k in names) growth_factors[[k]] <- as.character()

# let's get also all the individual identifiers in the table:
growth_factors[i,1] <- as.character(unique(uni$PlantID))
growth_factors[i,2] <- as.character(unique(uni$Genotype))
growth_factors[i,3] <- as.character(unique(uni$Treatment))
growth_factors[i,4] <- as.character(unique(uni$plate.no))

# let's calculate the model:  
Area_mm2 <- uni$LR.number
time_d <- uni$Day
model_C <- lm(Area_mm2~ time_d)
# add model parts into the main table with growth factors
growth_factors[i,5] <- as.numeric(model_C$coefficients[[1]])
growth_factors[i,6] <- as.numeric(model_C$coefficients[[2]])

# calculate predicted Area values for this specific sample:
timevalues <- unique(time_d)
timevalues
## [1] 0 1 2 3
Area.pred <- predict(model_C,list(Time=timevalues))
Area.pred
##    1    2    3    4 
## -0.5  1.0  2.5  4.0
uni$MRL_pred <- Area.pred
done <- uni

growth_factors2 <- growth_factors
growth_factors2
for(i in 2:length(unique(new_data$PlantID))){
  the_one <- unique(new_data$PlantID)[i]
  uni <- subset(new_data, new_data$PlantID == the_one)

# let's get also all the individual identifiers in the table:
  growth_factors[i,1] <- as.character(unique(uni$PlantID))
  growth_factors[i,2] <- as.character(unique(uni$Genotype))
  growth_factors[i,3] <- as.character(unique(uni$Treatment))
  growth_factors[i,4] <- as.character(unique(uni$plate.no))

# let's calculate the model:  
  Area_mm2 <- uni$LR.number
  time_d <- uni$Day
  model_C <- lm(Area_mm2~ time_d)
# add model parts into the main table with growth factors
  growth_factors[i,5] <- as.numeric(model_C$coefficients[[1]])
  growth_factors[i,6] <- as.numeric(model_C$coefficients[[2]])

# calculate predicted Area values for this specific sample:
  timevalues <- unique(time_d)
  timevalues

  Area.pred <- predict(model_C,list(Time=timevalues))
  Area.pred

  uni$MRL_pred <- Area.pred
  done <- rbind(done,uni)
  growth_factors

}
## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading
dim(growth_factors2)
## [1] 1 6
head(growth_factors)
dim(done)
## [1] 567  21
dim(growth_factors)
## [1] 150   6
growth_factors
growth_factors <- na.omit(growth_factors)
write.csv(growth_factors, "noLR_GR.csv", row.names = FALSE)

# plotting:

#Let's have an order in which we would like to plot things first:
unique(growth_factors$Genotype)
## [1] "NIL14.8.1" "M042.4"    "NIL14.8.2" "NIL17.2.4" "M234.1"    "M218.3"   
## [7] "S.lyco"    "NIL17.2.2" "M050.11"
growth_factors$Genotype <- factor(growth_factors$Genotype,levels = c("S.lyco", "NIL14.8.1", "NIL14.8.2", "NIL17.2.2", "NIL17.2.4", "M042.4", "M050.11", "M218.3", "M234.1"))
growth_factors$intercept <- as.numeric(growth_factors$intercept)
growth_factors$delta <- as.numeric(growth_factors$delta)

my_box_plot <- ggplot(data = growth_factors, mapping = aes(x = Genotype, y = intercept, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + facet_wrap(~ Treatment)
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("GR intercept")

my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("LRno_GR_all_Intercept.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2
my_box_plot <- ggplot(data = growth_factors, mapping = aes(x = Genotype, y = delta, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + facet_wrap(~ Treatment)
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("GR delta (# LR / MR / day)")
my_box_plot <- my_box_plot + ggtitle("Increase in Lateral Root Number")
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("LRno_GR_all_Delta.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2

And we can do the same for average Lateral Root Length:

i=1

the_one <- unique(new_data$PlantID)[i]
uni <- subset(new_data, new_data$PlantID == the_one)

names <- c(text="Plant.ID", "Genotype", "Treatment", "plate.no", "intercept", "delta")
growth_factors <- data.frame()
growth_factors
for (k in names) growth_factors[[k]] <- as.character()

# let's get also all the individual identifiers in the table:
growth_factors[i,1] <- as.character(unique(uni$PlantID))
growth_factors[i,2] <- as.character(unique(uni$Genotype))
growth_factors[i,3] <- as.character(unique(uni$Treatment))
growth_factors[i,4] <- as.character(unique(uni$plate.no))

# let's calculate the model:  
Area_mm2 <- uni$aLRL
time_d <- uni$Day
model_C <- lm(Area_mm2~ time_d)
# add model parts into the main table with growth factors
growth_factors[i,5] <- as.numeric(model_C$coefficients[[1]])
growth_factors[i,6] <- as.numeric(model_C$coefficients[[2]])

# calculate predicted Area values for this specific sample:
timevalues <- unique(time_d)
timevalues
## [1] 0 1 2 3
Area.pred <- predict(model_C,list(Time=timevalues))
Area.pred
##          1          2          3          4 
## -0.0343125  0.1459255  0.3261635  0.5064015
uni$MRL_pred <- Area.pred
done <- uni

growth_factors2 <- growth_factors
growth_factors2
for(i in 2:length(unique(new_data$PlantID))){
  the_one <- unique(new_data$PlantID)[i]
  uni <- subset(new_data, new_data$PlantID == the_one)

# let's get also all the individual identifiers in the table:
  growth_factors[i,1] <- as.character(unique(uni$PlantID))
  growth_factors[i,2] <- as.character(unique(uni$Genotype))
  growth_factors[i,3] <- as.character(unique(uni$Treatment))
  growth_factors[i,4] <- as.character(unique(uni$plate.no))

# let's calculate the model:  
  Area_mm2 <- uni$aLRL
  time_d <- uni$Day
  model_C <- lm(Area_mm2~ time_d)
# add model parts into the main table with growth factors
  growth_factors[i,5] <- as.numeric(model_C$coefficients[[1]])
  growth_factors[i,6] <- as.numeric(model_C$coefficients[[2]])

# calculate predicted Area values for this specific sample:
  timevalues <- unique(time_d)
  timevalues

  Area.pred <- predict(model_C,list(Time=timevalues))
  Area.pred

  uni$MRL_pred <- Area.pred
  done <- rbind(done,uni)
  growth_factors

}
## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading

## Warning in predict.lm(model_C, list(Time = timevalues)): prediction from a
## rank-deficient fit may be misleading
dim(growth_factors2)
## [1] 1 6
head(growth_factors)
dim(done)
## [1] 567  21
dim(growth_factors)
## [1] 150   6
growth_factors
growth_factors <- na.omit(growth_factors)
write.csv(growth_factors, "aLRL_GR.csv", row.names = FALSE)

# plotting:

#Let's have an order in which we would like to plot things first:
unique(growth_factors$Genotype)
## [1] "NIL14.8.1" "M042.4"    "NIL14.8.2" "NIL17.2.4" "M234.1"    "M218.3"   
## [7] "S.lyco"    "NIL17.2.2" "M050.11"
growth_factors$Genotype <- factor(growth_factors$Genotype,levels = c("S.lyco", "NIL14.8.1", "NIL14.8.2", "NIL17.2.2", "NIL17.2.4", "M042.4", "M050.11", "M218.3", "M234.1"))
growth_factors$intercept <- as.numeric(growth_factors$intercept)
growth_factors$delta <- as.numeric(growth_factors$delta)

my_box_plot <- ggplot(data = growth_factors, mapping = aes(x = Genotype, y = intercept, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + facet_wrap(~ Treatment)
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("GR intercept")

my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("aLRL_GR_all_Intercept.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2
my_box_plot <- ggplot(data = growth_factors, mapping = aes(x = Genotype, y = delta, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + facet_wrap(~ Treatment)
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("GR delta (cm LR / day)")
my_box_plot <- my_box_plot + ggtitle("average Lateral Root Growth")
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("aLRL_GR_all_Delta.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2

And then we would also like to have a look at the relative effect of salt for each of these lines. Thus… we need to calculate the average for control conditions of each MR, LRno and aLRL GR and preferably put them all into one dataframe:

aLRL <- read.csv("aLRL_GR.csv")
MR <- read.csv("Main_Root_GR.csv")
noLR <- read.csv("noLR_GR.csv")

colnames(MR)[6] <- "GR_MR"
MR <- MR[,c(1:4,6)]
colnames(aLRL)[6] <- "GR_aLRL"
aLRL <- aLRL[,c(1:4,6)]
colnames(noLR)[6] <- "GR_noLR"
noLR <- noLR[,c(1:4,6)]

library(reshape)
## 
## Attaching package: 'reshape'
## The following object is masked from 'package:plotly':
## 
##     rename
all_GR <- merge(MR, aLRL, by=c("Plant.ID", "Genotype", "Treatment", "plate.no"))
all_GR <- merge(all_GR, noLR, by=c("Plant.ID", "Genotype", "Treatment", "plate.no"))
head(all_GR)
sum_GR <- summaryBy(GR_MR + GR_aLRL + GR_noLR ~ Treatment + Genotype, data = all_GR)
head(sum_GR)
c_GR <- subset(sum_GR, sum_GR$Treatment == "1/4MS")
c_GR <- c_GR[,2:5]

rel_GR<- merge(all_GR, c_GR, by=c("Genotype"))
rel_GR
rel_GR$Rel_MR <- rel_GR$GR_MR/rel_GR$GR_MR.mean
rel_GR$Rel_aLRL <- rel_GR$GR_aLRL/rel_GR$GR_aLRL.mean
rel_GR$Rel_noLR <- rel_GR$GR_noLR/rel_GR$GR_noLR.mean
rel_GR
rel_GR_salt <- subset(rel_GR, rel_GR$Treatment == "100 mM NaCl")
colnames(rel_GR_salt)
##  [1] "Genotype"     "Plant.ID"     "Treatment"    "plate.no"    
##  [5] "GR_MR"        "GR_aLRL"      "GR_noLR"      "GR_MR.mean"  
##  [9] "GR_aLRL.mean" "GR_noLR.mean" "Rel_MR"       "Rel_aLRL"    
## [13] "Rel_noLR"
rel_GR_salt <- rel_GR_salt[,c(1,11:13)]

head(rel_GR_salt)
# ok so from here we can already do some visualizations:

rel_GR_salt$Genotype <- factor(rel_GR_salt$Genotype,levels = c("S.lyco", "NIL14.8.1", "NIL14.8.2", "NIL17.2.2", "NIL17.2.4", "M042.4", "M050.11", "M218.3", "M234.1"))

my_box_plot <- ggplot(data = rel_GR_salt, mapping = aes(x = Genotype, y = Rel_MR, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("fraction of control")
my_box_plot <- my_box_plot + ggtitle("Relative Main Root Growth")
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("Relative_MR_GR_all_Delta.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2
my_box_plot <- ggplot(data = rel_GR_salt, mapping = aes(x = Genotype, y = Rel_aLRL, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("fraction of control")
my_box_plot <- my_box_plot + ggtitle("Relative average Lateral Root Growth")
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("Relative_aLRL_GR_all_Delta.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2
my_box_plot <- ggplot(data = rel_GR_salt, mapping = aes(x = Genotype, y = Rel_noLR, colour = Genotype)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + theme(legend.position = "none")
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("fraction of control")
my_box_plot <- my_box_plot + ggtitle("Relative Lateral Root Number")
my_box_plot <- my_box_plot + theme(axis.text.x = element_text(angle=70, vjust=0.5))
my_box_plot <- my_box_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "S.lyco") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("royalblue", "coral3", "chocolate", "darkseagreen", "darkseagreen3", "tomato4", "tomato3", "tomato2", "tomato1"))
my_box_plot

quartz.save("Relative_noLR_GR_all_Delta.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2
# and now let's go for the finale - and compare the relative effect on RSA:

colnames(rel_GR_salt)[2] <- "Rel. MR GR"
colnames(rel_GR_salt)[3] <- "Rel. aLRL GR"
colnames(rel_GR_salt)[4] <- "Rel. LR# GR"
melt_salt <- melt(rel_GR_salt, id="Genotype")

my_comparisons <- list(c("Rel. MR GR","Rel. aLRL GR"), c("Rel. aLRL GR", "Rel. LR# GR"), c("Rel. MR GR","Rel. LR# GR"))

my_box_plot <- ggplot(data = melt_salt, mapping = aes(x = variable, y = value, colour = variable)) 
my_box_plot <- my_box_plot + geom_beeswarm(alpha=0.6, priority = "density")
my_box_plot <- my_box_plot + stat_summary(fun.y=mean, geom="point", shape=95, size=10, color="black", fill="black")
my_box_plot <- my_box_plot + facet_wrap(~ Genotype, ncol = 9)
my_box_plot <- my_box_plot + theme(axis.text.x = element_blank(),
                                   axis.ticks.x = element_blank(),
                                   axis.title.x = element_blank(),
                                   legend.position = "bottom",
                                   legend.title = element_blank(),
                                   strip.text.x = element_text(size = 10, face = "bold"))
my_box_plot <- my_box_plot + xlab("")
my_box_plot <- my_box_plot + ylab("fraction of control")
my_box_plot <- my_box_plot + ggtitle("Relative Changes in Root System Architecture")
my_box_plot <- my_box_plot + stat_compare_means(comparisons = my_comparisons, label = "p.signif") 
my_box_plot <- my_box_plot + scale_color_manual(values=c("deepskyblue3", "gold", "orchid3"))
my_box_plot

quartz.save("Relative_ALL_GR_all.pdf", type= "pdf", dpi=200)
## quartz_off_screen 
##                 2