Optimal coloring (Yemini et al. 2021) on the atlas of dandiset 000715¶
Paper. Yemini et al. (2021), Cell 184:272-288, doi:10.1016/j.cell.2020.12.012; claims written from the accepted author manuscript (PMC10494711) with Methods S1 and S2.
Data. Neighbors come from the statistical atlas of neuron positions fit on the ten heads of
DANDI:000715, version 0.241009.1514, the
same atlas and adjacency as ms2-neuropal-color-violations in reproduce_yemini_2021.ipynb (the
version is chosen there). The reporters and their expression come from data.mat of the paper's
Optimal_Coloring code (github.com/Eviatar/Optimal_Coloring at 990c597), an outside input: the
claim is data: with_inputs. The adjacency stored in the same data.mat, the paper's own input
for its 235 head and tail neurons, is run as a reference.
| Claim | Status here |
|---|---|
| fig7de-optimal-coloring-violations | tested, head atlas (with the outside reporter table) |
The dandiset's other claims are in reproduce_yemini_2021.ipynb.
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import neuropal as npl
pd.set_option("display.width", 200)
started = time.time()
worms = npl.load_all()
oc = npl.load_optimal_coloring()
print(f"loaded {len(worms)} files from dandiset {npl.DANDISET_ID} version {npl.DANDISET_VERSION} and "
f"Optimal_Coloring data.mat in {time.time() - started:.0f} s")
loaded 10 files from dandiset 000715 version 0.241009.1514 and Optimal_Coloring data.mat in 15 s
Data check¶
feats = [npl.worm_features(w) for w in worms]
print("annotated neurons per worm:", {w["subject"]: len(f[0]) for w, f in zip(worms, feats)})
names, mu, sigma, atlas_info = npl.fit_atlas(feats, seed=0)
print(f"head atlas: {len(names)} neurons, worms per neuron {atlas_info['counts'].min()} to {atlas_info['counts'].max()}, "
f"positions {mu[:, :3].min():.1f} to {mu[:, :3].max():.1f} um")
BA = oc["BrainAtlas"]
print(f"\ndata.mat BrainAtlas table: {BA.shape[0]} neurons x {BA.shape[1] - 3} reporter columns (after 3 descriptive columns "
f"{list(BA.columns[:3])}); values {sorted(np.unique(BA.iloc[:, 3:].to_numpy()).tolist())}")
print(f"data.mat neurons to color: {len(oc['neurons'])}; all in the table: {set(oc['neurons']) <= set(BA.Neuron.astype(str))}")
R_ref, rep_ref = npl.reporter_matrix(BA, oc["neurons"])
print(f"reporters with a distinct pattern over these {len(oc['neurons'])} neurons: {R_ref.shape[1]} "
f"(Methods S2 section 4.1: 1192 reporters, unique = 654)")
A_ref = np.asarray(oc["neuron_prob_distances"]).astype(int)
print(f"data.mat adjacency: {A_ref.shape}, values {np.unique(A_ref).tolist()}, symmetric {bool((A_ref == A_ref.T).all())}, "
f"diagonal entries set {int(np.diag(A_ref).sum())}, adjacent ordered pairs {int(A_ref.sum())}")
annotated neurons per worm: {'11_YAaLR': 193, '20_YAaLR': 193, '38_YAaDV': 190, '55_YAaDV': 192, '56_YAaDV': 191, '62_YAaLR': 186, '64_YAaDV': 189, '7_YAaLR': 188, '70_YAaLR': 193, '76_YAaDV': 188}
head atlas: 198 neurons, worms per neuron 1 to 10, positions 5.4 to 187.1 um
data.mat BrainAtlas table: 302 neurons x 1189 reporter columns (after 3 descriptive columns ['NeuronClass', 'Neuron', 'NeuronType']); values [0.0, 1.0]
data.mat neurons to color: 235; all in the table: True
reporters with a distinct pattern over these 235 neurons: 654 (Methods S2 section 4.1: 1192 reporters, unique = 654)
data.mat adjacency: (235, 235), values [0, 1], symmetric False, diagonal entries set 0, adjacent ordered pairs 5655
# The neighbors of ms2-neuropal-color-violations: atlas neurons that data.mat colors, adjacency from
# computeProbAdjacency on the head atlas.
shared = [n for n in names if n in oc["neurons"]]
ia = [names.index(n) for n in shared]
A_head = npl.prob_adjacency(mu[ia, :3], sigma[ia][:, :3, :3])
np.fill_diagonal(A_head, 0)
R_head, rep_head = npl.reporter_matrix(BA, shared)
print(f"head atlas neurons in data.mat: {len(shared)} of {len(names)}; left out: {[n for n in names if n not in oc['neurons']]}")
print(f"head adjacency: {int(A_head.sum())} adjacent ordered pairs, symmetric {bool((A_head == A_head.T).all())}, "
f"neighbors per neuron {A_head.sum(1).min()} to {A_head.sum(1).max()} (median {int(np.median(A_head.sum(1)))})")
print(f"reporters with a distinct pattern over the {len(shared)} head neurons: {R_head.shape[1]}")
j = [oc["neurons"].index(n) for n in shared]
print(f"agreement with data.mat's adjacency on the same neurons: {100 * (A_head == A_ref[np.ix_(j, j)]).mean():.1f}% of entries")
head atlas neurons in data.mat: 192 of 198; left out: ['AMSOL', 'AMSOR', 'DA2', 'DB3', 'HMC', 'VD3'] head adjacency: 4065 adjacent ordered pairs, symmetric False, neighbors per neuron 0 to 70 (median 16) reporters with a distinct pattern over the 192 head neurons: 543 agreement with data.mat's adjacency on the same neurons: 96.1% of entries
fig7de-optimal-coloring-violations¶
The optimal-coloring algorithm, using WormBase reporters, finds NeuroPAL alternatives with 1.7% color violations with 39 reporters and 10.5% when restricted to 3 reporters.
Criterion: adjacency as in ms2-neuropal-color-violations; the Brain Atlas reporters of data.mat;
3 colors; margin 1/3; optimal_color_solver.m; the best of 10,000 runs at sparsity 100 for the
roughly 40-reporter case and the best of 7,500 runs for the 3-reporter case, at the sparsity whose
solutions use 3 reporters after the code's Otsu threshold. Violations are the code's mean shortfall
below the margin over neighboring pairs (find_optimal_colors.m), which decides the claim; the
Methods share of neighboring pairs below the margin is shown as a variant.
Details taken from the paper and its code:
- Three intensities is the margin of 1/3: "1/3 as our minimum color margin (corresponding to discriminating 3 levels of color intensity)" (Methods S2 section 4.2). The solver has no other intensity setting.
- 300 iterations: optimal_color_parameters.m (
iterations = 300) and Methods S2 section 4.2 ("about 25 seconds for 300 iterations in total"). - The solution kept from a set of runs is the one with the fewest violations: "chosen as those with the minimal color violations, among 10,000 simulations" and "among 7,500 simulations" (Methods S2 Figure 1 legend). The figure's labels give the reporter counts of the chosen solutions (Reporters=39, Reporters=3). The best run over all 10,000 is the main reading of the first value, and the best run keeping exactly 39 reporters a variant; for the second value the best run keeping exactly 3 reporters is the main reading, and the best over all runs a variant.
- Each solution's colors are those of find_optimal_colors.m: reporters whose summed weight exceeds the Otsu threshold of X are kept, and R(:, kept) X(kept, :) is scaled by each channel's maximum.
Runs are drawn in batches of 250 from one generator per part, with fixed seeds.
K, MARGIN, ITERS, BATCH = 3, 1 / 3, 300, 250
def run_many(A, R, sparsity, runs, rng):
out = []
for b in range(0, runs, BATCH):
X = npl.optimal_color_solver(A, R, K, MARGIN, sparsity, ITERS, min(BATCH, runs - b), rng)
out.append(npl.solution_violations(X, R, A, MARGIN))
t = pd.DataFrame(np.concatenate(out), columns=["reporters", "share", "shortfall"])
t["reporters"] = t.reporters.astype(int)
t[["share", "shortfall"]] *= 100
return t
def best(t, label, reporters=None):
s = t if reporters is None else t[t.reporters == reporters]
if len(s) == 0:
return dict(selection=label, runs=0)
r = s.loc[s.shortfall.idxmin()]
return dict(selection=label, runs=len(s), reporters_of_best=int(r.reporters), shortfall_pct=r.shortfall,
methods_share_pct=r.share, lowest_share_pct_any_run=s.share.min())
t0 = time.time()
r40 = run_many(A_head, R_head, 100, 10_000, np.random.default_rng(7140))
print(f"sparsity 100, head atlas: {len(r40)} runs in {time.time() - t0:.0f} s")
print("reporters kept per run: median", int(r40.reporters.median()), ", range", r40.reporters.min(), "to", r40.reporters.max())
print(f"shortfall over runs: median {r40.shortfall.median():.2f}%, 5th to 95th percentile "
f"{r40.shortfall.quantile(0.05):.2f} to {r40.shortfall.quantile(0.95):.2f}%, range {r40.shortfall.min():.2f} to {r40.shortfall.max():.2f}%")
tab40 = pd.DataFrame([best(r40, "best of all runs (main)"), best(r40, "best run with 39 reporters", 39),
best(r40, "best run with 40 reporters", 40)])
print(tab40.round(3).to_string(index=False))
sparsity 100, head atlas: 10000 runs in 637 s
reporters kept per run: median 30 , range 12 to 54
shortfall over runs: median 2.76%, 5th to 95th percentile 2.16 to 3.59%, range 1.72 to 9.27%
selection runs reporters_of_best shortfall_pct methods_share_pct lowest_share_pct_any_run
best of all runs (main) 10000 27 1.717 12.964 10.923
best run with 39 reporters 205 39 1.859 14.367 14.145
best run with 40 reporters 129 40 1.979 15.301 13.850
grid = [250, 500, 1000, 1500, 2000, 3000]
rng = np.random.default_rng(7141)
scan = {s: run_many(A_head, R_head, s, 100, rng) for s in grid}
rows = [dict(sparsity=s, runs=len(t), share_of_runs_with_3_reporters=(t.reporters == 3).mean(),
median_reporters=int(t.reporters.median()), reporter_range=f"{t.reporters.min()} to {t.reporters.max()}",
best_3_reporter_shortfall_pct=t[t.reporters == 3].shortfall.min() if (t.reporters == 3).any() else np.nan)
for s, t in scan.items()]
scan_tab = pd.DataFrame(rows)
print(scan_tab.round(3).to_string(index=False))
sp3 = int(scan_tab.sparsity[scan_tab.share_of_runs_with_3_reporters.idxmax()])
print(f"\nsparsity whose runs most often keep 3 reporters: {sp3}")
sparsity runs share_of_runs_with_3_reporters median_reporters reporter_range best_3_reporter_shortfall_pct
250 100 0.00 15 8 to 24 NaN
500 100 0.00 8 4 to 13 NaN
1000 100 0.47 4 3 to 8 9.897
1500 100 0.95 3 2 to 5 9.897
2000 100 0.94 3 2 to 3 9.897
3000 100 0.89 3 2 to 3 15.178
sparsity whose runs most often keep 3 reporters: 1500
t0 = time.time()
r3 = run_many(A_head, R_head, sp3, 7_500, np.random.default_rng(7142))
print(f"sparsity {sp3}, head atlas: {len(r3)} runs in {time.time() - t0:.0f} s; reporters kept:",
r3.reporters.value_counts().sort_index().to_dict())
s3 = r3[r3.reporters == 3].shortfall
print(f"shortfall of 3-reporter runs: median {s3.median():.2f}%, 5th to 95th percentile {s3.quantile(0.05):.2f} to {s3.quantile(0.95):.2f}%")
tab3 = pd.DataFrame([best(r3, "best run with 3 reporters (main)", 3), best(r3, "best of all runs")])
print(tab3.round(3).to_string(index=False))
sparsity 1500, head atlas: 7500 runs in 491 s; reporters kept: {2: 20, 3: 7106, 4: 324, 5: 34, 6: 11, 7: 2, 8: 3}
shortfall of 3-reporter runs: median 14.21%, 5th to 95th percentile 9.90 to 20.69%
selection runs reporters_of_best shortfall_pct methods_share_pct lowest_share_pct_any_run
best run with 3 reporters (main) 7106 3 6.945 20.836 20.836
best of all runs 7500 7 4.935 18.204 18.204
# The paper's 7,500 simulations span 2 to 6 colors, 1 to 3 intensities and 2 to 71 reporters
# (Methods S2 Figure 6 legend), so its 3-color, 3-intensity, 3-reporter solution may be the best of
# fewer runs. The best of random subsets of the 3-reporter runs shows how the best depends on the
# number of runs (1,000 subsets per size).
rng = np.random.default_rng(7145)
v = s3.to_numpy()
rows = []
for n in (10, 30, 100, 300, 1000, 3000, len(v)):
m = np.array([v[rng.choice(len(v), n, replace=False)].min() for _ in range(1000)])
rows.append(dict(runs=n, median_best_pct=np.median(m), p5_best_pct=np.percentile(m, 5), p95_best_pct=np.percentile(m, 95),
share_of_subsets_within_9_45_to_11_55=np.mean((m >= 9.45) & (m <= 11.55))))
print(pd.DataFrame(rows).round(3).to_string(index=False))
print(f"\n3-reporter runs at or below 11.55%: {(v <= 11.55).mean():.3f}; at or below 10.5%: {(v <= 10.5).mean():.3f}; "
f"below 9.45%: {(v < 9.45).mean():.4f} ({int((v < 9.45).sum())} runs); distinct values below 9.45%: {np.unique(v[v < 9.45].round(3)).tolist()}")
runs median_best_pct p5_best_pct p95_best_pct share_of_subsets_within_9_45_to_11_55 10 9.897 9.717 9.897 0.951 30 9.897 8.380 9.897 0.878 100 9.717 7.831 9.897 0.657 300 8.380 6.945 9.717 0.255 1000 7.831 6.945 8.380 0.011 3000 6.945 6.945 8.069 0.000 7106 6.945 6.945 6.945 0.000 3-reporter runs at or below 11.55%: 0.395; at or below 10.5%: 0.360; below 9.45%: 0.0042 (30 runs); distinct values below 9.45%: [6.945, 7.733, 7.831, 8.069, 8.38]
Reference: the same runs on the paper's own adjacency in data.mat (235 head and tail neurons),
1,000 runs at each sparsity to keep the notebook within its time limit.
t0 = time.time()
ref40 = run_many(A_ref, R_ref, 100, 1_000, np.random.default_rng(7143))
ref3 = run_many(A_ref, R_ref, sp3, 1_000, np.random.default_rng(7144))
print(f"data.mat adjacency: {len(ref40) + len(ref3)} runs in {time.time() - t0:.0f} s")
print(f"sparsity 100: reporters kept median {int(ref40.reporters.median())}, range {ref40.reporters.min()} to {ref40.reporters.max()}; "
f"shortfall median {ref40.shortfall.median():.2f}%")
print(f"sparsity {sp3}: reporters kept", ref3.reporters.value_counts().sort_index().to_dict())
print(pd.DataFrame([best(ref40, "sparsity 100, best of all runs"), best(ref40, "sparsity 100, best run with 39 reporters", 39),
best(ref3, f"sparsity {sp3}, best run with 3 reporters", 3), best(ref3, f"sparsity {sp3}, best of all runs")]
).round(3).to_string(index=False))
data.mat adjacency: 2000 runs in 195 s
sparsity 100: reporters kept median 34, range 19 to 53; shortfall median 2.67%
sparsity 1500: reporters kept {2: 3, 3: 949, 4: 16, 5: 15, 6: 16, 7: 1}
selection runs reporters_of_best shortfall_pct methods_share_pct lowest_share_pct_any_run
sparsity 100, best of all runs 1000 31 1.737 12.042 12.042
sparsity 100, best run with 39 reporters 29 39 1.993 16.994 15.296
sparsity 1500, best run with 3 reporters 949 3 8.712 26.136 26.136
sparsity 1500, best of all runs 1000 5 5.199 15.915 15.915
fig, axes = plt.subplots(1, 2, figsize=(10, 3.2))
g = r40.groupby("reporters").shortfall
ax = axes[0]
ax.plot(g.min().index, g.min().values, "o-", ms=3, label="best run")
ax.plot(g.median().index, g.median().values, "o-", ms=3, color="0.6", label="median run")
ax.axhline(1.7, color="C3", lw=0.8, ls="--", label="paper, 1.7%")
ax.set_xlabel("reporters kept"); ax.set_ylabel("violations, mean shortfall (%)"); ax.set_title("sparsity 100, 10,000 runs, head atlas")
ax.legend(fontsize=8)
ax = axes[1]
ax.hist(s3, bins=40, color="0.4")
ax.axvline(10.5, color="C3", lw=0.8, ls="--", label="paper, 10.5%")
ax.set_xlabel("violations, mean shortfall (%)"); ax.set_ylabel("runs"); ax.set_title(f"3-reporter runs at sparsity {sp3}")
ax.legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"left: reporters {g.min().index.min()} to {g.min().index.max()}, best {g.min().min():.2f} to {g.min().max():.2f}%, "
f"median {g.median().min():.2f} to {g.median().max():.2f}%; right: {len(s3)} runs, {s3.min():.2f} to {s3.max():.2f}%; "
f"no point outside the axes")
print(f"notebook time {(time.time() - started) / 60:.1f} min")
left: reporters 12 to 54, best 1.72 to 5.57%, median 1.99 to 5.98%; right: 7106 runs, 6.95 to 20.69%; no point outside the axes notebook time 22.9 min
Result: blocked: the 3-reporter value depends on the sparsity, which the paper does not give. On the dandiset's head atlas, the best of 10,000 runs at sparsity 100 has a mean shortfall inside 1.7 +/- 0.17% (1.71 to 1.72% in the runs made for this attempt, keeping 27 to 30 reporters), and the best run keeping exactly 39 reporters is also inside (1.83 to 1.86%); the values printed above are this run's, and they differ in the second decimal between machines with the same seeds. For 3 reporters the best of the 7,500 runs at sparsity 1500 (7,106 of them keep 3 reporters) has 6.95%, below the tolerance of 9.45 to 11.55%. The 3-reporter solutions take few distinct values; 9.90% is their 5th percentile and only 30 runs (0.4%) fall below 9.45%, so the best of 10 runs lies within the tolerance in 95% of random subsets, the best of 100 runs in 66% and the best of 1,000 runs in 1%. On the paper's own adjacency in data.mat (1,000 runs each) the best is 1.74 to 1.81% at sparsity 100 and 8.71% for 3 reporters. The Methods share of neighboring pairs below the margin, a variant, is 13.65% and 20.84% for the two best solutions on the head atlas, far above both expected values.
The criterion reads the sparsity for 3 reporters as the sparsity whose solutions use 3 reporters,
which fits 1500, 2000 and 3000 of the scan equally well. optimal_coloring_sparsity.ipynb runs the
others. With the criterion's 7,500 runs the best 3-reporter value is 8.38% at 2000, below the
tolerance, and 9.90% at 3000, inside it (6,253 runs keep 3 reporters, none below 9.45%); at 1000,
where about 56% of runs keep 3 reporters, the best of 1,000 runs is 7.83%. The outcome would be
partially reproduced at 1500 and 2000 and reproduced at 3000. No sparsity is better supported than
the others, so the claim is blocked (methods_underspecified) with every result shown, as for the
other claims of this attempt whose readings disagree. The first record gave partially reproduced
at 1500; the change was proposed after verification round 2 and agreed in round 3. Whether the
paper's 3-reporter value is the best of 7,500 runs of this one setting, as the criterion reads it,
or of the smaller part of its 7,500 simulations that had this setting, is asked of a reader.
Summary¶
| Claim | Outcome |
|---|---|
| fig7de-optimal-coloring-violations | blocked: the 3-reporter value is below the tolerance at sparsities 1500 and 2000 and inside it at 3000, and the paper does not give the sparsity (head atlas, with the outside reporter table) |
How the analysis differs from the paper's: the neighbors come from the head atlas of the dandiset
(192 neurons), since its tails are not there, where the paper used 235 head and tail neurons; the
solver is a Python port of optimal_color_solver.m with its own random starts; the sparsity for 3
reporters is not given, so 1500, the sparsity of the scan whose runs most often keep 3 reporters,
is run here and 1000, 2000 and 3000 in optimal_coloring_sparsity.ipynb; the runs on the paper's
own adjacency are 1,000 per setting, not 10,000 and 7,500; sparsity 1000 has 1,000 runs. The
solver's results at sparsity 100 differ in the second decimal between machines (the best was 1.71
to 1.72% in the runs made for this attempt), though seeds are fixed.
What was found about the dataset is reported in reproduce_yemini_2021.ipynb; this notebook uses
only the dandiset's annotated positions, through the atlas.