Skills › Scientific Agent Skills › flowkit › analyze_gates.py
analyze_gates.py
204 lines · skills/flowkit/scripts/analyze_gates.py at 5f2bc03
1#!/usr/bin/env python3
2"""Apply a supplied GatingML strategy or FlowJo workspace to explicit FCS files."""
4from __future__ import annotations
6import argparse
7import hashlib
8import json
9from importlib.metadata import version
10from pathlib import Path
13def file_record(path: Path) -> dict:
14 with path.open("rb") as handle:
15 digest = hashlib.file_digest(handle, "sha256").hexdigest()
16 return {"path": str(path.resolve()), "sha256": digest}
19def add_denominators(report, event_counts: dict[str, int]):
20 """Attach population identities and denominators, including quadrant children."""
21 import pandas as pd
23 report = report.copy()
24 population_paths = []
25 counts = {}
26 for row in report.itertuples(index=False):
27 # FlowKit reports a quadrant at its owning QuadrantGate's path, with
28 # that gate stored separately in quadrant_parent. Its children use
29 # the full strategy path, so include the owner when building the LUT.
30 path = tuple(row.gate_path)
31 if pd.notna(row.quadrant_parent):
32 path += (row.quadrant_parent,)
33 path += (row.gate_name,)
34 key = (row.sample_id, path)
35 if key in counts:
36 raise ValueError(f"Duplicate population in analysis report: {key!r}")
37 counts[key] = int(row.count)
38 population_paths.append(path)
40 parent_counts = []
41 for row in report.itertuples(index=False):
42 path = tuple(row.gate_path)
43 if path == ("root",):
44 parent_count = event_counts[row.sample_id]
45 else:
46 key = (row.sample_id, path)
47 if key not in counts:
48 raise ValueError(f"Missing parent population in analysis report: {key!r}")
49 parent_count = counts[key]
50 if not 0 <= row.count <= parent_count <= event_counts[row.sample_id]:
51 raise ValueError(f"Inconsistent gate/parent/sample counts for {row.gate_name!r}")
52 parent_counts.append(parent_count)
54 report["population_path"] = population_paths
55 report["sample_event_count"] = report["sample_id"].map(event_counts)
56 report["parent_event_count"] = parent_counts
57 report["relative_percent_defined"] = report["parent_event_count"] > 0
58 # FlowKit 1.3.2 returns zero for ordinary gates below an empty parent and
59 # NaN for quadrants. Neither is a defined percentage; CSV writes blank.
60 report.loc[~report["relative_percent_defined"], "relative_percent"] = float("nan")
61 return report
64def analyze(
65 *,
66 fcs_paths: list[Path],
67 definition: Path,
68 mode: str,
69 output_dir: Path,
70 group: str | None = None,
71 filename_as_id: bool = False,
72) -> dict:
73 """Analyze all supplied samples; reject incomplete workspace groups and ID collisions."""
74 import flowkit as fk
76 if mode not in {"gatingml", "workspace"}:
77 raise ValueError("mode must be gatingml or workspace")
78 if (mode == "workspace") != (group is not None):
79 raise ValueError("A group is required for workspace mode and only valid there")
80 if output_dir.exists():
81 raise FileExistsError(f"Output directory already exists: {output_dir}")
82 if not fcs_paths:
83 raise ValueError("Supply at least one FCS file")
84 for path in [definition, *fcs_paths]:
85 if not path.is_file():
86 raise ValueError(f"Expected a local file: {path}")
88 # Hash the actual inputs, not paths or metadata supplied by a workspace.
89 definition_record = file_record(definition)
90 samples = []
91 records = []
92 seen_ids = set()
93 for path in fcs_paths:
94 record = file_record(path)
95 sample = fk.Sample(
96 path,
97 filename_as_id=filename_as_id,
98 use_flowjo_labels=(mode == "workspace"),
99 )
100 if sample.id in seen_ids:
101 raise ValueError(f"Duplicate sample ID: {sample.id!r}; check FCS $FIL values")
102 if sample.event_count == 0:
103 raise ValueError(f"Sample has no events: {sample.id!r}")
104 seen_ids.add(sample.id)
105 samples.append(sample)
106 record.update(sample_id=sample.id, event_count=sample.event_count)
107 records.append(record)
109 if mode == "gatingml":
110 strategy = fk.parse_gating_xml(str(definition))
111 if not strategy.get_gate_ids():
112 raise ValueError("GatingML strategy has no gates")
113 analysis = fk.Session(strategy, fcs_samples=samples)
114 analysis.analyze_samples(use_mp=False, cache_events=False)
115 report = analysis.get_analysis_report()
116 else:
117 analysis = fk.Workspace(
118 str(definition),
119 fcs_samples=samples,
120 # Retain missing-sample gate metadata so completeness can be checked.
121 load_missing_file_data=True,
122 find_fcs_files_from_wsp=False,
123 )
124 expected = set(analysis.get_sample_ids(group_name=group, loaded_only=False))
125 if not expected:
126 raise ValueError(f"Workspace group is empty: {group!r}")
127 if expected != seen_ids:
128 raise ValueError(
129 "FCS inputs must exactly match the selected workspace group; "
130 f"missing IDs={sorted(expected - seen_ids)}, "
131 f"unexpected IDs={sorted(seen_ids - expected)}"
132 )
133 for sample_id in sorted(expected):
134 if not analysis.get_gate_ids(sample_id):
135 raise ValueError(f"Workspace sample has no gates: {sample_id!r}")
136 analysis.analyze_samples(group_name=group, use_mp=False, cache_events=False)
137 report = analysis.get_analysis_report(group_name=group)
139 if report.empty or set(report["sample_id"]) != seen_ids:
140 raise ValueError("Analysis did not produce gate results for every input sample")
141 report = add_denominators(report, {sample.id: sample.event_count for sample in samples})
142 # Preserve path components without inventing an additional separator.
143 for column in ("gate_path", "population_path"):
144 report[column] = report[column].map(lambda path: json.dumps(list(path)))
145 manifest = {
146 "mode": mode,
147 "group": group,
148 "filename_as_id": filename_as_id,
149 "definition": definition_record,
150 "samples": records,
151 "versions": {
152 package: version(package)
153 for package in ("flowkit", "flowio", "flowutils", "numpy", "pandas")
154 },
155 "report": "gate_report.csv",
156 "percentages": {
157 "absolute_percent": "100 * gate count / total sample events",
158 "relative_percent": "100 * gate count / parent population count",
159 "empty_parent": "relative_percent is blank; relative_percent_defined is false",
160 },
161 "population_identity": "sample_id plus population_path (JSON array including gate name and quadrant owner)",
162 "event_processing": "FlowKit FCS preprocessing; compensation and transforms from gate definitions",
163 "use_mp": False,
164 }
165 output_dir.mkdir(parents=True, exist_ok=False)
166 report.to_csv(output_dir / "gate_report.csv", index=False)
167 (output_dir / "provenance.json").write_text(
168 json.dumps(manifest, indent=2, allow_nan=False) + "\n", encoding="utf-8"
169 )
170 return manifest
173def main(argv: list[str] | None = None) -> int:
174 parser = argparse.ArgumentParser(description=__doc__)
175 definitions = parser.add_mutually_exclusive_group(required=True)
176 definitions.add_argument("--gatingml", type=Path, help="GatingML 2.0 XML strategy")
177 definitions.add_argument("--workspace", type=Path, help="FlowJo 10 .wsp file")
178 parser.add_argument("--fcs", nargs="+", type=Path, required=True, help="Explicit local FCS files")
179 parser.add_argument("--group", help="Workspace group; all its FCS files must be supplied")
180 parser.add_argument("--output-dir", type=Path, required=True, help="New results directory")
181 parser.add_argument(
182 "--filename-as-id", action="store_true",
183 help="Use current file basenames instead of FCS $FIL metadata for sample IDs",
184 )
185 args = parser.parse_args(argv)
186 if bool(args.workspace) != (args.group is not None):
187 parser.error("--workspace requires --group; --group is not valid with --gatingml")
188 try:
189 manifest = analyze(
190 fcs_paths=args.fcs,
191 definition=args.workspace or args.gatingml,
192 mode="workspace" if args.workspace else "gatingml",
193 output_dir=args.output_dir,
194 group=args.group,
195 filename_as_id=args.filename_as_id,
196 )
197 except (ValueError, KeyError, OSError) as exc:
198 parser.exit(2, f"Analysis failed: {exc}\n")
199 print(f"Analyzed {len(manifest['samples'])} samples: {args.output_dir / 'gate_report.csv'}")
200 return 0
203if __name__ == "__main__":
204 raise SystemExit(main())