Compare Groups: Multi-Group Differential Analysis Engine

Performs differential analysis across N groups. For 2 groups, uses the Wilcoxon-Mann-Whitney test. For 3+ groups, uses Kruskal-Wallis for global significance followed by pairwise Wilcoxon post-hoc tests with Benjamini-Hochberg (BH/FDR) correction.

Usage

av_compare_groups(
  obj,
  group_by,
  compare_by,
  abs_score = TRUE,
  min_n = 5,
  max_items = NULL,
  parallel = FALSE
)

Arguments

obj

An AlphaVarSet object.

group_by

Column name defining the groups to compare (e.g., "tissue"). Supports 2 or more groups.

compare_by

Column name defining the items to compare (e.g., "gene_name"). Must be a categorical column (character or factor).

abs_score

Logical. If TRUE (default), uses absolute scores.

min_n

Minimum sample size per group per item (default 5). Items with fewer observations in any group are excluded with a warning.

max_items

Optional integer. If set and the data contains more than max_items unique items, only the top max_items by within-item variance are tested. Variance is the natural prefilter for "differential signal between groups" because low-spread items cannot drive a significant group difference regardless of group structure. NULL (default) tests all items.

parallel

Logical. Request that the item-level tests run via furrr::future_map(). The request is honoured only when all four conditions hold: furrr is installed, a forking future::plan() is active, it has more than one worker, and at least 750 items pass the min_n guard and therefore actually run a test. Otherwise the call falls back to sequential lapply() and reports the reason once per session. The conditions are measurements, not caution. The mapped closure captures the whole per-item split of the data, so a serialising plan (multisession, cluster) resends the full dataset to every worker and lost at every size tested, aborting above roughly 700 items at the 500 MiB default of future.globals.maxSize. With no plan set, future evaluates sequentially and adds bookkeeping for no benefit. The one backend that pays is future::plan(future::multicore), which shares the split list copy-on-write; it needs a Unix-alike and is unavailable on Windows and, by default, inside RStudio. The gate counts testable items rather than mapped items because items dropped at min_n cost almost nothing to map: splitting the example data by track_name maps 1881 items of which none are testable, and forking those ran 18 times slower than FALSE. Even on eight workers the measured gain at the largest testable set this data reaches is 1.1x to 1.4x, so FALSE (default) is the right setting for typical AlphaGenome exports. Results are identical either way.

Value

The updated AlphaVarSet object with results stored in obj$stats$comparisons[[comparison_name]] as a list containing:

method

Character. e.g., "Wilcoxon-Mann-Whitney (BH)" or "Kruskal-Wallis + Pairwise Wilcoxon (BH)".

global_p

Numeric. For multi-item comparisons with 3+ groups, this is the minimum p-value across all item-specific Kruskal-Wallis tests (not a single omnibus test across items). NA_real_ for 2-group comparisons.

results_tidy

A tibble with columns: item, group_1, group_2, n_1, n_2, median_1, median_2, median_diff, p_val, p_adj, effect_size, effect_type, method, significance_stars. n_1 / n_2 are the sample sizes used in each pairwise test. median_diff is median_1 - median_2 (renamed from diff to avoid masking base::diff()). The effect_size column contains Rank-Biserial Correlation (r) for each pairwise Wilcoxon comparison (bounded from -1 to 1). Positive r indicates group_1 is stochastically larger than group_2 (i.e., group_1 has higher ranks).

params

List. The parameters used for this comparison: group_by, compare_by, groups, abs_score, min_n.

dropped_groups

A named list (keyed by item) of character vectors listing groups that were excluded due to n < min_n while at least 2 valid groups remained. Always empty for 2-group comparisons (an item is fully skipped if either group falls below min_n).

skipped_items

Character vector of items removed entirely because too few groups had n >= min_n to run any test. Mirrors the console warning so it can be inspected post-hoc.

Details

This function implements the multi-group differential analysis engine for comparing variant effects across biological or technical variables.

Statistical Strategy:

• 2 groups: Wilcoxon-Mann-Whitney (non-parametric, robust to outliers). No valid omnibus test; global_p is set to NA_real_.

• 3+ groups: Kruskal-Wallis for omnibus test, followed by pairwise Wilcoxon tests as post-hoc comparisons.

• Correction: Benjamini-Hochberg FDR across all pairwise tests.

• Effect Sizes: Rank-Biserial Correlation (\(r = 1 - 2U / (n_1 \cdot n_2)\)) for each pairwise Wilcoxon test. Positive r indicates group_1 is stochastically larger than group_2. Per-item Epsilon-squared (\(\varepsilon^2 = H / (N - 1)\)) is stored in results_tidy for Kruskal-Wallis items.

• Scores: Uses the active score via get_sc(obj). Switch with av_set_active_score(obj, "raw") before calling.

Examples

example_file <- system.file("extdata", "bcl11a_rs1427407_example.csv", package = "AlphaVaR")
obj <- av_import_csv(example_file)
#> Reading bcl11a_rs1427407_example.csv ...

# The two alleles at chr2:60490908, tested modality by modality.
obj <- av_compare_groups(obj, group_by = "variant_id", compare_by = "output_type")
#> Warning: 1 items skipped due to insufficient sample size (n < 5 ): SPLICE_SITES 
obj$stats$comparisons[["variant_id_by_output_type"]]$results_tidy
#> # A tibble: 10 × 16
#>    item       group_1 group_2   n_1   n_2 median_1 median_2 median_diff    p_val
#>    <chr>      <chr>   <chr>   <int> <int>    <dbl>    <dbl>       <dbl>    <dbl>
#>  1 SPLICE_SI… chr2:6… chr2:6…   160   160    0.474    0.268      0.206  0       
#>  2 ATAC       chr2:6… chr2:6…   160   160    0.892    0.725      0.167  9.70e-13
#>  3 CAGE       chr2:6… chr2:6…   160   160    0.768    0.611      0.157  3.94e- 4
#>  4 CHIP_TF    chr2:6… chr2:6…   164   164    0.539    0.415      0.123  4.72e- 3
#>  5 RNA_SEQ    chr2:6… chr2:6…   160   160    0.313    0.251      0.0628 9.25e- 2
#>  6 CONTACT_M… chr2:6… chr2:6…    28    28    0.555    0.523      0.0315 1.95e- 1
#>  7 DNASE      chr2:6… chr2:6…   160   160    0.563    0.509      0.0545 2.61e- 1
#>  8 SPLICE_JU… chr2:6… chr2:6…   160   160    0.264    0.294     -0.0301 3.06e- 1
#>  9 CHIP_HIST… chr2:6… chr2:6…   160   160    0.661    0.622      0.0398 3.51e- 1
#> 10 PROCAP     chr2:6… chr2:6…    24    24    0.596    0.488      0.108  4.15e- 1
#> # ℹ 7 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> #   epsilon_sq <dbl>, kw_p <dbl>, p_adj <dbl>, significance_stars <chr>

# Five biosample types: Kruskal-Wallis per modality plus pairwise post-hoc.
obj <- av_compare_groups(obj, group_by = "biosample_type", compare_by = "output_type")
#> Warning: Dropped 4 of 2676 rows with NA in 'biosample_type', 'output_type', or 'quantile_score'.
#> Warning: 2 items skipped due to insufficient sample size (n < 5 ): CONTACT_MAPS, SPLICE_SITES 
res <- obj$stats$comparisons[["biosample_type_by_output_type"]]
res$method
#> [1] "Kruskal-Wallis + Pairwise Wilcoxon (BH)"
res$global_p
#> [1] 0.0002449143
head(res$results_tidy)
#> # A tibble: 6 × 16
#>   item    group_1      group_2   n_1   n_2 median_1 median_2 median_diff   p_val
#>   <chr>   <chr>        <chr>   <int> <int>    <dbl>    <dbl>       <dbl>   <dbl>
#> 1 ATAC    cell_line    tissue    198    97    0.852    0.668      0.184  7.54e-5
#> 2 DNASE   primary_cell tissue     87   108    0.692    0.459      0.233  2.21e-4
#> 3 CAGE    primary_cell tissue    187   108    0.614    0.761     -0.147  5.63e-3
#> 4 DNASE   cell_line    primar…   108    87    0.529    0.692     -0.163  1.03e-2
#> 5 ATAC    primary_cell tissue     22    97    0.774    0.668      0.106  1.98e-2
#> 6 RNA_SEQ primary_cell tissue     78   167    0.318    0.261      0.0569 2.14e-2
#> # ℹ 7 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> #   kw_p <dbl>, epsilon_sq <dbl>, p_adj <dbl>, significance_stars <chr>

# max_items prefilters by within-item variance instead of testing all 23 genes.
obj <- av_compare_groups(obj, group_by = "variant_id", compare_by = "gene_name",
                         max_items = 5)
#> Warning: Dropped 1712 of 2676 rows with NA in 'variant_id', 'gene_name', or 'quantile_score'.
#> Prefilter: keeping top 5 of 23 items by within-item variance.
#> Warning: 2 items skipped due to insufficient sample size (n < 5 ): IFITM3P9, NONOP2 
obj$stats$comparisons[["variant_id_by_gene_name"]]$results_tidy
#> # A tibble: 3 × 16
#>   item           group_1 group_2   n_1   n_2 median_1 median_2 median_diff p_val
#>   <chr>          <chr>   <chr>   <int> <int>    <dbl>    <dbl>       <dbl> <dbl>
#> 1 ENSG000002675… chr2:6… chr2:6…     5     9    0.756    0.537      0.219  0.110
#> 2 RNU4-51P       chr2:6… chr2:6…     8     8    0.548    0.395      0.154  0.713
#> 3 PAPOLG         chr2:6… chr2:6…     8     7    0.793    0.834     -0.0408 0.562
#> # ℹ 7 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> #   epsilon_sq <dbl>, kw_p <dbl>, p_adj <dbl>, significance_stars <chr>