Reproducing Yemini et al. (2021), NeuroPAL, from dandiset 000715¶
Paper. Yemini E, Lin A, Nejatbakhsh A, Varol E, Sun R, Mena GE, Samuel ADT, Paninski L, Venkatachalam V, Hobert O (2021). NeuroPAL: A Multicolor Atlas for Whole-Brain Neuronal Identification in C. elegans. Cell 184:272-288. doi:10.1016/j.cell.2020.12.012. The claims file was written from the accepted author manuscript in PubMed Central (PMC10494711) with Methods S1 and S2; the version of record was not read.
Dandiset. DANDI:000715, version 0.241009.1514 (doi:10.48324/dandi.000715/0.241009.1514). The paper cites no DANDI version (its data statement points to Zenodo). There are four published versions (0.240424.1847, 0.240611.1953, 0.240614.1942, 0.241009.1514) and the draft; all list the same ten assets with the same asset and blob identifiers and sizes, so their content is identical. The most recent published version is pinned.
What the dandiset holds. Ten NWB files, each one worm head: a five-channel volume and the hand-annotated neurons (one centre voxel and one name per neuron). The data check below shows that these are the heads of the paper's ten atlas worms (otIs669), whose head and tail volumes are also on Zenodo (doi:10.5281/zenodo.3906530); the tails are not in the dandiset. The files' own metadata say otherwise (strain OH16230, a spinning disk confocal, a microfluidic chip).
Outside inputs. The NeuroPAL colors of each neuron, from data.mat of the paper's
Optimal_Coloring code (github.com/Eviatar/Optimal_Coloring at 990c597), for
ms2-neuropal-color-violations (and the reporter table of the same file for fig7de, in
optimal_coloring.ipynb). Two head volumes of the Zenodo record, read only to check what the
dandiset's files are, and the hand annotations of all ten worms (*_ID.mat), compared with the
dandiset's. The ganglion of each neuron, from the CELL_ID neuron tables (Data/NeuroPAL/herm_data.mat
at d9f5df9), for the head part of fig7bc-ventral-ganglion-lowest.
| Claim | Status here |
|---|---|
| fig2cd-variability-differs | tested, head only |
| fig2c-ollr-more-variable-than-olll | tested |
| fig7bc-identification-accuracy | tested, head only |
| fig7bc-ventral-ganglion-lowest | blocked: tails not in the dandiset; head part reported |
| fig7b-eight-labels-above-90 | tested |
| ms1-uncertain-first-faster | tested, head only |
| ms1-position-only-accuracy | tested, head only |
| ms1-detection-accuracy | tested, head only |
| ms1-calibration | tested, head only |
| ms2-neuropal-color-violations | tested, head atlas (with outside colors) |
| fig7de-optimal-coloring-violations | tested in optimal_coloring.ipynb, head atlas (with outside reporters) |
| other claims (figures 1, 3 to 6, S1, S2) | not tested on the dandiset; see the record |
import time
import io
import warnings
import numpy as np
import pandas as pd
import requests
import scipy.io as sio
from scipy import stats
import matplotlib.pyplot as plt
import neuropal as npl
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 30)
pd.set_option("display.max_rows", 300)
warnings.filterwarnings("ignore", category=RuntimeWarning)
started = time.time()
worms = npl.load_all()
print(f"loaded {len(worms)} files from dandiset {npl.DANDISET_ID} version {npl.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s")
loaded 10 files from dandiset 000715 version 0.241009.1514 in 14 s
Data check¶
info = pd.DataFrame([{
"subject": w["subject"], "location": w["location"], "strain": w["strain"],
"device": w["device"][0], "chip": w["subject_description"], "born": w["date_of_birth"][:10],
"session": w["session_start"][:10], "shape": w["shape"], "grid_um": tuple(np.round(w["grid"], 4)),
"rois": len(w["labels"]), "voxels_per_roi": w["voxels_per_roi"].tolist(),
"named": sum(1 for l in w["labels"] if l), "unnamed": sum(1 for l in w["labels"] if not l),
"duplicate_names": len([l for l in set(w["labels"]) if l and w["labels"].count(l) > 1]),
} for w in worms])
print(info.to_string(index=False))
print("\nChannel order in every file:", worms[0]["channels"], "; RGBW indices:", worms[0]["rgbw"].tolist())
print("Paper (Methods S1, MICCAI atlas paper): 10 young-adult otIs669 (OH15262) worms, head and tail, "
"Zeiss LSM 880, about 0.2 x 0.2 x 0.8 um voxels, 195 head neurons annotated per worm.")
subject location strain device chip born session shape grid_um rois voxels_per_roi named unnamed duplicate_names 11_YAaLR head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (1011, 284, 42, 5) (0.2076, 0.2076, 0.75) 194 [1] 193 1 0 20_YAaLR head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (944, 309, 41, 5) (0.2076, 0.2076, 0.75) 194 [1] 193 1 0 38_YAaDV head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (921, 232, 38, 5) (0.2076, 0.2076, 0.75) 196 [1] 190 6 0 55_YAaDV head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (905, 218, 39, 5) (0.2078, 0.2078, 0.75) 192 [1] 192 0 0 56_YAaDV head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (935, 222, 37, 5) (0.2078, 0.2078, 0.75) 196 [1] 191 5 0 62_YAaLR head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (773, 239, 41, 5) (0.2076, 0.2076, 0.75) 189 [1] 186 3 0 64_YAaDV head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (831, 213, 41, 5) (0.2078, 0.2078, 0.75) 192 [1] 189 3 0 7_YAaLR head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (880, 264, 29, 5) (0.2077, 0.2077, 1.0215) 192 [1] 188 4 0 70_YAaLR head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (789, 219, 41, 5) (0.2076, 0.2076, 0.75) 193 [1] 193 0 0 76_YAaDV head OH16230 Spinning disk confocal NeuroPAL worm in microfluidic chip 2021-01-05 2021-01-07 (802, 226, 39, 5) (0.2077, 0.2077, 0.75) 193 [1] 188 5 0 Channel order in every file: ['mTagBFP2', 'Tag RFP-T', 'GFP-GCaMP', 'CyOFP1', 'mNeptune 2.5'] ; RGBW indices: [4, 3, 0, 1] Paper (Methods S1, MICCAI atlas paper): 10 young-adult otIs669 (OH15262) worms, head and tail, Zeiss LSM 880, about 0.2 x 0.2 x 0.8 um voxels, 195 head neurons annotated per worm.
Every file is a head (location), with one voxel per annotated neuron. The files say strain
OH16230 (the NeuroPAL;GCaMP6s strain of the whole-brain recordings), a spinning disk confocal, a
worm "in microfluidic chip", birth on 2021-01-05 and imaging on 2021-01-07, the same for all ten.
The paper's atlas worms are otIs669 (OH15262), imaged on a Zeiss LSM 880, and the paper was
published in January 2021. The subject identifiers (7_YAaLR, 38_YAaDV, ...) are those of the ten
otIs669 worms in the Zenodo record. The next cell compares two volumes with the Zenodo copies.
import fsspec, zipfile, h5py, remfile
ZENODO_ZIP = "https://zenodo.org/records/3906530/files/NeuroPAL%20data.zip?download=1"
zf = zipfile.ZipFile(fsspec.open(ZENODO_ZIP, block_size=2**22).open())
members = {m.split("/")[-1][:-4]: m for m in zf.namelist() if m.endswith(".mat") and "YAa" in m and "_ID" not in m}
print("Zenodo head volumes (otIs669 Young Adults):", sorted(members))
print("Dandiset subjects:", sorted(w["subject"] for w in worms))
print("Same set of identifiers:", sorted(members) == sorted(w["subject"] for w in worms))
for subject in ["7_YAaLR", "76_YAaDV"]:
m = sio.loadmat(io.BytesIO(zf.read(members[subject])), squeeze_me=True, struct_as_record=False,
variable_names=["data", "info"])
zen = m["data"]
w = next(w for w in worms if w["subject"] == subject)
asset = next(a for a in npl.list_assets() if subject in a["path"].replace("-", "_"))
with h5py.File(remfile.File(asset["url"]), "r") as f:
dan = f["acquisition/NeuroPALImageRaw/data"][()]
same = zen.shape == tuple(np.array(dan.shape)[[1, 0, 2, 3]]) and np.array_equal(np.transpose(zen, (1, 0, 2, 3)), dan)
print(f"{subject}: Zenodo {zen.shape}, dandiset {dan.shape}; identical voxel values after swapping x and y: {same}; "
f"Zenodo scale {np.round(m['info'].scale, 4)} against dandiset grid {np.round(w['grid'], 4)}; "
f"Zenodo RGBW (1-based) {m['info'].RGBW}, DIC/GFP channel (1-based) {m['info'].DIC}")
del zen, dan
# The annotations of all ten worms against the Zenodo *_ID.mat files (names and positions); reading
# the other eight volumes in full would take about 100 s each, so they are compared by annotations.
ids = {m.split("/")[-1][:-7]: m for m in zf.namelist() if m.endswith("_ID.mat") and "YAa" in m}
rows = []
for w in worms:
zlab, zpos = npl.zenodo_annotations(zf, ids[w["subject"]])
zvox = zpos[:, [1, 0, 2]] - 1 # MATLAB (y, x, z), 1-based, to the dandiset's (x, y, z), 0-based
same_names = sorted(zlab) == sorted(w["labels"])
order = [zlab.index(l) for l in w["labels"] if zlab.count(l) == 1]
stored = np.array([w["vox"][i] for i, l in enumerate(w["labels"]) if zlab.count(l) == 1])
diff = np.abs(zvox[order] - stored)
rows.append(dict(subject=w["subject"], view=w["subject"][-2:], zenodo_neurons=len(zlab), dandiset_neurons=len(w["labels"]),
same_names=same_names, compared_by_unique_name=len(order), max_abs_diff_voxels=round(float(diff.max()), 3),
equal_to_rounded_zenodo=bool((np.round(zvox[order]) == stored).all())))
print("\nAnnotations against the Zenodo *_ID.mat files:")
print(pd.DataFrame(rows).to_string(index=False))
Zenodo head volumes (otIs669 Young Adults): ['11_YAaLR', '20_YAaLR', '38_YAaDV', '55_YAaDV', '56_YAaDV', '62_YAaLR', '64_YAaDV', '70_YAaLR', '76_YAaDV', '7_YAaLR'] Dandiset subjects: ['11_YAaLR', '20_YAaLR', '38_YAaDV', '55_YAaDV', '56_YAaDV', '62_YAaLR', '64_YAaDV', '70_YAaLR', '76_YAaDV', '7_YAaLR'] Same set of identifiers: True
7_YAaLR: Zenodo (264, 880, 29, 5), dandiset (880, 264, 29, 5); identical voxel values after swapping x and y: True; Zenodo scale [0.2077 0.2077 1.0215] against dandiset grid [0.2077 0.2077 1.0215]; Zenodo RGBW (1-based) [5 4 1 2], DIC/GFP channel (1-based) 3
76_YAaDV: Zenodo (226, 802, 39, 5), dandiset (802, 226, 39, 5); identical voxel values after swapping x and y: True; Zenodo scale [0.2077 0.2077 0.75 ] against dandiset grid [0.2077 0.2077 0.75 ]; Zenodo RGBW (1-based) [5 4 1 2], DIC/GFP channel (1-based) 3
Annotations against the Zenodo *_ID.mat files: subject view zenodo_neurons dandiset_neurons same_names compared_by_unique_name max_abs_diff_voxels equal_to_rounded_zenodo 11_YAaLR LR 194 194 True 194 0.499 True 20_YAaLR LR 194 194 True 194 0.500 True 38_YAaDV DV 196 196 True 190 0.498 True 55_YAaDV DV 192 192 True 192 0.495 True 56_YAaDV DV 196 196 True 191 0.497 True 62_YAaLR LR 189 189 True 186 0.499 True 64_YAaDV DV 192 192 True 189 0.498 True 7_YAaLR LR 192 192 True 188 0.499 True 70_YAaLR LR 193 193 True 193 0.500 True 76_YAaDV DV 193 193 True 188 0.498 True
The two volumes checked are identical, voxel for voxel, to the Zenodo volumes of the otIs669 atlas
worms. For all ten worms the annotations hold the same names as the Zenodo *_ID.mat files, and
every stored voxel equals the Zenodo position, made 0-based and rounded (largest difference 0.5
voxel), comparing the neurons whose name occurs once. So the dandiset holds the heads of the paper's ten
atlas worms, and its strain, device, subject description and dates are wrong. The Zenodo files
mark the third channel (1-based) as the DIC (transmitted light) channel, which the dandiset calls
"GFP-GCaMP"; the next cell shows its values.
ch = pd.DataFrame([dict(subject=w["subject"], **{f"mean_{c}": round(float(v)) for c, v in zip(w["channels"], w["channel_mean"])},
**{f"max_{c}": int(v) for c, v in zip(w["channels"], w["channel_range"][1])}) for w in worms])
print(ch.to_string(index=False))
print("\nPanneuronal TagRFP-T (z-scored, 3x3x1 patch) at the stored voxel and shifted by -1 or +1 voxel on every axis:")
print(pd.DataFrame([dict(subject=w["subject"], **{f"offset_{k}": round(v, 2) for k, v in w["white_at_offset"].items()}) for w in worms]).to_string(index=False))
counts = pd.Series([l for w in worms for l in set(w["labels"]) if l]).value_counts()
print(f"\n{counts.size} distinct neuron names; held by all 10 worms: {(counts == 10).sum()}; by 8 or 9: {counts.between(8, 9).sum()}; "
f"by fewer than 8: {(counts < 8).sum()} ({', '.join(sorted(counts[counts < 8].index))})")
pos = pd.DataFrame([dict(subject=w["subject"], **{f"{a}_range_um": f"{w['pos'][:, i].min():.1f}-{w['pos'][:, i].max():.1f}" for i, a in enumerate('xyz')},
**{f"{c}_range": f"{w['rgbw_colors'][:, i].min():.2f}-{w['rgbw_colors'][:, i].max():.2f}" for i, c in enumerate('RGBW')}) for w in worms])
print("\nNeuron positions (um = voxel x grid) and colors (median z-score in the patch):")
print(pos.to_string(index=False))
sig = [(w["shape"], round(float(w["channel_mean"].sum()), 1), len(w["labels"])) for w in worms]
print("\nPairs of files with the same volume shape, channel means and neuron count (possible duplicates):",
[(worms[i]["subject"], worms[j]["subject"]) for i in range(10) for j in range(i + 1, 10) if sig[i] == sig[j]])
subject mean_mTagBFP2 mean_Tag RFP-T mean_GFP-GCaMP mean_CyOFP1 mean_mNeptune 2.5 max_mTagBFP2 max_Tag RFP-T max_GFP-GCaMP max_CyOFP1 max_mNeptune 2.5 11_YAaLR 702 1377 21313 624 312 65535 65535 51836 65535 65535 20_YAaLR 671 1414 19737 680 312 65535 65535 45460 65535 65535 38_YAaDV 893 2017 20592 880 468 65535 65535 47158 65535 65535 55_YAaDV 1016 1946 19504 763 485 65535 65535 51385 65535 65535 56_YAaDV 812 1810 20008 722 443 65535 65535 46540 65535 65535 62_YAaLR 864 1659 20673 632 416 65535 65535 48528 65535 65535 64_YAaDV 919 1911 19217 724 500 65535 65535 52398 65535 65535 7_YAaLR 737 1735 21015 769 393 65535 65535 50359 65535 65535 70_YAaLR 782 1773 20823 690 451 65535 65535 44093 65535 65535 76_YAaDV 878 1969 17302 789 476 65535 65535 41927 65535 65535 Panneuronal TagRFP-T (z-scored, 3x3x1 patch) at the stored voxel and shifted by -1 or +1 voxel on every axis: subject offset_-1 offset_0 offset_1 11_YAaLR 8.32 9.29 9.26 20_YAaLR 9.36 9.46 8.35 38_YAaDV 6.68 7.50 7.46 55_YAaDV 6.39 7.22 7.24 56_YAaDV 6.80 7.62 7.62 62_YAaLR 6.99 7.88 7.83 64_YAaDV 6.62 7.43 7.34 7_YAaLR 7.58 8.00 6.54 70_YAaLR 7.71 7.78 6.89 76_YAaDV 6.19 7.01 7.01 198 distinct neuron names; held by all 10 worms: 168; by 8 or 9: 23; by fewer than 8: 7 (AIYL, AS2, DA2, DB3, HMC, SIBDL, VD3) Neuron positions (um = voxel x grid) and colors (median z-score in the patch): subject x_range_um y_range_um z_range_um R_range G_range B_range W_range 11_YAaLR 18.1-203.4 9.8-48.4 3.8-24.8 -0.11-25.66 -0.15-35.76 -0.21-20.63 -0.04-12.44 20_YAaLR 12.0-191.6 13.1-61.6 5.2-24.0 -0.09-24.77 -0.12-32.11 -0.19-19.62 0.16-12.03 38_YAaDV 5.4-189.9 3.1-38.6 3.8-22.5 -0.13-20.70 -0.14-26.72 -0.21-15.64 0.23-9.53 55_YAaDV 13.3-185.3 6.9-37.8 8.2-27.8 -0.11-18.78 -0.13-28.60 -0.23-14.97 0.37-9.47 56_YAaDV 6.0-193.2 7.5-35.7 3.8-21.8 -0.13-20.68 -0.15-32.59 -0.21-17.16 0.09-10.20 62_YAaLR 7.7-158.0 9.3-43.2 4.5-24.8 -0.11-21.31 -0.03-36.56 -0.23-18.51 0.19-10.87 64_YAaDV 15.4-170.8 7.9-38.2 4.5-24.8 -0.13-18.75 -0.15-28.52 -0.21-15.40 0.14-9.53 7_YAaLR 8.9-178.2 18.5-49.0 4.1-24.5 -0.12-22.09 -0.17-29.52 -0.20-18.68 0.01-10.60 70_YAaLR 18.7-161.3 12.9-40.9 4.5-24.8 -0.13-20.19 -0.14-32.16 -0.19-17.02 0.02-10.08 76_YAaDV 10.0-162.7 8.1-40.9 3.8-23.2 -0.14-20.16 -0.17-28.49 -0.21-16.15 -0.07-9.51 Pairs of files with the same volume shape, channel means and neuron count (possible duplicates): []
The third channel has a mean near 20,000 in every file, far above the fluorescence channels, as expected of transmitted light, not of GCaMP. The panneuronal channel is brighter at the stored voxel than one voxel lower on every axis and about as bright as one voxel higher, so the voxel indices are read as 0-based centres, which the comparison with the Zenodo annotations above confirms. Positions span about 150 to 200 um along x (the head's length), as expected. No two files share shape, intensities and neuron count, so no recording is deposited twice. These are static volumes with no time axis or trial table, so there is no clock to check. About 3% of the stored neurons per worm have no name (empty label); they are left out. The names cover 198 neurons, of which 168 are named in all ten worms; the paper reports 195 head neurons annotated in each worm.
The statistical atlas (MS1 Algorithm 4)¶
The atlas is fit on the named neurons: position (um) and the R, G, B landmark colors (mNeptune2.5,
CyOFP1, mTagBFP2) read from the z-scored volume. Algorithm 4 as written alternates an affine
regression of the atlas means on each worm (equations 4.15 and 4.16) with the moments of the
aligned worms (equations 4.13 and 4.14) until the change is below a tolerance. On these data it
does not converge: the atlas is defined only up to an affine map, and the regression shrinks it
along its minor axes at every alternation (table below). The paper does not say how this was
avoided. The atlas here keeps the gauge of the initial worm (after each update the atlas is mapped
back, by the affine map that best fits its means to the initial worm's neurons) and stops after
npl.ATLAS_ITER = 20 alternations, where the position block has settled. This is an assumption;
the results are shown under other choices.
feats = [npl.worm_features(w) for w in worms]
def atlas_summary(names, mu, sigma):
coeff, m = npl._pca(mu[:, :3])
P = (mu[:, :3] - m) @ coeff
sd = np.sqrt(np.einsum("ai,nab,bi->ni", coeff, sigma[:, :3, :3], coeff))
csd = np.sqrt(np.diagonal(sigma[:, 3:, 3:], axis1=1, axis2=2))
return P.max(0) - P.min(0), np.median(sd, 0), np.median(csd, 0)
rows = []
for gauge in (False, True):
for it in (1, 5, 20, 100, 300):
names, mu, sigma, inf = npl.fit_atlas(feats, seed=0, max_iter=it, gauge=gauge)
ext, sd, csd = atlas_summary(names, mu, sigma)
rows.append(dict(gauge=gauge, iterations=it, converged=inf["converged"], change_mu=round(inf["dmu"], 4),
change_sigma=round(inf["dsig"], 4), extent_axes_um=np.round(ext, 1).tolist(),
median_sd_axes_um=np.round(sd, 2).tolist(), median_color_sd=np.round(csd, 2).tolist()))
print(pd.DataFrame(rows).to_string(index=False))
names, mu, sigma, atlas_info = npl.fit_atlas(feats, seed=0)
print(f"\nAtlas used: {len(names)} neurons, initial worm {worms[atlas_info['init_worm']]['subject']} (seed 0), "
f"{atlas_info['iterations']} alternations, gauge fixed; worms per neuron: min {atlas_info['counts'].min()}, "
f"median {int(np.median(atlas_info['counts']))}")
gauge iterations converged change_mu change_sigma extent_axes_um median_sd_axes_um median_color_sd False 1 False 89.8304 491.5575 [156.9, 22.5, 15.6] [2.19, 1.52, 1.06] [0.99, 1.54, 1.18] False 5 False 3.4915 8.8537 [156.4, 20.4, 14.4] [2.02, 1.38, 1.01] [0.92, 1.22, 1.0] False 20 False 2.0528 4.3402 [155.8, 15.4, 13.5] [2.06, 0.94, 0.92] [0.92, 0.96, 0.96] False 100 False 0.7055 1.1074 [150.0, 3.4, 0.9] [1.83, 0.24, 0.12] [0.88, 0.32, 0.75] False 300 False 0.3700 0.6165 [127.3, 0.0, 0.0] [1.55, 0.1, 0.1] [0.43, 0.16, 0.17] True 1 False 86.5017 568.6476 [161.3, 25.4, 16.9] [2.23, 1.69, 1.14] [1.14, 1.83, 1.43] True 5 False 0.9632 7.8928 [163.9, 25.4, 16.8] [2.11, 1.75, 1.15] [1.05, 1.74, 1.25] True 20 False 0.4846 5.6585 [166.6, 25.4, 16.8] [2.19, 1.75, 1.19] [1.06, 1.75, 1.28] True 100 False 0.9294 50.6137 [168.3, 25.5, 16.7] [2.21, 1.8, 1.2] [1.44, 2.36, 1.34] True 300 False 0.2896 54.7370 [168.4, 25.5, 16.7] [2.19, 1.8, 1.2] [2.29, 6.18, 3.29] Atlas used: 198 neurons, initial worm 70_YAaLR (seed 0), 20 alternations, gauge fixed; worms per neuron: min 1, median 10
Without the gauge the extent of the atlas along its second and third axes falls from about 22 and 16 um to 0 by 300 alternations; with it the position block holds (extents within a few percent from 5 to 300 alternations) while the color covariances keep growing slowly. The tolerance of 1e-6 is never reached.
fig2cd-variability-differs¶
In a probabilistic atlas built from 10 young-adult NeuroPAL worms, the positional variability of cell bodies differs across neuron types.
Criterion: per neuron, the standard deviation of its positional covariance along each principal axis of the atlas means; reproduced when along at least one axis the largest is more than twice the smallest, within the head and within the tail. Only heads are in the dandiset.
def axis_sd(mu, sigma):
coeff, m = npl._pca(mu[:, :3])
return np.sqrt(np.einsum("ai,nab,bi->ni", coeff, sigma[:, :3, :3], coeff)), coeff
def sd_table(mu, sigma, counts=None):
sd, _ = axis_sd(mu, sigma)
keep = np.ones(len(sd), bool) if counts is None else counts >= 8
s = sd[keep]
return s, pd.DataFrame({"axis": ["1 (major, anterior-posterior)", "2", "3"], "min_sd_um": s.min(0).round(3),
"median_sd_um": np.median(s, 0).round(3), "max_sd_um": s.max(0).round(3),
"max_over_min": (s.max(0) / s.min(0)).round(2),
"p90_over_p10": (np.percentile(s, 90, 0) / np.percentile(s, 10, 0)).round(2)})
s_all, t_all = sd_table(mu, sigma)
print(f"All {len(names)} atlas neurons (the ridge of {npl.RIDGE} um^2 bounds each SD below by {np.sqrt(npl.RIDGE):.2f} um):")
print(t_all.to_string(index=False))
s8, t8 = sd_table(mu, sigma, atlas_info["counts"])
print(f"\nNeurons named in at least 8 worms ({len(s8)}):")
print(t8.to_string(index=False))
print("\nVariants (neurons in at least 8 worms), largest over smallest SD per axis:")
for label, kw in [("initial worm 11_YAaLR", dict(init=0)), ("initial worm 38_YAaDV", dict(init=2)), ("5 alternations", dict(max_iter=5)),
("100 alternations", dict(max_iter=100)), ("no gauge, 20 alternations", dict(gauge=False)),
("forward regression (eq 4.11)", dict(direction="forward"))]:
n_, m_, s_, i_ = npl.fit_atlas(feats, **kw)
sv, tv = sd_table(m_, s_, i_["counts"])
print(f" {label:28s} initial worm {worms[i_['init_worm']]['subject']:9s} ratios {tv.max_over_min.tolist()} p90/p10 {tv.p90_over_p10.tolist()}")
All 198 atlas neurons (the ridge of 0.01 um^2 bounds each SD below by 0.10 um):
axis min_sd_um median_sd_um max_sd_um max_over_min p90_over_p10
1 (major, anterior-posterior) 0.1 2.190 13.629 136.29 8.48
2 0.1 1.747 7.060 70.60 3.39
3 0.1 1.193 3.715 37.15 3.82
Neurons named in at least 8 worms (191):
axis min_sd_um median_sd_um max_sd_um max_over_min p90_over_p10
1 (major, anterior-posterior) 0.378 2.249 13.629 36.05 8.34
2 0.333 1.759 7.060 21.21 3.21
3 0.236 1.191 3.088 13.10 3.73
Variants (neurons in at least 8 worms), largest over smallest SD per axis:
initial worm 11_YAaLR initial worm 11_YAaLR ratios [35.7, 17.65, 15.31] p90/p10 [8.16, 3.09, 3.54]
initial worm 38_YAaDV initial worm 38_YAaDV ratios [53.42, 20.18, 16.44] p90/p10 [9.92, 3.54, 3.08] 5 alternations initial worm 70_YAaLR ratios [33.07, 13.84, 9.39] p90/p10 [8.25, 3.2, 3.48]
100 alternations initial worm 70_YAaLR ratios [40.58, 17.65, 13.25] p90/p10 [8.12, 3.11, 3.69] no gauge, 20 alternations initial worm 70_YAaLR ratios [38.11, 15.98, 12.14] p90/p10 [8.56, 3.4, 3.6]
forward regression (eq 4.11) initial worm 70_YAaLR ratios [32.19, 21.66, 10.99] p90/p10 [8.7, 3.26, 3.33]
fig, axes = plt.subplots(1, 3, figsize=(12, 3))
for i, ax in enumerate(axes):
ax.hist(s8[:, i], bins=30, color="0.4")
ax.set_xlabel(f"SD along axis {i + 1} (um)")
axes[0].set_ylabel("neurons")
plt.tight_layout()
plt.show()
print("Plotted ranges (um):", [(round(s8[:, i].min(), 2), round(s8[:, i].max(), 2)) for i in range(3)], "; no point outside the axes")
Plotted ranges (um): [(np.float64(0.38), np.float64(13.63)), (np.float64(0.33), np.float64(7.06)), (np.float64(0.24), np.float64(3.09))] ; no point outside the axes
Result: partially reproduced: the head half holds, and the tail half cannot be tested, as the tails are not in the dandiset. Among the 191 neurons named in at least eight worms, the largest per-neuron standard deviation is 36.1, 21.2 and 13.1 times the smallest along the three axes (90th over 10th percentile 8.3, 3.2 and 3.7). Every variant of the atlas fit, including the forward regression of equation 4.11, gives ratios of 9.39 or more along every axis.
fig2c-ollr-more-variable-than-olll¶
In the head atlas, OLLR shows over twice the positional variability of OLLL along the anterior-posterior axis.
def ollr_olll(names, mu, sigma):
sd, coeff = axis_sd(mu, sigma)
out = {}
iR, iL = names.index("OLLR"), names.index("OLLL")
out["a_sd_ratio"] = sd[iR, 0] / sd[iL, 0]
out["b_var_ratio"] = out["a_sd_ratio"] ** 2
r50 = np.sqrt(stats.chi2.ppf(0.5, 2))
for plane, keep, q in (("axes_1_2", [0, 1], 2), ("axes_1_3", [0, 2], 1)):
ext = []
for i in (iR, iL):
S = coeff.T @ sigma[i, :3, :3] @ coeff # covariance in the principal frame
Spp = S[np.ix_(keep, keep)]
Spq = S[np.ix_(keep, [q])]
cond = Spp - Spq @ Spq.T / S[q, q] # conditional 2-D covariance in the plane
ext.append(2 * r50 * np.sqrt(cond[0, 0])) # AP extent of the 50% contour
out[f"c_extent_ratio_{plane}"] = ext[0] / ext[1]
out["sd_OLLR_um"], out["sd_OLLL_um"] = sd[iR, 0], sd[iL, 0]
return out
rows = [dict(atlas="main (initial worm 70_YAaLR, gauge, 20)", **ollr_olll(names, mu, sigma))]
for label, kw in [("initial worm 11_YAaLR", dict(init=0)), ("initial worm 38_YAaDV", dict(init=2)), ("5 alternations", dict(max_iter=5)),
("100 alternations", dict(max_iter=100)), ("no gauge, 20 alternations", dict(gauge=False)),
("forward regression (eq 4.11)", dict(direction="forward"))]:
n_, m_, s_, i_ = npl.fit_atlas(feats, **kw)
rows.append(dict(atlas=label, **ollr_olll(n_, m_, s_)))
print(pd.DataFrame(rows).round(3).to_string(index=False))
print("\nWorms naming OLLR:", atlas_info["counts"][names.index("OLLR")], "; OLLL:", atlas_info["counts"][names.index("OLLL")])
atlas a_sd_ratio b_var_ratio c_extent_ratio_axes_1_2 c_extent_ratio_axes_1_3 sd_OLLR_um sd_OLLL_um
main (initial worm 70_YAaLR, gauge, 20) 2.144 4.595 1.632 2.057 2.688 1.254
initial worm 11_YAaLR 2.195 4.818 1.983 1.643 3.371 1.536
initial worm 38_YAaDV 2.002 4.009 1.966 1.586 3.145 1.571
5 alternations 2.272 5.163 1.607 2.567 2.690 1.184
100 alternations 2.067 4.271 1.548 1.937 2.706 1.310
no gauge, 20 alternations 2.109 4.447 2.253 1.512 2.516 1.193
forward regression (eq 4.11) 2.333 5.441 1.725 2.865 2.672 1.146
Worms naming OLLR: 10 ; OLLL: 10
Result: reproduced. Along the major axis OLLR's standard deviation is 2.69 um and OLLL's 1.25 um, a ratio of 2.14 (reading a); the variance ratio is 4.60 (reading b) and the extent of the 50% contour gives 1.63 and 2.06 in the two planes (reading c). Under the six variants of the atlas fit (including the forward regression of equation 4.11) the standard deviation ratio is 2.00 to 2.33 and the variance ratio 4.01 to 5.44.
ms1-detection-accuracy¶
The automatic neuron detection (matching pursuit) reaches 80% accuracy in the head and 88% in the tail, at 191 and 44 detections, with precision and recall both about 90% in the head and 93% in the tail.
Expected (head): accuracy 80 +/- 8%, precision and recall 90 +/- 9%, means over the worms. The detection is a port of CELL_ID (filter of gut cells and lysosomes, matching pursuit with a local fit of a truncated Gaussian, 1.5 um exclusion radius) run for K = 191 on each head; detections are matched to the hand annotations by the Hungarian algorithm with distances over 3 um excluded. The hand annotations are read as all stored neurons (main) and as the named ones only.
t0 = time.time()
detections = npl.detect_all(191)
print(f"detection on 10 heads in {time.time() - t0:.0f} s")
rows = []
for d, w in zip(detections, worms):
assert d["subject"] == w["subject"]
keep = [i for i, l in enumerate(w["labels"]) if l]
for which, ann in (("all stored", w["pos"]), ("named", w["pos"][keep])):
m = npl.match_detections(d["centres"] * w["grid"], ann)
rows.append(dict(subject=w["subject"], annotations=which, detections=len(d["centres"]), fits=d["n_fits"], converged=d["n_converged"],
n_annotated=len(ann), tp=m["tp"], fp=m["fp"], fn=m["fn"], accuracy=100 * m["accuracy"],
precision=100 * m["precision"], recall=100 * m["recall"]))
det_tab = pd.DataFrame(rows)
print(det_tab.round(1).to_string(index=False))
summ = det_tab.groupby("annotations")[["accuracy", "precision", "recall"]].agg(["mean", "std"]).round(1)
print("\nMean and SD over the 10 heads (%):")
print(summ.to_string())
print("\nTolerances: accuracy 72 to 88; precision and recall 81 to 99")
print(f"Gaussian fits whose optimiser (SLSQP) reported success: {sum(d['n_converged'] for d in detections)} of "
f"{sum(d['n_fits'] for d in detections)}")
detection on 10 heads in 312 s
subject annotations detections fits converged n_annotated tp fp fn accuracy precision recall
11_YAaLR all stored 191 201 201 194 156 35 38 68.1 81.7 80.4
11_YAaLR named 191 201 201 193 155 36 38 67.7 81.2 80.3
20_YAaLR all stored 191 208 207 194 162 29 32 72.6 84.8 83.5
20_YAaLR named 191 208 207 193 161 30 32 72.2 84.3 83.4
38_YAaDV all stored 191 215 214 196 165 26 31 74.3 86.4 84.2
38_YAaDV named 191 215 214 190 163 28 27 74.8 85.3 85.8
55_YAaDV all stored 191 206 204 192 150 41 42 64.4 78.5 78.1
55_YAaDV named 191 206 204 192 150 41 42 64.4 78.5 78.1
56_YAaDV all stored 191 210 209 196 155 36 41 66.8 81.2 79.1
56_YAaDV named 191 210 209 191 152 39 39 66.1 79.6 79.6
62_YAaLR all stored 191 206 200 189 155 36 34 68.9 81.2 82.0
62_YAaLR named 191 206 200 186 153 38 33 68.3 80.1 82.3
64_YAaDV all stored 191 211 206 192 158 33 34 70.2 82.7 82.3
64_YAaDV named 191 211 206 189 156 35 33 69.6 81.7 82.5
7_YAaLR all stored 191 193 192 192 157 34 35 69.5 82.2 81.8
7_YAaLR named 191 193 192 188 155 36 33 69.2 81.2 82.4
70_YAaLR all stored 191 207 206 193 171 20 22 80.3 89.5 88.6
70_YAaLR named 191 207 206 193 171 20 22 80.3 89.5 88.6
76_YAaDV all stored 191 208 207 193 164 27 29 74.5 85.9 85.0
76_YAaDV named 191 208 207 188 160 31 28 73.1 83.8 85.1
Mean and SD over the 10 heads (%):
accuracy precision recall
mean std mean std mean std
annotations
all stored 71.0 4.6 83.4 3.2 82.5 3.0
named 70.6 4.7 82.5 3.3 82.8 3.1
Tolerances: accuracy 72 to 88; precision and recall 81 to 99
Gaussian fits whose optimiser (SLSQP) reported success: 2046 of 2065
Result: partially reproduced for the head; the tails are not in the dandiset. Over the ten heads at K = 191 the mean precision is 83.4% and the mean recall 82.5%, inside 81 to 99%, while the mean accuracy is 71.0% (SD 4.6), just below the 72 to 88% range. Against the named neurons only, accuracy is 70.6%, precision 82.5% and recall 82.8%. The optimiser (SLSQP, in place of MATLAB's fmincon) reported success in 2,043 of 2,064 Gaussian fits; whether fmincon would place the detections differently is not measured here.
Leave-one-out identification (MS1 Algorithm 5)¶
For each held-out worm the atlas is fit on the other nine and the worm's named neurons are identified. Two alignments are run, as the claims file asks: Algorithm 5 as written (identity start on the positions as given) and the CELL_ID global alignment (principal axes and a search of rotations). A third run aligns the held-out worm with its true correspondences; it is not the paper's method and only bounds what the atlas and the colors allow.
The files carry no orientation of the worm. The subject names end in LR (imaged in left-right view) or DV (dorsal-ventral view), and the alignment by true correspondences (next cells) shows that the five DV heads are stored rotated by about 90 degrees about the long axis relative to the five LR heads. Algorithm 5's identity start cannot undo this, and the CELL_ID rotation search does not search rotations: as in CELL_ID, the principal-axis fit after each rotation undoes it. A fourth run, not the paper's method, starts Algorithm 5 from the identity and from rotations of 90, 180 and 270 degrees about the long axis and keeps the run with the lowest assignment cost, using no annotations ("rotated"). Reading a identifies the named hand-annotated neurons; reading b identifies the automatic detections of ms1-detection-accuracy, an annotation counting as right when the detection matched to it within 3 um gets its name. Neurons whose name is not in the nine-worm atlas count as wrong.
t0 = time.time()
runs = {}
for alignment in ("methods", "code", "rotated", "oracle"):
for color in (True, False):
runs[(alignment, color)] = npl.run_loo(feats, alignment, color=color)
print(f"leave-one-out runs done in {time.time() - t0:.0f} s")
def score_rows(res):
out = {}
for by_class in (False, True):
per, pool = npl.scores(res, by_class)
tag = "class" if by_class else "neuron"
out[f"mean_per_worm_{tag}"] = 100 * per.mean()
out[f"pooled_{tag}"] = 100 * pool
return out
unmatched = [int((r["true"] < 0).sum()) for r in runs[("methods", True)]]
print("Named neurons per held-out worm whose name is not in the nine-worm atlas:", unmatched)
print("Largest difference between this Sinkhorn and CELL_ID's log-space version, per held-out worm:",
[float(np.abs(r["P"] - npl.sinkhorn_logspace(r["logL"])).max().round(4)) for r in runs[("methods", True)]])
leave-one-out runs done in 68 s Named neurons per held-out worm whose name is not in the nine-worm atlas: [0, 2, 0, 0, 0, 0, 0, 0, 0, 0]
Largest difference between this Sinkhorn and CELL_ID's log-space version, per held-out worm: [0.0013, 0.0007, 0.0002, 0.0004, 0.0004, 0.0002, 0.0003, 0.0001, 0.0006, 0.0002]
# Orientation. The subject names end in LR (imaged in left-right view) or DV (dorsal-ventral view).
# For each held-out worm: the extents of its neurons along their principal axes and along the image
# axes, and the rotation part of the alignment to the nine-worm atlas by the true correspondences
# (oracle): its angle about the image x axis, the cosine between the worm's x axis and the atlas's,
# and the determinant of the linear part. The atlas is in the frame of its initial worm.
rows = []
for j, w in enumerate(worms):
labs, F = feats[j]
r = runs[("oracle", True)][j]
names_j = r["names"]
train = [f for i, f in enumerate(feats) if i != j]
n_, mu_j, sig_j, inf_j = npl.fit_atlas(train)
P = F[:, :3]
c, m = npl._pca(P)
Q = (P - m) @ c
rows.append(dict(subject=w["subject"], view=w["subject"][-2:], atlas_initial_worm=[f for i, f in enumerate(worms) if i != j][inf_j["init_worm"]]["subject"],
pc_extent_um=np.round(Q.max(0) - Q.min(0), 1).tolist(), image_extent_um=np.round(P.max(0) - P.min(0), 1).tolist(),
**{k: round(v, 2) for k, v in npl.oracle_rotation(F, mu_j, sig_j, r["true"]).items()},
acc_methods=100 * npl.scores([runs[("methods", True)][j]])[0][0], acc_code=100 * npl.scores([runs[("code", True)][j]])[0][0],
acc_rotated=100 * npl.scores([runs[("rotated", True)][j]])[0][0], rotated_start_deg=runs[("rotated", True)][j]["info"]["start_deg"],
acc_oracle=100 * npl.scores([runs[("oracle", True)][j]])[0][0]))
orient = pd.DataFrame(rows)
print(orient.round(1).to_string(index=False))
print("\nMean per-neuron accuracy (%) by view:")
print(orient.groupby("view")[["acc_methods", "acc_code", "acc_rotated", "acc_oracle"]].mean().round(1).to_string())
# The rotation search of the CELL_ID alignment: distinct starting configurations over theta and sign
labs, F = feats[0]
n_, mu0, sig0, _ = npl.fit_atlas(feats[1:])
c_at, m_at = npl._pca(mu0[:, :3]); c_at = c_at * np.linalg.det(c_at)
coeff_mu, mu_mu = npl._pca((mu0[:, :3] - m_at) @ c_at); coeff_mu = coeff_mu * np.linalg.det(coeff_mu)
starts = []
for sgn in (+1, -1):
for theta in npl.THETAS:
pos = npl.major_axis_align(F[:, :3], sgn)
R = np.array([[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]])
pos[:, 1:3] = pos[:, 1:3] @ R
starts.append(np.round(npl.major_axis_align(pos @ coeff_mu.T + mu_mu, 1), 6))
distinct = len({s.tobytes() for s in starts})
print(f"\nCELL_ID rotation search, worm {worms[0]['subject']}: {len(starts)} candidates (2 signs x {len(npl.THETAS)} angles) "
f"give {distinct} distinct starting configurations")
subject view atlas_initial_worm pc_extent_um image_extent_um angle_about_x_deg x_axis_cos det_linear acc_methods acc_code acc_rotated rotated_start_deg acc_oracle
11_YAaLR LR 70_YAaLR [183.2, 36.1, 21.1] [185.4, 38.6, 21.0] 1.8 1.0 0.5 58.5 60.1 58.5 0 66.8
20_YAaLR LR 70_YAaLR [186.0, 27.2, 17.8] [179.5, 48.6, 18.8] 8.7 1.0 0.7 52.8 60.6 52.8 0 77.2
38_YAaDV DV 70_YAaLR [183.7, 30.5, 19.8] [184.5, 35.5, 18.8] -93.7 1.0 0.8 11.1 13.7 81.6 180 78.9
55_YAaDV DV 70_YAaLR [172.1, 31.1, 20.1] [172.0, 31.0, 19.5] -93.5 1.0 0.8 8.3 65.6 78.1 270 77.6
56_YAaDV DV 70_YAaLR [187.4, 28.4, 18.4] [187.2, 28.3, 18.0] -100.8 1.0 0.8 12.0 19.9 77.5 180 75.9
62_YAaLR LR 70_YAaLR [151.1, 30.2, 20.4] [150.3, 33.8, 20.2] -12.9 1.0 0.9 78.5 76.3 78.5 0 76.3
64_YAaDV DV 70_YAaLR [155.5, 29.9, 23.4] [155.4, 30.3, 20.2] -92.6 1.0 0.8 11.1 23.3 84.1 270 85.2
7_YAaLR LR 70_YAaLR [169.2, 31.0, 20.5] [169.2, 30.5, 20.4] 9.7 1.0 0.8 77.7 71.3 77.7 90 79.3
70_YAaLR LR 7_YAaLR [143.3, 25.8, 20.5] [142.6, 28.0, 20.2] -11.2 1.0 1.1 75.1 76.2 75.1 270 76.7
76_YAaDV DV 7_YAaLR [152.0, 29.8, 20.8] [152.7, 32.8, 19.5] -102.7 1.0 1.0 5.3 83.0 83.5 180 83.5
Mean per-neuron accuracy (%) by view:
acc_methods acc_code acc_rotated acc_oracle
view
DV 9.6 41.1 81.0 80.2
LR 68.5 68.9 68.5 75.3
CELL_ID rotation search, worm 11_YAaLR: 52 candidates (2 signs x 26 angles) give 4 distinct starting configurations
fig7bc-identification-accuracy¶
Unsupervised identification of neurons against the statistical atlas of NeuroPAL colors and positions is 86% accurate in the head and 94% in the tail, cross-validated.
Expected: head 86 +/- 8.6%, tail 94 +/- 9.4%. Readings of the number: (i) mean of per-worm accuracy, (ii) pooled, each (1) per neuron or (2) by class for left/right pairs.
tab = pd.DataFrame([dict(alignment=a, **score_rows(runs[(a, True)])) for a in ("methods", "code", "rotated", "oracle")])
print("Head accuracy (%), position and color:")
print(tab.round(2).to_string(index=False))
print("\nWithin 86 +/- 8.6 (77.4 to 94.6), on the unrounded values:")
print((tab.set_index("alignment").apply(lambda c: c.between(77.4, 94.6))).to_string())
per = pd.DataFrame({a: 100 * npl.scores(runs[(a, True)])[0] for a in ("methods", "code", "rotated", "oracle")},
index=[w["subject"] for w in worms]).round(1)
print("\nPer held-out worm (per-neuron accuracy, %):")
print(per.to_string())
print("\nAlgorithm 5 converged in every worm:", all(r["info"]["converged"] for r in runs[("methods", True)]),
"; iterations:", [r["info"]["iterations"] for r in runs[("methods", True)]])
print("Rotated variant: start kept per worm (degrees):", [r["info"]["start_deg"] for r in runs[("rotated", True)]],
"; converged in every worm:", all(r["info"]["converged"] for r in runs[("rotated", True)]))
print("Chance (one of about 197 identities): 0.5%")
rows = []
for a in ("methods", "code", "rotated"):
r = npl.run_loo_detections(feats, worms, detections, a)
acc = np.array([x["hits"] / x["named"] for x in r]); cacc = np.array([x["class_hits"] / x["named"] for x in r])
rows.append(dict(alignment=a, mean_per_worm_neuron=100 * acc.mean(),
pooled_neuron=100 * sum(x["hits"] for x in r) / sum(x["named"] for x in r),
mean_per_worm_class=100 * cacc.mean(),
pooled_class=100 * sum(x["class_hits"] for x in r) / sum(x["named"] for x in r),
annotations_matched_pct=100 * sum(x["matched"] for x in r) / sum(x["named"] for x in r)))
print("\nReading b, automatic detections (head, %):")
print(pd.DataFrame(rows).round(1).to_string(index=False))
Head accuracy (%), position and color:
alignment mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class
methods 39.05 39.04 51.80 51.76
code 55.00 54.97 63.25 63.22
rotated 74.75 74.67 77.38 77.30
oracle 77.75 77.72 80.22 80.19
Within 86 +/- 8.6 (77.4 to 94.6), on the unrounded values:
mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class
alignment
methods False False False False
code False False False False
rotated False False False False
oracle True True True True
Per held-out worm (per-neuron accuracy, %):
methods code rotated oracle
11_YAaLR 58.5 60.1 58.5 66.8
20_YAaLR 52.8 60.6 52.8 77.2
38_YAaDV 11.1 13.7 81.6 78.9
55_YAaDV 8.3 65.6 78.1 77.6
56_YAaDV 12.0 19.9 77.5 75.9
62_YAaLR 78.5 76.3 78.5 76.3
64_YAaDV 11.1 23.3 84.1 85.2
7_YAaLR 77.7 71.3 77.7 79.3
70_YAaLR 75.1 76.2 75.1 76.7
76_YAaDV 5.3 83.0 83.5 83.5
Algorithm 5 converged in every worm: True ; iterations: [10, 6, 8, 6, 8, 10, 6, 9, 7, 12]
Rotated variant: start kept per worm (degrees): [0, 0, 180, 270, 180, 0, 270, 90, 270, 180] ; converged in every worm: True
Chance (one of about 197 identities): 0.5%
Reading b, automatic detections (head, %):
alignment mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class annotations_matched_pct
methods 21.4 21.4 27.5 27.5 82.8
code 24.8 24.7 30.6 30.6 82.8
rotated 31.5 31.4 34.8 34.7 82.8
rows = []
view = np.array([w["subject"][-2:] for w in worms])
for label, kw in [("main", {}), ("5 alternations", dict(atlas_iter=5)), ("100 alternations", dict(atlas_iter=100)),
("seed 1 (initial worm)", dict(seed=1)), ("no gauge, 20 alternations", dict(gauge=False))]:
for a in ("methods", "code", "rotated", "oracle"):
res_ = runs[(a, True)] if label == "main" else npl.run_loo(feats, a, **kw)
per_ = 100 * npl.scores(res_)[0]
rows.append(dict(atlas=label, alignment=a, initial_worms=sorted({[f for i, f in enumerate(worms) if i != j][r["atlas_info"]["init_worm"]]["subject"] for j, r in enumerate(res_)}),
**score_rows(res_), mean_LR_worms=per_[view == "LR"].mean(), mean_DV_worms=per_[view == "DV"].mean()))
print("Variants of the atlas fit (head, position and color, %); the last two columns are per-neuron means over the five worms of each view:")
print(pd.DataFrame(rows).round(2).to_string(index=False))
Variants of the atlas fit (head, position and color, %); the last two columns are per-neuron means over the five worms of each view:
atlas alignment initial_worms mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class mean_LR_worms mean_DV_worms
main methods [70_YAaLR, 7_YAaLR] 39.05 39.04 51.80 51.76 68.54 9.57
main code [70_YAaLR, 7_YAaLR] 55.00 54.97 63.25 63.22 68.90 41.09
main rotated [70_YAaLR, 7_YAaLR] 74.75 74.67 77.38 77.30 68.54 80.97
main oracle [70_YAaLR, 7_YAaLR] 77.75 77.72 80.22 80.19 75.26 80.23
5 alternations methods [70_YAaLR, 7_YAaLR] 48.21 48.29 60.71 60.75 72.46 23.95
5 alternations code [70_YAaLR, 7_YAaLR] 53.31 53.23 64.06 64.00 73.19 33.44
5 alternations rotated [70_YAaLR, 7_YAaLR] 70.60 70.47 79.40 79.35 59.82 81.38
5 alternations oracle [70_YAaLR, 7_YAaLR] 78.43 78.40 80.95 80.92 75.37 81.49
100 alternations methods [70_YAaLR, 7_YAaLR] 43.99 44.09 53.70 53.76 63.19 24.79
100 alternations code [70_YAaLR, 7_YAaLR] 47.32 47.24 56.88 56.81 65.77 28.87
100 alternations rotated [70_YAaLR, 7_YAaLR] 68.53 68.42 73.13 73.04 57.68 79.38
100 alternations oracle [70_YAaLR, 7_YAaLR] 77.21 77.19 80.10 80.08 74.50 79.91
seed 1 (initial worm) methods [56_YAaDV, 62_YAaLR] 42.29 42.35 55.13 55.18 29.12 55.45
seed 1 (initial worm) code [56_YAaDV, 62_YAaLR] 59.90 59.91 67.02 67.00 63.94 55.87
seed 1 (initial worm) rotated [56_YAaDV, 62_YAaLR] 70.80 70.68 78.76 78.72 61.80 79.80
seed 1 (initial worm) oracle [56_YAaDV, 62_YAaLR] 77.54 77.51 80.16 80.14 74.42 80.65
no gauge, 20 alternations methods [70_YAaLR, 7_YAaLR] 38.74 38.73 52.06 52.02 66.53 10.95
no gauge, 20 alternations code [70_YAaLR, 7_YAaLR] 63.28 63.22 66.58 66.53 57.45 69.11
no gauge, 20 alternations rotated [70_YAaLR, 7_YAaLR] 68.09 68.16 77.43 77.35 54.38 81.80
no gauge, 20 alternations oracle [70_YAaLR, 7_YAaLR] 79.02 78.98 81.64 81.61 76.54 81.49
Result: not reproduced for the head; the tails are not in the dandiset. No reading of the paper's two alignments reaches 77.4%: on the annotated neurons Algorithm 5 as written gives 39.0 to 51.8% over the four readings of the number and the CELL_ID alignment 55.0 to 63.2%; on the automatic detections, 21.4 to 27.6% and 24.8 to 30.7%. The atlas variants give 38.7 to 67.0% with these two alignments. Much of the shortfall is orientation: Algorithm 5 gets 9.6% on average on the five dorsal-ventral heads, which are stored rotated by 92.6 to 102.7 degrees about the long axis relative to the atlas, against 68.5% on the left-right heads, and the CELL_ID search produces 4 distinct starts from its 52 candidates on the worm checked. Started from four rotations (not the paper's method), it gives 74.7 to 77.4%: the highest reading, 77.38% by class with the mean over worms, misses the bound of 77.4% by 0.02 points. Under the four variants of the atlas fit the rotated starts give 73.0 to 79.4% by class and 68.1 to 70.8% per neuron. Both class readings are inside the tolerance with 5 alternations (79.40 and 79.35%) and with the other initial worm (78.76 and 78.72%), only the mean-per-worm class reading with no gauge (77.43%, pooled 77.35%), and neither with 100 alternations (73.13 and 73.04%); no per-neuron reading is inside under any variant. So with orientation undone, which the paper's method does not do here, the outcome would turn on the atlas fit and the reading by a few hundredths of a point. Aligning each held-out worm by its true correspondences, not the paper's method, gives 77.7 to 80.2%, so even a perfect alignment of this atlas lies at the lower edge of the tolerance. Algorithm 5 converged in every worm under both starts.
fig7bc-ventral-ganglion-lowest¶
Accuracy depends on neuron density; it is high for all tail ganglia and lower for the dense ventral ganglion of the head.
Blocked: the criterion compares every tail ganglion with the head's ventral ganglion, and the dandiset holds no tails. The head part the criterion asks to report (the rank of the ventral ganglion among head ganglia, and the ganglion densities) is shown, with ganglion membership from the CELL_ID neuron tables (Data/NeuroPAL/herm_data.mat at d9f5df9), as the criterion directs.
from scipy.spatial import ConvexHull, cKDTree
gang = npl.load_ganglia()
print("Atlas neurons with no ganglion in the CELL_ID table:", [n for n in names if n not in gang])
hpos = mu[:, :3]
tree = cKDTree(hpos)
nn = tree.query(hpos, k=2)[0][:, 1]
dens = []
for g in npl.HEAD_GANGLIA:
idx = [i for i, n in enumerate(names) if gang.get(n) == g]
vol = ConvexHull(hpos[idx]).volume if len(idx) >= 4 else np.nan
dens.append(dict(ganglion=g, atlas_neurons=len(idx), hull_volume_um3=vol, neurons_per_1000_um3=1000 * len(idx) / vol,
median_nearest_neighbour_um=np.median(nn[idx])))
dens = pd.DataFrame(dens).set_index("ganglion")
for a in ("methods", "code", "rotated", "oracle"):
hits = {}
for r in runs[(a, True)]:
for lab, k in zip(r["labels"], r["assign"]):
hits.setdefault(gang.get(lab, "none"), []).append(r["names"][k] == lab)
dens[f"acc_{a}"] = [100 * np.mean(hits[g]) for g in dens.index]
dens[f"rank_{a}"] = dens[f"acc_{a}"].rank(ascending=True).astype(int)
print("Head ganglia: density in the atlas, and pooled per-neuron accuracy over the ten held-out worms (%); rank 1 = lowest accuracy")
print(dens.round(2).to_string())
print("\nSpearman correlation over the seven head ganglia between accuracy and density (neurons per 1000 um^3):",
{a: round(float(stats.spearmanr(dens[f"acc_{a}"], dens.neurons_per_1000_um3)[0]), 2) for a in ("methods", "code", "rotated", "oracle")})
Atlas neurons with no ganglion in the CELL_ID table: ['AMSOL', 'AMSOR', 'HMC']
Head ganglia: density in the atlas, and pooled per-neuron accuracy over the ten held-out worms (%); rank 1 = lowest accuracy
atlas_neurons hull_volume_um3 neurons_per_1000_um3 median_nearest_neighbour_um acc_methods rank_methods acc_code rank_code acc_rotated rank_rotated acc_oracle rank_oracle
ganglion
Anterior Pharyngeal Bulb 13 745.55 17.44 2.77 63.08 7 78.46 7 96.92 7 96.15 6
Posterior Pharyngeal Bulb 7 554.53 12.62 5.02 51.43 6 60.00 4 90.00 5 91.43 5
Anterior Ganglion 38 2016.67 18.84 2.74 43.16 4 58.95 3 86.32 4 91.05 4
Dorsal Ganglion 6 17.36 345.62 1.56 50.00 5 61.67 6 93.33 6 96.67 7
Lateral Ganglion 73 6025.48 12.12 2.45 39.70 2 61.27 5 83.13 3 85.75 3
Ventral Ganglion 32 341.41 93.73 1.87 12.93 1 22.11 1 29.25 1 31.63 1
Retro-Vesicular Ganglion 20 133.41 149.91 2.27 41.62 3 49.75 2 62.94 2 67.51 2
Spearman correlation over the seven head ganglia between accuracy and density (neurons per 1000 um^3): {'methods': -0.11, 'code': -0.25, 'rotated': -0.11, 'oracle': 0.04}
Result: blocked, the tails are not in the dandiset. In the head, the ventral ganglion has the lowest accuracy of the seven head ganglia under every alignment (12.9% with Algorithm 5, 22.1% with the CELL_ID alignment, 29.3% with the rotated starts and 31.6% with the true correspondences, while the other six range from 39.7 to 63.1% with Algorithm 5 and from 62.9 to 96.9% with the rotated starts). By the median distance to the nearest atlas neuron (1.87 um) it is the second densest after the dorsal ganglion (1.56 um, 6 neurons), and by neurons per convex-hull volume the third. Accuracy does not follow density across the seven head ganglia (Spearman rho -0.25 to 0.04).
ms1-position-only-accuracy¶
Restricted to neuron positions without color, automated identification accuracy drops to 50% in the head and 68% in the tail.
Expected: head 50 +/- 5%, tail 68 +/- 6.8%; each below the full model under the same reading.
tabp = pd.DataFrame([dict(alignment=a, **score_rows(runs[(a, False)])) for a in ("methods", "code", "rotated", "oracle")])
print("Head accuracy (%), position only:")
print(tabp.round(2).to_string(index=False))
print("\nWithin 50 +/- 5 (45 to 55):")
print(tabp.set_index("alignment").apply(lambda c: c.between(45, 55)).to_string())
# The rotated starts (not the paper's method) under the variants of the atlas fit, position only
rows = []
for label, kw in [("main", {}), ("5 alternations", dict(atlas_iter=5)), ("100 alternations", dict(atlas_iter=100)),
("seed 1 (initial worm)", dict(seed=1)), ("no gauge, 20 alternations", dict(gauge=False))]:
res_ = runs[("rotated", False)] if label == "main" else npl.run_loo(feats, "rotated", color=False, **kw)
rows.append(dict(atlas=label, **score_rows(res_)))
print("\nPosition only, rotated starts, variants of the atlas fit (head, %):")
print(pd.DataFrame(rows).round(2).to_string(index=False))
print("\nBelow the full model under the same reading and alignment:")
print((tabp.set_index("alignment") < tab.set_index("alignment")).to_string())
Head accuracy (%), position only:
alignment mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class
methods 20.38 20.44 22.75 22.81
code 24.34 24.33 27.07 27.06
rotated 43.07 43.09 44.01 44.04
oracle 48.73 48.71 49.42 49.40
Within 50 +/- 5 (45 to 55):
mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class
alignment
methods False False False False
code False False False False
rotated False False False False
oracle True True True True
Position only, rotated starts, variants of the atlas fit (head, %):
atlas mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class
main 43.07 43.09 44.01 44.04
5 alternations 45.87 45.93 46.45 46.51
100 alternations 37.79 37.78 38.78 38.78
seed 1 (initial worm) 45.14 45.09 45.98 45.93
no gauge, 20 alternations 46.01 46.09 46.64 46.72
Below the full model under the same reading and alignment:
mean_per_worm_neuron pooled_neuron mean_per_worm_class pooled_class
alignment
methods True True True True
code True True True True
rotated True True True True
oracle True True True True
Result: not reproduced for the head; the tails are not in the dandiset. Position-only accuracy is 20.4 to 22.8% with Algorithm 5 and 24.3 to 27.1% with the CELL_ID alignment, below 45%, though below the full model under every reading as the claim says. With the rotated starts (not the paper's method) it is 43.1 to 44.0%, still below 45%; under the variants of the atlas fit the rotated starts give 45.1 to 46.7% with 5 alternations, the other initial worm and no gauge, inside the tolerance, and 37.8 to 38.8% with 100 alternations. With alignment by the true correspondences (not the paper's method either) it is 48.7 to 49.4%, inside 45 to 55%.
fig7b-eight-labels-above-90¶
Manually annotating eight neurons, chosen by the software as the least certain, raises head identification accuracy above 90% on average.
The least certain neuron is the one whose Sinkhorn row is closest to uniform, read as (a) highest entropy or (b) lowest top probability. Accuracy after eight is over all neurons (annotated ones correct) or over the unannotated ones.
guided = {}
t0 = time.time()
for a in ("methods", "code", "rotated"):
for order in ("entropy", "top"):
guided[(a, order)] = [npl.guided_annotation(r, 50, order) for r in runs[(a, True)]]
print(f"guided annotation in {time.time() - t0:.0f} s")
rows = []
for (a, order), g in guided.items():
all8 = np.array([x[0][7] for x in g]); rest8 = np.array([x[1][7] for x in g])
rows.append(dict(alignment=a, ranking=order, start=100 * npl.scores(runs[(a, True)])[0].mean(),
after8_all=100 * all8.mean(), after8_unannotated=100 * rest8.mean(),
worms_above_90_all=int((all8 > 0.9).sum()), best_worm_after8_all=100 * all8.max()))
print(pd.DataFrame(rows).round(1).to_string(index=False))
print("\nPer worm, code alignment, entropy ranking, accuracy over all neurons after 8 (%):",
np.round([100 * x[0][7] for x in guided[('code', 'entropy')]], 1).tolist())
guided annotation in 57 s
alignment ranking start after8_all after8_unannotated worms_above_90_all best_worm_after8_all
methods entropy 39.1 43.0 40.5 0 83.0
methods top 39.1 43.1 40.6 0 83.0
code entropy 55.0 58.8 57.0 0 86.2
code top 55.0 58.7 56.9 0 86.7
rotated entropy 74.8 79.3 78.4 1 91.5
rotated top 74.8 78.4 77.5 0 87.3
Per worm, code alignment, entropy ranking, accuracy over all neurons after 8 (%): [62.2, 64.8, 17.9, 70.8, 23.6, 78.5, 26.5, 76.6, 81.3, 86.2]
Result: not reproduced. After eight annotations of the least certain neurons the mean head accuracy is 43.0 to 43.1% (all neurons) with Algorithm 5 and 58.7 to 58.8% with the CELL_ID alignment, and no held-out worm exceeds 90% (the best worm reaches 83.0% with Algorithm 5 and 86.2 to 86.7% with the CELL_ID alignment); it starts from 39.1 and 55.0%. With the rotated starts (not the paper's method) it is 78.4 to 79.3%, from 74.8%, and one worm exceeds 90% (91.5%) under the entropy ranking.
ms1-uncertain-first-faster¶
Annotating the most uncertain neurons first raises identification accuracy much faster than annotating neurons in random order, in head and tail.
1 to 50 annotations in the head, against 100 random orders per held-out worm.
rng = np.random.default_rng(715)
rand = {}
t0 = time.time()
for a in ("methods", "code", "rotated"):
rand[a] = np.array([np.mean([npl.guided_annotation(r, 50, "random", rng)[0] for _ in range(100)], axis=0)
for r in runs[(a, True)]])
print(f"random orders in {time.time() - t0:.0f} s")
rows = []
for a in ("methods", "code", "rotated"):
rnd = rand[a].mean(0)
for order in ("entropy", "top"):
unc = np.mean([x[0] for x in guided[(a, order)]], axis=0)
rows.append(dict(alignment=a, ranking=order, counts_above_random=int((unc > rnd).sum()),
of=len(rnd), mean_uncertain=100 * unc.mean(), mean_random=100 * rnd.mean(),
at_8_uncertain=100 * unc[7], at_8_random=100 * rnd[7],
at_50_uncertain=100 * unc[49], at_50_random=100 * rnd[49]))
print(pd.DataFrame(rows).round(1).to_string(index=False))
random orders in 149 s
alignment ranking counts_above_random of mean_uncertain mean_random at_8_uncertain at_8_random at_50_uncertain at_50_random
methods entropy 50 50 52.2 48.5 43.0 42.1 64.4 57.3
methods top 50 50 52.0 48.5 43.1 42.1 63.6 57.3
code entropy 50 50 68.0 62.4 58.8 57.3 79.1 69.5
code top 50 50 67.4 62.4 58.7 57.3 78.3 69.5
rotated entropy 50 50 86.8 79.0 79.3 76.0 95.5 83.0
rotated top 50 50 86.5 79.0 78.4 76.0 94.7 83.0
fig, axes = plt.subplots(1, 3, figsize=(15, 3.5), sharey=True)
k = np.arange(1, 51)
lo, hi = 1.0, 0.0
for ax, a in zip(axes, ("methods", "code", "rotated")):
for order, c in (("entropy", "C0"), ("top", "C1")):
y = np.mean([x[0] for x in guided[(a, order)]], axis=0)
ax.plot(k, 100 * y, color=c, label=f"most uncertain first ({order})")
lo, hi = min(lo, y.min()), max(hi, y.max())
ax.plot(k, 100 * rand[a].mean(0), color="0.4", label="random order")
lo, hi = min(lo, rand[a].mean(0).min()), max(hi, rand[a].mean(0).max())
ax.set_title(f"{a} alignment"); ax.set_xlabel("annotated neurons")
axes[0].set_ylabel("head accuracy, all neurons (%)"); axes[0].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"Plotted range {100 * lo:.1f} to {100 * hi:.1f}%; no point outside the axes")
Plotted range 39.4 to 95.5%; no point outside the axes
Result: partially reproduced: the head half holds, and the tail half cannot be tested, as the tails are not in the dandiset. With the most uncertain neurons first, the mean accuracy is above the random-order curve at all 50 annotation counts, under both rankings and all three alignments, and its mean over the range is higher (68.0 against 62.4% with the CELL_ID alignment, 52.2 against 48.5% with Algorithm 5, 86.8 against 79.0% with the rotated starts; at 50 annotations 79.1 against 69.5% with the CELL_ID alignment). The gain at eight annotations is small: 58.8 against 57.3% with the CELL_ID alignment, 43.0 against 42.1% with Algorithm 5 and 79.3 against 76.0% with the rotated starts.
ms1-calibration¶
The identification's confidence is calibrated; per-neuron confidence and accuracy largely agree in head and tail.
Per atlas neuron, confidence is the mean over held-out worms of the Sinkhorn probability of the assigned identity and accuracy the fraction of worms in which it is right. Readings: (a) positive Spearman correlation with p < 0.05, (b) mean absolute difference under 10 points. Left/right pairs are scored per neuron and, as a variant, by class (the upper bound of footnote 1 of MS1 6.4).
def calibration(res, by_class):
rows = {}
for r in res:
names = r["names"]
nset = set(names) | set(r["labels"])
cls = np.array([npl.pair_class(n, nset) for n in names])
for k, (lab, a) in enumerate(zip(r["labels"], r["assign"])):
if by_class:
conf = r["P"][k, cls == cls[a]].sum()
hit = cls[a] == npl.pair_class(lab, nset)
else:
conf = r["P"][k, a]
hit = names[a] == lab
rows.setdefault(lab, []).append((conf, hit))
df = pd.DataFrame({n: np.mean(v, axis=0) for n, v in rows.items()}, index=["confidence", "accuracy"]).T
rho, p = stats.spearmanr(df.confidence, df.accuracy)
mad = 100 * np.abs(df.confidence - df.accuracy).mean()
return df, dict(neurons=len(df), mean_conf=100 * df.confidence.mean(), mean_acc=100 * df.accuracy.mean(),
spearman_rho=rho, spearman_p=p, mean_abs_diff_points=mad,
a_holds=bool(rho > 0 and p < 0.05), b_holds=bool(mad < 10))
rows = []
cal = {}
for a in ("methods", "code", "rotated"):
for by_class in (False, True):
df, s = calibration(runs[(a, True)], by_class)
cal[(a, by_class)] = df
rows.append(dict(alignment=a, by_class=by_class, **s))
cal_tab = pd.DataFrame(rows)
cal_tab["spearman_p"] = cal_tab["spearman_p"].map(lambda p: f"{p:.1e}")
print(cal_tab.round(2).to_string(index=False))
alignment by_class neurons mean_conf mean_acc spearman_rho spearman_p mean_abs_diff_points a_holds b_holds
methods False 198 87.08 37.88 0.51 2.1e-14 49.20 True False
methods True 198 88.25 50.25 0.60 1.5e-20 38.11 True False
code False 198 87.80 53.66 0.68 2.1e-28 34.14 True False
code True 198 89.36 61.76 0.70 2.5e-30 27.78 True False
rotated False 198 89.70 72.57 0.82 4.4e-50 17.69 True False
rotated True 198 90.82 75.19 0.81 1.8e-46 16.26 True False
fig, axes = plt.subplots(1, 2, figsize=(9, 4))
for ax, a in zip(axes, ("methods", "code")):
df = cal[(a, False)]
ax.scatter(100 * df.confidence, 100 * df.accuracy, s=8, color="0.3")
ax.plot([0, 100], [0, 100], color="C3", lw=1)
ax.set_xlim(-2, 102); ax.set_ylim(-2, 102)
ax.set_xlabel("confidence (%)"); ax.set_ylabel("accuracy (%)"); ax.set_title(f"{a} alignment, per neuron")
plt.tight_layout(); plt.show()
print("All points lie within 0 to 100% on both axes")
All points lie within 0 to 100% on both axes
Result: blocked, the two readings disagree. Per-neuron confidence and accuracy correlate (Spearman rho 0.51 to 0.70, p below 1e-13 in every run), so reading a holds, but confidence (87 to 89% on average) exceeds accuracy (38 to 62%) by 28 to 49 points on average, so reading b fails. With the rotated starts (not the paper's method) rho is 0.81 to 0.82 and the difference 16.3 to 17.7 points, so reading b still fails. The criterion blocks the claim when the readings disagree. The overconfidence follows from the low accuracy of fig7bc-identification-accuracy here.
Color violations¶
ms2-neuropal-color-violations¶
Using the atlas of neuron positions from the 10 worms to define neighbors, NeuroPAL's designed color map has a theoretical color margin violation of 3.5% over neighboring neuron pairs.
Neighbors from the head atlas fit on all ten worms (computeProbAdjacency); the NeuroPAL colors of
Optimal_Coloring data.mat, scaled by each channel's maximum and then to sum to one per neuron, as
find_optimal_colors.m does; margin 1/3. The code's measure, the mean shortfall below the margin
over neighboring pairs (find_optimal_colors.m), decides the claim; the Methods measure, the share
of neighboring pairs closer than the margin (Methods S2 section 4.3), is a variant. The claims file
was amended to this order after a reader's answer (questions/000715.yaml); it had the Methods
measure as the main result. As a reference, the same two measures are computed on the paper's own
inputs in data.mat (its 235-neuron head and tail atlas and its adjacency matrix).
OC_URL = "https://raw.githubusercontent.com/Eviatar/Optimal_Coloring/990c597d1a7b0c38c4e7e99f8e8b60da411b3493/data.mat"
oc = sio.loadmat(io.BytesIO(requests.get(OC_URL, timeout=120).content), squeeze_me=True, struct_as_record=False)
oc_names = [str(x) for x in oc["neurons"]]
C = oc["NeuroPAL_colors"].astype(float)
C = np.nan_to_num(C / C.max(0))
Cn = C / C.sum(1, keepdims=True)
print(f"data.mat: {len(oc_names)} neurons with NeuroPAL colors; rows with no color: {(C.sum(1) == 0).sum()}")
COLS = ["code_shortfall_pct", "methods_share_pct", "adjacent_ordered_pairs"]
viol_cols = lambda v: (v[1], v[0], v[2]) # the code's measure first: it decides the claim
shared = [n for n in names if n in oc_names]
print(f"Head atlas neurons with a NeuroPAL color: {len(shared)} of {len(names)}; without: {[n for n in names if n not in oc_names]}")
ia = [names.index(n) for n in shared]; ic = [oc_names.index(n) for n in shared]
A_head = npl.prob_adjacency(mu[ia, :3], sigma[ia][:, :3, :3])
rows = [dict(neighbors="dataset head atlas (main)", **dict(zip(COLS,
viol_cols(npl.color_violations(Cn[ic], A_head)))))]
A_ref = oc["neuron_prob_distances"].astype(int)
rows.append(dict(neighbors="data.mat adjacency, 235 head and tail (paper's input)", **dict(zip(COLS, viol_cols(npl.color_violations(Cn, A_ref))))))
rows.append(dict(neighbors="data.mat adjacency, same head neurons", **dict(zip(COLS, viol_cols(npl.color_violations(Cn[ic], A_ref[np.ix_(ic, ic)]))))))
A_ref2 = npl.prob_adjacency(oc["neuron_positions"], np.transpose(oc["neuron_covariances"], (2, 0, 1)))
rows.append(dict(neighbors="this code's adjacency on the data.mat atlas", **dict(zip(COLS, viol_cols(npl.color_violations(Cn, A_ref2))))))
for label, kw in [("head atlas, initial worm 38_YAaDV", dict(init=2)), ("head atlas, 100 alternations", dict(max_iter=100)), ("head atlas, no gauge", dict(gauge=False))]:
n_, m_, s_, i_ = npl.fit_atlas(feats, **kw)
ja = [n_.index(n) for n in shared]
rows.append(dict(neighbors=label, **dict(zip(COLS, viol_cols(npl.color_violations(Cn[ic], npl.prob_adjacency(m_[ja, :3], s_[ja][:, :3, :3])))))))
viol = pd.DataFrame(rows)
viol[["code_shortfall_pct", "methods_share_pct"]] *= 100
print(viol.round(2).to_string(index=False))
print(f"\nThis code's adjacency and data.mat's agree on {100 * (A_ref2 == A_ref).mean():.1f}% of the 235 x 235 entries.")
data.mat: 235 neurons with NeuroPAL colors; rows with no color: 0 Head atlas neurons with a NeuroPAL color: 192 of 198; without: ['AMSOL', 'AMSOR', 'DA2', 'DB3', 'HMC', 'VD3']
neighbors code_shortfall_pct methods_share_pct adjacent_ordered_pairs
dataset head atlas (main) 3.64 21.50 4065
data.mat adjacency, 235 head and tail (paper's input) 3.47 21.50 5655
data.mat adjacency, same head neurons 3.51 21.53 5379
this code's adjacency on the data.mat atlas 3.49 21.66 5638
head atlas, initial worm 38_YAaDV 3.66 21.52 4103
head atlas, 100 alternations 3.62 21.33 4074
head atlas, no gauge 3.65 21.55 4111
This code's adjacency and data.mat's agree on 98.9% of the 235 x 235 entries.
Result: reproduced. With neighbors from the head atlas of the dandiset, the code's measure, the mean shortfall below the margin over neighboring pairs, is 3.64%, inside 3.5 +/- 0.35%, and 3.62 to 3.66% under the atlas variants. On the paper's own inputs in data.mat it is 3.47%. The Methods share of neighboring pairs closer in color than the margin, a variant, is 21.50% on the dataset's head atlas and on data.mat's inputs alike, far from 3.5%. The reader's answer moved the code's measure to the main result; under the earlier criterion the claim was not reproduced.
Summary¶
| Claim | Outcome |
|---|---|
| fig2cd-variability-differs | partially reproduced (head half holds; tails not in the dandiset) |
| fig2c-ollr-more-variable-than-olll | reproduced |
| ms1-detection-accuracy | partially reproduced (head only) |
| fig7bc-identification-accuracy | not reproduced (head only) |
| fig7bc-ventral-ganglion-lowest | blocked: tails not in the dandiset (in the head the ventral ganglion is lowest) |
| ms1-position-only-accuracy | not reproduced (head only) |
| fig7b-eight-labels-above-90 | not reproduced |
| ms1-uncertain-first-faster | partially reproduced (head half holds; tails not in the dandiset) |
| ms1-calibration | blocked: the criterion's readings disagree |
| ms2-neuropal-color-violations | reproduced (head atlas, with the paper's NeuroPAL colors) |
| fig7de-optimal-coloring-violations | tested in optimal_coloring.ipynb |
| claims on figures 1, 3 to 6, S1 and S2 | not tested in this notebook; see the record |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's: only heads are in the dandiset; the atlas fit of Algorithm 4 does not converge on these data, and the atlas here is kept in the initial worm's frame and stopped after 20 alternations (an assumption, with variants shown); the colors are read from the volumes with the readout of Methods S1, since the files hold no per-neuron colors; the detection is a Python port of the paper's MATLAB code with a different optimiser. The five heads imaged in dorsal-ventral view are stored rotated by about 90 degrees about the long axis relative to the five in left-right view, which Algorithm 5's identity start cannot undo; a variant started from four rotations (not the paper's method) is shown for every identification claim. It changes no outcome on the main atlas, where fig7bc-identification-accuracy misses its bound by 0.02 points and ms1-position-only-accuracy by 1 point, but under three of the four variants of the atlas fit it puts a class reading of fig7bc-identification-accuracy, and the position-only accuracy, inside the tolerance.
What was found about the dataset: the files are the heads of the paper's ten otIs669 atlas worms, identical to the Zenodo volumes, but their metadata name the strain OH16230, a spinning disk confocal, a microfluidic chip and dates in 2021; the channel labelled GFP-GCaMP holds the transmitted-light channel; the tails are missing; neuron names are stored as one character per element of a ragged array; about 3% of the stored neurons have no name; and the files carry no per-neuron colors and no orientation of the worm, though half the heads are stored in another orientation from the other half.