Compare Against Reference: Many-to-One Differential Analysis
Source:R/05_inferential.R
av_compare_against_reference.RdTests all other groups against a single reference element using pairwise Wilcoxon tests with Benjamini-Hochberg (BH) correction. This is the statistically appropriate approach for targeted biological hypotheses (e.g., "Is each region type enriched compared to intergenic?").
Usage
av_compare_against_reference(
obj,
group_by,
compare_by,
reference_level,
abs_score = TRUE,
min_n = 5
)Arguments
- obj
An
AlphaVarSetobject.- group_by
Column name defining the groups (e.g., "region_type").
- compare_by
Column name defining the items (e.g., "gene_name"). Must be a categorical column (character or factor).
- reference_level
Character. The reference group label to test against (e.g., "intergenic"). Must exist in the
group_bycolumn.- abs_score
Logical. If TRUE (default), uses absolute scores.
- min_n
Minimum sample size per group per item (default 5).
Value
The updated AlphaVarSet object with results stored in
obj$stats$comparisons[[comparison_name]] (named
"<group_by>_vs_<reference_level>_by_<compare_by>"). Structure
matches av_compare_groups() output (method, global_p, results_tidy,
params, dropped_groups).
Details
Why Many-to-One instead of All-vs-All? When the biological question targets enrichment relative to a neutral baseline (e.g., intergenic), testing all \(k(k-1)/2\) pairs inflates the multiple testing burden unnecessarily. Many-to-One tests only \(k-1\) comparisons per item, maximizing statistical power under BH correction.
Effect Size: Rank-Biserial Correlation (r). Positive r indicates the non-reference group is stochastically larger than the reference.
Examples
example_file <- system.file("extdata", "bcl11a_rs1427407_example.csv", package = "AlphaVaR")
obj <- av_import_csv(example_file)
#> Reading bcl11a_rs1427407_example.csv ...
# Cell lines as the baseline: k - 1 tests per modality instead of all pairs.
obj <- av_compare_against_reference(
obj,
group_by = "biosample_type",
compare_by = "output_type",
reference_level = "cell_line"
)
#> Warning: Dropped 4 of 2676 rows with NA in 'biosample_type', 'output_type', or 'quantile_score'.
#> Warning: 1 items skipped (reference group n < 5 ): SPLICE_SITES
res <- obj$stats$comparisons[["biosample_type_vs_cell_line_by_output_type"]]
res$method
#> [1] "Many-to-One Wilcoxon vs cell_line (BH)"
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 tissue cell_l… 97 198 0.668 0.852 -0.184 7.54e-5
#> 2 DNASE primar… cell_l… 87 108 0.692 0.529 0.163 1.03e-2
#> 3 CHIP_HISTONE in_vit… cell_l… 34 115 0.533 0.655 -0.123 6.53e-2
#> 4 DNASE tissue cell_l… 108 108 0.459 0.529 -0.0698 6.33e-2
#> 5 PROCAP primar… cell_l… 8 40 0.397 0.596 -0.199 1.28e-1
#> 6 CHIP_TF tissue cell_l… 21 286 0.569 0.468 0.101 1.13e-1
#> # ℹ 7 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> # epsilon_sq <dbl>, kw_p <dbl>, p_adj <dbl>, significance_stars <chr>