diagFDR quick start: scope, calibration, stability and entrapment on a simulated list

This vignette runs end to end without external data. It uses the ground-truth simulator of the diagFDR paper (Supplementary Note S6) so that every diagnostic can be compared with the realised false discovery proportion, which is known exactly in simulation and unknown on real data.

library(diagFDR)
sim <- sim_tdc_example          # 20,000 units, Gaussian nulls, correct specification
sim
#> <dfdr_tbl>
#>   unit:    simulated_unit
#>   scope:   global
#>   q_from:  tdc_counting
#>   n:       20000
#> # A tibble: 20,000 × 6
#>    id       is_decoy      q           pep run   score
#>    <chr>    <lgl>     <dbl>         <dbl> <chr> <dbl>
#>  1 u0000001 TRUE     0.352  NA            <NA>  0.391
#>  2 u0000002 FALSE    0       0.0000000448 <NA>  6.11 
#>  3 u0000003 FALSE    0       0.000197     <NA>  4.31 
#>  4 u0000004 TRUE     0.238  NA            <NA>  0.935
#>  5 u0000005 FALSE    0.190   0.998        <NA>  1.14 
#>  6 u0000006 TRUE     0.211  NA            <NA>  1.06 
#>  7 u0000007 TRUE     0.116  NA            <NA>  1.50 
#>  8 u0000008 TRUE     0.0324 NA            <NA>  2.13 
#>  9 u0000009 FALSE    0.359   1.000        <NA>  0.344
#> 10 u0000010 FALSE    0.330   1.000        <NA>  0.500
#> # ℹ 19,990 more rows

1. Headline support and the FDP-CV benchmark

An FDR claim at a threshold \(\alpha\) should come with the accepted target and decoy counts, \(T_\alpha\) and \(D_\alpha\), the recomputed decoy ratio \(D_\alpha/T_\alpha\), the decoy-support index \(D_\alpha^{-1/2}\) and the FDP-CV benchmark \(\sqrt2\,D_\alpha^{-1/2}\) (Proposition 1).

h <- dfdr_headline(sim, alpha = 0.01)
h[, c("alpha", "T_alpha", "D_alpha", "FDR_hat", "support_index", "CV_benchmark")]
#> # A tibble: 1 × 6
#>   alpha T_alpha D_alpha FDR_hat support_index CV_benchmark
#>   <dbl>   <int>   <int>   <dbl>         <dbl>        <dbl>
#> 1  0.01    6565      65 0.00990         0.124        0.175

In simulation the realised FDP is known:

dfdr_realised_fdp(sim, alpha = c(0.05, 0.01, 0.002))
#>   alpha T_alpha D_alpha V_alpha decoy_ratio         FDP
#> 1 0.050    7005     350     395 0.049964311 0.056388294
#> 2 0.010    6565      65      74 0.009900990 0.011271896
#> 3 0.002    6308      12      14 0.001902346 0.002219404

The benchmark says how much the realised FDP is expected to vary across experiments with this decoy support. dfdr_check_cv_benchmark() checks it by Monte Carlo (the paper uses 1000 replicates of \(2\times10^5\) units; here a small run suffices):

dfdr_check_cv_benchmark(B = 30, n = 20000, alpha = 0.01)[, c("alpha", "CV_FDP", "mean_D", "kappa")]
#>   alpha    CV_FDP mean_D    kappa
#> 1  0.01 0.1890189   65.4 1.528541

kappa is the ratio of the observed CV to \(D_\alpha^{-1/2}\); Proposition 1 predicts \(\sqrt2\approx1.41\) under independence.

2. Stability across thresholds

alphas <- c(0.0005, 0.001, 0.002, 0.005, 0.01, 0.02, 0.05)
dfdr_curve_stability(sim, alphas = alphas)[, c("alpha", "T_alpha", "D_alpha", "support_index")]
#> # A tibble: 7 × 4
#>    alpha T_alpha D_alpha support_index
#>    <dbl>   <int>   <int>         <dbl>
#> 1 0.0005    6065       3        0.577 
#> 2 0.001     6186       6        0.408 
#> 3 0.002     6308      12        0.289 
#> 4 0.005     6461      32        0.177 
#> 5 0.01      6565      65        0.124 
#> 6 0.02      6705     134        0.0864
#> 7 0.05      7005     350        0.0535
dfdr_local_tail(sim, alphas = alphas, win_rel = 0.2)[, c("alpha", "D_alpha_win")]
#> # A tibble: 7 × 2
#>    alpha D_alpha_win
#>    <dbl>       <int>
#> 1 0.0005           0
#> 2 0.001            1
#> 3 0.002            3
#> 4 0.005            6
#> 5 0.01            14
#> 6 0.02            28
#> 7 0.05            75
dfdr_elasticity(sim, alphas = alphas, eps = 0.2)[, c("alpha", "jaccard")]
#> # A tibble: 7 × 2
#>    alpha jaccard
#>    <dbl>   <dbl>
#> 1 0.0005   1    
#> 2 0.001    0.995
#> 3 0.002    0.997
#> 4 0.005    0.998
#> 5 0.01     0.997
#> 6 0.02     0.992
#> 7 0.05     0.989

3. Calibration: PEP concordance and equal chance

The simulated list carries the true PEP of every target, so the IPE should be close to zero under the target–decoy comparator \(D_b/T_b\):

dfdr_ipe(sim)$IPE
#> [1] 0.002067698
dfdr_equal_chance_qbands(sim)$pooled
#> # A tibble: 1 × 11
#>   qmax_export low_lo low_hi N_test N_D_test pi_D_hat effect_abs ci95_lo ci95_hi
#>         <dbl>  <dbl>  <dbl>  <int>    <int>    <dbl>      <dbl>   <dbl>   <dbl>
#> 1           1    0.2    0.5   9973     4968    0.498    0.00186   0.488   0.508
#> # ℹ 2 more variables: p_value_binom <dbl>, pass_minN <lgl>

4. What a misspecified decoy model looks like

Decoys that are under-dispersed (scenario S2 of the paper) inflate the realised FDP while the reported decoy ratio stays at \(\alpha\):

bad <- dfdr_simulate_tdc(n = 20000, decoy_scale = 0.85, seed = 2)
dfdr_realised_fdp(bad, alpha = 0.01)
#>   alpha T_alpha D_alpha V_alpha decoy_ratio        FDP
#> 1  0.01    6877      68     194 0.009888033 0.02820998

This is why external validation (entrapment) is a distinct requirement.

5. Label propagation at protein-group level

When several exported rows collapse into one reported unit and the rows do not all carry the same target/decoy label, the rule that assigns the unit’s label must be stated. The four rules are compared by dfdr_label_rule_audit():

set.seed(1)
rep <- data.frame(Protein.Group = rep(paste0("P", 1:2000), each = 3),
                  Decoy = rbinom(6000, 1, 0.12),
                  Global.PG.Q.Value = c(runif(6000, 0, 0.3)))
label_is_intrinsic(rep, "Protein.Group")
#> $intrinsic
#> [1] FALSE
#> 
#> $n_mixed
#> [1] 678
#> 
#> $n_keys
#> [1] 2000
dfdr_label_rule_audit(rep, "Protein.Group", "Global.PG.Q.Value", alpha = 0.05)
#>          rule n_units T_alpha D_alpha decoy_ratio support_index CV_benchmark
#> 1       split    2678     740     125  0.16891892    0.08944272    0.1264911
#> 2 winning_row    2000     728     111  0.15247253    0.09491580    0.1342312
#> 3   any_decoy    2000     544     295  0.54227941    0.05822225    0.0823387
#> 4    majority    2000     807      32  0.03965304    0.17677670    0.2500000
#>   n_mixed
#> 1     678
#> 2     678
#> 3     678
#> 4     678

Picked-protein competition is the procedural alternative that needs no rule:

pg <- diann_global_proteingroup(rep, rule = "split", two_pass = FALSE)
pk <- dfdr_picked_protein(pg, alpha = 0.05)
c(T = pk$T_alpha, D = pk$D_alpha, ratio = pk$decoy_ratio, pairs = pk$n_pairs)
#>     T     D ratio pairs 
#>    19     0     0   678

6. External validation by entrapment

Entrapment matches are accepted identifications known to be false because they map only to a proteome absent from the sample. On a real DIA-NN entrapment search the ids come from diann_entrapment_ids(rep, accessions) with the accessions read by read_fasta_accessions(); here we label some incorrect targets as entrapment to show the output.

truth <- attr(sim, "truth")
ent_ids <- sim$id[!sim$is_decoy & truth$is_null & runif(nrow(sim)) < 0.3]
dfdr_entrapment(sim, ent_ids, alpha = c(0.05, 0.01), r = 0.43)
#>   alpha T_alpha D_alpha decoy_ratio E_alpha     fdp_raw          lo          hi
#> 1  0.05    7005     350  0.04996431     122 0.017416131 0.014606662 0.020754597
#> 2  0.01    6565      65  0.00990099      26 0.003960396 0.002704189 0.005796771
#>      r fdp_corrected fdp_corrected_target_only
#> 1 0.43    0.05791876               0.041220533
#> 2 0.43    0.01317062               0.009246845

fdp_raw is \(E_\alpha/T_\alpha\) with a Wilson interval; fdp_corrected is the size-corrected estimate of Wen et al. (2025), which needs the entrapment-to-target size ratio r (dfdr_entrapment_size_ratio() computes it from the FASTA files at protein and peptide level).

Two scope choices that overlap almost entirely can still differ sharply where they disagree. dfdr_entrapment_discordant() computes the entrapment rate on the shared core and on the discordant set:

a1 <- dfdr_accepted(sim, 0.01)
a2 <- dfdr_accepted(sim, 0.03)      # a looser list standing in for a mis-scoped one
d  <- dfdr_entrapment_discordant(a1, a2, ent_ids)
d
#>          set    n  E         fdp          lo          hi
#> 1       core 6565 26 0.003960396 0.002704189 0.005796771
#> 2      added  235 40 0.170212766 0.127560688 0.223473255
#> 3    removed    0  0          NA          NA          NA
#> 4 discordant  235 40 0.170212766 0.127560688 0.223473255
#> 5      list1 6565 26 0.003960396 0.002704189 0.005796771
#> 6      list2 6800 66 0.009705882 0.007636771 0.012328635
attr(d, "ratio"); attr(d, "jaccard")
#> [1] 42.97872
#> [1] 0.9654412

7. Report

dfdr_run_all() runs the standard suite on a named list of dfdr_tbl objects and dfdr_render_report() writes the HTML report with alerts tagged by their evidential basis; see the DIA-NN, Spectronaut, MaxQuant, mokapot and mzIdentML vignettes for the adapters.