code.qmd

---
title: "Harvest index is a key trait for screening drought-tolerant potato genotypes (_Solanum tuberosum_)"
subtitle: "Journal of Crop Science and Biotechnology"
format:
  html:
    toc: true
    toc-location: left
    number-sections: true
    self-contained: true
    output-file: "ESM_1"
editor_options: 
  chunk_output_type: console
execute: 
  echo: true
  message: false
  warning: false
---

# Setup

```{r}
#| label: setup

source("https://inkaverse.com/setup.r")
library(gtsummary)
library(factoextra)
library(corrplot)
library(emmeans)
library(multcomp)

devtools::session_info()
```

# Googlesheets connect

```{r}
url <- "https://docs.google.com/spreadsheets/d/1dfgpmCKdPmxRHozrZp0iE_xMGsvKTIcztDpMWYSEGaY"
# browseURL(url)
gs <- as_sheets_id(url)
```

# Trait rename

```{r}
geno <- c("CIP720088" = "G01"
          , "CIP392797.22" = "G02"
          , "CIP397077.16" = "G03"
          , "CIP398192.213" = "G04"
          , "CIP398180.612" = "G05"
          , "CIP398208.704" = "G06"
          , "CIP398098.119" = "G07"
          , "CIP398190.89" = "G08"
          , "CIP398192.592" = "G09"
          , "CIP398201.510" = "G10"
          , "CIP398203.244" = "G11"
          , "CIP398203.5" = "G12"
          , "CIP398208.219" = "G13"
          , "CIP398208.33" = "G14"
          , "CIP398208.620" = "G15")

traits <- c("SPAD_29" = "Chlorophyll concentration (SPAD) at 29 dap"
            , "SPAD_59" = "Chlorophyll concentration (SPAD) at 59 dap"
            , "SPAD_76" = "Chlorophyll concentration (SPAD) at 76 dap"
            , "SPAD_83" = "Chlorophyll concentration (SPAD) at 83 dap"
            , "HGT" = "Plant height (cm)"
            , "RWC" = "Relative water content (%)"
            , "LOP" = "Leaf osmotic potential (MPa)"
            , "LDW" = "Leaf dry weight (g)"
            , "SDW" = "Stem dry weight (g)"
            , "RDW" = "Root dry weight (g)"
            , "TDW" = "Tuber dry weight (g)"
            , "NTUB" = "Tuber number (N°)" 
            , "TRS" = "Total transpiration (mL)"
            , "LFA" = "Leaf area (cm^2^)"
            , "RTL" = "Root length (cm)"
            , "TDB" = "Total dry biomass (g)"
            , "HI"  = "Harvest index (HI)" 
            , "SLA" = "Specific leaf area (cm^2^ g^-1^)"
            , "RCC" = "Relative chlorophyll content (RCC)"
            , "WUE_B" = "Biomass water use efficiency (g L^-1^)"
            , "WUE_T" = "Tuber water use efficiency (g L^-1^)"
            )
```

# Import data

```{r}
fb <- gs %>% 
  range_read("fb") %>% 
  rename_with("tolower") %>%
  rename_with(~gsub("\\s+|\\.", "_", .)) %>% 
  mutate(treat = case_when(treat %in% "wellwater" ~ "WW"
                           , treat %in% "drydown" ~ "WD")) %>% 
  dplyr::select(block, treat, genotype,
         spad_29 = spad_29dap, 
         spad_59 = spad_59dap, 
         spad_76 = spad_76dap, 
         spad_83 = spad_83dap,
         hgt = hgt_86dap,
         rwc = rwc_84dap,
         lop = op_84dap,
         ldw = leafdw,
         sdw = stemdw,
         rdw = rootdw,
         tdw = tubdw,
         ntub,
         trs = ttrns,
         lfa = la,
         rtl = rlg
         ) %>% 
  mutate(tdb = (ldw+sdw+rdw+tdw),
         hi = tdw/(ldw+sdw+rdw+tdw),
         sla = lfa/ldw,
         rcc = (spad_83/lfa),
         wue_b = (ldw+sdw+rdw+tdw)/trs,
         wue_t = tdw/trs
         ) %>% 
  mutate(gnt = recode(genotype, !!!geno), .after = genotype) %>% 
  dplyr::select(genotype, gnt, treat, block, everything()) %>% 
  mutate(across(where(is.character), as.factor)) %>% 
  rename_with(., toupper, where(is.numeric)) 

# str(fb)

fb %>% web_table()
```

# Abbreviations

```{r}
gs %>% 
  range_read("var") %>% 
  filter(include == "x") %>%
  dplyr::select(Variable, Abbreviation) %>%
  kable()
```

# Table 1

```{r}
tab <- gs %>% 
  range_read("gnt") %>% 
  dplyr::select(Number, Genotypes, Adaptability, "Growing period" = GPL, "Heat tolerance" = Heat, "Dry matter (%)")
  
# tab %>% sheet_write(gs, data = .,  sheet = "tab1")

tab %>% web_table()
```

# Table 2

```{r}
tab <- fb %>%
  dplyr::select(!c(block, gnt, genotype)) %>% 
  tbl_summary(
    by = treat, 
    statistic = list(all_continuous() ~ "{mean} ± {sd}"),
    missing = "no") %>% 
  add_p() %>% 
  bold_labels() %>% 
  as_tibble() %>% 
  mutate(`**Characteristic**` = recode(
    `**Characteristic**`
    , "__SPAD_29__" = "Chlorophyll concentration (SPAD) at 29 dap"
    , "__SPAD_59__" = "Chlorophyll concentration (SPAD) at 59 dap"
    , "__SPAD_76__" = "Chlorophyll concentration (SPAD) at 76 dap"
    , "__SPAD_83__" = "Chlorophyll concentration (SPAD) at 83 dap"
    , "__HGT__" = "Plant height (cm)"
    , "__RWC__" = "Relative water content (%)"
    , "__LOP__" = "Leaf osmotic potential (MPa)"
    , "__LDW__" = "Leaf dry weight (g)"
    , "__SDW__" = "Stem dry weight (g)"
    , "__RDW__" = "Root dry weight (g)"
    , "__TDW__" = "Tuber dry weight (g)"
    , "__NTUB__" = "Tuber number (N°)" 
    , "__TRS__" = "Total transpiration (mL)"
    , "__LFA__" = "Leaf area (cm^2^)"
    , "__RTL__" = "Root length (cm)"
    , "__TDB__" = "Total dry biomass (g)"
    , "__HI__"  = "Harvest index (HI)" 
    , "__SLA__" = "Specific leaf area (cm^2^ g^-1^)"
    , "__RCC__" = "Relative chlorophyll content (RCC)"
    , "__WUE_B__" = "Biomass water use efficiency (g L^-1^)"
    , "__WUE_T__" = "Tuber water use efficiency (g L^-1^)"
    )) %>% 
  rename(Variable = `**Characteristic**`
         , `Water deficit` = `**WD**, N = 75`
         , `Well-Watered`= `**WW**, N = 75`
         , `p-value` = `**p-value**`) 

# tab %>% sheet_write(gs, data = .,  sheet = "tab2")

tab %>% web_table()
```

# Figure 1

```{r}
fts <- gs %>% 
  range_read("ftsw") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = fts, -ID, -Genotype, -Treatment)

model <- fts %>% 
  inti::mean_comparison(response = "fts"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt1 <- model$comparison %>% 
  mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(
  type = "line", color = F
  , x = "day"
  , y = "fts"
  , group = "Treatment"
  , ylab = "Fraction of transpirable soil water"
  , xlab =  "Days"
  , glab = "Treatment"
  , legend = "none"
  , ylimits = c(0, 1.08, 0.1)
  , error =  "ste"
  ) +
  scale_x_continuous(breaks = c(0:100)*2
                     , limits = c(34, 88))

trns <- gs %>% 
  range_read("trans") %>%
  filter(Treatment != "preharvest") %>%
  tidyr::gather(key = day, value = trans, -ID, -Genotype, -Treatment) %>%
  filter(day != "TOTAL") %>%
  drop_na()


model <- trns %>% 
  inti::mean_comparison(response = "trans"
                        , model_factors = "Treatment*day"
                        , c("Treatment", "day"))

plt2 <- model$comparison %>% 
    mutate(Treatment = case_when(
    Treatment == "drydown" ~ "Water Deficit (WD)"
    , Treatment == "wellwater" ~ "Well Watered (WW)"
  )) %>% 
  mutate(across(day, as.numeric)) %>% 
  plot_smr(type = "line", color = F
            , x = "day"
            , y = "trans"
            , group = "Treatment"
            , ylab = "Transpiration (mL)"
            , xlab =  "Days"
            , glab = "Irrigation"
            , ylimits = c(0, 700, 100)
            , error = "ste"
            ) +
  scale_x_continuous(breaks = c(0:100)*2
                     , limits = c(34, 88)
                     ) +
  theme(axis.title.x=element_blank(),
        axis.text.x=element_blank(),
        axis.ticks.x=element_blank())

plot <- ggdraw(xlim = c(0, 0.5), ylim = c(0,0.9)) +
  draw_plot(plt2, width = 0.5, height = 0.43 , x = 0.0, y = 0.47 ) +
  draw_plot(plt1, width = 0.496, height = 0.5, x = 0.004, y = 0.0 ) +
          draw_plot_label(
            label = c("a", "b")
            , size = c(16, 15)
            , x = c(0.465, 0.465)
            , y = c(0.79, 0.49)
            )

plot %>% 
  ggsave2(plot = .
          , "files/Fig1.pdf", width = 18, height = 15, units = "cm")

"files/Fig1.pdf" %>% include_graphics()

plot %>% 
  ggsave2(plot = .
          , "files/Fig1.eps", width = 18, height = 15, units = "cm")
```

# Figure 2

```{r}
# tdw

model <- fb %>% 
  lm(TDW ~ block + gnt*treat, data = .)
  
anova(model)

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcTDW <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
  mutate(across(emmean, ~round(., 2))) %>% 
  unite(TDW, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

tdw <- mc %>% 
  plot_smr(type = "line"
          , color = F
          , x = "gnt"
          , xlab = ""
          , y = "emmean"
          , ylab = "Tuber dry weight (g)"
          , group = "treat"
          , glab = "Treatment"
          , legend = "none"
          , error = "SE"
          , ylimits = c(0, 100, 20)
          , sig = ".group"
          )

# rcc

model <- fb %>% 
  lm(RCC ~ block + gnt*treat, data = .)
  
anova(model)

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcRCC <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
  mutate(across(emmean, ~round(., 2))) %>% 
  unite(RCC, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

rcc <- mc %>% 
  plot_smr(type = "line"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "Relative chlorophyll content"
        , ylimits = c(0, 0.1, 0.02)
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , error = "SE"
        , sig = ".group"
        ) 

# hi

model <- fb %>% 
  lm(HI ~ block + gnt*treat, data = .)
  
anova(model)

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcHI <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
  mutate(across(emmean, ~round(., 2))) %>% 
  unite(HI, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

hi <- mc %>% 
  plot_smr(type = "line"
        , color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "Harvest Index"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 1, 0.2)
        , error = "SE"
        , sig = ".group"
        ) 

# wue_t

model <- fb %>% 
  lm(WUE_T ~ block + gnt*treat, data = .)
  
anova(model)

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcWUE_T <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
    mutate(across(emmean, ~round(., 2))) %>% 
  unite(WUE_T, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

wue_t <- mc %>% 
  plot_smr(type = "line", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "Tuber water use efficiency (g/L)"
        , group = "treat"
        , glab = "Treatment"
        , legend = "none"
        , ylimits = c(0, 10, 2)
        , error = "SE"
        , sig = ".group"
        ) 

# spad_29

model <- fb %>% 
  lm(SPAD_29 ~ block + gnt*treat, data = .)
  
anova(model)

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcSPAD_29 <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
    mutate(across(emmean, ~round(., 2))) %>% 
  unite(SPAD_29, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

spad_29 <- mc %>% 
    mutate(treat  = case_when(
    treat  == "WD" ~ "Water Deficit (WD)"
    , treat  == "WW" ~ "Well Watered (WW)"
    )) %>%
  plot_smr(type = "line"
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "SPAD at 29 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend = "top"
        , ylimits = c(30, 70, 5)
        , error = "SE"
        , color = F
        , sig = ".group"
        )

# spad_83

model <- fb %>% 
  lm(SPAD_83 ~ block + gnt*treat, data = .)
  
anova(model)

mc <- emmeans(model, ~ treat | gnt) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws))

mcSPAD_83 <- emmeans(model, ~ gnt | treat) %>% 
  cld(Letters = letters, reversed = T) %>% 
  mutate(across(.group, trimws)) %>% 
    mutate(across(emmean, ~round(., 2))) %>% 
  unite(SPAD_83, c(emmean, .group), sep = " ") %>% 
  dplyr::select(1:3)

spad_83 <- mc %>% 
  plot_smr(type = "line", color = F
        , x = "gnt"
        , xlab =  ""
        , y = "emmean"
        , ylab = "SPAD at 83 dap"
        , group = "treat"
        , glab = "Treatment"
        , legend= "none"
        , ylimits = c(30, 70, 5)
        , error = "SE"
        , sig = ".group"
        ) 

plot <- ggdraw(xlim = c(0.0, 1), ylim = c(0.0, 1)) +
  draw_plot(spad_29 , width = 0.49, height = 0.35, x = 0.0, y = 0.6) +
  draw_plot(spad_83, width = 0.49, height = 0.3, x = 0.5, y = 0.6) +
  draw_plot(rcc,  width = 0.49, height = 0.3, x = 0.0, y = 0.3) +
  draw_plot(tdw, width = 0.49, height = 0.3, x = 0.5, y = 0.3) +
  draw_plot(hi , width = 0.49, height = 0.3, x = 0.0, y = 0.0) +
  draw_plot(wue_t, width = 0.49, height = 0.3, x = 0.5, y = 0.0) +
          draw_plot_label(
            label = c("a", "b", "c", "d", "e", "f"),
            x = c(0.06, 0.56, 0.06, 0.56, 0.06, 0.56),
            y = c(0.89, 0.89, 0.59, 0.59, 0.29, 0.29))

plot %>% 
  ggsave(plot = .
         , "files/Fig2.pdf"
         , units = "cm", width = 28, height = 25)

"files/Fig2.pdf" %>% include_graphics()

plot %>% 
  ggsave(plot = .
         , "files/Fig2.eps"
         , units = "cm", width = 28, height = 25)


geno4treat <- merge(mcSPAD_29, mcSPAD_83) %>% 
  merge(., mcRCC) %>% 
  merge(., mcTDW) %>%
  merge(., mcHI) %>%
  merge(., mcWUE_T) %>% 
  arrange(desc(treat)) %>% 
  rename(Genotype = gnt, Treatment = treat)
  
# geno4treat %>% sheet_write(gs, data = .,  sheet = "tabS1")
```

# Rename variables in PCA

```{r}

vars <- names(fb)

model <- 5:length(vars) %>% map(function(x) {
  
    fb %>% 
      H2cal(trait = vars[x]
            , gen.name = "genotype"
            , rep.n = 4
            , fixed.model = "0 + (1|block) + genotype"
            , random.model = "1 + (1|block) + (1|genotype)"
            , emmeans = T
            , plot_diag = T
            , outliers.rm = T
            )
})

tabsmr <- 1:length(model) %>% map(function(x) {
      model[[x]][["tabsmr"]] 
      }) %>% 
  bind_rows()

varnames <- tabsmr %>% 
  dplyr::select(trait, repeatability) %>% 
  mutate(name = paste0("(",round(repeatability, 2), ") ", trait)) %>% 
  dplyr::select(!repeatability) %>% 
  deframe()
```

# Figure 3

```{r potatos}
fig <- magick::image_read("files/potatos.png") %>% 
  grid::rasterGrob() %>%
  ggsave("files/Fig3.pdf", plot =  ., width = 13, height = 13)

fig <- magick::image_read("files/potatos.png") %>% 
  grid::rasterGrob() %>%
  ggsave("files/Fig3.eps", plot =  ., width = 13, height = 13)

"files/Fig3.pdf" %>% include_graphics()
```

# Figure 4

```{r}
library(gridGraphics)

mvdata <- fb %>% 
  group_by(genotype, treat) %>% 
  summarise(across(where(is.numeric), ~mean(.x, na.rm = T))) %>% 
  ungroup() %>% 
  unite(trait, c(genotype, treat), sep = "-") 

plot <- ~{
  
  mvdata %>% 
    dplyr::select(where(is.numeric)) %>% 
    dplyr::select(!ends_with(c("29", "59", "76"))) %>% 
    pairs.panels(x = .
                 , hist.col="red"
                 , pch = 21
                 , stars = TRUE
                 , scale = FALSE
                 , lm = TRUE
                 , cex.axis = 0.6
                 ) 
  }

list(plot)  %>%  
  plot_grid(plotlist = .
            , nrow = 1
            ) %>% 
  ggsave2("files/Fig4.pdf", plot =  ., width = 30, height = 30, units = "cm")

"files/Fig4.pdf" %>% include_graphics()

list(plot)  %>%  
  plot_grid(plotlist = .
            , nrow = 1
            ) %>% 
  ggsave("files/Fig4.eps", plot =  ., width = 30, height = 30, units = "cm")
```

# Figure 5

```{r}
mv <- mvdata %>% 
  pivot_longer(!trait) %>% 
  mutate(name = recode(name, !!!varnames)) %>% 
  pivot_wider() %>% 
  column_to_rownames("trait") %>% 
  PCA(scale.unit = T, graph = F)

hcpc <- HCPC(mv, graph = F)

tree <- ~ {
  
  par(cex=0.7)
  
  plot.HCPC(x = hcpc
          , choice = "map"
          , legend = list(x = "topright"
                          , cex = 0.7
                          , inset = 0.001
                          , box.lty=0
                          )
          , draw.tree = F
          )
}

pca <- mv %>% 
  plot.PCA(choix = "var", autoLab = "auto", cex = 0.6)

plot <- plot_grid(pca, tree
                  , nrow = 1, labels = "auto")
plot %>% 
  ggsave2(plot = .
          , "files/Fig5.pdf"
          , units = "cm", width = 30, height = 15)

"files/Fig5.pdf" %>% include_graphics()

plot %>% 
  ggsave2(plot = .
          , "files/Fig5.eps"
          , units = "cm", width = 30, height = 15)
```

# Supplementary information

# Supplementary Table 1

```{r}
h2plot <- tabsmr %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  dplyr::select(trait, mean, std, min, max, V.g, V.e, repeatability) %>% 
  mutate(across(where(is.numeric), ~ round(.x, digits = 2))) %>% 
  mutate(trait = recode(trait, !!!traits)) %>% 
  rename(Trait = trait, H2 = repeatability)
  
h2plot  %>% web_table()

# h2plot %>% sheet_write(gs, data = .,  sheet = "tabS1")
```

# Supplementary Figure 1

```{r}
var <- get_pca_var(mv)

tmp <- tempfile(fileext = ".png")
ppi <- 300
png(tmp, width=3.8*ppi, height=10*ppi, res=ppi)
corrplot(var$cor, 
         method="number",
         tl.col="black", 
         tl.srt=45,)
graphics.off()

pt1 <- png::readPNG(tmp) %>%
  grid::rasterGrob(interpolate = TRUE)

pt2 <- fviz_eig(mv, 
                addlabels=TRUE,
                hjust = 0.05,
                barfill="white",
                barcolor ="darkblue",
                linecolor ="red") + 
  ylim(0, 50) + 
  labs(
    title = "PCA - percentage of explained variances",
    y = "Variance (%)") +
  theme_minimal()

pt3 <- fviz_contrib(mv,
                     choice = "var", 
                     axes = 1, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 1 - variables contribution") 

pt4 <- fviz_contrib(mv,
                     choice = "var", 
                     axes = 2, 
                     top = 10,
                     fill="white",
                     color ="darkblue",
                     sort.val = "desc") +
  ylim(0, 18) + 
  labs(title = "Dim 2 - variables contribution") 

plot <- ggdraw(xlim = c(0.0, 1.0), ylim = c(0, 1.0))+
  draw_plot(pt1,  width = 0.4, height = 0.99, x = 0.6, y = 0.0) +  
  draw_plot(pt2,  width = 0.6, height = 0.34, x = 0.03, y = 0.66) +
  draw_plot(pt3, width = 0.6, height = 0.34, x = 0.03, y = 0.33) + 
  draw_plot(pt4, width = 0.6, height = 0.34, x = 0.03, y = 0.0) +
        draw_plot_label(
    label = c("a", "b", "c", "d"),
    x = c(0.005, 0.005, 0.005, 0.65),
    y = c(0.999, 0.67, 0.34, 0.999))

ggsave2(plot = plot, "files/FigS1.pdf", height = 25, width = 30, units = "cm")

"files/FigS1.pdf" %>% include_graphics()
```

# Convert figures to pdf to jpg

```{r}
library(magick)

figs <- list.files(path = "files/", pattern = "eps", full.names = T)

convert <- 1:length(figs) %>% map(\(x) {
  
  figname <- figs[x] %>% basename() %>% gsub("eps", "jpeg", .)
  
  path2xport <- dirname(figs[x])
  
  figs[x] %>% 
    image_read() %>% 
    image_write(file.path(path2xport, figname))
})
```