## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
library(pepVet)
data("aa_properties", package = "pepVet")

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

## ----paths--------------------------------------------------------------------
bsa_path <- system.file("extdata", "P02769.fasta", package = "pepVet")
h3_path <- system.file("extdata", "P68431.fasta", package = "pepVet")
proteome_path <- system.file(
  "extdata", "small_proteome_50_proteins.fasta",
  package = "pepVet"
)

## ----digest-basic-------------------------------------------------------------
digest_protein(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)

## ----digest-mc----------------------------------------------------------------
# mc=0: strict cleavage only
digest_protein("AKRTPK", enzyme = "trypsin", missed_cleavages = 0L)

# mc=1: also include once-skipped joins
digest_protein("AKRTPK", enzyme = "trypsin", missed_cleavages = 1L)

## ----cleavage-annotations-----------------------------------------------------
annotate_cleavage_sites(bsa_path, enzyme = "trypsin")

## ----score-basic--------------------------------------------------------------
digest_result <- digest_protein(bsa_path,
  enzyme = "trypsin",
  missed_cleavages = 1L
)
score_peptides(digest_result)

## ----preset-standard----------------------------------------------------------
pepvet_preset("standard")

## ----preset-evaluate----------------------------------------------------------
syn_path <- system.file("extdata", "P37840_isoforms.fasta", package = "pepVet")
syn_proteome <- digest_protein(syn_path, enzyme = "trypsin")
targeted_preset <- pepvet_preset("targeted")

do.call(
  evaluate_digest,
  c(list(sequence = syn_path, enzyme = "trypsin", proteome = syn_proteome), targeted_preset)
)$scores

## ----score-proteome, eval=FALSE-----------------------------------------------
# # Digest the proteome background first
# proteome_digest <- digest_protein(proteome_path, enzyme = "trypsin")
# 
# # Score BSA in the context of that proteome
# bsa_digest <- digest_protein(bsa_path, enzyme = "trypsin")
# score_peptides(bsa_digest, proteome = proteome_digest)

## ----evaluate-----------------------------------------------------------------
ev <- evaluate_digest(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)
names(ev)

## ----evaluate-scores----------------------------------------------------------
ev$scores

## ----evaluate-params----------------------------------------------------------
ev$params

## ----evaluate-cleavage-efficiency---------------------------------------------
ev_eff <- evaluate_digest(
  bsa_path,
  enzyme = "trypsin",
  missed_cleavages = 1L,
  include_cleavage_efficiency = TRUE
)

ev_eff$peptides
ev_eff$scores[, c("protein_id", "n_high_efficiency_sites", "n_low_efficiency_sites")]

## ----compare------------------------------------------------------------------
comp <- compare_digests(
  bsa_path,
  enzymes = c(
    "trypsin", "lysc", "glutamyl endopeptidase",
    "asp-n endopeptidase", "chymotrypsin-high"
  ),
  missed_cleavages = 1L
)
comp

## ----recommend----------------------------------------------------------------
recommend_enzyme(
  bsa_path,
  enzymes = c(
    "trypsin", "lysc", "glutamyl endopeptidase",
    "asp-n endopeptidase", "chymotrypsin-high"
  ),
  missed_cleavages = 1L
)

## ----batch, eval=FALSE--------------------------------------------------------
# batch <- batch_evaluate(proteome_path,
#   enzyme = "trypsin",
#   missed_cleavages = 1L
# )
# 
# # Number of proteins evaluated
# nrow(batch)
# 
# # Score and verdict for the first few proteins
# batch[, c("protein_id", "composite_score", "verdict")]

## ----proteome-onboarding, eval=FALSE------------------------------------------
# proteome_digest <- digest_protein(proteome_path, enzyme = "trypsin")
# batch <- batch_evaluate(proteome_path,
#   enzyme = "trypsin",
#   proteome = proteome_digest
# )

## ----summarize-batch, eval=FALSE----------------------------------------------
# summary <- summarize_batch(batch)
# 
# # Verdict distribution (Good / Moderate / Poor and their percentages)
# summary$verdict_counts
# 
# # Composite score distribution
# summary$score_distribution
# 
# # Per-component mean scores: lowest value is the weakest dimension
# summary$component_summary
# 
# # Proteins in the bottom 10% by composite score
# summary$problem_proteins
# 
# # Moderate / Poor proteins where hydrophobicity or short-protein flags are set:
# # likely candidates for switching to a less specific enzyme
# summary$enzyme_switch_candidates

## ----triage, eval=FALSE-------------------------------------------------------
# triaged <- triage_proteins(batch)
# 
# # Count of each action
# table(triaged$action)
# 
# # Proteins that should be tried with a different enzyme
# triaged[
#   triaged$action == "try_other_enzyme",
#   c("protein_id", "verdict", "composite_score")
# ]

## ----export-skyline-----------------------------------------------------------
peps <- digest_protein(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)

# Skyline transition list: one row per valid peptide per charge state
# Columns: Protein, Peptide Sequence, Precursor Charge, Precursor Mz
export_peptide_list(peps, format = "skyline", charges = 2:3)

## ----export-generic-----------------------------------------------------------
# Generic annotated table: all peptide columns plus gravy, pI, and valid flag
export_peptide_list(peps, format = "generic")

## ----export-fasta-------------------------------------------------------------
# FASTA character vector for valid peptides only
# Each header: >protein_id|start-end
head(export_peptide_list(peps, format = "fasta"))

## ----export-file, eval=FALSE--------------------------------------------------
# export_peptide_list(peps, format = "skyline", file = "bsa_transitions.csv")
# export_peptide_list(peps, format = "fasta", file = "bsa_peptides.fasta")

## ----pepvet-check-------------------------------------------------------------
result <- pepvet_check(bsa_path, enzyme = "trypsin")

## ----pepvet-check-pipe, eval=FALSE--------------------------------------------
# result <- pepvet_check(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)
# result$scores
# result$peptides

## ----report-single------------------------------------------------------------
digest_report(ev)

## ----report-compare-----------------------------------------------------------
digest_report(comp)

## ----pipe, eval=FALSE---------------------------------------------------------
# # Compare enzymes and print a report
# comp <- compare_digests(bsa_path,
#   enzymes = c("trypsin", "lysc", "glutamyl endopeptidase")
# )
# digest_report(comp)
# 
# # Get the top model rank directly from the sequence
# winner <- recommend_enzyme(bsa_path,
#   enzymes = c("trypsin", "lysc", "glutamyl endopeptidase")
# )

## ----h3-trypsin---------------------------------------------------------------
ev_h3_trypsin <- evaluate_digest(h3_path, enzyme = "trypsin")
ev_h3_trypsin$scores$verdict

## ----h3-compare---------------------------------------------------------------
compare_digests(
  h3_path,
  enzymes = c(
    "trypsin", "lysc", "glutamyl endopeptidase",
    "asp-n endopeptidase"
  )
)

## ----aa-properties------------------------------------------------------------
aa_properties

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