---
title: "Hanlon_2023_Figure_4"
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 4A

msr <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_1")
msr$DAG <- as.character(msr$DAG)
msr$ppm <- as.numeric(msr$ppm)
msr2 <- msr %>% 
  filter(Depth1 >= 80 & Depth1 < 140)

p1 <- ggplot(msr2) +
  geom_point(aes(y = ppm, x = Depth1, color = DAG)) +
  theme_mh() +
  geom_smooth(aes(y = ppm, x = Depth1, color = DAG, fill = DAG)) +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3")) + 
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3")) + 
  coord_flip() +
  scale_x_reverse() +
  xlab("Column depth (cm)") +
  theme(legend.position = "top")
p1
ggsave("Figure4A.png", width = 2.5, height = 3)

# Figure 4B

crm <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_2")
rs <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_3")

crm$DAG <- as.character(crm$DAG)

rs$DAG <- as.character(rs$DAG)
rs$Depth <- factor(rs$Depth, levels = c("85_90", "90_95", "95_100", "100_105", "105_110", "110_120", "120_160"), 
                   labels = c("85-90", "90-95", "95-100", "100-105", "105-110", "110-120", "120-160"))

rs2 <- rs %>% 
  filter(root_class == "total_cm") %>% 
  group_by(Genotype, DAG, Depth, Depth1) %>% 
  summarise(SE = stderror(length),
            mean_length = mean(length))
p2 <- ggplot(data = rs2) +
  geom_point(mapping = aes(x = Depth1, y = mean_length, color = Genotype), alpha = 0.5, size = 2) +
  geom_line(mapping = aes(x = Depth1, y = mean_length, color = Genotype), show.legend = F) + 
  geom_errorbar(mapping = aes(x = Depth1, y = mean_length, ymin = mean_length - SE, 
                    ymax = mean_length + SE, color = Genotype), 
                alpha = 0.5, width = 0, show.legend = F) + 
  facet_wrap(.~DAG, ncol = 3, scales = "free_x") +
  scale_x_reverse() +
  coord_flip() +
  scale_color_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701"),
                     name = "Genotype") +
  scale_fill_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701"),
                    guide = "none") +
  xlab("Column Depth (cm)") +
  ylab("Total Root Length (cm)") +
  theme_mh() +
  theme(legend.position = "none",
        strip.text = element_blank())
p3 <- ggplot(crm) +
  geom_boxplot(aes(x = Genotype, y = cr_mass, fill = Genotype), alpha = 0.5) +
  facet_wrap(.~DAG, ncol = 3, scales = "free_x") +
  scale_fill_manual(values = c("#7F3C8D", "#11A579", "#3969AC", "#F2B701")) +
  ylab("Root mass (g)") +
  theme_mh() +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -10),
        legend.key.size = unit(0.8, 'lines'),
        axis.title.x = element_blank(),
        axis.text.x = element_blank(),
        axis.ticks.x = element_blank())
p3 + p2 + 
  plot_layout(nrow = 2, heights = c(2, 3))

ggsave("Figure4B.png", width = 5, height = 3)
```

# Figure 4C


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

sr <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_5") 
xrf_final <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_4_6") 

sr$DAG <- as.character(sr$DAG)

rs2 <- root_length %>% 
  pivot_longer(total_cm:thick, names_to = "trait", values_to = "value") %>% 
  pivot_wider(names_from = c(trait, Depth1, Depth2), values_from = value) %>% 
  group_by(Column) %>% 
  summarise(total_shallow = sum(total_cm_85_90, total_cm_90_95, na.rm = TRUE),
            total_at = sum(total_cm_95_100, total_cm_100_105, total_cm_105_110, na.rm = TRUE),
            total_deep = sum(total_cm_110_120, total_cm_120_160, na.rm = TRUE),
            fine_shallow = sum(fine_85_90, fine_90_95, na.rm = TRUE),
            fine_at = sum(fine_95_100, fine_100_105, fine_105_110, na.rm = TRUE),
            fine_deep = sum(fine_110_120, fine_120_160, na.rm = TRUE),
            moderate_shallow = sum(moderate_85_90, moderate_90_95, na.rm = TRUE),
            moderate_at = sum(moderate_95_100, moderate_100_105, moderate_105_110, na.rm = TRUE),
            moderate_deep = sum(moderate_110_120, moderate_120_160, na.rm = TRUE),
            thick_shallow = sum(thick_85_90, thick_90_95, na.rm = TRUE),
            thick_at = sum(thick_95_100, thick_100_105, thick_105_110, na.rm = TRUE),
            thick_deep = sum(thick_110_120, thick_120_160, na.rm = TRUE)) %>% 
  left_join(crm, .) %>% 
  dplyr::select(Column, Genotype, DAG, cr_mass, total_shallow:thick_deep) %>% 
  pivot_longer(cr_mass:thick_deep, names_to = "root_trait", values_to = "value") %>% 
  left_join(., sr)

xic <- rs2 %>% 
  ungroup() %>% 
  group_by(root_trait, Sr, DAG) %>% 
  mutate(r2 = round((cor(value, count))^2, 2)) %>% 
  dplyr::select(root_trait, Sr, DAG, r2) %>% 
  distinct()

xic2 <- rs2 %>% 
  ungroup() %>% 
  group_by(root_trait, Sr) %>% 
  mutate(r2 = round((cor(value, count))^2, 2)) %>% 
  dplyr::select(root_trait, Sr, r2) %>% 
  distinct()


# Leaves 8-12

sr2 <-  xrf_final %>% 
  dplyr::filter(Element == "Sr") %>% 
  dplyr::filter(leaf %in% c("L8", "L9", "L10", "L11", "L12")) %>% 
  left_join(., crm %>% dplyr::select(Column, DAG, Genotype)) %>% 
  group_by(Column, Genotype, DAG) %>% 
  summarise(plant_total_sr = sum(Net),
            plant_mean_sr = mean(Net)) %>% 
  pivot_longer(plant_total_sr:plant_mean_sr, names_to = "Sr", values_to = "count") %>% 
  drop_na()

sr2$DAG <- as.character(sr2$DAG)

rs2 <- root_length %>% 
  pivot_longer(total_cm:thick, names_to = "trait", values_to = "value") %>% 
  pivot_wider(names_from = c(trait, Depth1, Depth2), values_from = value) %>% 
  group_by(Column) %>% 
  summarise(total_shallow = sum(total_cm_85_90, total_cm_90_95, na.rm = TRUE),
            total_at = sum(total_cm_95_100, total_cm_100_105, total_cm_105_110, na.rm = TRUE),
            total_deep = sum(total_cm_110_120, total_cm_120_160, na.rm = TRUE),
            fine_shallow = sum(fine_85_90, fine_90_95, na.rm = TRUE),
            fine_at = sum(fine_95_100, fine_100_105, fine_105_110, na.rm = TRUE),
            fine_deep = sum(fine_110_120, fine_120_160, na.rm = TRUE),
            moderate_shallow = sum(moderate_85_90, moderate_90_95, na.rm = TRUE),
            moderate_at = sum(moderate_95_100, moderate_100_105, moderate_105_110, na.rm = TRUE),
            moderate_deep = sum(moderate_110_120, moderate_120_160, na.rm = TRUE),
            thick_shallow = sum(thick_85_90, thick_90_95, na.rm = TRUE),
            thick_at = sum(thick_95_100, thick_100_105, thick_105_110, na.rm = TRUE),
            thick_deep = sum(thick_110_120, thick_120_160, na.rm = TRUE)) %>% 
  left_join(crm, .) %>% 
  dplyr::select(Column, Genotype, DAG, cr_mass, total_shallow:thick_deep) %>% 
  pivot_longer(cr_mass:thick_deep, names_to = "root_trait", values_to = "value") %>% 
  left_join(., sr2) %>% 
  drop_na()


xic <- rs2 %>% 
  ungroup() %>% 
  group_by(root_trait, Sr, DAG) %>% 
  dplyr::summarize(cos = stats::cor.test(value, count)$estimate,
                   p = stats::cor.test(value, count)$p.value) %>%
  mutate(r2 = round(cos^2, 2)) %>% 
  dplyr::select(root_trait, Sr, DAG, r2, p) %>% 
  distinct()

xicb <- rs2 %>% 
  ungroup() %>% 
  group_by(root_trait, Sr) %>% 
  dplyr::summarize(cos = stats::cor.test(value, count)$estimate,
                   p = stats::cor.test(value, count)$p.value) %>%
  mutate(r2 = round(cos^2, 2)) %>% 
  dplyr::select(root_trait, Sr, r2, p) %>% 
  distinct()

# just leaf 10

sr2 <-  xrf_final %>% 
  dplyr::filter(Element == "Sr") %>% 
  dplyr::filter(leaf %in% c("L10")) %>% 
  left_join(., crm %>% dplyr::select(Column, DAG, Genotype)) %>% 
  group_by(Column, Genotype, DAG) %>% 
  summarise(plant_total_sr = sum(Net),
            plant_mean_sr = mean(Net)) %>% 
  pivot_longer(plant_total_sr:plant_mean_sr, names_to = "Sr", values_to = "count") %>% 
  drop_na()

sr2$DAG <- as.character(sr2$DAG)

rs2 <- root_length %>% 
  pivot_longer(total_cm:thick, names_to = "trait", values_to = "value") %>% 
  pivot_wider(names_from = c(trait, Depth1, Depth2), values_from = value) %>% 
  group_by(Column) %>% 
  summarise(total_shallow = sum(total_cm_85_90, total_cm_90_95, na.rm = TRUE),
            total_at = sum(total_cm_95_100, total_cm_100_105, total_cm_105_110, na.rm = TRUE),
            total_deep = sum(total_cm_110_120, total_cm_120_160, na.rm = TRUE),
            fine_shallow = sum(fine_85_90, fine_90_95, na.rm = TRUE),
            fine_at = sum(fine_95_100, fine_100_105, fine_105_110, na.rm = TRUE),
            fine_deep = sum(fine_110_120, fine_120_160, na.rm = TRUE),
            moderate_shallow = sum(moderate_85_90, moderate_90_95, na.rm = TRUE),
            moderate_at = sum(moderate_95_100, moderate_100_105, moderate_105_110, na.rm = TRUE),
            moderate_deep = sum(moderate_110_120, moderate_120_160, na.rm = TRUE),
            thick_shallow = sum(thick_85_90, thick_90_95, na.rm = TRUE),
            thick_at = sum(thick_95_100, thick_100_105, thick_105_110, na.rm = TRUE),
            thick_deep = sum(thick_110_120, thick_120_160, na.rm = TRUE)) %>% 
  left_join(crm, .) %>% 
  dplyr::select(Column, Genotype, DAG, cr_mass, total_shallow:thick_deep) %>% 
  pivot_longer(cr_mass:thick_deep, names_to = "root_trait", values_to = "value") %>% 
  left_join(., sr2) %>% 
  drop_na()

xic2 <- rs2 %>% 
  ungroup() %>% 
  group_by(root_trait, Sr, DAG) %>% 
  dplyr::summarize(cos = stats::cor.test(value, count)$estimate,
                   p = stats::cor.test(value, count)$p.value) %>%
  mutate(r2 = round(cos^2, 2)) %>% 
  dplyr::select(root_trait, Sr, DAG, r2, p) %>% 
  distinct()

xicb2 <- rs2 %>% 
  ungroup() %>% 
  group_by(root_trait, Sr) %>% 
  dplyr::summarize(cos = stats::cor.test(value, count)$estimate,
                   p = stats::cor.test(value, count)$p.value) %>%
  mutate(r2 = round(cos^2, 2)) %>% 
  dplyr::select(root_trait, Sr, r2, p) %>% 
  distinct()

# combine and plot

xicbb <- xicb %>% mutate(DAG = "all")
xicbb2 <- xicb2 %>% mutate(DAG = "all")
xx <- xic %>%
  bind_rows(., xicbb)  %>% 
  mutate(leaf = "L8_12") 

xx2 <- xic2 %>% 
  bind_rows(., xicbb2) %>% 
  mutate(leaf = "L10")

xx3 <- bind_rows(xx, xx2)

xx4 <- xx3 %>% 
  filter(Sr == "plant_total_sr") %>% 
  filter(root_trait %in% c("cr_mass", "total_shallow", "total_at", "total_deep", "fine_shallow", "fine_at", "fine_deep")) %>% 
  mutate(depth_class = case_when(root_trait %in% c("cr_mass", "total_shallow", "fine_shallow") ~ "Shallow",
                                 root_trait %in% c("total_at", "fine_at") ~ "At Sr",
                                 TRUE ~ "Deep"))
# Make plot - have l8-12 and l10 for all, then split up days 
xa <- xx4 %>% filter(DAG == "all")

xa$leaf <- factor(xa$leaf, levels = c("L10", "L8_12"), labels = c("L10", "L8-12"))
xa$depth_class <- factor(xa$depth_class, levels = c("Shallow", "At Sr", "Deep"))
xa$root_trait <- factor(xa$root_trait, 
                        levels = c("cr_mass", "total_shallow", "fine_shallow",
                                   "total_at", "fine_at",
                                   "total_deep", "fine_deep"),
                        labels = c("Shallow mass", "Total length", "Fine length",
                                   "Total length", "Fine length",
                                   "Total length", "Fine length"))

x1 <- xa %>% 
   mutate(sig = case_when(p < 0.05 ~ "sig1",
                         p >= 0.05 & p < 0.1 ~ "sig2")) %>% 
  mutate(r3 = case_when(sig %in% c("sig1", "sig2") ~ r2)) %>% 
  ggplot(aes(x = leaf, y = root_trait, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(limits = c(0,1), expand = c(0,0), 
                     name = expression(R^2)) +
  geom_text(aes(label = round(r3, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  facet_wrap(depth_class~.,  scales = "free_y", ncol = 1, strip.position = "left") +
  theme_mh() +
  scale_x_discrete(position = "top") +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        strip.placement = "outside",
        legend.position = "top",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -15),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = margin(l = -0),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2))) 

xa2 <- xx4 %>% filter(DAG != "all") %>% filter(leaf == "L8_12")

xa2$depth_class <- factor(xa2$depth_class, levels = c("Shallow", "At Sr", "Deep"))
xa2$root_trait <- factor(xa2$root_trait, 
                        levels = c("cr_mass", "total_shallow", "fine_shallow",
                                   "total_at", "fine_at",
                                   "total_deep", "fine_deep"),
                        labels = c("Shallow mass", "Total length", "Fine length",
                                   "Total length", "Fine length",
                                   "Total length", "Fine length"))

x2 <- xa2 %>%  
   mutate(sig = case_when(p < 0.05 ~ "sig1",
                         p >= 0.05 & p < 0.1 ~ "sig2")) %>% 
  mutate(r3 = case_when(sig %in% c("sig1", "sig2") ~ r2)) %>% 
  ggplot(aes(x = DAG, y = root_trait, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
    scale_fill_viridis(limits = c(0,1), expand = c(0,0), 
                     name = expression(R^2)) +
  geom_text(aes(label = round(r3, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  facet_wrap(depth_class~.,  scales = "free_y", ncol = 1) +
  theme_mh() +
  scale_x_discrete(position = "top") +
  theme(axis.title.y = element_blank(),
        axis.title.x = element_blank(),
        axis.text.y = element_blank(),
        strip.text = element_blank(), 
        strip.placement = "outside",
        legend.position = "top",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -15),
        axis.ticks.x = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = margin(l = -0),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2))) 

x1 + x2 +
  plot_layout(nrow = 1, widths = c(2, 4), guides = "collect") & theme(legend.position = 'top')

ggsave("Figure4C.png", width = 4, height = 3)
```



