NADIA 0.99.3
In most proteomics experiments, the proteins that truly change between
conditions are unknown. The comparisons in vignette("choosing-methods")
therefore rely on indirect criteria. A spike-in experiment provides a known
answer instead: proteins from one or more foreign organisms are mixed into a
constant background at predefined ratios. We know in advance which species
should change, the expected direction and size of each change, and that the
background proteins should remain unchanged.
benchmarking_proteomics() compares a differential-expression result with this
known experimental design. Given the DE table, a protein-to-species annotation,
and the expected fold changes, it classifies proteins as TP, FP, TN, or FN;
calculates performance metrics; and produces diagnostic tables and plots. This
vignette explains those classification rules and shows how to interpret the
resulting benchmark.
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("NADIA")
This vignette also uses highcharter, which is optional;
BiocManager::install("NADIA", dependencies = TRUE) installs it together
with NADIA.
library(NADIA)
data(nadia_dia)
The example dataset originates from a three-proteome spike-in experiment. Each sample contains a constant 80% Homo sapiens (HeLa) background, while Escherichia coli and Saccharomyces cerevisiae proteins make up the remaining 20% in reciprocal proportions. Across the complete experimental design, the E. coli fraction increases and the yeast fraction decreases from condition A to condition D. The four condition mixtures were prepared from the same HeLa stock and then divided into four aliquots, each of which was carried independently through the complete sample-preparation workflow and LC–ESI–MS/MS analysis. These are therefore termed pseudo-biological replicates: they are not independent biological specimens, but they capture more experimental variation than repeated injections of the same prepared sample.
Figure 1: Complete three-proteome spike-in design underlying the benchmark dataset
HeLa remains constant at 80%, while E. coli increases and yeast decreases across conditions. The dataset distributed with NADIA retains conditions A, B, and D and 2,000 protein groups.
The reduced dataset distributed with NADIA retains conditions A, B, and D: 12 samples and 2,000 protein groups. These groups were selected at random with a fixed seed rather than by abundance or data completeness, preserving a realistic missing-value structure. Condition D contains nominally 0% yeast, which produces the strongest biological missingness in the experiment and makes this condition particularly informative for assessing whether a pipeline recovers the expected biological changes.
The processed nadia_dia object identifies proteins through
PG.ProteinGroups, but it deliberately keeps the experimental species
annotation in a separate file. The following expandable example shows how that
file is read, how its organism column is converted into the species_df input
required by benchmarking_proteomics(), and a few representative rows from
each species.
spikein <- read.delim(system.file("extdata", "nadia_dia_spikein.tsv.gz",
package = "NADIA"))
species_df <- data.frame(Protein.IDs = spikein$PG.ProteinGroups,
Species = spikein$PG.OrganismId)
species_example <- do.call(rbind, lapply(c("HUMAN", "ECOLI", "YEAST"),
function(species) head(species_df[species_df$Species == species, ], 2)))
rownames(species_example) <- NULL
species_example
#> Protein.IDs Species
#> 1 A0A024RBG1 HUMAN
#> 2 A0A140T897;P02769 HUMAN
#> 3 P00370 ECOLI
#> 4 P00550 ECOLI
#> 5 O13535 YEAST
#> 6 O13563 YEAST
table(species_df$Species)
#>
#> ECOLI HUMAN YEAST
#> 338 1256 406
Each protein in species_df must map to exactly one species. NADIA collapses
exact duplicate rows with a message, rejects proteins assigned to conflicting
species, and verifies that joining the annotation does not change the number of
DE rows. These checks prevent duplicated mappings from silently inflating the
confusion matrix and every metric derived from it.
expected_values says what each species should do in each comparison, as a
log2 fold change. This is the specification of the experiment, taken from how
much of each organism was mixed in — not something derived from the data.
NADIA accepts either a base R data.frame or a tibble. The first option below
is executed in this vignette; the commented alternative creates the same
expected object with tibble::tribble().
# Option 1: base R data.frame
expected <- data.frame(
Comparison = rep(c("B-A", "D-A", "D-B"), each = 2),
Species = rep(c("ECOLI", "YEAST"), times = 3),
expected_logFC = c(1, -0.58, 2, -3.3, 1, -2.72))
# Option 2: equivalent row-wise tibble
# expected <- tibble::tribble(
# ~Comparison, ~Species, ~expected_logFC,
# "B-A", "ECOLI", 1,
# "B-A", "YEAST", -0.58,
# "D-A", "ECOLI", 2,
# "D-A", "YEAST", -3.3,
# "D-B", "ECOLI", 1,
# "D-B", "YEAST", -2.72
# )
expected
The resulting table contains one expected change for each spike-in species and comparison:
#> Comparison Species expected_logFC
#> 1 B-A ECOLI 1.00
#> 2 B-A YEAST -0.58
#> 3 D-A ECOLI 2.00
#> 4 D-A YEAST -3.30
#> 5 D-B ECOLI 1.00
#> 6 D-B YEAST -2.72
Condition D requires a practical convention. Its nominal yeast proportion is
0%, so the exact yeast fold change relative to a non-zero condition would be
zero on the ratio scale and -Inf after log2 transformation. Because the
benchmark requires a finite expected value, NADIA represents this complete
depletion with a finite floor. For D-A, the chosen value is approximately
-3.3, corresponding to a tenfold decrease (2^-3.3 is approximately 0.10).
Using the same floor relative to condition B gives approximately -2.72 for
D-B. These values are an operational scoring convention; they do not imply
that yeast was intentionally added to condition D.
Because this table defines the biological truth, NADIA validates it before
scoring. expected_logFC must be numeric, finite, and non-zero.
Keep species annotation and expected changes separate. Every protein that
may enter the benchmark, including the human background, must have its species
recorded in species_df. By contrast, expected_values contains only the
species expected to change in each comparison. Human proteins are therefore
present in species_df but absent from expected_values, so NADIA treats them
as background: a significant human protein is an FP and a non-significant one
is a TN. The same rule applies to any other species omitted from
expected_values; background species should not be added with
expected_logFC = 0.
The species annotation and the expected changes now provide the experimental
truth against which the pipeline can be scored. First,
process_proteomics() runs the complete analysis and produces a
differential-expression table. benchmarking_proteomics() then compares every
protein in that table with the declared truth, assigns TP, FP, TN, or FN, and
summarizes how well the expected changes were recovered.
benchmarking_proteomics() returns a named list containing numerical tables and
their corresponding diagnostic plots. Its principal inputs and controls are
summarized below.
Data and experimental truth
de_res: the differential-expression table, including protein IDs, log2 fold changes, p-values, and comparisons.
expected_values: the expected Comparison, Species, and expected_logFC combinations.
species_df: the protein-to-species mapping. It is needed unless de_res already contains a Species column.
Classification settings
alpha: significance threshold; the default is 0.05.
lfc_thr: minimum absolute log2 fold change; the default is 0.
p_col: p-value column used for classification; the default is adj.P.Val.
Selection and output
comparisons and assay: optional filters selecting the results to score.
output_dir: export directory; NULL keeps the results in memory without writing files.
verbose and species_colors: control progress messages and plot colors.
The diagram groups the outputs by the question they answer. The following sections inspect these results individually.
Figure 2: Inputs, classification steps, and output families produced by benchmarking_proteomics()
In this example, res$DEPs_results is the differential-expression table,
expected is the truth table defined above, and species_df supplies the
protein-to-species annotation. Setting output_dir = NULL keeps all outputs in
the returned bench object.
res <- process_proteomics(nadia_dia, verbose = FALSE)
bench <- benchmarking_proteomics(res$DEPs_results, expected,
species_df = species_df,
output_dir = NULL, verbose = FALSE)
The first summary to inspect is metrics_table. Each row represents one
comparison and the following code selects the principal classification and
ranking indicators.
metric_cols <- c("Comparison", "Sensitivity", "Specificity", "Precision",
"F1", "MCC", "AUC")
knitr::kable(bench$metrics_table[, metric_cols], digits = 3,
caption = "Classification and ranking metrics by comparison.")
| Comparison | Sensitivity | Specificity | Precision | F1 | MCC | AUC |
|---|---|---|---|---|---|---|
| B-A | 0.631 | 0.979 | 0.946 | 0.757 | 0.682 | 0.832 |
| D-A | 0.933 | 0.952 | 0.920 | 0.926 | 0.882 | 0.978 |
| D-B | 0.847 | 0.972 | 0.947 | 0.894 | 0.840 | 0.949 |
All six indicators are interpreted so that better performance gives a larger value, except that MCC can also be negative:
The reported metrics follow from how each protein is matched to the experimental truth and to the result returned by the differential-expression analysis.
Two independent things go into every classification.
truth — what the protein is. It is 1 if the protein belongs to a
spike-in species with a declared expectation, 0 if it is background. This comes
from species_df and expected_values, that is, from how the experiment was
mixed. It is a property of the protein and no analysis can change it.
predicted — what the pipeline said. It is 1 only if the protein is
significant (adj.P.Val <= alpha, with alpha = 0.05 by default) and its
logFC has the sign the experiment specified. Significance alone is not enough.
Crossing them gives the four cells:
truth = 1 (spike-in) |
truth = 0 (background) |
|
|---|---|---|
predicted = 1 |
TP — detected, right direction | FP — background wrongly flagged |
predicted = 0 |
FN — missed, or found with the wrong sign | TN — correctly left alone |
The same rules in biological terms are:
| Protein scenario | Analysis result | Class |
|---|---|---|
| Significant human protein | A change that should not exist is called significant | FP |
| Non-significant human protein | The background is correctly left unchanged | TN |
| Significant spike-in, correct sign | The expected change is recovered | TP |
| Non-significant spike-in | A real expected change is missed | FN |
| Significant spike-in, wrong sign | The expected direction is not recovered | FN |
The FN label alone does not distinguish a spike-in protein that was not
significant from one that was significant in the wrong direction. To preserve
this useful diagnostic information, classified_df includes two additional
columns without changing the classification metrics. is_significant records
whether the protein passed both the adjusted p-value and fold-change thresholds,
irrespective of direction. direction_error identifies the subset of FNs that
passed those thresholds but had the opposite sign. The following code retrieves
these cases; there are nine in this analysis:
wrong_sign <- subset(bench$classified_df, direction_error)
wrong_sign <- merge(wrong_sign, expected,
by = c("Comparison", "Species"))
nrow(wrong_sign)
table(wrong_sign$Comparison, wrong_sign$Species)
wrong_sign[, c("Protein.IDs", "Comparison", "Species", "logFC",
"expected_logFC", "adj.P.Val", "is_significant",
"direction_error", "classification")]
#> [1] 9
#>
#> YEAST
#> B-A 1
#> D-A 4
#> D-B 4
#> Protein.IDs Comparison Species logFC expected_logFC adj.P.Val
#> 1 Q04231 B-A YEAST 6.848494 -0.58 3.623551e-05
#> 2 P25574 D-A YEAST 1.430882 -3.30 1.897355e-03
#> 3 P33202 D-A YEAST 1.606620 -3.30 1.010409e-03
#> 4 Q04231 D-A YEAST 4.645648 -3.30 4.141794e-04
#> 5 Q06408 D-A YEAST 1.340623 -3.30 3.315274e-05
#> 6 P25574 D-B YEAST 1.928146 -2.72 2.375836e-04
#> 7 P33202 D-B YEAST 1.629626 -2.72 1.097067e-03
#> 8 Q03835 D-B YEAST 0.642464 -2.72 3.524624e-02
#> 9 Q06408 D-B YEAST 1.566569 -2.72 1.037085e-05
#> is_significant direction_error classification
#> 1 TRUE TRUE FN
#> 2 TRUE TRUE FN
#> 3 TRUE TRUE FN
#> 4 TRUE TRUE FN
#> 5 TRUE TRUE FN
#> 6 TRUE TRUE FN
#> 7 TRUE TRUE FN
#> 8 TRUE TRUE FN
#> 9 TRUE TRUE FN
All nine are yeast proteins that rose significantly even though the spike-in
design specified a decrease. They occur across all three comparisons, whose
expected yeast changes are -0.58, -3.3, and -2.72. Every one is an FN and
retains truth = 1.
NADIA also provides gg_signif_bars, which displays significant proteins by
species, comparison, and observed direction. Because the plot uses
is_significant rather than predicted, these nine FNs remain visible in the
UP facet:
bench$gg_signif_bars
metrics_table summarizes each comparison with one value per performance
indicator. A confusion plot instead displays the four underlying classification
outcomes—TP, FP, TN, and FN—and therefore shows whether errors arise from
missed spike-in proteins or from significant calls in the unchanged background.
The general purpose of these plots is to turn the summary benchmark scores into a diagnostic view of pipeline performance. They show whether the analysis recovers the expected biological changes while keeping the background stable, allow comparisons to be assessed side by side, and identify the species or error type responsible for weaker performance. This makes it possible to interpret why a metric is high or low rather than treating it as an isolated number.
The species-level plot separates these outcomes by organism within each comparison. Each stacked bar represents all proteins from one species and shows the percentage assigned to each applicable class. This is the most informative view for a spike-in experiment because it reveals whether recovery differs between E. coli, yeast, and the human background.
bench$gg_confusion_by_species
The dependence on effect size is immediately visible. In B-A, 87.9% of the
E. coli proteins are recovered as TP, whereas only 42.5% of the yeast proteins
are recovered. The expected yeast change in this comparison is a log2 fold
change of -0.58, equivalent to a 1.5-fold difference and therefore relatively
difficult to detect. In D-A, where the expected yeast log2 fold change is
-3.3, recovery rises to 90.9%. The benchmark is thus evaluating the pipeline
against explicitly defined effect sizes, and smaller effects are expected to be
more difficult to recover.
The plot is convenient for comparing patterns, whereas the corresponding table
provides the exact percentage for every class. It also reports N, the number
of proteins evaluated for each species and comparison.
species_cols <- c("Comparison", "Species", "N", "TP_pct", "FP_pct",
"FN_pct", "TN_pct")
knitr::kable(bench$confusion_by_species[, species_cols], digits = 1,
caption = "Classification percentages within each species.")
| Comparison | Species | N | TP_pct | FP_pct | FN_pct | TN_pct |
|---|---|---|---|---|---|---|
| B-A | HUMAN | 1254 | 0.0 | 2.2 | 0.0 | 97.8 |
| B-A | YEAST | 405 | 42.5 | 0.0 | 57.5 | 0.0 |
| B-A | ECOLI | 338 | 87.9 | 0.0 | 12.1 | 0.0 |
| D-A | HUMAN | 1254 | 0.0 | 4.8 | 0.0 | 95.2 |
| D-A | YEAST | 405 | 90.9 | 0.0 | 9.1 | 0.0 |
| D-A | ECOLI | 338 | 96.2 | 0.0 | 3.8 | 0.0 |
| D-B | HUMAN | 1254 | 0.0 | 2.8 | 0.0 | 97.2 |
| D-B | YEAST | 405 | 86.9 | 0.0 | 13.1 | 0.0 |
| D-B | ECOLI | 338 | 82.0 | 0.0 | 18.0 | 0.0 |
The zero-valued columns follow directly from the role of each organism in the
experimental design. Human proteins form the unchanged background, so they can
only be TN when correctly left non-significant or FP when called significant;
their TP and FN percentages must therefore be zero. E. coli and yeast are the
expected-changing species, so their proteins can only be TP when the expected
change is recovered or FN when it is missed or detected in the wrong direction;
their FP and TN percentages must be zero. Within each spike-in species, TP_pct
is consequently its recovery rate and FN_pct is its missed or incorrectly
directed fraction.
For a more compact summary, the overall plot pools the three species within each comparison. This makes the total balance of correct and incorrect classifications easy to compare. However, it hides which spike-in species contributed the TP and FN proteins, and its percentages use all scored proteins as the denominator rather than the class-specific denominators used for sensitivity and specificity.
bench$gg_confusion_overall
Correct classifications dominate the overall composition: TP and TN together
represent 84.9% of the scored proteins in B-A, 94.5% in D-A, and 92.5% in
D-B. The larger FN segment in B-A (13.7%) is consistent with the weak yeast
effect in that comparison. FP proteins account for only 1.4%, 3.0%, and 1.8% of
all proteins in B-A, D-A, and D-B, respectively. These are percentages of
the complete benchmark and should be interpreted as the composition of each
bar; class-specific performance is reported separately in metrics_table.
The figure emphasizes relative proportions. The following table complements it with the exact TP, FP, TN, and FN counts and the resulting accuracy for each comparison.
overall_cols <- c("Comparison", "N", "TP", "FP", "FN", "TN", "Accuracy")
knitr::kable(bench$confusion_overall[, overall_cols], digits = 3,
caption = "Classification counts across all scored proteins.")
| Comparison | N | TP | FP | FN | TN | Accuracy |
|---|---|---|---|---|---|---|
| B-A | 1997 | 469 | 27 | 274 | 1227 | 0.849 |
| D-A | 1997 | 693 | 60 | 50 | 1194 | 0.945 |
| D-B | 1997 | 629 | 35 | 114 | 1219 | 0.925 |
All three comparisons evaluate the same 1,997 proteins. B-A contains 469 TP
and 274 FN, giving the lowest accuracy (0.849). Recovery improves markedly in
D-A, with 693 TP and only 50 FN, and this comparison achieves the highest
accuracy (0.945) despite having 60 FP. D-B is intermediate, with 629 TP, 114
FN, and an accuracy of 0.925. The table is included because the exact counts are
the quantities from which the summary metrics are calculated and cannot be read
precisely from the percentage-based plot alone.
The metrics heatmap provides a compact comparison of pipeline performance across all three contrasts. Columns represent comparisons and rows show six metrics: AUC, sensitivity, specificity, precision, F1, and accuracy. Each cell contains the numerical value of the metric, while the common red-to-green scale highlights lower and higher performance, respectively. Because all six metrics range from 0 to 1 and larger values indicate better performance, the heatmap makes both consistently strong results and comparison-specific weaknesses easy to identify.
bench$gg_heatmap
The heatmap identifies D-A as the strongest comparison overall. It has the
highest AUC (0.978), sensitivity (0.933), F1 (0.926), and accuracy (0.945),
showing that the expected changes are both well ranked and successfully
recovered. D-B also performs well, but its lower sensitivity (0.847) indicates
that it misses more spike-in proteins than D-A; its high specificity (0.972)
and precision (0.947) show that background control and the reliability of its
positive calls remain strong. B-A is the most difficult comparison. Although
its specificity (0.979) and precision (0.946) are high, its sensitivity falls to
0.631, which also lowers F1 to 0.757 and AUC to 0.832. Together with the
species-level results in Section 7.1, this pattern shows that the main weakness
in B-A is failure to recover the small expected yeast change, rather than an
excess of significant calls among human background proteins.
Spike-in benchmarks are often imbalanced: here, 1,254 of 1,997 proteins are unchanged human background. A method that called nothing significant would therefore achieve 62.8% accuracy, showing how accuracy alone can conceal poor recovery of the expected changes.
The OpDEA framework developed by Peng et al. (2024) evaluates proteomics differential-expression workflows using nMCC, G-mean, and pAUC at FPR limits of 1%, 5%, and 10%. nMCC and G-mean assess classification across spike-in and background proteins, whereas pAUC measures ranking in the low-FPR region. Together, they reduce the influence of class imbalance and reveal poor spike-in recovery that strong background performance might otherwise conceal.
NADIA reports these five values for every comparison in opdea_metrics:
knitr::kable(bench$opdea_metrics, digits = 3,
caption = "OpDEA metrics by comparison.")
| Comparison | nMCC | G_mean | pAUC_001 | pAUC_005 | pAUC_010 |
|---|---|---|---|---|---|
| B-A | 0.841 | 0.786 | 0.768 | 0.807 | 0.821 |
| D-A | 0.941 | 0.942 | 0.912 | 0.951 | 0.961 |
| D-B | 0.920 | 0.907 | 0.868 | 0.912 | 0.925 |
[0, 1]. It uses
TP, FP, TN, and FN in a balanced measure that is robust to class imbalance.
A value of 1 is perfect, 0.5 indicates no association, and values below 0.5
indicate disagreement with the expected classes.sqrt(Sensitivity * Specificity). It is high only when the pipeline both
recovers spike-in proteins and preserves the background, and becomes zero if
either component is zero.pAUC_001, pAUC_005, and pAUC_010. The McClish-corrected scale ranges from
0.5 (random) to 1 (perfect) and focuses the evaluation on stringent low-FPR
regions.All five indicators identify D-A as the strongest comparison. Its nMCC and
G-mean are both approximately 0.94, and it has the highest pAUC across all FPR
limits (0.912–0.961), indicating balanced classification and strong ranking
even in the most stringent low-FPR region.
D-B ranks second, with nMCC = 0.920, G-mean = 0.907, and pAUC values of
0.868–0.925. B-A performs worst: its lower G-mean (0.786) and pAUC
(0.768–0.821) reflect reduced sensitivity to the small yeast effect. The
consistent pattern across the three metric families indicates that the main
limitation in B-A is recovery of the expected spike-in changes, rather than
control of the human background.
A receiver operating characteristic (ROC) curve evaluates how well a continuous
score separates two known classes across every possible threshold. In this
benchmark, spike-in proteins are the positive class, background proteins are the
negative class, and the ranking score is derived from adj.P.Val. Each point on
the curve combines the true positive rate (TPR, or sensitivity) with the false
positive rate (FPR, or 1 - Specificity) obtained at one threshold.
Curves closer to the upper-left corner indicate better separation: they recover more spike-in proteins while misclassifying fewer background proteins. The area under the curve (AUC) summarizes the complete curve from 0 to 1. An AUC of 1 represents perfect ranking, 0.5 corresponds to random ranking, and values below 0.5 indicate that the classes are ranked in the opposite order.
The first plot shows the complete FPR and TPR ranges from 0 to 1.
bench$gg_roc
The full AUC is highest for D-A (0.978), followed by D-B (0.949) and B-A
(0.832). This agrees with the previous results: the larger expected changes are
ranked more effectively, whereas the small yeast effect makes B-A more
difficult. The legend preserves the input order (B-A, D-A, D-B); it is
not sorted by AUC, so performance should be read from the displayed values.
The full AUC averages performance over all possible FPRs, including very large rates that would rarely be accepted in a proteomics differential-expression analysis. Most analyses use stringent adjusted p-value thresholds and operate near the left edge of the ROC curve, rather than accepting large fractions of significant background proteins. Consequently, full AUC values can obscure differences in the practically relevant region, and a ranking based on the full curve can shrink, widen, or even reverse when only low FPRs are considered.
To examine this operationally relevant part of the curve, NADIA restricts the
x-axis to FPR <= 0.10. The zoom is not a different ROC analysis; it enlarges
the first 10% of the same curves so that differences near stringent operating
points can be seen.
bench$gg_roc_zoom
This figure differs from the full curve in two ways. First, the legend reports
the McClish-corrected pAUC over FPR <= 0.05, although the axis extends to 10%
to provide context around that interval. Second, the legend is deliberately
sorted by decreasing pAUC_005: D-A (0.951), D-B (0.912), and B-A
(0.807). These values focus on a more relevant operating range and confirm that
D-A provides the strongest ranking under a low false-positive-rate constraint.
Both AUC and pAUC use adj.P.Val to rank spike-in and background proteins; they
do not consider the expected fold-change direction. A significant spike-in with
the wrong sign can therefore rank highly in a ROC analysis while remaining an
FN in the direction-aware confusion matrix. ROC and confusion metrics should be
interpreted together.
An FPR range is not numerically equivalent to an adjusted p-value threshold.
adj.P.Val <= 0.05 defines which proteins the analysis calls significant,
whereas the resulting FPR is the proportion of background proteins that are
called positive. The following table links the two empirically by reporting the
observed operating point at the threshold used in this vignette. TPR is the
proportion of spike-in proteins called positive at the same threshold.
op <- do.call(rbind, lapply(unique(bench$classified_df$Comparison), function(comp) {
d <- subset(bench$classified_df, Comparison == comp)
called_positive <- !is.na(d$adj.P.Val) & d$adj.P.Val <= 0.05
data.frame(
Comparison = comp,
FPR_at_adjP005 = sum(called_positive & d$truth == 0) / sum(d$truth == 0),
TPR_at_adjP005 = sum(called_positive & d$truth == 1) / sum(d$truth == 1))
}))
knitr::kable(op, digits = 4,
col.names = c("Comparison", "FPR", "TPR"),
caption = "Observed ROC operating points at adj.P.Val <= 0.05.")
| Comparison | FPR | TPR |
|---|---|---|
| B-A | 0.0215 | 0.6326 |
| D-A | 0.0478 | 0.9381 |
| D-B | 0.0279 | 0.8520 |
At adj.P.Val <= 0.05, the observed FPR is 2.2% for B-A, 4.8% for D-A,
and 2.8% for D-B; the corresponding ROC TPRs are 63.3%, 93.8%, and 85.2%.
These TPRs are slightly higher than the direction-aware sensitivity values in
metrics_table because ROC analysis counts a significant spike-in as positive
even when its observed fold-change direction is incorrect. All
three points therefore lie within both the displayed 10% zoom and the 5% region
summarized by pAUC_005. This correspondence is specific to these results and
must not be interpreted as a general identity between adj.P.Val <= 0.05 and
FPR <= 0.05.
The table also explains why the zoom is more informative for this analysis: it
resolves the region containing the actual operating points, whereas most of the
full curve lies far beyond them. A difference in full AUC may arise mainly at
high FPRs, while a difference in pAUC_005 reflects performance under a
stringent false-positive-rate budget. The next subsection shows how restricting
the FPR range can change the comparison between analytical pipelines.
The search engine determines which proteins and quantitative evidence enter the
downstream differential-expression analysis. To measure its influence, the same
twelve injections were processed with Spectronaut and DIA-NN and evaluated
against the same expected spike-in ratios, with the downstream benchmark kept
fixed. With the full reports, that design isolates the search-engine pipeline.
The example reports shipped with NADIA, however, are trimmed to 2,000 proteins
sampled separately for each engine (only a few hundred are shared), so the
results below describe these two example datasets and do not isolate a
search-engine effect: part of any difference comes from the proteins retained.
This comparison also provides a
practical example of why both views are needed: the ranking suggested by the
full ROC range can reverse when the low-FPR region is examined.
The DIA-NN organism annotation is extracted from Protein.Names, after which
benchmarking_proteomics() can evaluate its results in exactly the same way as
the Spectronaut results.
Full-range ROC comparison. We begin with the AUC over the complete ROC range. This provides one global measure of how well each engine ranks spike-in proteins ahead of background proteins across all possible thresholds.
diann <- preprocess_diann(
system.file("extdata", "nadia_diann_report.tsv.gz", package = "NADIA"),
condition_order = c("A", "B", "D"), verbose = FALSE)
species_dn <- data.frame(
Protein.IDs = diann$protein_id$PG.ProteinGroups,
Species = toupper(sub(".*_", "",
sub(";.*", "", diann$protein_id$PG.ProteinNames))))
bench_dn <- benchmarking_proteomics(
process_proteomics(diann, verbose = FALSE)$DEPs_results,
expected, species_df = species_dn, output_dir = NULL, verbose = FALSE)
engine_auc <- rbind(
cbind(Engine = "Spectronaut",
bench$metrics_table[, c("Comparison", "TP", "FP", "TN", "FN", "AUC")]),
cbind(Engine = "DIA-NN",
bench_dn$metrics_table[, c("Comparison", "TP", "FP", "TN", "FN", "AUC")]))
engine_auc$N <- rowSums(engine_auc[, c("TP", "FP", "TN", "FN")])
engine_auc <- engine_auc[, c("Engine", "Comparison", "N", "AUC")]
knitr::kable(engine_auc, digits = 3,
caption = "Full ROC AUC by search engine and comparison.")
| Engine | Comparison | N | AUC |
|---|---|---|---|
| Spectronaut | B-A | 1997 | 0.832 |
| Spectronaut | D-A | 1997 | 0.978 |
| Spectronaut | D-B | 1997 | 0.949 |
| DIA-NN | B-A | 1965 | 0.909 |
| DIA-NN | D-A | 1965 | 0.995 |
| DIA-NN | D-B | 1965 | 0.967 |
DIA-NN has the higher full AUC in every comparison: 0.909 versus 0.832 in
B-A, 0.995 versus 0.978 in D-A, and 0.967 versus 0.949 in D-B. Thus, when
the entire ROC range is considered, DIA-NN gives the stronger overall ranking in
these example data. Both engines nevertheless show the same biological pattern, performing
best in D-A and worst in the more difficult B-A comparison.
The N column also shows that Spectronaut and DIA-NN contribute 1,997 and 1,965
proteins, respectively. AUC is based on rates, which makes it more comparable
than raw TP or FP counts across protein sets of different sizes, but it does not
remove differences in how difficult each set is to classify. The benchmark
remains conditional on the proteins each example retains; identification
coverage is a separate question.
Low-FPR comparison. Full AUC includes thresholds with large false-positive rates that would rarely be accepted in practice. We therefore repeat the comparison using pAUC at FPR limits of 1%, 5%, and 10%, focusing progressively on the region in which a proteomics analysis is more likely to operate.
engine_pauc <- cbind(
Engine = rep(c("Spectronaut", "DIA-NN"), each = 3),
rbind(bench$opdea_metrics, bench_dn$opdea_metrics)[
, c("Comparison", "pAUC_001", "pAUC_005", "pAUC_010")])
knitr::kable(engine_pauc, digits = 3,
caption = "Partial AUC by search engine and FPR range.")
| Engine | Comparison | pAUC_001 | pAUC_005 | pAUC_010 |
|---|---|---|---|---|
| Spectronaut | B-A | 0.768 | 0.807 | 0.821 |
| Spectronaut | D-A | 0.912 | 0.951 | 0.961 |
| Spectronaut | D-B | 0.868 | 0.912 | 0.925 |
| DIA-NN | B-A | 0.688 | 0.841 | 0.870 |
| DIA-NN | D-A | 0.855 | 0.961 | 0.978 |
| DIA-NN | D-B | 0.778 | 0.916 | 0.939 |
The result now depends on the FPR limit. At 1%, Spectronaut has the higher pAUC
in all three comparisons: 0.768 versus 0.688 in B-A, 0.912 versus 0.855 in
D-A, and 0.868 versus 0.778 in D-B. At 5% and 10%, the order reverses and
DIA-NN has the higher pAUC in every comparison. In these data, its greater
recovery becomes advantageous once the analysis allows a slightly larger
false-positive budget.
Choosing a search engine therefore requires an explicit decision about the acceptable FPR. A 1% limit means tolerating approximately one false positive for every 100 truly unchanged background proteins; it does not mean one false discovery per 100 reported proteins. If that strict budget is required, Spectronaut performs better here. If limits of 5% or 10% are acceptable, DIA-NN performs better. The relevant operating range, rather than a single global AUC, should determine which engine is preferred for the study.
A volcano plot displays the magnitude and statistical evidence of every
protein-level change. The x-axis shows the log2 fold change, with negative and
positive values indicating changes in opposite directions, while the y-axis
shows -log10(adj.P.Val). Proteins with stronger statistical evidence therefore
appear higher in the plot.
benchmarking_proteomics() returns one interactive Highcharts volcano plot per
comparison in hc_volcano_list. The interactive version allows individual
proteins to be examined by hovering over the points and species to be shown or
hidden from the legend.
The following code retrieves the plot for B-A:
bench$hc_volcano_list[["B-A"]]
Coloring proteins by species turns the volcano plot into a diagnostic view of
the spike-in experiment. Human background proteins are shown in gray and
should remain concentrated around log2FC = 0; yeast proteins are blue and are
expected near -0.58; and E. coli proteins are orange and are expected near
1.00. The blue and orange vertical dashed lines mark these expected spike-in
log2 fold changes, whereas the solid vertical line at zero represents no
change.
The subtitle summarizes the significant proteins for each species (restricted
to the expected direction for spike-ins). MED is the median observed log2
fold change, MAD is its median absolute deviation, and RCV is the robust
coefficient of variation (MAD / abs(MED) * 100); n is the number of
proteins summarized.
Comparing MED with the expected line assesses fold-change recovery, while
lower MAD and RCV values indicate a tighter distribution.
The horizontal dashed line marks alpha = 0.05
(adj.P.Val <= 0.05). Significant proteins have stronger colors, whereas
non-significant proteins appear lighter. In B-A, most human proteins remain
near zero and the spike-in species cluster around their expected values.
Deviations flag proteins for review. Spike-ins with the correct direction but
unusually large fold changes are often imputed or presence–absence cases,
because one condition has little or no measured signal.
After scoring one pipeline, the next question is whether another combination
of preprocessing methods recovers the known changes more reliably.
benchmarking_multiple() brings the results of separate
benchmarking_proteomics() runs into a common comparison. It does not rerun
the differential-expression analyses; it validates their benchmark summaries,
ranks the pipelines, and generates comparative plots.
The function accepts either combined tables already held in memory or a parent directory containing the files exported by each benchmark. This section begins with the in-memory workflow; the directory-based workflow is introduced in Section 12.3.
Benchmark inputs
opdea_combined: combined OpDEA table containing Assay, Comparison, nMCC, G_mean, and the three pAUC columns. This is the required input for an in-memory comparison.
confusion_combined: optional combined TP, FP, TN, and FN counts used to generate comparative confusion plots.
classified_combined: optional protein-level classified results used to draw multi-pipeline ROC curves.
bench_metrics_combined: optional classification metrics used for the extended 11-metric ranking.
results_dir: parent directory containing exported benchmark folders. It is used when opdea_combined is not supplied.
File discovery and method names
pattern: regular expression identifying OpDEA files; the default matches benchmark_opdea_metrics.tsv.
method_names: optional pipeline names for imported files. If omitted, names are taken from their parent folders.
recursive: search within subdirectories; the default is TRUE.
strip_prefix: prefix removed from imported method names; the default removes benchmark_.
Ranking and plots
metrics: metrics used for the standard ranking; by default these are nMCC, G_mean, pAUC_001, pAUC_005, and pAUC_010.
extended_metrics: metrics used when an extended ranking is requested through bench_metrics_combined.
p_col: p-value column used for multi-pipeline ROC curves; the default is adj.P.Val.
plots: “all” or any selection of ranking_heatmap_mean, ranking_heatmap_median, ranking_bars_mean, ranking_bars_median, metrics_heatmap_mean, metrics_heatmap_median, metrics_comparison, extended_ranking_bars, and extended_ranking_heatmap. Confusion and ROC plots are generated automatically when their input tables are available.
Messages and exports
verbose: print progress messages; the default is TRUE.
output_dir: optional directory for the combined tables and plots; NULL keeps them in memory.
export_plots and export_tables: control PNG and TSV export; both default to TRUE when output_dir is provided.
plot_width, plot_height, and plot_dpi: dimensions and resolution of exported plots; defaults are 12 inches, 8 inches, and 150 dpi.
The example tests three pipelines chosen to separate the effects of normalization and imputation:
| Pipeline | Normalization | Imputation |
|---|---|---|
cycloess_combo |
cycloess | Two-stage combo: Impseqrob for MAR values and min for MNAR values |
quantile_combo |
quantile | The same two-stage combo imputation |
cycloess_min |
cycloess | Minimum-value imputation applied as a single method |
Comparing the first two pipelines isolates the normalization choice, whereas
comparing cycloess_combo with cycloess_min isolates the imputation strategy.
Each pipeline is processed and scored separately; only the OpDEA and confusion
tables needed below are then combined and passed to benchmarking_multiple().
pipelines <- list(
cycloess_combo = list(norm_method = "cycloess", imp_method = "combo"),
quantile_combo = list(norm_method = "quantile", imp_method = "combo"),
cycloess_min = list(norm_method = "cycloess", imp_method = "min"))
invisible(capture.output(
scored <- lapply(names(pipelines), function(p) {
processed <- do.call(
process_proteomics,
c(list(nadia_dia, verbose = FALSE), pipelines[[p]]))
benchmark <- benchmarking_proteomics(
processed$DEPs_results, expected,
species_df = species_df,
output_dir = NULL, verbose = FALSE)
list(
opdea = cbind(Assay = p, benchmark$opdea_metrics),
confusion = cbind(Assay = p, benchmark$metrics_table))
})))
opdea <- do.call(rbind, lapply(scored, `[[`, "opdea"))
confusion <- do.call(rbind, lapply(scored, `[[`, "confusion"))
multi_bench <- benchmarking_multiple(
opdea_combined = opdea,
confusion_combined = confusion,
plots = "ranking_bars_mean",
verbose = FALSE)
The classification counts should be inspected before reducing performance to a single rank. They show whether a pipeline loses expected spike-ins as false negatives or incorrectly calls background proteins as false positives. Because all nine rows contain the same 1,997 proteins, their counts are directly comparable.
count_cols <- c("Assay", "Comparison", "TP", "FP", "TN", "FN")
knitr::kable(
multi_bench$confusion_combined[, count_cols],
caption = "Classification counts by pipeline and comparison.")
| Assay | Comparison | TP | FP | TN | FN |
|---|---|---|---|---|---|
| cycloess_combo | B-A | 469 | 27 | 1227 | 274 |
| cycloess_combo | D-A | 693 | 60 | 1194 | 50 |
| cycloess_combo | D-B | 629 | 35 | 1219 | 114 |
| quantile_combo | B-A | 494 | 35 | 1219 | 249 |
| quantile_combo | D-A | 709 | 917 | 337 | 34 |
| quantile_combo | D-B | 628 | 934 | 320 | 115 |
| cycloess_min | B-A | 400 | 27 | 1227 | 343 |
| cycloess_min | D-A | 659 | 48 | 1206 | 84 |
| cycloess_min | D-B | 569 | 28 | 1226 | 174 |
The stacked plot presents the same counts as one bar per pipeline and comparison, making the balance between correct classifications and errors easier to compare. Labels are hidden for segments smaller than 4% of a bar to avoid overlap; the table retains every exact count.
confusion_plot <- multi_bench$gg_confusion_stacked
bar_total <- ave(
confusion_plot$data$Count,
confusion_plot$data$Assay,
confusion_plot$data$Comparison,
FUN = sum)
confusion_plot$data$label_text[
confusion_plot$data$Count < 0.04 * bar_total] <- ""
confusion_plot
cycloess_combo provides the best balance: it keeps FP counts low while
recovering more TP than cycloess_min in every comparison. This is why it later
ranks first even though it does not minimize every error count individually.
quantile_combo recovers slightly more TP in B-A and D-A, but its FP counts
rise to 917 and 934 in D-A and D-B, respectively, indicating poor control
of the background. In contrast, cycloess_min is the most conservative
pipeline: it produces the fewest or joint-fewest FP, but also the most FN in all
three comparisons. The plot therefore reveals the error pattern behind the
ranking rather than only identifying a winner.
The confusion counts describe each type of error but do not provide one overall
ordering. benchmarking_multiple() therefore averages each OpDEA metric across
comparisons, ranks the pipelines from the highest metric value to the lowest,
and calculates rank_final as the mean of the five metric-specific ranks. A
lower final rank indicates better and more consistent performance.
multi_bench$mean_ranking
#> Assay rank_nMCC rank_G_mean rank_pAUC_001 rank_pAUC_005
#> 1 cycloess_combo 1 1 1 1
#> 2 cycloess_min 2 2 2 2
#> 3 quantile_combo 3 3 3 3
#> rank_pAUC_010 rank_final
#> 1 1 1
#> 2 2 2
#> 3 3 3
The result is unanimous: cycloess_combo ranks first for nMCC, G_mean, and
all three pAUC limits. cycloess_min ranks second and quantile_combo third.
The ordering reflects the trade-offs observed above: the first pipeline balances
spike-in recovery and background control, the second loses sensitivity, and the
third is strongly penalized for its false positives in D-A and D-B.
The final ranks can also be displayed as horizontal bars. This view is useful when many pipelines are compared because the methods are ordered directly and the best result remains visually identifiable.
multi_bench$gg_ranking_bars_mean
The green bar for cycloess_combo has the lowest possible final rank of 1,
followed by cycloess_min at 2 and quantile_combo at 3. The result agrees with
vignette("choosing-methods"), where cycloess was also preferred using proxy
criteria such as replicate consistency and condition separation. The two
approaches may not always select the same pipeline. When they differ, the
spike-in result is more informative because it compares each pipeline with
changes known in advance, whereas the proxy criteria evaluate only indirect
signs of good performance.
The in-memory workflow is convenient when every pipeline can be run in one R
session. A directory-based workflow is more suitable for long analyses, work
performed at different times, or pipelines implemented outside NADIA. Each
pipeline is first evaluated with benchmarking_proteomics() and exported to
its own subdirectory; benchmarking_multiple() then discovers and combines
those benchmark files through results_dir.
External differential-expression results must be benchmarked first.
benchmarking_multiple() does not score a raw MSstats or other external result
table directly. Adapt that table to the input expected by
benchmarking_proteomics(), run the individual benchmark with an
output_dir, and only then include its exported folder in the multiple
comparison.
An external differential-expression table requires one row per protein and comparison with the following fields:
| NADIA field | Required content | MSstats ComparisonResult |
|---|---|---|
Protein.IDs |
Protein identifier matching species_df |
Protein |
logFC |
Numeric log2 fold change with the same contrast direction as expected |
log2FC |
Comparison |
Comparison label matching expected$Comparison exactly |
Label |
adj.P.Val |
Numeric adjusted p-value used for scoring | adj.pvalue |
P.Value |
Optional raw p-value | pvalue |
Species can instead be included directly as a Species column, but using the
same species_df for every pipeline prevents annotation differences from
confounding the comparison. The output of
MSstats::groupComparison()
can be adapted as follows:
# msstats_result is the list returned by MSstats::groupComparison()
msstats_source <- msstats_result$ComparisonResult
msstats_de <- data.frame(
Protein.IDs = msstats_source$Protein,
logFC = msstats_source$log2FC,
Comparison = msstats_source$Label,
P.Value = msstats_source$pvalue,
adj.P.Val = msstats_source$adj.pvalue,
stringsAsFactors = FALSE)
Check the comparison labels and fold-change direction before scoring: for
example, B-A must represent the same contrast in msstats_de and expected.
After adaptation, both NADIA and external results can be exported under one
parent directory.
results_root <- "benchmark_results"
dir.create(results_root, recursive = TRUE, showWarnings = FALSE)
# Export one individual benchmark for each NADIA pipeline
for (p in names(pipelines)) {
processed <- do.call(
process_proteomics,
c(list(nadia_dia, verbose = FALSE), pipelines[[p]]))
benchmarking_proteomics(
processed$DEPs_results, expected,
species_df = species_df,
output_dir = file.path(results_root, p),
verbose = FALSE)
}
# Export the benchmark of the adapted external pipeline
benchmarking_proteomics(
msstats_de, expected,
species_df = species_df,
output_dir = file.path(results_root, "MSstats"),
verbose = FALSE)
# Import, compare, and optionally export all individual benchmarks
multi_from_files <- benchmarking_multiple(
results_dir = results_root,
output_dir = file.path(results_root, "multiple_summary"),
verbose = FALSE)
At minimum, every method folder must contain
benchmark_opdea_metrics.tsv. The additional files written automatically by
benchmarking_proteomics() extend the comparison:
benchmark_confusion_overall.tsv adds the confusion plots,
benchmark_classified.tsv adds multi-pipeline ROC curves, and
benchmark_metrics.tsv enables the extended ranking. Folder names become the
pipeline names shown in the tables and figures. All methods should use the same
experimental truth, thresholds, comparison directions, and protein universe;
if the detected protein sets differ, retain and report N alongside the
performance metrics.
A spike-in benchmark provides direct evidence of how well a pipeline recovers predefined abundance changes while preserving an unchanged background. It is therefore a strong method for comparing pipelines under the experimental conditions represented by the benchmark.
However, the highest-ranked pipeline is not guaranteed to be optimal for every
biological study. Spike-in proteins are added at controlled amounts and are
usually detected across conditions, whereas real proteins may be absent from
one condition, follow a missing-not-at-random (MNAR) pattern, or change by
amounts outside the range tested here. These situations can respond differently
to normalization, imputation, and differential-expression methods. Because
each normalization–imputation combination can produce a different downstream
result, vignette("choosing-methods") explains how to compare complete
preprocessing strategies when the true changes are unknown, while
vignette("missing-values") examines missingness and its treatment in greater
detail.
The ranking should therefore be interpreted as strong evidence about analytical performance, together with an assessment of how closely the benchmark reflects the biological experiment of interest.
Peng H, Wang H, Kong W, Li J, Goh WWB (2024). Optimizing differential expression analysis for proteomics data via high-performing rules and ensemble inference. Nature Communications 15, 3922. doi:10.1038/s41467-024-47899-w.
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] NADIA_0.99.3 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_2.0.0
#> [3] shape_1.4.6.1 magrittr_2.0.5
#> [5] magick_2.9.1 farver_2.1.2
#> [7] nloptr_2.2.1 rmarkdown_2.32
#> [9] GlobalOptions_0.1.4 vctrs_0.7.3
#> [11] memoise_2.0.1 Cairo_1.7-0
#> [13] minqa_1.2.8 janitor_2.2.1
#> [15] BiocBaseUtils_1.15.1 htmltools_0.5.9
#> [17] S4Arrays_1.13.2 curl_8.1.0
#> [19] broom_1.0.13 pROC_1.19.1
#> [21] SparseArray_1.13.4 TTR_0.24.4
#> [23] sass_0.4.10 bslib_0.12.0
#> [25] htmlwidgets_1.6.4 zoo_1.9-1
#> [27] lubridate_1.9.5 cachem_1.1.0
#> [29] lifecycle_1.0.5 iterators_1.0.14
#> [31] pkgconfig_2.0.3 Matrix_1.7-6
#> [33] R6_2.6.1 fastmap_1.2.0
#> [35] rbibutils_2.4.1 MatrixGenerics_1.25.0
#> [37] snakecase_0.11.1 clue_0.3-68
#> [39] digest_0.6.39 colorspace_2.1-4
#> [41] tidyHeatmap_1.13.1 AnnotationDbi_1.75.2
#> [43] patchwork_1.3.2 S4Vectors_0.51.10
#> [45] RSQLite_3.53.3 rrcovNA_0.5-3
#> [47] GenomicRanges_1.65.4 labeling_0.4.3
#> [49] timechange_0.4.0 mgcv_1.9-5
#> [51] httr_1.4.9 abind_1.4-8
#> [53] compiler_4.6.1 bit64_4.8.6
#> [55] withr_3.0.3 doParallel_1.0.17
#> [57] S7_0.2.2 backports_1.5.1
#> [59] BiocParallel_1.47.0 DBI_1.3.0
#> [61] viridis_0.6.5 dendextend_1.19.1
#> [63] MASS_7.3-66 DelayedArray_0.39.8
#> [65] rjson_0.2.23 tools_4.6.1
#> [67] rrcov_1.7-7 otel_0.3.0
#> [69] quantmod_0.4.29 glue_1.8.1
#> [71] nlme_3.1-172 grid_4.6.1
#> [73] cluster_2.1.8.3 generics_0.1.4
#> [75] sva_3.61.0 gtable_0.3.6
#> [77] tidyr_1.3.2 data.table_1.18.6.1
#> [79] XVector_0.53.0 BiocGenerics_0.59.13
#> [81] foreach_1.5.2 pillar_1.11.1
#> [83] stringr_1.6.0 limma_3.99.4
#> [85] genefilter_1.95.0 logging_0.10-111
#> [87] robustbase_0.99-7 circlize_0.4.18
#> [89] splines_4.6.1 dplyr_1.2.1
#> [91] lattice_0.23-1 survival_3.8-12
#> [93] bit_4.6.0 annotate_1.91.1
#> [95] tidyselect_1.2.1 ComplexHeatmap_2.29.0
#> [97] Biostrings_2.81.9 knitr_1.52
#> [99] BERT_1.9.0 reformulas_0.4.4
#> [101] gridExtra_2.3.1 bookdown_0.48
#> [103] IRanges_2.47.5 Seqinfo_1.3.2
#> [105] edgeR_4.99.12 SummarizedExperiment_1.43.0
#> [107] stats4_4.6.1 xfun_0.62
#> [109] Biobase_2.73.2 statmod_1.5.2
#> [111] matrixStats_1.5.0 DEoptimR_1.2-2
#> [113] stringi_1.8.9 yaml_2.3.12
#> [115] boot_1.3-32 evaluate_1.0.5
#> [117] codetools_0.2-20 MsCoreUtils_1.25.4
#> [119] tibble_3.3.1 BiocManager_1.30.27
#> [121] cli_3.6.6 xtable_1.8-8
#> [123] Rdpack_2.6.6 jquerylib_0.1.4
#> [125] dichromat_2.0-1 Rcpp_1.1.2
#> [127] norm_1.0-11.1 png_0.1-9
#> [129] XML_3.99-0.25 parallel_4.6.1
#> [131] blob_1.3.0 ggplot2_4.0.3
#> [133] assertthat_0.2.1 lme4_2.1-0
#> [135] viridisLite_0.4.3 mvtnorm_1.4-2
#> [137] rlist_0.4.6.2 scales_1.4.0
#> [139] xts_0.14.3 pcaPP_2.0-7
#> [141] purrr_1.2.2 highcharter_0.9.5
#> [143] crayon_1.5.3 GetoptLong_1.1.1
#> [145] rlang_1.3.0 KEGGREST_1.53.6