---
title: "Hanlon_2023_Figure_6"
author: "Meredith Hanlon"
date: "2023-04-11"
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")
    )
}


```

# Figure 6 A



```{r}
soil <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_5_1")

s2 <- soil %>% 
  filter(Element == "Sr" & Line == "K12") %>% 
  filter(Core_1 == "5" & Core_2 == "8") %>% 
  mutate(d = as.numeric(Depth) *10) %>% 
  filter(week %in% c(6, 8, 10))

s2$d <- as.character(s2$d)
s2$week <- factor(s2$week, levels = c(6, 8, 10),
                  labels = c("6wk", "8wk", "10wk"))

p1 <- ggplot(s2) +
  geom_boxplot(aes(x = Net, y = d, fill = week)) +
  theme_mh() +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3")) + 
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3")) + 
  scale_y_discrete(limits=rev) +
  ylab("Soil depth (cm)") +
  xlab("XRF Counts") +
  theme(legend.position = "top")
p1

ggsave("Figure5A.png", width = 2, height = 3, units = "in")
```

# Figure 5 B-D

```{r}

cl <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_5_2")

cl$leaf_no <- as.numeric(cl$leaf_no)

lsr1 <- cl %>% 
  dplyr::select(Core_Position, week, Genotype, Subplot, plant, plant_position, leaf_no, leaf_rel_ear) %>% 
  group_by(Core_Position, plant, week, Genotype, Subplot, plant_position) %>% 
  slice_min(leaf_no, n = 3) %>% 
  mutate(leaf_group = "old")

lsr2 <- cl %>% 
  dplyr::select(Core_Position, week, Genotype, Subplot, plant, plant_position,leaf_no, leaf_rel_ear) %>% 
  group_by(Core_Position, plant, week, Genotype, Subplot, plant_position) %>% 
  slice_max(leaf_no, n = 3) %>% 
  mutate(leaf_group = "young")

lsr3 <- cl %>% 
  dplyr::select(Core_Position, week, Genotype, Subplot, plant, plant_position,leaf_no, leaf_rel_ear) %>% 
  group_by(Core_Position, plant, week, Genotype, Subplot, plant_position) %>% 
  filter(leaf_rel_ear %in% c(-1, 0, 1)) %>% 
  mutate(leaf_group = "ear")

lsr1$Subplot <- as.character(lsr1$Subplot)
lsr2$Subplot <- as.character(lsr2$Subplot)
lsr3$Subplot <- as.character(lsr3$Subplot)

lsr <- bind_rows(lsr1, lsr2) %>% bind_rows(., lsr3) %>% distinct()

cl$Subplot <- as.character(cl$Subplot)

lsrp <- lsr %>% 
  dplyr::select(-Core_Position) %>% 
  ungroup() %>% 
  dplyr::select(week, Genotype, Subplot, plant, plant_position, leaf_group, leaf_no) %>% 
  left_join(cl %>% ungroup() %>% dplyr::select(week, Genotype, Subplot, plant_position, leaf_no, Element, leaf_rel_ear, Net)) %>% 
  filter(Element == "Sr") %>% 
  dplyr::select(-leaf_no) %>% 
  distinct()

lsrp$leaf_group = factor(lsrp$leaf_group, levels = c("old", "ear", "young"))
lsrp$week = factor(lsrp$week, levels = c("6wk", "8wk", "10wk"))

lsrp$plant_position <- factor(lsrp$plant_position, 
                             levels = c("At Sr injection", "In row", "Out of row"),
                             labels = c("At", "In", "Out"))

stat_test <- lsrp %>% 
  group_by(leaf_group, week) %>% 
  dunn_test(Net ~ plant_position)

st <- stat_test %>% add_y_position() %>% 
  filter(p.adj.signif != "ns")

lsrp %>% 
  ggplot() +
  geom_boxplot(aes(x = plant_position, y = Net, fill = plant_position)) +
  facet_grid(leaf_group~ week) +
  theme_mh() +
  scale_fill_manual(values = c( "#7F3C8D", "#11A579", "#3969AC"), guide = "none") +
  ylab("XRF Counts") +
  xlab("Plant position relative to Sr injection") +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -10),
        legend.key.size = unit(0.8, 'lines')) +
  stat_pvalue_manual(st) +
  scale_y_continuous(expand = expansion(mult = c(0.00, 0.1)))

ggsave("Figure5B.png", width = 5, height = 3, units = "in")

cl2 <- cl %>% ungroup() %>% 
  dplyr::select(week, Genotype, root_trait, root_length) 

cl2$root_trait <- factor(cl2$root_trait, levels = c( "total_length", "shallow20",  "deep30", "deep40", "mid", "d2s"),
                        labels = c("Total length", "<20 cm", ">30 cm", ">40 cm", "20<x<40 cm", "Deep to shallow ratio"))

stat_test <- cl2 %>% 
  group_by(root_trait) %>% 
  dunn_test(root_length ~ week)

st <- stat_test %>% 
  add_xy_position(x = "week", scales = "free", step.increase = 0.2) %>% 
  filter(p.adj.signif != "ns")

cl2 %>% 
  ggplot(mapping = aes(x = week, y = root_length)) +
  geom_boxplot(aes(fill = week)) +
  facet_wrap(.~root_trait, scales = "free_y", ncol = 3) +
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3")) +
  ylab("Root Length (cm)") +
  theme_mh() +
  theme(legend.position = "none",
        axis.title.x = element_blank()) +
    stat_pvalue_manual(st) +
  scale_y_continuous(expand = expansion(mult = c(0.00, 0.1)))

ggsave("Figure5C.png", width = 7, height = 2.5, units = "in")


cl$root_trait <- factor(cl$root_trait, levels = c( "total_length", "shallow20",  "deep30", "deep40", "mid", "d2s"),
                        labels = c("Total length", "<20 cm", ">30 cm", ">40 cm", "20<x<40 cm", "Deep to shallow ratio"))

rla <- cl %>%
  dplyr::select(week, Genotype, Subplot, root_trait, root_length) %>% 
  group_by(week, Genotype, Subplot, root_trait) %>% 
  summarise(root_length_mean = mean(root_length, na.rm = TRUE))

rla$Subplot <- as.character(rla$Subplot)

# get leaf groups

lsr1 <- cl %>% 
  filter(plant %in% c("p5", "p8")) %>% 
  unite(id, week, Genotype, Subplot, plant, leaf_no, leaf_rel_ear, sep = "_") %>% 
  separate(id, into = c("week", "Genotype", "Subplot", "plant", "leaf_no", 'leaf_rel_ear'), sep = "_") %>% 
  dplyr::select(week, Genotype, Subplot, plant, leaf_no, leaf_rel_ear) %>% 
  slice_min(leaf_no, n = 3) %>% 
  mutate(leaf_group = "old")
lsr2 <- cl %>% 
  filter(plant %in% c("p5", "p8")) %>% 
  dplyr::select(week, Genotype, Subplot, plant, leaf_no, leaf_rel_ear) %>% 
  slice_max(leaf_no, n = 3) %>% 
  mutate(leaf_group = "young")
lsr3 <- cl %>% 
  filter(plant %in% c("p5", "p8")) %>% 
  dplyr::select(week, Genotype, Subplot, plant, leaf_no, leaf_rel_ear) %>% 
  filter(leaf_rel_ear %in% c(-1, 0, 1)) %>% 
  mutate(leaf_group = "ear")
lsr2$Subplot <- as.character(lsr2$Subplot)
lsr3$Subplot <- as.character(lsr3$Subplot)
lsr2$leaf_no <- as.character(lsr2$leaf_no)
lsr3$leaf_no <- as.character(lsr3$leaf_no)
lsr2$leaf_rel_ear<- as.character(lsr2$leaf_rel_ear)
lsr3$leaf_rel_ear <- as.character(lsr3$leaf_rel_ear)
lsr <- bind_rows(lsr1, lsr2) %>% bind_rows(., lsr3)

rla$Subplot <- as.character(rla$Subplot)
cl$Subplot <- as.character(cl$Subplot)
cl$leaf_no <- as.character(cl$leaf_no)
cl$leaf_rel_ear <- as.character(cl$leaf_rel_ear)

lsr4 <- lsr %>% 
  ungroup() %>% 
  left_join(cl %>% ungroup() %>% dplyr::select(week, Genotype, Subplot, leaf_no, Element, leaf_rel_ear, Net)) %>% 
  filter(Element == "Sr") %>% 
  group_by(week, Genotype, Subplot, leaf_group) %>% 
  # remove extreme outliers
  summarise(sr_mean = mean(Net, na.rm = TRUE)) %>% 
  left_join(rla) 


lsr4$week <- factor(lsr4$week, 
                    levels = c("6wk", "8wk", "10wk"))

lsr4 %>% 
  filter(leaf_group == "ear") %>% 
  ggplot() +
  geom_smooth(aes(x = sr_mean, y = root_length_mean), method = "lm", se = FALSE, color = "grey40") +
  geom_point(aes(x = sr_mean, y = root_length_mean, color = week)) +
  theme_mh() +
  facet_wrap(root_trait~., scales = "free", ncol = 3) +
  stat_cor(aes(x = sr_mean, y = root_length_mean, label = paste(..rr.label.., ..p.label.., sep = "~`,`~")), method = "pearson") +
  theme(legend.position = "top") +
  theme_mh() + 
  xlab("Sr Counts") +
  ylab("Root trait measurement") +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3")) +
  theme(legend.position = "none",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -15),
        axis.ticks.x = element_blank(),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2)))

ggsave("Figure5D.png", width = 7, height = 3, units = "in")
```