Skip to contents

Sweeps every element of group_by, testing each element's score distribution against the pooled rest (the all-element generalization of av_compare_target). Reports a tidy table ranked by shift, with BH correction across elements and, for >= 3 elements, a Kruskal-Wallis omnibus. One-vs-rest is a marginal screen ("which element stands out vs the baseline?"); for a specific pairwise contrast use av_compare_groups.

Usage

av_calc_shift(
  obj,
  group_by = "output_type",
  abs_score = TRUE,
  omnibus = TRUE,
  min_n = 5L,
  balance_by = NULL,
  per = "variant_id"
)

Arguments

obj

An AlphaVarSet object.

group_by

Character. Column whose elements are swept (default "output_type").

abs_score

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

omnibus

Logical. If TRUE (default), adds a Kruskal-Wallis omnibus across all elements when there are >= 3 of them.

min_n

Minimum n per arm (default 5). Elements below it are dropped.

balance_by

Character or NULL. When set, collapses each arm to one balanced value per unit per (two-step weighted mean over balance_by) before each element's test. Must differ from per and group_by.

per

Character. Unit (replication) unit for balance_by; default "variant_id". When a per unit straddles several elements of group_by (a per-track column), the sweep uses paired Wilcoxon signed-rank tests plus a Friedman omnibus blocked by per; otherwise Mann-Whitney plus a balanced Kruskal-Wallis omnibus. Ignored when balance_by is NULL.

Value

The updated AlphaVarSet with results at obj$stats$comparisons[[paste0(group_by, "_vs_rest")]]: a list with method, global_p (Kruskal-Wallis or Friedman p, or NA_real_ when no omnibus), global_effect / global_effect_type (epsilon-squared or Kendall's W), results_tidy (one row per kept element, schema as in av_compare_target), params, dropped_groups.

See also

av_compare_target for a single element.

Examples

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

obj <- av_calc_shift(obj, group_by = "output_type")
#> Warning: 1 element(s) dropped (n < 5 in one arm): SPLICE_SITES
res <- obj$stats$comparisons[["output_type_vs_rest"]]
res$method
#> [1] "Kruskal-Wallis Omnibus + Wilcoxon One-vs-Rest (BH)"
res$global_p
#> [1] 4.609982e-130
res$results_tidy
#> # A tibble: 10 × 12
#>    group   n_group n_rest median_group median_rest median_diff    p_val    p_adj
#>    <chr>     <int>  <int>        <dbl>       <dbl>       <dbl>    <dbl>    <dbl>
#>  1 ATAC        320   2356        0.814       0.454     0.360   0        0       
#>  2 CAGE        320   2356        0.696       0.465     0.231   0        0       
#>  3 SPLICE…     320   2356        0.277       0.532    -0.254   7.46e-43 2.49e-42
#>  4 RNA_SEQ     320   2356        0.270       0.517    -0.247   2.29e-26 5.72e-26
#>  5 SPLICE…     320   2356        0.336       0.517    -0.181   2.71e-14 5.43e-14
#>  6 CHIP_H…     320   2356        0.640       0.465     0.175   5.02e-10 8.37e-10
#>  7 DNASE       320   2356        0.529       0.483     0.0461  3.05e- 1 4.36e- 1
#>  8 CONTAC…      56   2620        0.535       0.488     0.0461  5.30e- 1 6.63e- 1
#>  9 CHIP_TF     328   2348        0.491       0.494    -0.00288 7.21e- 1 8.01e- 1
#> 10 PROCAP       48   2628        0.565       0.494     0.0705  8.50e- 1 8.50e- 1
#> # ℹ 4 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> #   significance_stars <chr>

# Tracks without a biosample label cannot be blocked by biosample.
obj_bs <- av_subset(obj, !is.na(biosample_name))
obj_bs <- av_calc_shift(obj_bs, group_by = "output_type",
                        balance_by = "variant_id", per = "biosample_name")
#> Warning: Friedman omnibus: dropped 620 of 620 'biosample_name' units with incomplete element coverage.
obj_bs$stats$comparisons[["output_type_vs_rest"]]$results_tidy
#> # A tibble: 10 × 12
#>    group   n_group n_rest median_group median_rest median_diff    p_val    p_adj
#>    <chr>     <int>  <int>        <dbl>       <dbl>       <dbl>    <dbl>    <dbl>
#>  1 SPLICE…     177    177        0.294       0.467    -0.178   1.07e-17 1.07e-16
#>  2 CHIP_H…     126    126        0.623       0.415     0.192   1.67e-13 8.37e-13
#>  3 ATAC         69     69        0.686       0.416     0.289   2.05e-11 6.83e-11
#>  4 RNA_SEQ     172    172        0.298       0.469    -0.146   8.05e- 6 2.01e- 5
#>  5 CAGE         78     78        0.607       0.425     0.187   2.57e- 5 5.14e- 5
#>  6 SPLICE…     202    202        0.365       0.441    -0.0787  1.67e- 4 2.79e- 4
#>  7 CONTAC…       8      8        0.491       0.412     0.129   5.87e- 2 8.39e- 2
#>  8 DNASE       149    149        0.535       0.475     0.0296  2.35e- 1 2.93e- 1
#>  9 PROCAP        6      6        0.511       0.482    -0.00843 8.34e- 1 9.27e- 1
#> 10 CHIP_TF      45     45        0.506       0.497     0.0223  9.64e- 1 9.64e- 1
#> # ℹ 4 more variables: effect_size <dbl>, effect_type <chr>, method <chr>,
#> #   significance_stars <chr>