## Griffiths M & York LM (2019) Targeting root nutrient uptake kinetics for increasing plant productivity

#################
## User setup  ##
#################
# 1) Set R working directory to folder with kinetics meta-analysis csv data
# 2) Install or load following packages
library(tidyverse) #ggplot2, purrr, tibble, dplyr, tidyr, stringr, readr, forcats

##############################
## Metaanalysis Regressions ##
##############################
dat <- read_csv(file=paste("Kinetics_metaanalysis_figure_subset.csv"), na = c("NA", "na", "n.a.", ""))

dat <- mutate(dat,
              LUR.HUR = (LUR / HUR)
			  )

dat <- mutate(dat,
              HUR.LUR = (HUR / LUR)
			 )

dat <- mutate(dat,
              Imax.Km = (Imax / Km)
			  )

dat$Geno <- as.factor(dat$Geno)

dat <- dat %>% filter(Nutrient=="Nitrate")

###########################
## Metaanalysis Boxplots ##
###########################
tempdf <- dat %>% group_by(Geno) %>% summarise(avgdep = median(Imax, na.rm = TRUE))
tempdf <- tempdf %>% arrange(avgdep)
dat$Geno <- factor(dat$Geno, levels = unique(tempdf$Geno), ordered = TRUE);
ggplot(dat, aes(x=Geno, y=Imax)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(shape=Geno), position=position_jitter(0.05)) +
  theme_bw() + 
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,axis.text.x = element_text(angle=45,vjust=0.5)
  ) + 
  xlab(bquote("")) +
  ylab(bquote("Imax (µmol" ~ g^-1 ~ h^-1*")"))
ggsave(file=paste("Nitrate_Imax_boxplot.png"), width=7, height=3, dpi=300)
ggsave(file=paste("Nitrate_Imax_boxplot.pdf"), width=7, height=3, dpi=300, useDingbats=FALSE)

ggplot(dat, aes(x=Geno, y=Km)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(shape=Geno), position=position_jitter(0.05)) +
  theme_bw() + 
  theme(
    plot.background = element_blank()
    ,panel.grid.major = element_blank()
    ,panel.grid.minor = element_blank()
    ,axis.text.x = element_text(angle=45,vjust=0.5)
  ) + 
  xlab(bquote("")) +
  ylab(bquote("Km (µM)"))
ggsave(file=paste("Nitrate_Km_boxplot.png"), width=7, height=3, dpi=300)
ggsave(file=paste("Nitrate_Km_boxplot.pdf"), width=7, height=3, dpi=300, useDingbats=FALSE)

#################
## Regressions ##
#################
## part 1 define model
ggplotRegression  <- function(dat, xvar, yvar){
  fml <- paste(yvar, "~", xvar)
  fit <- lm(fml, dat)
  ggplot(fit$model, aes_string(x = names(fit$model)[2], y = names(fit$model)[1])) + 
    geom_point(data = dat, aes(shape=Geno)) +
    stat_smooth(method = "lm", col = "red") +
    theme_bw() + 
    theme(
      plot.background = element_blank()
      ,panel.grid.major = element_blank()
      ,panel.grid.minor = element_blank()
    ) +
    
    geom_text(aes(x = 6, y = max(dat$Imax*0.95)), hjust = 0, 
              label = paste("R2 = ",signif(summary(fit)$r.squared, 5),
                            " Intercept =",signif(fit$coef[[1]],5)
              )) +
    geom_text(aes(x = 6, y = max(dat$Imax*0.90)), hjust = 0, 
              label = paste("Slope =",signif(fit$coef[[2]], 5),
                            " P =",signif(summary(fit)$coef[2,4], 5)
              ))
}

## part 2 plot graph
p <- ggplotRegression(dat, "HUR", "Imax")
p + scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  xlab(bquote("Highest uptake rate reported (µmol" ~ g^-1 ~ h^-1*")")) +
  ylab(bquote("Imax (µmol" ~ g^-1 ~ h^-1*")"))
ggsave(file=paste("Imax_vs_HUR_all.png"), width=4.75, height=2.75, dpi=300)
ggsave(file=paste("Imax_vs_HUR_all.pdf"), width=4.75, height=2.75, dpi=300, useDingbats=FALSE)

p <- ggplotRegression(dat, "Km", "Imax")
p + scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  xlab(bquote("Km (µM)")) + 
  ylab(bquote("Imax (µmol" ~ g^-1 ~ h^-1*")"))
ggsave(file=paste("Imax_vs_Km_all.png"), width=4.75, height=2.75, dpi=300)
ggsave(file=paste("Imax_vs_Km_all.pdf"), width=4.75, height=2.75, dpi=300, useDingbats=FALSE)

filtereddata <- dat %>% filter(Range_high<5000)
p <- ggplotRegression(filtereddata, "Range_high", "Imax")
p + scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  xlab(bquote("Highest concentration tested (µM)")) +
  ylab(bquote("Imax (µmol" ~ g^-1 ~ h^-1*")"))
ggsave(file=paste("Imax_vs_conc_all.png"), width=4.75, height=2.75, dpi=300)
ggsave(file=paste("Imax_vs_conc_all.pdf"), width=4.75, height=2.75, dpi=300, useDingbats=FALSE)

p <- ggplotRegression(dat, "at.50uM", "Imax")
p + scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0, 0)) +
  xlab(bquote("Uptake rate at 50µM (µmol" ~ g^-1 ~ h^-1*")")) +
  ylab(bquote("Imax (µmol" ~ g^-1 ~ h^-1*")"))
ggsave(file=paste("Imax_vs_50uM_all.png"), width=4.75, height=2.75, dpi=300)
ggsave(file=paste("Imax_vs_50uM_all.pdf"), width=4.75, height=2.75, dpi=300, useDingbats=FALSE)

############################
##       Pie chart        ##
## Species/nutrient count ##
############################
dat <- read_csv(file=paste("Kinetics_metaanalysis_csv.csv"), na = c("NA", "na", "n.a.", ""))

uval <- unique(dat[c("Reference", "Species", "Nutrient")])
uval <- uval %>% group_by(Species, Nutrient) %>% summarise(Count = n())
uval <- uval %>% arrange(desc(Species)) %>% mutate(lab.ypos = cumsum(Count) - 0.5*Count)
uval

mycols <- c('#bababa'
,'#d73027'
,'#fc8d59'
,'#2ca25f'
,'#1c9099'
,'#91bfdb'
,'#4575b4'
)

lvl0 <- tibble(Species = "Parent", Count = 0, level = 0, fill = NA)
lvl1 <- uval %>%
  group_by(Species) %>%
  summarise(Count = sum(Count)) %>%
  ungroup() %>%
  mutate(level = 1) %>%
  mutate(fill = Species)
lvl2 <- uval %>%
  select(Species = Nutrient, Count, fill = Species) %>%
  mutate(level = 2)

bind_rows(lvl0, lvl1, lvl2) %>%
  mutate(Species = as.factor(Species) %>% fct_reorder2(fill, Count)) %>%
  arrange(fill, Species) %>%
  mutate(level = as.factor(level)) %>%
  ggplot(aes(x = level, y = Count, fill = fill, alpha = level)) +
  geom_col(width = 1, color = "white", position = position_stack()) +
  #geom_col(width = 1, color = "white", size = 0.25, position = position_stack()) +
  geom_text(aes(label = Species), position = position_stack(vjust = 0.5), color = "white") +
  geom_text(aes(label = Count), position = position_stack(vjust = 0.6), color = "white") +
  coord_polar(theta = "y") +
  scale_alpha_manual(values = c("0" = 0, "1" = 1, "2" = 0.7), guide = F) +
  scale_x_discrete(breaks = NULL) +
  scale_y_continuous(breaks = NULL) +
  scale_fill_manual(values = mycols, na.translate = F) +
  labs(x = NULL, y = NULL) +
  labs(fill='Species') +
  theme_void()

ggsave(file=paste("Metaanalysis_species_nutrient_pies.png"), width=8, height=8, dpi=300)
ggsave(file=paste("Metaanalysis_species_nutrient_pies.pdf"), width=8, height=8, dpi=300)

#####################
##    Pie chart    ##
## Plant age count ##
#####################
dat <- read_csv(file=paste("Kinetics_metaanalysis_csv.csv"), na = c("NA", "na", "n.a.", ""))

uval <- unique(dat[c("Reference", "Plant_age_d")])
uval <- mutate(uval, calc = (Plant_age_d-1))
uval <- mutate(uval, bin = floor(calc/7))
uval <- uval %>% group_by(bin) %>% summarise(Count = n())
uval <- na.omit(uval)
uval <- uval %>% arrange(desc(bin)) %>% mutate(lab.ypos = cumsum(Count) - 0.5*Count)
uval$bin <- as.factor(uval$bin)

mycols <- c('#9e0142',
            '#d53e4f',
            '#f46d43',
            '#fdae61',
            '#fee08b',
            '#e6f598',
            '#abdda4',
            '#66c2a5',
            '#3288bd')

ggplot(uval, aes(x = 2, y = Count, fill = bin)) +
  geom_bar(stat = "identity", color = "white") +
  coord_polar(theta = "y", start = 0) + 
  geom_text(aes(y = lab.ypos, label = Count), color = "white")+
  scale_fill_manual(values = mycols) +
  theme_void() +
  xlim(0.5, 2.5) +
  labs(fill='Plant age') 

ggsave(file=paste("Metaanalysis_plantage.png"), width=4.75, height=3.1968, dpi=300)
ggsave(file=paste("Metaanalysis_plantage.pdf"), width=4.75, height=3.1968, dpi=300)

########################
##      Pie chart     ##
## Pretreatment count ##
########################
dat <- read_csv(file=paste("Kinetics_metaanalysis_csv.csv"), na = c("NA", "na", "n.a.", ""))

uval <- unique(dat[c("Reference", "Plant_pretreatment")])
uval <- uval %>% group_by(Plant_pretreatment) %>% summarise(Count = n())
uval <- uval %>% arrange(desc(Plant_pretreatment)) %>% mutate(lab.ypos = cumsum(Count) - 0.5*Count)
uval

mycols <- c('#d73027'
,'#fc8d59'
,'#fee090'
,'#bababa'
)

ggplot(uval, aes(x = 2, y = Count, fill = Plant_pretreatment)) +
  geom_bar(stat = "identity", color = "white") +
  coord_polar(theta = "y", start = 0) + 
  geom_text(aes(y = lab.ypos, label = Count), color = "white")+
  scale_fill_manual(values = mycols) +
  theme_void() +
  xlim(0.5, 2.5) +
  labs(fill='Plant pretreatment') 

ggsave(file=paste("Metaanalysis_pretreatment_pie.png"), width=4.75, height=3.1968, dpi=300)
ggsave(file=paste("Metaanalysis_pretreatment_pie.pdf"), width=4.75, height=3.1968, dpi=300)


#####################
##    Pie chart    ##
## Root type count ##
#####################
dat <- read_csv(file=paste("Kinetics_metaanalysis_csv.csv"), na = c("NA", "na", "n.a.", ""))

uval <- unique(dat[c("Reference", "Root_sample_type")])
uval <- uval %>% group_by(Root_sample_type) %>% summarise(Count = n())
uval <- uval %>% arrange(desc(Root_sample_type)) %>% mutate(lab.ypos = cumsum(Count) - 0.5*Count)
uval

mycols <- c('#bababa'
,'#4575b4'
)

ggplot(uval, aes(x = 2, y = Count, fill = Root_sample_type)) +
  geom_bar(stat = "identity", color = "white") +
  coord_polar(theta = "y", start = 0) + 
  geom_text(aes(y = lab.ypos, label = Count), color = "white")+
  scale_fill_manual(values = mycols) +
  theme_void() +
  xlim(0.5, 2.5) +
  labs(fill='Root sample')

ggsave(file=paste("Metaanalysis_roottype.png"), width=4.75, height=3.1968, dpi=300)
ggsave(file=paste("Metaanalysis_roottype.pdf"), width=4.75, height=3.1968, dpi=300)