---
title: "pepVet Compared to Other Tools"
description: >
  A descriptive comparison of pepVet, MS-Digest, ExPASy PeptideMass, Protein Cleaver, ProteaseGuru, and PeptideRanger.
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
vignette: >
  %\VignetteIndexEntry{pepVet Compared to Other Tools}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

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

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

Six tools in this comparison simulate or evaluate proteolytic digests. The analysis records peptide overlap for one BSA trypsin configuration, then describes scoring, filtering, ranking, structural context, and programmatic access.

This vignette is package documentation, not paper evidence, an experimental validation, or a performance benchmark. It compares retained outputs and documented interfaces. The local artifacts do not retain service versions or retrieval dates for MS-Digest and PeptideMass, and the assessed Protein Cleaver source revision was not recorded. The capability entries therefore do not certify current-release equivalence.

## Comparison basis

The overlap analysis records agreement and disagreement among three retained outputs before the capability review separates the tools by purpose and interface.

### Common ground

The first analysis compares peptide sets from MS-Digest, ExPASy PeptideMass, and pepVet for BSA with trypsin and one missed cleavage.

```{r overlap-table}
knitr::kable(overlap, caption = "BSA + trypsin MC=1 peptide overlap")
```

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

The three outputs share 116 peptides. pepVet returns 154 because it reports all fragments from the cleavage engine before scoring filters. The assessed MS-Digest and ExPASy outputs contain 130 and 151 peptides. On peptides shared by each pair, the mass correlations are 0.999983 for MS-Digest and ExPASy, 0.999821 for pepVet and MS-Digest, and 1.000000 for pepVet and ExPASy. These values describe the assessed configuration and do not establish equivalence across other inputs or tool versions.

### How tools differ

The capability matrix records the assessed interfaces and outputs.

```{r cap-table}
knitr::kable(caps, caption = "Tool capability matrix")
```

MS-Digest (Baker, Clauser, and Burlingame, UCSF) returns peptide masses and coordinates for a selected digest rule. Its [official instructions](https://prospector.ucsf.edu/prospector/html/instruct/digestman.htm) describe multiple-protein input and hydrophobicity indices, and the [automation documentation](https://prospector.ucsf.edu/prospector/html/misc/invar.htm) describes command-line and XML interfaces. The retained output does not establish which service version produced it. MS-Digest does not supply a pepVet-style composite digest score or score-band label.

PeptideMass (Wilkins et al. 1997) accepts a UniProt accession and returns neutral average peptide masses together with the source protein's pI. It can use signal-peptide and PTM annotations from UniProt. It does not score or rank digests.

Protein Cleaver (Koulouras and Xu 2025) is an R Shiny application built with `cleaver`. It flags peptides as identifiable or non-identifiable using configurable length and mass windows, ranks proteases using maximum coverage and identifiable-peptide count, and includes sequence and 3D structure mapping. Its [official paper](https://doi.org/10.3389/fbinf.2025.1576317) and [repository documentation](https://github.com/gkoulouras/ProteinCleaver), checked 2026-07-14, do not document retention-time prediction. The unrecorded source revision assessed for this comparison used the same cleavage engine as pepVet and a different mass table, including 71.079 Da for alanine compared with pepVet's 71.037 Da. The assessed revision had no GRAVY filter, composite score, verdict, or workflow presets.

ProteaseGuru (Miller et al. 2021) is a desktop application for comparing multiple proteases. It reports SSRCalc hydrophobicity, coverage maps, and shared or exclusive peptides. The assessed documentation did not provide a pepVet-style composite score, workflow presets, or an R API.

PeptideRanger (Riley et al. 2023) is an R package that uses a random forest to prioritize peptides for MS analysis. It produces an RF score from 0 to 1 per peptide. In this comparison, pepVet summarizes digest-level criteria and PeptideRanger supplies a peptide-level suitability score. We do not treat either score as a calibrated probability.

## Worked comparisons

The following analyses compare digest-level scores, workflow settings, and peptide-level classifications without treating the tools as interchangeable.

### pepVet scoring in practice

The benchmark contains five proteins and five enzymes. The pepVet composite combines five components: length, sequence coverage, peptide count, hydrophobicity, and charge richness.

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

Highest composite score within the five-enzyme grid:

- BSA: trypsin at 0.885 (`Good`).
- H3.1: Lys-C at 0.769 (`Good`).
- BACE1: chymotrypsin-high at 0.769 (`Good`).
- Lysozyme C: trypsin at 0.825 (`Good`).
- Ubiquitin fusion: Lys-C at 0.885 (`Good`).

Two pairs receive `Poor`: BACE1 with Lys-C at 0.389 and Lysozyme C with Glu-C at 0.269. These labels follow the current pepVet settings and are not experimental outcomes.

### Workflow preset effects

Six presets record a length range, GRAVY range, weight vector, and pI flag. Each artifact row applies the complete preset. Because some presets weight `S_unique`, all rows use the same supplied background made from the BSA, H3.1, and BACE1 tryptic digests. Uniqueness in this comparison means absence from those three supplied digests, not uniqueness in a species or search database. The figure reports length-valid counts only. GRAVY affects `S_hydro`, not the definition of a length-valid peptide.

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

Within these rows, the targeted preset has the narrowest length range and the FFPE preset has the broadest. Composite differences reflect the complete configurations, not only the displayed length-valid counts. Verdicts are heuristic score bands under those configurations, not validated outcomes.

### Pipeline with PeptideRanger

pepVet ranks enzyme-digest combinations under an explicit scoring model. PeptideRanger supplies a peptide-level random-forest score. The comparison below groups peptides by a combined length-and-GRAVY window. That grouping is narrower than pepVet's definition of a valid peptide, which depends on the active length range.

The analysis calculates the difference between mean PeptideRanger scores inside and outside the selected window across 25 protein-enzyme combinations. A positive value means that the inside-window group has the higher mean score in that combination. It does not validate either model.

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

The score difference is positive in `r sum(pr_data$mean_score_difference > 0, na.rm = TRUE)` of `r sum(!is.na(pr_data$mean_score_difference))` computable combinations. A missing value means that one comparison group is empty. Small window-pass groups can make the difference unstable, so this descriptive contrast is not evidence of improved experimental detection.

A two-stage pipeline run looks like this:

This code passes pepVet's length-valid peptides to PeptideRanger. It does not reproduce the narrower combined length-and-GRAVY grouping used in the descriptive artifact above.

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

### Binary versus graded: pepVet and Protein Cleaver

The comparison applies a simulated Protein Cleaver length-and-mass rule to the same peptides used by pepVet. Protein Cleaver uses a binary flag, while pepVet calculates digest-level components and a composite score.

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

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

Under the simulated rule, every peptide inside the selected pepVet window is also PC-identifiable. The Protein Cleaver rule accepts additional peptides because it uses a different length range and a mass window rather than a GRAVY window. The exact result depends on the simulated settings and does not replace a run of the external application.

## Choosing an output type

pepVet can be used for enzyme ranking and inspection of digest-level score components. PeptideRanger can add a peptide-level suitability score, while Protein Cleaver can add structural context when a PDB file is available. MS-Digest and ExPASy provide established peptide-mass calculations. The appropriate tool depends on the question and on whether the user needs a digest ranking, peptide prioritization, structural annotation, or mass calculation.

## References

Koulouras G, Xu Y. Protein cleaver: an interactive web interface for in silico prediction and systematic annotation of protein digestion-derived peptides. *Frontiers in Bioinformatics*. 2025, 5:1576317. doi: [10.3389/fbinf.2025.1576317](https://doi.org/10.3389/fbinf.2025.1576317).

Miller RM, Ibrahim K, Smith LM. ProteaseGuru: A Tool for Protease Selection in Bottom-Up Proteomics. *Journal of Proteome Research*. 2021, 20(4):1936-1942. doi: [10.1021/acs.jproteome.0c00954](https://doi.org/10.1021/acs.jproteome.0c00954).

Riley RM, Spencer Miko SE, Morin RD, Morin GB, Negri GL. PeptideRanger: An R Package to Optimize Synthetic Peptide Selection for Mass Spectrometry Applications. *Journal of Proteome Research*. 2023, 22(2):526-531. doi: [10.1021/acs.jproteome.2c00538](https://doi.org/10.1021/acs.jproteome.2c00538).

Wilkins MR, Lindskog I, Gasteiger E, Bairoch A, Sanchez JC, Hochstrasser DF, Appel RD. Detailed peptide characterization using PEPTIDEMASS, a World-Wide-Web-accessible tool. *Electrophoresis*. 1997, 18(3-4):403-408. doi: [10.1002/elps.1150180314](https://doi.org/10.1002/elps.1150180314).

## Session info

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