## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()
set.seed(123)

## -----------------------------------------------------------------------------
library(ggplot2)
library(dplyr)
library(MetaboDynamics)

# Simulate data resembling experimental observations
mu_mean <- 12; sd_mean <- 5
prior_mean_abundance <- c(mu_mean, sd_mean)
prior_sd_abundance <- 2

n_metabolites <- 300
mu <- rnorm(n_metabolites, mean = prior_mean_abundance[1], sd = prior_mean_abundance[2])
sigma <- rexp(n_metabolites, rate = 1/prior_sd_abundance)

n_observations <- 3
abundances <- data.frame(
  metabolite = rep(seq_len(n_metabolites), each = n_observations),
  abundance = NA
)

for (i in seq_len(n_metabolites)){
  abundances[abundances$metabolite == i, ]$abundance <- rnorm(3, mean = mu[i], sd = sigma[i])
}
abundances$abundance <- exp(abundances$abundance)

# Visualize: Black = Data, Red = Prior
ggplot(abundances, aes(x = log(abundance))) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(aes(x = rnorm(nrow(abundances), mean = 12, sd = 10)), 
               color = "red", linewidth = 1) +
  labs(title = "Distribution of log-transformed metabolite abundances",
       subtitle = "Black: Observed Data | Red: Prior N(12, 10)")

## -----------------------------------------------------------------------------
sds <- abundances %>%
  group_by(metabolite) %>%
  summarise(sd_log_abundance = sd(log(abundance)))

# Visualize: Black = Data SDs, Red = Prior
ggplot(sds, aes(x = sd_log_abundance)) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(aes(x = rexp(nrow(sds), rate = 1/4)), 
               color = "red", linewidth = 1) +
  labs(title = "Distribution of metabolite specific standard deviations",
       subtitle = "Black: Observed SDs | Red: Prior Exp(1/4)")

## -----------------------------------------------------------------------------
# Simulate counts
counts <- data.frame(
  observation = 1:20, 
  counts = as.integer(rexp(20, 1e-7))
)

# Visualize: Black = Data, Red = Prior
ggplot(counts, aes(x = counts)) +
  geom_density(fill = "grey80", alpha = 0.5) +
  geom_density(data = data.frame(counts = rexp(1e5, 1/5e7)), 
               aes(x = counts), color = "red", linewidth = 1) +
  labs(title = "Distribution of cell counts",
       subtitle = "Black: Observed Counts | Red: Prior Exp(1/5e7)")

## ----eval = FALSE-------------------------------------------------------------
# # # Prepare data (ensure at least two time points per condition)
# # abundances$condition <- "A"
# # abundances <- rbind(abundances, abundances)
# # abundances$time <- rep(c(1, 2), each = nrow(abundances) / 2)
# #
# # counts$condition <- "A"
# # counts <- rbind(counts, counts)
# # counts$time <- rep(c(1, 2), each = nrow(counts) / 2)
# #
# # # Fit model with custom priors
# # fit <- fit_dynamics_model(
# #   model = "raw_plus_counts",
# #   model_option = "sd_per_condition",
# #   data = abundances,
# #   counts = counts,
# #   scaled_measurement = "abundance",
# #
# #   # Prior for mean of metabolite specific abundances (Normal)
# #   prior_mean_abundance = c(12, 10),
# #
# #   # Prior for metabolite specific standard deviation (Exponential mean)
# #   prior_sd_abundance = 4,
# #
# #   # Prior for cell counts (Exponential mean)
# #   prior_counts = 5e7
# # )