This R-notebook describes the 1st experiment where the metabolites extracted from cowpea leaves (14 days after drought / mock treatment application) were used to enrich agar plates with 4 days old Arabidopsis seedlings (Col-0) exposed to osmotic stress (-0.25 MPa >> 250 g PEG8000 / 1 L of liquid 1/2 MS used for 24h infusions). The control plates (noPEG) were not infused with a liquid MS, but rather kept at 4C for o/n during the time of infusion of PEG plates.

The treatments are as follow: - mock - an empty ependorf tube that underwent speed-vac and sonicator bath treatments as negative control - E02 - an empty metabolite extraction - 100 uL of original extraction dissolved in 1 ml of MQ water after speed-vac - CP01 - a metabolite extraction from MMJ22-06CP01 sample (UCR779 accession grown in Control conditions) - 100 uL of original extraction dissolved in 1 ml of MQ water after speed-vac - CP06 - a metabolite extraction from MMJ22-06CP06 sample (UCR779 accession grown in Drought conditions - 10% soil water-holding capacity) - 100 uL of original extraction dissolved in 1 ml of MQ water after speed-vac

Unfortunately - due to light conditions during the day time - the automated root tracing did not work - and thus we will only analyze the last timepoint for this experiment (for now).

Load and prepare the data

list.files(pattern="20230327")
## [1] "20230327_MMJCP1_plate1.csv" "20230327_MMJCP1_plate2.csv"
## [3] "20230327_MMJCP1_plate3.csv" "20230327_MMJCP1_plate4.csv"
## [5] "20230327_MMJCP2_plate1.csv" "20230327_MMJCP2_plate2.csv"
## [7] "20230327_MMJCP2_plate3.csv" "20230327_MMJCP2_plate4.csv"
data1 <- read.csv("20230327_MMJCP1_plate1.csv")
data2 <- read.csv("20230327_MMJCP1_plate2.csv")
data3 <- read.csv("20230327_MMJCP1_plate3.csv")
data4 <- read.csv("20230327_MMJCP1_plate4.csv")
data5 <- read.csv("20230327_MMJCP2_plate1.csv")
data6 <- read.csv("20230327_MMJCP2_plate2.csv")
data7 <- read.csv("20230327_MMJCP2_plate3.csv")
data8 <- read.csv("20230327_MMJCP2_plate4.csv")
head(data1)

Let’s add the individual treatments for each plate:

data1$treatment1 <- "E"
data1$treatment2 <- "noPEG"
data2$treatment1 <- "CP01"
data2$treatment2 <- "noPEG"
data3$treatment1 <- "CP01"
data3$treatment2 <- "PEG"
data4$treatment1 <- "E"
data4$treatment2 <- "PEG"
data5$treatment1 <- "mock"
data5$treatment2 <- "PEG"
data6$treatment1 <- "mock"
data6$treatment2 <- "noPEG"
data7$treatment1 <- "CP06"
data7$treatment2 <- "noPEG"
data8$treatment1 <- "CP06"
data8$treatment2 <- "PEG"

data <- rbind(data1, data2)
data <- rbind(data, data3)
data <- rbind(data, data4)
data <- rbind(data, data5)
data <- rbind(data, data6)
data <- rbind(data, data7)
data <- rbind(data, data8)
data

let’s clean up the data to contain only interesting collumns:

colnames(data)
##  [1] "image"                 "root_name"             "root"                 
##  [4] "length"                "surface"               "volume"               
##  [7] "diameter"              "root_order"            "root_ontology"        
## [10] "parent_name"           "parent"                "insertion_position"   
## [13] "insertion_angle"       "n_child"               "child_density"        
## [16] "first_child"           "insertion_first_child" "last_child"           
## [19] "insertion_last_child"  "treatment1"            "treatment2"
data.new <- data[,c(1:4, 20:21)]
data.new

plot data

library(ggplot2)
library(ggpubr)
library(ggbeeswarm)
library(ggsci)

MRL_graph <- ggplot(data = data.new, mapping = aes(x = treatment1, y = length, color = treatment1)) 
MRL_graph <- MRL_graph + geom_beeswarm(alpha=0.6, priority = "density")
MRL_graph <- MRL_graph + facet_wrap(~ treatment2, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
MRL_graph <- MRL_graph + stat_compare_means(aes(group = treatment1), method = "aov")
MRL_graph <- MRL_graph + theme_bw() + theme(legend.position = "none") + ylab("root length (cm)") + xlab("") + color_palette("aaas")
MRL_graph

Now - let’s calculate Tukey HSD:

library(stringr)
library(multcompView)

MapC <- subset(data.new, data.new$treatment2 == "noPEG")
MapD <- subset(data.new, data.new$treatment2 == "PEG")

Output <- TukeyHSD(aov(length ~ treatment1, data = MapC))
P7 = Output$treatment1[,'p adj']
stat.test<- multcompLetters(P7)
stat.test
## CP06    E mock CP01 
##  "a" "ab"  "c"  "b"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
test$treatment2 <- "noPEG"
colnames(test)[1] <- "Tukey"
testC <- test 

Output <- TukeyHSD(aov(length ~ treatment1, data = MapD))
P7 = Output$treatment1[,'p adj']
stat.test<- multcompLetters(P7)
stat.test
## CP06    E mock CP01 
##  "a" "ab"  "b"  "b"
test <- as.data.frame(stat.test$Letters)
test$group1 <- rownames(test)
test$group2 <- rownames(test)
test$treatment2 <- "PEG"
colnames(test)[1] <- "Tukey"
testD <- test

test <- rbind(testC, testD)
test
MRL_graph <- ggplot(data = data.new, mapping = aes(x = treatment1, y = length, color = treatment1)) 
MRL_graph <- MRL_graph + geom_beeswarm(alpha=0.6, priority = "density")
MRL_graph <- MRL_graph + facet_wrap(~ treatment2, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
MRL_graph <- MRL_graph + stat_compare_means(aes(group = treatment1), method = "aov")
MRL_graph <- MRL_graph + theme_bw() + theme(legend.position = "none") + ylab("root length (cm)") + xlab("") + color_palette("aaas")
MRL_graph <- MRL_graph + stat_pvalue_manual(test, label = "Tukey", y.position = 3)
MRL_graph