Delta Analysis in PhIP-seq

Overview

The Delta framework answers two related but distinct questions about peptide-level prevalence:

Where do peptides shift? — deltaplot() and deltaplot_interactive() place each peptide at its pooled prevalence (x) and prevalence shift Δ (y), providing a global view of which features move between groups.
Is there a statistically significant global shift? — compute_delta() aggregates per-peptide z-scores into a single Stouffer-type statistic and assesses significance via subject-label permutation. forestplot() and forestplot_interactive() then display the top and bottom features ranked by that statistic.

This vignette walks through the complete Delta workflow in phiper:

Visualising raw prevalence shift with deltaplot() and deltaplot_interactive()
Running the permutation test with compute_delta()
Displaying results with forestplot() and forestplot_interactive()
Plotting prevalence along proteins with profile_plot()

Setup

library(phiper)

Load the bundled example dataset. It contains two patient groups (A, B) measured at two timepoints (T1, T2) across 1 000 simulated peptides.

pd <- load_example_data()
pd
#> ── <phip_data> ───────────────────────────────────────────────────────────────── 
#> 
#> counts (first 5 rows): 
#> # A tibble: 5 × 9
#>   sample_id subject_id group timepoint peptide_id     exist counts_control
#>   <chr>     <chr>      <chr> <chr>     <chr>          <int>          <int>
#> 1 A_T1_1    A_1        A     T1        agilent_102300     0              7
#> 2 A_T1_1    A_1        A     T1        agilent_102431     0             25
#> 3 A_T1_1    A_1        A     T1        agilent_103785     0              4
#> 4 A_T1_1    A_1        A     T1        agilent_103792     0             35
#> 5 A_T1_1    A_1        A     T1        agilent_10384      0              0
#> # ℹ 2 more variables: counts_hits <int>, fold_change <dbl>
#> 
#> table size: 80,000 rows x 9 columns
#> 
#> peptide library preview (first 5 rows): 
#> # A tibble: 5 × 8
#>   peptide_id Fullname                    species genus family order class common
#>   <chr>      <chr>                       <chr>   <chr> <chr>  <chr> <chr> <chr> 
#> 1 agilent_1  Chromodomain-helicase-DNA-… Homo s… Homo  Homin… Prim… Mamm… Human 
#> 2 agilent_2  integral membrane protein   Mycoba… Myco… Mycob… Myco… Acti… <NA>  
#> 3 agilent_3  hypothetical protein (6/16… Mycoba… Myco… Mycob… Myco… Acti… <NA>  
#> 4 agilent_4  envelope protein (5/8) & a… Orthof… Orth… Flavi… Amar… Flas… JEV   
#> 5 agilent_5  Myosin-7 & beta-myosin hea… Homo s… Homo  Homin… Prim… Mamm… Human 
#> ... plus 40 more columns
#> 
#> library size: 357,190 rows x 48 columns
#> 
#> meta flags: 
#>   con:            <duckdb_connection>
#>   peptide_libraries: combined
#>   longitudinal:   TRUE
#>   exist:          TRUE
#>   fold_change:    TRUE
#>   raw_counts:     FALSE
#>   extra_cols:     group, counts_control, counts_hits
#>   peptide_con:    <duckdb_connection>
#>   materialise_table: TRUE
#>   finalizer_env:  <environment>
#>   full_cross:     FALSE

Note. The example data are entirely simulated and have no biological meaning. They exist solely to demonstrate the API.

The Delta plots consume prevalence tables produced by compute_pop() (see the POP vignette for details). Compute a group-level prevalence table first:

pop_group <- compute_pop(
  pd,
  rank_cols  = "peptide_id",
  group_cols = "group"
)
#> [13:58:51] INFO  compute_pop
#> [13:58:51] INFO  compute_pop
#>                    - ranks : peptide_id
#>                    - group_cols: group
#>                    - exist_col : exist
#>                    - pop_k_min : 1
#>                    - paired : FALSE
#> [13:58:51] INFO  ranks resolved
#>                    - available: peptide_id
#> [13:58:52] INFO  computing cohort sizes and validating binary group_cols
#> [13:58:52] INFO  computing presence per sample via k-of-n rule
#> [13:58:52] INFO  counting present samples per feature (pop, unpaired)
#> [13:58:52] INFO  building pairwise comparisons
#> [13:58:53] OK    materialized; computing Fisher p-values
#>                    - table: ph_pop_20261019_135852
#> [13:58:53] OK    done (compute_pop, unpaired)
#>                    - rows : 1950
#>                    - ranks : peptide_id
#>                    - k_min : 1
#> [13:58:53] OK    compute_pop - done
#>                  -> elapsed: 1.863s
head(pop_group)
#>         rank        feature group_col group1 n1 N1 prop1 percent1 group2 n2 N2
#> 1 peptide_id agilent_100642     group      A  0 38     0        0      B  8 42
#> 2 peptide_id agilent_100997     group      A  0 38     0        0      B  6 42
#> 3 peptide_id  agilent_10133     group      A  0 38     0        0      B  5 42
#> 4 peptide_id agilent_101516     group      A  0 38     0        0      B 17 42
#> 5 peptide_id agilent_101615     group      A  0 38     0        0      B  9 42
#> 6 peptide_id agilent_101998     group      A  0 38     0        0      B 11 42
#>       prop2 percent2      ratio delta_ratio        p_raw n_peptides
#> 1 0.1904762 19.04762 0.06907895   -6.238095 5.758809e-03          1
#> 2 0.1428571 14.28571 0.09210526   -4.428571 2.664380e-02          1
#> 3 0.1190476 11.90476 0.11052632   -3.523810 5.626494e-02          1
#> 4 0.4047619 40.47619 0.03250774  -14.380952 2.792824e-06          1
#> 5 0.2142857 21.42857 0.06140351   -7.142857 2.625712e-03          1
#> 6 0.2619048 26.19048 0.05023923   -8.952381 5.233874e-04          1

For the unpaired group comparison with compute_delta(), each subject must appear at most once per group. Because pd contains two timepoints per subject, restrict to a single timepoint first:

pd_t1 <- pd |> dplyr::filter(timepoint == "T1") |> dplyr::collect()

Visualising prevalence shift

`deltaplot()` — static ggplot2

deltaplot() plots each peptide at:

x-axis — pooled prevalence (prop1 + prop2) / 2
y-axis — prevalence shift prop2 − prop1

A dashed horizontal line marks Δ = 0. Arrows and labels indicate the direction of enrichment. A GAM smooth (controlled by add_smooth) summarises the trend across pooled prevalence.

The only required input is a data frame containing group1, group2, prop1, and prop2 — which is exactly what compute_pop() returns.

deltaplot(
  pop_group,
  group_pair_values = c("A", "B"),
  group_labels      = c("Group A", "Group B")
)
#> [13:58:54] INFO  Preparing delta prevalence plot.

Use group_pair_values to select a specific contrast when the table contains more than one pair. group_labels provides display names for the axis annotations.

Adjusting the smooth

The smooth uses a GAM with basis dimension smooth_k. Lower values produce smoother curves; set add_smooth = FALSE to suppress it.

deltaplot(
  pop_group,
  group_pair_values  = c("A", "B"),
  group_labels       = c("Group A", "Group B"),
  smooth_k           = 3,
  point_alpha        = 0.15,
  point_size         = 0.4
)
#> [13:58:54] INFO  Preparing delta prevalence plot.

Jitter and appearance

deltaplot(
  pop_group,
  group_pair_values   = c("A", "B"),
  group_labels        = c("Group A", "Group B"),
  point_jitter_width  = 0.01,
  point_jitter_height = 0.01,
  point_alpha         = 0.20,
  arrow_color         = "steelblue"
)
#> [13:58:54] INFO  Preparing delta prevalence plot.

`deltaplot_interactive()` — plotly

The interactive version mirrors the static API. Hovering over a point shows the feature identifier, group-level counts and percentages, and p-value.

deltaplot_interactive(
  pop_group,
  group_pair_values = c("A", "B"),
  group_labels      = c("Group A", "Group B")
)

Disable the smooth if you have few peptides or want a cleaner display:

deltaplot_interactive(
  pop_group,
  group_pair_values = c("A", "B"),
  group_labels      = c("Group A", "Group B"),
  add_smooth        = FALSE,
  point_size        = 5,
  point_alpha       = 0.5
)

Testing global shift: `compute_delta()`

compute_delta() performs a permutation test for a global (antigen-wide) shift in prevalence between two groups. For each (rank, feature) stratum it computes a per-peptide z-score, aggregates them via a Stouffer statistic, and assesses significance by permuting subject labels.

Unpaired design — group comparison

set.seed(42)
delta_group <- compute_delta(
  pd_t1,
  rank_cols      = "peptide_id",
  group_cols     = "group",
  B_permutations = 2000L
)

delta_group
#> # A tibble: 1,950 × 19
#>    rank       feature        group_col group1 group2 design   n_subjects_paired
#>    <chr>      <chr>          <chr>     <chr>  <chr>  <chr>                <int>
#>  1 peptide_id agilent_100642 group     A      B      unpaired                NA
#>  2 peptide_id agilent_100997 group     A      B      unpaired                NA
#>  3 peptide_id agilent_10133  group     A      B      unpaired                NA
#>  4 peptide_id agilent_101516 group     A      B      unpaired                NA
#>  5 peptide_id agilent_101615 group     A      B      unpaired                NA
#>  6 peptide_id agilent_101998 group     A      B      unpaired                NA
#>  7 peptide_id agilent_102297 group     A      B      unpaired                NA
#>  8 peptide_id agilent_102300 group     A      B      unpaired                NA
#>  9 peptide_id agilent_102431 group     A      B      unpaired                NA
#> 10 peptide_id agilent_102806 group     A      B      unpaired                NA
#> # ℹ 1,940 more rows
#> # ℹ 12 more variables: n_peptides_used <int>, m_eff <dbl>, T_obs <dbl>,
#> #   T_null_mean <dbl>, T_null_sd <dbl>, T_obs_stand <dbl>, Z_from_p <dbl>,
#> #   p_perm <dbl>, b <int>, max_delta <dbl>, frac_delta_pos <dbl>,
#> #   frac_delta_pos_w <dbl>

What’s in the output?

Column	Description
`rank`	Rank at which the feature is defined
`feature`	Feature identifier
`group_col`	Name of the grouping column
`group1`, `group2`	The two groups compared
`design`	`"paired"` or `"unpaired"`
`n_subjects_paired`	Paired subjects used (or `NA` for unpaired)
`n_peptides_used`	Peptides contributing to the test
`m_eff`	Effective number of peptides after weighting
`T_obs`	Observed Stouffer-type test statistic
`p_perm`	Two-sided permutation p-value
`b`	Number of permutations actually used
`max_delta`	Maximum absolute peptide-level prevalence difference
`frac_delta_pos`	Unweighted fraction of positive peptide-level deltas
`frac_delta_pos_w`	Weighted fraction of positive peptide-level deltas

Choosing the per-peptide statistic

stat_mode controls how each peptide’s z-score is computed:

`stat_mode`	Description
`"srlr"` (default)	Signed root likelihood ratio
`"diff"`	Wald-type z: Δ / SE
`"asin"`	Arcsine-transformed z (variance-stabilised)
`"score"`	Pooled score z (2×2 score test)
`"mcnemar"`	Signed McNemar z — paired designs only
`"srlr_paired"`	Signed root deviance — paired designs only

# Arcsine-transformed statistic
delta_asin <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "group",
  stat_mode      = "asin",
  B_permutations = 2000L
)

Aggregation and stratification

aggregate_stat controls how peptide-level z-scores are combined:

"stouffer" (default) — weighted Stouffer combination
"maxmean" — mean of the dominant direction
"af" — absolute fraction

strat_bins bins peptides by pooled prevalence before aggregation. The default bins c(0.002, 0.005, ..., 0.50) stratify across the prevalence range so that low- and high-prevalence peptides are treated separately. Set strat_bins = 0 for a single global statistic without stratification.

# No stratification
delta_nostrat <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "group",
  strat_bins     = 0,
  B_permutations = 2000L
)

# Custom bins
delta_custom <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "group",
  strat_bins     = c(0.05, 0.25, 0.50),
  B_permutations = 2000L
)

Weighting

weight_mode controls peptide weights in the Stouffer statistic:

"equal" (default) — all peptides weight 1
"se_invvar" — inverse-SE weights (downweights noisy peptides)
"n_eff_sqrt" — square root of expected positives (emphasises common peptides)

delta_invvar <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "group",
  weight_mode    = "se_invvar",
  B_permutations = 2000L
)

Paired design — timepoint comparison

When the same subjects are measured at two timepoints, set paired_by to the subject-linking column and stat_mode to a paired statistic. The function then performs sign-flipping permutations instead of label shuffling.

set.seed(42)
delta_paired <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "timepoint",
  paired_by      = "subject_id",
  stat_mode      = "srlr_paired",
  B_permutations = 2000L
)
delta_paired
#> # A tibble: 1,945 × 19
#>    rank       feature        group_col group1 group2 design n_subjects_paired
#>    <chr>      <chr>          <chr>     <chr>  <chr>  <chr>              <int>
#>  1 peptide_id agilent_100642 timepoint T1     T2     paired                37
#>  2 peptide_id agilent_100997 timepoint T1     T2     paired                37
#>  3 peptide_id agilent_10133  timepoint T1     T2     paired                37
#>  4 peptide_id agilent_101516 timepoint T1     T2     paired                37
#>  5 peptide_id agilent_101615 timepoint T1     T2     paired                37
#>  6 peptide_id agilent_101998 timepoint T1     T2     paired                37
#>  7 peptide_id agilent_102297 timepoint T1     T2     paired                37
#>  8 peptide_id agilent_102300 timepoint T1     T2     paired                37
#>  9 peptide_id agilent_102431 timepoint T1     T2     paired                37
#> 10 peptide_id agilent_102806 timepoint T1     T2     paired                37
#> # ℹ 1,935 more rows
#> # ℹ 12 more variables: n_peptides_used <int>, m_eff <dbl>, T_obs <dbl>,
#> #   T_null_mean <dbl>, T_null_sd <dbl>, T_obs_stand <dbl>, Z_from_p <dbl>,
#> #   p_perm <dbl>, b <int>, max_delta <dbl>, frac_delta_pos <dbl>,
#> #   frac_delta_pos_w <dbl>

Forest plots

Forest plots display the most extreme positive and negative features for a chosen rank, sorted by the selected test statistic.

`forestplot()` — static ggplot2

forestplot(
  delta_group,
  rank_of_interest  = "peptide_id",
  statistic_to_plot = "T",
  n_neg_each        = 10,
  n_pos_each        = 10,
  left_label        = "More in A",
  right_label       = "More in B"
)$plot

statistic_to_plot selects what is plotted on the x-axis:

Value	Description
`"T"`	Raw `T_obs` from `compute_delta()`
`"T_stand"`	Permutation-standardised T (requires `T_obs_stand` column)
`"Z_from_p"`	Signed Z derived from permutation p-values (requires `Z_from_p`)

Diverging colours

Set use_diverging_colors = TRUE to shade bars from blue (negative) to red (positive):

forestplot(
  delta_group,
  rank_of_interest     = "peptide_id",
  statistic_to_plot    = "T",
  n_neg_each           = 10,
  n_pos_each           = 10,
  left_label           = "More in A",
  right_label          = "More in B",
  use_diverging_colors = TRUE
)$plot

Filtering to significant features

Pass a column name to filter_significant to restrict the plot to features below a given significance threshold. With p_perm this retains only permutation-significant features.

forestplot(
  delta_group,
  rank_of_interest   = "peptide_id",
  statistic_to_plot  = "T",
  filter_significant = "p_perm",
  sig_level          = 0.05,
  left_label         = "More in A",
  right_label        = "More in B"
)$plot

`forestplot_interactive()` — plotly

The interactive version mirrors the static API and returns a plotly widget. Hovering over a bar shows the feature identifier and statistic value.

forestplot_interactive(
  delta_group,
  rank_of_interest     = "peptide_id",
  statistic_to_plot    = "T",
  n_neg_each           = 10,
  n_pos_each           = 10,
  left_label           = "More in A",
  right_label          = "More in B",
  use_diverging_colors = TRUE
)$plot

Profiles along a protein: `profile_plot()`

forestplot() shows which features shift; profile_plot() shows where inside a protein the reactivity sits. Each library peptide of a protein is drawn as a segment over the amino acids it covers, at its prevalence in each group, with one panel per protein. The positions come from the phiperio peptide library the peptide IDs belong to. Peptides of the protein that are missing from the compute_pop() table, because they are present in no sample, are drawn at 0%.

pop_t1 <- compute_pop(pd_t1, rank_cols = "peptide_id", group_cols = "group")
#> [14:02:16] INFO  compute_pop
#> [14:02:16] INFO  compute_pop
#>                    - ranks : peptide_id
#>                    - group_cols: group
#>                    - exist_col : exist
#>                    - pop_k_min : 1
#>                    - paired : FALSE
#> [14:02:16] INFO  ranks resolved
#>                    - available: peptide_id
#> [14:02:16] INFO  computing cohort sizes and validating binary group_cols
#> [14:02:16] INFO  computing presence per sample via k-of-n rule
#> [14:02:16] INFO  counting present samples per feature (pop, unpaired)
#> [14:02:16] INFO  building pairwise comparisons
#> [14:02:17] OK    materialized; computing Fisher p-values
#>                    - table: ph_pop_20261019_140216
#> [14:02:17] OK    done (compute_pop, unpaired)
#>                    - rows : 1950
#>                    - ranks : peptide_id
#>                    - k_min : 1
#> [14:02:17] OK    compute_pop - done
#>                  -> elapsed: 1.282s

proteins <- c("protseq_9788", "protseq_25936")
profile_plot(pop_t1, proteins = proteins, ncol = 2)
#> [14:02:17] INFO  Retrieving peptide metadata into DuckDB cache
#>                  -> get_peptide_library(library = combined, force_refresh =
#>                     FALSE)
#> [14:02:17] INFO  Opened DuckDB connection
#>                    - cache dir:
#>                      /github/home/.cache/R/phiperio/peptide_meta/phip_cache.duckdb
#>                    - tables: peptide_meta_combined
#> [14:02:17] OK    Using cached peptide_meta_combined (fast path)
#> [14:02:17] OK    Retrieving peptide metadata into DuckDB cache - done
#>                  -> elapsed: 0.019s

With show_overlap = TRUE (the default), a shaded band marks the longest run of overlapping hit peptides per protein, group pair and direction, in the colour of the group in which they are more common. Hit peptides have a BH-adjusted p_raw of at most alpha (default 0.05), corrected over all peptides of the pair. The simulated example data has no overlapping hit peptides, so no band appears here; the examples in ?profile_plot show one.

A table with several group pairs gives one line per group. Here the prevalence by patient group and by timepoint, over both timepoints:

pop_two <- compute_pop(
  pd,
  rank_cols  = "peptide_id",
  group_cols = c("group", "timepoint")
)
#> [14:02:17] INFO  compute_pop
#> [14:02:17] INFO  compute_pop
#>                    - ranks : peptide_id
#>                    - group_cols: group, timepoint
#>                    - exist_col : exist
#>                    - pop_k_min : 1
#>                    - paired : FALSE
#> [14:02:17] INFO  ranks resolved
#>                    - available: peptide_id
#> [14:02:17] INFO  computing cohort sizes and validating binary group_cols
#> [14:02:18] INFO  computing presence per sample via k-of-n rule
#> [14:02:18] INFO  counting present samples per feature (pop, unpaired)
#> [14:02:18] INFO  building pairwise comparisons
#> [14:02:19] OK    materialized; computing Fisher p-values
#>                    - table: ph_pop_20261019_140218
#> [14:02:20] OK    done (compute_pop, unpaired)
#>                    - rows : 3900
#>                    - ranks : peptide_id
#>                    - k_min : 1
#> [14:02:20] OK    compute_pop - done
#>                  -> elapsed: 2.389s
profile_plot(pop_two, proteins = proteins, ncol = 2)
#> [14:02:20] INFO  Retrieving peptide metadata into DuckDB cache
#>                  -> get_peptide_library(library = combined, force_refresh =
#>                     FALSE)
#> [14:02:20] INFO  Opened DuckDB connection
#>                    - cache dir:
#>                      /github/home/.cache/R/phiperio/peptide_meta/phip_cache.duckdb
#>                    - tables: peptide_meta_combined
#> [14:02:20] OK    Using cached peptide_meta_combined (fast path)
#> [14:02:20] OK    Retrieving peptide metadata into DuckDB cache - done
#>                  -> elapsed: 0.018s

Names of proteins and groups are used as panel titles and legend labels. delta_tbl adds the statistic and permutation p-value of compute_delta(), run at the protein level, to the panel titles:

set.seed(42)
delta_protein <- compute_delta(
  pd_t1,
  rank_cols         = "protein_seq_id",
  group_cols        = "group",
  rank_feature_keep = list(protein_seq_id = proteins),
  B_permutations    = 2000L
)
#> [14:02:20] INFO  Retrieving peptide metadata into DuckDB cache
#>                  -> get_peptide_library(library = combined, force_refresh =
#>                     FALSE)
#> [14:02:20] INFO  Opened DuckDB connection
#>                    - cache dir:
#>                      /github/home/.cache/R/phiperio/peptide_meta/phip_cache.duckdb
#>                    - tables: peptide_meta_combined
#> [14:02:20] OK    Using cached peptide_meta_combined (fast path)
#> [14:02:20] OK    Retrieving peptide metadata into DuckDB cache - done
#>                  -> elapsed: 0.018s

profile_plot(
  pop_t1,
  proteins  = c("IgA1 protease" = "protseq_9788", "Desmoglein-3" = "protseq_25936"),
  groups    = c("Group A" = "A", "Group B" = "B"),
  delta_tbl = delta_protein,
  ncol      = 2
)
#> [14:02:20] INFO  Retrieving peptide metadata into DuckDB cache
#>                  -> get_peptide_library(library = combined, force_refresh =
#>                     FALSE)
#> [14:02:20] INFO  Opened DuckDB connection
#>                    - cache dir:
#>                      /github/home/.cache/R/phiperio/peptide_meta/phip_cache.duckdb
#>                    - tables: peptide_meta_combined
#> [14:02:20] OK    Using cached peptide_meta_combined (fast path)
#> [14:02:20] OK    Retrieving peptide metadata into DuckDB cache - done
#>                  -> elapsed: 0.018s

Putting it all together

A typical Delta analysis runs in four steps:

# 1. Compute prevalence (needed for deltaplot)
pop <- compute_pop(
  pd,
  rank_cols  = "peptide_id",
  group_cols = "group"
)

# 2. Visualise raw shift
deltaplot(pop, group_pair_values = c("A", "B"),
          group_labels = c("Group A", "Group B"))
deltaplot_interactive(pop, group_pair_values = c("A", "B"),
                      group_labels = c("Group A", "Group B"))

# 3. Permutation test (unpaired)
set.seed(42)
delta <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "group",
  B_permutations = 2000L
)

# Optional: paired test across timepoints
delta_paired <- compute_delta(
  pd,
  rank_cols      = "peptide_id",
  group_cols     = "timepoint",
  paired_by      = "subject_id",
  stat_mode      = "srlr_paired",
  B_permutations = 2000L
)

# 4. Forest plots
forestplot(delta, rank_of_interest = "peptide_id",
           left_label = "More in A", right_label = "More in B")$plot
forestplot_interactive(delta, rank_of_interest = "peptide_id",
                       left_label = "More in A",
                       right_label = "More in B")$plot

Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 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.32.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8    
#>  [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: Etc/UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] dplyr_1.2.1    phiper_0.4.23  rmarkdown_2.32
#> 
#> loaded via a namespace (and not attached):
#>  [1] future_1.80.0       sass_0.4.10         utf8_1.2.6         
#>  [4] generics_0.1.4      tidyr_1.3.2         lattice_0.23-1     
#>  [7] listenv_1.1.0       digest_0.6.39       magrittr_2.0.5     
#> [10] evaluate_1.0.5      grid_4.6.1          RColorBrewer_1.1-3 
#> [13] fastmap_1.2.0       blob_1.3.0          Matrix_1.7-6       
#> [16] jsonlite_2.0.0      DBI_1.3.0           mgcv_1.9-4         
#> [19] purrr_1.2.2         scales_1.4.0        permute_0.9-10     
#> [22] codetools_0.2-20    jquerylib_0.1.4     duckdb_1.5.6       
#> [25] cli_3.6.6           rlang_1.3.0         chk_0.11.0         
#> [28] dbplyr_2.6.0        phiperio_0.6.3      future.apply_1.21.0
#> [31] parallelly_1.48.0   splines_4.6.1       withr_3.0.3        
#> [34] cachem_1.1.0        yaml_2.3.12         vegan_2.7-6        
#> [37] otel_0.2.0          Rtsne_0.17          parallel_4.6.1     
#> [40] tools_4.6.1         ggsignif_0.6.4      ggplot2_4.0.3      
#> [43] globals_0.20.0      buildtools_1.0.0    vctrs_0.7.3        
#> [46] R6_2.6.1            lifecycle_1.0.5     MASS_7.3-66        
#> [49] cluster_2.1.8.3     pkgconfig_2.0.3     RcppParallel_6.2.1 
#> [52] pillar_1.11.1       bslib_0.12.0        gtable_0.3.6       
#> [55] glue_1.8.1          Rcpp_1.1.2          xfun_0.62          
#> [58] tibble_3.3.1        tidyselect_1.2.1    sys_3.4.3          
#> [61] knitr_1.52          farver_2.1.2        nlme_3.1-171       
#> [64] htmltools_0.5.9     parallelDist_0.2.7  maketools_1.3.2    
#> [67] labeling_0.4.3      compiler_4.6.1      S7_0.2.2