Skip to content

Commit 977059c

Browse files
Merge branch 'dev'
2 parents 56efd93 + 1f27924 commit 977059c

2 files changed

Lines changed: 74 additions & 188 deletions

File tree

Lines changed: 57 additions & 178 deletions
Original file line numberDiff line numberDiff line change
@@ -1,106 +1,60 @@
11
#!/usr/bin/env python3
22
"""
3-
Run project-specific two-group tests for mean IF cluster data.
3+
Run project-specific two-group tests for mean intensities from cluster data.
44
5-
This is a simpler companion to the Psi/Iso 2x2 ANOVA script for cases where
6-
only two conditions are present, e.g.:
5+
Prereqs:
6+
- `cstats_mean_IF` to generate mean intensity CSVs for each sample.
77
8-
AwS = Awake Saline
9-
AnS = Anesthetized Saline
10-
11-
Input convention:
12-
- each input CSV is one sample
8+
Inputs:
9+
- CSVs in the working directory (one per sample) or specified with -i
1310
- the condition is encoded as the first underscore-delimited field in the filename
14-
e.g., AwS_sample01_cluster_mean_IF.csv -> condition AwS
15-
- each CSV contains one row per cluster
11+
e.g., AwS_sample01_cluster_mean_intensity.csv -> condition AwS
12+
- each CSV contains one row per cluster (columns: cluster_ID, mean_intensity, etc.)
1613
1714
For each cluster, outputs:
1815
- group means, SD, SEM, and n for the two conditions
1916
- difference: second condition minus first condition
2017
default: AnS - AwS
21-
- Welch t-test by default, or equal-variance t-test with --equal_var
22-
- Holm-Sidak and Benjamini-Hochberg FDR correction across clusters
18+
- Student's t-test (equal-variance)
2319
- higher-mean group
2420
- Cohen's d and Hedges' g, using condition order:
2521
positive = second condition > first condition
2622
27-
Example:
28-
mean_IF_2group_ttest.py \
29-
-i 'cluster_mean_IF_AnS_v_AwS_vox_p_tstat1_q0.05_rev_cluster_index_LH/*.csv' \
30-
-o mean_IF_AnS_vs_AwS_ttest_results.csv
31-
32-
Equivalent explicit conditions:
33-
mean_IF_2group_ttest.py \
34-
-i '*.csv' \
35-
--conditions AwS AnS \
36-
-o mean_IF_AnS_minus_AwS_ttest_results.csv
23+
Usage:
24+
------
25+
./mean_IF_2group_ttest.py -c AwS AnS [OPTIONS]
3726
"""
3827

39-
import argparse
40-
from pathlib import Path
41-
import glob as globlib
42-
4328
import numpy as np
4429
import pandas as pd
30+
from pathlib import Path
31+
from rich import print
32+
from rich.traceback import install
4533
from scipy.stats import ttest_ind
46-
from statsmodels.stats.multitest import multipletests
47-
4834

49-
DEFAULT_CONDITIONS = ["AwS", "AnS"]
35+
from unravel.core.help_formatter import RichArgumentParser, SuppressMetavar, SM
36+
from unravel.core.config import Configuration
37+
from unravel.core.utils import log_command, match_files, verbose_end_msg, verbose_start_msg
5038

5139

5240
def parse_args():
53-
p = argparse.ArgumentParser(
54-
description="Project-specific two-group t-tests for mean_IF cluster CSVs."
55-
)
56-
p.add_argument(
57-
"-i", "--input_dir", dest="inputs", nargs="*", default=None,
58-
help=(
59-
"Input directory, CSV file, or quoted glob pattern. May be repeated as "
60-
"multiple values after -i. Default: current directory."
61-
),
62-
)
63-
p.add_argument(
64-
"-g", "--glob", default="*.csv",
65-
help="Input CSV glob used when an input is a directory. Default: *.csv",
66-
)
67-
p.add_argument(
68-
"--conditions", nargs=2, default=DEFAULT_CONDITIONS, metavar=("REF", "TEST"),
69-
help=(
70-
"Two condition prefixes to compare. Difference/effect size is TEST - REF. "
71-
"Default: AwS AnS, so diff = AnS - AwS."
72-
),
73-
)
74-
p.add_argument(
75-
"-vcol", "--value_col", default="mean_IF_intensity",
76-
help="Dependent-variable column. Default: mean_IF_intensity",
77-
)
78-
p.add_argument(
79-
"-ccol", "--cluster_col", default="cluster_ID",
80-
help="Cluster ID column. Default: cluster_ID",
81-
)
82-
p.add_argument(
83-
"-scol", "--sample_col", default="sample",
84-
help="Sample column. If absent, filename stem is used. Default: sample",
85-
)
86-
p.add_argument(
87-
"-o", "--out", default="mean_IF_2group_ttest_results.csv",
88-
help="Output CSV path. Default: mean_IF_2group_ttest_results.csv",
89-
)
90-
p.add_argument(
91-
"--raw_out", default="mean_IF_2group_ttest_raw_data.csv",
92-
help="Raw long-format output CSV path. Default: mean_IF_2group_ttest_raw_data.csv",
93-
)
94-
p.add_argument(
95-
"--equal_var", action="store_true", default=False,
96-
help="Use equal-variance t-tests. Default: Welch t-tests.",
97-
)
98-
p.add_argument(
99-
"--min_n_per_group", type=int, default=2,
100-
help="Minimum n per condition needed for the t-test. Default: 2",
101-
)
102-
return p.parse_args()
41+
parser = RichArgumentParser(formatter_class=SuppressMetavar, add_help=False, docstring=__doc__)
42+
43+
reqs = parser.add_argument_group('Required arguments')
44+
reqs.add_argument('-c', '--conditions', help='Conditions. Each condition must match first word of CSV filenames (underscore-separated).', required=True, nargs='*', action=SM)
45+
46+
opts = parser.add_argument_group('Optional arguments')
47+
opts.add_argument('-i', '--input', help="Path(s) or glob pattern(s) for input CSV files. Default: '*.csv'", default='*.csv', nargs='*', action=SM)
48+
opts.add_argument("-o", "--out", default="out/mean_intensity_2group_ttest_results.csv", help="Output CSV path. Default: out/mean_intensity_2group_ttest_results.csv")
49+
opts.add_argument("-r", "--raw_out", default="out/mean_intensity_2group_ttest_raw_data.csv", help="Raw long-format output CSV path. Default: out/mean_intensity_2group_ttest_raw_data.csv")
50+
opts.add_argument("-vc", "--value_col", default="mean_intensity", help="Dependent-variable column. Default: mean_intensity")
51+
opts.add_argument("-cc", "--cluster_col", default="cluster_ID", help="Cluster ID column. Default: cluster_ID")
52+
opts.add_argument("-sc", "--sample_col", default="sample", help="Sample column. If absent, filename stem is used. Default: sample")
10353

54+
general = parser.add_argument_group('General arguments')
55+
general.add_argument('-v', '--verbose', help='Verbose output.', action='store_true')
56+
57+
return parser.parse_args()
10458

10559
def p_to_sig(p):
10660
if pd.isna(p):
@@ -128,58 +82,7 @@ def safe_float(x):
12882
def infer_condition(filename):
12983
return Path(filename).name.split("_")[0]
13084

131-
132-
def expand_input_csvs(inputs, glob_pattern):
133-
"""Expand directories, individual CSV files, and quoted glob patterns."""
134-
if not inputs:
135-
inputs = ["."]
136-
137-
csvs = []
138-
missing_inputs = []
139-
140-
for item in inputs:
141-
item = str(item)
142-
143-
# Allows commands like:
144-
# -i 'cluster_mean_IF_dir/*.csv'
145-
if globlib.has_magic(item):
146-
matches = [Path(m) for m in globlib.glob(item, recursive=True)]
147-
csvs.extend([
148-
m for m in matches
149-
if m.is_file() and m.name.lower().endswith(".csv")
150-
])
151-
continue
152-
153-
path = Path(item).expanduser()
154-
155-
if path.is_dir():
156-
csvs.extend([
157-
m for m in path.glob(glob_pattern)
158-
if m.is_file() and m.name.lower().endswith(".csv")
159-
])
160-
elif path.is_file():
161-
if path.name.lower().endswith(".csv"):
162-
csvs.append(path)
163-
else:
164-
missing_inputs.append(item)
165-
166-
# De-duplicate while keeping deterministic order.
167-
csvs = sorted({str(p): p for p in csvs}.values(), key=lambda p: str(p))
168-
169-
if not csvs:
170-
msg = (
171-
"No input CSVs found. "
172-
f"inputs={inputs!r}; directory glob={glob_pattern!r}"
173-
)
174-
if missing_inputs:
175-
msg += f"; missing inputs={missing_inputs!r}"
176-
raise FileNotFoundError(msg)
177-
178-
return csvs
179-
180-
181-
def load_mean_if_data(inputs, glob_pattern, conditions, value_col, cluster_col, sample_col):
182-
csvs = expand_input_csvs(inputs, glob_pattern)
85+
def load_mean_intensity_data(csvs, conditions, value_col, cluster_col, sample_col):
18386
conditions = [str(c) for c in conditions]
18487
condition_set = set(conditions)
18588

@@ -259,7 +162,7 @@ def higher_by_mean(mean_ref, mean_test, ref_cond, test_cond):
259162
return test_cond if mean_test > mean_ref else ref_cond
260163

261164

262-
def cohen_d_and_hedges_g(vals_ref, vals_test):
165+
def hedges_g(vals_ref, vals_test):
263166
"""
264167
Standardized mean difference using pooled SD.
265168
@@ -289,18 +192,17 @@ def cohen_d_and_hedges_g(vals_ref, vals_test):
289192
# Small-sample correction.
290193
correction = 1 - (3 / (4 * (n_ref + n_test) - 9)) if (n_ref + n_test) > 2 else np.nan
291194
g = d * correction if not pd.isna(correction) else np.nan
292-
return safe_float(d), safe_float(g)
195+
return safe_float(g)
293196

294197

295-
def run_ttest(vals_ref, vals_test, equal_var, min_n_per_group):
198+
def run_ttest(vals_ref, vals_test, min_n_per_group=2):
296199
if len(vals_ref) < min_n_per_group or len(vals_test) < min_n_per_group:
297200
return np.nan, np.nan
298201

299202
try:
300203
res = ttest_ind(
301204
vals_test,
302205
vals_ref,
303-
equal_var=equal_var,
304206
alternative="two-sided",
305207
nan_policy="omit",
306208
)
@@ -309,38 +211,7 @@ def run_ttest(vals_ref, vals_test, equal_var, min_n_per_group):
309211
return np.nan, np.nan
310212

311213

312-
def add_multiple_testing_columns(results):
313-
"""Add Holm-Sidak and BH-FDR adjusted p-values across clusters."""
314-
results = results.copy()
315-
316-
results["p_holm_sidak"] = np.nan
317-
results["sig_holm_sidak"] = ""
318-
results["reject_holm_sidak_0.05"] = False
319-
320-
results["p_fdr_bh"] = np.nan
321-
results["sig_fdr_bh"] = ""
322-
results["reject_fdr_bh_0.05"] = False
323-
324-
pvals = results["p"].to_numpy(dtype=float)
325-
valid = ~pd.isna(pvals)
326-
327-
if valid.sum() > 0:
328-
reject_hs, p_hs, _, _ = multipletests(pvals[valid], alpha=0.05, method="holm-sidak")
329-
reject_bh, p_bh, _, _ = multipletests(pvals[valid], alpha=0.05, method="fdr_bh")
330-
331-
valid_idx = results.index[valid]
332-
results.loc[valid_idx, "p_holm_sidak"] = p_hs
333-
results.loc[valid_idx, "reject_holm_sidak_0.05"] = reject_hs
334-
results.loc[valid_idx, "sig_holm_sidak"] = [p_to_sig(p) for p in p_hs]
335-
336-
results.loc[valid_idx, "p_fdr_bh"] = p_bh
337-
results.loc[valid_idx, "reject_fdr_bh_0.05"] = reject_bh
338-
results.loc[valid_idx, "sig_fdr_bh"] = [p_to_sig(p) for p in p_bh]
339-
340-
return results
341-
342-
343-
def build_results(data, conditions, value_col, cluster_col, equal_var, min_n_per_group):
214+
def build_results(data, conditions, value_col, cluster_col, min_n_per_group=2):
344215
ref_cond, test_cond = [str(c) for c in conditions]
345216
diff_col = f"{test_cond}_minus_{ref_cond}"
346217

@@ -364,11 +235,10 @@ def build_results(data, conditions, value_col, cluster_col, equal_var, min_n_per
364235
t_stat, p = run_ttest(
365236
vals_ref=vals_ref,
366237
vals_test=vals_test,
367-
equal_var=equal_var,
368238
min_n_per_group=min_n_per_group,
369239
)
370240

371-
cohen_d, hedges_g = cohen_d_and_hedges_g(vals_ref, vals_test)
241+
hedges_g_val = hedges_g(vals_ref, vals_test)
372242

373243
rows.append({
374244
"cluster_ID": cluster_id,
@@ -385,9 +255,8 @@ def build_results(data, conditions, value_col, cluster_col, equal_var, min_n_per
385255
"t_stat": t_stat,
386256
"p": p,
387257
"sig": p_to_sig(p),
388-
"cohen_d": cohen_d,
389-
"hedges_g": hedges_g,
390-
"test": "Student_ttest" if equal_var else "Welch_ttest",
258+
"hedges_g": hedges_g_val,
259+
"test": "Student_ttest",
391260
"effect_direction": (
392261
f"{test_cond}>{ref_cond}" if not pd.isna(diff) and diff > 0
393262
else f"{test_cond}<{ref_cond}" if not pd.isna(diff) and diff < 0
@@ -397,18 +266,28 @@ def build_results(data, conditions, value_col, cluster_col, equal_var, min_n_per
397266
})
398267

399268
results = pd.DataFrame(rows)
400-
results = add_multiple_testing_columns(results)
401269

402270
sort_cols = ["p", "cluster_ID"]
403271
return results.sort_values(sort_cols, na_position="last")
404272

405273

274+
@log_command
406275
def main():
276+
install()
407277
args = parse_args()
278+
Configuration.verbose = args.verbose
279+
verbose_start_msg()
280+
281+
csvs = match_files(args.input)
408282

409-
data = load_mean_if_data(
410-
inputs=args.inputs,
411-
glob_pattern=args.glob,
283+
# Filter CSVs to only those that match the requested conditions based on filename prefix
284+
csvs = [
285+
f for f in csvs
286+
if any(f.name.startswith(f'{group}_') for group in args.conditions)
287+
]
288+
289+
data = load_mean_intensity_data(
290+
csvs=csvs,
412291
conditions=args.conditions,
413292
value_col=args.value_col,
414293
cluster_col=args.cluster_col,
@@ -424,8 +303,6 @@ def main():
424303
conditions=args.conditions,
425304
value_col=args.value_col,
426305
cluster_col=args.cluster_col,
427-
equal_var=args.equal_var,
428-
min_n_per_group=args.min_n_per_group,
429306
)
430307

431308
out = Path(args.out)
@@ -438,6 +315,8 @@ def main():
438315
print(f"Comparison: {test_cond} - {ref_cond}")
439316
print(results.head(20).to_string(index=False))
440317

318+
verbose_end_msg()
319+
441320

442321
if __name__ == "__main__":
443322
main()

0 commit comments

Comments
 (0)