---
title: "Supplemental Figure 9"
author: "Meredith Hanlon"
date: "2023-04-13"
output: html_document
---
```{r global_options, echo = FALSE}
knitr::opts_chunk$set(warning=FALSE, message=FALSE, echo=FALSE, error = TRUE)
```

```{r}
library(tidyverse)
library(rstatix)
library(ggpubr)
library(readxl)
library(stats)
library(viridis)
library(patchwork)
library(ggpmisc)
library(ggthemes)
library(ggrepel)
library(FactoMineR)
library(factoextra)
library(corrplot)
library(mclust)
library(broom)
library(RColorBrewer)
library(ggh4x)

stderror <- function(x) sd(x, na.rm = TRUE) / sqrt(length(na.omit(x)))

theme_mh <- function(){ 
    font <- "sans"   #assign font family up front
    theme(
      panel.grid.major = element_blank(),
      panel.grid.minor = element_blank(),
      panel.background = element_blank(),
      panel.border = element_rect(color = "black", fill = NA),
      axis.ticks = element_line(color = "black"),   
      axis.line = element_line(color = "black"),
      strip.background = element_blank(),
      legend.key = element_rect(fill = NA),
      text = element_text(family = font, size = 10, color = "black"),
      axis.text = element_text(family = font, size = 8, color = "black")
    )
}


```

# Supplemental Figure 9

```{r}
bm <- read_excel("Hanlon_2023_data.xlsx", sheet = "SFigure_9_1")
root_length <- read_excel("Hanlon_2023_data.xlsx", sheet = "SFigure_9_2")
xrf <- read_excel("Hanlon_2023_data.xlsx", sheet = "SFigure_9_3")
```


# Biomass

```{r}
bm2 <- bm %>% 
  pivot_longer(Leaf:Total, names_to = "tissue", values_to = "weight") %>% 
  dplyr::filter(tissue == "Total") %>% 
  filter(treatment %in% c("c", "t1", "t2", "t3", "t6", "t7"))

bm2$treatment2 <- factor(bm2$treatment,
                   levels = c("c", "t1", "t2", "t3", "t6", "t7"),
                   labels = c("None", "0.44 mM Sr", "1.3 mM Sr", "2.6 mM Sr",
                              "0.44 mM Rb", "1.3 mM Rb"))

stats <- bm2 %>% 
  dunn_test(weight ~ treatment2)


bm2$treatment2 <- factor(bm2$treatment,
                   levels = c("c", "t1", "t2", "t3", "t6", "t7"),
                   labels = c("None", "0.44 mM Sr", "1.3 mM Sr", "2.6 mM Sr",
                              "0.44 mM Rb", "1.3 mM Rb"))
bm3 <- bm2 %>% 
  mutate(Trmt = case_when(treatment == "t3" ~ "2.6 mM",
                          treatment %in% c("t1", "t6") ~ "0.44 mM",
                          treatment %in% c("t2", "t7") ~ "1.3 mM")) %>% 
  mutate(chem = case_when(treatment %in% c("t1", "t2", "t3") ~ "Sr",
                          treatment %in% c("t6", "t7") ~ "Rb",
                          treatment == "c" ~ "None")) %>% 
  mutate(Trmt = replace_na(Trmt, "")) %>% 
  as.data.frame()

p1 <- ggplot(bm3, aes(x = interaction(Trmt, chem, sep = "&"), y = weight)) +
  geom_boxplot() +
  ggbeeswarm::geom_beeswarm(aes(color = treatment)) +
  theme_mh() +
  ylab("Shoot dry biomass (g)") +
  theme(legend.position = "none",
        axis.title.x = element_blank()) +
  scale_color_tableau() +
  guides(x = guide_axis_nested(delim = "&"))
p1

ggsave("Supplemental_Figure_9A.png", width = 4, height = 2)
```

```{r}
rl2 <- root_length 

rl2$treatment2 <- factor(rl2$treatment,
                   levels = c("None", "0.44 mM Sr", "1.3 mM Sr", "2.6 mM Sr",
                              "0.44 mM Rb", "1.3 mM Rb"))

stats <- rl2 %>% 
  group_by(root_class) %>% 
  dunn_test(length ~ treatment2)


rl3 <- rl2 %>% 
  mutate(Trmt = case_when(treatment == "2.6 mM Sr" ~ "2.6 mM",
                          treatment %in% c("0.44 mM Rb", "0.44 mM Sr") ~ "0.44 mM",
                          treatment %in% c("1.3 mM Sr", "1.3 mM Rb") ~ "1.3 mM")) %>% 
  mutate(chem = case_when(treatment %in% c("0.44 mM Sr", "1.3 mM Sr", "2.6 mM Sr") ~ "Sr",
                          treatment %in% c("0.44 mM Rb", "1.3 mM Rb") ~ "Rb",
                          treatment == "None" ~ "None")) %>% 
  mutate(Trmt = replace_na(Trmt, "")) %>% 
  as.data.frame()

p1 <- ggplot(rl3) +
  geom_boxplot(aes(x = interaction(Trmt, chem, sep = "&"), y = length)) +
  ggbeeswarm::geom_beeswarm(aes(x = interaction(Trmt, chem, sep = "&"), y = length, color = root_class)) +
  theme_mh() +
  facet_wrap(.~root_class, scales = "free_y") +
  theme(legend.position = "none",
        axis.title.x = element_blank()) +
  scale_color_tableau() +
  guides(x = guide_axis_nested(delim = "&"))
p1

ggsave("Supplemental_Figure_9B.png", width = 7.5, height = 3)
```

```{r}
x3 <- xrf %>% dplyr::filter(Element == "Sr") %>% 
  dplyr::filter(treatment %in% c("c", "t1", "t2", "t3")) %>%
  mutate(experiment_group = "Sr")

x4 <- xrf %>% dplyr::filter(Element == "Rb") %>% 
  dplyr::filter(treatment %in% c("c", "t6", "t7")) %>% 
  mutate(experiment_group = "Rb")
  
x5 <- bind_rows(x3, x4)
x5$DAP <- as.character(x5$DAP)

x5$treatment <- factor(x5$treatment,
                   levels = c("c", "t1", "t2", "t3", "t6", "t7"),
                   labels = c("None", "0.44 mM Sr", "1.3 mM Sr", "2.6 mM Sr",
                              "0.44 mM Rb", "1.3 mM Rb"))

stats <- x5 %>% 
  ungroup() %>% 
  group_by(Element, experiment_group, DAP) %>% 
  t_test(Net ~ treatment, ref.group = "None") %>% 
  filter(p.adj.signif != "ns")  %>%  
  add_x_position() %>% 
  mutate(y.position = c(15500))


p1 <- ggplot(x5) +
  geom_boxplot(aes(x = treatment, y = Net, fill = DAP)) +
  theme_mh () +
  theme(legend.position = "none",
        axis.title.x = element_blank()) +
  scale_fill_tableau() +
  facet_grid2(vars(experiment_group), vars(DAP),remove_labels = "y",
              scales = "free", independent = "x") +
  stat_pvalue_manual(stats)+
  scale_y_continuous(expand = c(0.2, 0))
p1


ggsave("Supplemental_Figure_9C.png", width = 7.5, height = 3)

```
