Saving a strategy and applying it to a batch (Gating-ML 2.0)#
The reason to model a gating strategy as an object that does not belong to any one sample: one strategy, ninety files. This notebook saves a strategy, reads it back, and runs it across a batch — and then does the thing that actually matters, which is comparing the resulting frequencies across samples and noticing where the strategy stopped fitting.
Continues from gating.
# No data of your own is needed. ov.datasets.flow_demo_fcs() writes a SIMULATED
# FCS 3.1 file into the current working directory: 59,600 events, nine known
# populations, a real $SPILLOVER keyword, and a detector noise floor that pushes
# dim events below zero — which is why a log display axis will not do.
import numpy as np, pathlib
import omicverse as ov
ov.style() # omicverse's plotting defaults — call this once, up front
path = ov.datasets.flow_demo_fcs()
print(path)
🔬 Starting plot initialization...
🧬 Detecting GPU devices…
✅ Apple Silicon MPS detected
• [MPS] Apple Silicon GPU - Metal Performance Shaders available
____ _ _ __
/ __ \____ ___ (_)___| | / /__ _____________
/ / / / __ `__ \/ / ___/ | / / _ \/ ___/ ___/ _ \
/ /_/ / / / / / / / /__ | |/ / __/ / (__ ) __/
\____/_/ /_/ /_/_/\___/ |___/\___/_/ /____/\___/
🔖 Version: 2.2.4 📚 Tutorials: https://omicverse.readthedocs.io/
✅ plot_set complete.
flow_demo.fcs
import omicverse as ov
import matplotlib.pyplot as plt
import numpy as np
adata = ov.datasets.flow_demo() # = read_fcs(flow_demo_fcs()) + obs['population']
ov.flow.compensate(adata)
lin = ov.flow.Linear(t=262144.0)
lg = ov.flow.Logicle(t=262144.0, m=4.5, w=1.0, a=0.0)
adata
AnnData object with n_obs × n_vars = 59600 × 8
obs: 'sample', 'population'
var: 'n', 'channel', 'marker', 'PnB', 'PnE', 'PnG', 'PnR'
uns: 'meta', 'fcs', 'flow_demo', 'flow'
layers: 'uncompensated'
def at(v):
"""A data value expressed on the logicle scale — how a threshold is really
written. Typing 0.42 instead of at(600) is how gates drift."""
return float(lg.apply(np.array([float(v)]))[0])
cells = ov.flow.PolygonGate(
name="Cells", dims=("FSC-A", "SSC-A"),
transforms={"FSC-A": lin, "SSC-A": lin},
vertices=np.array([[0.12,0.01],[0.12,0.17],[0.40,0.32],
[0.68,0.32],[0.68,0.09],[0.38,0.01]]))
singlets = ov.flow.PolygonGate(
name="Singlets", dims=("FSC-A", "FSC-H"),
transforms={"FSC-A": lin, "FSC-H": lin},
vertices=np.array([[0.05,0.03],[0.05,0.10],[0.75,0.86],[0.75,0.69]]))
live = ov.flow.RectangleGate(
name="Live", dims=("Viability",), transforms={"Viability": lg},
bounds=((None, at(1200)),))
cd3 = ov.flow.RectangleGate(
name="CD3+", dims=("CD3",), transforms={"CD3": lg},
bounds=((at(600), None),))
quad = ov.flow.QuadrantGate(
name="CD4/CD8", dims=("CD4", "CD8"), transforms={"CD4": lg, "CD8": lg},
dividers=(at(600), at(600)),
quadrant_names=("DN", "CD4+CD8-", "CD4-CD8+", "DP"))
gs = ov.flow.GatingStrategy("T cell panel")
gs.add_gate(cells)
gs.add_gate(singlets, parent="Cells")
gs.add_gate(live, parent="Singlets")
gs.add_gate(cd3, parent="Live")
gs.add_gate(quad, parent="CD3+")
print(gs.tree())
root
└ Cells
└ Singlets
└ Live
└ CD3+
└ CD4/CD8
└ DN
└ CD4+CD8-
└ CD4-CD8+
└ DP
Round-tripping through a dict#
to_dict is JSON-able and is the seam a gate editor works against — a UI renders
and edits these dicts and never needs the Python objects.
import json
d = gs.to_dict()
print(json.dumps(d["gates"][3], indent=2)[:520], "...")
{
"kind": "rectangle",
"name": "CD3+",
"dims": [
"CD3"
],
"transforms": {
"CD3": {
"type": "logicle",
"t": 262144.0,
"w": 1.0,
"m": 4.5,
"a": 0.0
}
},
"bounds": [
[
0.34689636224578013,
null
]
],
"parent": "Live"
} ...
restored = ov.flow.GatingStrategy.from_dict(d)
print(restored)
print(restored.tree())
GatingStrategy('T cell panel', 5 gates)
root
└ Cells
└ Singlets
└ Live
└ CD3+
└ CD4/CD8
└ DN
└ CD4+CD8-
└ CD4-CD8+
└ DP
The round trip has to preserve the transforms, not just the geometry, or a restored gate lands somewhere else entirely. Check it rather than assume it:
a = gs.apply(adata.copy())
b = restored.apply(adata.copy())
identical = all(np.array_equal(a.masks[k], b.masks[k]) for k in a.masks)
print(f"populations: {sorted(a.masks) == sorted(b.masks)}")
print(f"every mask identical after the round trip: {identical}")
populations: True
every mask identical after the round trip: True
Gating-ML 2.0#
The dict is convenient but omicverse-specific. Gating-ML 2.0 is the ISAC interchange standard, and it is what FlowJo, Cytobank and flowCore read.
ov.flow.write_gatingml(gs, "strategy.xml")
xml = pathlib.Path("strategy.xml").read_text()
print(f"{len(xml):,} bytes")
print(xml[:700], "...")
4,757 bytes
<?xml version='1.0' encoding='utf-8'?>
<ns0:Gating-ML xmlns:ns0="http://www.isac-net.org/std/Gating-ML/v2.0/gating" xmlns:ns1="http://www.isac-net.org/std/Gating-ML/v2.0/transformations" xmlns:ns2="http://www.isac-net.org/std/Gating-ML/v2.0/datatypes" xmlns:gating="http://www.isac-net.org/std/Gating-ML/v2.0/gating" xmlns:transforms="http://www.isac-net.org/std/Gating-ML/v2.0/transformations" xmlns:data-type="http://www.isac-net.org/std/Gating-ML/v2.0/datatypes">
<ns1:transformation ns1:id="xf_FSC_negA_linear">
<ns1:flin ns1:T="262144.0" ns1:A="0.0" />
</ns1:transformation>
<ns1:transformation ns1:id="xf_SSC_negA_linear">
<ns1:flin ns1:T="262144.0" ns1:A="0.0" />
</ns1:transfo ...
from_file = ov.flow.read_gatingml("strategy.xml")
c = from_file.apply(adata.copy())
print(from_file)
print("masks identical to the original:",
all(np.array_equal(a.masks[k], c.masks[k]) for k in a.masks))
GatingStrategy('strategy', 5 gates)
masks identical to the original: True
One strategy, several samples#
Now the point of all of it. Three files with different biology — the second has lost most of its CD4 compartment, the third is a lower-quality run — gated by the same strategy object.
All three are synthetic: they are re-syntheses of the same simulated panel with a population thinned, not real patients. The frequencies below are therefore a demonstration of the mechanism and nothing else.
def variant(seed, thin=None):
"""A second sample: re-simulate with a different seed, optionally thinning
a population to mimic a different patient."""
ad = ov.datasets.flow_demo(seed=seed)
ov.flow.compensate(ad)
if thin:
rng = np.random.default_rng(seed + 100)
keep = np.ones(ad.n_obs, bool)
for pop, drop in thin.items():
idx = np.flatnonzero((ad.obs["population"] == pop).to_numpy())
keep[rng.choice(idx, int(len(idx) * drop), replace=False)] = False
ad = ad[keep].copy()
return ad
batch = {
"healthy": adata,
"CD4-depleted": variant(1, thin={"CD4 T": 0.85}),
"smaller run": variant(2, thin={"CD4 T": 0.2, "CD8 T": 0.2}),
}
{k: v.n_obs for k, v in batch.items()}
{'healthy': 59600, 'CD4-depleted': 45150, 'smaller run': 54200}
import pandas as pd
rows = []
for name, sample in batch.items():
r = gs.apply(sample.copy())
s = r.stats().set_index("population")["freq_parent"]
rows.append(s.rename(name))
freqs = pd.concat(rows, axis=1)
freqs.style.format("{:.1%}")
| healthy | CD4-depleted | smaller run | |
|---|---|---|---|
| population | |||
| CD3+ | 69.5% | 53.3% | 65.0% |
| CD4+CD8- | 58.8% | 17.7% | 57.9% |
| CD4-CD8+ | 34.6% | 69.2% | 34.0% |
| CD4/CD8 | 100.0% | 100.0% | 100.0% |
| Cells | 81.5% | 75.6% | 79.7% |
| DN | 4.8% | 9.7% | 6.0% |
| DP | 1.7% | 3.5% | 2.1% |
| Live | 92.8% | 89.4% | 91.8% |
| Singlets | 92.2% | 88.9% | 91.2% |
Easier to read as a picture — and this is the figure a batch actually produces, one bar group per population, one colour per sample:
fig, ax = plt.subplots(figsize=(8.6, 3.2))
order = ["Cells", "Singlets", "Live", "CD3+", "CD4+CD8-", "CD4-CD8+", "DN", "DP"]
sub = freqs.loc[order]
x, w = np.arange(len(order)), 0.26
for i, col in enumerate(sub.columns):
ax.bar(x + (i - 1) * w, sub[col], w, label=col)
ax.set_xticks(x, labels=order, rotation=30, ha="right")
ax.set_ylabel("frequency of parent")
ax.set_title("one strategy, three samples", fontsize=9)
ax.legend(fontsize=8, frameon=False)
plt.tight_layout()
The CD4 depletion is unmissable in CD4+CD8- — and note that
freq_parent is what makes the comparison legitimate, since the three files do
not have the same number of events.
Side by side, on the same axes, with the same gate:
fig, axes = plt.subplots(1, 3, figsize=(14, 4.2))
for ax, (name, sample) in zip(axes, batch.items()):
r = gs.apply(sample.copy())
ov.flow.biaxial(sample, "CD4", "CD8", gates=[quad], result=r,
population="CD3+", ax=ax, max_events=30000,
title=name)
plt.tight_layout()
Diffing two strategies#
Because a strategy is data, two of them can be compared — which is the other thing that is hard to do by hand. Here a colleague has tightened the viability threshold, three nodes above the populations anyone will report:
theirs = ov.flow.GatingStrategy.from_dict(gs.to_dict())
theirs.name = "T cell panel (revised)"
theirs.gates["Live"].bounds = ((None, at(300)),)
mine = gs.apply(adata.copy()).stats().set_index("population")["count"]
yours = theirs.apply(adata.copy()).stats().set_index("population")["count"]
diff = pd.DataFrame({"original": mine, "revised": yours})
diff["delta"] = diff["revised"] - diff["original"]
diff["pct"] = (diff["delta"] / diff["original"]).map("{:+.1%}".format)
diff
original revised delta pct
population
CD3+ 28895 27954 -941 -3.3%
CD4+CD8- 16999 16453 -546 -3.2%
CD4-CD8+ 9998 9674 -324 -3.2%
CD4/CD8 28895 27954 -941 -3.3%
Cells 48587 48587 0 +0.0%
DN 1398 1343 -55 -3.9%
DP 500 484 -16 -3.2%
Live 41574 40085 -1489 -3.6%
Singlets 44785 44785 0 +0.0%
One threshold, moved at a node nobody reports, and every population below it changes. This is the argument for keeping the strategy as an object that can be diffed rather than as a sequence of clicks: the number that moved is not the one that was edited.
What is next#
Everything so far is manual and auditable. The last notebook adds unsupervised clustering — and is explicit about what it does not replace.