## ----include = FALSE----------------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width = 7,
  fig.height = 5
)

## ----setup, message = FALSE---------------------------------------------------
library(pepdiff)
library(dplyr)

## ----generate-data------------------------------------------------------------
set.seed(42)

n_peptides <- 30
n_reps <- 10  # Good power with 10 reps

peptides <- paste0("PEP_", sprintf("%03d", 1:n_peptides))
genes <- paste0("GENE_", LETTERS[((1:n_peptides - 1) %% 26) + 1])

# 10 peptides with real effects (3-fold up or down)
diff_peptides <- peptides[1:10]
effect_directions <- c(rep(3, 5), rep(0.33, 5))

sim_data <- expand.grid(
  peptide = peptides,
  treatment = c("ctrl", "trt"),
  bio_rep = 1:n_reps,
  stringsAsFactors = FALSE
) %>%
  mutate(
    gene_id = genes[match(peptide, peptides)],
    base = rep(rgamma(n_peptides, shape = 5, rate = 0.5), each = 2 * n_reps),
    effect = case_when(
      peptide %in% diff_peptides & treatment == "trt" ~
        effect_directions[match(peptide, diff_peptides)],
      TRUE ~ 1
    ),
    value = rgamma(n(), shape = 15, rate = 15 / (base * effect))
  ) %>%
  select(peptide, gene_id, treatment, bio_rep, value)

# Import
temp_file <- tempfile(fileext = ".csv")
write.csv(sim_data, temp_file, row.names = FALSE)

dat <- read_pepdiff(
  temp_file,
  id = "peptide",
  gene = "gene_id",
  value = "value",
  factors = "treatment",
  replicate = "bio_rep"
)

dat

## ----run-all-tests, message = FALSE, results = 'hide'-------------------------
res_wilcox <- compare(dat, compare = "treatment", ref = "ctrl",
                       method = "pairwise", test = "wilcoxon")
res_boot <- compare(dat, compare = "treatment", ref = "ctrl",
                     method = "pairwise", test = "bootstrap_t")
res_bayes <- compare(dat, compare = "treatment", ref = "ctrl",
                      method = "pairwise", test = "bayes_t")
res_rp <- compare(dat, compare = "treatment", ref = "ctrl",
                   method = "pairwise", test = "rankprod")

## ----wilcoxon-----------------------------------------------------------------
res_wilcox

## ----bootstrap----------------------------------------------------------------
res_boot

## ----bayes--------------------------------------------------------------------
res_bayes

## ----bayes-strict, eval = FALSE-----------------------------------------------
# # Require strong evidence (BF > 10)
# res_bayes_strict <- compare(dat, compare = "treatment", ref = "ctrl",
#                              method = "pairwise", test = "bayes_t",
#                              bf_threshold = 10)

## ----rankprod-----------------------------------------------------------------
res_rp

## ----compare-results----------------------------------------------------------
comparison <- tibble(
  peptide = res_wilcox$results$peptide,
  truly_diff = peptide %in% diff_peptides,
  wilcoxon = res_wilcox$results$significant,
  bootstrap = res_boot$results$significant,
  bayes = res_bayes$results$significant,
  rankprod = res_rp$results$significant
)

# How many significant calls per test?
comparison %>%
  summarise(
    across(wilcoxon:rankprod, sum)
  )

## ----agreement----------------------------------------------------------------
# Agreement on truly differential peptides
comparison %>%
  filter(truly_diff) %>%
  summarise(
    detected_by_all = sum(wilcoxon & bootstrap & bayes & rankprod),
    detected_by_any = sum(wilcoxon | bootstrap | bayes | rankprod),
    missed_by_all = sum(!wilcoxon & !bootstrap & !bayes & !rankprod)
  )

## ----fdr-example--------------------------------------------------------------
# Raw p-values vs FDR-adjusted
res_wilcox$results %>%
  select(peptide, p_value, fdr, significant) %>%
  arrange(p_value) %>%
  head(10)

## ----plot-results, fig.height = 6---------------------------------------------
plot(res_wilcox)

## ----plot-bayes, fig.height = 6-----------------------------------------------
plot(res_bayes)

## ----cleanup, include = FALSE-------------------------------------------------
unlink(temp_file)