Formatting co-culture growth curves from Synergy H1

Author

Shane Hogle

Published

July 22, 2025

Abstract
The ancestral (streptomycin sensitive) and evolved (streptomycin resistant) forms of HAMBI_1287 and HAMBI_1977 were grown on adjacent 0.22 um membrane-separated wells in a Verasys co-culture plate. Note that due to the unique format of the plate rows A and H do not contain any data. Measurements were collected with a synergy H1 multimode microplate reader.

1 Setup

1.1 Libraries

Show/hide code
library(here)
library(tidyverse)
library(readxl)
library(stringr)
library(lubridate)
library(fs)
library(ggforce)
library(slider)
source(here::here("R", "utils_gcurves.R"))

1.2 Global variables

Show/hide code
data_raw <- here::here("_data_raw", "coculture_plate")
data <- here::here("data", "coculture_plate")

# make processed data directory if it doesn't exist
fs::dir_create(data)

2 Tidying growth curves from Synergy H1 multimode microplate reader

2.1 Read growth curves

Show/hide code
plate01 <- readxl::read_xlsx(here::here(data_raw, "20250714_init", "rawdata_versarys.xlsx"), sheet = 3, skip = 1) %>% 
  # set interval start to be first cell and make all intervals relative to that
  # use time_length to just create an hours variable of type numeric
  mutate(seconds = lubridate::time_length(lubridate::interval(Time[1], Time), unit = "second")) %>% 
  tidyr::pivot_longer(c(-seconds, -Time), names_to = "well", values_to = "OD600") %>%
  mutate(hours = lubridate::time_length(seconds, unit = "hours")) %>% 
  # converting the well format so it matches the samplesheet
  mutate(well = paste0(str_extract(well, "^[A-H]"), str_pad(str_extract(well, "\\d+"), width = 2, pad = "0", side = "left"))) %>% 
  dplyr::select(seconds, hours, well, OD600) %>% 
  # create a plate variable for later combining
  mutate(plate_name = "plate01")

3 Format growth curves

3.1 Read sample metadata

Show/hide code
samplesheet01 <- readxl::read_xlsx(here::here(data_raw, "20250714_init", "samplesheet.xlsx")) %>% 
  mutate(strainID = paste0(strain, str_to_upper(str_sub(evolution, 1, 1)))) %>% 
  mutate(plate_name = "plate01")

3.2 Join with metadata to remove ununsed samples

Show/hide code
coculture_gcurves_sm <-plate01 %>% 
  left_join(samplesheet01, by = join_by(well, plate_name)) %>% 
  dplyr::group_by(plate_name, well) %>% 
  dplyr::mutate(OD600_rollmean = slider::slide_dbl(OD600, mean, .before = 2, .after = 2)) %>% 
  ungroup() %>% 
  relocate(OD600_rollmean, .after = "OD600") %>% 
  filter(str_detect(well, "^A|^H", negate = TRUE))

readr::write_tsv(coculture_gcurves_sm, here::here(data, "coculture_gcurves_smooth.tsv"))

4 Inspect growth curves

4.1 plate01 (Ancestral HAMBI_1287 ANC)

This plate contains three replicates for the ancestral form of HAMBI_1287

Figure 1: Growth curves for the first coculture plate. X-axis is time in hours (48 hour incubation). Y axis is the absorbance scaled for each well. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed. Note that rows A and H are meaningless because of the construction of the plate. Columns 2 and 10 were not inoculated with bacteria.

5 Growth curve statistics

Show/hide code
library("growthrates")
Loading required package: lattice
Loading required package: deSolve
Show/hide code
library("DescTools")

Using the tool growthrates to estimate mu_max. I have found this works a lot better the gcplyr and is more convenient than using another tool outside of R. Nonparametric estimate growth rates by spline is very fast. Fitting to a model takes more time resources. Generally it is best to try multiple approaches and to visualize/check the data to make sure it makes sense.

Show/hide code
coculture_gcurves_sm <- coculture_gcurves_sm %>% 
  # make uniq id
  mutate(id = paste0(plate_name, "|", well))

5.1 Spline based estiamte

Smoothing splines are a quick method to estimate maximum growth. The method is called nonparametric, because the growth rate is directly estimated from the smoothed data without being restricted to a specific model formula.

From growthrates documentation:

The method was inspired by an algorithm of Kahm et al. (2010), with different settings and assumptions. In the moment, spline fitting is always done with log-transformed data, assuming exponential growth at the time point of the maximum of the first derivative of the spline fit. All the hard work is done by function smooth.spline from package stats, that is highly user configurable. Normally, smoothness is automatically determined via cross-validation. This works well in many cases, whereas manual adjustment is required otherwise, e.g. by setting spar to a fixed value [0, 1] that also disables cross-validation.

5.1.1 Fit

Show/hide code
set.seed(45278)
many_spline <- growthrates::all_splines(OD600_rollmean ~ hours | id, data = coculture_gcurves_sm, spar = 0.5)

readr::write_rds(many_spline, here::here(data, "coculture_spline_fits"))

5.1.2 Results

Show/hide code
many_spline_res <- growthrates::results(many_spline)

5.1.3 Predictions

Show/hide code
many_spline_xy <- purrr::map(many_spline@fits, \(x) data.frame(x = x@xy[1], y = x@xy[2])) %>% 
  purrr::list_rbind(names_to = "id") 

many_spline_fitted <- purrr::map(many_spline@fits, \(x) data.frame(x@FUN(x@obs$time, x@par))) %>% 
  purrr::list_rbind(names_to = "id") %>% 
  dplyr::rename(hours = time, predicted = y) %>% 
  dplyr::left_join(coculture_gcurves_sm, by = dplyr::join_by(id, hours)) %>% 
  dplyr::group_by(id) %>% 
  # this step makes sure we don't plot fits that go outside the range of the data
  dplyr::mutate(predicted = dplyr::if_else(dplyr::between(predicted, min(OD600_rollmean), max(OD600_rollmean)), predicted, NA_real_)) %>% 
  dplyr::ungroup()

5.1.4 Plot

5.1.4.1 Plate01

Figure 2: As in Figure 1. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from the first derivative of a smoothing spline. Red dot is hours and OD600 at which maximum growth rate is reached.

5.2 AUC

Calculates AUC using DescTools package

Show/hide code
many_auc_res <- coculture_gcurves_sm %>% 
  dplyr::summarize(auc = DescTools::AUC(hours, OD600_rollmean),
            max_od = max(OD600_rollmean),
            min_od = min(OD600_rollmean),
            .by = id) %>% 
  dplyr::left_join(dplyr::distinct(dplyr::select(coculture_gcurves_sm, plate_name:id)), by = join_by(id)) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(auc, max_od, min_od, .after="strainID")

6 Write all output

Show/hide code
readr::write_tsv(many_auc_res, here::here(data, "coculture_gcurve_auc_results.tsv")) 
Show/hide code
many_spline_res %>% 
  dplyr::left_join(dplyr::distinct(dplyr::select(coculture_gcurves_sm, plate_name:id)), by = join_by(id)) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(y0:r2, .after="strainID") %>% 
  readr::write_tsv(here::here(data, "coculture_gcurve_spline_results.tsv"))