analysis.Rmd

---
title: "Analysis and visualization"
date: "Last knitted on `r format(Sys.Date(), '%d %b %Y')`"
author: "Toby Halamka, Sebastian Kopf"
output:
  html_document: 
    df_print: paged
    number_sections: no
    toc: yes
    toc_float: true
    toc_depth: 3
    code_folding: show
editor_options:
  chunk_output_type: console
---

# Libraries & Scripts

```{r setup, include = TRUE, message=FALSE}
# load libraries
library(tidyverse) # dplyr, tidyr, ggplot2
library(cowplot) # multi panel plots
library(ggrepel) # plot annotation placement
library(latex2exp) # formatting of plot labels
library(readxl) # read excel data
library(openxlsx) # write tables in excel format

# load scripts
source(file.path("scripts", "growth_functions.R"))
source(file.path("scripts", "plotting_functions.R"))
source(file.path("scripts", "table_functions.R"))

# Global Knitting options
knitr::opts_chunk$set(
  echo = TRUE, # switch to FALSE to avoid code blocks from displaying
  dev = c("png", "pdf" , "postscript"), fig.keep = "all",
  dev.args = list(pdf = list(encoding = "WinAnsi", useDingbats = FALSE)),
  fig.path = file.path("figures", "")
)

# output folders
if (!dir.exists("figures")) dir.create("figures")
if (!dir.exists("tables")) dir.create("tables")
```

# Data Loading & Calculations

## Metadata

```{r}
experiments <- read_excel(file.path("data", "metadata.xlsx")) %>%
  # units adjustment for paper guidelines (space before all units)
  mutate(
    `C source` = str_replace_all(`C source`, "g/L", " g/L")
  )
```

## Growth rates

```{r}
# load data
growth_curves <- read_excel(file.path("data", "growth_data.xlsx"))

# convert to tidy format
growth_curves_tidy <- growth_curves %>%
  # list what to include
  pivot_longer(matches("^rep\\d"), names_to = "replicate", values_to = "OD") %>% 
  # mark death phases
  mark_death_phase(time = time.hours, N = OD, group_by = c(bug_ID, exp_ID, replicate))

# calculate growth rates
growth_rates <- 
  growth_curves_tidy %>%
  filter(!death_phase) %>% 
  estimate_growth_curve_parameters(
    time = time.hours, 
    N = OD,
    group_by = c(bug_ID, exp_ID, replicate)
  ) %>%
  # growth rates from 1/hr to 1/d
  mutate(r.1_d = r * 24) %>%
  # add metadata
  left_join(experiments, by = c("bug_ID", "exp_ID"))
```

## Lipid abundances

```{r}
# load data
fames <- read_excel(file.path("data", "fames_data.xlsx"))
tetraethers <- read_excel(file.path("data", "tetraether_data.xlsx"))

# convert to tidy format
fames_tidy <- fames %>%
  pivot_longer(-c(bug_ID, exp_ID, replicate), names_to = "compound", values_to="amount.ug") %>%
  mutate(compound = str_remove(compound, fixed(".ug")))
tetraethers_tidy <- tetraethers %>%
  pivot_longer(-c(bug_ID, exp_ID, replicate), names_to = "compound", values_to="amount.ng") %>%
  mutate(
    compound = str_remove(compound, fixed(".ng")),
    amount.ug = amount.ng / 1000
  )

# calculate percentage membrane composition
lipids <- 
  bind_rows(
    mutate(fames_tidy, category = "FAMEs"),
    mutate(tetraethers_tidy, category = "tetraethers")
  ) %>%
  group_by(bug_ID, exp_ID, replicate) %>%
  mutate(rel_amount = amount.ug/sum(amount.ug, na.rm = TRUE)) %>%
  ungroup() %>%
  # add metadata
  left_join(experiments, by = c("bug_ID", "exp_ID"))

# replicate averages
lipids_means <-
  lipids %>%
  group_by(bug_ID, exp_ID, category, compound) %>%
  summarize(
    n = sum(!is.na(rel_amount)),
    rel_amount_mean = if (n > 0) mean(rel_amount, na.rm = TRUE) else NA_real_,
    rel_amount_min = if (n > 0) rel_amount_mean - sd(rel_amount, na.rm = TRUE) else NA_real_,
    rel_amount_max = if (n > 0) rel_amount_mean + sd(rel_amount, na.rm = TRUE) else NA_real_,
    .groups = "drop"
  ) %>%
  # add metadata
  left_join(experiments, by = c("bug_ID", "exp_ID"))
```

## MS2 spectra

```{r}
# load data
ms2 <- tibble(compound = c("GDGT-1a", "GTGT-1a", "GDGT-1cisomer")) %>%
  mutate(data = map(compound, ~read_excel("data/ms2_spectra.xlsx", sheet = .x)))

# prepare ms2 data for plotting
ms2_w_peaks <- 
  ms2 %>%
  mutate(compound = as_factor(compound)) %>%
  unnest(cols = data) %>% 
  arrange(compound, Mass) %>%
  # calculate relative intensity
  group_by(compound) %>%
  mutate(rel_intensity = Intensity/max(Intensity)) %>% 
  # find peaks
  mutate(
    signal = rel_intensity > 1e-6,
    peak_nr = cumsum(c(0, diff(signal)) > 0)
  ) %>% 
  # identify peak apex intensity
  group_by(compound, peak_nr) %>%
  mutate(
    peak_rel_intensity = ifelse(peak_nr > 0, max(rel_intensity), NA_real_),
    peak_mass = Mass[rel_intensity == peak_rel_intensity[1]][1]
  ) %>%
  ungroup()
```

## Chromatograms

```{r}
# load tetraether chromatographic data
chroms <- 
  tibble(
    compound = c("GTGT-1a", "GDGT-1cisomer", "GDGT-1a"),
    mz = c(1024, 1018, 1022),
    data = map(
      compound, 
      ~readxl::read_excel("data/chroms.xlsx", sheet = .x, col_types = c("numeric", "numeric", "numeric"))
    )
  ) %>%
  unnest(cols = data)

# normalize chromatograms & make tidy
chroms_norm_tidy <-
  chroms %>%
  group_by(compound) %>%
  mutate(
    stdnormalized = Intensity_std/max(Intensity_std, na.rm = TRUE),
    samplenormalized = Intensity_sample/max(Intensity_sample, na.rm = TRUE) 
  ) %>%
  ungroup() %>%
  pivot_longer(cols = matches("normalized"), names_to = "type") %>%
  filter(!is.na(value))
```


# Visualization

## Figure 1: growth rates & lipid abundances

### Setup

```{r}
# base_plot
C_source_levels <- str_replace(unique(sort(growth_rates$`C source`)), "30 g", "\n30 g")

base_plot <- 
  ggplot() +
  aes(x = as.numeric(`C source`)) +
  geom_vline(xintercept = 1.5, linetype = 2, color = "grey50") +
  scale_x_continuous(
    limits = c(0.5, 2.5), expand = c(0, 0), 
    breaks = c(1,2), labels = function(x) C_source_levels[x]
  ) + 
  theme_figure(grid = TRUE, text_size = 20, axis_text_size = 14) +
  theme(
    panel.grid.major.x = element_blank(),
    strip.background.y = element_rect(fill = "grey50"),
    strip.text.y = element_text(color = "white"),
    axis.ticks.x = element_blank()
  ) +
  labs(x = NULL, y = NULL, fill = NULL, shape = NULL)
```


### Growth rates

```{r "fig1_growth_rates", fig.width=8, fig.height=4, warning=FALSE}
# growth rates
gr_plot_data <- 
  growth_rates %>%
  filter(bug_ID == "e.agg" & (near(`% O2`, 1) | near(`% O2`, 21))) %>%
  group_by(`% O2`, `C source`) %>%
  summarize(
    n = sum(r),
    r.1_d_mean = mean(r.1_d),
    r.1_d_sd = sd(r.1_d),
    .groups = "drop"
  ) %>%
  arrange(desc(`% O2`), `C source`) %>%
  mutate(
    `C source` = factor(str_replace(`C source`, "30 g", "\n30 g"), levels = C_source_levels),
    category = "A: growth rates",
    O2 = sprintf("$%.0f\\,%%\\,O_2$", `% O2`) %>% as_factor(),
  )

gr_plot <- 
  base_plot %+% gr_plot_data %+%
  aes(y = r.1_d_mean, shape = `C source`, ymin = r.1_d_mean - r.1_d_sd, ymax = r.1_d_mean + r.1_d_sd) + 
  geom_errorbar(width = 0) + 
  scale_y_continuous(NULL, labels = function(x) paste0(x, " / day")) +
  facet_grid(category ~ O2, labeller = latex_labeller) +
  theme(legend.key.height = unit(1, "cm")) + 
  labs(shape = NULL) +
  coord_cartesian(ylim = c(0.4, 1.0))

gr_plot + geom_point(size = 4, color = "black")
```

### Tetraether abundances

```{r "fig1_tetraether_abundances", fig.width=8, fig.height=4, warning=FALSE}
ab_plot_data <- 
  lipids_means %>%
  arrange(desc(`% O2`), desc(category), `C source`) %>%
  filter(
    near(`% O2`, 1) | near(`% O2`, 21),
    !str_detect(`C source`, "solid"),
    str_detect(compound, "br") | compound %in% c("iC15:0", "C16:0", "iDA")
  ) %>%
  mutate(
    `C source` = factor(str_replace(`C source`, "30 g", "\n30 g"), levels = C_source_levels),
    O2 = sprintf("$%.0f\\,%%\\,O_2$", `% O2`) %>% as_factor(),
    compound = factor(
      compound,
      levels = c("brGTGT_1a", "brGDGT_1c_isomer", "brGDGT_1a", "iC15:0", "C16:0", "iDA"),
      labels = c("brGTGT Ia", "brGDGT Ic isomer", "brGDGT Ia", "iso-C15:0", "C16:0", "iso-diabolic acid")
    )
  )

brgdgt_ab_plot <- 
  base_plot %+% 
  mutate(filter(ab_plot_data, category == "tetraethers"), category = "B: tetraethers") %+%
  aes(y = rel_amount_mean, fill = compound) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9, preserve = "single")) +
  geom_errorbar(
    mapping  = aes(ymin = rel_amount_min, ymax = rel_amount_max), 
    position = position_dodge(width = 0.9, preserve = "single"), width = 0, color = "black"
  ) +
  facet_grid(category ~ O2, scales = "free_y", labeller = latex_labeller, drop = FALSE) +
  scale_y_continuous(
    NULL, breaks = c(0.005, 0.015, 0.025), expand = expansion(mult = c(0, 0.03)), 
    labels = function(x) paste0(100*x, " %")) +
  scale_fill_manual(values = c("#009E73", "#0072B2", "#E69F00"))

brgdgt_ab_plot
```

### FAME abundances

```{r "fig1_FAME_abundances", fig.width=8, fig.height=4, warning=FALSE}
fames_ab_plot <- 
  base_plot %+% 
  mutate(filter(ab_plot_data, category == "FAMEs"), category = "C: FAMEs") %+%
  aes(y = rel_amount_mean, fill = compound) +
  geom_bar(stat = "identity", position = position_dodge(width = 0.9, preserve = "single")) +
  geom_errorbar(
     mapping  = aes(ymin = rel_amount_min, ymax = rel_amount_max), 
     position = position_dodge(width = 0.9, preserve = "single"), width = 0, color = "black"
  ) +
  facet_grid(category ~ O2, scales = "free_y", labeller = latex_labeller, switch = "x") +
  scale_y_continuous(
    NULL, breaks = c(0.1, 0.3, 0.5, 0.7), expand = expansion(mult = c(0, 0.03)), 
    labels = function(x) paste0(100*x, " %")) +
  scale_fill_manual(values = c("#FCDE05", "#FF5473", "#A696FF"))

fames_ab_plot
```


### Combined

```{r "fig1_combined", fig.width=8, fig.height=8, warning=FALSE}
# big version
gr_plot_top <- gr_plot + 
  theme(
    axis.text.x = element_blank(), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_blank(),
    legend.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold"),
    strip.text = element_text(face = "bold")
  )

brgdgt_ab_plot_middle <- 
  brgdgt_ab_plot + 
  theme(
    axis.text.x = element_blank(), 
    axis.title.x = element_blank(), 
    axis.ticks.x = element_blank(),
    strip.background.x = element_blank(), 
    strip.text.x = element_blank(),
    legend.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold")
  )

fames_ab_plot_bottom <-
  fames_ab_plot +
  theme(
    strip.background.x = element_blank(), 
    strip.text.x = element_blank(),
    legend.title = element_text(face = "bold"),
    axis.text = element_text(face = "bold")
  )

plot_grid(
  gr_plot_top + geom_point(size = 4, color = "black") + 
    labs(shape = "a: growth rates") + 
    theme(strip.background.y = element_blank(), strip.text.y = element_blank()), 
  brgdgt_ab_plot_middle + 
    labs(fill = "b: tetraethers") + 
    theme(strip.background.y = element_blank(), strip.text.y = element_blank()), 
  fames_ab_plot_bottom + 
    labs(fill = "c: FAMEs") + 
    theme(strip.background.y = element_blank(), strip.text.y = element_blank()), 
  align = "v", ncol = 1, axis = "lr", rel_heights = c(1.15, 0.95, 1.2)
)
```

```{r "fig1_combined_small_no_legend", fig.width=2.36, fig.height=3.15, warning=FALSE}
# small version
small_theme <- 
  theme(
    legend.position = "none", 
    text = element_text(size = 8),
    axis.text = element_text(size = 8),
    strip.text = element_text(size = 8),
    plot.margin = margin(t = 0, b = 0),
    panel.spacing.x = unit(1, "pt")
  )

plot_grid(
  gr_plot_top + small_theme + geom_point(size = 2, color = "black") + 
    theme(plot.margin = margin(t = 1)), 
  brgdgt_ab_plot_middle + small_theme, 
  fames_ab_plot_bottom + small_theme + theme(
    axis.text.x = element_text(size = 6)
  ), 
  align = "v", ncol = 1, axis = "lr", rel_heights = c(1.15, 0.95, 1.2)
)
```

## Figure 2: tetraether structures

### Spectra 

```{r "fig2_spectra_without_labels", fig.width=12, fig.height=12}
# define colors
ms_colors <- 
  tribble(
    ~compound,       ~type, ~color,
    "GDGT-1a",       "label", "#E69F00",
    "GDGT-1cisomer", "label", "#0072B2",
    "GTGT-1a",       "label", "#009E73",
    "GDGT-1a",       "formula", "black",
    "GDGT-1cisomer", "formula", "black",
    "GTGT-1a",       "formula", "black"
  ) %>%
  mutate(compound = as_factor(compound))

# add colors to ms2 and pick overall mass range
ms2_w_colors <- ms2_w_peaks %>%
  # keep order from the colors table
  mutate(compound = factor(compound, levels = levels(ms_colors$compound))) %>%
  # add colors 
  left_join(filter(ms_colors, type == "label"), by = "compound") %>%
  # select overall mass range to show
  filter(Mass >= 290, Mass <= 1050) 

# adjust to change the max shown on the y axis
max_y <- 0.18

# spectra plot
spectra_plot_without_labels <- 
  ms2_w_colors %>%
  ggplot() +
  # lines
  geom_line(aes(x = Mass, y = rel_intensity, color = color), size = 1.5) +
  # scales
  scale_color_identity(drop = FALSE) +
  scale_x_continuous(expand = c(0, 0), breaks = (0:10)*100) +
  scale_y_continuous(labels = function(x) paste0(100*x, " %")) +
  coord_cartesian(ylim = c(0.01, max_y)) +
  # facet
  facet_grid(compound ~ .) +
  # theme
  theme_figure(grid = FALSE, legend = FALSE) +
  expand_limits(x = 1050) + 
  # labels
  labs(x = "mass", y = "relative intensity")

spectra_plot_without_labels
```

```{r "fig2_spectra", fig.width=12, fig.height=12}
# load peaks to highlight
peaks <- readxl::read_excel("data/ms2_fragments.xlsx") %>%
  filter(!is.na(include) & include) %>%
  mutate(order = row_number()) %>%
  pivot_longer(cols = c(label, formula), names_to = "type", values_to = "label") %>%
  filter(!is.na(label), nchar(label) > 0) %>%
  mutate(
    peak_join_mass = sprintf("%.1f", peak),
    label = str_replace_all(label, "-", "\\\\,-"),
    label_tex = 
      ifelse(
        type == "formula",
        sprintf("$%s:\\,%s$", peak_join_mass, label),
        sprintf("$\\[%s\\]^+$", label)
      ) %>% latex2exp::TeX() %>% as.character()
  )

# process peaks to generate peak labels
peak_labels <- ms2_w_peaks %>%
  mutate(peak_join_mass = sprintf("%.1f", peak_mass)) %>%
  left_join(peaks, by = "peak_join_mass") %>% 
  filter(peak_nr > 0, !is.na(label_compound), label_compound == compound) %>%
  arrange(desc(type), desc(order)) %>%
  select(compound = label_compound, type, molecular_ion, 
         peak_join_mass, peak_mass, peak_rel_intensity, label_tex) %>%
  # calculate base position for label above top peak
  group_by(compound, peak_join_mass) %>%
  mutate(peak_rel_intensity = max(peak_rel_intensity), peak_mass = mean(peak_mass)) %>%
  unique() %>%
  mutate(y_offset = 0:(n() - 1)) %>%
  ungroup() %>%
  # keep order from the colors table
  mutate(compound = factor(compound, levels = levels(ms_colors$compound))) %>%
  # add colors
  left_join(ms_colors, by = c("compound", "type"))

# make plot with labels
spectra_plot <- spectra_plot_without_labels +
   # labels
  geom_text_repel(
    data = function(df) filter(peak_labels, compound %in% unique(df$compound)),
    mapping = 
      aes(
        x = peak_mass,
        y = ifelse(molecular_ion, max_y - 0.01, peak_rel_intensity), # + y_offset * 0.01, 
        hjust = ifelse(peak_mass > 900, 1, 0),
        vjust = ifelse(molecular_ion, 1, 0),
        label = label_tex, color = color
      ), 
    size = 4, parse = TRUE, nudge_y = 0.1, 
    seed = 42, min.segment.length = 0, force = 4
  ) 

spectra_plot
```

### Chromatograms

```{r "fig2_chromatograms", fig.width=6, fig.height=6}
# define colors
chroms_colors <- 
  tibble(
    # the order here determines panel and legend order
    compound = c("GDGT-1a", "GTGT-1a", "GDGT-1cisomer"),
    color = c("#E69F00", "#009E73", "#0072B2")
  ) %>% 
  crossing(type = c("samplenormalized", "stdnormalized")) %>%
  mutate(
    compound_type = paste(compound, type),
    color = ifelse(type == "stdnormalized", "thistle4", color) %>% as_factor(),
    compound = as_factor(compound)
  )

# plot
chroms_plot <- chroms_norm_tidy %>%
  # use as_factor to copy order in panels and legends
  mutate(compound = factor(compound, levels = levels(chroms_colors$compound))) %>%
  arrange(compound) %>%
  mutate(
    mz = as_factor(sprintf("m/z = %.0f", mz)),
    compound_type = paste(compound, type)
  ) %>%
  left_join(select(chroms_colors, compound_type, color), by = "compound_type") %>%
  ggplot() +
  aes(x = Time, y = value, color = color) +
  # data
  geom_line(size = 1.05) +
  # x acis breaks
  scale_x_continuous(breaks = c(30, 35, 40, 45, 50, 55)) + 
  # y axis starts at 0 but goes slightly above the max
  scale_y_continuous(expand = expansion(mult = c(0, 0.1))) +
  # manual color scale - every other is a standard
  scale_color_identity() +
  # plot all in one using a facet grid
  facet_grid(mz ~ .) +
  # always use this for easy defaults
  theme_figure(grid = FALSE) +
  # additionally disable the y axis ticks and legend
  theme(
    axis.text.y = element_blank(),
    axis.ticks.y = element_blank(),
    legend.position = "none"
  ) +
  # all labels in one
  labs(y = "Intensity (normalized)", x = "Time (minutes)")

chroms_plot
```

### Combined

```{r "fig2_combined_without_labels", fig.width=14, fig.height=12, warning=FALSE}
plot_grid(
  spectra_plot_without_labels + 
    theme(strip.background = element_blank(), strip.text = element_blank()) +
    labs(y = NULL),
  chroms_plot, 
  align = "h", ncol = 2, axis = "bt", rel_widths = c(4, 3)
)
```

```{r "fig2_combined", fig.width=14, fig.height=12, warning=FALSE}
plot_grid(
  spectra_plot + 
    theme(strip.background = element_blank(), strip.text = element_blank()) +
    labs(y = NULL),
  chroms_plot, 
  align = "h", ncol = 2, axis = "bt", rel_widths = c(4, 3)
)
```

## Figure S1: growth curves

```{r "figS1_growth_curves", fig.width = 10, fig.height = 6.5}
ggplot() + 
  aes(x = time.hours/24, y = OD, color = replicate) +
  geom_line(
    data = generate_logistic_curve(growth_rates, time = time.hours, N = OD) %>%
       mutate(panel = sprintf("%s: %.0f %% O2\n%s", bug_ID, `% O2`, `C source`) %>% as_factor())
  ) +
  geom_point(
    data = growth_curves_tidy %>%
      left_join(experiments, by = c("bug_ID", "exp_ID")) %>%
      mutate(panel = sprintf("%s: %.0f %% O2\n%s", bug_ID, `% O2`, `C source`) %>% as_factor())
  ) + 
  scale_color_brewer(palette = "Set1") +
  facet_wrap(~panel, scales = "free") +
  theme_figure(text_size = 14) +
  labs(x = "Time [days]")
```

# Tables

## Table S1: lipid data

```{r}
# lipids data table
lipids_table <- 
  lipids %>%
  arrange(desc(`% O2`), `C source`, replicate) %>%
   mutate(
    header = paste("%", compound),
    rel_amount.percent = 100 * rel_amount
  ) %>%
  select(organism = bug_ID, `% O2`, `C source`, replicate, header, rel_amount.percent) %>%
  pivot_wider(names_from = header, values_from = rel_amount.percent) %>%
  select(
    organism, `% O2`:replicate, 
    `% C8:0`, `% C14:0`, `% iso-C15:0` = `% iC15:0`, `% C15:0`, 
    `% C16:1`, `% C16:0`, `% iso-C17:0` = `% iC17:0`, `% C17:1`, 
    `% C17:0`, `% C18:1`, `% C18:0`, `% C20:0`, `% C22:0`, 
    `% squalene`, `% iso-diabolic acid` = `% iDA`,
    `% brGDGT Ia` = `% brGDGT_1a`, `% brGDGT Ic isomer` = `% brGDGT_1c_isomer`, 
    `% brGTGT Ia` = `% brGTGT_1a`, everything()
  ) %>%
  mutate(`% O2` = sprintf("%.0f", `% O2`))
lipids_table %>% export_to_excel(file = "tables/table_S1.xlsx")
lipids_table %>% knitr::kable()
```

## Table S2: growth data

```{r, warning=FALSE}
# growth rates data table
gr_table <- growth_rates %>%
  arrange(desc(`% O2`), `C source`, replicate) %>%
  select(organism = bug_ID, `% O2`, `C source`, replicate, `growth rate [1/d]` = r.1_d, `K [OD600]` = K) %>%
  mutate(`% O2` = sprintf("%.0f", `% O2`))
gr_table %>% export_to_excel(file = "tables/table_S2.xlsx")
gr_table %>% knitr::kable()
```