In maomlab/BayesPharma: Tools for Bayesian Analysis of Non-Linear Pharmacology Models

# inspired by https://www.jumpingrivers.com/blog/knitr-default-options-settings-hooks/
knitr::opts_chunk$set(
  echo = FALSE,
  cache = TRUE,
  cache.path = "cache/apply_model_tQ_nSmase2_kinetics/",
  fig.path = "apply_model_tQ_nSmase2_kinetics_files/",
  fig.retina = 2, # Control using dpi
  fig.width = 6,  # generated images
  fig.height = 5, # generated images
  fig.pos = "t",  # pdf mode
  fig.align = "center",
  dpi = if (knitr::is_latex_output()) 72 else 300, 
  out.width = "100%")

suppressWarnings(suppressMessages(library(tidyverse)))
library(BayesPharma)
library(GeomIndicator)
suppressWarnings(suppressMessages(library(tidybayes)))
library(deSolve)
library(tictoc)

data_path <- "~/Google Drive/My Drive/momeara_lab/Drug Discovery/nSMase2/Pharmacology Analysis/raw_data/08262022_DKswitch.cpds.xlsx"

data <- readxl::read_excel(
  path = data_path,
  sheet = "Data 08-28-22-151019_kinetic",
  skip = 8) |>
  dplyr::filter(!is.na(Dose_uM)) |>
  tidyr::fill(Temp) |>
  tidyr::pivot_longer(
    cols = c(-Dose_uM, -Time_min, -Temp),
    names_to = "Drug",
    values_to = "RFU")

data <- data |>
  dplyr::left_join(
    data |>
      dplyr::distinct(Drug, Dose_uM) |>
      dplyr::mutate(series_index = dplyr::row_number()),
    by = c("Drug", "Dose_uM"))

data <- data |>
  dplyr::mutate(
    P = RFU / 3000,
    ET = 1,
    ST = Dose_uM,
    time = Time_min / 60)

ggplot2::ggplot(data = data) + 
  ggplot2::theme_bw() +
  ggplot2::geom_line(
    mapping = ggplot2::aes(
      x = time,
      y = P,
      group = log10(Dose_uM),
      color = Dose_uM)) + 
  # ggplot2::geom_smooth(
  #   data = data |>
  #     dplyr::filter(
  #       treatment != "Blank",
  #       time_s <= 200),
  #   mapping = ggplot2::aes(
  #     x = time_s,
  #     y = product_concentration * 0.7,
  #     group = paste0(substrate_concentration_mg_ml, "_", replica)),
  #   method = "lm",
  #   formula = y ~ x,
  #   color = "orange") +
  ggplot2::facet_wrap(
    facets = dplyr::vars(Drug)) +
  ggplot2::theme(legend.position = "bottom")

data_single <- tQ_model_generate(
  time = seq(0.00, 3, by=.05),
  kcat = 3,
  kM = 5,
  ET = 10,
  ST = 10) |>
  as.data.frame() |>
  dplyr::rename(P_true = 2) |>
  dplyr::mutate(
    P = rnorm(dplyr::n(), P_true, 0.5), # add some observational noise
    ST = 10, ET = 10)

data_339 <- data |> dplyr::filter(Drug == 339)
plot_lineweaver_burk <- function(data){

  plot_data <- data_339 |>
    dplyr::filter(Dose_uM > 0) |>
    dplyr::group_by(Drug, Dose_uM) |>
    dplyr::summarize(
      Vmax_inverse = 1/max(P),
      S_inverse = 1/Dose_uM[1],
      .groups = "drop")

  linear_fit <- lm(
    data = plot_data |> dplyr::filter(Vmax_inverse > 0.8),
    formula = Vmax_inverse ~ S_inverse)

  # Derive formula for x-intercept for a line
  # The formula for a line is y = mx + b so if y = 0 then x = -b/m
  x_intercept <- -coef(linear_fit)[1] / coef(linear_fit)[2]
  Km_data = data.frame(
    S_inverse = x_intercept,
    Vmax_inverse = 0,
    label = paste0(
    "-1/Km = ", signif(x_intercept, 2), "\n",
    "Km = ", signif(-1/x_intercept, 2)))

  line_data <- tibble::tibble(
    S_inverse = seq(x_intercept, max(plot_data$S_inverse), length.out = 100))
  line_data$Vmax_inverse <- predict(linear_fit, newdata = line_data)


  plot <- ggplot2::ggplot(data = plot_data) +
    ggplot2::theme_bw() +
    ggplot2::geom_vline(xintercept = 0) +
    ggplot2::geom_hline(yintercept = 0) +
    ggplot2::geom_point(
      mapping = ggplot2::aes(
        x = S_inverse,
        y = Vmax_inverse),
      color = "blue") +
    ggplot2::geom_point(
      data = Km_data,
      mapping = ggplot2::aes(x=S_inverse, y=Vmax_inverse),
      size = 1.3,
      color = "orange") +
    ggrepel::geom_label_repel(
      data = Km_data,
      mapping = ggplot2::aes(x = S_inverse, y = Vmax_inverse, label = label)) +
    ggplot2::geom_line(
      data = line_data,
      mapping = ggplot2::aes(x = S_inverse, y = Vmax_inverse),
      color = "red") +
    ggplot2::facet_wrap(facets = dplyr::vars(Drug), scales = "free_y")
}

scale_estimates <- data |>
  dplyr::group_by(Drug, ET) |>
  dplyr::summarize(
    Vmax = max(RFU),
    Kcat = max(RFU) / ET[1],
    .groups = "drop")

# n_series  subsample adapt_delta runtime comment
# 1 NA  .999  12.6  
# 2 NA  .999  58.3  
# 2 NA  .99   28.5  Seems more reliable
# 2 NA  .9    14.5  Some fraction of the chains fail to sample at all
# 3 40  .999  302.3 Only two chains ran at a time? 
# 3 40  .9    681.6 Sat on zero samples for a long time and then flew
# 4 40  .9    665.4 Sat on zero samples for a long time and then flew
# 6 NA  .9    653.0 Sat on zero samples for a long time and then flew
# all NA  .9  20595.8
# all NA  .0  13351.1 Tighter prior for logkM = normal(-3, 2)

runtime <- system.time({
  model <- brms::brm(
      data = data |>
        dplyr::filter(Drug == 339) |>
        dplyr::filter(ST > 0) |>
        dplyr::group_by(series_index) |>
        dplyr::arrange(series_index, time),
      formula = brms::brmsformula(
        P ~ tQ_multiple(series_index, time, kcat, 10^logkM, ET, ST),
        kcat ~ 1,
        logkM ~ 1,
        nl = TRUE,
        loop = FALSE),
      stanvars = BayesPharma::tQ_stanvar(),
      prior = c(
        brms::prior(prior = gamma(9, 100), lb = 0, nlpar = "kcat"),
        brms::prior(prior = normal(-3, 2), nlpar = "logkM")),
      init = function() list(b_kcat = .09, b_logkM = -3, sigma = 1),
      chains = 4,
      iter = 4000,
      cores = 4,
      #control = list(adapt_delta = 0.9),
      #backend = "cmdstanr",
      algorithm = "meanfield")
})

model_prior <- model |>
  update(sample_prior = "only")

draws_prior <- model_prior |>
  tidybayes::tidy_draws() |>
  tidybayes::gather_variables() |>
  dplyr::filter(.variable |> stringr::str_detect("b_(kcat|logkM)")) |>
  tidyr::separate(
    .variable,
    into = c("class", "parameter", "treatment"), sep = "_") |>
  dplyr::mutate(treatment = treatment |> stringr::str_replace("treatment", ""))

draws_posterior <- model |>
  tidybayes::tidy_draws() |>
  tidybayes::gather_variables() |>
  dplyr::filter(.variable |> stringr::str_detect("b_(kcat|logkM)")) |>
  tidyr::separate(
    .variable,
    into = c("class", "parameter", "treatment"), sep = "_") |>
  dplyr::mutate(treatment = treatment |> stringr::str_replace("treatment", ""))

draws_prior_pairs <- draws_prior |>
  tidyr::pivot_wider(
    id_cols = c(".chain", ".iteration", ".draw", "treatment"),
    names_from = "parameter",
    values_from = ".value")

draws_posterior_pairs <- draws_posterior |>
  tidyr::pivot_wider(
    id_cols = c(".chain", ".iteration", ".draw", "treatment"),
    names_from = "parameter",
    values_from = ".value")

ggplot2::ggplot() +
  ggplot2::theme_bw() +
  ggplot2::geom_point(
    data = draws_prior_pairs,
    mapping = ggplot2::aes(
      x = log10(kcat),
      y = logkM),
    color = "orange",
    size = 0.8,
    shape = 16,
    alpha = .6) +
  ggplot2::geom_point(
    data = draws_posterior_pairs,
    mapping = ggplot2::aes(
      x = log10(kcat),
      y = logkM),
    color = "blue",
    size = .8,
    shape = 16,
    alpha = .3) +
  ggplot2::scale_x_continuous("Catalytic constant: log(kcat)") +
  ggplot2::scale_y_continuous("Michaelis constant: log(kM)") +
  ggplot2::facet_wrap(
    facets = dplyr::vars(treatment),
    nrow = 1,
    scales = "free")

#BayesPharma::plot_posterior_draws(
#  model = model)

maomlab/BayesPharma documentation built on Aug. 24, 2024, 8:45 a.m.