---
title: "Getting Started with pepVet"
description: >
  A walkthrough of pepVet's core pipeline from a single protein digest to multi-enzyme comparison, batch evaluation, and peptide export.
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{Getting Started with pepVet}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

pepVet evaluates proteolytic digests before acquisition. A protein sequence and enzyme produce peptide coordinates, score components, a composite label, and the resolved scoring parameters. The [tool comparison](pepVet-tool-comparison.html) describes differences between pepVet and five other digestion tools.

## The core pipeline

```
Protein sequence
      |
      v
    digest_protein()    ->  tibble of peptides (one row per peptide)
      |
      v
    score_peptides()    ->  one-row tibble of component scores + verdict
      |
      v
    evaluate_digest()   ->  convenience wrapper: digest + score + params in one call
      |
      v
    pepvet_check()      ->  evaluate + print report in one interactive call
      |
      v
    compare_digests()   ->  run evaluate_digest across N enzymes, rank by composite
      |
      v
    recommend_enzyme()  ->  return the name of the top-ranked enzyme
      |
      v
    batch_evaluate()    ->  run evaluate_digest across all proteins in a FASTA
      |                     returns a flat tibble: one row per protein
      v
    summarize_batch()   ->  proteome-level verdict distribution and score stats
      |
      v
    triage_proteins()   ->  per-protein action labels
      |
      v
    export_peptide_list()  ->  export valid peptides to Skyline, generic CSV, or FASTA
      |
      v
    digest_report()     ->  compact console output for any of the above
```

Every output is programmatic. `digest_protein()` returns a tibble with one row per peptide. `score_peptides()` returns a tibble with one row per protein. `compare_digests()` returns a tibble with one row per enzyme. `evaluate_digest()` returns a named list with the score tibble, the peptide tibble, and the resolved parameters.

## Protein input

Every pepVet function that takes a protein sequence accepts the same input shapes. You can pass:

- a raw character sequence
- a named character vector (one entry per protein)
- a `Biostrings::AAString` or `AAStringSet`
- a path to a FASTA file (single- or multi-entry)

The package ships fixed FASTA fixtures in `inst/extdata/` for these examples.

```{r 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"
)
```

## Step 1: Digest a protein

`digest_protein()` cleaves a protein sequence using any cleaver-compatible enzyme rule and returns a tibble with one row per peptide fragment.

```{r digest-basic}
digest_protein(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)
```

The output columns are:

| Column | Description |
|---|---|
| `protein_id` | Full FASTA header or supplied name |
| `peptide` | Amino acid sequence of the fragment |
| `start` | 1-based start position in the parent protein |
| `end` | 1-based end position (inclusive) |
| `length` | Peptide length in residues |
| `missed_cleavages` | Number of internal cleavage sites skipped |

`digest_protein()` expands missed cleavages explicitly and retains exact start and end coordinates for each joined product.

```{r 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)
```

The validation layer checks enzyme names against a registry of 40 cleaver-compatible rules, normalises case, and trims whitespace. Common examples include `trypsin`, `lysc`, `glutamyl endopeptidase`, `asp-n endopeptidase`, `chymotrypsin-high`, and `thermolysin`. The input `" Trypsin "` therefore resolves to trypsin.

### Step 1a: Inspect cleavage efficiency

`annotate_cleavage_sites()` classifies local trypsin-family motifs using the package's sequence rules.

```{r cleavage-annotations}
annotate_cleavage_sites(bsa_path, enzyme = "trypsin")
```

pepVet currently labels tryptic sites as:

- `low` for proline-blocked (`KP` / `RP`) and acidic P1' (`KD`, `KE`, `RD`, `RE`) motifs
- `medium` for adjacent-basic sites such as `KK`, `KR`, and `RR`
- `high` for other tryptic sites

These labels are sequence-local annotations. They do not model structure, PTMs, extended subsite preferences, cleavage probability, or abundance.

## Step 2: Score the peptide set

`score_peptides()` summarises a digest tibble into per-protein component scores and a weighted composite verdict.

```{r score-basic}
digest_result <- digest_protein(bsa_path,
  enzyme = "trypsin",
  missed_cleavages = 1L
)
score_peptides(digest_result)
```

Only peptides inside the active `length_range` count as valid. The default window is 7 to 25 residues. Fragments outside that window still appear in the digest tibble, but they do not contribute to the component numerators. See the [scoring model](scoring-model.html) article for the formula of each component.

Composite score verdict thresholds are `Good` at `>= 0.65`, `Moderate` at `>= 0.40`, and `Poor` below `0.40`. These strict, conservative boundaries are package design choices for triage. They are not calibrated probabilities, and the current PeptideAtlas analysis did not reliably calibrate them. Treat the labels as prompts for review, not predictions of experimental success.

### Scoring model

pepVet is not a peptide detectability predictor. Existing ML-based tools learn peptide-level detection patterns from experimental MS data. pepVet does something narrower: it ranks the digest quality of a protein for a chosen enzyme-workflow combination using explicit physicochemical criteria.

- ML detectability tools learn nonlinear feature interactions from observed detections.
- pepVet uses a rule-based multi-criteria decision analysis with linear weights.
- ML detectability tools output per-peptide predictions or suitability scores, depending on the model.
- pepVet outputs per-protein digest-quality rankings.

pepVet scores describe which digest-level components raise or lower a ranking. They do not estimate the observation probability of a specific peptide.

The [scoring model](scoring-model.html) separates package design choices, expert priors, physical constants, and the limited PeptideAtlas evaluation. The [workflow presets](workflow-presets.html) article lists the resolved settings for each preset.

### Known limitations

- Verdict thresholds are heuristic labels, not calibrated classifications.
- Weight vectors are expert priors, not empirically trained coefficients.
- The model omits PTMs, chemical labels, and digestion chemistry beyond the current rules.
- GRAVY and peptide-length windows assume conventional C18 reversed-phase LC with ESI.
- `S_charge` reflects extra internal basic-residue richness, not whether a peptide can ionize at all.
- Proteome-aware uniqueness is only as relevant as the proteome digest you supply.

### Scope

pepVet scores are interpretable rankings, not calibrated probabilities. The composite score has no unit. It ranks proteins by digest quality within a given enzyme-workflow combination.

Cross-workflow comparisons are not valid when the resolved scoring configuration changes. `score_peptides()` records a `preset_used` column. The function reports a preset name when the resolved GRAVY window, peptide-length window, weights, and pI setting exactly match a shipped preset. It marks other configurations as `"custom"`.

### Step 2a: Apply a workflow preset

Workflow presets give you a starting parameter set for common experiment types. Each preset returns a list with `gravy_range`, `length_range`, and `weights`, so you can pass the result into `evaluate_digest()` or `score_peptides()` directly.

```{r preset-standard}
pepvet_preset("standard")
```

```{r 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
```

Six presets ship with pepVet. Each bundles a length range, a GRAVY window, and scoring weights for a common experiment type. See the [workflow presets](workflow-presets.html) article for the full reference.

Presets with non-zero `S_unique` weights require a proteome digest. pepVet rejects those presets when `proteome` is missing because the component is undefined without a comparison background.

### Step 2b: Add proteome-aware uniqueness

When you supply a `proteome` digest, `score_peptides()` adds a sixth component, `S_unique`. It measures the fraction of valid peptides that appear in exactly one protein in the supplied background. Shared sequences can map to more than one protein, so the relevance of this component depends on the background.

```{r 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)
```

Proteome-aware scoring switches the default weight set to give `S_unique` a 20% share: `c(S_length = 0.160, S_coverage = 0.279, S_count = 0.181, S_hydro = 0.110, S_charge = 0.070, S_unique = 0.200)`.

## Step 3: Evaluate in one call

`evaluate_digest()` wraps `digest_protein()` and `score_peptides()` into a single call and returns a named list.

```{r evaluate}
ev <- evaluate_digest(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)
names(ev)
```

```{r evaluate-scores}
ev$scores
```

```{r evaluate-params}
ev$params
```

The `$params` element records the resolved enzyme name, missed cleavages, full protein IDs, and `preset_used` value.

When you need peptide-level cleavage-risk context, request it explicitly. The flag adds a `cleavage_efficiency` column to the peptide table. The protein-level efficiency counts (`n_high_efficiency_sites`, `n_low_efficiency_sites`) are always present in the scores table regardless of this flag.

```{r 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")]
```

The protein-level counts are informational and do not contribute to the score. They summarize the sequence-local labels returned by the annotation rules.

## Step 4: Compare enzymes

`compare_digests()` runs `evaluate_digest()` across a vector of enzyme names for a single protein and returns a tibble sorted from best to worst composite score.

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

`compare_digests()` defines one row per enzyme for one protein and therefore rejects multi-protein input.

## Step 5: Get the top model rank

`recommend_enzyme()` runs `compare_digests()` internally and returns the name of the highest-scoring enzyme as a character string. If multiple enzymes tie on the composite score, the function returns all tied names in alphabetical order. Despite the function's established name, this is a model ranking under the supplied settings, not an experimental recommendation.

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

## Step 6: Batch across a proteome

`batch_evaluate()` runs `evaluate_digest()` independently for every protein in a multi-FASTA file and returns a **flat tibble** with one row per protein. Columns include `protein_id`, `protein_length`, all component scores, `composite_score`, `verdict`, count fields, and four sequence-level difficulty flags.

```{r 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")]
```

`batch_evaluate()` evaluates each protein independently. Valid-peptide counts and hydrophobicity flags follow the active `length_range` and `gravy_range`. Short-protein and low-complexity flags remain fixed sequence-level heuristics. The returned tibble serves as input to `summarize_batch()` and `triage_proteins()`.

For proteome-aware uniqueness scoring, supply a background proteome digest:

```{r proteome-onboarding, eval=FALSE}
proteome_digest <- digest_protein(proteome_path, enzyme = "trypsin")
batch <- batch_evaluate(proteome_path,
  enzyme = "trypsin",
  proteome = proteome_digest
)
```

### Step 6a: Summarize a batch

`summarize_batch()` computes proteome-level aggregate statistics from the tibble returned by `batch_evaluate()`. It returns a named list with five elements.

```{r 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
```

`enzyme_switch_candidates` is a heuristic derived from sequence-level difficulty flags, not from running alternative enzymes. Use `compare_digests()` to assess a specific alternative enzyme under the active scoring settings.

### Step 6b: Triage proteins

`triage_proteins()` appends a deterministic `action` label derived from the verdict and difficulty flags.

```{r 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")
]
```

The action rules are:

- `proceed`: the verdict is `Good`.
- `skip`: a non-Good row has no valid peptides or has the low-complexity flag.
- `try_other_enzyme`: a remaining row has the hydrophobicity flag, the short-protein flag, or a `Poor` verdict.
- `consider_alternative`: any remaining `Moderate` row.

These labels do not run another enzyme or predict that the suggested action will improve an experiment.

### Step 6c: Export a peptide list

`export_peptide_list()` filters valid peptides from any `digest_protein()` tibble and writes them in a format compatible with downstream tools.

```{r 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)
```

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

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

Write directly to a file by passing a path to the `file` argument:

```{r export-file, eval=FALSE}
export_peptide_list(peps, format = "skyline", file = "bsa_transitions.csv")
export_peptide_list(peps, format = "fasta", file = "bsa_peptides.fasta")
```

## Step 7: Quick interactive check

`pepvet_check()` combines `evaluate_digest()` and `digest_report()` in one interactive call. It prints a console report and returns the result invisibly.

```{r pepvet-check}
result <- pepvet_check(bsa_path, enzyme = "trypsin")
```

The returned value is the full `evaluate_digest()` result list, so you can pipe it into any downstream step:

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

## Step 8: Report to the console

`digest_report()` prints an ASCII-safe summary for an `evaluate_digest()` result or a ranked table for a `compare_digests()` result. Long protein identifiers wrap, and comparison tables simplify when the terminal is narrow. The function returns its input invisibly.

```{r report-single}
digest_report(ev)
```

```{r report-compare}
digest_report(comp)
```

`digest_report()` returns `invisible(x)`, so a script can print the report and retain the object:

```{r 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")
)
```

## Additional examples and references

The remaining sections connect the core workflow to a challenging protein, the shipped amino acid data, and the package's narrower role among proteomics tools.

### A challenging protein: Histone H3.1

Histone H3.1 is a package example in which strict tryptic digestion produces many short fragments from its lysine- and arginine-rich N-terminal tail. Many of those fragments fall below pepVet's conservative 7-residue boundary. This is a scoring result under the active settings, not a universal detection threshold.

```{r h3-trypsin}
ev_h3_trypsin <- evaluate_digest(h3_path, enzyme = "trypsin")
ev_h3_trypsin$scores$verdict
```

In this in-silico comparison, LysC recognises only lysine and cleaves less frequently on H3.1, producing longer peptides with a higher score under the default settings:

```{r h3-compare}
compare_digests(
  h3_path,
  enzymes = c(
    "trypsin", "lysc", "glutamyl endopeptidase",
    "asp-n endopeptidase"
  )
)
```

### Amino acid reference data

The package ships `aa_properties`, a 22-row tibble containing the 20 standard amino acids plus selenocysteine and pyrrolysine.

```{r aa-properties}
aa_properties
```

The `kyte_doolittle` column contains the Kyte-Doolittle (1982) hydrophobicity scale used in `S_hydro`. The `is_basic` flag supports `S_charge`, and the `pKa_side_chain` column holds reference pKa values for the seven ionisable residues (C, D, E, H, K, R, Y) used in `calculate_pI()`.

### How pepVet differs from common tools

pepVet covers a narrow part of the workflow. It is not a spectral library builder or a peptide detectability predictor. It compares enzyme-workflow combinations using explicit digest-level criteria before acquisition. This describes the package scope rather than an experimental performance claim.

The [tool-comparison vignette](pepVet-tool-comparison.html) compares pepVet against five tools on a panel of five proteins and five enzymes. It covers baseline peptide overlap, capability differences across 16 dimensions, the scoring model in practice, workflow preset effects, pipeline integration with PeptideRanger, and a binary versus graded classification comparison with Protein Cleaver.

### Related articles

- [Scoring model](scoring-model.html): formulas, weights, thresholds, and limitations.
- [Enzyme selection](enzyme-selection.html): worked comparisons across several proteins.
- [Workflow presets](workflow-presets.html): resolved settings for the six shipped presets.
- [Visualisation](visualisation.html): plotting functions with examples.
- [Weight sensitivity](weight-sensitivity.html): Dirichlet Monte Carlo perturbation of the weights.
- [Score diagnostics](score-diagnostics.html): VIF, PCA, and ablation summaries.
- [Tool comparison](pepVet-tool-comparison.html): descriptive comparison with five other tools.

## Session info

```{r session-info}
sessionInfo()
```