## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(
  echo = TRUE,
  warning = FALSE,
  message = FALSE
)

## ----installation, eval=FALSE-------------------------------------------------
# if (!requireNamespace("remotes", quietly = TRUE)) {
#   install.packages("remotes")
# }
# 
# remotes::install_github("MohammadDarbalaei/barmixR")

## ----load-package-------------------------------------------------------------
library(barmixR)

## -----------------------------------------------------------------------------
library(MGLM)
library(tidyverse)


set.seed(1001)

# --- Reduced size for vignette ---
k <- 4              # number of clones (cell lines)
n_reps <- 3         # fewer replicates
n_reads <- 2000     # fewer sequencing reads

# Only control + one treatment
conditions <- c("DMSO", "TreatA")

# Helper function
generate_profile <- function(pattern, k){
  breaks <- floor(seq(0, k, length.out = length(pattern) + 1))
  prof <- numeric(k)
  for(i in seq_along(pattern)){
    idx <- (breaks[i] + 1):breaks[i + 1]
    prof[idx] <- pattern[i]
  }
  prof + runif(k, 0, 0.05)
}

# Define profiles (simplified)
alpha_profiles <- list(
  DMSO   = rep(1/k, k),
  TreatA = generate_profile(c(2, 0.5), k)
)

# Simulate counts
all_data <- list()

for(cond in conditions){

  alpha <- alpha_profiles[[cond]]
  alpha_scaled <- alpha / sum(alpha) * 20  # smaller concentration

  mat <- rdirmn(n_reps, n_reads, alpha_scaled)

  rownames(mat) <- paste0(cond, "_rep", 1:n_reps)
  colnames(mat) <- paste0("clone", 1:k)

  all_data[[cond]] <- mat
}

data_count <- do.call(rbind, all_data)

## -----------------------------------------------------------------------------
print(data_count[1:2,1:2])

## -----------------------------------------------------------------------------
get_lognorm_params <- function(mean, sd){
  sigma2 <- log(1 + (sd / mean)^2)
  mu <- log(mean) - sigma2 / 2
  sigma <- sqrt(sigma2)
  list(meanlog = mu, sdlog = sigma)
}

params_t0 <- get_lognorm_params(100, 40)
params_t14 <- get_lognorm_params(300, 150)

sim_vols <- function(params, seed){
  set.seed(seed)
  rlnorm(n_reps, params$meanlog, params$sdlog)
}

DMSO_t0  <- sim_vols(params_t0, 1)
DMSO_t14 <- sim_vols(params_t14, 101)

TreatA_t14 <- sim_vols(get_lognorm_params(200, 120), 102)

# Combine
V_t0  <- round(c(DMSO_t0, DMSO_t0), 1)
V_t14 <- round(c(DMSO_t14, TreatA_t14), 1)

## -----------------------------------------------------------------------------
condition_count <- factor(
  sub("_rep[0-9]+$","",rownames(data_count)),
  levels=conditions
)

cell_line <- factor(paste0("clone",1:k))

data <- list(
  data_count=data_count,
  condition_count=condition_count,
  V=V_t14,
  VT0=V_t0,
  condition_v=condition_count,
  cell_line=cell_line
)

time_d <- 14

## -----------------------------------------------------------------------------
print(data)

## ----model_fit, cache=TRUE----------------------------------------------------
# Reduced number of posterior sampling iterations
# to keep the vignette computationally lightweight.
# For real analyses, increase iterations and chains
# to ensure proper convergence and reliable inference.

fit <- barmixRQTR(
  data = data,
  time_d = time_d,
  control = list(
    chains = 1,
    iter_count = 100,
    iter_V = 100,
    cores = 1
  ),
  in_vivo = TRUE
)

model <- fit

## ----fig.width=14, fig.height=8, out.width="100%"-----------------------------
ppc_barcode <- ppcBarcodes(model)
print(ppc_barcode$ppc_plot)

## ----fig.width=18, fig.height=8, out.width="100%"-----------------------------
sampled_fraction <- ppc_barcode$sampled_fraction

fraction_ratio_results <- fractionRatio(model, sampled_fraction)

fraction_ratio_results$plot_ratio_fraction

## ----fig.width=18, fig.height=8, out.width="100%"-----------------------------
ppc_population_results <- ppcPopulation(model)

ppc_population_results$ppc_plot

## -----------------------------------------------------------------------------
sampled_population <- ppc_population_results$sampled_population

## ----fig.width=18, fig.height=8, out.width="100%"-----------------------------
population_ratio_results <- populationRatio(model, sampled_population)

population_ratio_results$plot_ratio_population

## ----fig.width=18, fig.height=7, out.width="100%"-----------------------------
resistance <- QTRresistance(
  model,
  fraction_ratio_results$li_sam_ratio_relative,
  population_ratio_results$li_sam_ratio_V
)

print(resistance$treatment_resistance)

## ----fig.width=8, fig.height=6, out.width="100%"------------------------------
heatmap <- QTRheatmap(
  model,
  resistance$summary_table
)

print(heatmap$heatmap_within)

## ----fig.width=6, fig.height=6, out.width="100%"------------------------------
print(heatmap$heatmap_treatment)

## ----fig.width=18, fig.height=3, out.width="100%"-----------------------------
cell_order <- paste0("clone", 1:4)

decision <- QTRDecision(
  model = model,
  summary_table = resistance$summary_table,
  cell_order = cell_order
  # legend_text = "DMSO = control condition"
)

print(decision$rank_plot)


## ----session-info-------------------------------------------------------------
sessionInfo()