Experimental description:
The seeds were sterilized in 50% bleach for 10 min with gentle shaking, then rinsed 5 times with deionized water, and incubated at 4 °C overnight. Then seeds were germinated in ¼ MS + sucrose (with vitamins) control plate for 4 complete days in growth chamber #35. 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 20 h light /4 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 in treatment plates.
getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20211222_ER_3_tomato_accessions_RSA_2nd_replicate"
setwd("C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20211222_ER_3_tomato_accessions_RSA_2nd_replicate/")
list.files(pattern = ".csv")
## [1] "2022jan_ER_growth_factors.csv" "ER_d5.csv"
## [3] "ER_d6.csv" "ER_d7.csv"
## [5] "ER_d8.csv" "ER_d9.csv"
my_data_d5<- read.csv("ER_d5.csv")
my_data_d6<- read.csv("ER_d6.csv")
my_data_d7<- read.csv("ER_d7.csv")
my_data_d8<- read.csv("ER_d8.csv")
my_data_d9<- read.csv("ER_d9.csv")
head(my_data_d5)
head(my_data_d6)
head(my_data_d9)
library(readr)
library(dplyr)
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
my_data_d5<-my_data_d5%>%mutate(day=5,my_data_d5 )
head(my_data_d5)
my_data_d6<-my_data_d6%>%mutate(day=6,my_data_d6 )
head(my_data_d6)
my_data_d7<-my_data_d7%>%mutate(day=7,my_data_d7 )
head(my_data_d7)
my_data_d8<-my_data_d8%>%mutate(day=8,my_data_d8 )
head(my_data_d8)
my_data_d9<-my_data_d9%>%mutate(day=9,my_data_d9 )
head(my_data_d9)
now I fuse days
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
colnames(Main_root)
## [1] "image" "root_name" "root"
## [4] "length" "vector_length" "surface"
## [7] "volume" "direction" "diameter"
## [10] "root_order" "root_ontology" "parent_name"
## [13] "parent" "insertion_position" "insertion_angle"
## [16] "n_child" "child_density" "first_child"
## [19] "insertion_first_child" "last_child" "insertion_last_child"
## [22] "day"
MR_data<-Main_root[, c(1:4, 10, 16,17, 19:22)]
MR_data
Lateral_root <- subset(all_data, all_data$root_order != 0)
Lateral_root
unique(Lateral_root$root_order)
## [1] 1 2
colnames(Lateral_root)
## [1] "image" "root_name" "root"
## [4] "length" "vector_length" "surface"
## [7] "volume" "direction" "diameter"
## [10] "root_order" "root_ontology" "parent_name"
## [13] "parent" "insertion_position" "insertion_angle"
## [16] "n_child" "child_density" "first_child"
## [19] "insertion_first_child" "last_child" "insertion_last_child"
## [22] "day"
MR_data2 <- subset(MR_data, MR_data$n_child >0)
MR_data2
I have started to confuse now….by [1], etc…..
temporary <- subset(Lateral_root, Lateral_root$parent == MR_data2$root[1])
temporary
temporary <- subset(temporary, temporary$day == MR_data2$day[3])
temporary
again what is dim here tries to do????
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.1701659
LR_number
## [1] 1
MR_data2$LRL <- 0
MR_data2$LRno <- 0
MR_data2
MR_data2$LRL[1] <- total_LRL
MR_data2$LRno[1] <- LR_number
MR_data2
MR_data_noChild <- subset(MR_data, MR_data$n_child == 0)
MR_data_noChild
length(MR_data2$root)
## [1] 182
for(i in 2:182){
temporary <- subset(Lateral_root, Lateral_root$parent == MR_data2$root[i])
temporary <- subset(temporary, temporary$day == MR_data2$day[i])
total_LRL <- sum(temporary$length)
LR_number <- dim(temporary)[1]
MR_data2$LRL[i] <- total_LRL
MR_data2$LRno[i] <- LR_number
}
MR_data2$check <- MR_data2$n_child - MR_data2$LRno
MR_data2
unique(MR_data2$check)
## [1] 0
#let's remove check column from Child:
MR_data_Child2 <- MR_data2[,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
MR_all$root_name[1]
## [1] " pl3_c_a"
text <- strsplit(x = MR_all$root_name[1], split = "_")
text
## [[1]]
## [1] " pl3" "c" "a"
plate <- text[[1]][1]
plate <- gsub(" ", "", plate)
plate
## [1] "pl3"
cond <- text[[1]][2]
cond
## [1] "c"
genotype <- text[[1]][3]
genotype
## [1] "a"
dim(MR_all)
## [1] 240 13
for(i in 1:240){
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
length(unique(MR_all$root_name))
## [1] 49
unique(MR_all$day)
## [1] 5 6 7 8 9
MR_all$day <- as.numeric(as.character(MR_all$day))
unique(MR_all$day)
## [1] 5 6 7 8 9
240/49
## [1] 4.897959
library(ggplot2)
library(ggpubr)
unique(MR_all$condition)
## [1] "c" "s" "0"
histogram_TRS_ER <- ggdensity(MR_all, x = "TRS",
add = "mean", rug = TRUE, facet.by = "day",
color = "condition", fill = "condition")
#palette = c("#00AFBB", "#E7B800"))
histogram_TRS_ER
## Warning: Groups with fewer than two data points have been dropped.
## Groups with fewer than two data points have been dropped.
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
## Warning in max(ids, na.rm = TRUE): no non-missing arguments to max; returning
## -Inf
we have a problem - since we have three conditions - rather than two expected ones - lets check what are the roots with “0” as condition (should be “c”):
MR_all2 <- subset(MR_all, MR_all$condition != "0")
MR_funky <- subset(MR_all, MR_all$condition == "0")
MR_funky
# wait - but they also have no genotype - let's check images
temp <- subset(MR_all2, MR_all2$image == "EricsTomatoLines_202102015011.rsml")
temp
temp2 <- subset(MR_all2, MR_all2$image == "EricsTomatoLines_202102017025.rsml")
temp2
conclusion >> the two roots with condition “0” are actually lateral roots that were not classified as such (whatever reason) - and thus we can ignore them and remove them from the dataset - continuing with MR_all2 onwards
MR_all2$TRS <- as.numeric(MR_all2$TRS)
unique(MR_all2$condition)
## [1] "c" "s"
histogram_TRS_ER <- ggdensity(MR_all2, x = "TRS",
add = "mean", rug = TRUE, facet.by = "day",
color = "condition", fill = "condition",
palette = c("#00AFBB", "#E7B800"))
histogram_TRS_ER
histogram_LRno_ER <- ggdensity(MR_all2, x = "LRno",
add = "mean", rug = TRUE, facet.by = "day",
color = "condition", fill = "condition",
palette = c("#00AFBB", "#E7B800"))
histogram_LRno_ER
pdf("histogram.TRS.ER.pdf")
plot(histogram_TRS_ER)
# if plotting multiple graphs - this command is extremely important
dev.off()
## png
## 2
pdf("histogram.LRno.ER.pdf")
plot(histogram_LRno_ER)
# 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_ER, histogram_LRno_ER, ncol=2,
align = "hv", labels=c("AUTO"),
label_size = 24)
dev.off()
## png
## 2
unique(MR_all2$genotype)
## [1] "a" "b" "c"
TRS_lgraph <- ggplot(data=MR_all2, aes(x= genotype, y=TRS, color = condition))
TRS_lgraph <- TRS_lgraph + geom_boxplot()
TRS_lgraph <- TRS_lgraph + facet_grid(day ~ condition, scales = "free") + scale_color_manual(values=c("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)
better_TRS_graph <- ggplot(data=MR_all2, aes(x= genotype, y=TRS, color = genotype))
better_TRS_graph <- better_TRS_graph + geom_beeswarm(alpha=0.6, priority = "density")
better_TRS_graph <- better_TRS_graph + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
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
#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_all2)
my_graph <- ggplot(data=MR_all2, 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)
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_all2, MR_all2$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)
There is one pant that didnot grow but then stops to grow (pl6_c_a), so let;s remove them:
geno_a$root_name <- gsub(" ", "", geno_a$root_name)
bye_a <- c("pl6_c_a")
OK_geno_a <- subset(geno_a, !(geno_a$root_name %in% bye_a))
unique(OK_geno_a$root_name)
## [1] "pl3_c_a" "pl4_c_a" "pl4_s_a" "pl8_c_a" "pl1_c_a" "pl2_s_a" "pl2_c_a"
## [8] "pl5_s_a" "pl5_c_a" "pl7_s_a" "pl7_c_a" "pl1_s_a" "pl3_s_a" "pl6_s_a"
## [15] "pl8_s_a"
my_graph <- ggplot(data=OK_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 (cm)") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph
my_graph <- ggplot(data=OK_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=OK_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=OK_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
so let’s move to the next genotype:
geno_b <- subset(MR_all2, MR_all2$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)
Most of the above look good except plant “pl2_c_b” and “pl4_c_b” and “pl1_s_b”so let’s remove them:
geno_b$root_name <- gsub(" ", "", geno_b$root_name)
bye_b <- c("pl2_c_b","pl4_c_b","pl1_s_b")
OK_geno_b <- subset(geno_b, !(geno_b$root_name %in% bye_b))
unique(OK_geno_b$root_name)
## [1] "pl2_s_b" "pl3_s_b" "pl1_c_b" "pl3_c_b" "pl4_s_b" "pl5_c_b" "pl6_s_b"
## [8] "pl6_c_b" "pl7_c_b" "pl8_s_b" "pl8_c_b" "pl5_s_b" "pl7_s_b"
my_graph <- ggplot(data=OK_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 (cm)") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph
my_graph <- ggplot(data=OK_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=OK_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=OK_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_all2, MR_all2$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
library(plotly)
ggplotly(my_graph)
There is one pant that didnot grow at all (pl8_s_c), and one plant that grows very odd (pl5_s_c), so let;s remove them:
geno_c$root_name <- gsub(" ", "", geno_c$root_name)
bye_c <- c("pl8_s_c", "pl5_s_c")
OK_geno_c <- subset(geno_c, !(geno_c$root_name %in% bye_c))
unique(OK_geno_c$root_name)
## [1] "pl1_c_c" "pl2_c_c" "pl3_c_c" "pl4_c_c" "pl5_c_c" "pl6_c_c" "pl7_c_c"
## [8] "pl8_c_c" "pl3_s_c" "pl4_s_c" "pl7_s_c" "pl1_s_c" "pl6_s_c" "pl2_s_c"
my_graph <- ggplot(data=OK_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 (cm)") + xlab("time (day)") + ggtitle("Main Root Length") + theme(legend.position='none')
my_graph
my_graph <- ggplot(data=OK_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=OK_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=OK_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
MR_OK <- rbind(OK_geno_a, OK_geno_b)
MR_OK <- rbind(MR_OK, OK_geno_c)
MR_time_graph <- ggplot(data=MR_OK, 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("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 <- MR_time_graph + scale_y_continuous(name = " length (cm)", breaks = c(0:10))
MR_time_graph <- MR_time_graph + scale_color_manual(values=c("blue", "dark orange", "red"))
MR_time_graph
LRL_time_graph <- ggplot(data=MR_OK, aes(x= day, y=LRL, group = root_name, color = genotype))
LRL_time_graph <- LRL_time_graph + geom_line(alpha = 0.2)
LRL_time_graph <- LRL_time_graph + facet_grid(~ condition)
LRL_time_graph <- LRL_time_graph + ylab("length (cm)") + xlab("time (days after germination)") + ggtitle("Lateral Root Length")
LRL_time_graph <- LRL_time_graph + stat_summary(fun.y=mean, aes(group= genotype), size=0.7, geom="line", linetype = "dashed")
LRL_time_graph
LRno_time_graph <- ggplot(data=MR_OK, 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_OK, 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 52 rows containing non-finite values (`stat_summary()`).
## Warning: Removed 52 rows containing missing values (`geom_line()`).
TRS_time_graph <- ggplot(data=MR_OK, 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 <- TRS_time_graph + scale_color_manual(values=c("blue", "dark orange", "red"))
TRS_time_graph
### 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_OK, MR_OK$root_name == unique(MR_OK$root_name)[1])
temp2 <- temp1[order(temp1$day),]
temp2
# 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
temp_MR2 <- subset(temp_MR, temp_MR$MRdouble == "no")
temp_MR2 <- subset(temp_MR2, temp_MR2$day < 8)
temp_MR2
plot(temp_MR2$length~ temp_MR2$day)
# let's add the regression line to this graph
abline(lm(temp_MR2$length~ temp_MR2$day))
MR_model <- lm(temp_MR2$length~ temp_MR2$day)
MR_model
##
## Call:
## lm(formula = temp_MR2$length ~ temp_MR2$day)
##
## Coefficients:
## (Intercept) temp_MR2$day
## -6.317 2.163
summary(MR_model)
##
## Call:
## lm(formula = temp_MR2$length ~ temp_MR2$day)
##
## Residuals:
## 1 2 3
## -0.002641 0.005281 -0.002641
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -6.316808 0.027696 -228.1 0.00279 **
## temp_MR2$day 2.162885 0.004574 472.9 0.00135 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.006468 on 1 degrees of freedom
## Multiple R-squared: 1, Adjusted R-squared: 1
## F-statistic: 2.236e+05 on 1 and 1 DF, p-value: 0.001346
MR_growth_rate <- MR_model$coefficients[[2]]
# Let's remove all NA for LRno and aLRL in temp2
LR_temp <- temp2[,c("day", "LRno", "aLRL")]
LR_temp2 <- na.omit(LR_temp)
# Let's start with LRno
plot(LR_temp2$LRno ~ LR_temp2$day)
abline(lm(LR_temp2$LRno ~ LR_temp2$day))
LRno_model <- lm(LR_temp2$LRno ~ LR_temp2$day)
LRno_model
##
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
##
## Coefficients:
## (Intercept) LR_temp2$day
## -22.4 4.6
summary(LRno_model)
##
## Call:
## lm(formula = LR_temp2$LRno ~ LR_temp2$day)
##
## Residuals:
## 1 2 3 4 5
## 0.4 -2.2 2.2 0.6 -1.0
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -22.400 4.364 -5.134 0.01432 *
## LR_temp2$day 4.600 0.611 7.529 0.00486 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.932 on 3 degrees of freedom
## Multiple R-squared: 0.9497, Adjusted R-squared: 0.933
## F-statistic: 56.68 on 1 and 3 DF, p-value: 0.004858
LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase
## [1] 4.6
# 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
## -1.4744 0.3326
summary(aLRL_model)
##
## Call:
## lm(formula = LR_temp2$aLRL ~ LR_temp2$day)
##
## Residuals:
## 1 2 3 4 5
## -0.01834 0.20510 -0.24279 -0.05636 0.11239
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.47440 0.44628 -3.304 0.0456 *
## LR_temp2$day 0.33258 0.06249 5.322 0.0130 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1976 on 3 degrees of freedom
## Multiple R-squared: 0.9042, Adjusted R-squared: 0.8723
## F-statistic: 28.32 on 1 and 3 DF, p-value: 0.01296
aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth
## [1] 0.3325816
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] <- temp2$genotype[1]
growth_factors[1,3] <- 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
length(unique(MR_OK$root_name))
## [1] 42
# 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:42){
temp1 <- subset(MR_OK, MR_OK$root_name == unique(MR_OK$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 <- subset(temp_MR2, temp_MR2$day < 8)
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] <- temp2$genotype[1]
growth_factors[e,3] <- 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
write.csv(growth_factors, "2022jan_ER_growth_factors.csv", row.names = FALSE)
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 = "genotype", fill="genotype", color="genotype",
facet.by = c("cond"), ncol=2,
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") + stat_compare_means(method="t.test", ref.group = "a",
label = "p.signif", hide.ns = T)
MR.delta_p_geno <- MR.delta_p_geno + ggtitle("Main Root Growth")
MR.delta_p_geno <- MR.delta_p_geno + scale_color_manual(values=c("blue", "dark orange", "red"))
MR.delta_p_geno
LRno.delta_p_geno <- ggerrorplot(growth_factors, y = "LRno.delta", x = "genotype", fill="genotype", color="genotype",
facet.by = c("cond"), ncol=2,
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") + stat_compare_means(method="t.test", ref.group = "a",
label = "p.signif", hide.ns = T)
LRno.delta_p_geno <- LRno.delta_p_geno + ggtitle("Lateral Root Number Increase")
LRno.delta_p_geno
aLRL.delta_p_geno <- ggerrorplot(growth_factors, y = "aLRL.delta", x = "genotype", fill="genotype", color="genotype",
facet.by = c("cond"), ncol=2,
desc_stat = "mean_sd", add = "jitter",
add.params = list(color = "darkgray"),
xlab="Genotype", ylab="Growth Rate (cm / day)")
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "a",
label = "p.signif", hide.ns = T)
aLRL.delta_p_geno <- aLRL.delta_p_geno + ggtitle("average Lateral Root Growth")
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)
##
## Attaching package: 'doBy'
## The following object is masked from 'package:dplyr':
##
## order_by
avg_growth <- summaryBy(data = growth_factors, MR.delta + aLRL.delta + LRno.delta ~ genotype + cond)
avg_growth
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
STI_growth_factors <- merge(growth_factors, 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)
STI <- subset(STI_growth_factors, STI_growth_factors$cond == "s")
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_compare_means(label = "p.signif", method = "t.test", ref.group = "a")
MR.STI_plot <- MR.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
MR.STI_plot
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("Fraction of control") + ggtitle("Salt Tolerance Index based on average Lateral Root Growth")
aLRL.STI_plot <- aLRL.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "a")
aLRL.STI_plot <- aLRL.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
aLRL.STI_plot
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("Fraction of control") + ggtitle("Salt Tolerance Index based on average Increase in Lateral Root Number")
LRno.STI_plot <- LRno.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "a")
LRno.STI_plot <- LRno.STI_plot + scale_color_manual(values=c("royalblue", "coral3", "tomato4"))
LRno.STI_plot
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
library(cowplot)
plot_grid(MR_time_graph, LRL_time_graph, LRno_time_graph, labels = c("AUTO"), ncol = 2)
pdf("20220114_ER_3accessions_RSA_analysis_Fig1.pdf", width = 13, height = 10)
plot_grid(MR_time_graph, LRL_time_graph, LRno_time_graph, labels = c("AUTO"), ncol = 2)
dev.off()
## png
## 2
library(cowplot)
plot_grid(MR.delta_p_geno, aLRL.delta_p_geno, LRno.delta_p_geno, labels = c("AUTO"), ncol = 2)
pdf("20220114_ER_3accessions_RSA_analysis_Fig2.pdf", width = 13, height = 10)
plot_grid(MR.delta_p_geno, aLRL.delta_p_geno, LRno.delta_p_geno, labels = c("AUTO"), ncol = 2)
dev.off()
## png
## 2
library(cowplot)
plot_grid(MR.STI_plot, aLRL.STI_plot, LRno.STI_plot, labels = c("AUTO"), ncol = 2)
pdf("20220114_ER_3accessions_RSA_analysis_Fig3.pdf", width = 13, height = 10)
plot_grid(MR.STI_plot, aLRL.STI_plot, LRno.STI_plot, labels = c("AUTO"), ncol = 2)
dev.off()
## png
## 2
pdf("20230128_ER_3accessions_TRS_Fig4.pdf", width = 13, height = 10)
plot_grid(TRS_time_graph, aLRL.STI_plot, LRno.STI_plot, labels = c("AUTO"), ncol = 2)
dev.off()
## png
## 2
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).
The preview shows you a rendered HTML copy of the contents of the editor. Consequently, unlike Knit, Preview does not run any R code chunks. Instead, the output of the chunk when it was last run in the editor is displayed.