Skip to content

How to compare specific pairwise contrasts for lfc and heatmap (without reference group)? 6 treatment groups, to compare 3 pairs #320

@suminayeoh

Description

@suminayeoh

Hello ANCOM-BC experts - I’d appreciate advice on how to parameterize ANCOM-BC2 so pairwise contrasts for all my requested comparisons show up reproducibly (I’m seeing single-index columns referencing one baseline and missing the two-index pair columns I expect).

Short experimental design

Subset: ps_Chap3_DA_ITS_AT = samples where (Arrival_time == "CA" & Month == "12M") | (Arrival_time == "LA" & Month == "24M") and Treatment %in% c("K","M","KM").

N: 30 samples (6 Treat_AT groups × 5 each).

I am trying to study within-treatment arrival-time comparisons (eg. K treatment CA concurrent-arrival vs K treatment late-arrival). Intially I tried to run Treatment * Arrival_time + Block but model failed. So I combined Treatment & Arrival into a variable and ran Treat_AT + Block instead:
Treat_AT = paste(Treatment, Arrival_time, sep = "_") with enforced levels: K_CA, K_LA, KM_CA, KM_LA, M_CA, M_LA.
Block is Block 1 to 5 (was supposed to be covariate as Block were found to be significant in beta diversity analysis)

Exact ANCOM-BC2 call / parameters (what I used)

res <- ancombc2(
data = ps_Chap3_DA_ITS_AT,
tax_level = <NULL or "Phylum"/"Family"/"Genus">,
fix_formula = "Treat_AT + Block",
rand_formula = NULL,
group = "Treat_AT",
p_adj_method = "BH",
prv_cut = 0.10,
lib_cut = 1000,
s0_perc = 0.05,
pseudo_sens = TRUE,
struc_zero = TRUE,
neg_lb = TRUE,
dunnet = FALSE,
alpha = 0.05,
n_cl = 1,
iter_control = list(tol = 1e-2, max_iter = 20, verbose = TRUE),
em_control = list(tol = 1e-5, max_iter = 100),
lme_control = lme4::lmerControl(),
global = TRUE,
pairwise = TRUE
)

Contrasts I specifically want (within-treatment arrival-time comparisons)
K_CA vs K_LA
M_CA vs M_LA
KM_CA vs KM_LA

(Under my enforced ordering these map to Treat_AT1 vs Treat_AT2, Treat_AT5 vs Treat_AT6, Treat_AT3 vs Treat_AT4.)

Problem / question (brief)
res$res_pair shows lfc_Treat_AT1..lfc_Treat_AT5 and pairwise columns like lfc_Treat_AT2_Treat_AT1, but no Treat_AT6 token (so the M_CA vs M_LA pairwise column such as q_Treat_AT6_Treat_AT5 is missing). I did not set dunnet = TRUE or an explicit reference manually; I forced the factor levels in phyloseq before running.

Questions

Is it expected ANCOM-BC2 parameterizes with a single-reference index even when pairwise = TRUE?
Would releveling Treat_AT (so a different reference) force explicit two-index pairwise columns for all contrasts? Example I might try:
sample_data(ps_Chap3_DA_ITS_AT)$Treat_AT <- relevel(sample_data(ps_Chap3_DA_ITS_AT)$Treat_AT, ref="K_CA")

re-run
Alternatives: should I set group = NULL and compute contrasts externally, or change group usage? Which guarantees explicit two-index pairwise columns?

Any quick snippet to detect which reference/parameterization ANCOM used so I can re-run only affected models?

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions