This is RSA for the genome paper in which I try to test how does H2O2 and asc treatments impact the RSA of the two accessions under salt stress.

getwd()
## [1] "C:/Users/Julkowska Lab/Desktop/R codes by Maryam/20241212_RSA_LA1511_M248_ASC_H2O2_treatment"
list.files(pattern = ".csv")
## [1] "2024Dec_2acc_ASC_H2O2_growth_factors.csv"
## [2] "D5-ROS.csv"                              
## [3] "D9-ROS.csv"
d5<-read.csv("D5-ROS.csv")
d9<-read.csv("D9-ROS.csv")
head(d5)
d5$day <- 5
d9$day <- 9
head(d9)
all_data <- rbind(d5, d9)
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"
MR_data<-Main_root[, c(1:4, 10, 14:17, 19:22)]
MR_data
Lateral_root <- subset(all_data, all_data$root_order != 0)
Lateral_root
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
temporary <- subset(Lateral_root, Lateral_root$parent == MR_data2$root[1])
temporary
temporary <- subset(temporary, temporary$day == MR_data2$day[3])
temporary
dim(temporary)
## [1]  7 22
total_LRL <- sum(temporary$length)
angle_LR.avg <- mean(temporary$insertion_angle)
angle_LR.sd <- sd(temporary$insertion_angle)
LR_number <- dim(temporary)[1]
total_LRL
## [1] 3.725968
angle_LR.avg
## [1] 73.12841
angle_LR.sd
## [1] 29.20661
LR_number
## [1] 7
MR_data2$LRL <- 0
MR_data2$LRno <- 0
MR_data2$LRangle.avg <- 0
MR_data2$LRangle.sd <- 0
MR_data2
MR_data2$LRL[1] <- total_LRL
MR_data2$LRno[1] <- LR_number
MR_data2$LRangle.avg[1] <- angle_LR.avg
MR_data2$LRangle.sd[1] <- angle_LR.sd
MR_data2
MR_data_noChild <- subset(MR_data, MR_data$n_child == 0)
MR_data_noChild
length(MR_data2$root)
## [1] 119
for(i in 2:119){
  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)
  angle_LR.avg <- mean(temporary$insertion_angle)
  angle_LR.sd <- sd(temporary$insertion_angle)
  LR_number <- dim(temporary)[1]
  MR_data2$LRL[i] <- total_LRL
  MR_data2$LRno[i] <- LR_number
  MR_data2$LRangle.avg[i] <- angle_LR.avg
  MR_data2$LRangle.sd[i] <- angle_LR.sd
    
}

MR_data2$check <- MR_data2$n_child - MR_data2$LRno
MR_data2
unique(MR_data2$check)
## [1] 0
colnames(MR_data2)
##  [1] "image"                 "root_name"             "root"                 
##  [4] "length"                "root_order"            "insertion_position"   
##  [7] "insertion_angle"       "n_child"               "child_density"        
## [10] "insertion_first_child" "last_child"            "insertion_last_child" 
## [13] "day"                   "LRL"                   "LRno"                 
## [16] "LRangle.avg"           "LRangle.sd"            "check"
MR_data_Child2 <- MR_data2[,1:17]
MR_data_Child2
MR_data_noChild
MR_data_noChild$LRL <- 0
MR_data_noChild$LRno <- 0
MR_data_noChild$LRangle.avg <- 0
MR_data_noChild$LRangle.sd <- 0
MR_all <-rbind(MR_data_Child2, MR_data_noChild)
MR_all
dim(MR_data_Child2)
## [1] 119  17
dim(MR_data_noChild)
## [1] 130  17
dim(MR_all)
## [1] 249  17
MR_all$root_name[50]
## [1] " PL1-C-ROS-M248"
text <- strsplit(x = MR_all$root_name[50], split = "-")
text
## [[1]]
## [1] " PL1" "C"    "ROS"  "M248"
plate <- text[[1]][1]
plate <- gsub(" ", "", plate)
plate
## [1] "PL1"
cond <- text[[1]][2:3]
cond
## [1] "C"   "ROS"
genotype <- text[[1]][4]
genotype
## [1] "M248"
dim(MR_all)
## [1] 249  17
for(i in 1:249){
  text <- strsplit(x = MR_all$root_name[i], split = "-")
  plate <- text[[1]][1]
  cond <- paste(text[[1]][2], text[[1]][3], sep="_")
  genotype <- text[[1]][4]
   
  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

for(i in 1:nrow(MR_all)){
  if(MR_all$LRno[i] == 0){
    MR_all$aLRL[i] <- 0
  } else {
  MR_all$aLRL[i] <- MR_all$LRL[i]/ MR_all$LRno[i]  
  }
}


MR_all$MRpLRL <- MR_all$length / MR_all$LRL
MR_all
length(unique(MR_all$root_name))
## [1] 123
unique(MR_all$day)
## [1] 9 5
MR_all$day <- as.numeric(as.character(MR_all$day))
unique(MR_all$day)
## [1] 9 5
249/123
## [1] 2.02439
unique(MR_all$root_name)
##   [1] " PL1-C-ASC-M248"   " PL2-C-ASC-M248"   " PL3-C-ASC-M248"  
##   [4] " PL4-C-ASC-M248"   " PL5-C-ASC-M248"   " PL6-C-ASC-M248"  
##   [7] " PL7-C-ASC-M248"   " PL8-C-ASC-M248"   " PL9-C-ASC-M248"  
##  [10] " PL10-C-ASC-M248"  " PL1-C-ASC-LA"     " PL2-C-ASC-LA"    
##  [13] " PL3-C-ASC-LA"     " PL4-C-ASC-LA"     " PL5-C-ASC-LA"    
##  [16] " PL6-C-ASC-LA"     " PL7-C-ASC-LA"     " PL8-C-ASC-LA"    
##  [19] " PL9-C-ASC-LA"     " PL10-C-ASC-LA"    " PL1-S-ASC-M248"  
##  [22] " PL2-S-ASC-M248"   " PL3-S-ASC-M248"   " PL4-S-ASC-M248"  
##  [25] " PL5-S-ASC-M248"   " PL6-S-ASC-M248"   " PL7-S-ASC-M248"  
##  [28] " PL8-S-ASC-M248"   " PL9-S-ASC-M248"   " PL10-S-ASC-M248" 
##  [31] " PL1-S-ASC-LA"     " PL2-S-ASC-LA"     " PL3-S-ASC-LA"    
##  [34] " PL4-S-ASC-LA"     " PL5-S-ASC-LA"     " PL6-S-ASC-LA"    
##  [37] " PL7-S-ASC-LA"     " PL8-S-ASC-LA"     " PL9-S-ASC-LA"    
##  [40] " PL10-S-ASC-LA"    " PL10-C-ROS-M248"  " PL9-C-ROS-M248"  
##  [43] " PL8-C-ROS-M248"   " PL7-C-ROS-M248"   " PL6-C-ROS-M248"  
##  [46] " PL5-C-ROS-M248"   " PL4-C-ROS-M248"   " PL3-C-ROS-M248"  
##  [49] " PL2-C-ROS-M248"   " PL1-C-ROS-M248"   " PL11-C-ROS-LA"   
##  [52] " PL10-C-ROS-LA"    " PL9-C-ROS-LA"     " PL8-C-ROS-LA"    
##  [55] " PL7-C-ROS-LA"     " PL6-C-ROS-LA"     " PL5-C-ROS-LA"    
##  [58] " PL4-C-ROS-LA"     " PL3-C-ROS-LA"     " PL2-C-ROS-LA"    
##  [61] " PL1-C-ROS-LA"     " PL10-S-ROS-LA"    " PL9-S-ROS-LA"    
##  [64] " PL8-S-ROS-LA"     " PL7-S-ROS-LA"     " PL5-S-ROS-LA"    
##  [67] " PL4-S-ROS-LA"     " PL3-S-ROS-LA"     " PL1-S-ROS-LA"    
##  [70] " PL11-S-ROS-M248"  " PL10-S-ROS-M248"  " PL9-S-ROS-M248"  
##  [73] " PL8-S-ROS-M248"   " PL7-S-ROS-M248"   " PL6-S-ROS-M248"  
##  [76] " PL5-S-ROS-M248"   " PL4-S-ROS-M248"   " PL3-S-ROS-M248"  
##  [79] " PL2-S-ROS-M248"   " PL1-S-ROS-M248"   " PL10-C-ONLY-LA"  
##  [82] " PL9-C-ONLY-LA"    " PL8-C-ONLY-LA"    " PL7-C-ONLY-LA"   
##  [85] " PL6-C-ONLY-LA"    " PL5-C-ONLY-LA"    " PL4-C-ONLY-LA"   
##  [88] " PL3-C-ONLY-LA"    " PL2-C-ONLY-LA"    " PL1-C-ONLY-LA"   
##  [91] " PL10-C-ONLY-M248" " PL9-C-ONLY-M248"  " PL8-C-ONLY-M248" 
##  [94] " PL7-C-ONLY-M248"  " PL6-C-ONLY-M248"  " PL5-C-ONLY-M248" 
##  [97] " PL4-C-ONLY-M248"  " PL3-C-ONLY-M248"  " PL2-C-ONLY-M248" 
## [100] " PL1-C-ONLY-M248"  " PL10-S-ONLY-M248" " PL9-S-ONLY-M248" 
## [103] " PL8-S-ONLY-M248"  " PL7-S-ONLY-M248"  " PL6-S-ONLY-M248" 
## [106] " PL5-S-ONLY-M248"  " PL4-S-ONLY-M248"  " PL3-S-ONLY-M248" 
## [109] " PL2-S-ONLY-M248"  " PL1-S-ONLY-M248"  " PL10-S-ONLY-LA"  
## [112] " PL9-S-ONLY-LA"    " PL8-S-ONLY-LA"    " PL7-S-ONLY-LA"   
## [115] " PL6-S-ONLY-LA"    " PL5-S-ONLY-LA"    " PL4-S-ONLY-LA"   
## [118] " PL3-S-ONLY-LA"    " PL2-S-ONLY-LA"    " PL6-S-ROS-LA"    
## [121] " PL2-S-ROS-LA"     " PL1-S-ONLY-LA"    " root_0"
unique(MR_all$genotype)
## [1] "M248" "LA"   NA

##it seems that I have one main root that does not have name, and I have NA in genotype as well…….not sure I got ride of NA in my genotype?????????/

good_stuff <- c("M248", "LA")
funk <- subset(MR_all, !(MR_all$genotype %in% good_stuff))
funk
unique(MR_all$genotype)
## [1] "M248" "LA"   NA
MR_all <- subset(MR_all, (MR_all$genotype %in% good_stuff))
MR_d9 <- subset(MR_all, MR_all$day == 5)
MR_d5 <- subset(MR_all, MR_all$day == 9)

suspish <- subset(MR_d9, !(MR_d9$root_name %in% MR_d5$root_name))
suspish
MR_clean <- subset(MR_all, (MR_all$root_name %in% MR_d5$root_name))

MR_d9 <- subset(MR_clean, MR_clean$day == 5)
MR_d5 <- subset(MR_clean, MR_clean$day == 9)

suspish <- subset(MR_d9, !(MR_d9$root_name %in% MR_d5$root_name))
suspish
suspish <- subset(MR_d9, !(MR_d5$root_name %in% MR_d9$root_name))
suspish
unique(MR_all$cond)
## [1] "C_ASC"  "S_ASC"  "C_ROS"  "S_ROS"  "C_ONLY" "S_ONLY"
library(ggplot2)
library(ggpubr)
MR_all$condition<- factor(MR_all$condition, levels=c("C_ONLY", "C_ASC", "C_ROS",  "S_ONLY", "S_ASC" , "S_ROS" ))


histogram_TRS_all <- ggdensity(MR_all, x = "TRS",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition")
                           

histogram_TRS_all

pdf("histogram.TRS.all.treatment.pdf")
plot(histogram_TRS_all)
# if plotting multiple graphs - this command is extremely important 
dev.off()
## png 
##   2
histogram_LRno_all <- ggdensity(MR_all, x = "LRno",
                           add = "mean", rug = TRUE, facet.by = "day",
                           color = "condition", fill = "condition")
                           

histogram_LRno_all

pdf("histogram.LRno.all.treatment.pdf")
plot(histogram_LRno_all)
# 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_Histogram.pdf", height = 15, width = 12)
plot_grid(histogram_TRS_all, histogram_LRno_all, ncol=2,
          align = "hv", labels=c("AUTO"), 
          label_size = 24)
dev.off()
## png 
##   2
unique(MR_all$genotype)
## [1] "M248" "LA"
TRS_lgraph <- ggplot(data=MR_all, aes(x= genotype, y=TRS, color = condition)) 
TRS_lgraph <- TRS_lgraph + geom_boxplot()
TRS_lgraph <- TRS_lgraph + facet_grid(day ~ condition, scales = "free") + scale_color_manual(values=c("blue","turquoise3", "cyan","red", "maroon3", "dark orange"))
TRS_lgraph <- TRS_lgraph + ylab("Total root size (cm)") + xlab("Genotype") + theme(legend.position='none')
TRS_lgraph

#I would like to have a different test, tukey within each accession

unique(MR_all$day)
## [1] 9 5
TRS_data <- subset(MR_all, MR_all$day == 9)
TRS_data
TRS_graph_d9 <- ggerrorplot(TRS_data, y = "TRS", x = "condition", fill="condition", color="condition", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "genotype", palette = "jco",
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Total Root Size (cm)") 
TRS_graph_d9 <- TRS_graph_d9 + rremove("legend") + stat_compare_means(method="t.test", ref.group = "C_ONLY", 
                                                              label = "p.signif", hide.ns = T)
TRS_graph_d9 <- TRS_graph_d9 + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
TRS_graph_d9

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
my_graph <- ggplot(data=MR_all, aes(x= day, y=length, group = root_name, color = condition)) 
my_graph <- my_graph + geom_line(alpha = 0.2) 
my_graph <- my_graph + facet_grid(~ genotype) 
my_graph <- my_graph + ylab("Main Root Length (cm)") + xlab("time (day)") 
my_graph

#IT Seems odd…why replicates are either cluster below average or above average? in theory, there should be variaotion across the scale.

my_graph <- ggerrorplot(MR_all, y = "length", x = "condition", fill="cond", color="condition", 
                        desc_stat = "mean_sd", add = "jitter", facet.by = "genotype", palette = "jco",
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Main Root Length (cm)") 
my_graph <- my_graph + rremove("legend") + stat_compare_means(method="t.test", ref.group = "C_ONLY", 
                                                              label = "p.signif", hide.ns = T)
my_graph <- my_graph + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
my_graph

#I skipped curating data…and jumped into growth factor calculations

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

MR_model <- lm(temp_MR$length~ temp_MR$day)
MR_model
## 
## Call:
## lm(formula = temp_MR$length ~ temp_MR$day)
## 
## Coefficients:
## (Intercept)  temp_MR$day  
##      -4.128        1.154
MR_growth_rate <- MR_model$coefficients[[2]]
MR_growth_rate
## [1] 1.154348
#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  
##        -8.75          1.75
LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))
LRno_increase
## [1] 1.75
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.6654        0.1331
aLRL_growth <- as.numeric(as.character(aLRL_model$coefficients[[2]]))
aLRL_growth
## [1] 0.1330703
temp2
names <- c(text="root_name", "genotype", "cond.final", "MR.delta", "LRno.delta", "aLRL.delta")
growth_factors <- data.frame()

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

growth_factors[1,1] <- temp2$root_name[1]
growth_factors[1,2] <- as.character(temp2$genotype[1])
growth_factors[1,3] <- as.character(temp2$condition[1])
if(MR_growth_rate < 0){
  MR_growth_rate <- 0}
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$root_name))
## [1] 122

growth factor calculations - loop

# 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:122){

  temp1 <- subset(MR_all, MR_all$root_name == unique(MR_all$root_name)[e])
  temp2 <- temp1[order(temp1$day),]
  temp2

############ MR calculations
  temp_MR <- temp2[,c("day", "length")]
  if(dim(temp_MR)[1] > 1){
    MR_model <- lm(temp_MR$length~ temp_MR$day)
    MR_growth_rate <- MR_model$coefficients[[2]]} else{MR_growth_rate <- 0}
  if(MR_growth_rate < 0){MR_growth_rate <- 0}
    
  ############ LRno calculations
  LR_temp <- temp2[,c("day", "LRno", "aLRL")]
  LR_temp$aLRL <- as.numeric(as.character(LR_temp$aLRL))
  #LR_temp2 <- na.omit(LR_temp)
  dim(LR_temp)
  
  ####################### safety precaution to calculate LR growth rate only for the plants that have LR at #least for two days: 
  if(dim(LR_temp)[1] > 1){
  LRno_model <- lm(LR_temp$LRno ~ LR_temp$day)
  LRno_increase <- as.numeric(as.character(LRno_model$coefficients[[2]]))

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


write.csv(growth_factors, "2024Dec_2acc_ASC_H2O2_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 = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Main root growth rate (cm / day)") 
MR.delta_p_geno <- MR.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "C_ONLY", 
                                                              label = "p.signif", hide.ns = T)
MR.delta_p_geno

GF_M248<- subset(growth_factors, growth_factors$genotype == "M248")
GF_LA<- subset(growth_factors, growth_factors$genotype == "LA")


aov(MR.delta ~ cond.final, data = GF_M248)
## Call:
##    aov(formula = MR.delta ~ cond.final, data = GF_M248)
## 
## Terms:
##                 cond.final Residuals
## Sum of Squares   0.5399151 2.5174095
## Deg. of Freedom          5        55
## 
## Residual standard error: 0.2139418
## Estimated effects may be unbalanced
Output <- TukeyHSD(aov(MR.delta ~ cond.final, data = GF_M248))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = MR.delta ~ cond.final, data = GF_M248)
## 
## $cond.final
##                       diff         lwr         upr     p adj
## C_ASC-C_ONLY   0.041306852 -0.24118475  0.32379846 0.9979929
## C_ROS-C_ONLY   0.177217030 -0.10527458  0.45970864 0.4419330
## S_ONLY-C_ONLY -0.027633365 -0.31012497  0.25485824 0.9997140
## S_ASC-C_ONLY  -0.138982428 -0.42147403  0.14350918 0.6950823
## S_ROS-C_ONLY   0.048911770 -0.22708491  0.32490845 0.9950164
## C_ROS-C_ASC    0.135910178 -0.14658143  0.41840178 0.7145828
## S_ONLY-C_ASC  -0.068940217 -0.35143182  0.21355139 0.9786104
## S_ASC-C_ASC   -0.180289280 -0.46278089  0.10220233 0.4224852
## S_ROS-C_ASC    0.007604918 -0.26839176  0.28360160 0.9999995
## S_ONLY-C_ROS  -0.204850395 -0.48734200  0.07764121 0.2821409
## S_ASC-C_ROS   -0.316199458 -0.59869106 -0.03370785 0.0197707
## S_ROS-C_ROS   -0.128305260 -0.40430194  0.14769142 0.7429283
## S_ASC-S_ONLY  -0.111349063 -0.39384067  0.17114254 0.8518535
## S_ROS-S_ONLY   0.076545135 -0.19945154  0.35254181 0.9628154
## S_ROS-S_ASC    0.187894198 -0.08810248  0.46389088 0.3500509
Pval = Output$cond.final[,'p adj']
Pval
##  C_ASC-C_ONLY  C_ROS-C_ONLY S_ONLY-C_ONLY  S_ASC-C_ONLY  S_ROS-C_ONLY 
##    0.99799293    0.44193296    0.99971400    0.69508231    0.99501644 
##   C_ROS-C_ASC  S_ONLY-C_ASC   S_ASC-C_ASC   S_ROS-C_ASC  S_ONLY-C_ROS 
##    0.71458283    0.97861036    0.42248523    0.99999947    0.28214088 
##   S_ASC-C_ROS   S_ROS-C_ROS  S_ASC-S_ONLY  S_ROS-S_ONLY   S_ROS-S_ASC 
##    0.01977066    0.74292826    0.85185346    0.96281536    0.35005094
library(multcompView)
## Warning: package 'multcompView' was built under R version 4.3.2
stat.test<- multcompLetters(Pval)
stat.test
##  C_ASC  C_ROS S_ONLY  S_ASC  S_ROS C_ONLY 
##   "ab"    "a"   "ab"    "b"   "ab"   "ab"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$genotype <- "M248"
test
# Now same for LA:
Output <- TukeyHSD(aov(MR.delta ~ cond.final, data = GF_LA))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = MR.delta ~ cond.final, data = GF_LA)
## 
## $cond.final
##                      diff        lwr         upr     p adj
## C_ASC-C_ONLY   0.10954643 -0.3873125  0.60640540 0.9863958
## C_ROS-C_ONLY   0.20493502 -0.2805004  0.69037042 0.8119676
## S_ONLY-C_ONLY -0.35187647 -0.8487354  0.14498250 0.3070870
## S_ASC-C_ONLY  -0.39998594 -0.8968449  0.09687303 0.1824065
## S_ROS-C_ONLY  -0.22147760 -0.7183366  0.27538137 0.7749364
## C_ROS-C_ASC    0.09538859 -0.3900468  0.58082399 0.9919505
## S_ONLY-C_ASC  -0.46142290 -0.9582819  0.03543607 0.0831067
## S_ASC-C_ASC   -0.50953237 -1.0063913 -0.01267340 0.0413264
## S_ROS-C_ASC   -0.33102403 -0.8278830  0.16583494 0.3740810
## S_ONLY-C_ROS  -0.55681148 -1.0422469 -0.07137609 0.0157392
## S_ASC-C_ROS   -0.60492096 -1.0903564 -0.11948556 0.0067275
## S_ROS-C_ROS   -0.42641262 -0.9118480  0.05902278 0.1161603
## S_ASC-S_ONLY  -0.04810948 -0.5449684  0.44874949 0.9997279
## S_ROS-S_ONLY   0.13039887 -0.3664601  0.62725783 0.9706351
## S_ROS-S_ASC    0.17850834 -0.3183506  0.67536731 0.8946026
Pval = Output$cond.final[,'p adj']
Pval
##  C_ASC-C_ONLY  C_ROS-C_ONLY S_ONLY-C_ONLY  S_ASC-C_ONLY  S_ROS-C_ONLY 
##   0.986395791   0.811967554   0.307087044   0.182406541   0.774936388 
##   C_ROS-C_ASC  S_ONLY-C_ASC   S_ASC-C_ASC   S_ROS-C_ASC  S_ONLY-C_ROS 
##   0.991950459   0.083106704   0.041326379   0.374081039   0.015739204 
##   S_ASC-C_ROS   S_ROS-C_ROS  S_ASC-S_ONLY  S_ROS-S_ONLY   S_ROS-S_ASC 
##   0.006727503   0.116160311   0.999727947   0.970635060   0.894602572
stat.test<- multcompLetters(Pval)
stat.test
##  C_ASC  C_ROS S_ONLY  S_ASC  S_ROS C_ONLY 
##   "ab"    "a"   "bc"    "c"  "abc"  "abc"
test2 <- as.data.frame(stat.test$Letters)
test2$group1 <- rownames(test)
test2$group2 <- rownames(test)
colnames(test2)[1] <- "Tukey"
test2
test2$genotype <- "LA"
test2
test3 <- rbind(test, test2)
test3
MR.delta_p_geno <- ggerrorplot(growth_factors, y = "MR.delta", x = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Main root growth rate (cm / day)") 
MR.delta_p_geno <- MR.delta_p_geno + rremove("legend") +stat_pvalue_manual(test3, label = "Tukey", y.position = 2)  + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
MR.delta_p_geno

Output <- TukeyHSD(aov(LRno.delta ~ cond.final, data = GF_M248))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = LRno.delta ~ cond.final, data = GF_M248)
## 
## $cond.final
##                diff        lwr         upr     p adj
## C_ASC-C_ONLY  -0.45 -1.2270986  0.32709856 0.5314348
## C_ROS-C_ONLY  -0.40 -1.1770986  0.37709856 0.6532531
## S_ONLY-C_ONLY -1.90 -2.6770986 -1.12290144 0.0000000
## S_ASC-C_ONLY  -1.30 -2.0770986 -0.52290144 0.0001085
## S_ROS-C_ONLY  -1.80 -2.5592318 -1.04076816 0.0000001
## C_ROS-C_ASC    0.05 -0.7270986  0.82709856 0.9999638
## S_ONLY-C_ASC  -1.45 -2.2270986 -0.67290144 0.0000142
## S_ASC-C_ASC   -0.85 -1.6270986 -0.07290144 0.0242916
## S_ROS-C_ASC   -1.35 -2.1092318 -0.59076816 0.0000361
## S_ONLY-C_ROS  -1.50 -2.2770986 -0.72290144 0.0000071
## S_ASC-C_ROS   -0.90 -1.6770986 -0.12290144 0.0143453
## S_ROS-C_ROS   -1.40 -2.1592318 -0.64076816 0.0000179
## S_ASC-S_ONLY   0.60 -0.1770986  1.37709856 0.2199215
## S_ROS-S_ONLY   0.10 -0.6592318  0.85923184 0.9987846
## S_ROS-S_ASC   -0.50 -1.2592318  0.25923184 0.3870697
Pval = Output$cond.final[,'p adj']
Pval
##  C_ASC-C_ONLY  C_ROS-C_ONLY S_ONLY-C_ONLY  S_ASC-C_ONLY  S_ROS-C_ONLY 
##  5.314348e-01  6.532531e-01  2.445801e-08  1.084645e-04  5.578008e-08 
##   C_ROS-C_ASC  S_ONLY-C_ASC   S_ASC-C_ASC   S_ROS-C_ASC  S_ONLY-C_ROS 
##  9.999638e-01  1.419284e-05  2.429165e-02  3.611648e-05  7.102623e-06 
##   S_ASC-C_ROS   S_ROS-C_ROS  S_ASC-S_ONLY  S_ROS-S_ONLY   S_ROS-S_ASC 
##  1.434529e-02  1.794707e-05  2.199215e-01  9.987846e-01  3.870697e-01
library(multcompView)
stat.test<- multcompLetters(Pval)
stat.test
##  C_ASC  C_ROS S_ONLY  S_ASC  S_ROS C_ONLY 
##    "a"    "a"    "b"    "b"    "b"    "a"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$genotype <- "M248"
test
# Now same for LA:
Output <- TukeyHSD(aov(LRno.delta ~ cond.final, data = GF_LA))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = LRno.delta ~ cond.final, data = GF_LA)
## 
## $cond.final
##                     diff         lwr        upr     p adj
## C_ASC-C_ONLY   1.1500000  0.38015512  1.9198449 0.0006616
## C_ROS-C_ONLY   0.7795455  0.02740052  1.5316904 0.0380372
## S_ONLY-C_ONLY -0.9250000 -1.69484488 -0.1551551 0.0099253
## S_ASC-C_ONLY  -0.2750000 -1.04484488  0.4948449 0.8968704
## S_ROS-C_ONLY  -1.2000000 -1.96984488 -0.4301551 0.0003470
## C_ROS-C_ASC   -0.3704545 -1.12259948  0.3816904 0.6940992
## S_ONLY-C_ASC  -2.0750000 -2.84484488 -1.3051551 0.0000000
## S_ASC-C_ASC   -1.4250000 -2.19484488 -0.6551551 0.0000166
## S_ROS-C_ASC   -2.3500000 -3.11984488 -1.5801551 0.0000000
## S_ONLY-C_ROS  -1.7045455 -2.45669039 -0.9524005 0.0000002
## S_ASC-C_ROS   -1.0545455 -1.80669039 -0.3024005 0.0016072
## S_ROS-C_ROS   -1.9795455 -2.73169039 -1.2274005 0.0000000
## S_ASC-S_ONLY   0.6500000 -0.11984488  1.4198449 0.1440740
## S_ROS-S_ONLY  -0.2750000 -1.04484488  0.4948449 0.8968704
## S_ROS-S_ASC   -0.9250000 -1.69484488 -0.1551551 0.0099253
Pval = Output$cond.final[,'p adj']
Pval
##  C_ASC-C_ONLY  C_ROS-C_ONLY S_ONLY-C_ONLY  S_ASC-C_ONLY  S_ROS-C_ONLY 
##  6.615807e-04  3.803720e-02  9.925329e-03  8.968704e-01  3.470202e-04 
##   C_ROS-C_ASC  S_ONLY-C_ASC   S_ASC-C_ASC   S_ROS-C_ASC  S_ONLY-C_ROS 
##  6.940992e-01  1.521435e-09  1.664460e-05  3.374545e-11  1.780586e-07 
##   S_ASC-C_ROS   S_ROS-C_ROS  S_ASC-S_ONLY  S_ROS-S_ONLY   S_ROS-S_ASC 
##  1.607239e-03  3.070617e-09  1.440740e-01  8.968704e-01  9.925329e-03
stat.test<- multcompLetters(Pval)
stat.test
##  C_ASC  C_ROS S_ONLY  S_ASC  S_ROS C_ONLY 
##    "a"    "a"   "bc"   "bd"    "c"    "d"
test2 <- as.data.frame(stat.test$Letters)
test2$group1 <- rownames(test)
test2$group2 <- rownames(test)
colnames(test2)[1] <- "Tukey"
test2
test2$genotype <- "LA"
test2
test3 <- rbind(test, test2)
test3
LRno.delta_p_geno <- ggerrorplot(growth_factors, y = "LRno.delta", x = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Lateral root number increase rate (# LR / day)") 
LRno.delta_p_geno <- LRno.delta_p_geno + rremove("legend") 
LRno.delta_p_geno <- LRno.delta_p_geno  +stat_pvalue_manual(test3, label = "Tukey", y.position = 4.5)  + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
LRno.delta_p_geno

Output <- TukeyHSD(aov(aLRL.delta ~ cond.final, data = GF_M248))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = aLRL.delta ~ cond.final, data = GF_M248)
## 
## $cond.final
##                       diff          lwr         upr     p adj
## C_ASC-C_ONLY  -0.032477349 -0.118379268 0.053424570 0.8725688
## C_ROS-C_ONLY  -0.062154650 -0.148056569 0.023747269 0.2844431
## S_ONLY-C_ONLY  0.029595554 -0.056306365 0.115497473 0.9102075
## S_ASC-C_ONLY  -0.041281202 -0.127183121 0.044620717 0.7155593
## S_ROS-C_ONLY  -0.047777648 -0.131704546 0.036149250 0.5500172
## C_ROS-C_ASC   -0.029677301 -0.115579220 0.056224618 0.9092453
## S_ONLY-C_ASC   0.062072903 -0.023829016 0.147974822 0.2858162
## S_ASC-C_ASC   -0.008803853 -0.094705771 0.077098066 0.9996407
## S_ROS-C_ASC   -0.015300299 -0.099227197 0.068626600 0.9943127
## S_ONLY-C_ROS   0.091750204  0.005848285 0.177652123 0.0297792
## S_ASC-C_ROS    0.020873448 -0.065028470 0.106775367 0.9790134
## S_ROS-C_ROS    0.014377002 -0.069549896 0.098303901 0.9957495
## S_ASC-S_ONLY  -0.070876755 -0.156778674 0.015025163 0.1619844
## S_ROS-S_ONLY  -0.077373201 -0.161300099 0.006553697 0.0870401
## S_ROS-S_ASC   -0.006496446 -0.090423344 0.077430452 0.9999096
Pval = Output$cond.final[,'p adj']
Pval
##  C_ASC-C_ONLY  C_ROS-C_ONLY S_ONLY-C_ONLY  S_ASC-C_ONLY  S_ROS-C_ONLY 
##    0.87256877    0.28444314    0.91020750    0.71555926    0.55001717 
##   C_ROS-C_ASC  S_ONLY-C_ASC   S_ASC-C_ASC   S_ROS-C_ASC  S_ONLY-C_ROS 
##    0.90924534    0.28581623    0.99964071    0.99431272    0.02977920 
##   S_ASC-C_ROS   S_ROS-C_ROS  S_ASC-S_ONLY  S_ROS-S_ONLY   S_ROS-S_ASC 
##    0.97901341    0.99574954    0.16198443    0.08704007    0.99990958
library(multcompView)
stat.test<- multcompLetters(Pval)
stat.test
##  C_ASC  C_ROS S_ONLY  S_ASC  S_ROS C_ONLY 
##   "ab"    "a"    "b"   "ab"   "ab"   "ab"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$genotype <- "M248"
test
# Now same for LA:
Output <- TukeyHSD(aov(aLRL.delta ~ cond.final, data = GF_LA))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = aLRL.delta ~ cond.final, data = GF_LA)
## 
## $cond.final
##                       diff         lwr           upr     p adj
## C_ASC-C_ONLY   0.005062847 -0.06738493  0.0775106281 0.9999455
## C_ROS-C_ONLY  -0.002259133 -0.07304123  0.0685229587 0.9999989
## S_ONLY-C_ONLY -0.095351881 -0.16779966 -0.0229040998 0.0035839
## S_ASC-C_ONLY  -0.048105522 -0.12055330  0.0243422594 0.3778352
## S_ROS-C_ONLY  -0.072029741 -0.14447752  0.0004180403 0.0521695
## C_ROS-C_ASC   -0.007321981 -0.07810407  0.0634601116 0.9996240
## S_ONLY-C_ASC  -0.100414728 -0.17286251 -0.0279669469 0.0018708
## S_ASC-C_ASC   -0.053168369 -0.12561615  0.0192794123 0.2698001
## S_ROS-C_ASC   -0.077092588 -0.14954037 -0.0046448068 0.0307181
## S_ONLY-C_ROS  -0.093092747 -0.16387484 -0.0223106553 0.0036150
## S_ASC-C_ROS   -0.045846388 -0.11662848  0.0249357040 0.4057833
## S_ROS-C_ROS   -0.069770607 -0.14055270  0.0010114849 0.0555224
## S_ASC-S_ONLY   0.047246359 -0.02520142  0.1196941403 0.3980905
## S_ROS-S_ONLY   0.023322140 -0.04912564  0.0957699212 0.9312493
## S_ROS-S_ASC   -0.023924219 -0.09637200  0.0485235619 0.9239201
Pval = Output$cond.final[,'p adj']
Pval
##  C_ASC-C_ONLY  C_ROS-C_ONLY S_ONLY-C_ONLY  S_ASC-C_ONLY  S_ROS-C_ONLY 
##   0.999945454   0.999998894   0.003583881   0.377835231   0.052169543 
##   C_ROS-C_ASC  S_ONLY-C_ASC   S_ASC-C_ASC   S_ROS-C_ASC  S_ONLY-C_ROS 
##   0.999624009   0.001870759   0.269800075   0.030718115   0.003614993 
##   S_ASC-C_ROS   S_ROS-C_ROS  S_ASC-S_ONLY  S_ROS-S_ONLY   S_ROS-S_ASC 
##   0.405783298   0.055522411   0.398090511   0.931249336   0.923920086
stat.test<- multcompLetters(Pval)
stat.test
##  C_ASC  C_ROS S_ONLY  S_ASC  S_ROS C_ONLY 
##    "a"   "ab"    "c"  "abc"   "bc"   "ab"
test2 <- as.data.frame(stat.test$Letters)
test2$group1 <- rownames(test)
test2$group2 <- rownames(test)
colnames(test2)[1] <- "Tukey"
test2
test2$genotype <- "LA"
test2
test3 <- rbind(test, test2)
test3
aLRL.delta_p_geno <- ggerrorplot(growth_factors, y = "aLRL.delta", x = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Average lateral root growth rate (cm / day)") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") 
aLRL.delta_p_geno <- aLRL.delta_p_geno  +stat_pvalue_manual(test3, label = "Tukey", y.position = 0.5)  + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
aLRL.delta_p_geno

#this is the growth rate with ttest

aLRL.delta_p_geno <- ggerrorplot(growth_factors, y = "aLRL.delta", x = "cond.final", fill="cond.final", color="cond.final", 
                        facet.by = c("genotype"), ncol=2,
                        desc_stat = "mean_sd", add = "jitter", 
                        add.params = list(color = "darkgray"),
                        xlab="", ylab="Average lateral root growth rate (cm / day)") 
aLRL.delta_p_geno <- aLRL.delta_p_geno + rremove("legend") + stat_compare_means(method="t.test", ref.group = "C_ONLY", 
                                                              label = "p.signif", hide.ns = T)
aLRL.delta_p_geno <- aLRL.delta_p_geno #+ ggtitle("average Lateral Root Growth")
aLRL.delta_p_geno

### growth factor calculations - Salt Stress Tolerance Indexes calculations. in theory, I have 3 STI. for c-only vs s-only, for c-asc vs s-asc, and for c-ros vs s-ros

library(doBy)
unique(growth_factors$cond.final)
## [1] C_ASC  S_ASC  C_ROS  S_ROS  C_ONLY S_ONLY
## Levels: C_ONLY C_ASC C_ROS S_ONLY S_ASC S_ROS

#need to check wheather I am doing it right # growth factor for Only

avg_growth <- summaryBy(data = growth_factors, MR.delta + aLRL.delta + LRno.delta ~ genotype + cond.final + ONLY)
avg_growth
control_treatments <- c("C_ONLY", "C_ASC", "C_ROS")
avg_growth5 <- subset(avg_growth, avg_growth$cond.final %in% control_treatments)
colnames(avg_growth5) <- gsub(".mean", ".control", colnames(avg_growth5))
avg_growth5$ONLY <- avg_growth5$cond.final
avg_growth5$ONLY <- gsub("C_", "", avg_growth5$ONLY)

avg_growth5 <- avg_growth5[,c(1,3:6)]
avg_growth5
growth_factors5 <- growth_factors
growth_factors5$ONLY <- growth_factors5$cond.final
growth_factors5$ONLY <- gsub("C_", "", growth_factors5$ONLY)
growth_factors5$ONLY <- gsub("S_", "", growth_factors5$ONLY)
unique(growth_factors5$ONLY)
## [1] "ASC"  "ROS"  "ONLY"
unique(growth_factors5$cond.final)
## [1] C_ASC  S_ASC  C_ROS  S_ROS  C_ONLY S_ONLY
## Levels: C_ONLY C_ASC C_ROS S_ONLY S_ASC S_ROS
growth_factors5
growth_factors6 <- merge(growth_factors5, avg_growth5, by = c("genotype", "ONLY"))
growth_factors6
growth_factors6$MR.STI <- growth_factors6$MR.delta / growth_factors6$MR.delta.control
growth_factors6$aLRL.STI <- growth_factors6$aLRL.delta / growth_factors6$aLRL.delta.control
growth_factors6$LRno.STI <- growth_factors6$LRno.delta / growth_factors6$LRno.delta.control
head(growth_factors6)
growth_factors6$ONLY <- factor(growth_factors6$ONLY, levels=c("ONLY", "ASC", "ROS"))
growth_factors6
#install.packages("ggbeeswarm")
library(ggbeeswarm)
## Warning: package 'ggbeeswarm' was built under R version 4.3.3
MR.STI_plot <- ggplot(data = growth_factors6, mapping = aes(x = ONLY, y = MR.STI, colour = ONLY)) 
MR.STI_plot <- MR.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
MR.STI_plot <- MR.STI_plot + facet_grid(~ genotype)
MR.STI_plot <- MR.STI_plot + stat_summary(fun=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("Salt tolerance index based on main root growth (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 = "ONLY", hide.ns = T) 
MR.STI_plot <- MR.STI_plot # + scale_color_manual(values=c("royalblue", "coral3"))
MR.STI_plot

aLRL.STI_plot <- ggplot(data = growth_factors6, mapping = aes(x = ONLY, y = aLRL.STI, colour = ONLY)) 
aLRL.STI_plot <- aLRL.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
aLRL.STI_plot <- aLRL.STI_plot + facet_grid(~ genotype)
aLRL.STI_plot <- aLRL.STI_plot + stat_summary(fun=mean, geom="point", shape=95, size=10, color="black", fill="black")
aLRL.STI_plot <- aLRL.STI_plot + theme(legend.position = "none")
aLRL.STI_plot <- aLRL.STI_plot + xlab("")
aLRL.STI_plot <- aLRL.STI_plot + ylab("Salt tolerance index based on average lateral root growth (Fraction of control)")
aLRL.STI_plot <- aLRL.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "ONLY", hide.ns = T) 
aLRL.STI_plot <- aLRL.STI_plot #+ scale_color_manual(values=c("royalblue", "coral3"))
aLRL.STI_plot

LRno.STI_plot <- ggplot(data = growth_factors6, mapping = aes(x = ONLY, y = LRno.STI, colour = ONLY)) 
LRno.STI_plot <- LRno.STI_plot + geom_beeswarm(alpha=0.6, priority = "density")
LRno.STI_plot <- LRno.STI_plot + facet_grid(~ genotype)
LRno.STI_plot <- LRno.STI_plot + stat_summary(fun=mean, geom="point", shape=95, size=10, color="black", fill="black")
LRno.STI_plot <- LRno.STI_plot + theme(legend.position = "none")
LRno.STI_plot <- LRno.STI_plot + xlab("")
LRno.STI_plot <- LRno.STI_plot + ylab("Salt tolerance index based on average increase in lateral root number (Fraction of control)")
LRno.STI_plot <- LRno.STI_plot + stat_compare_means(label = "p.signif", method = "t.test", ref.group = "ONLY", hide.ns = T) 
LRno.STI_plot <- LRno.STI_plot #+ scale_color_manual(values=c("royalblue", "coral3"))
LRno.STI_plot

Now - lets have a look what happens if we compare STIs of individual RSA components:

growth_factors6
GF6 <- growth_factors6[,c(1:4, 11:13)]
control <- c("C_ROS", "C_ASC", "C_ONLY")
GF6 <- subset(GF6, !(GF6$cond.final %in% control))
GF6
library(reshape2)
mGF6 <- melt(GF6, id=c("genotype", "ONLY", "root_name", "cond.final"))
mGF6
STI_plot <- ggplot(data = mGF6, mapping = aes(x = variable, y = value, colour = ONLY, group = root_name)) 
STI_plot <- STI_plot + geom_line(alpha=0.6)
STI_plot <- STI_plot + facet_grid(ONLY ~ genotype)
STI_plot <- STI_plot# + stat_summary(fun=mean, geom="point", shape=95, size=10, color="black", fill="black")
STI_plot <- STI_plot + theme(legend.position = "none")
STI_plot <- STI_plot + xlab("")
STI_plot <- STI_plot + ylab("Salt tolerance index (Fraction of control)")
STI_plot

STI_plot2 <- ggplot(data = mGF6, mapping = aes(x = variable, y = value, colour = ONLY)) 
STI_plot2 <- STI_plot2 + geom_beeswarm(alpha=0.6, priority = "density")
STI_plot2 <- STI_plot2 + facet_grid(ONLY ~ genotype)
STI_plot2 <- STI_plot2 + stat_summary(fun=mean, geom="point", shape=95, size=10, color="black", fill="black")
STI_plot2 <- STI_plot2 + theme(legend.position = "none")
STI_plot2 <- STI_plot2 + xlab("")
STI_plot2 <- STI_plot2 + ylab("Salt tolerance index (Fraction of control)")
STI_plot2

now - we would like to compare between treatments - how does ROS/ASC compare to ONLY treatment. For that - we need to merge ONLY with trait:

mGF6$id <- paste(mGF6$ONLY, mGF6$variable, sep="_")
mGF6
mGF6$id <- factor(mGF6$id, levels=c("ONLY_MR.STI", "ONLY_aLRL.STI", "ONLY_LRno.STI", "ASC_MR.STI", "ASC_aLRL.STI", "ASC_LRno.STI", "ROS_MR.STI", "ROS_aLRL.STI", "ROS_LRno.STI"))
mGF6
STI_plot3 <- ggplot(data = mGF6, mapping = aes(x = id, y = value, colour = ONLY)) 
STI_plot3 <- STI_plot3 + geom_beeswarm(alpha=0.6, priority = "density")
STI_plot3 <- STI_plot3 + facet_grid(~ genotype)
STI_plot3 <- STI_plot3 + stat_summary(fun=mean, geom="point", shape=95, size=10, color="black", fill="black")
STI_plot3 <- STI_plot3 + theme(legend.position = "none")
STI_plot3 <- STI_plot3 + xlab("")
STI_plot3 <- STI_plot3 + ylab("Salt tolerance index (Fraction of control)") + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))
STI_plot3

now - let’s add Tukey leters to it:

# M248
mGF6_M248 <- subset(mGF6, mGF6$genotype == "M248")
mGF6_LA <- subset(mGF6, mGF6$genotype == "LA")


Output <- TukeyHSD(aov(value ~ id, data = mGF6_M248))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = value ~ id, data = mGF6_M248)
## 
## $id
##                                     diff          lwr         upr     p adj
## ONLY_aLRL.STI-ONLY_MR.STI    0.174007549 -0.229531124  0.57754622 0.9047560
## ONLY_LRno.STI-ONLY_MR.STI   -0.535812960 -0.939351633 -0.13227429 0.0018766
## ASC_MR.STI-ONLY_MR.STI      -0.205983044 -0.609521717  0.19755563 0.7882492
## ASC_aLRL.STI-ONLY_MR.STI    -0.006814419 -0.410353092  0.39672425 1.0000000
## ASC_LRno.STI-ONLY_MR.STI    -0.258300998 -0.661839671  0.14523768 0.5219149
## ROS_MR.STI-ONLY_MR.STI      -0.107700873 -0.501961555  0.28655981 0.9939604
## ROS_aLRL.STI-ONLY_MR.STI     0.130591771 -0.263668911  0.52485245 0.9789165
## ROS_LRno.STI-ONLY_MR.STI    -0.442814005 -0.837074686 -0.04855332 0.0161294
## ONLY_LRno.STI-ONLY_aLRL.STI -0.709820509 -1.113359182 -0.30628184 0.0000092
## ASC_MR.STI-ONLY_aLRL.STI    -0.379990593 -0.783529266  0.02354808 0.0812164
## ASC_aLRL.STI-ONLY_aLRL.STI  -0.180821968 -0.584360641  0.22271671 0.8842749
## ASC_LRno.STI-ONLY_aLRL.STI  -0.432308547 -0.835847220 -0.02876987 0.0264117
## ROS_MR.STI-ONLY_aLRL.STI    -0.281708422 -0.675969104  0.11255226 0.3688612
## ROS_aLRL.STI-ONLY_aLRL.STI  -0.043415778 -0.437676460  0.35084490 0.9999927
## ROS_LRno.STI-ONLY_aLRL.STI  -0.616821553 -1.011082235 -0.22256087 0.0001124
## ASC_MR.STI-ONLY_LRno.STI     0.329829915 -0.073708758  0.73336859 0.2013311
## ASC_aLRL.STI-ONLY_LRno.STI   0.528998540  0.125459867  0.93253721 0.0022661
## ASC_LRno.STI-ONLY_LRno.STI   0.277511962 -0.126026711  0.68105063 0.4218530
## ROS_MR.STI-ONLY_LRno.STI     0.428112087  0.033851405  0.82237277 0.0230495
## ROS_aLRL.STI-ONLY_LRno.STI   0.666404730  0.272144049  1.06066541 0.0000226
## ROS_LRno.STI-ONLY_LRno.STI   0.092998955 -0.301261727  0.48725964 0.9978057
## ASC_aLRL.STI-ASC_MR.STI      0.199168625 -0.204370048  0.60270730 0.8172429
## ASC_LRno.STI-ASC_MR.STI     -0.052317954 -0.455856627  0.35122072 0.9999743
## ROS_MR.STI-ASC_MR.STI        0.098282171 -0.295978510  0.49254285 0.9967709
## ROS_aLRL.STI-ASC_MR.STI      0.336574815 -0.057685867  0.73083550 0.1573301
## ROS_LRno.STI-ASC_MR.STI     -0.236830960 -0.631091642  0.15742972 0.6070785
## ASC_LRno.STI-ASC_aLRL.STI   -0.251486579 -0.655025252  0.15205209 0.5584057
## ROS_MR.STI-ASC_aLRL.STI     -0.100886454 -0.495147135  0.29337423 0.9961311
## ROS_aLRL.STI-ASC_aLRL.STI    0.137406190 -0.256854492  0.53166687 0.9712150
## ROS_LRno.STI-ASC_aLRL.STI   -0.435999585 -0.830260267 -0.04173890 0.0190601
## ROS_MR.STI-ASC_LRno.STI      0.150600125 -0.243660557  0.54486081 0.9507075
## ROS_aLRL.STI-ASC_LRno.STI    0.388892769 -0.005367913  0.78315345 0.0561667
## ROS_LRno.STI-ASC_LRno.STI   -0.184513007 -0.578773688  0.20974768 0.8570087
## ROS_aLRL.STI-ROS_MR.STI      0.238292644 -0.146466385  0.62305167 0.5667957
## ROS_LRno.STI-ROS_MR.STI     -0.335113132 -0.719872160  0.04964590 0.1391359
## ROS_LRno.STI-ROS_aLRL.STI   -0.573405775 -0.958164803 -0.18864675 0.0002804
Pval = Output$id[,'p adj']
Pval
##   ONLY_aLRL.STI-ONLY_MR.STI   ONLY_LRno.STI-ONLY_MR.STI 
##                9.047560e-01                1.876576e-03 
##      ASC_MR.STI-ONLY_MR.STI    ASC_aLRL.STI-ONLY_MR.STI 
##                7.882492e-01                1.000000e+00 
##    ASC_LRno.STI-ONLY_MR.STI      ROS_MR.STI-ONLY_MR.STI 
##                5.219149e-01                9.939604e-01 
##    ROS_aLRL.STI-ONLY_MR.STI    ROS_LRno.STI-ONLY_MR.STI 
##                9.789165e-01                1.612938e-02 
## ONLY_LRno.STI-ONLY_aLRL.STI    ASC_MR.STI-ONLY_aLRL.STI 
##                9.154416e-06                8.121635e-02 
##  ASC_aLRL.STI-ONLY_aLRL.STI  ASC_LRno.STI-ONLY_aLRL.STI 
##                8.842749e-01                2.641174e-02 
##    ROS_MR.STI-ONLY_aLRL.STI  ROS_aLRL.STI-ONLY_aLRL.STI 
##                3.688612e-01                9.999927e-01 
##  ROS_LRno.STI-ONLY_aLRL.STI    ASC_MR.STI-ONLY_LRno.STI 
##                1.123996e-04                2.013311e-01 
##  ASC_aLRL.STI-ONLY_LRno.STI  ASC_LRno.STI-ONLY_LRno.STI 
##                2.266075e-03                4.218530e-01 
##    ROS_MR.STI-ONLY_LRno.STI  ROS_aLRL.STI-ONLY_LRno.STI 
##                2.304950e-02                2.258298e-05 
##  ROS_LRno.STI-ONLY_LRno.STI     ASC_aLRL.STI-ASC_MR.STI 
##                9.978057e-01                8.172429e-01 
##     ASC_LRno.STI-ASC_MR.STI       ROS_MR.STI-ASC_MR.STI 
##                9.999743e-01                9.967709e-01 
##     ROS_aLRL.STI-ASC_MR.STI     ROS_LRno.STI-ASC_MR.STI 
##                1.573301e-01                6.070785e-01 
##   ASC_LRno.STI-ASC_aLRL.STI     ROS_MR.STI-ASC_aLRL.STI 
##                5.584057e-01                9.961311e-01 
##   ROS_aLRL.STI-ASC_aLRL.STI   ROS_LRno.STI-ASC_aLRL.STI 
##                9.712150e-01                1.906010e-02 
##     ROS_MR.STI-ASC_LRno.STI   ROS_aLRL.STI-ASC_LRno.STI 
##                9.507075e-01                5.616669e-02 
##   ROS_LRno.STI-ASC_LRno.STI     ROS_aLRL.STI-ROS_MR.STI 
##                8.570087e-01                5.667957e-01 
##     ROS_LRno.STI-ROS_MR.STI   ROS_LRno.STI-ROS_aLRL.STI 
##                1.391359e-01                2.803870e-04
stat.test<- multcompLetters(Pval)
stat.test
## ONLY_aLRL.STI ONLY_LRno.STI    ASC_MR.STI  ASC_aLRL.STI  ASC_LRno.STI 
##           "a"           "b"        "abcd"          "ac"         "bcd" 
##    ROS_MR.STI  ROS_aLRL.STI  ROS_LRno.STI   ONLY_MR.STI 
##         "acd"          "ac"          "bd"          "ac"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
colnames(test)[1] <- "Tukey"
test
test$genotype <- "M248"
test
# Now same for LA:
Output <- TukeyHSD(aov(value ~ id, data = mGF6_LA))
Output
##   Tukey multiple comparisons of means
##     95% family-wise confidence level
## 
## Fit: aov(formula = value ~ id, data = mGF6_LA)
## 
## $id
##                                    diff          lwr          upr     p adj
## ONLY_aLRL.STI-ONLY_MR.STI   -0.34318173 -0.684039953 -0.002323515 0.0471461
## ONLY_LRno.STI-ONLY_MR.STI   -0.26642746 -0.607285682  0.074430756 0.2517883
## ASC_MR.STI-ONLY_MR.STI      -0.09423934 -0.435097562  0.246618876 0.9933884
## ASC_aLRL.STI-ONLY_MR.STI    -0.05358145 -0.394439666  0.287276772 0.9998869
## ASC_LRno.STI-ONLY_MR.STI    -0.21861344 -0.559471655  0.122244783 0.5180389
## ROS_MR.STI-ONLY_MR.STI      -0.01111592 -0.351974139  0.329742299 1.0000000
## ROS_aLRL.STI-ONLY_MR.STI    -0.18139623 -0.522254453  0.159461985 0.7470630
## ROS_LRno.STI-ONLY_MR.STI    -0.52067014 -0.861528357 -0.179811919 0.0001828
## ONLY_LRno.STI-ONLY_aLRL.STI  0.07675427 -0.264103948  0.417612490 0.9984040
## ASC_MR.STI-ONLY_aLRL.STI     0.24894239 -0.091915828  0.589800610 0.3383353
## ASC_aLRL.STI-ONLY_aLRL.STI   0.28960029 -0.051257932  0.630458506 0.1615114
## ASC_LRno.STI-ONLY_aLRL.STI   0.12456830 -0.216289921  0.465426517 0.9615277
## ROS_MR.STI-ONLY_aLRL.STI     0.33206581 -0.008792405  0.672924033 0.0621938
## ROS_aLRL.STI-ONLY_aLRL.STI   0.16178550 -0.179072719  0.502643719 0.8463534
## ROS_LRno.STI-ONLY_aLRL.STI  -0.17748840 -0.518346623  0.163369815 0.7686390
## ASC_MR.STI-ONLY_LRno.STI     0.17218812 -0.168670100  0.513046338 0.7965768
## ASC_aLRL.STI-ONLY_LRno.STI   0.21284602 -0.128012204  0.553704234 0.5545617
## ASC_LRno.STI-ONLY_LRno.STI   0.04781403 -0.293044192  0.388672246 0.9999524
## ROS_MR.STI-ONLY_LRno.STI     0.25531154 -0.085546677  0.596169762 0.3050587
## ROS_aLRL.STI-ONLY_LRno.STI   0.08503123 -0.255826990  0.425889448 0.9967216
## ROS_LRno.STI-ONLY_LRno.STI  -0.25424268 -0.595100895  0.086615544 0.3105082
## ASC_aLRL.STI-ASC_MR.STI      0.04065790 -0.300200323  0.381516115 0.9999864
## ASC_LRno.STI-ASC_MR.STI     -0.12437409 -0.465232311  0.216484127 0.9618777
## ROS_MR.STI-ASC_MR.STI        0.08312342 -0.257734796  0.423981642 0.9972002
## ROS_aLRL.STI-ASC_MR.STI     -0.08715689 -0.428015110  0.253701328 0.9961120
## ROS_LRno.STI-ASC_MR.STI     -0.42643079 -0.767289014 -0.085572576 0.0043940
## ASC_LRno.STI-ASC_aLRL.STI   -0.16503199 -0.505890207  0.175826231 0.8315858
## ROS_MR.STI-ASC_aLRL.STI      0.04246553 -0.298392692  0.383323746 0.9999809
## ROS_aLRL.STI-ASC_aLRL.STI   -0.12781479 -0.468673006  0.213043432 0.9553297
## ROS_LRno.STI-ASC_aLRL.STI   -0.46708869 -0.807946910 -0.126230472 0.0011760
## ROS_MR.STI-ASC_LRno.STI      0.20749752 -0.133360704  0.548355734 0.5885226
## ROS_aLRL.STI-ASC_LRno.STI    0.03721720 -0.303641018  0.378075420 0.9999931
## ROS_LRno.STI-ASC_LRno.STI   -0.30205670 -0.642914922  0.038801516 0.1242955
## ROS_aLRL.STI-ROS_MR.STI     -0.17028031 -0.511138533  0.170577905 0.8062276
## ROS_LRno.STI-ROS_MR.STI     -0.50955422 -0.850412437 -0.168695999 0.0002718
## ROS_LRno.STI-ROS_aLRL.STI   -0.33927390 -0.680132123  0.001584315 0.0520308
Pval = Output$id[,'p adj']
Pval
##   ONLY_aLRL.STI-ONLY_MR.STI   ONLY_LRno.STI-ONLY_MR.STI 
##                0.0471461421                0.2517882981 
##      ASC_MR.STI-ONLY_MR.STI    ASC_aLRL.STI-ONLY_MR.STI 
##                0.9933884122                0.9998869112 
##    ASC_LRno.STI-ONLY_MR.STI      ROS_MR.STI-ONLY_MR.STI 
##                0.5180389044                0.9999999995 
##    ROS_aLRL.STI-ONLY_MR.STI    ROS_LRno.STI-ONLY_MR.STI 
##                0.7470630234                0.0001827977 
## ONLY_LRno.STI-ONLY_aLRL.STI    ASC_MR.STI-ONLY_aLRL.STI 
##                0.9984040220                0.3383352862 
##  ASC_aLRL.STI-ONLY_aLRL.STI  ASC_LRno.STI-ONLY_aLRL.STI 
##                0.1615113739                0.9615276647 
##    ROS_MR.STI-ONLY_aLRL.STI  ROS_aLRL.STI-ONLY_aLRL.STI 
##                0.0621938231                0.8463534042 
##  ROS_LRno.STI-ONLY_aLRL.STI    ASC_MR.STI-ONLY_LRno.STI 
##                0.7686390143                0.7965767866 
##  ASC_aLRL.STI-ONLY_LRno.STI  ASC_LRno.STI-ONLY_LRno.STI 
##                0.5545616758                0.9999524429 
##    ROS_MR.STI-ONLY_LRno.STI  ROS_aLRL.STI-ONLY_LRno.STI 
##                0.3050587379                0.9967216015 
##  ROS_LRno.STI-ONLY_LRno.STI     ASC_aLRL.STI-ASC_MR.STI 
##                0.3105081744                0.9999863506 
##     ASC_LRno.STI-ASC_MR.STI       ROS_MR.STI-ASC_MR.STI 
##                0.9618776790                0.9972002320 
##     ROS_aLRL.STI-ASC_MR.STI     ROS_LRno.STI-ASC_MR.STI 
##                0.9961119834                0.0043939971 
##   ASC_LRno.STI-ASC_aLRL.STI     ROS_MR.STI-ASC_aLRL.STI 
##                0.8315858094                0.9999808923 
##   ROS_aLRL.STI-ASC_aLRL.STI   ROS_LRno.STI-ASC_aLRL.STI 
##                0.9553296753                0.0011759969 
##     ROS_MR.STI-ASC_LRno.STI   ROS_aLRL.STI-ASC_LRno.STI 
##                0.5885225841                0.9999931309 
##   ROS_LRno.STI-ASC_LRno.STI     ROS_aLRL.STI-ROS_MR.STI 
##                0.1242954745                0.8062276180 
##     ROS_LRno.STI-ROS_MR.STI   ROS_LRno.STI-ROS_aLRL.STI 
##                0.0002717524                0.0520307806
stat.test<- multcompLetters(Pval)
stat.test
## ONLY_aLRL.STI ONLY_LRno.STI    ASC_MR.STI  ASC_aLRL.STI  ASC_LRno.STI 
##          "ab"         "abc"          "ac"          "ac"         "abc" 
##    ROS_MR.STI  ROS_aLRL.STI  ROS_LRno.STI   ONLY_MR.STI 
##          "ac"         "abc"           "b"           "c"
test2 <- as.data.frame(stat.test$Letters)
test2$group1 <- rownames(test)
test2$group2 <- rownames(test)
colnames(test2)[1] <- "Tukey"
test2
test2$genotype <- "LA"
test2
test3 <- rbind(test, test2)
test3
STI_plot3 <- STI_plot3 + stat_pvalue_manual(test3, label = "Tukey", y.position = 2) 
STI_plot3