Formatting bioscreen growth curves

Author

Shane Hogle

Published

May 24, 2025

1 Introduction

Before the main experiment we grew each ancestral and evolved form of each species in 100% R2A media at a range of streptomycin concentrations. We measured the optical density of the cultures in a high-throughput platereader (Bioscreen) over 48 hours. From this data we can estimate the growth rate and the carrying capacity of each species in the different streptomycin concentrations.

We also grew all the ancestral and evolved forms of each species for 48 hours in 100% R2A media. We filtered the spent media of each species and evolved form, collected the filtrate, and grew each species/evo form in the filtrate of every other species. Again, we measured the optical density of the cultures grown on these filtrates in a high-throughput platereader (Bioscreen) over 48 hours.

In this notebook we will read the output from the bioscreen plate reader and format it for later plotting and analysis.

2 Setup

2.1 Libraries

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

2.2 Global variables

Show/hide code
data_raw <- here::here("_data_raw", "monocultures", "20240328_bioscreen")
data <- here::here("data", "monocultures", "20240328_bioscreen")

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

2.3 Functions

For plotting results

Show/hide code
plotplate <- function(df, dfxy, unsmoothed=TRUE, predicted=FALSE, plate, rows, cols, page){
  dffilt <- dplyr::filter(df, plate_name == {{ plate }})
  xyfilt <- if (!is.null(dfxy)){ left_join(dfxy, distinct(dffilt, id, bioscreen_well, plate_name), by = join_by(id)) %>% 
      drop_na()}
  
  ggplot(dffilt, aes(x = hours)) +
    list(
      ggplot2::geom_line(aes(y=OD600_rollmean), color = "blue"), 
      if (unsmoothed) {ggplot2::geom_line(aes(y=OD600), color = "orange", lty = 2)},
      if (predicted) {ggplot2::geom_line(aes(y=predicted), color = "orange")}, 
      if (!is.null(dfxy)) {ggplot2::geom_point(data = xyfilt, aes(x = x, y = y), color = "red", size = 2)},
      ggplot2::labs(x = "Hours", y = "OD600"), 
      ggplot2::scale_x_continuous(breaks = seq(0, 48, 12), labels = seq(0, 48, 12)), 
      ggforce::facet_wrap_paginate(~ bioscreen_well, nrow = rows, ncol = cols, page = page, scales = "free_y"), 
      ggplot2::theme(axis.text = element_text(size = 5))
    )
}

3 Tidying growth curves from bioscreen

3.1 Deal with encoding issues

First we need to know what encoding the file has

Show/hide code
readr::guess_encoding(
  here::here(data_raw,
    "bioscreen_strep_01.txt"
  )
)

Bioscreen files are in UTF-16E.

3.1.1 List files

Show/hide code
files <- fs::dir_ls(
  path = data_raw,
  all = FALSE,
  recurse = TRUE,
  type = "file",
  glob = "*.txt",
  regexp = NULL,
  invert = FALSE,
  fail = TRUE
)

files_named <- rlang::set_names(files, nm = purrr::map_chr(files, fs::path_file) %>% str_remove(".txt"))

3.1.2 Read file using base R

Must use read.delim as it is the only tool that I could find that handled the UTF-16LE correctly

Show/hide code
raw_df <- purrr::map(
  files_named,
  read.delim,
  # this is the encoding we deduced above
  fileEncoding = "UTF-16LE",
  skipNul = TRUE,
  # telling it there is no header for colulmn names
  header = FALSE,
  sep = " ",
  colClasses = "character",
  strip.white = TRUE,
  na.strings = "",
  # skip the first 3 lines because
  skip = 3,
  quote = "",
  allowEscapes = FALSE,
  comment.char = "",
) %>% purrr::list_rbind(names_to = "plate_name")

3.1.3 Fix locale issues

Show/hide code
raw_df_decimal <- readr::type_convert(
  raw_df,
  locale = readr::locale(decimal_mark = ",", grouping_mark = "."),
  cols(
    .default = col_double(),
    plate_name = col_character(),
    V1 = col_time(format = ""),
    V203 = col_logical()
  )
) %>%
  tibble::as_tibble()

3.1.4 Fix weird time encoding

The V1 variable has time encoded as Hr:Min:Sec which we want to convert to just seconds

The rounding step is necessary because the bioscreen actually doesn’t output consistent intervals. Sometimes it is 00:30:06 and other times it is 00:30:05. This becomes a problem later on when trying to combine multiple runs at once

Show/hide code
round_any <- function(x, accuracy, f=round){
  f(x/ accuracy) * accuracy
}

raw_df_decimal_sec <- dplyr::mutate(raw_df_decimal, 
                                    V1 = round_any(lubridate::period_to_seconds(lubridate::hms(V1)), 100)
                                    )

3.1.5 Remove extraneous columns.

Some of the columns have no useful data whatsoever. We need to get rid of those.

Show/hide code
dplyr::select(raw_df_decimal_sec, tidyselect::last_col(offset = 1), tidyselect::last_col())
Show/hide code
dplyr::select(raw_df_decimal_sec, tidyselect::num_range("V", 1:4))

The last columns is V203 and has nothing. Column V2 is just blank. Now we’ll just keep the columns we need.

Show/hide code
gcurves_slurped_fmt <- raw_df_decimal_sec %>% 
  dplyr::select(plate_name, V1, V3:V202) %>% 
  dplyr::mutate(hours = lubridate::time_length(V1, unit = "hours")) %>%
  tidyr::pivot_longer(c(-plate_name, -hours, -V1), names_to = "well", values_to = "OD600") %>%
  # here we remove the V part and subtract 2 so that bioscreen well numbers will go from 1:100 per plate
  dplyr::mutate(bioscreen_well = as.numeric(str_remove(well, "V"))-2) %>% 
  dplyr::rename(seconds = V1) 

4 Format growth curves

4.1 Read metadata

Streptomycin and filtrate conditions

Show/hide code
gcurves_slurped_fmt_md <- readr::read_tsv(
  here::here(data_raw, "plate_conditions.tsv"),
  col_names = TRUE,
  cols(
    plate_name = col_character(),
    plate_number = col_double(),
    bioscreen_well = col_double(),
    replicate = col_double(),
    sp_hist = col_character(),
    strep_conc = col_double(),
    sp_filtrate = col_character()
  )
)

4.2 Join with metadata to remove ununsed samples

Show/hide code
gcurves <- dplyr::left_join(
  gcurves_slurped_fmt_md,
  gcurves_slurped_fmt,
  by = dplyr::join_by(plate_name, bioscreen_well)
)

One thing I’ve realized is that many methods for inferring growth rates struggle when the density of observations is too high (e.g., one measurement every 5 minutes). In reality I’ve found that taking one measurement every 15 minutes is sufficient. Here we thin it out so that measurements are in 20 minute intervals. This seems to improve the fitting procedure a lot without much of a cost.

Show/hide code
gcurves_thin <- gcurves %>% 
  # 1200 is 20 minutes so by ensuring modulo = 0 we include only time points
  # 0, 20, 40, 60 minutes and so on...
  dplyr::filter(seconds %% 1200 == 0)

Now we’ll do some smoothing to reduce the “jaggedness” of the curves a bit. We use the slider package with a 5 point rolling mean for each focal observation we take the mean including the focal point and two points before and after.

Show/hide code
gcurves_thin_sm <- gcurves_thin %>% 
  dplyr::group_by(plate_name, bioscreen_well) %>% 
  dplyr::mutate(OD600_rollmean = slider::slide_dbl(OD600, mean, .before = 2, .after = 2)) %>% 
  ungroup()

readr::write_tsv(gcurves_thin_sm, here::here(data, "gcurves_formatted_thinned.tsv"))

5 Inspect growth curves

5.1 bioscreen_strep_01

Figure 1: Bioscreen plate “bioscreen_strep_01” wells 1-100. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

Figure 2: Bioscreen plate “bioscreen_strep_01” wells 101-200. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

5.2 bioscreen_strep_02

Figure 3: Bioscreen plate “bioscreen_strep_02” wells 1-70. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

Figure 4: Bioscreen plate “bioscreen_strep_02” wells 71-140. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

5.3 bioscreen_strep_03

Figure 5: Bioscreen plate “bioscreen_strep_03” wells 1-70 Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

Figure 6: Bioscreen plate “bioscreen_strep_03” wells 71-140. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

5.4 bioscreen_pairwise_filtrates_01

Figure 7: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 1-80 Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

Figure 8: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 81-160 Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

5.5 bioscreen_pairwise_filtrates_02

Figure 9: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 1-80. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

Figure 10: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 81-160. Blue line is smoothed with a moving average window of 9 points. Orange is non-smoothed

5.6 Conclusions

Growth curves all look reasonable. Can proceed with the analysis.

6 Growth curve statistics

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

Attaching package: 'DescTools'
The following object is masked _by_ '.GlobalEnv':

    %nin%

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
gcurves_thin_sm <- gcurves_thin_sm %>% 
  # make uniq id
  mutate(id = paste0(plate_name, "|", bioscreen_well))

6.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.

6.1.1 Fit

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

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

6.1.2 Results

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

6.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") %>% 
  dplyr::mutate(id = stringr::str_remove_all(id, "bioscreen_pairwise_|bioscreen_"))

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(gcurves_thin_sm, by = dplyr::join_by(id, hours)) %>% 
  dplyr::mutate(id = stringr::str_remove_all(id, "bioscreen_pairwise_|bioscreen_")) %>%
  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()

6.1.4 Plot

6.1.4.1 bioscreen_strep_01

Figure 11: Bioscreen plate “bioscreen_strep_01” wells 1-100. 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.

Figure 12: Bioscreen plate “bioscreen_strep_01” wells 101-200. 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.

6.1.4.2 bioscreen_strep_02

Figure 13: Bioscreen plate “bioscreen_strep_02” wells 1-70. 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.

Figure 14: Bioscreen plate “bioscreen_strep_02” wells 71-140. 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.

6.1.4.3 bioscreen_strep_03

Figure 15: Bioscreen plate “bioscreen_strep_03” wells 1-70. 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.

Figure 16: Bioscreen plate “bioscreen_strep_03” wells 71-140. 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.

6.1.4.4 bioscreen_pairwise_filtrates_01

Figure 17: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate. Red dot is hours and OD600 at which maximum growth rate is reached.

Figure 18: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate. Red dot is hours and OD600 at which maximum growth rate is reached.

6.1.4.5 bioscreen_pairwise_filtrates_02

Figure 19: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate. Red dot is hours and OD600 at which maximum growth rate is reached.

Figure 20: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate. Red dot is hours and OD600 at which maximum growth rate is reached.

6.2 Linear model based estimate

This method determine maximum growth rates from the log-linear part of a growth curve using a heuristic approach similar to the “growth rates made easy” method of Hall et al. (2013).

From growthrates documentation:

The algorithm works as follows:

1. Fit linear regressions to all subsets of h consecutive data points. If for example h = 5, fit a linear regression to points 1 . . . 5, 2 . . . 6, 3. . . 7 and so on. The method seeks the highest rate of exponential growth, so the dependent variable is of course log-transformed.

2. Find the subset with the highest slope bmax and include also the data points of adjacent subsets that have a slope of at least quota · bmax, e.g. all data sets that have at least 95% of the maximum slope.

3. Fit a new linear model to the extended data window identified in step 2.

6.2.1 Fit

Show/hide code
set.seed(45278)
many_linear <- growthrates::all_easylinear(OD600_rollmean ~ hours | id, data = gcurves_thin_sm)

readr::write_rds(many_linear, here::here(data, "linear_fits"))

6.2.2 Results

Show/hide code
many_linear_res <- growthrates::results(many_linear)

6.2.3 Predictions

Show/hide code
many_linear_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(gcurves_thin_sm, by = dplyr::join_by(id, hours)) %>% 
  dplyr::mutate(id = stringr::str_remove_all(id, "bioscreen_pairwise_|bioscreen_")) %>%
  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()

6.2.4 Plot

6.2.4.1 bioscreen_strep_01

Figure 21: Bioscreen plate “bioscreen_strep_01” wells 1-100. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

Figure 22: Bioscreen plate “bioscreen_strep_01” wells 101-200. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

6.2.4.2 bioscreen_strep_02

Figure 23: Bioscreen plate “bioscreen_strep_02” wells 1-70. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

Figure 24: Bioscreen plate “bioscreen_strep_02” wells 71-140. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

6.2.4.3 bioscreen_strep_03

Figure 25: Bioscreen plate “bioscreen_strep_03” wells 1-70. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

Figure 26: Bioscreen plate “bioscreen_strep_03” wells 71-140. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

6.2.4.4 bioscreen_pairwise_filtrates_01

Figure 27: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

Figure 28: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

6.2.4.5 bioscreen_pairwise_filtrates_02

Figure 29: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

Figure 30: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is slope of max predicted growth rate from a piecewise linear model fit.

6.3 Baranyi parametric fit estimate

This uses the growth model of Baranyi and Roberts (1995) written as analytical solution of the system of differential equations.

From growthrates documentation:

The version of the equation used in this package has the following form:

\(A = time + \frac{1}{mumax} \times \log((−\mu_{max} \times time) +exp(−h_{0})−exp(−\mu_{max} \times time−h0))\)

\(\log(y) = \log(y_{0}) + \mu_{max} \times A − log(1 + \frac{\exp(\mu_{max} \times A) − 1}{\exp(\log(K) − \log(y_{0})})\)

6.3.1 Fit

Show/hide code
set.seed(783622)

# starting values and limits for the baranyi model
p   <- c(y0 = 0.03, mumax = .1, K = 0.1, h0 = 1)
lower   <- c(y0 = 0.001, mumax = 1e-2, K = 0.005, h0 = 0)
upper   <- c(y0 = 0.1,   mumax = 1,    K = 0.9,   h0 = 10)

many_baranyi <- growthrates::all_growthmodels(
                   OD600_rollmean ~ growthrates::grow_baranyi(hours, parms) | id,
                   data = gcurves_thin_sm,
                   p = p, lower = lower, upper = upper,
                   transform = "log", ncores = 8)

readr::write_rds(many_baranyi, here::here(data, "barayani_fits"))

6.3.2 Results

Show/hide code
many_baranyi_res <- growthrates::results(many_baranyi)

6.3.3 Predictions

Show/hide code
many_baranyi_fitted <- purrr::map(many_baranyi@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(gcurves_thin_sm, by = dplyr::join_by(id, hours)) %>% 
  dplyr::mutate(id = stringr::str_remove_all(id, "bioscreen_pairwise_|bioscreen_"))

6.3.4 Plot

6.3.4.1 bioscreen_strep_01

Figure 31: Bioscreen plate “bioscreen_strep_01” wells 1-100. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

Figure 32: Bioscreen plate “bioscreen_strep_01” wells 101-200. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

6.3.4.2 bioscreen_strep_02

Figure 33: Bioscreen plate “bioscreen_strep_02” wells 1-70. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

Figure 34: Bioscreen plate “bioscreen_strep_02” wells 71-140. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

6.3.4.3 bioscreen_strep_03

Figure 35: Bioscreen plate “bioscreen_strep_03” wells 1-70. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

Figure 36: Bioscreen plate “bioscreen_strep_03” wells 71-140. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

6.3.4.4 bioscreen_pairwise_filtrates_01

Figure 37: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

Figure 38: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

6.3.4.5 bioscreen_pairwise_filtrates_02

Figure 39: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

Figure 40: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Baranyi 1995 model.

6.4 Huang parametric fit estimate

Huangs growth model (ref 1, ref 2) written as analytical solution of the differential equations.

From growthrates documentation:

The version of the equation used in this package has the following form:

\(B = time + \frac{1}{alpha} \times \log(\frac{1 + exp(−\alpha \times (time - \lambda))}{1 + exp(\alpha \times \lambda})\)

\(\log(y) = \log(y_{0}) + \log(K) − \log(y_{0} + (K − y_{0}) \times exp(−\mu_{max} \times B))\)

In contrast to the original publication, all parameters related to population abundance (y, y0, K) are given as untransformed values. They are not log-transformed. In general, using log-transformed parameters would indeed be a good idea to avoid the need of constained optimization, but tests showed that box-constrained optimization worked resonably well. Therefore, handling of optionally log-transformed parameters was removed from the package to avoid confusion. If you want to discuss this, please let me know.

6.4.1 Fit

Show/hide code
set.seed(123784)

# starting values and limits for the huang model
p   <- c(y0 = 0.03, mumax = .1, K = 0.1, alpha = 1.5, lambda = 3)
lower   <- c(y0 = 0.001, mumax = 1e-2, K = 0.005, alpha = 0.1, lambda = 0.03)
upper   <- c(y0 = 0.1,   mumax = 1,    K = 0.9,   alpha = 15, lambda = 15)

many_huang <- growthrates::all_growthmodels(
                   OD600_rollmean ~ growthrates::grow_huang(hours, parms) | id,
                   data = gcurves_thin_sm,
                   p = p, lower = lower, upper = upper,
                   transform = "log", ncores = 8)

readr::write_rds(many_huang, here::here(data, "huang_fits"))

6.4.2 Results

Show/hide code
many_huang_res <- growthrates::results(many_huang)

6.4.3 Predictions

Show/hide code
many_huang_fitted <- purrr::map(many_huang@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(gcurves_thin_sm, by = dplyr::join_by(id, hours)) %>% 
  dplyr::mutate(id = stringr::str_remove_all(id, "bioscreen_pairwise_|bioscreen_"))

6.4.4 Plot

6.4.4.1 bioscreen_strep_01

Figure 41: Bioscreen plate “bioscreen_strep_01” wells 1-100. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

Figure 42: Bioscreen plate “bioscreen_strep_01” wells 101-200. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

6.4.4.2 bioscreen_strep_02

Figure 43: Bioscreen plate “bioscreen_strep_02” wells 1-70. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

Figure 44: Bioscreen plate “bioscreen_strep_02” wells 71-140. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huangv model.

6.4.4.3 bioscreen_strep_03

Figure 45: Bioscreen plate “bioscreen_strep_03” wells 1-70. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

Figure 46: Bioscreen plate “bioscreen_strep_03” wells 71-140. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

6.4.4.4 bioscreen_pairwise_filtrates_01

Figure 47: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

Figure 48: Bioscreen plate “bioscreen_pairwise_filtrates_01” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

6.4.4.5 bioscreen_pairwise_filtrates_02

Figure 49: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 1-80. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

Figure 50: Bioscreen plate “bioscreen_pairwise_filtrates_02” wells 81-160. Blue line is smoothed with a moving average window of 5 points. Orange is fitted to the Huang model.

6.5 AUC

Calculates AUC using DescTools package

Show/hide code
many_auc_res <- gcurves_thin_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(gcurves_thin_sm, id, plate_name:sp_filtrate))) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(auc, max_od, min_od, .after="sp_filtrate")
Joining with `by = join_by(id)`

7 Write all output

Show/hide code
readr::write_tsv(many_auc_res, here::here(data, "gcurve_auc_results.tsv")) 
Show/hide code
many_spline_res %>% 
  dplyr::left_join(dplyr::distinct(dplyr::select(gcurves_thin_sm, id, plate_name:sp_filtrate)), by = join_by(id)) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(y0:r2, .after="sp_filtrate") %>% 
  readr::write_tsv(here::here(data, "gcurve_spline_results.tsv"))
Show/hide code
many_linear_res %>% 
  dplyr::left_join(dplyr::distinct(dplyr::select(gcurves_thin_sm, id, plate_name:sp_filtrate)), by = join_by(id)) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(y0:r2, .after="sp_filtrate") %>% 
  readr::write_tsv(here::here(data, "gcurve_linear_results.tsv"))
Show/hide code
many_baranyi_res %>% 
  dplyr::left_join(dplyr::distinct(dplyr::select(gcurves_thin_sm, id, plate_name:sp_filtrate)), by = join_by(id)) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(y0:r2, .after="sp_filtrate") %>% 
  readr::write_tsv(here::here(data, "gcurve_baranyi_results.tsv"))
Show/hide code
many_huang_res %>% 
  dplyr::left_join(dplyr::distinct(dplyr::select(gcurves_thin_sm, id, plate_name:sp_filtrate)), by = join_by(id)) %>% 
  dplyr::select(-id) %>% 
  dplyr::relocate(y0:r2, .after="sp_filtrate") %>% 
  readr::write_tsv(here::here(data, "gcurve_huang_results.tsv"))