Skip to content

Refactor group analysis CLI: replace contrast_column/values with formula-based model + pairwise contrasts #147

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Draft
Copilot wants to merge 4 commits into
base: main
Choose a base branch
from
Draft

Refactor group analysis CLI: replace contrast_column/values with formula-based model + pairwise contrasts #147

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
eb47be8
Initial plan
Copilot Apr 5, 2026
3e4a0f9
Refactor group analysis CLI: replace contrast_column/values with mode…
Copilot Apr 5, 2026
28f6bca
Address code review: improve strata naming, add file existence check,…
Copilot Apr 5, 2026
27d929a
Address code review: move imports to top of Snakefile, rename private…
Copilot Apr 5, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
15 changes: 10 additions & 5 deletions spimquant/config/snakebids.yml
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -250,15 +250,20 @@ parse_args:
action: store
nargs: '+'

--contrast_column:
help: "Column name in participants.tsv to use for defining group contrasts (e.g., 'treatment', 'genotype'). Required for group-level statistical analysis."
--model:
help: "Statistical model formula in patsy/statsmodels format, with 'metric' as a placeholder for the response variable (e.g. 'metric ~ C(treatment) * C(genotype) + age'). Required for group-level analysis."
default: null
type: str

--contrast_values:
help: "Two group values for contrast comparison (e.g., 'control' 'drug'). Used with --contrast_column for statistical testing. Provide exactly 2 values."
--pairwise:
help: "Factor(s) in participants.tsv for which all pairwise comparisons should be computed. Can be specified multiple times (e.g. '--pairwise treatment --pairwise genotype')."
default: null
nargs: 2
nargs: '+'

--within:
help: "Factors defining strata within which pairwise contrasts are computed (e.g. '--within genotype sex'). All combinations of these factor levels will be used as separate strata."
default: null
nargs: '+'


#--- workflow specific configuration --
Expand Down
158 changes: 80 additions & 78 deletions spimquant/workflow/Snakefile
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,7 @@
import os
from itertools import combinations as _combinations, product as _product

import pandas as _pd
from zarrnii import ZarrNii
from snakemake.utils import format
from snakebids import bids, generate_inputs, get_wildcard_constraints, set_bids_spec
Expand Down Expand Up @@ -141,19 +143,80 @@ for seg in config["crop_atlas_segs"]:
)
else:
crop_atlas_segs.append(seg)
# Validate that contrast arguments are provided when using group analysis level
# Validate that model/pairwise arguments are provided when using group analysis level
# and generate pairwise contrast labels from participants.tsv at planning time
pairwise_contrast_labels = []
pairwise_contrast_info = {} # label -> dict with factor, levelA, levelB, strata

if config.get("analysis_level") == "group":
if config.get("contrast_column") is None or config.get("contrast_values") is None:
if not config.get("model"):
raise ValueError(
"When using group analysis level, --model must be specified "
"(e.g. 'metric ~ C(treatment) + age')."
)
if not config.get("pairwise"):
raise ValueError(
"When using group analysis level, both --contrast_column and "
"--contrast_values must be specified for filtering group data by contrasts."
"When using group analysis level, --pairwise must be specified "
"(e.g. '--pairwise treatment')."
)

_participants_tsv = os.path.join(config["bids_dir"], "participants.tsv")
if not os.path.exists(_participants_tsv):
raise ValueError(
f"participants.tsv not found at '{_participants_tsv}'. "
"This file is required for group-level analysis."
)
_participants_df = _pd.read_csv(_participants_tsv, sep="\t")

_pairwise_factors = config["pairwise"]
_within_factors = config.get("within") or []

for _factor in _pairwise_factors:
if _factor not in _participants_df.columns:
raise ValueError(
f"Pairwise factor '{_factor}' not found in participants.tsv. "
f"Available columns: {list(_participants_df.columns)}"
)
_levels = sorted(str(v) for v in _participants_df[_factor].dropna().unique())

if _within_factors:
_within_level_lists = []
for _wf in _within_factors:
if _wf not in _participants_df.columns:
raise ValueError(
f"Within factor '{_wf}' not found in participants.tsv. "
f"Available columns: {list(_participants_df.columns)}"
)
_within_level_lists.append(
sorted(str(v) for v in _participants_df[_wf].dropna().unique())
)
for _lA, _lB in _combinations(_levels, 2):
for _stratum in _product(*_within_level_lists):
_strata_dict = dict(zip(_within_factors, _stratum))
_strata_str = "+".join(f"{f}-{v}" for f, v in _strata_dict.items())
_label = f"{_factor}+{_lA}vs{_lB}+{_strata_str}"
pairwise_contrast_labels.append(_label)
pairwise_contrast_info[_label] = {
"factor": _factor,
"levelA": _lA,
"levelB": _lB,
"strata": _strata_dict,
}
else:
for _lA, _lB in _combinations(_levels, 2):
_label = f"{_factor}+{_lA}vs{_lB}"
pairwise_contrast_labels.append(_label)
pairwise_contrast_info[_label] = {
"factor": _factor,
"levelA": _lA,
"levelB": _lB,
"strata": {},
}


wildcard_constraints:
stain="[a-zA-Z0-9]+",
contrast_column="[a-zA-Z0-9_]+",
contrast_value="[a-zA-Z0-9_]+",
pairwise_contrast="[a-zA-Z0-9+_-]+",


rule all_templatereg_deform_zooms:
Expand Down Expand Up @@ -511,27 +574,32 @@ rule all_imaris_crops:
rule all_group_stats:
"""Target rule for group-level statistical analysis."""
input:
# Per-contrast groupstats TSV and PNG
expand(
bids(
root=root,
datatype="group",
seg="{seg}",
from_="{template}",
desc="{desc}",
contrast="{pairwise_contrast}",
suffix="groupstats.{ext}",
),
seg=atlas_segs,
desc=config["seg_method"],
template=config["template"],
pairwise_contrast=pairwise_contrast_labels,
ext=["png", "tsv"],
),
# Per-contrast NIfTI stat maps for each seg metric
expand(
bids(
root=root,
datatype="group",
seg="{seg}",
space="{template}",
desc="{desc}",
contrast="{pairwise_contrast}",
metric="{stain}+{metric}",
suffix="{stat}.nii.gz",
),
Expand All @@ -541,7 +609,9 @@ rule all_group_stats:
stain=stains_for_seg,
metric=config["seg_metrics"],
stat=config["stats_maps"],
pairwise_contrast=pairwise_contrast_labels,
),
# All-subjects count maps (not contrast-specific)
expand(
bids(
root=root,
Expand All @@ -556,55 +626,19 @@ rule all_group_stats:
level=range(4),
stain=stains_for_seg,
),
# Add contrast-filtered outputs
expand(
bids(
root=root,
datatype="group",
level="{level}",
space="{template}",
desc="{desc}",
contrast="{contrast_column}+{contrast_value}",
suffix="{stain}+count.nii.gz",
),
desc=config["seg_method"],
template=config["template"],
level=range(4),
stain=stains_for_seg,
contrast_column=config["contrast_column"],
contrast_value=config["contrast_values"],
),
# Add group-averaged segstats maps by contrast
expand(
bids(
root=root,
datatype="group",
seg="{seg}",
space="{template}",
desc="{desc}",
contrast="{contrast_column}+{contrast_value}",
metric="{stain}+{metric}",
suffix="groupavg.nii.gz",
),
seg=atlas_segs,
desc=config["seg_method"],
template=config["template"],
stain=stains_for_seg,
metric=config["seg_metrics"],
contrast_column=config["contrast_column"],
contrast_value=config["contrast_values"],
),


rule all_group_stats_coloc:
input:
# Per-contrast NIfTI stat maps for coloc metrics
expand(
bids(
root=root,
datatype="group",
seg="{seg}",
space="{template}",
desc="{desc}",
contrast="{pairwise_contrast}",
metric="coloc+{metric}",
suffix="{stat}.nii.gz",
),
Expand All @@ -613,7 +647,9 @@ rule all_group_stats_coloc:
template=config["template"],
metric=config["coloc_seg_metrics"],
stat=config["stats_maps"],
pairwise_contrast=pairwise_contrast_labels,
),
# All-subjects colocalization count maps (not contrast-specific)
expand(
bids(
root=root,
Expand All @@ -627,40 +663,6 @@ rule all_group_stats_coloc:
template=config["template"],
level=range(4),
),
expand(
bids(
root=root,
datatype="group",
space="{template}",
level="{level}",
desc="{desc}",
contrast="{contrast_column}+{contrast_value}",
suffix="coloccount.nii.gz",
),
desc=config["seg_method"],
template=config["template"],
level=range(4),
contrast_column=config["contrast_column"],
contrast_value=config["contrast_values"],
),
expand(
bids(
root=root,
datatype="group",
seg="{seg}",
space="{template}",
desc="{desc}",
contrast="{contrast_column}+{contrast_value}",
metric="coloc+{metric}",
suffix="groupavg.nii.gz",
),
seg=atlas_segs,
desc=config["seg_method"],
template=config["template"],
metric=config["coloc_seg_metrics"],
contrast_column=config["contrast_column"],
contrast_value=config["contrast_values"],
),


rule all_qc:
Expand Down
Loading