---
title: "Figure_1_Hanlon_2023"
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(broom)
library(tidyverse)
library(rstatix)
library(ggpubr)
library(readxl)
library(stats)
library(viridis)
library(patchwork)

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 1A 

root depth by genotype for the 2020 XRF data

```{r, fig.width = 2.5, fig.height = 4}
d <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1A")
d2 <- d %>% 
  dplyr::rename(Genotype = genotype, Rep = rep, Week = time) %>% 
  mutate(Genotype = case_when(Genotype == "99d" ~ "DKC49 x 94RIB",
                              Genotype == "85d" ~ "DKC35 x 88RIB",
                              TRUE ~ Genotype))


# Change the names for the hybrids - 85d = DKC35 X 88RIB, 99d = DKC49 X 94RIB

cs <- d2 %>% 
  dplyr::select(Genotype, Rep, Week, seg_length_seg_1:seg_length_seg_6) %>% 
  pivot_longer(seg_length_seg_1:seg_length_seg_6, names_to = "trait", values_to = "root_length") %>% 
  separate(trait, into=c("x","y","z","Segment")) %>% 
  dplyr::select(-x, -y, -z) %>% 
  group_by(Genotype, Week, Segment) %>% 
  summarise(SE = stderror(root_length),
            root_length = mean(root_length))


# get d95 and d50
cs2 <- d2 %>%   
  dplyr::select(Genotype, Rep, Week, d95) %>% 
  pivot_longer(d95, names_to = "trait", values_to = "dtr") %>%
  group_by(Genotype, Week, trait) %>% 
  summarise(SE = stderror(dtr),
            dtrs = mean(dtr))

cs$Segment <- 10*as.numeric(cs$Segment)
cs$Genotype <- factor(cs$Genotype, levels = c("B73", "PHB47", "DKC35 x 88RIB", "DKC49 x 94RIB"))
cs$Week <- factor(cs$Week, levels = c("6w", "8w", "10w"))

cs2$Genotype <- factor(cs2$Genotype, levels = c("B73", "PHB47", "DKC35 x 88RIB", "DKC49 x 94RIB"))
cs2$Week <- factor(cs2$Week, levels = c("6w", "8w", "10w"))

ggplot() +
  geom_vline(data = cs2, mapping = aes(xintercept = dtrs, color = Week), linetype = 2, show.legend = F) +
  geom_rect(data = cs2, aes(xmin = dtrs-SE, xmax = dtrs+SE, ymin = 0, 
                            ymax = 800, fill = Week), alpha = 0.07) +
  geom_point(data = cs, mapping = aes(x = Segment, y = root_length, color = Week), alpha = 0.5, size = 2) +
  geom_line(data = cs, mapping = aes(x = Segment, y = root_length, color = Week), show.legend = F) + 
  geom_errorbar(data = cs, mapping = aes(x = Segment, y = root_length, ymin = root_length - SE, 
                    ymax = root_length + SE, color = Week), 
                alpha = 0.5, width = 0, show.legend = F) + 
  scale_x_reverse() +
  facet_wrap(.~Genotype, ncol = 2) +
  coord_flip() +
  scale_color_manual(values = c("#1B9E77", "#D95F02", "#7570B3"),
                     name = "Week") +
  scale_fill_manual(values = c("#1B9E77", "#D95F02", "#7570B3"),
                    guide = "none") +
  xlab("Core Segment Depth (cm)") +
  ylab("Root Length (cm)") +
  theme_mh() +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -10),
        legend.key.size = unit(0, 'lines'))

ggsave("Figure1A.png", width = 2.50, height = 4, units = "in")
```

# Figure 1B

Compare XRF and ICP


```{r, fig.width = 5, fig.height = 4}
xio <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1B")

xio3 <- xio %>% 
  pivot_longer(cols = c(ICP, XRF), names_to = "method", values_to = "elemental_value") %>% 
  group_by(method) %>% 
  mutate(scale_ev = scale(elemental_value, center = FALSE)) %>% 
  filter(!Element %in% c("Ni", "Pb")) %>% 
  unite(el.m, Element, method, sep = ", ", remove = FALSE)
xio3$Week <- factor(xio3$Week, levels = c("6w", "8w", "10w"))
xio3$Genotype <- factor(xio3$Genotype, levels = c("B73", "PHB47", "DKC35 x 88RIB", "DKC49 x 94RIB"))

# scale all xrf by 1000
# only show select elements because of space constraints

xio3 %>% 
  mutate(elemental_value = case_when(method == "XRF" ~ elemental_value/1000,
                                     TRUE ~ elemental_value)) %>% 
  filter(Element %in% c("Sr", "K")) %>% 
  ggplot() +
  geom_boxplot(aes(x = Genotype, y = elemental_value, color = Week, fill = method),
               outlier.size = 0.5, lwd = 0.2) +
  facet_wrap(.~el.m, ncol = 2, scales = "free_y" ) +
  scale_color_manual(values = c("#9ECAE1", "#4292C6", "#084594")) +
  scale_fill_manual(values = c("white", "grey80"), guide = "none")+
  theme_mh() +
  theme(legend.position = "top") +
  ylab("Counts (XRF) or ppm (ICP)") +
  theme_mh() +
  theme(legend.position = "top",
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -5),
        legend.key.size = unit(0.5, 'lines'),
        axis.text.x = element_text(angle = 45, hjust = 1),
        panel.spacing = unit(0, "lines"))

ggsave("Figure1B.png", width = 5, height = 4, units = "in")
```

# Figure 1C

Correlation of XRF, ICP and root traits

```{r}
d <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1A")
xio <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1B")
ic <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1C")
lx <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1C_2")

xic <- xio %>% 
  dplyr::select(Element, ICP, XRF) %>% 
  mutate(r2 = round((cor(XRF, ICP))^2, 2)) %>% 
  dplyr::select(Element, r2) %>% 
  mutate(method = "XRF_ICP_corr")



d2 <- d %>% 
  dplyr::rename(Genotype = genotype, Rep = rep, Week = time) %>% 
  mutate(Genotype = case_when(Genotype == "99d" ~ "DKC49 x 94RIB",
                              Genotype == "85d" ~ "DKC35 x 88RIB",
                              TRUE ~ Genotype))
d2$Genotype <- factor(d2$Genotype, levels = c("B73", "PHB47", "DKC35 x 88RIB", "DKC49 x 94RIB"))
d2$Week <- factor(d2$Week, levels = c("6w", "8w", "10w"))

d3 <- d2 %>% 
  mutate(shallow_20 = seg_length_seg_1 + seg_length_seg_2,
         deep_30 = seg_length_seg_4 + seg_length_seg_5 + seg_length_seg_6,
         deep_40 = seg_length_seg_5 + seg_length_seg_6,
         percent_shallow_20 = (seg_length_seg_1 + seg_length_seg_2)/total_root_length * 100,
         percent_deep_30 = (seg_length_seg_4 + seg_length_seg_5 + seg_length_seg_6)/total_root_length * 100,
         percent_deep_40 = (seg_length_seg_5 + seg_length_seg_6)/total_root_length * 100,
         ) %>% 
  dplyr::select(Genotype, Rep, Week, total_root_length, shallow_20, deep_30, deep_40, percent_shallow_20, percent_deep_30, percent_deep_40, d50, d70, d95)

# first get total cors for all leaves, then do each leaf group individually, pivot, and combine

cor_all_icp <- ic %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>%   
  group_by(root_trait, Element) %>% 
  dplyr::summarize(cos = stats::cor.test(ICP, root_value)$estimate,
                   p_all = stats::cor.test(ICP, root_value)$p.value) %>%
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  rename(r2_all = r2)

cor_leaves_icp <- ic %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element, leaf) %>% 
  dplyr::summarize(cos = stats::cor.test(ICP, root_value)$estimate,
                   p = stats::cor.test(ICP, root_value)$p.value) %>%
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  pivot_wider(names_from = leaf, values_from = c(r2, p))

cor_icp <- full_join(cor_all_icp, cor_leaves_icp) %>% 
  mutate(method = "ICP") %>% 
  pivot_longer(cols = c(p_all:p_young), names_to = c("test"), values_to = c("value")) %>% 
  separate(test, into = c("test", "p")) %>% 
  pivot_wider(names_from = p, values_from = value)

# xrf cors

cor_all_xrf <- lx %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element) %>% 
  dplyr::summarize(cos = stats::cor.test(value, root_value)$estimate,
                   p_all = stats::cor.test(value, root_value)$p.value) %>% 
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  rename(r2_all = r2)

cor_leaves_xrf <- lx %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element, leaf) %>% 
    dplyr::summarize(cos = stats::cor.test(value, root_value)$estimate,
                   p = stats::cor.test(value, root_value)$p.value) %>%
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  pivot_wider(names_from = leaf, values_from = c(r2, p))

cor_xrf <- full_join(cor_all_xrf, cor_leaves_xrf) %>% 
  mutate(method = "XRF") %>% 
  pivot_longer(cols = c(p_all:p_young), names_to = c("test"), values_to = c("value")) %>% 
  separate(test, into = c("test", "p")) %>% 
  pivot_wider(names_from = p, values_from = value)

inds <- full_join(cor_icp, cor_xrf) %>% 
  pivot_longer(all:young, names_to = "Leaves", values_to = "val") %>% 
  pivot_wider(names_from = test, values_from = val)

sup1 <- inds %>% 
  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 = method, y = Element, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(name = "r2", limits = c(0,1)) +
  geom_text(aes(label = round(r3, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  theme_mh() +
  facet_wrap(.~root_trait+Leaves)

ggsave("Supplemental Figure 1.png", width = 10, height = 30, units = "in")


# plot averages

d2 <- d %>% 
  dplyr::rename(Genotype = genotype, Rep = rep, Week = time) %>% 
  mutate(Genotype = case_when(Genotype == "99d" ~ "DKC49 x 94RIB",
                              Genotype == "85d" ~ "DKC35 x 88RIB",
                              TRUE ~ Genotype))
d2$Genotype <- factor(d2$Genotype, levels = c("B73", "PHB47", "DKC35 x 88RIB", "DKC49 x 94RIB"))
d2$Week <- factor(d2$Week, levels = c("6w", "8w", "10w"))

d3 <- d2 %>% 
  mutate(shallow_20 = seg_length_seg_1 + seg_length_seg_2,
         deep_30 = seg_length_seg_4 + seg_length_seg_5 + seg_length_seg_6,
         deep_40 = seg_length_seg_5 + seg_length_seg_6,
         percent_shallow_20 = (seg_length_seg_1 + seg_length_seg_2)/total_root_length * 100,
         percent_deep_30 = (seg_length_seg_4 + seg_length_seg_5 + seg_length_seg_6)/total_root_length * 100,
         percent_deep_40 = (seg_length_seg_5 + seg_length_seg_6)/total_root_length * 100,
         ) %>% 
  dplyr::select(Genotype, Rep, Week, total_root_length, shallow_20, deep_30, deep_40, percent_shallow_20, percent_deep_30, percent_deep_40, d50, d70, d95) %>% 
  ungroup() %>% 
  group_by(Genotype, Week) %>% 
  summarise(across(where(is.numeric), ~mean(.,na.rm = TRUE)))

ic <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1C_3")

cor_all_icp <- ic %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>%   
  group_by(root_trait, Element) %>% 
  dplyr::summarize(cos = stats::cor.test(ICP, root_value)$estimate,
                   p_all = stats::cor.test(ICP, root_value)$p.value) %>%
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  rename(r2_all = r2)

cor_leaves_icp <- ic %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element, leaf) %>% 
  dplyr::summarize(cos = stats::cor.test(ICP, root_value)$estimate,
                   p = stats::cor.test(ICP, root_value)$p.value) %>%
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  pivot_wider(names_from = leaf, values_from = c(r2, p))

cor_icp <- full_join(cor_all_icp, cor_leaves_icp) %>% 
  mutate(method = "ICP") %>% 
  pivot_longer(cols = c(p_all:p_young), names_to = c("test"), values_to = c("value")) %>% 
  separate(test, into = c("test", "p")) %>% 
  pivot_wider(names_from = p, values_from = value)

lx <- read_excel("Hanlon_2023_data.xlsx", sheet = "Figure_1C_4")

cor_all_xrf <- lx %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element) %>% 
  dplyr::summarize(cos = stats::cor.test(value, root_value)$estimate,
                   p_all = stats::cor.test(value, root_value)$p.value) %>% 
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  rename(r2_all = r2)

cor_leaves_xrf <- lx %>%
  pivot_longer(cols = total_root_length:d95, names_to = 'root_trait', values_to = "root_value") %>% 
  drop_na() %>% 
  group_by(root_trait, Element, leaf) %>% 
    dplyr::summarize(cos = stats::cor.test(value, root_value)$estimate,
                   p = stats::cor.test(value, root_value)$p.value) %>%
  mutate(r2 = round((cos)^2, 2)) %>% 
  dplyr::select(-cos) %>% 
  pivot_wider(names_from = leaf, values_from = c(r2, p))

cor_xrf <- full_join(cor_all_xrf, cor_leaves_xrf) %>% 
  mutate(method = "XRF") %>% 
  pivot_longer(cols = c(p_all:p_young), names_to = c("test"), values_to = c("value")) %>% 
  separate(test, into = c("test", "p")) %>% 
  pivot_wider(names_from = p, values_from = value)

inds <- full_join(cor_icp, cor_xrf) %>% 
  pivot_longer(all:young, names_to = "Leaves", values_to = "val") %>% 
  pivot_wider(names_from = test, values_from = val)


# just patchwork the xrf icp corrs

xic <- xio %>% 
  dplyr::select(Element, ICP, XRF) %>% 
  group_by(Element) %>% 
  dplyr::summarize(cos = stats::cor.test(ICP, XRF)$estimate,
                   p = stats::cor.test(ICP, XRF)$p.value) %>%
  mutate(r2 = cos^2, 2) %>% 
  dplyr::select(Element, r2, p) %>% 
  mutate(method = "XRF_ICP_corr") %>% 
  distinct()

sup2 <- inds %>% 
  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 = method, y = Element, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(name = "r2", limits = c(0,1)) +
  geom_text(aes(label = round(r3, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  theme_mh() +
  facet_wrap(.~root_trait+Leaves)

ggsave("Supplemental Figure 2.png", width = 10, height = 30, units = "in")

# only plot ear leaf
inds <- inds %>% 
  filter(root_trait %in% c("total_root_length", "shallow_20", "deep_30", "deep_40", "d50", "d70", "d95"))

inds$root_trait <- factor(inds$root_trait, 
                          levels = c("total_root_length", "shallow_20", "deep_30", "deep_40", "d50", "d70", "d95"),
                          labels = c("Total", "<20cm",">30cm", ">40cm", "d50", "d70", "d95"))

cp <- inds %>% 
  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)) %>% 
  filter(Leaves == "ear") %>% 
  ggplot(aes(x = method, y = Element, 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") +
  theme_mh() +
  xlab(NULL) +
  facet_wrap(.~root_trait, ncol = 10) +
  theme(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.text.y = element_blank(),
        axis.title.y = element_blank(),
        axis.ticks.y = element_blank(),
        plot.margin = margin(l = -10),
        legend.key.height = unit(2, "mm"),
        legend.title = element_text(margin = margin(b = 10, t = -2))) 



xicp <- xic %>% 
  left_join((inds %>% ungroup() %>% dplyr::select(Element) %>% distinct), .) %>% 
  mutate(method = "XRF\nICP") %>% 
  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)) %>% 
  left_join((inds %>% ungroup() %>% dplyr::select(Element) %>% distinct), .) %>% 
  mutate(method = "XRF\nICP") %>% 
  ggplot(aes(x = method, y = Element, fill = r2)) +
  geom_tile(color = "white", linewidth = 0.1) +
  scale_fill_viridis(name = "r2", limits = c(0,1), expand = c(0,0), na.value = "white") +
  geom_text(aes(label = round(r3, 2), color = sig),  size = 3) +
  scale_color_manual(values = c("black", "white", NA), guide = "none") +
  theme_mh() +
  xlab(NULL) +
  facet_grid(.~method) +
  theme(legend.position = "none",
        panel.spacing.x = unit(0, "lines"),
        legend.margin = margin(0, 0, 0, 0),
        legend.box.margin = margin(b = -5),
        axis.ticks.x = element_blank(),
        axis.text.x = element_blank(),
        plot.margin = margin(r = -10))

xicp + cp +
  plot_layout(widths = c(0.5, 6))

ggsave("Figure1C.png", width = 8, height = 4, units = "in")
```

