Evaluating pepVet Parameters Against PeptideAtlas
Source:vignettes/articles/peptideatlas-concordance.Rmd
peptideatlas-concordance.RmdpepVet’s default scoring parameters (length range
[7, 25], GRAVY range [-1.0, 0.6], verdict
thresholds 0.65 / 0.40) are
literature-informed priors and package design choices. The package uses
intentionally strict and conservative length and GRAVY windows. Before
this work, no analysis had evaluated their relationship with peptide
observations using pepVet outputs.
This article uses a January 2026 human PeptideAtlas build with 5.3 million unique observed peptide sequences.
Approach
We sampled 500 proteins from the UniProt human reference proteome, limited to sequences at least 50 residues long with standard amino acids. We digested each sequence with trypsin and one missed cleavage, then checked whether each theoretical peptide appeared in the PeptideAtlas sequence list.
The analysis has three parts:
Observation-rate difference: for each protein, compute the observation rate inside and outside the selected length-and-GRAVY window. Report the mean difference across proteins with a bootstrap interval.
Grid search: vary length and GRAVY boundaries across 900 specified combinations and identify the highest enrichment within that grid.
Threshold assessment: fit logistic regression for protein-level observation status and compute Youden’s J across a threshold grid. Inspect the class balance before interpreting the selected threshold.
In this article, “inside-window” means that a theoretical peptide
meets both the selected length and GRAVY criteria. It is not the
package’s length-valid definition, which depends only on the active
length range. The retained artifact columns use the earlier
valid and invalid names for these two
combined-window groups.
Results
The three analyses report the default-window association, the bounded parameter grid, and the protein-level threshold assessment separately.
Enrichment
Using the default parameters, the mean observation rates across proteins are:
valid_rate <- mean(pp$detection_rate_valid, na.rm = TRUE)
invalid_rate <- mean(pp$detection_rate_invalid, na.rm = TRUE)
enrichment <- valid_rate - invalid_rate| Metric | Value |
|---|---|
| Observation rate (inside selected window) | 57.3% |
| Observation rate (outside selected window) | 19.3% |
| Mean observation-rate difference | 38.0 percentage points |
| Observation-rate ratio | 3.0x |
The bootstrap 95% interval for the mean difference is [36.1, 40.0] percentage points. The interval is above zero for this sample and resampling procedure.
Within these sampled human tryptic digests, peptides inside the default windows had a higher mean PeptideAtlas observation rate than peptides outside them. This is concordance evidence for the sampled data. It does not estimate a causal effect of the filters or establish performance for other enzymes, organisms, instruments, or acquisition methods.
Signal decomposition
The full script calculates the observation-rate difference for each filter alone.
This analysis is done at the peptide level in the full script. In
this sample, the length filter [7, 25] alone has a 47.5
percentage-point difference, while the GRAVY filter
[-1.0, 0.6] alone has an 18.9 percentage-point difference.
The combined difference is 38.0 percentage points because the selected
sets overlap.
In this decomposition, the length filter had the larger individual association with observation status. The analysis does not establish a general ordering of determinants of LC-MS/MS detection.
GRAVY min finding
23.8% of observed peptides have GRAVY below -1.0. These are more hydrophilic peptides that the default window excludes. Lowering the GRAVY minimum to -1.5 raises enrichment from 38.0 to 42.1 percentage points while changing the inside-window observation rate from 57.3% to 57.0%. This was the largest individual change among the tested GRAVY-floor values.
The existing DIA preset includes wider GRAVY limits of [-1.0, 0.8]. The package default remains unchanged because this analysis covers one human tryptic sample.
Grid search
The full grid search evaluated 900 combinations of length and GRAVY boundaries.
top <- head(gs, 5)[, c("length_min", "length_max", "gravy_min", "gravy_max",
"detection_rate_valid", "detection_rate_invalid", "enrichment")]
knitr::kable(top, digits = 4)| length_min | length_max | gravy_min | gravy_max | detection_rate_valid | detection_rate_invalid | enrichment |
|---|---|---|---|---|---|---|
| 7 | 35 | -1.5 | 1.0 | 0.5416 | 0.0906 | 0.4511 |
| 7 | 30 | -1.5 | 1.0 | 0.5562 | 0.1056 | 0.4506 |
| 7 | 40 | -1.5 | 1.0 | 0.5303 | 0.0835 | 0.4467 |
| 8 | 30 | -1.5 | 1.0 | 0.5699 | 0.1245 | 0.4454 |
| 7 | 30 | -1.5 | 0.8 | 0.5583 | 0.1150 | 0.4434 |
The highest enrichment among the 900 evaluated combinations occurs at
length [7, 35] and GRAVY [-1.5, 1.0]: 45.1%,
or 7.1 percentage points above the default result (bootstrap 95%
interval for the difference [6.2, 7.9] percentage points). This is a
bounded grid-search result, not an estimate of universally optimal
boundaries.
The highest-enrichment boundaries in this grid are wider than the defaults in two ways:
- GRAVY: wider on both ends (-1.5 to 1.0 vs -1.0 to 0.6). The lower end contributes more to the improvement.
- Length: the upper boundary is 35 rather than 25. Within the tested grid, widening this boundary alone lowers enrichment, while the combination with the wider GRAVY range has the highest observed enrichment.
The wider selected range reflects a different tradeoff from the package defaults. pepVet retains strict, conservative defaults for triage rather than claiming that they are estimated detection boundaries. Users can inspect the DIA preset or custom ranges when broader inclusion is appropriate for their workflow.
Threshold assessment
The following logistic regression relates composite score to protein-level observation status in this sample:
# Re-fit from pre-computed data to display coefficients
logit_fit <- glm(I(n_observed > 0) ~ composite_score,
data = pp, family = binomial)
knitr::kable(summary(logit_fit)$coefficients, digits = 4)| Estimate | Std. Error | z value | Pr(>|z|) | |
|---|---|---|---|---|
| (Intercept) | 0.7001 | 0.7795 | 0.8981 | 0.3691 |
| composite_score | 2.6432 | 1.1706 | 2.2580 | 0.0239 |
In this sample, the fitted model associates composite score with binary protein-level observation status (coefficient test p = 0.024). A 0.1 increase in composite score corresponds to an estimated 0.26 increase in log-odds under this model. The analysis describes the sampled data and does not provide calibrated probabilities for new data.
Youden’s J across the threshold grid:
best_idx <- which.max(tc$youden_j)
best_th <- tc$threshold[best_idx]
knitr::kable(data.frame(
Threshold = best_th,
Sensitivity = tc$sensitivity[best_idx],
Specificity = tc$specificity[best_idx],
Youden_J = tc$youden_j[best_idx],
Default = 0.65
), digits = 3)| Threshold | Sensitivity | Specificity | Youden_J | Default |
|---|---|---|---|---|
| 0.875 | 0.944 | 0.231 | 0.174 | 0.65 |
The threshold with the highest Youden J on this grid is 0.875, compared with the default 0.65. However, 461 of 500 proteins (92.2%) have at least one observed peptide. With only 39 negative examples, the specificity and selected threshold are unstable. The maximum Youden J is 0.174.
This analysis does not provide a reliable calibration of the verdict thresholds. The package retains 0.65 and 0.40 as heuristic defaults.
Interpretation
The result supports a limited statement about the sampled human tryptic data. It does not calibrate pepVet for other datasets or workflows.
What this analysis supports
| Parameter | Status | Evidence |
|---|---|---|
| Length range [7, 25] | Concordant in this sample | 47.5 percentage-point difference alone. The highest-enrichment grid range was wider. |
| GRAVY range [-1.0, 0.6] | Concordant in this sample | 18.9 percentage-point difference alone. A lower floor increased the difference in the tested grid. |
| Verdict thresholds 0.65 / 0.40 | Retained, not validated | Class imbalance prevents reliable calibration |
| Expert-prior weight vector | Not tested | A separate sensitivity analysis measures stability only under its specified perturbation model |
Limitations
Human + trypsin only. The PeptideAtlas build is human-specific and trypsin is the only enzyme tested. Transfer to other organisms or enzymes was not evaluated.
No abundance weighting. Proteins that are not expressed in the samples that contributed to PeptideAtlas will appear as “unobserved” regardless of digest quality. PaxDB abundance scores were not integrated.
Binary observation status. PeptideAtlas records presence or absence, not detection confidence or intensity. The analysis treats a peptide observed once and a peptide observed in hundreds of experiments identically.
Protein-level threshold analysis has 39 negative examples. With 92.2% of proteins having at least one observed peptide, specificity estimates and the selected threshold are unstable.
Reproducibility
The full analysis script is at
data-raw/peptideatlas-concordance.R. It requires the
PeptideAtlas build FASTA (168 MB) and UniProt human reference proteome
(14 MB), which are not shipped with the package. Pre-computed results
are in inst/extdata/peptideatlas-concordance/.
Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] pepVet_0.99.0
#>
#> loaded via a namespace (and not attached):
#> [1] vctrs_0.7.3 crayon_1.5.3 cli_3.6.6
#> [4] knitr_1.51 rlang_1.3.0 xfun_0.60
#> [7] otel_0.2.0 generics_0.1.4 textshaping_1.0.5
#> [10] jsonlite_2.0.0 glue_1.8.1 S4Vectors_0.51.5
#> [13] Biostrings_2.81.5 htmltools_0.5.9 stats4_4.6.1
#> [16] ragg_1.5.2 sass_0.4.10 cleaver_1.51.0
#> [19] rmarkdown_2.31 Seqinfo_1.3.0 tibble_3.3.1
#> [22] evaluate_1.0.5 jquerylib_0.1.4 fastmap_1.2.0
#> [25] IRanges_2.47.2 yaml_2.3.12 lifecycle_1.0.5
#> [28] compiler_4.6.1 fs_2.1.0 pkgconfig_2.0.3
#> [31] XVector_0.53.0 systemfonts_1.3.2 digest_0.6.39
#> [34] R6_2.6.1 pillar_1.11.1 magrittr_2.0.5
#> [37] bslib_0.11.0 tools_4.6.1 pkgdown_2.2.1
#> [40] BiocGenerics_0.59.10 cachem_1.1.0 desc_1.4.3