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