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

## -----------------------------------------------------------------------------
library(ggplot2)
library(dplyr)
library(MetaboDynamics)
mu_mean <- 12      # Mean of the prior for metabolite mean
sd_mean <- 5       # SD of the prior for metabolite mean
prior_mean_abundance <- c(mu_mean,sd_mean)
prior_sd_abundance <- 2  # Rate parameter for exponential prior on SD

# Step 1: Sample metabolite-specific parameters from priors
n_metabolites <- 300
mu <- rnorm(n_metabolites, mean = prior_mean_abundance[1], sd = prior_mean_abundance[2])        # μ ~ N(12, 5)
sigma <- rexp(n_metabolites, rate = 1/prior_sd_abundance)      # σ ~ Exp(1/2), so mean = 2

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

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

abundances[abundances$metabolite==i,]$abundance <- abundance
}

abundances$abundance <- exp(abundances$abundance)

ggplot(abundances,aes(x=log(abundance)))+
  geom_density()+
  ggtitle("Distribution of log-transformed metabolite abundances")

## -----------------------------------------------------------------------------
ggplot(abundances,aes(x=log(abundance)))+
  geom_density()+
  geom_density(aes(x=rnorm(nrow(abundances),mean = 12,sd = 10)),col="red")+
  ggtitle("Distribution of log-transformed metabolite abundances",
          "red line = N(12,10)")

## -----------------------------------------------------------------------------
ggplot(abundances,aes(x=log(abundance),col=as.factor(metabolite)))+
  geom_density()+
  guides(col="none")+
  ggtitle("Distribution of log-transformed metabolite abundances","colour = metabolite")

## -----------------------------------------------------------------------------
sds <- abundances%>%group_by(metabolite)%>%
  mutate(sd_metabolite=sd(log(abundance)))%>%select(metabolite,sd_metabolite)%>%distinct()

ggplot(sds,aes(x=sd_metabolite))+
  geom_density()+
  ggtitle("Distribution of metabolite specific standard deviations")
  

## -----------------------------------------------------------------------------
ggplot(sds,aes(x=sd_metabolite))+
  geom_density()+
  geom_density(aes(x=rexp(nrow(sds),1/4)),col="red")+
  ggtitle("Distribution of metabolite specific standard deviations","red line: prior for standard deviation")

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

ggplot(counts,aes(x=counts))+
  geom_density()+  # visualize distribution
  ggtitle("Distribution of observed cell counts")

## -----------------------------------------------------------------------------
ggplot(counts)+
  geom_density(aes(x=counts))+ # observed data
  geom_density(data=as.data.frame(cbind(observations=1:1e5,counts=rexp(1e5,1/1e7))),
               aes(x=counts),col="red")+
  ggtitle("Distribution of cell counts","red line = prior for cell counts")

## -----------------------------------------------------------------------------
ggplot(counts)+
  geom_density(aes(x=counts))+ # observed data
  geom_density(data=as.data.frame(cbind(observations=1:1e5,counts=rexp(1e5,1/5e7))),
               aes(x=counts),col="red")+
  ggtitle("Distribution of cell counts","red line = prior for cell counts")

## -----------------------------------------------------------------------------
# we need at least two time points in one condition
abundances$condition <- "A"
abundances <- rbind(abundances,abundances)
abundances$time <- rep(c(1,2),each=nrow(abundances)/2)

counts$condition <- "A"
coutns <- rbind(counts,counts)
counts$time <- rep(c(1,2),each=nrow(counts)/2)


# out commented to reduce run-time of vignettes
# 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
#     prior_mean_abundance = c(12,10), # N(12,10)
# 
#     # prior for metabolite specific standard deviation
#     prior_sd_abundance = 4, # exponential(1/4)
# 
#     # prior for cell counts
#     prior_counts = 5e7, # exponential(1/5e7),
#   )