## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(pepVet)
has_ggplot2 <- requireNamespace("ggplot2", quietly = TRUE)

## ----data, include = FALSE----------------------------------------------------
base_dir <- system.file("extdata", "comparison-data", package = "pepVet")
overlap <- read.csv(file.path(base_dir, "sectionA-overlap.csv"))
caps <- read.csv(file.path(base_dir, "sectionB-capabilities.csv"))
scores <- read.csv(file.path(base_dir, "sectionC-scores-all.csv"))
presets <- read.csv(file.path(base_dir, "sectionC2-presets.csv"))
pr_data <- read.csv(file.path(base_dir, "sectionD-peptideranger.csv"))
pc_data <- read.csv(file.path(base_dir, "sectionE-protein-cleaver.csv"))

enzyme_labels <- c(
  "trypsin" = "Trypsin",
  "lysc" = "Lys-C",
  "chymotrypsin-high" = "Chymotrypsin",
  "glutamyl endopeptidase" = "Glu-C",
  "asp-n endopeptidase" = "Asp-N"
)
enzyme_fill <- c(
  "trypsin" = "#3B7A9E",
  "lysc" = "#6BA292",
  "chymotrypsin-high" = "#D4A76A",
  "glutamyl endopeptidase" = "#C46A6A",
  "asp-n endopeptidase" = "#8B7EB5"
)
verdict_fill <- c(
  "Good" = "#3A8C5F",
  "Moderate" = "#D4A76A",
  "Poor" = "#C46A6A"
)

## ----install, eval = FALSE----------------------------------------------------
# if (!requireNamespace("BiocManager", quietly = TRUE)) {
#   install.packages("BiocManager")
# }
# BiocManager::install("pepVet")

## ----overlap-table------------------------------------------------------------
knitr::kable(overlap, caption = "BSA + trypsin MC=1 peptide overlap")

## ----overlap-fig, eval=has_ggplot2, fig.width=8, fig.height=5, fig.cap="Peptide overlap among MS-Digest, ExPASy PeptideMass, and pepVet for BSA trypsin MC=1. The lighter segment marks peptides outside the 116 shared by all three tools.", fig.alt="Stacked bars showing 154 peptides from pepVet, 130 from MS-Digest, and 151 from ExPASy, with 116 common to all three."----
ggplot2::ggplot(overlap, ggplot2::aes(x = Tool)) +
  ggplot2::geom_col(
    ggplot2::aes(y = N_unique),
    fill = "#AFC6D5", color = "black", width = 0.65
  ) +
  ggplot2::geom_col(
    ggplot2::aes(y = Shared_all_three),
    fill = "#3B7A9E", color = "black", width = 0.65
  ) +
  ggplot2::geom_text(
    ggplot2::aes(y = N_unique, label = N_unique),
    vjust = -0.35
  ) +
  ggplot2::scale_y_continuous(
    limits = c(0, max(overlap$N_unique) * 1.12),
    expand = c(0, 0)
  ) +
  ggplot2::labs(
    title = "Peptide overlap: BSA trypsin MC=1",
    subtitle = "Dark bars show the 116 peptides shared by all three tools",
    x = NULL, y = "Number of peptides"
  ) +
  ggplot2::theme_minimal(base_size = 10) +
  ggplot2::theme(panel.grid.major.x = ggplot2::element_blank())

## ----cap-table----------------------------------------------------------------
knitr::kable(caps, caption = "Tool capability matrix")

## ----enzyme-fig, eval=has_ggplot2, fig.width=8, fig.height=5, fig.cap="Enzyme comparison across five proteins. Dashed lines mark the Moderate and Good verdict thresholds. Each bar is one enzyme. pepVet's composite score is a weighted sum of 5 components.", fig.alt="Grouped bar chart of composite scores by protein and enzyme with dashed verdict thresholds."----
ggplot2::ggplot(
  scores,
  ggplot2::aes(x = protein, y = composite, fill = enzyme)
) +
  ggplot2::geom_col(
    position = "dodge", color = "black", alpha = 0.85
  ) +
  ggplot2::geom_hline(
    yintercept = c(0.40, 0.65),
    linetype = "dashed", linewidth = 0.5,
    color = c("#D4A76A", "#3A8C5F")
  ) +
  ggplot2::scale_fill_manual(
    values = enzyme_fill,
    labels = enzyme_labels[names(enzyme_fill)]
  ) +
  ggplot2::scale_y_continuous(
    limits = c(0, 1), expand = c(0, 0.02)
  ) +
  ggplot2::labs(
    title = "Enzyme comparison across proteins",
    y = "Composite score", x = NULL, fill = "Enzyme"
  ) +
  ggplot2::theme_minimal(base_size = 11) +
  ggplot2::theme(legend.position = "bottom")

## ----preset-fig, eval=has_ggplot2, fig.width=10, fig.height=5, fig.cap="Workflow preset effects on length-valid peptide count. Grey bars show all theoretical digest products and coloured bars show peptides inside each preset's length range.", fig.alt="Faceted bar chart of total and length-valid peptide counts across 6 presets for BSA, H3, and BACE1."----
ggplot2::ggplot(presets, ggplot2::aes(x = preset)) +
  ggplot2::geom_col(
    ggplot2::aes(y = n_total),
    fill = "#D8DDE6", alpha = 0.7, width = 0.85
  ) +
  ggplot2::geom_col(
    ggplot2::aes(y = n_length_valid, fill = verdict),
    color = "black", alpha = 0.85, width = 0.85
  ) +
  ggplot2::geom_text(
    ggplot2::aes(y = n_length_valid, label = n_length_valid),
    vjust = -0.3, size = 3.2
  ) +
  ggplot2::facet_wrap(~protein, nrow = 1) +
  ggplot2::scale_fill_manual(values = verdict_fill) +
  ggplot2::labs(
    title = "Preset effects on length-valid peptide count",
    x = NULL, y = "Number of peptides", fill = "Verdict"
  ) +
  ggplot2::theme_minimal(base_size = 10) +
  ggplot2::theme(
    axis.text.x = ggplot2::element_text(
      angle = 35, hjust = 1, size = 8
    ),
    legend.position = "bottom",
    strip.background = ggplot2::element_rect(
      fill = "#F0F0F0", color = NA
    ),
    strip.text = ggplot2::element_text(face = "bold")
  )

## ----pr-fig, eval=has_ggplot2, fig.width=9, fig.height=6, fig.cap="PeptideRanger mean-score difference for peptides inside versus outside the selected length-and-GRAVY window. Positive values mean that the inside-window group has the higher mean score.", fig.alt="Horizontal bar chart of PeptideRanger mean-score differences for computable protein-enzyme combinations, colored by sign."----
pr_plot <- pr_data[!is.na(pr_data$mean_score_difference), ]
pr_plot$sign <- ifelse(
  pr_plot$mean_score_difference >= 0, "positive", "negative"
)
pr_plot$enzyme_label <- enzyme_labels[pr_plot$enzyme]
enzyme_order <- tapply(
  pr_plot$mean_score_difference, pr_plot$enzyme_label, mean, na.rm = TRUE
)
pr_plot$enzyme_label <- factor(
  pr_plot$enzyme_label, levels = names(sort(enzyme_order))
)

ggplot2::ggplot(
  pr_plot,
  ggplot2::aes(x = mean_score_difference, y = enzyme_label, fill = sign)
) +
  ggplot2::geom_col(color = "black", alpha = 0.85, width = 0.7) +
  ggplot2::geom_vline(
    xintercept = 0, linewidth = 0.5, color = "grey50"
  ) +
  ggplot2::scale_fill_manual(values = c(
    "positive" = "#3A8C5F", "negative" = "#C46A6A"
  )) +
  ggplot2::facet_grid(protein ~ ., scales = "free_y", space = "free_y") +
  ggplot2::labs(
    title = "PeptideRanger score difference by selected window",
    x = "Mean score difference", y = NULL, fill = NULL
  ) +
  ggplot2::theme_minimal(base_size = 10) +
  ggplot2::theme(
    legend.position = "none",
    panel.grid.major.y = ggplot2::element_blank(),
    strip.text.y = ggplot2::element_text(
      angle = 0, hjust = 1, face = "bold", size = 9
    ),
    strip.background = ggplot2::element_rect(
      fill = "#F0F0F0", color = NA
    )
  )

## ----pr-code, eval = requireNamespace("PeptideRanger", quietly = TRUE)--------
# preset <- pepvet_preset("standard")
# ev <- do.call(
#   evaluate_digest,
#   c(
#     list(
#       sequence = system.file("extdata", "P02769.fasta", package = "pepVet"),
#       enzyme = "trypsin",
#       missed_cleavages = 1L
#     ),
#     preset
#   )
# )
# length_valid <- ev$peptides
# length_valid <- length_valid[
#   length_valid$length >= preset$length_range[1] &
#     length_valid$length <= preset$length_range[2],
# ]
# pr_predict <- getExportedValue("PeptideRanger", "peptide_predictions")
# pr_model <- getExportedValue("PeptideRanger", "RFmodel_ProteomicsDB")
# pr_res <- pr_predict(
#   length_valid$peptide,
#   prediction_model = pr_model
# )
# head(pr_res)

## ----pc-table-----------------------------------------------------------------
knitr::kable(
  pc_data,
  caption = "Protein Cleaver identifiability vs selected pepVet window"
)

## ----pc-fig, eval=has_ggplot2, fig.width=8, fig.height=5, fig.cap="Overlap between the selected pepVet length-and-GRAVY window and the simulated Protein Cleaver identifiability rule. 'Both filters' passes both definitions.", fig.alt="Bar chart comparing peptides passing both filters with peptides passing only the simulated Protein Cleaver rule for BSA, H3, and BACE1."----
pc_plot <- rbind(
  data.frame(
    protein = pc_data$protein,
    classifier = "Both filters",
    count = pc_data$both_pass
  ),
  data.frame(
    protein = pc_data$protein,
    classifier = "PC identifiable only",
    count = pc_data$pc_only
  )
)

ggplot2::ggplot(
  pc_plot,
  ggplot2::aes(x = protein, y = count, fill = classifier)
) +
  ggplot2::geom_col(
    position = "dodge", color = "black", alpha = 0.85, width = 0.6
  ) +
  ggplot2::geom_text(
    ggplot2::aes(label = count),
    position = ggplot2::position_dodge(width = 0.6),
    vjust = -0.3, size = 3.5
  ) +
  ggplot2::scale_fill_manual(values = c(
    "Both filters" = "#3A8C5F",
    "PC identifiable only" = "#C46A6A"
  )) +
  ggplot2::labs(
    title = "pepVet vs Protein Cleaver",
    x = NULL, y = "Number of peptides", fill = "Classification"
  ) +
  ggplot2::theme_minimal(base_size = 10) +
  ggplot2::theme(legend.position = "bottom")

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