This experiment is describing the analysis of the root cross section. The roots were collected in the lab of Prof. Avat Shekoofa (Uni. of Tennessee), and the root samples were sent to Julkowska lab (BTI). Magda did hand cross sections and stained them with toluidine blue to visualize the xylem vessels. The individual images were stitched into an ortho-mosaic using Photoshop. The images were subsequently quantified in ImageJ for the area of the a) root, b) secondary xylem, c) xylem vessel members and the d) number of xylem vessel members (as characterized here)

All of the image data is available at Zenodo repository

#Data import and wrestling

list.files(pattern="csv")
## [1] "Avat_cotton_collaborationn - Coding.csv"      
## [2] "Avat_cotton_collaborationn - ImageJ Clean.csv"
data <- read.csv("Avat_cotton_collaborationn - ImageJ Clean.csv")
coding <- read.csv("Avat_cotton_collaborationn - Coding.csv")

data
coding

OK - let’s keep only the relevant tissues that were measured in all samples:

unique(data$Trait)
## [1] "Root Thickness" "Xylem"          "ind.xyl"        "cum.xyl"       
## [5] "vessel.no"
keep <- c("Root Thickness", "Xylem", "cum.xyl", "vessel.no" )
data1 <- subset(data, data$Trait %in% keep)
data1

cool - now I have to make sure to change my image name into actual sample number:

data1$image <- gsub(".png", "", data1$image)
data1$image <- gsub("Sample", "", data1$image)
data1$image <- gsub("a", "", data1$image)
data1$image <- gsub("b", "", data1$image)
data1$image <- gsub("c", "", data1$image)
unique(data1$image)
##  [1] "01" "02" "03" "04" "05" "06" "07" "08" "09" "10" "11" "12" "13" "14" "15"
## [16] "16" "17" "18" "19" "20" "21" "22" "23" "24" "25" "26" "27" "28" "29" "30"
## [31] "31" "32"
data2 <- data1
data2$image <- as.numeric(as.character(data2$image))
data2
library(reshape2)
data3 <- dcast(data2, image ~ Trait)
## Using Value as value column: use value.var to override.
data3
colnames(coding)[1] <- "image"
data4 <- merge(coding, data3, by = "image", all = TRUE)
data4

OK - now we got almost everything - except that I d like to calculate some ratios too

colnames(data4)[7] <- "Root.area"
colnames(data4)[9] <- "Xylem.area"
colnames(data4)[6] <- "Xylem.vessel.area"
data4$Xyl.p.root <- data4$Xylem.area / data4$Root.area
data4$XylVessel.p.root <- data4$Xylem.vessel.area / data4$Root.area
data4$XylVessel.p.Xylem <- data4$Xylem.vessel.area / data4$Xylem.area
data4$XylVessel.mean.area <- data4$Xylem.vessel.area / data4$vessel.no
data4

Data plotting:

OK - now we got everything organized - let’s start plotting the data!

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

data4$VPD..kPa. <- factor(data4$VPD..kPa., levels = c("2.1", "1.1"))

XylVA_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = Xylem.vessel.area, color = Cultivar)) 
XylVA_graph <- XylVA_graph + geom_beeswarm(alpha=0.6, priority = "density")
XylVA_graph <- XylVA_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
## Warning: The `fun.y` argument of `stat_summary()` is deprecated as of ggplot2 3.3.0.
## ℹ Please use the `fun` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
XylVA_graph <- XylVA_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
XylVA_graph <- XylVA_graph + theme_bw() + theme(legend.position = "none") + ylab("Xylem vessel area (pix)") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
XylVA_graph

XylVA_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = Xylem.vessel.area, color = VPD..kPa.)) 
XylVA_graph2 <- XylVA_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
XylVA_graph2 <- XylVA_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
XylVA_graph2 <- XylVA_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
XylVA_graph2 <- XylVA_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Xylem vessel area (pix)") + xlab("") + color_palette("aaas")
XylVA_graph2

data4
Root_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = Root.area, color = Cultivar)) 
Root_graph <- Root_graph + geom_beeswarm(alpha=0.6, priority = "density")
Root_graph <- Root_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Root_graph <- Root_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
Root_graph <- Root_graph + theme_bw() + theme(legend.position = "none") + ylab("Root area (pix)") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
Root_graph

Root_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = Xylem.vessel.area, color = VPD..kPa.)) 
Root_graph2 <- Root_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
Root_graph2 <- Root_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Root_graph2 <- Root_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
Root_graph2 <- Root_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Root area (pix)") + xlab("") + color_palette("aaas")
Root_graph2

Vessel.no_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = vessel.no, color = Cultivar)) 
Vessel.no_graph <- Vessel.no_graph + geom_beeswarm(alpha=0.6, priority = "density")
Vessel.no_graph <- Vessel.no_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Vessel.no_graph <- Vessel.no_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
Vessel.no_graph <- Vessel.no_graph + theme_bw() + theme(legend.position = "none") + ylab("Xylem vessel number") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
Vessel.no_graph

Vessel.no_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = vessel.no, color = VPD..kPa.)) 
Vessel.no_graph2 <- Vessel.no_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
Vessel.no_graph2 <- Vessel.no_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Vessel.no_graph2 <- Vessel.no_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
Vessel.no_graph2 <- Vessel.no_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Xylem vessel number") + xlab("") + color_palette("aaas")
Vessel.no_graph2

Xyl.a_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = Xylem.area, color = Cultivar)) 
Xyl.a_graph <- Xyl.a_graph + geom_beeswarm(alpha=0.6, priority = "density")
Xyl.a_graph <- Xyl.a_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Xyl.a_graph <- Xyl.a_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
Xyl.a_graph <- Xyl.a_graph + theme_bw() + theme(legend.position = "none") + ylab("Xylem area (pix)") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
Xyl.a_graph

Xyl.a_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = vessel.no, color = VPD..kPa.)) 
Xyl.a_graph2 <- Xyl.a_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
Xyl.a_graph2 <- Xyl.a_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Xyl.a_graph2 <- Xyl.a_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
Xyl.a_graph2 <- Xyl.a_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Xylem area (pix)") + xlab("") + color_palette("aaas")
Xyl.a_graph2

Xyl.p.root_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = Xyl.p.root, color = Cultivar)) 
Xyl.p.root_graph <- Xyl.p.root_graph + geom_beeswarm(alpha=0.6, priority = "density")
Xyl.p.root_graph <- Xyl.p.root_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Xyl.p.root_graph <- Xyl.p.root_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
Xyl.p.root_graph <- Xyl.p.root_graph + theme_bw() + theme(legend.position = "none") + ylab("Xylem per root area (ratio)") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
Xyl.p.root_graph

Xyl.p.root_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = Xyl.p.root, color = VPD..kPa.)) 
Xyl.p.root_graph2 <- Xyl.p.root_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
Xyl.p.root_graph2 <- Xyl.p.root_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
Xyl.p.root_graph2 <- Xyl.p.root_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
Xyl.p.root_graph2 <- Xyl.p.root_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Xylem per root area (ratio)") + xlab("") + color_palette("aaas")
Xyl.p.root_graph2

XylV.p.Xyl_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = XylVessel.p.Xylem, color = Cultivar)) 
XylV.p.Xyl_graph <- XylV.p.Xyl_graph + geom_beeswarm(alpha=0.6, priority = "density")
XylV.p.Xyl_graph <- XylV.p.Xyl_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
XylV.p.Xyl_graph <- XylV.p.Xyl_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
XylV.p.Xyl_graph <- XylV.p.Xyl_graph + theme_bw() + theme(legend.position = "none") + ylab("Xylem Vessel per Xylem area (ratio)") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
XylV.p.Xyl_graph

XylV.p.Xyl_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = XylVessel.p.Xylem, color = VPD..kPa.)) 
XylV.p.Xyl_graph2 <- XylV.p.Xyl_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
XylV.p.Xyl_graph2 <- XylV.p.Xyl_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
XylV.p.Xyl_graph2 <- XylV.p.Xyl_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
XylV.p.Xyl_graph2 <- XylV.p.Xyl_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Xylem Vessel per Xylem area (ratio)") + xlab("") + color_palette("aaas")
XylV.p.Xyl_graph2

XylV.meanArea_graph <- ggplot(data = data4, mapping = aes(x = Cultivar, y = XylVessel.mean.area, color = Cultivar)) 
XylV.meanArea_graph <- XylV.meanArea_graph + geom_beeswarm(alpha=0.6, priority = "density")
XylV.meanArea_graph <- XylV.meanArea_graph + facet_wrap(~ VPD..kPa., ncol = 2) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
XylV.meanArea_graph <- XylV.meanArea_graph + stat_compare_means(aes(group = Cultivar), method = "aov")
XylV.meanArea_graph <- XylV.meanArea_graph + theme_bw() + theme(legend.position = "none") + ylab("Xylem Vessel mean area (pix)") + xlab("") + color_palette("jco") + theme(axis.text.x = element_text(angle = 90))
XylV.meanArea_graph

XylV.meanArea_graph2 <- ggplot(data = data4, mapping = aes(x = VPD..kPa., y = XylVessel.mean.area, color = VPD..kPa.)) 
XylV.meanArea_graph2 <- XylV.meanArea_graph2 + geom_beeswarm(alpha=0.6, priority = "density")
XylV.meanArea_graph2 <- XylV.meanArea_graph2 + facet_wrap(~ Cultivar, ncol = 4) + stat_summary(fun.y=mean, geom="point", shape=95, size=6, color="black", fill="black")
XylV.meanArea_graph2 <- XylV.meanArea_graph2 + stat_compare_means(aes(group = VPD..kPa.), method = "aov")
XylV.meanArea_graph2 <- XylV.meanArea_graph2 + theme_bw() + theme(legend.position = "none") + ylab("Xylem Vessel mean area (pix)") + xlab("") + color_palette("aaas")
XylV.meanArea_graph2

OK - let’s make them into final pdf’s to share with Avat and her team:

library(cowplot)
## 
## Attaching package: 'cowplot'
## The following object is masked from 'package:ggpubr':
## 
##     get_legend
pdf("Graphs_traits_raw_btw_genotype_comparison.pdf", width = 10, height = 10)
plot_grid(Root_graph, Xyl.a_graph, XylVA_graph, Vessel.no_graph, labels = "AUTO", cols = 2)
## Warning in plot_grid(Root_graph, Xyl.a_graph, XylVA_graph, Vessel.no_graph, :
## Argument 'cols' is deprecated. Use 'ncol' instead.
dev.off()
## quartz_off_screen 
##                 2
pdf("Graphs_traits_raw_btw_treatment_comparison.pdf", width = 15, height = 10)
plot_grid(Root_graph2, Xyl.a_graph2, XylVA_graph2, Vessel.no_graph2, labels = "AUTO", cols = 2)
## Warning in plot_grid(Root_graph2, Xyl.a_graph2, XylVA_graph2,
## Vessel.no_graph2, : Argument 'cols' is deprecated. Use 'ncol' instead.
dev.off()
## quartz_off_screen 
##                 2
pdf("Graphs_traits_derived_btw_genotype_comparison.pdf", width = 5, height = 15)
plot_grid(Xyl.p.root_graph, XylV.p.Xyl_graph, XylV.meanArea_graph, labels = "AUTO", cols = 1)
## Warning in plot_grid(Xyl.p.root_graph, XylV.p.Xyl_graph, XylV.meanArea_graph, :
## Argument 'cols' is deprecated. Use 'ncol' instead.
dev.off()
## quartz_off_screen 
##                 2
pdf("Graphs_traits_derived_btw_treatment_comparison.pdf", width = 7.5, height = 15)
plot_grid(Xyl.p.root_graph2, XylV.p.Xyl_graph2, XylV.meanArea_graph2, labels = "AUTO", cols = 1)
## Warning in plot_grid(Xyl.p.root_graph2, XylV.p.Xyl_graph2,
## XylV.meanArea_graph2, : Argument 'cols' is deprecated. Use 'ncol' instead.
dev.off()
## quartz_off_screen 
##                 2