this is the RSA experiment for the 5 Arabidopsis genotypes.4 d germination in 1/2 ms, at d5, transferred to the +/- 0, 75, 125 salt, and scanned 5 times, for everyother day. root and shoots harvested after two weeks of stress. a.col, b.duf,c.wrky,d/e.2xko.

getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20231120_RSA_duf_wrky_2xko"
list.files(pattern = ".csv")
## [1] "20231120_duf_wrky_2xko_growth_factors.csv"
## [2] "d11.csv"                                  
## [3] "d13.csv"                                  
## [4] "d5.csv"                                   
## [5] "d7.csv"                                   
## [6] "d9.csv"
d5<-read.csv("d5.csv")
d7<-read.csv("d7.csv")
d9<-read.csv("d9.csv")
d11<-read.csv("d11.csv")
d13<-read.csv("d13.csv")
head(d5)
d5$day <- 5
d7$day <- 7
d9$day <- 9
d11$day <- 11
d13$day <- 13

head(d13)
all_data <- rbind(d5, d7)
all_data <- rbind(all_data, d9)
all_data <- rbind(all_data,d11)
all_data <- rbind(all_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] " c33ea405-4f9d-480f-84f0-cd20721cfe26"
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.2559965
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
dim(MR_data_3)
## [1] 584  12
for(i in 1:584){
  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
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]  5  7  9 11 13
?strsplit()
## starting httpd help server ... done
MR_all$root_name[12]
## [1] " pl1-125-b"
text <- strsplit(x = MR_all$root_name[12], split = "-")
text
## [[1]]
## [1] " pl1" "125"  "b"
plate <- text[[1]][1]
plate
## [1] " pl1"
cond <- text[[1]][2]
cond
## [1] "125"
genotype <- text[[1]][3]
genotype
## [1] "b"

OK - now let’s loop it:

dim(MR_all)
## [1] 916  12
for(i in 1:916){
  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] "a" "c" "e" "d" "b" NA
MR_all.nona <- na.omit(MR_all)
unique(MR_all.nona$genotype)
## [1] "a" "c" "e" "d" "b"
#get genotype information as it should be:
#MR_all.nona$genotype <- gsub("c1", "Col-0", MR_all.nona$genotype)
#MR_all.nona$genotype <- gsub("c2", "Col-0", MR_all.nona$genotype)
#MR_all.nona$genotype <- gsub("duf1", "duf-5", MR_all.nona$genotype)
#MR_all.nona$genotype <- gsub("duf2", "duf-5", MR_all.nona$genotype)
#unique(MR_all.nona$genotype)

unique(MR_all.nona$condition)
## [1] "125" "75"  "0"
#MR_all.nona$condition<- gsub("c","0", MR_all.nona$condition)
dim(MR_all.nona)
## [1] 858  15
unique(MR_all$day)
## [1]  5  7  9 11 13
MR_all$condition <- as.numeric(MR_all$condition)
MR_all
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
MR_all.nona$condition <- as.numeric(MR_all.nona$condition)
MR_all.nona
library(ggplot2)
library(ggpubr)

MR_all <- MR_all.nona
unique(MR_all.nona$condition)
## [1] 125  75   0
MR_all.nona$condition <- as.factor(MR_all.nona$condition)


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

pdf("histogram.TRS.pdf")
plot(histogram_TRS)
# 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_LRL, ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## png 
##   2
unique(MR_all$genotype)
## [1] "a" "c" "e" "d" "b"
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 = "a") 
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: 14, 5, 2, 8, 11
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 3, 9, 15, 12, 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: 11, 8, 5, 14, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 3, 9, 15, 6, 12
## 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: 11, 8, 5, 14, 2
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 3, 15, 9, 12, 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: 2, 14, 11, 8, 5
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 3, 15, 12, 9, 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: 11, 8, 5, 2, 14
## Warning: Computation failed in `stat_compare_means()`
## Caused by error:
## ! Can't find specified reference group: 1. Allowed values include one of: 15, 3, 12, 9, 6

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_a <- subset(MR_all.nona, MR_all.nona$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)
geno_b <- subset(MR_all.nona, MR_all.nona$genotype == "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)
geno_c <- subset(MR_all.nona, MR_all.nona$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)
geno_d <- subset(MR_all.nona, MR_all.nona$genotype == "d")

my_graph <- ggplot(data=geno_d, 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)
geno_e <- subset(MR_all.nona, MR_all.nona$genotype == "e")

my_graph <- ggplot(data=geno_e, 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)
MR_all_geno <- rbind(geno_a, geno_b,geno_c, geno_d,geno_e)

MR_time_graph <- ggplot(data=MR_all_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", "plum", "rosybrown1", "hotpink","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_all_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 + scale_color_manual(values= c("blue", "plum", "rosybrown1", "hotpink","red"))

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_all_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 + scale_color_manual(values= c("blue", "plum", "rosybrown1", "hotpink","red"))

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_all_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 + scale_color_manual(values= c("blue", "plum", "rosybrown1", "hotpink","red"))
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_all_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 + scale_color_manual(values= c("blue", "plum", "rosybrown1", "hotpink","red"))

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 312 rows containing non-finite values (`stat_summary()`).
## Warning: Removed 300 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_all_geno, MR_all_geno$root_name == unique(MR_all_geno$root_name)[1])
temp1
temp2 <- temp1[order(temp1$day),]
temp2
temp2$condition[1]
## [1] 125
## 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)

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  
##     -0.7286       0.3021
summary(MR_model)
## 
## Call:
## lm(formula = temp_MR$length ~ temp_MR$day)
## 
## Residuals:
##         1         2         3         4         5 
##  0.097690 -0.055932 -0.110891 -0.001185  0.070317 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.72859    0.14895  -4.891 0.016345 *  
## temp_MR$day  0.30206    0.01579  19.131 0.000312 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.09986 on 3 degrees of freedom
## Multiple R-squared:  0.9919, Adjusted R-squared:  0.9892 
## F-statistic:   366 on 1 and 3 DF,  p-value: 0.0003119
MR_model$coefficients[[2]]
## [1] 0.302063
MR_growth_rate <- MR_model$coefficients[[2]]
# 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  
##        -4.25          0.85
summary(LRno_model)
## 
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
## 
## Residuals:
##    1    2    3    4    5 
##  1.0 -0.7 -1.4  0.9  0.2 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)  
## (Intercept)   -4.2500     1.7858   -2.38   0.0976 .
## LR_temp2$day   0.8500     0.1893    4.49   0.0206 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.197 on 3 degrees of freedom
## Multiple R-squared:  0.8705, Adjusted R-squared:  0.8273 
## F-statistic: 20.16 on 1 and 3 DF,  p-value: 0.02061
LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase
## [1] 0.85

###why my average lateral root length is negative!!!!Well, because they are not growing anymore.

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.28170      -0.01171
summary(aLRL_model)
## 
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
## 
## Residuals:
##        1        2        3        4        5 
##  0.03287  0.07014 -0.13474 -0.07241  0.10414 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)
## (Intercept)   0.28170    0.17276   1.631    0.201
## LR_temp2$day -0.01171    0.01831  -0.640    0.568
## 
## Residual standard error: 0.1158 on 3 degrees of freedom
## Multiple R-squared:   0.12,  Adjusted R-squared:  -0.1733 
## F-statistic: 0.4092 on 1 and 3 DF,  p-value: 0.5679
aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth
## [1] -0.01171444
temp2$condition[1]
## [1] 125
## Levels: 0 75 125
temp2$condition[1][1]
## [1] 125
## 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
length(unique(MR_all_geno$root_name))
## [1] 179
for(e in 1:179){
  temp1 <- subset(MR_all_geno, MR_all_geno$root_name == unique(MR_all_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
growth_factors <- subset(growth_factors, growth_factors$MR.delta > 0)
growth_factorsb <- subset(growth_factors, growth_factors$aLRL.delta > 0)
dim(growth_factors)
## [1] 173   6
dim(growth_factorsb)
## [1] 145   6
growth_factors
unique(growth_factors$condition)
## [1] "3" "2" "1"
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, "20231120_duf_wrky_2xko_growth_factors.csv", row.names = FALSE)
### Visualizing the growth factors:

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


growth_factors2$MR.delta <- as.numeric(as.character(growth_factors2$MR.delta))
growth_factors2$aLRL.delta <- as.numeric(as.character(growth_factors2$aLRL.delta))
growth_factors2$LRno.delta <- as.numeric(as.character(growth_factors2$LRno.delta))

growth_factors2$condition <- factor(growth_factors2$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 = "a", 
                                                              label = "p.signif")
MR.delta_p_geno <- MR.delta_p_geno + ggtitle("Main Root Growth")
MR.delta_p_geno

growth_factors2 <- subset(growth_factors2, growth_factors2$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 = "a", 
                                                              label = "p.signif")
LRno.delta_p_geno <- LRno.delta_p_geno + ggtitle("Lateral Root Increase")
LRno.delta_p_geno

growth_factors2 <- subset(growth_factors2, growth_factors2$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 = "a", 
                                                              label = "p.signif")
aLRL.delta_p_geno2

library(doBy)
avg_growth <- summaryBy(data = growth_factors2, MR.delta + aLRL.delta + LRno.delta ~ genotype + condition)
avg_growth
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
STI_growth_factors <- merge(growth_factors2, avg_growth_C, id="genotype")
STI_growth_factors
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)
odd <- subset(STI_growth_factors, STI_growth_factors$aLRL.STI < 0)
odd
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
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 = "a")
STI125_plot <- STI125_plot + scale_color_manual(values= c("blue", "plum", "rosybrown1", "hotpink","red"))
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 = "a") 
STI75_plot <- STI75_plot + scale_color_manual(values= c("blue", "plum", "rosybrown1", "hotpink","red"))
STI75_plot