Weight Sensitivity Analysis • pepVet

Analysis scope

pepVet’s scoring weights (S_length = 0.200, S_coverage = 0.348, S_count = 0.226, S_hydro = 0.138, S_charge = 0.088) determine how each component contributes to the composite. sensitivity_analysis() perturbs these weights and records how often the verdict or enzyme ranking changes. It asks whether the result is stable within the specified perturbation model.

In this article, “stable” means that the verdict is unchanged in more than 95% of perturbations. This operational label does not show that the weights or verdict are experimentally correct. A less stable result indicates that the score lies close enough to a decision boundary for modest changes in weighting to alter the label.

Single-protein examples

The BSA and Histone H3.1 examples contrast a stable result with one that changes near a verdict boundary under the same perturbation model.

A stable case: BSA with trypsin

bsa_path <- system.file("extdata", "P02769.fasta", package = "pepVet")
res <- evaluate_digest(bsa_path, enzyme = "trypsin", missed_cleavages = 1L)
withr::local_seed(20260714)
sens <- sensitivity_analysis(res, n_iter = 5000L)

The default composite for BSA trypsin is 0.885, above the 0.65 Good threshold. The simulation records the verdict in each weight draw:

sens$summary$verdict_pct
#>     Good Moderate     Poor 
#>        1        0        0

All 5,000 draws return Good. The central 95% simulation interval for the composite is 0.851 to 0.916. This result establishes stability only for the configured Dirichlet distribution and random seed.

A borderline case: Histone H3.1 with trypsin

h3_path <- system.file("extdata", "P68431.fasta", package = "pepVet")
res_h3 <- evaluate_digest(h3_path, enzyme = "trypsin", missed_cleavages = 1L)
withr::local_seed(20260714)
sens_h3 <- sensitivity_analysis(res_h3, n_iter = 5000L)

Histone H3.1 is lysine- and arginine-rich. Tryptic digestion produces many short products, and the default composite sits near the 0.65 boundary. The sensitivity analysis measures how often the label changes under the specified weight perturbations:

sens_h3$summary$verdict_pct
#>     Good Moderate     Poor 
#>   0.0716   0.9284   0.0000

The verdict changes between Good and Moderate across the sampled weights. In this deterministic render, 7.2% of draws are Good and 92.8% are Moderate. The result is sensitive to weighting near the 0.65 boundary. Compare other enzymes if that sensitivity affects the decision, but evaluate them under the same scoring configuration.

Visualising single-protein verdict distributions

p <- plot_weight_sensitivity(sens_h3)
p

The density spans both sides of the 0.65 boundary. In the BSA simulation, all sampled composites remain above that boundary:

p_bsa <- plot_weight_sensitivity(sens)
p_bsa

Batch analysis: 50 proteins across 10 enzymes

The batch example applies the same perturbation to the package’s 50-protein documentation fixture and 10 enzymes. It demonstrates batch behavior without treating this small fixture as a proteome-wide validation set.

small_path <- system.file("extdata", "small_proteome_50_proteins.fasta",
  package = "pepVet")
enzymes <- c("trypsin", "lysc", "chymotrypsin-high", "asp-n endopeptidase",
  "glutamyl endopeptidase", "arg-c proteinase", "thermolysin",
  "pepsin", "staphylococcal peptidase i", "proteinase k")
batch_comp <- batch_compare_enzymes(small_path, enzymes = enzymes)
#> Scoring 50 proteins against 10 enzymes.
#> Warning: Protein "sp|P08246|ELNE_HUMAN Neutrophil elastase OS=Homo sapiens OX=9606
#> GN=ELANE PE=1 SV=1" has no cleavage sites for "lysc". S_count and composite
#> score set to 0; verdict set to Poor.
#> Warning: Protein "sp|P31358|CD52_HUMAN CAMPATH-1 antigen OS=Homo sapiens OX=9606 GN=CD52
#> PE=1 SV=1" has no cleavage sites for "glutamyl endopeptidase". S_count and
#> composite score set to 0; verdict set to Poor.
#> Warning: Protein "sp|P31358|CD52_HUMAN CAMPATH-1 antigen OS=Homo sapiens OX=9606 GN=CD52
#> PE=1 SV=1" has no cleavage sites for "staphylococcal peptidase i". S_count and
#> composite score set to 0; verdict set to Poor.
withr::local_seed(20260714)
sens <- sensitivity_analysis(batch_comp, n_iter = 2000L)

Overall instability across the supplied fixture:

sens$summary$total_instability
#> [1] 0.081841

This means 8.2% of all protein-enzyme-iteration triples produced a different verdict than the default weights would give.

Pairs with fewer than 5% verdict changes

mean(sens$per_protein$verdict_instability < 0.05)
#> [1] 0.668

67% of protein-enzyme pairs have a verdict that changes in fewer than 5% of iterations. These verdicts are stable under the tested perturbations, which is narrower than evidence of experimental validity.

Scores near verdict boundaries

The rows with the largest instability have mean composites near the 0.65 or 0.40 thresholds in this simulation.

worst <- head(sens$per_protein[order(sens$per_protein$verdict_instability,
  decreasing = TRUE), ], 10)
worst[, c("protein_id", "verdict_instability", "composite_mean",
  "composite_lo", "composite_hi")]
#> # A tibble: 10 × 5
#>    protein_id       verdict_instability composite_mean composite_lo composite_hi
#>    <chr>                          <dbl>          <dbl>        <dbl>        <dbl>
#>  1 sp|Q92956|TNR14…               0.498          0.400        0.350        0.449
#>  2 sp|P13667|PDIA4…               0.485          0.649        0.617        0.682
#>  3 sp|P08246|ELNE_…               0.481          0.401        0.353        0.449
#>  4 sp|P15529|MCP_H…               0.476          0.651        0.598        0.701
#>  5 sp|P06400|RB_HU…               0.475          0.652        0.600        0.706
#>  6 sp|P01889|HLAB_…               0.468          0.652        0.600        0.703
#>  7 sp|O75475|PSIP1…               0.468          0.403        0.348        0.463
#>  8 sp|Q6ZWK4|RHEX_…               0.462          0.648        0.603        0.698
#>  9 sp|O95273|CCDB1…               0.462          0.652        0.608        0.692
#> 10 sp|P08174|DAF_H…               0.461          0.653        0.598        0.708

The displayed rows have simulation intervals that cross a verdict boundary, matching the pattern in the H3.1 example.

Instability by enzyme

Mean verdict instability is the fraction of iterations in which the verdict differs from the default-weight verdict. A value of 0.05 means 5% of iterations changed the label.

instability_by_enzyme <- sort(tapply(sens$per_protein$verdict_instability,
  sens$per_protein$enzyme, mean))
instability_by_enzyme
#>                    trypsin          chymotrypsin-high 
#>                    0.04016                    0.04491 
#>               proteinase k                       lysc 
#>                    0.05671                    0.08269 
#>                thermolysin           arg-c proteinase 
#>                    0.08275                    0.09271 
#>                     pepsin        asp-n endopeptidase 
#>                    0.09371                    0.09635 
#>     glutamyl endopeptidase staphylococcal peptidase i 
#>                    0.11250                    0.11592

For trypsin, a mean of 4% means that this fraction of sampled iterations changed the verdict, averaged across the 50 proteins. The displayed means describe this fixture and perturbation model. They do not isolate why enzymes differ.

Visualising the batch results

p_batch <- plot_weight_sensitivity(sens)
p_batch

Each panel shows per-protein verdict instability for one enzyme. Values near zero changed label rarely under the sampled weights. Larger values identify protein-enzyme pairs that depend more on the weight vector.

Interpretation limits

The analysis identifies verdicts and rankings that depend on the tested weights. It does not test the weights or verdict thresholds against observed peptide detection. Neither high nor low stability establishes experimental success or failure.

Technical notes

The perturbation is a package-chosen Dirichlet distribution with concentration parameters alpha_i = nu * w0_i, centered on the resolved reference weights. At the default nu = 63, a reference weight of 0.2 has standard deviation 0.05 and an exact marginal 95% interval of approximately 0.112 to 0.306. Dirichlet draws are non-negative and sum to 1, so they preserve the weight simplex. A reference weight of zero remains zero in every draw. The cited simulation-based multicriteria literature motivates exploring uncertain weight spaces, but it does not empirically validate this distribution or concentration for proteomics.

For batch inputs, the analysis processes 100 score rows per chunk by default. It classifies simulated scores without constructing a character verdict matrix and aggregates optional squared weight-score associations within each chunk.

Proteins with no cleavage sites for a given enzyme receive composite = 0 and verdict = “Poor” as a hard override. This ensures undigestible proteins remain negative controls under all weight perturbations.

References

Butler J, Jia J, Dyer J. Simulation techniques for the sensitivity analysis of multi-criteria decision models. European Journal of Operational Research. 1997. 103:531-546.

Lahdelma R, Hokkanen J, Salminen P. SMAA: Stochastic multiobjective acceptability analysis. European Journal of Operational Research. 1998. 106(1):137-143.

Mazurek J, Strzalka D. On the Monte Carlo weights in multiple criteria decision analysis. PLOS ONE. 2022. 17(10):e0268950.

Tervonen T, Figueira JR. A survey on stochastic multicriteria acceptability analysis methods. Journal of Multi-Criteria Decision Analysis. 2008. 15(1-2):1-14.

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] gtable_0.3.6         jsonlite_2.0.0       compiler_4.6.1      
#>  [4] crayon_1.5.3         Biostrings_2.81.5    jquerylib_0.1.4     
#>  [7] scales_1.4.0         systemfonts_1.3.2    IRanges_2.47.2      
#> [10] Seqinfo_1.3.0        textshaping_1.0.5    yaml_2.3.12         
#> [13] fastmap_1.2.0        ggplot2_4.0.3        R6_2.6.1            
#> [16] XVector_0.53.0       cleaver_1.51.0       labeling_0.4.3      
#> [19] generics_0.1.4       knitr_1.51           BiocGenerics_0.59.10
#> [22] tibble_3.3.1         desc_1.4.3           RColorBrewer_1.1-3  
#> [25] bslib_0.11.0         pillar_1.11.1        rlang_1.3.0         
#> [28] utf8_1.2.6           cachem_1.1.0         xfun_0.60           
#> [31] S7_0.2.2             fs_2.1.0             sass_0.4.10         
#> [34] otel_0.2.0           cli_3.6.6            pkgdown_2.2.1       
#> [37] withr_3.0.3          magrittr_2.0.5       grid_4.6.1          
#> [40] digest_0.6.39        lifecycle_1.0.5      S4Vectors_0.51.5    
#> [43] vctrs_0.7.3          evaluate_1.0.5       glue_1.8.1          
#> [46] farver_2.1.2         ragg_1.5.2           stats4_4.6.1        
#> [49] rmarkdown_2.31       tools_4.6.1          pkgconfig_2.0.3     
#> [52] htmltools_0.5.9