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

## ----toy----------------------------------------------------------------------
library(diagFDR)

set.seed(2)

n <- 6000
toy <- data.frame(
  id = paste0("run1||scan", seq_len(n)),
  is_decoy = sample(c(FALSE, TRUE), n, replace = TRUE, prob = c(0.98, 0.02)),
  q = pmin(1, runif(n)^3),                                  # many small q-values
  pep = pmin(1, pmax(0, runif(n)^3 + rnorm(n, sd = 0.02))), # correlated toy PEP
  run = "run1",
  score = rnorm(n)
)

x <- as_dfdr_tbl(
  toy,
  unit = "psm",
  scope = "global",
  q_source = "toy mokapot",
  q_max_export = 1
)

diag <- dfdr_run_all(
  xs = list(mokapot = x),
  alpha_main = 0.01,
  alphas = c(5e-4, 1e-3, 2e-3, 5e-3, 1e-2, 2e-2, 5e-2, 1e-1, 2e-1),
  low_conf = c(0.2, 0.5)
)

## ----headline-----------------------------------------------------------------
diag$tables$headline

## ----stability-plots----------------------------------------------------------
diag$plots$dalpha
diag$plots$cv

## ----boundary-stability-------------------------------------------------------
diag$plots$dwin
diag$plots$elasticity

## ----equal-chance-------------------------------------------------------------
diag$tables$equal_chance_pooled
diag$plots$equal_chance__mokapot

## ----pep----------------------------------------------------------------------
diag$tables$pep_IPE
diag$plots$pep_reliability__mokapot

## ----sumpep-------------------------------------------------------------------
# sumpep is a list-of-tibbles (one per list)
diag$tables$sumpep$mokapot

## ----real-mokapot, eval=FALSE-------------------------------------------------
# library(diagFDR)
# 
# # Read mokapot outputs (targets + decoys)
# raw <- read_mokapot_psms(
#   target_path = "path/to/your_output.mokapot.psms.txt",
#   decoy_path  = "path/to/your_output.mokapot.decoy.psms.txt"
# )
# 
# # Construct competed winners (1 per run+spectrum_id; max mokapot score)
# x <- mokapot_competed_universe(
#   raw,
#   id_mode = "runxid",
#   unit = "psm",
#   scope = "global",
#   q_source = "mokapot q-value",
#   q_max_export = 1
# )
# 
# # Run diagnostics
# diag <- dfdr_run_all(
#   xs = list(mokapot = x),
#   alpha_main = 0.01,
#   alphas = c(5e-4, 1e-3, 2e-3, 5e-3, 1e-2, 2e-2, 5e-2, 1e-1, 2e-1),
#   low_conf = c(0.2, 0.5)
# )
# 
# # Write outputs to disk (tables + plots + summary + manifest + README)
# dfdr_write_report(
#   diag,
#   out_dir = "diagFDR_mokapot_out",
#   formats = c("csv", "png", "manifest", "readme", "summary")
# )
# 
# # Render a single HTML report (requires rmarkdown in Suggests)
# dfdr_render_report(diag, out_dir = "diagFDR_mokapot_out")