import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import numpy as np
import pandas as pd
import scipy.stats as ss
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
from IPython.display import display
import scala

pd.set_option("display.width", 200, "display.max_columns", 40, "display.max_rows", 200)
T0 = time.time()
print("dandiset", scala.DANDISET, "version", scala.VERSION)

dandiset 000008 version 0.211014.0809

vmeta = scala.version_metadata()
assets = scala.list_assets()
assets_0808 = scala.list_assets("0.211014.0808")
same = set(zip(assets.path, assets["size"])) == set(zip(assets_0808.path, assets_0808["size"]))
print("assets in 0.211014.0809:", len(assets), "| in 0.211014.0808:", len(assets_0808),
      "| same paths and sizes:", same)
print("archive summary:", {k: vmeta["assetsSummary"][k] for k in ["numberOfFiles", "numberOfSubjects", "numberOfCells", "numberOfBytes"]})
t = time.time()
meta, feats, stim, binned, errors = scala.load_all(assets, workers=8)
print(f"read and extracted {len(meta)} files in {(time.time() - t) / 60:.1f} min; errors: {len(errors)}")
display(errors)

assets in 0.211014.0809: 1328 | in 0.211014.0808: 1328 | same paths and sizes: True
archive summary: {'numberOfFiles': 1328, 'numberOfSubjects': 266, 'numberOfCells': 1328, 'numberOfBytes': 11922334254}

read and extracted 1328 files in 15.4 min; errors: 0

print("files:", len(meta), "| distinct cell ids:", meta.cell_id.nunique(), "| distinct session ids:", meta.session_id.nunique(),
      "| distinct mice:", meta["subject.subject_id"].nunique(), "| distinct slices:", meta.slice_id.nunique())
print("cells per mouse:", meta.groupby("subject.subject_id").size().describe()[["min", "50%", "max"]].to_dict())
print("recording dates per mouse (max):", meta.groupby("subject.subject_id")["session_start_time"].nunique().max())
print("\nsampling rate (Hz):", meta.rate.astype(str).value_counts().to_dict())
print("trace length (samples):", meta.lengths.astype(str).value_counts().to_dict())
st_nan = meta.starting_times.map(lambda d: d["nan"])
st_val = meta.starting_times.map(lambda d: len(d["values"]))
print("series with starting_time NaN:", int(st_nan.sum()), "of", int(meta.n_sweeps.sum()), "| files with any numeric starting_time:", int((st_val > 0).sum()))
print("voltage unit as stored:", meta.voltage_units.astype(str).value_counts().to_dict(),
      "| conversion:", meta.conversion.astype(str).value_counts().to_dict())
print(f"voltage values: min {meta.vmin.min():.1f}, max {meta.vmax.max():.1f}; median pre-stimulus value per file "
      f"{meta.v_baseline.quantile(0.05):.1f} to {meta.v_baseline.quantile(0.95):.1f} (5th to 95th percentile)")
print("stimulus unit as stored:", stim.unit.value_counts().to_dict(), "| stimulus rate:", stim.rate.value_counts().to_dict(),
      "| stimulus length (samples):", stim.length.value_counts().to_dict())
print(f"stimulus amplitudes: {stim.amplitude.min():.0f} to {stim.amplitude.max():.0f}; distinct levels per series:",
      stim.levels.value_counts().to_dict())
print("stimulus onset sample:", stim.onset.value_counts().to_dict(), "| offset sample:", stim.offset.value_counts().to_dict())
print("sweeps per file:", meta.n_sweeps.value_counts().sort_index().to_dict())
print("sweeps without a stimulus series per file:", (meta.n_sweeps - stim.groupby("path").size().reindex(meta.index)).value_counts().to_dict())

files: 1328 | distinct cell ids: 1328 | distinct session ids: 1328 | distinct mice: 266 | distinct slices: 1260
cells per mouse: {'min': 1.0, '50%': 5.0, 'max': 12.0}
recording dates per mouse (max): 1

sampling rate (Hz): {'25000.0': 1328}
trace length (samples): {'[25000]': 790, '[22500]': 538}
series with starting_time NaN: 104409 of 104409 | files with any numeric starting_time: 0
voltage unit as stored: {"['volts']": 1328} | conversion: {'[1.0]': 1328}
voltage values: min -1025.4, max 1025.7; median pre-stimulus value per file -78.5 to -54.7 (5th to 95th percentile)
stimulus unit as stored: {'amperes': 103081} | stimulus rate: {25000.0: 103081} | stimulus length (samples): {22500: 103081}
stimulus amplitudes: -320 to 1780; distinct levels per series: {2: 103081}
stimulus onset sample: {2500: 103081} | offset sample: {17500: 103081}
sweeps per file: {24: 1, 30: 2, 35: 1, 38: 3, 40: 5, 41: 2, 42: 3, 43: 1, 45: 6, 47: 2, 48: 3, 50: 1, 51: 1, 52: 1, 55: 1, 56: 1, 57: 2, 58: 3, 59: 1, 60: 3, 61: 2, 62: 5, 63: 4, 65: 2, 66: 1, 67: 3, 68: 2, 69: 2, 71: 2, 72: 3, 73: 1, 74: 1, 76: 1, 78: 3, 79: 4, 80: 1248, 100: 1}
sweeps without a stimulus series per file: {1: 1328}

# Metadata fields: presence and values
fields = [c for c in meta.columns if c in ["RNA_family", "RNA_type", "RNA_type_confidence", "RNA_type_top3", "ALM_VISp_top3",
          "cortical_thickness_4x", "cortical_thickness_um", "soma_depth_um", "inferred_layer", "targeted_layer", "exon_reads",
          "intron_reads", "intergenic_reads", "number_of_genes_detected", "sequencing_batch", "slice_id", "cell_id",
          "subject.age", "subject.cre", "subject.date_of_birth", "subject.genotype", "subject.sex", "subject.species",
          "subject.subject_id", "experimenter", "institution", "lab"]]
rows = []
for c in fields:
    v = meta[c]
    num = pd.to_numeric(v, errors="coerce")
    rows.append({"field": c, "non-missing": int(v.notna().sum()), "distinct": int(v.astype(str).nunique()),
                 "range": f"{num.min():.4g} to {num.max():.4g}" if num.notna().sum() > 0.5 * v.notna().sum() and v.notna().any() else "",
                 "example": str(v.dropna().iloc[0])[:40] if v.notna().any() else ""})
display(pd.DataFrame(rows))
print("RNA_family values:", meta.RNA_family.value_counts().to_dict())
added = {"has_soma_depth_um", "session_id", "experimenter", "institution", "lab", "experiment_description",
         "session_start_time", "session_description", "identifier", "rate", "lengths", "starting_times", "voltage_units",
         "conversion", "n_sweeps", "sweeps", "vmin", "vmax", "v_baseline", "hash_all", "hash_first", "cell_id"}
icephys_fields = sorted(c for c in meta.columns if not c.startswith("subject.") and c not in added - {"cell_id"})
subject_fields = sorted(c for c in meta.columns if c.startswith("subject."))
print("every DandiIcephysMetadata field in any file:", icephys_fields)
print("every subject field in any file:", subject_fields)
f0 = scala.open_asset(assets.asset_id.iloc[0])
gen = []
f0["general"].visit(gen.append)
print(f"every entry under 'general' in {assets.path.iloc[0]}:", gen)
print("fields whose name mentions hold, yield or cDNA:",
      [c for c in icephys_fields + subject_fields + gen if any(k in c.lower() for k in ["hold", "yield", "cdna"])])
print("layers (inferred):", meta.inferred_layer.value_counts().to_dict(), "| (targeted):", meta.targeted_layer.value_counts().to_dict())

RNA_family values: {'Pvalb': 289, 'Sst': 272, 'IT': 254, 'Vip': 153, 'CT': 105, 'low quality': 97, 'Lamp5': 91, 'PT': 49, 'Sncg': 13, 'NP': 5}
every DandiIcephysMetadata field in any file: ['ALM_VISp_top3', 'RNA_family', 'RNA_type', 'RNA_type_confidence', 'RNA_type_top3', 'cell_id', 'cortical_thickness_4x', 'cortical_thickness_um', 'exon_reads', 'inferred_layer', 'intergenic_reads', 'intron_reads', 'number_of_genes_detected', 'sequencing_batch', 'slice_id', 'soma_depth_um', 'targeted_layer']
every subject field in any file: ['subject.age', 'subject.cre', 'subject.date_of_birth', 'subject.genotype', 'subject.sex', 'subject.species', 'subject.subject_id']

every entry under 'general' in sub-mouse-AAYYT/sub-mouse-AAYYT_ses-20180420-sample-2_slice-20180420-slice-2_cell-20180420-sample-2_icephys.nwb: ['DandiIcephysMetadata', 'DandiIcephysMetadata/RNA_family', 'DandiIcephysMetadata/cell_id', 'DandiIcephysMetadata/cortical_thickness_4x', 'DandiIcephysMetadata/cortical_thickness_um', 'DandiIcephysMetadata/exon_reads', 'DandiIcephysMetadata/inferred_layer', 'DandiIcephysMetadata/intergenic_reads', 'DandiIcephysMetadata/intron_reads', 'DandiIcephysMetadata/number_of_genes_detected', 'DandiIcephysMetadata/sequencing_batch', 'DandiIcephysMetadata/slice_id', 'DandiIcephysMetadata/targeted_layer', 'devices', 'devices/device', 'experiment_description', 'experimenter', 'institution', 'intracellular_ephys', 'intracellular_ephys/elec0', 'intracellular_ephys/elec0/description', 'intracellular_ephys/sweep_table', 'intracellular_ephys/sweep_table/id', 'intracellular_ephys/sweep_table/series', 'intracellular_ephys/sweep_table/series_index', 'intracellular_ephys/sweep_table/sweep_number', 'lab', 'session_id', 'subject', 'subject/age', 'subject/cre', 'subject/date_of_birth', 'subject/genotype', 'subject/sex', 'subject/subject_id']
fields whose name mentions hold, yield or cDNA: []
layers (inferred): {'5': 585, '2/3': 434, '6': 258, '1': 51} | (targeted): {'5': 576, '2/3': 427, '6': 259, '1': 66}

d4 = pd.to_numeric(meta.cortical_thickness_4x, errors="coerce")
th = pd.to_numeric(meta.cortical_thickness_um, errors="coerce")
ratio = d4 / th
print(f"cells with both values: {ratio.notna().sum()}; 'cortical_thickness_4x' / 'cortical_thickness_um' ranges "
      f"{ratio.min():.3f} to {ratio.max():.3f} (all below 1: {bool((ratio.dropna() < 1).all())})")
print(f"'cortical_thickness_um' {th.min():.0f} to {th.max():.0f}; 'cortical_thickness_4x' {d4.min():.0f} to {d4.max():.0f}")
lay = scala.layer_from_depth(ratio.values, meta.targeted_layer.astype(str).values)
agree = (lay == meta.inferred_layer.astype(str).values)
print(f"layer from that ratio with the Nissl boundaries agrees with 'inferred_layer' for {agree[ratio.notna().values].sum()} "
      f"of {ratio.notna().sum()} cells with the ratio, and for {agree[ratio.isna().values].sum()} of {ratio.isna().sum()} "
      "cells without it (the targeted layer)")
print("'inferred_layer' equals 'targeted_layer' for", int((meta.inferred_layer.astype(str) == meta.targeted_layer.astype(str)).sum()), "cells")

cells with both values: 1283; 'cortical_thickness_4x' / 'cortical_thickness_um' ranges 0.031 to 0.977 (all below 1: True)
'cortical_thickness_um' 415 to 706; 'cortical_thickness_4x' 16 to 675
layer from that ratio with the Nissl boundaries agrees with 'inferred_layer' for 1264 of 1283 cells with the ratio, and for 45 of 45 cells without it (the targeted layer)
'inferred_layer' equals 'targeted_layer' for 1168 cells

dup_all = meta.hash_all.duplicated(keep=False).sum()
dup_first = meta.hash_first.duplicated(keep=False).sum()
print("files sharing an identical voltage matrix with another file:", dup_all, "| sharing an identical first sweep:", dup_first)
dups = scala.duplicate_recordings(meta)
display(dups.drop(columns="hash_all"))
dropped = list(dups.cell_id[~dups.kept])
print("cells whose recordings are left out so that each recording counts once:", dropped)
big = meta[(meta.vmax > 500) | (meta.vmin < -500)]
print("files with voltage values beyond +-500 mV:", list(big.cell_id), f"| range {big.vmin.min():.1f} to {big.vmax.max():.1f}")

# Timing check on the negative-current sweeps (10 ms bins, 0 to 0.9 s)
def window(b, t0, t1):
    return np.mean(b[int(round(t0 * 100)):int(round(t1 * 100))])

rows = []
stim_by = {p: s.set_index("sweep") for p, s in stim.groupby("path")}
for p, sweeps in binned.items():
    s = stim_by[p]
    amp, on_true, on_shift, off_true, off_shift = [], [], [], [], []
    for sw, b in sweeps.items():
        b = np.asarray(b)
        amp.append(s.loc[int(sw), "amplitude"])
        on_true.append(window(b, 0.10, 0.15) - window(b, 0.05, 0.10))
        on_shift.append(window(b, 0.05, 0.10) - window(b, 0.00, 0.05))
        off_true.append(window(b, 0.70, 0.75) - window(b, 0.65, 0.70))
        off_shift.append(window(b, 0.65, 0.70) - window(b, 0.60, 0.65))
    if len(amp) >= 3:
        amp = np.array(amp)
        rows.append({"path": p, **{k: 1000 * np.polyfit(amp, np.array(v), 1)[0] for k, v in
                     [("onset", on_true), ("onset_shifted_50ms", on_shift), ("offset", off_true), ("offset_shifted_50ms", off_shift)]}})
timing = pd.DataFrame(rows).set_index("path")
print(f"cells with at least 3 negative-current sweeps: {len(timing)}")
print("slope of voltage change on current (mV per nA), median and fraction above 20 mV/nA:")
display(pd.DataFrame({"median": timing.median(), "fraction > 20": (timing > 20).mean(), "fraction < -20": (timing < -20).mean()}).round(3))

files sharing an identical voltage matrix with another file: 6 | sharing an identical first sweep: 6

cells whose recordings are left out so that each recording counts once: ['20190731_sample_8', '20181116_sample_4', '20190829_sample_6']
files with voltage values beyond +-500 mV: ['20190625_sample_13'] | range -1025.4 to 1025.7

cells with at least 3 negative-current sweeps: 1328
slope of voltage change on current (mV per nA), median and fraction above 20 mV/nA:

tab = scala.paper_meta_table()
ids = meta.cell_id.values
print("table rows:", len(tab), "| dandiset cells:", len(ids), "| matched by cell id:", int(np.isin(ids, tab.index).sum()),
      "| table cells not in the dandiset:", list(tab.index[~tab.index.isin(ids)]))
m = meta.set_index("cell_id")
tm = tab.loc[m.index]

def norm(v):
    if v is None or (isinstance(v, float) and np.isnan(v)):
        return "nan"
    try:
        return repr(float(v))        # exact: numbers are compared as floats, without rounding
    except (TypeError, ValueError):
        return str(v).strip()

pairs = [("subject.subject_id", "Mouse"), ("slice_id", "Slice"), ("subject.genotype", "Mouse genotype"),
         ("subject.sex", "Mouse gender"), ("subject.cre", "Cre"), ("targeted_layer", "Targeted layer"),
         ("inferred_layer", "Inferred layer"), ("exon_reads", "Exon reads"), ("intron_reads", "Intron reads"),
         ("intergenic_reads", "Intergenic reads"), ("sequencing_batch", "Sequencing batch"),
         ("number_of_genes_detected", "Number of genes detected"), ("RNA_type_confidence", "RNA type confidence"),
         ("ALM_VISp_top3", "ALM/VISp top-3"), ("RNA_family", "RNA family"), ("RNA_type", "RNA type"),
         ("cortical_thickness_4x", "Soma depth (4x)"), ("cortical_thickness_um", "Cortical thickness (4x)"),
         ("cortical_thickness_4x", "Soma depth (µm)"), ("cortical_thickness_um", "Cortical thickness (µm)")]
rows = []
for a, b in pairs:
    x = m[a].map(norm)
    y = tm[b].map(norm)
    x_paper = m[a].map(lambda v: norm(scala.paper_name(v)))
    rows.append({"file field": a, "table column": b, "equal": int((x == y).sum()),
                 "equal in the paper's spelling": int((x_paper == y).sum()), "of": len(m)})
link = pd.DataFrame(rows)
display(link)
age = m["subject.age"].str.extract(r"P(\d+)D")[0].astype(int)
print("age in days equal to the table's 'Mouse age':", int((age == tm["Mouse age"]).sum()), "of", len(m))
print("\ntable 'Exclusion reasons':")
display(tab.loc[tab["Exclusion reasons"].notna(), ["RNA family", "RNA type", "Exclusion reasons"]])
excluded = set(tab.index[tab["Exclusion reasons"].notna()])
print("excluded cells with a t-type:", sum(1 for c in excluded if isinstance(tab.loc[c, "RNA type"], str)),
      "| of them in the dandiset:", sum(1 for c in excluded if isinstance(tab.loc[c, "RNA type"], str) and c in m.index))
print("missing cell:", tab.loc[~tab.index.isin(ids), ["Mouse", "RNA family", "RNA type", "Soma depth (µm)", "Traced"]].to_dict("index"))

table rows: 1329 | dandiset cells: 1328 | matched by cell id: 1328 | table cells not in the dandiset: ['20190211_sample_10']

age in days equal to the table's 'Mouse age': 1328 of 1328

table 'Exclusion reasons':

excluded cells with a t-type: 5 | of them in the dandiset: 5
missing cell: {'20190211_sample_10': {'Mouse': 'mouse_PYVMS', 'RNA family': 'CT', 'RNA type': 'L6 CT Gpr139', 'Soma depth (µm)': '1046.7', 'Traced': 'y'}}

pft = scala.paper_feature_table()
t_all = scala.cell_table(meta, feats, excluded)
common = t_all.index.intersection(pft.index)
print("published rows:", len(pft), "| dandiset cells:", len(t_all), "| matched:", len(common))
rows = []
for f in scala.FEATURES_29:
    x = t_all.loc[common, f].astype(float).values
    y = pft.loc[common, f].astype(float).values
    eq = np.isclose(x, y, rtol=1e-4, atol=1e-6) | (np.isnan(x) & np.isnan(y))
    both = np.isfinite(x) & np.isfinite(y)
    rows.append({"feature": f, "equal": int(eq.sum()), "of": len(common), "largest difference": np.nanmax(np.abs(x - y)) if both.any() else 0})
fl = pd.DataFrame(rows)
display(fl.round(4))
non_rheo = [f for f in scala.FEATURES_29 if f != "Rheobase (pA)"]
x = t_all.loc[common, non_rheo].astype(float).values
y = pft.loc[common, non_rheo].astype(float).values
eqcell = (np.isclose(x, y, rtol=1e-4, atol=1e-6) | (np.isnan(x) & np.isnan(y))).all(axis=1)
print(f"cells equal on all 28 features other than rheobase: {eqcell.sum()} of {len(common)}")
diff = pd.Index(common)[~eqcell]
display(t_all.loc[diff, ["mouse", "ttype", "Max number of APs", "Input resistance (MOhm)"]].join(
    pft.loc[diff, ["Max number of APs", "Input resistance (MOhm)"]], rsuffix=" (published)").join(meta.set_index("cell_id").loc[diff, ["n_sweeps"]]))
for cid in diff:
    xd = t_all.loc[cid, non_rheo].astype(float)
    yd = pft.loc[cid, non_rheo].astype(float)
    bad = [f for f in non_rheo if not (np.isclose(xd[f], yd[f], rtol=1e-4, atol=1e-6) or (np.isnan(xd[f]) and np.isnan(yd[f])))]
    print(f"{cid}: {len(bad)} features differ: {bad if len(bad) <= 6 else str(len(bad)) + ' features'}")
r = t_all.loc[common, "Rheobase (pA)"] - pft.loc[common, "Rheobase (pA)"]
print(f"rheobase: equal for {int(np.isclose(r, 0).sum())} cells; differences elsewhere {r[~np.isclose(r, 0)].min():.1f} to {r[~np.isclose(r, 0)].max():.1f} pA")

published rows: 1328 | dandiset cells: 1328 | matched: 1328

cells equal on all 28 features other than rheobase: 1320 of 1328

20190610_sample_3: 6 features differ: ['Spike frequency adaptation', 'ISI Fano factor', 'ISI coefficient of variation', 'AP Fano factor', 'AP coefficient of variation', 'Wildness']
20190821_sample_2: 1 features differ: ['Membrane time constant (ms)']
20190906_sample_6: 5 features differ: ['ISI Fano factor', 'ISI coefficient of variation', 'AP Fano factor', 'AP coefficient of variation', 'Wildness']
20190513_sample_1: 5 features differ: ['Spike frequency adaptation', 'ISI Fano factor', 'ISI coefficient of variation', 'AP Fano factor', 'AP coefficient of variation']
20190128_sample_3: 25 features differ: 25 features
20180215_sample_1: 22 features differ: 22 features
20190823_sample_5: 5 features differ: ['Max number of APs', 'ISI Fano factor', 'ISI coefficient of variation', 'AP Fano factor', 'AP coefficient of variation']
20190222_sample_5: 25 features differ: 25 features
rheobase: equal for 1287 cells; differences elsewhere -40.0 to 20.0 pA

# The rheobase uses RANSAC, which draws random subsets and has no seed in the paper's code. Re-running the
# cells whose rheobase differs with other seeds shows whether the published value is one RANSAC can give.
pc = scala.paper_code()
rdiff_all = list(pd.Index(common)[~np.isclose(r, 0)])
rdiff = rdiff_all[:8]   # the first eight, ten seeds each, to keep the run time down
hits = 0
for cid in rdiff:
    f = scala.open_asset(assets.set_index("path").loc[t_all.loc[cid, "path"], "asset_id"])
    cell = scala.read_cell(f)
    vals = []
    for seed in range(10):
        import warnings
        tm_, v_, cur_, i0_ = scala.arrange_as_code(cell)
        np.random.seed(seed)
        with warnings.catch_warnings(), np.errstate(all="ignore"):
            warnings.simplefilter("ignore")
            df_, dfr_ = pc.extract_spike_features(tm_, cur_, v_)
            vals.append(float(pc.get_cell_features(df_, dfr_, tm_, cur_, v_, i0_)["Rheobase (pA)"].values[0]))
        pc.plt.close("all")
    hits += np.any(np.isclose(vals, pft.loc[cid, "Rheobase (pA)"], rtol=1e-4))
print(f"cells with a different rheobase: {len(rdiff_all)}; of the first {len(rdiff)}, the published value is reached with one of seeds 0 to 9 for {hits}")

cells with a different rheobase: 41; of the first 8, the published value is reached with one of seeds 0 to 9 for 4

print("distinct recorded cells:", meta.cell_id.nunique(), "| files:", len(meta), "| distinct subjects:", meta["subject.subject_id"].nunique())
print("reading (a) 1,329 cells:", meta.cell_id.nunique() == 1329, "| reading (b) 1,328 cells:", meta.cell_id.nunique() == 1328,
      "| 266 mice:", meta["subject.subject_id"].nunique() == 266)
print("distinct recordings (identical voltage matrices counted once):", meta.hash_all.nunique())
print("the six cells of the three duplicated pairs, as their files describe them:")
display(meta.set_index("cell_id").loc[dups.cell_id, ["subject.subject_id", "slice_id", "RNA_type", "exon_reads",
                                                     "number_of_genes_detected", "cortical_thickness_4x"]])

distinct recorded cells: 1328 | files: 1328 | distinct subjects: 266
reading (a) 1,329 cells: False | reading (b) 1,328 cells: True | 266 mice: True
distinct recordings (identical voltage matrices counted once): 1325
the six cells of the three duplicated pairs, as their files describe them:

mm = meta.assign(age=meta["subject.age"].str.extract(r"P(\d+)D")[0].astype(int))
g = mm.groupby("subject.subject_id")
print("mice with more than one age value:", int((g.age.nunique() > 1).sum()))
ages = g.age.first()
q = {m_: np.percentile(ages, [25, 50, 75], method=m_) for m_ in ["linear", "lower", "higher", "nearest"]}
print(f"n = {len(ages)} mice; min {ages.min()}, max {ages.max()} days")
display(pd.DataFrame(q, index=["Q1", "median", "Q3"]))
print("per cell instead of per mouse (not the claim's reading):", np.percentile(mm.age, [0, 25, 50, 75, 100]))

mice with more than one age value: 0
n = 266 mice; min 35, max 245 days

per cell instead of per mouse (not the claim's reading): [ 35.  64.  75. 100. 245.]

gm = meta.groupby("subject.subject_id").agg(genotype=("subject.genotype", "first"), n_genotypes=("subject.genotype", "nunique"),
                                            cre=("subject.cre", lambda s: ", ".join(sorted(set(s)))))
print("mice with more than one genotype:", int((gm.n_genotypes > 1).sum()))
display(gm.groupby("genotype").agg(mice=("cre", "size"), cre_labels=("cre", lambda s: sorted(set(", ".join(s).split(", "))))))
paper = {"wild type": ("wt/wt", 27), "Viaat-Cre": ("Viaat-Cre/wt; Ai9/wt", 24), "Sst-Cre": ("Sst-IRES-Cre/wt; Ai9/wt", 75),
         "Vip-Cre": ("Vip-IRES-Cre/wt; Ai9/wt", 46), "Pvalb-Cre": ("Pvalb-Cre/wt; Ai9/wt", 76), "Npy-Cre": ("Npy-Cre/wt; Ai9/wt", 2),
         "Vipr2-Cre": ("Vipr2-IRES2-Cre-D/wt; Ai9/wt", 7), "Slc17a8-Cre": ("Slc17a8-IRES2-Cre-D/wt; Ai9/wt", 6),
         "Gnb4-Cre": ("Gnb4-IRES2-CreER2-D knock-in/wt; Ai9/wt", 1), "Slc17a8-iCre": ("Vglut3-Cre/wt; Ai9/wt", 2)}
counts = gm.genotype.value_counts()
cmp_ = pd.DataFrame([{"line (paper)": k, "genotype in the files": v[0], "paper": v[1], "dandiset": int(counts.get(v[0], 0))} for k, v in paper.items()])
cmp_["match"] = cmp_.paper == cmp_.dandiset
display(cmp_)
print("matching lines:", int(cmp_.match.sum()), "of 10 | total mice:", len(gm),
      "| genotypes in the files not in the paper's list:", sorted(set(counts.index) - {v[0] for v in paper.values()}))

mice with more than one genotype: 0

matching lines: 9 of 10 | total mice: 266 | genotypes in the files not in the paper's list: ['Scnn1a-Tg3-Cre/wt; Ai9/wt']

rows = []
for p, s in stim.groupby("path"):
    s = s.sort_values("sweep")
    steps = s.amplitude.values
    sw = s.sweep.values
    rows.append({"path": p, "first": steps[0], "max": steps.max(),
                 "dur_ok": bool(np.all(np.abs((s.offset - s.onset) / s.rate * 1000 - 600) <= 1)),
                 "first_ok": bool(np.isclose(steps[0], -200)),
                 "step_ok": bool(np.allclose(steps, steps[0] + 20 * (sw - sw[0]))),
                 "max_ok": bool(steps.max() <= 1380),
                 "repeat": bool((np.diff(steps) < 0).any())})
sp_ = pd.DataFrame(rows).set_index("path")
sp_["conforms"] = sp_[["dur_ok", "first_ok", "step_ok", "max_ok"]].all(axis=1)
print("cells:", len(sp_))
print(sp_[["dur_ok", "first_ok", "step_ok", "max_ok", "conforms", "repeat"]].mean().round(4).to_dict())
print("first step (pA):", sp_["first"].value_counts().sort_index().to_dict())
print("largest step (pA):", sp_["max"].value_counts().sort_index().to_dict())
first_c = sp_["first"].map(lambda x: f"{x:.0f}" if x in (-200, -100) else "other")
max_c = sp_["max"].map(lambda x: "up to 1380" if x <= 1380 else f"{x:.0f}" if x == 1480 else "other above 1380")
print("first step against largest step (cells):")
display(pd.crosstab(first_c, max_c, margins=True))
print("cells with a second series starting again at -200 pA:", int(sp_.repeat.sum()))
print("intervals between pulses: not computable, as no series has a start time (see the data check)")

cells: 1328
{'dur_ok': 1.0, 'first_ok': 0.8908, 'step_ok': 1.0, 'max_ok': 0.9157, 'conforms': 0.8901, 'repeat': 0.0}
first step (pA): {-320.0: 3, -300.0: 2, -280.0: 1, -260.0: 1, -240.0: 4, -220.0: 8, -200.0: 1183, -180.0: 19, -160.0: 7, -140.0: 6, -120.0: 2, -100.0: 90, -80.0: 2}
largest step (pA): {280.0: 1, 380.0: 1, 420.0: 1, 480.0: 1, 540.0: 1, 560.0: 1, 580.0: 6, 600.0: 2, 620.0: 2, 640.0: 3, 660.0: 1, 680.0: 3, 720.0: 2, 740.0: 3, 760.0: 1, 780.0: 1, 800.0: 2, 840.0: 1, 860.0: 1, 880.0: 1, 920.0: 1, 940.0: 1, 960.0: 1, 980.0: 4, 1000.0: 1, 1020.0: 6, 1040.0: 3, 1080.0: 2, 1100.0: 1, 1120.0: 3, 1140.0: 3, 1160.0: 1, 1180.0: 1, 1200.0: 2, 1220.0: 4, 1240.0: 1, 1260.0: 3, 1280.0: 2, 1300.0: 1, 1320.0: 1, 1340.0: 5, 1360.0: 8, 1380.0: 1126, 1400.0: 14, 1420.0: 4, 1440.0: 4, 1460.0: 1, 1480.0: 87, 1500.0: 1, 1780.0: 1}
first step against largest step (cells):

cells with a second series starting again at -200 pA: 0
intervals between pulses: not computable, as no series has a start time (see the data check)

depth4x = pd.to_numeric(meta.cortical_thickness_4x, errors="coerce")
print("files with the field 'soma_depth_um':", int(meta.has_soma_depth_um.sum()), "| with a value in it:",
      int(pd.to_numeric(meta.get("soma_depth_um", pd.Series(dtype=float)), errors="coerce").notna().sum()))
print("cells with a soma depth (the field 'cortical_thickness_4x', shown above to hold it):", int(depth4x.notna().sum()))
missing = tab.loc[~tab.index.isin(meta.cell_id)]
print("the cell missing from the dandiset has a depth in the paper's table:", bool(pd.to_numeric(missing["Soma depth (µm)"], errors="coerce").notna().all()))
print("paper's table, cells with a µm soma depth:", int(pd.to_numeric(tab["Soma depth (µm)"], errors="coerce").notna().sum()))

files with the field 'soma_depth_um': 43 | with a value in it: 0
cells with a soma depth (the field 'cortical_thickness_4x', shown above to hold it): 1283
the cell missing from the dandiset has a depth in the paper's table: True
paper's table, cells with a µm soma depth: 1284

ex_id = "20181112_sample_11"
ex = t_all.loc[ex_id]
print("file:", ex.path)
legend = [("Input resistance (MOhm)", 274.80, 27.48), ("Membrane time constant (ms)", 21.95, 2.195),
          ("ISI adaptation index", 1.26, 0.126), ("ISI average adaptation index", 1.15, 0.115),
          ("AP amplitude adaptation index", 0.91, 0.091), ("AP amplitude average adaptation index", 0.99, 0.099),
          ("AP threshold (mV)", -40.18, 4.018), ("AP amplitude (mV)", 81.17, 8.117), ("AP width (ms)", 0.80, 0.08),
          ("Afterhyperpolarization (mV)", -12.60, 1.26), ("Afterdepolarization (mV)", 0, 0),
          ("Upstroke-to-downstroke ratio", 1.62, 0.162), ("Latency (ms)", 69.28, 6.928), ("Rheobase (pA)", 20.44, 2.044),
          ("Max number of APs", 32, 3.2), ("ISI coefficient of variation", 0.27, 0.027), ("ISI Fano factor", 0.0014, 0.00014),
          ("AP coefficient of variation", 0.17, 0.017), ("AP Fano factor", 1.32, 0.132), ("Sag ratio", 1.17, 0.117),
          ("Sag time (s)", 0.26, 0.026), ("Sag area (mV*s)", 31.16, 3.116), ("Rebound (mV)", 17.84, 1.784)]
e4 = pd.DataFrame([{"feature": f, "legend": v, "tolerance": tol, "measured": ex[f], "published table": pft.loc[ex_id, f]}
                   for f, v, tol in legend])
e4["within"] = np.abs(e4.measured - e4.legend) <= e4.tolerance + 1e-12
display(e4.round(4))
print("within tolerance:", int(e4.within.sum()), "of", len(e4))

file: sub-mouse-ZWZEZ/sub-mouse-ZWZEZ_ses-20181112-sample-11_slice-20181112-slice-7_cell-20181112-sample-11_icephys.nwb

within tolerance: 23 of 23

t = scala.cell_table(meta, feats, excluded, drop=dropped)      # each recording once
t_paper = scala.cell_table(meta, feats, excluded)               # all cells, as the paper
TABLES = [("each recording once", t), ("all cells, as the paper", t_paper)]
for label, tt in TABLES:
    X = scala.features17(tt)
    h17 = X.notna().all(axis=1)
    typ = (tt.ttype != "") & ~tt.excluded
    print(f"{label}: cells with all 17 features {int(h17.sum())}; of these with a t-type and not excluded {int((h17 & typ).sum())}; "
          f"with a t-type, exclusion not applied (files alone) {int((h17 & (tt.ttype != '')).sum())}")
X17 = scala.features17(t)
have17 = X17.notna().all(axis=1)
typed = (t.ttype != "") & ~t.excluded
print("features undefined, by feature (all cells):", (~np.isfinite(scala.features17(t_paper))).sum()[lambda s: s > 0].to_dict())

each recording once: cells with all 17 features 1317; of these with a t-type and not excluded 1216; with a t-type, exclusion not applied (files alone) 1221
all cells, as the paper: cells with all 17 features 1320; of these with a t-type and not excluded 1219; with a t-type, exclusion not applied (files alone) 1224
features undefined, by feature (all cells): {'ISI adaptation index': 8, 'ISI coefficient of variation': 8, 'AP amplitude adaptation index': 8, 'AP coefficient of variation': 4}

feat3 = ["Max number of APs", "Sag ratio", "Rebound (mV)"]

def family_medians(tt):
    X = scala.features17(tt)
    d = tt[(tt.ttype != "") & ~tt.excluded & X.notna().all(axis=1)]
    med_a = d[d.family.isin(scala.FAMILIES_7)].groupby("family")[feat3].median()
    four = ["Pvalb", "Sst", "Vip", "Lamp5"]
    grp_b = d.family.where(d.family.isin(four), other="excitatory").where(d.family.isin(four + scala.EXCITATORY))
    med_b = d.groupby(grp_b)[feat3].median()
    return d, med_a, med_b

for label, tt in TABLES:
    d, med_a, med_b = family_medians(tt)
    if label == TABLES[0][0]:
        print("(a) seven families, medians:"); display(med_a.round(3))
        print("(b) five groups of Fig. 2, medians:"); display(med_b.round(3))
        for fam_, f in [("Pvalb", "Max number of APs"), ("Sst", "Sag ratio"), ("Sst", "Rebound (mV)")]:
            p = ss.mannwhitneyu(d.loc[d.family == fam_, f], d.loc[d.family != fam_, f], alternative="two-sided").pvalue
            print(f"for information, {fam_} against all other cells, {f}: Mann-Whitney p = {p:.2g}")
        print("cells:", len(d), "| by family:", d.family.value_counts().to_dict())
        print("cells come from", d.mouse.nunique(), "mice; Pvalb from", d[d.family == "Pvalb"].mouse.nunique(),
              "and Sst from", d[d.family == "Sst"].mouse.nunique())
    for name, med in [("(a)", med_a), ("(b)", med_b)]:
        print(f"{label} {name}: highest max APs {med['Max number of APs'].idxmax()}, highest sag ratio {med['Sag ratio'].idxmax()}, "
              f"highest rebound {med['Rebound (mV)'].idxmax()}")

(a) seven families, medians:

(b) five groups of Fig. 2, medians:

for information, Pvalb against all other cells, Max number of APs: Mann-Whitney p = 4.8e-138
for information, Sst against all other cells, Sag ratio: Mann-Whitney p = 4.5e-38
for information, Sst against all other cells, Rebound (mV): Mann-Whitney p = 3.1e-59
cells: 1216 | by family: {'Pvalb': 287, 'Sst': 271, 'IT': 246, 'Vip': 153, 'CT': 102, 'Lamp5': 90, 'ET': 49, 'Sncg': 13, 'NP': 5}
cells come from 264 mice; Pvalb from 116 and Sst from 83
each recording once (a): highest max APs Pvalb, highest sag ratio ET, highest rebound Sst
each recording once (b): highest max APs Pvalb, highest sag ratio Sst, highest rebound Sst
all cells, as the paper (a): highest max APs Pvalb, highest sag ratio ET, highest rebound Sst
all cells, as the paper (b): highest max APs Pvalb, highest sag ratio Sst, highest rebound Sst

# By mouse: the family medians over mice (one median per mouse and family first), for the three features
d, _, _ = family_medians(t)
per_mouse = d[d.family.isin(scala.FAMILIES_7)].groupby(["family", "mouse"])[feat3].median().groupby("family").median()
display(per_mouse.round(3))

def family_confusions(tt, show=False):
    Zr = scala.ephys_representation(scala.features17(tt))
    zz = tt.loc[Zr.index]
    out = {}
    for handling in ["methods", "code"]:
        lab = zz.family.values.astype(object).copy()
        if handling == "methods":
            lab[zz.excluded.values] = "low quality"
        C = scala.family_confusion(scala.knn_predictions(Zr.values, lab, k=10), lab)
        off = C[scala.FAMILIES_7].values.copy()
        np.fill_diagonal(off, np.nan)
        o = pd.DataFrame(off, index=scala.FAMILIES_7, columns=scala.FAMILIES_7)
        et_it, it_et = C.loc["ET", "IT"], C.loc["IT", "ET"]
        o.loc["ET", "IT"] = o.loc["IT", "ET"] = np.nan
        rows_n = int(pd.Series(lab).isin(scala.FAMILIES_7).sum())
        out[handling] = C
        if show:
            print(f"\n{handling} handling of the excluded cells: row percentages (n = cells counted)")
            display((C[scala.FAMILIES_7] * 100).round(1).join(C["n"]))
        print(f"{handling}: cells {len(Zr)}; ET to IT {et_it:.3f}, IT to ET {it_et:.3f}; largest other off-diagonal entry "
              f"{np.nanmax(o.values):.3f} ({o.stack().idxmax()}); other entries at 25% or more: {int((o.values >= 0.25).sum())}; "
              f"cells in the seven rows {rows_n}, not counted {rows_n - int(C.n.sum())}")
    return Zr, out

for label, tt in TABLES:
    print(f"--- {label}")
    Zr, conf_ = family_confusions(tt, show=(label == TABLES[0][0]))
    if label == TABLES[0][0]:
        Z, conf = Zr, conf_
zt = t.loc[Z.index]

--- each recording once

methods handling of the excluded cells: row percentages (n = cells counted)

methods: cells 1317; ET to IT 0.490, IT to ET 0.012; largest other off-diagonal entry 0.144 (('Lamp5', 'Sst')); other entries at 25% or more: 0; cells in the seven rows 1198, not counted 0

code handling of the excluded cells: row percentages (n = cells counted)

code: cells 1317; ET to IT 0.490, IT to ET 0.012; largest other off-diagonal entry 0.144 (('Lamp5', 'Sst')); other entries at 25% or more: 0; cells in the seven rows 1203, not counted 0
--- all cells, as the paper
methods: cells 1320; ET to IT 0.490, IT to ET 0.012; largest other off-diagonal entry 0.143 (('Lamp5', 'Sst')); other entries at 25% or more: 0; cells in the seven rows 1201, not counted 0
code: cells 1320; ET to IT 0.490, IT to ET 0.012; largest other off-diagonal entry 0.143 (('Lamp5', 'Sst')); other entries at 25% or more: 0; cells in the seven rows 1206, not counted 0

# Chance reference: the same classifier with family labels permuted across cells (seed fixed)
rng = np.random.default_rng(31)
lab = zt.family.values.astype(object).copy(); lab[zt.excluded.values] = "low quality"
perm = rng.permutation(lab)
Cp = scala.family_confusion(scala.knn_predictions(Z.values, perm, k=10), perm)
print("mean diagonal, true labels:", round(float(np.nanmean(np.diag(conf['methods'][scala.FAMILIES_7].values))), 3),
      "| permuted labels:", round(float(np.nanmean(np.diag(Cp[scala.FAMILIES_7].values))), 3))
pred_m = pd.Series(scala.knn_predictions(Z.values, lab, k=10), index=Z.index)
et = zt[(zt.family == "ET") & ~zt.excluded]
print("ET cells:", len(et), "from", et.mouse.nunique(), "mice; predicted IT:", int((pred_m[et.index] == "IT").sum()),
      "| mice with an ET cell predicted IT:", et[pred_m[et.index] == "IT"].mouse.nunique())

mean diagonal, true labels: 0.808 | permuted labels: 0.157

ET cells:

 49 from 28 mice; predicted IT: 24 | mice with an ET cell predicted IT: 16

target = {"Sag ratio", "Input resistance (MOhm)", "AP width (ms)", "Rebound (mV)"}
# Layer as the claims file's opening comment defines it (normalized depth cut at the Nissl boundaries, the targeted
# layer without a depth; scala.cell_table, column `layer`), from the files' 4x depths. Two variants: the files' own
# `inferred_layer`, and the same rule on the paper's µm depths (mini-atlas table), which the code uses.
um_nd = (pd.to_numeric(tab["Soma depth (µm)"], errors="coerce")
         / pd.to_numeric(tab["Cortical thickness (µm)"], errors="coerce").replace(0, np.nan))

def layer_defs(tt):
    return [("depth from the files", tt.layer), ("files' inferred_layer", tt.inferred_layer),
            ("depth from the paper's µm table", pd.Series(scala.layer_from_depth(um_nd.reindex(tt.index).values,
                                                                                  tt.targeted_layer.values), index=tt.index))]

lt = layer_defs(t_paper)
print("cells whose layer differs from the main definition:",
      {name: int((lay != lt[0][1]).sum()) for name, lay in lt[1:]},
      "| of them IT or ET:", {name: int(((lay != lt[0][1]) & t_paper.family.isin(["IT", "ET"])).sum()) for name, lay in lt[1:]})
ed6 = {}
for label, tt in TABLES:
    for lname, lay in layer_defs(tt):
        it5 = tt[(tt.family == "IT") & (lay == "5") & ~tt.excluded]
        et5 = tt[(tt.family == "ET") & (lay == "5") & ~tt.excluded]
        rows = []
        for f in scala.FEATURES_29:
            if f == "Rebound number of APs":
                continue
            x, y = it5[f].dropna(), et5[f].dropna()
            mwu = ss.mannwhitneyu(x, y, alternative="two-sided")
            rows.append({"feature": f, "U": mwu.statistic, "p": mwu.pvalue, "median IT": x.median(), "median ET": y.median(),
                         "n IT": len(x), "n ET": len(y)})
        mw = pd.DataFrame(rows).sort_values("p").reset_index(drop=True)
        if label == TABLES[0][0] and lname == "depth from the files":
            print("L5 IT cells:", len(it5), "from", it5.mouse.nunique(), "mice | L5 ET cells:", len(et5), "from", et5.mouse.nunique(), "mice")
            display(mw.head(12))
        kept = mw[~mw.feature.isin(["Upstroke-to-downstroke ratio", "Sag time (s)", "Sag area (mV*s)"])]
        top4 = list(kept.feature[:4])
        ed6[(label, lname)] = kept
        print(f"{label}; layer {lname} (L5 IT {len(it5)}, L5 ET {len(et5)}): four smallest p after the omissions {top4}; "
              f"of the claim's four: {len(target & set(top4))}; fourth p = {kept.p.iloc[3]:.2g}; "
              f"fifth: {kept.feature.iloc[4]} (p = {kept.p.iloc[4]:.2g})")

cells whose layer differs from the main definition: {"files' inferred_layer": 19, "depth from the paper's µm table": 19} | of them IT or ET: {"files' inferred_layer": 3, "depth from the paper's µm table": 3}

L5 IT cells: 111 from 61 mice | L5 ET cells: 46 from 26 mice

each recording once; layer depth from the files (L5 IT 111, L5 ET 46): four smallest p after the omissions ['Sag ratio', 'Input resistance (MOhm)', 'AP width (ms)', 'Rebound (mV)']; of the claim's four: 4; fourth p = 1.8e-05; fifth: Afterhyperpolarization (mV) (p = 8.7e-05)
each recording once; layer files' inferred_layer (L5 IT 111, L5 ET 45): four smallest p after the omissions ['Sag ratio', 'Input resistance (MOhm)', 'AP width (ms)', 'Rebound (mV)']; of the claim's four: 4; fourth p = 2.9e-05; fifth: Afterhyperpolarization (mV) (p = 7.7e-05)
each recording once; layer depth from the paper's µm table (L5 IT 111, L5 ET 45): four smallest p after the omissions ['Sag ratio', 'Input resistance (MOhm)', 'AP width (ms)', 'Rebound (mV)']; of the claim's four: 4; fourth p = 2.9e-05; fifth: Afterhyperpolarization (mV) (p = 7.7e-05)

all cells, as the paper; layer depth from the files (L5 IT 111, L5 ET 46): four smallest p after the omissions ['Sag ratio', 'Input resistance (MOhm)', 'AP width (ms)', 'Rebound (mV)']; of the claim's four: 4; fourth p = 1.8e-05; fifth: Afterhyperpolarization (mV) (p = 8.7e-05)

all cells, as the paper; layer files' inferred_layer (L5 IT 111, L5 ET 45): four smallest p after the omissions ['Sag ratio', 'Input resistance (MOhm)', 'AP width (ms)', 'Rebound (mV)']; of the claim's four: 4; fourth p = 2.9e-05; fifth: Afterhyperpolarization (mV) (p = 7.7e-05)

all cells, as the paper; layer depth from the paper's µm table (L5 IT 111, L5 ET 45): four smallest p after the omissions ['Sag ratio', 'Input resistance (MOhm)', 'AP width (ms)', 'Rebound (mV)']; of the claim's four: 4; fourth p = 2.9e-05; fifth: Afterhyperpolarization (mV) (p = 7.7e-05)

ty = t_paper[(t_paper.ttype != "") & ~t_paper.excluded]
print("cells with a t-type, not excluded:", len(ty), "| inhibitory:", int(ty.family.isin(scala.INHIBITORY).sum()),
      "| excitatory:", int(ty.family.isin(scala.EXCITATORY).sum()), "| t-types covered:", ty.ttype.nunique())
print("files alone (no exclusion):", int((t_paper.ttype != "").sum()), "cells,", t_paper.loc[t_paper.ttype != "", "ttype"].nunique(), "t-types")
mt = tab.loc[~tab.index.isin(meta.cell_id)]
mtype = mt["RNA type"].iloc[0]
print("the cell missing from the dandiset:", mt.index[0], "|", mtype, "| other cells of that t-type in the dandiset:",
      int((ty.ttype == mtype).sum()))
ptab = tab[tab["RNA type"].notna() & tab["Exclusion reasons"].isna()]
print("paper's table for comparison:", len(ptab), "cells,", ptab["RNA type"].nunique(), "t-types")

cells with a t-type, not excluded: 1226 | inhibitory: 817 | excitatory: 409 | t-types covered: 76
files alone (no exclusion): 1231 cells, 76 t-types
the cell missing from the dandiset: 20190211_sample_10 | L6 CT Gpr139 | other cells of that t-type in the dandiset: 0
paper's table for comparison: 1227 cells, 77 t-types

def name_layers(name):
    if name.startswith("L2/3"): return {"2/3"}
    if name.startswith("L4/5"): return {"2/3", "5"}
    if name.startswith("L5/6"): return {"5", "6"}
    if name.startswith("L5"): return {"5"}
    if name.startswith("L6"): return {"6"}
    return None

def layer_of(x):
    return "1" if x < 0.07 else "2/3" if x < 0.29 else "5" if x < 0.73 else "6"

tp = t_paper
typed_p = (tp.ttype != "") & ~tp.excluded
um = (pd.to_numeric(tab["Soma depth (µm)"], errors="coerce") / pd.to_numeric(tab["Cortical thickness (µm)"], errors="coerce").replace(0, np.nan)).reindex(tp.index)
print(f"difference between the two normalized depths: median {np.nanmedian(np.abs(tp.depth_norm - um)):.4f}, "
      f"{int((np.abs(tp.depth_norm - um) > 0.05).sum())} cells differ by more than 0.05")
fig1b = {}
for label, dn in [("files (4x)", tp.depth_norm), ("paper's table (µm)", um)]:
    rows = []
    exc = tp[typed_p & tp.family.isin(scala.EXCITATORY)]
    for name, g in exc.groupby("ttype"):
        v = dn[g.index].dropna()
        if len(v) < 3:
            continue
        L = layer_of(v.median())
        rows.append({"t-type": name, "n with depth": len(v), "median": round(v.median(), 3), "layer": L,
                     "name allows": "/".join(sorted(name_layers(name))), "agrees": L in name_layers(name)})
    fig1b[label] = pd.DataFrame(rows)
    print(label, ": agreeing t-types", int(fig1b[label].agrees.sum()), "of", len(fig1b[label]),
          f"({fig1b[label].agrees.mean():.0%})")
display(fig1b["files (4x)"].merge(fig1b["paper's table (µm)"][["t-type", "median", "agrees"]], on="t-type", suffixes=("", " (µm)")))
l54 = tp[typed_p & (tp.ttype == "L5 IT_4")]
print("L5 IT_4 cells: depths", l54.depth_norm.round(3).tolist(), "| mice", l54.mouse.tolist())

difference between the two normalized depths: median 0.0050, 14 cells differ by more than 0.05

files (4x) : agreeing t-types 18 of 19 (95%)
paper's table (µm) : agreeing t-types 18 of 19 (95%)

L5 IT_4 cells: depths [0.753, 0.758, 0.755, 0.835] | mice ['mouse_BICAU', 'mouse_EBFUB', 'mouse_WVNFF', 'mouse_ZKETK']

for label, dn in [("files (4x)", tp.depth_norm), ("paper's table (µm)", um)]:
    v = dn[tp[typed_p & (tp.ttype == "L4/5 IT_1")].index].dropna()
    q1, med, q3 = v.quantile([0.25, 0.5, 0.75])
    print(f"{label}: n = {len(v)}, median {med:.3f} (reading a: within 0.235 to 0.40: {0.235 <= med <= 0.40}), "
          f"IQR {q1:.3f} to {q3:.3f} (reading b: contains 0.29: {q1 <= 0.29 <= q3}); range {v.min():.3f} to {v.max():.3f}")
l451 = tp[typed_p & (tp.ttype == "L4/5 IT_1")]
print("L4/5 IT_1 cells come from", l451.mouse.nunique(), "mice, at most", l451.mouse.value_counts().max(), "per mouse; median of per-mouse medians",
      round(l451.groupby("mouse").depth_norm.median().median(), 3))

files (4x): n = 22, median 0.279 (reading a: within 0.235 to 0.40: True), IQR 0.259 to 0.306 (reading b: contains 0.29: True); range 0.161 to 0.370
paper's table (µm): n = 22, median 0.277 (reading a: within 0.235 to 0.40: True), IQR 0.252 to 0.303 (reading b: contains 0.29: True); range 0.164 to 0.366
L4/5 IT_1 cells come from 14 mice, at most 4 per mouse; median of per-mouse medians 0.287

import io
from sklearn.decomposition import PCA
mf = pd.read_csv(io.StringIO(scala._get(f"{scala.MINI_ATLAS}/data/m1_patchseq_morph_features.csv").text)).set_index("cell id")
mf = mf.iloc[:, :-2].astype(float)   # the code drops the last two columns (normalized layer, cell class)
zp = pd.read_csv(io.StringIO(scala._get(f"{scala.MINI_ATLAS}/data/m1_patchseq_morph_zprofiles.csv").text)).set_index("cell id").astype(float)
print("cells with morphometrics:", len(mf), "| of them in the dandiset:", int(mf.index.isin(t.index).sum()))

def morph_representation(tt):
    # As phenotype-tsne.ipynb at the pinned commit. Cell 10: cells with any morphometric value, less two cells
    # without all dendritic features. Cell 11: inhibitory or excitatory family and not excluded; per class, the
    # features present in that class, standardized (np.std, ddof 0), 20 PCs scaled by the s.d. of PC1, and 0.25
    # added to the excitatory block. Cell 12: z-profiles, PCs 2 to 5 of a 5-component PCA, scaled by the s.d. of
    # the first kept PC, 0.25 added to the excitatory block. Cell 13: the four blocks in a 48-column matrix,
    # zero outside a cell's own class.
    fam = tt.family.where(~tt.excluded, "low quality")
    q = [c for c in tt.index if c in mf.index and mf.loc[c].notna().any() and c not in ["20180820_sample_1", "20180921_sample_3"]]
    inh = [c for c in q if fam[c] in scala.INHIBITORY]
    exc = [c for c in q if fam[c] in scala.EXCITATORY]
    def block(cs):
        M = mf.loc[cs]
        M = M.loc[:, M.notna().sum() > 0]
        M = (M - M.mean()) / M.std(ddof=0)
        P = PCA(n_components=20).fit_transform(M.values); P /= np.std(P[:, 0])
        Zp = PCA(n_components=5).fit_transform(zp.loc[cs].values)[:, 1:]; Zp /= np.std(Zp[:, 0])
        return P, Zp
    Pi, Zi = block(inh)
    Pe, Ze = block(exc)
    Mx = np.zeros((len(inh) + len(exc), 48))
    Mx[:len(inh), 0:20] = Pi; Mx[len(inh):, 20:40] = Pe + .25; Mx[:len(inh), 40:44] = Zi; Mx[len(inh):, 44:48] = Ze + .25
    return pd.DataFrame(Mx, index=inh + exc), len(inh), len(exc)

def separability(tt, show=False):
    Zr = scala.ephys_representation(scala.features17(tt))
    M, ni, ne = morph_representation(tt)
    res_ = {}
    for handling in ["methods", "code"]:
        fam = tt.family.copy()
        if handling == "methods":
            fam[tt.excluded] = "low quality"
        both = [c for c in M.index if c in Zr.index]
        for ncomp in [16, 17]:
            # Cell 17: the morphological matrix next to the ephys representation (cell 7, all 17 PCs scaled by the
            # s.d. of PC1) times sqrt(2), for the cells with both. The Methods' "16-dimensional" is the 16-PC variant.
            comb = np.concatenate([M.loc[both].values, Zr.loc[both].values[:, :ncomp] * np.sqrt(2)], axis=1)
            out = {}
            for name, R, cs in [("ephys", Zr.values, Zr.index), ("morphology", M.values, M.index), ("combined", comb, both)]:
                lab = fam[list(cs)].values.astype(object)
                C = scala.family_confusion(scala.knn_predictions(R, lab, k=10), lab)
                dg = np.diag(C[scala.FAMILIES_7].values.astype(float))
                out[name] = {"cells": len(cs), "inhibitory mean diagonal": np.nanmean(dg[:4]), "seven-family mean diagonal": np.nanmean(dg)}
            r_ = pd.DataFrame(out).T
            res_[(handling, ncomp)] = r_
            if show:
                print(f"\nexcluded cells: {handling}; ephys block of the combined representation: {ncomp} PCs "
                      f"(morphology: {ni} inhibitory, {ne} excitatory cells)")
                display(r_.round(3))
            print(f"{handling}, {ncomp} PCs: inhibitory, morphology below ephys: "
                  f"{r_.loc['morphology', 'inhibitory mean diagonal'] < r_.loc['ephys', 'inhibitory mean diagonal']}; "
                  f"seven families, combined highest: {r_['seven-family mean diagonal'].idxmax() == 'combined'}")
    return res_

for label, tt in TABLES:
    print(f"--- {label}")
    morph_res = separability(tt, show=True)

cells with morphometrics: 646 | of them in the dandiset: 645
--- each recording once

excluded cells: methods; ephys block of the combined representation: 16 PCs (morphology: 367 inhibitory, 268 excitatory cells)

methods, 16 PCs: inhibitory, morphology below ephys: True; seven families, combined highest: True

excluded cells: methods; ephys block of the combined representation: 17 PCs (morphology: 367 inhibitory, 268 excitatory cells)

methods, 17 PCs: inhibitory, morphology below ephys: True; seven families, combined highest: True

excluded cells: code; ephys block of the combined representation: 16 PCs (morphology: 367 inhibitory, 268 excitatory cells)

code, 16 PCs: inhibitory, morphology below ephys: True; seven families, combined highest: True

excluded cells: code; ephys block of the combined representation: 17 PCs (morphology: 367 inhibitory, 268 excitatory cells)

code, 17 PCs: inhibitory, morphology below ephys: True; seven families, combined highest: True
--- all cells, as the paper

excluded cells: methods; ephys block of the combined representation: 16 PCs (morphology: 367 inhibitory, 268 excitatory cells)

methods, 16 PCs: inhibitory, morphology below ephys: True; seven families, combined highest: True

excluded cells: methods; ephys block of the combined representation: 17 PCs (morphology: 367 inhibitory, 268 excitatory cells)

methods, 17 PCs: inhibitory, morphology below ephys: True; seven families, combined highest: True

excluded cells: code; ephys block of the combined representation: 16 PCs (morphology: 367 inhibitory, 268 excitatory cells)

def type_confusion(tt, handling, show=False):
    Zr = scala.ephys_representation(scala.features17(tt))
    zz = tt.loc[Zr.index]
    lab = zz.family.values.astype(object).copy()
    ttl = zz.ttype.values.astype(object).copy()
    if handling == "methods":
        lab[zz.excluded.values] = "low quality"
        ttl[zz.excluded.values] = ""
    pred = scala.knn_predictions(Zr.values, lab, k=10)
    rows = []
    for name in sorted(set(ttl) - {""}):
        ind = ttl == name
        if ind.sum() < 10:
            continue
        p = pred[ind]
        p = p[np.isin(p, scala.FAMILIES_7)]
        own = zz.family.values[ind][0]
        row = {"t-type": name, "family": own, "cells": int(ind.sum()), "counted": len(p)}
        row.update({f: np.mean(p == f) for f in scala.FAMILIES_7})
        row["wrong"] = 1 - row.get(own, 0)
        rows.append(row)
    R5 = pd.DataFrame(rows).set_index("t-type")
    pooled = pred[pd.Series(ttl).str.startswith("Vip Mybpc1").values]
    pooled = pooled[np.isin(pooled, scala.FAMILIES_7)]
    if show:
        display(R5.round(2))
    ii = R5.loc["Sst Pvalb Calb2", "Pvalb"] if "Sst Pvalb Calb2" in R5.index else np.nan
    iii = R5.loc["Lamp5 Egln3_1", ["Vip", "Sst"]].sum() if "Lamp5 Egln3_1" in R5.index else np.nan
    iv_a = {n: round(R5.loc[n, "Sst"], 3) for n in R5.index if n.startswith("Vip Mybpc1")}
    print(f"{handling}: t-types with at least 10 cells {len(R5)}; (i) at most 25% wrong: {int((R5.wrong <= 0.25).sum())}; "
          f"(ii) Sst Pvalb Calb2 as Pvalb {ii:.3f}; (iii) Lamp5 Egln3_1 as Vip or Sst {iii:.3f}; "
          f"(iv) Vip Mybpc1 as Sst, (a) per t-type {iv_a}, (b) pooled ({len(pooled)} cells) {np.mean(pooled == 'Sst'):.3f}")
    return R5

for label, tt in TABLES:
    print(f"--- {label}")
    for handling in ["methods", "code"]:
        R5 = type_confusion(tt, handling, show=(label == TABLES[0][0] and handling == "methods"))

--- each recording once

methods: t-types with at least 10 cells 35; (i) at most 25% wrong: 28; (ii) Sst Pvalb Calb2 as Pvalb 0.385; (iii) Lamp5 Egln3_1 as Vip or Sst 0.900; (iv) Vip Mybpc1 as Sst, (a) per t-type {'Vip Mybpc1_2': np.float64(0.235), 'Vip Mybpc1_3': np.float64(0.375)}, (b) pooled (59 cells) 0.254
code: t-types with at least 10 cells 36; (i) at most 25% wrong: 29; (ii) Sst Pvalb Calb2 as Pvalb 0.385; (iii) Lamp5 Egln3_1 as Vip or Sst 0.900; (iv) Vip Mybpc1 as Sst, (a) per t-type {'Vip Mybpc1_2': np.float64(0.235), 'Vip Mybpc1_3': np.float64(0.375)}, (b) pooled (59 cells) 0.254
--- all cells, as the paper

methods: t-types with at least 10 cells 35; (i) at most 25% wrong: 28; (ii) Sst Pvalb Calb2 as Pvalb 0.385; (iii) Lamp5 Egln3_1 as Vip or Sst 0.900; (iv) Vip Mybpc1 as Sst, (a) per t-type {'Vip Mybpc1_2': np.float64(0.235), 'Vip Mybpc1_3': np.float64(0.375)}, (b) pooled (59 cells) 0.254

code: t-types with at least 10 cells 36; (i) at most 25% wrong: 29; (ii) Sst Pvalb Calb2 as Pvalb 0.385; (iii) Lamp5 Egln3_1 as Vip or Sst 0.900; (iv) Vip Mybpc1 as Sst, (a) per t-type {'Vip Mybpc1_2': np.float64(0.235), 'Vip Mybpc1_3': np.float64(0.375)}, (b) pooled (59 cells) 0.254

from sklearn.cluster import KMeans
nbmeta = scala._get(f"{scala.MINI_ATLAS}/code/phenotype-tsne.ipynb").json()["metadata"]
print("phenotype-tsne.ipynb ran under Python", nbmeta["language_info"]["version"])

def variability(tt, handling, n_init):
    Zr = scala.ephys_representation(scala.features17(tt))
    zz = tt.loc[Zr.index]
    Zv = Zr.values
    ttl = zz.ttype.values.astype(object).copy()
    if handling == "methods":
        ttl[zz.excluded.values] = ""
    names_ = [n for n in sorted(set(ttl) - {""}) if (ttl == n).sum() >= 10]
    nv = pd.Series({n: scala.normalized_variance(Zv, ttl == n) for n in names_})
    restrict = np.isin(ttl, names_)
    km = KMeans(n_clusters=len(names_), random_state=42, n_init=n_init).fit(Zv[restrict])
    cvar = [np.sum(np.var(Zv[restrict][km.labels_ == k], axis=0)) / np.sum(np.var(Zv, axis=0))
            for k in range(len(names_)) if np.sum(km.labels_ == k) >= 10]
    lo, hi = min(cvar), max(cvar)
    above = nv[nv > hi]
    print(f"{handling} membership, n_init={n_init}: {len(names_)} t-types, {len(cvar)} clusters with at least 10 cells, "
          f"band {lo:.3f} to {hi:.3f}; t-types above the band: {len(above)} ({len(above) / len(names_):.0%}); "
          f"Pvalb Vipr2_2 {nv.get('Pvalb Vipr2_2', np.nan):.3f} (above: {nv.get('Pvalb Vipr2_2', 0) > hi}); "
          f"Vip Mybpc1_2 {nv.get('Vip Mybpc1_2', np.nan):.3f} (above: {nv.get('Vip Mybpc1_2', 0) > hi})")
    return nv, ttl, Zv

fig5b = {}
for label, tt in TABLES:
    print(f"--- {label}")
    for handling in ["methods", "code"]:
        for n_init in [10, "auto"]:
            fig5b[(label, handling, str(n_init))] = variability(tt, handling, n_init)
print("\nnormalized variances (each recording once, methods membership):")
print(fig5b[(TABLES[0][0], "methods", "10")][0].sort_values(ascending=False).round(3).to_dict())

phenotype-tsne.ipynb ran under Python 3.7.9
--- each recording once
methods membership, n_init=10: 35 t-types, 34 clusters with at least 10 cells, band 0.173 to 0.468; t-types above the band: 15 (43%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)
methods membership, n_init=auto: 35 t-types, 33 clusters with at least 10 cells, band 0.178 to 0.542; t-types above the band: 8 (23%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)
code membership, n_init=10: 36 t-types, 35 clusters with at least 10 cells, band 0.188 to 0.549; t-types above the band: 8 (22%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)

code membership, n_init=auto: 36 t-types, 35 clusters with at least 10 cells, band 0.186 to 0.468; t-types above the band: 16 (44%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)
--- all cells, as the paper
methods membership, n_init=10: 35 t-types, 33 clusters with at least 10 cells, band 0.185 to 0.463; t-types above the band: 17 (49%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)
methods membership, n_init=auto: 35 t-types, 32 clusters with at least 10 cells, band 0.194 to 0.483; t-types above the band: 14 (40%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)
code membership, n_init=10: 36 t-types, 35 clusters with at least 10 cells, band 0.165 to 0.453; t-types above the band: 20 (56%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)

code membership, n_init=auto: 36 t-types, 35 clusters with at least 10 cells, band 0.221 to 0.538; t-types above the band: 10 (28%); Pvalb Vipr2_2 0.449 (above: False); Vip Mybpc1_2 0.692 (above: True)

normalized variances (each recording once, methods membership):
{'L5 ET_1': 0.722, 'Vip Mybpc1_2': 0.692, 'Sst Myh8_1': 0.648, 'Vip Mybpc1_3': 0.629, 'Sst Chodl': 0.614, 'L5 IT_1': 0.6, 'Lamp5 Pdlim5_2': 0.592, 'Vip Serpinf1_1': 0.544, 'Vip Sncg': 0.516, 'Pvalb Gpr149': 0.505, 'L4/5 IT_1': 0.504, 'Pvalb Kank4': 0.5, 'Lamp5 Egln3_1': 0.493, 'L2/3 IT_3': 0.489, 'L6 CT Pou3f2': 0.482, 'Sst Pvalb Calb2': 0.467, 'Sst Calb2': 0.463, 'Pvalb Egfem1': 0.456, 'Sst Hpse': 0.453, 'Pvalb Vipr2_2': 0.449, 'Lamp5 Slc35d3': 0.448, 'L4/5 IT_2': 0.448, 'Pvalb Il1rapl2': 0.442, 'Pvalb Calb1_1': 0.434, 'L6 IT_2': 0.431, 'Vip Chat_1': 0.424, 'Pvalb Reln': 0.402, 'Sst Crhr2_1': 0.381, 'Lamp5 Egln3_2': 0.334, 'Sst Crhr2_2': 0.334, 'L6 CT Cpa6': 0.315, 'Vip Gpc3': 0.288, 'Sst Th_3': 0.288, 'Sst Etv1': 0.287, 'Pvalb Calb1_2': 0.271}

import igraph as ig
import leidenalg as la
from scipy.stats import entropy
from sklearn.neighbors import kneighbors_graph

def leiden_entropies(tt, label, show=False):
    Zv = scala.ephys_representation(scala.features17(tt)).values
    A = kneighbors_graph(Zv, 10)
    src, dst = A.nonzero()
    G = ig.Graph(directed=True)
    G.add_vertices(A.shape[0])
    G.add_edges(list(zip(src, dst)))
    parts = np.zeros((len(Zv), 100), dtype=int)
    for i in range(100):
        parts[:, i] = la.find_partition(G, la.RBConfigurationVertexPartition, resolution_parameter=2, seed=i).membership
    ncl = np.array([len(np.unique(parts[:, i])) for i in range(100)])
    if show:
        print(f"Leiden clusters per run: median {np.median(ncl):.0f}, range {ncl.min()} to {ncl.max()}")
    out = {}
    for handling in ["methods", "code"]:
        nv, ttl, _ = fig5b[(label, handling, "10")]
        ent = {}
        for n in nv.index:
            ind = ttl == n
            E = []
            for i in range(100):
                np.random.seed(i)
                sub = np.random.choice(ind.sum(), size=10, replace=False)
                _, pr = np.unique(parts[:, i][ind][sub], return_counts=True)
                E.append(entropy(pr))
            ent[n] = np.median(E)
        ent = pd.Series(ent)
        rho = ss.spearmanr(nv[ent.index], ent)
        p1 = ss.spearmanr(nv[ent.index], ent, alternative="greater").pvalue
        print(f"{label}, {handling} membership: {len(ent)} t-types; Spearman rho between normalized variance and median entropy "
              f"{rho.statistic:.3f}, one-sided p {p1:.3f}")
        out[handling] = (nv, ent)
    return out

for label, tt in TABLES:
    out8 = leiden_entropies(tt, label, show=(label == TABLES[0][0]))
    if label == TABLES[0][0]:
        nv_m, ent_m = out8["methods"]
top_v = set(nv_m.sort_values(ascending=False).index[:10])
top_e = set(ent_m.sort_values(ascending=False).index[:10])
print("t-types among the ten most variable by both measures:", sorted(top_v & top_e))
print("median entropies:", ent_m.sort_values(ascending=False).round(3).to_dict())

Leiden clusters per run: median 17, range 15 to 20

each recording once, methods membership: 35 t-types; Spearman rho between normalized variance and median entropy 0.266, one-sided p 0.061

each recording once, code membership: 36 t-types; Spearman rho between normalized variance and median entropy 0.313, one-sided p 0.031

all cells, as the paper, methods membership: 35 t-types; Spearman rho between normalized variance and median entropy 0.166, one-sided p 0.171

all cells, as the paper, code membership: 36 t-types; Spearman rho between normalized variance and median entropy 0.222, one-sided p 0.097
t-types among the ten most variable by both measures: ['L5 ET_1', 'L5 IT_1', 'Sst Chodl', 'Vip Mybpc1_3']
median entropies: {'Vip Mybpc1_3': 1.67, 'Sst Chodl': 1.643, 'Sst Pvalb Calb2': 1.498, 'L5 ET_1': 1.363, 'Lamp5 Egln3_1': 1.28, 'Pvalb Egfem1': 1.28, 'L5 IT_1': 1.28, 'L4/5 IT_2': 1.254, 'L6 CT Pou3f2': 1.224, 'Sst Crhr2_2': 1.221, 'Sst Crhr2_1': 1.221, 'Pvalb Calb1_1': 1.194, 'Pvalb Vipr2_2': 1.168, 'Lamp5 Egln3_2': 1.168, 'L6 CT Cpa6': 1.129, 'Sst Calb2': 1.089, 'L4/5 IT_1': 1.089, 'Vip Mybpc1_2': 1.089, 'Pvalb Gpr149': 1.03, 'Pvalb Reln': 1.03, 'Pvalb Kank4': 1.03, 'Vip Serpinf1_1': 0.95, 'Pvalb Calb1_2': 0.943, 'Sst Etv1': 0.943, 'Pvalb Il1rapl2': 0.943, 'Sst Hpse': 0.94, 'Lamp5 Pdlim5_2': 0.898, 'Vip Chat_1': 0.802, 'Sst Myh8_1': 0.802, 'L6 IT_2': 0.802, 'Lamp5 Slc35d3': 0.639, 'Sst Th_3': 0.639, 'L2/3 IT_3': 0.625, 'Vip Gpc3': 0.611, 'Vip Sncg': 0.611}

y = pd.to_numeric(tab["Yield (pg/µl)"].replace("?", np.nan), errors="coerce")
h = pd.to_numeric(tab["Hold Time (min)"], errors="coerce")
ok = y.notna() & h.notna() & (y > 0)
r_all = ss.spearmanr(h[ok], np.log(y[ok]))
ok_d = ok & tab.index.isin(meta.cell_id)
r_d = ss.spearmanr(h[ok_d], np.log(y[ok_d]))
print(f"all table cells with both: n = {int(ok.sum())}, rho = {r_all.statistic:.4f}, p = {r_all.pvalue:.3f}")
print(f"dandiset cells with both: n = {int(ok_d.sum())}, rho = {r_d.statistic:.4f}, p = {r_d.pvalue:.3f}")
print(f"hold time {h.min():.0f} to {h.max():.0f} min; yield {y[ok].min():.2f} to {y[ok].max():.1f} pg/µl")

all table cells with both: n = 1320, rho = -0.0124, p = 0.653
dandiset cells with both: n = 1319, rho = -0.0123, p = 0.654
hold time 3 to 85 min; yield 13.30 to 2770.8 pg/µl

# morph-reconstructions: an attempt, from this run, to list the Brain Image Library folder of the reconstructions
bil = "https://download.brainimagelibrary.org/3a/88/3a88a7687ab66069/"
try:
    r = scala.requests.get(bil, timeout=30)
    print("status", r.status_code, "| bytes", len(r.content))
    print(r.text[:2000])
except scala.requests.RequestException as e:
    print(f"{time.strftime('%Y-%m-%d %H:%M UTC', time.gmtime())}: no answer within 30 s:", type(e).__name__)

2026-10-21 06:42 UTC: no answer within 30 s: ConnectTimeout

print(f"run time: {(time.time() - T0) / 60:.1f} min")

run time: 18.5 min

Claim	Data	Here	Tested or blocked
link-metadata-table	link	this notebook	tested
link-ephys-features	link	this notebook	tested
cohort-cells-mice	dataset	this notebook	tested
cohort-ages	dataset	this notebook	tested
cohort-genotypes	dataset	this notebook	tested
stimulus-protocol	dataset	this notebook	tested (intervals not in the files)
soma-depth-count	dataset	this notebook	tested
ed4-example-cell	dataset	this notebook	tested
ephys-feature-count	dataset	this notebook	tested
fig2-family-features	dataset	this notebook	tested
fig3d-families-ephys	dataset	this notebook	tested
ed6-it-et-room-features	dataset	this notebook	tested
ed6-physiological-temperature	elsewhere (DANDI:000035)	none	blocked: the paper gives no measure
fig1a-ttype-coverage	dataset	this notebook	tested
fig1b-depth-layers	dataset	this notebook	tested
l45-it1-boundary	dataset	this notebook	tested
fig4c-distance-correlation	with_inputs	transcriptomics notebook	tested
ed7-within-families	with_inputs	transcriptomics notebook	tested
fig4d-it-depth	with_inputs	transcriptomics notebook	tested
fig4a-vip-tau	with_inputs	transcriptomics notebook	tested
fig4b-sst-rebound	with_inputs	transcriptomics notebook	tested
fig4e-pvalb-axon-shape	elsewhere (Brain Image Library)	none	blocked: no measure in the paper
fig3d-morph-combined	with_inputs	this notebook	tested
fig5a-ttypes-between-families	dataset	this notebook	tested
fig5b-ttype-variability	dataset	this notebook	tested
ed8-leiden-entropy	dataset	this notebook	tested
fig5c-vip-mybpc1-rebound	with_inputs	transcriptomics notebook	tested
fig5d-sst-pvalb-calb2-firing	with_inputs	transcriptomics notebook	tested
rnaseq-depth	elsewhere	sequencing notebook	tested
ttype-qc-counts	elsewhere	sequencing notebook	tested
ed2-smartseq-similarity	elsewhere	sequencing notebook	tested
morph-reconstructions	elsewhere (Brain Image Library)	none	blocked: the archive did not answer
fig2-rrr	with_inputs	transcriptomics notebook	tested
ed5-rrr-ion-channels	with_inputs	transcriptomics notebook	tested
hold-time-yield	with_inputs	this notebook	tested

	subject.subject_id	slice_id	RNA_type	exon_reads	number_of_genes_detected	cortical_thickness_4x
cell_id
20190731_sample_5	mouse_XKIBL	20190731_slice_5	Pvalb Il1rapl2	583064	8250	156.6
20190731_sample_8	mouse_XKIBL	20190731_slice_8	Pvalb Il1rapl2	2231783	11092	217.1
20181116_sample_2	mouse_AZZAO	20181116_slice_2	Pvalb Il1rapl2	239157	11303	118.7
20181116_sample_4	mouse_AZZAO	20181116_slice_2	Pvalb Gpr149	1795312	9820	312.9
20190829_sample_2	mouse_AWOAY	20190829_slice_2	Lamp5 Slc35d3	2424634	8587	284.4
20190829_sample_6	mouse_AWOAY	20190829_slice_4	Lamp5 Slc35d3	1431990	8557	380.5

	feature	legend	tolerance	measured	published table	within
0	Input resistance (MOhm)	274.8000	27.4800	274.8041	274.8041	True
1	Membrane time constant (ms)	21.9500	2.1950	21.9533	21.9533	True
2	ISI adaptation index	1.2600	0.1260	1.2584	1.2584	True
3	ISI average adaptation index	1.1500	0.1150	1.1509	1.1509	True
4	AP amplitude adaptation index	0.9100	0.0910	0.9101	0.9101	True
5	AP amplitude average adaptation index	0.9900	0.0990	0.9853	0.9853	True
6	AP threshold (mV)	-40.1800	4.0180	-43.1840	-43.1840	True
7	AP amplitude (mV)	81.1700	8.1170	81.1725	81.1725	True
8	AP width (ms)	0.8000	0.0800	0.8000	0.8000	True
9	Afterhyperpolarization (mV)	-12.6000	1.2600	-12.6011	-12.6011	True
10	Afterdepolarization (mV)	0.0000	0.0000	0.0000	0.0000	True
11	Upstroke-to-downstroke ratio	1.6200	0.1620	1.6245	1.6245	True
12	Latency (ms)	69.2800	6.9280	69.2800	69.2800	True
13	Rheobase (pA)	20.4400	2.0440	20.4444	20.4444	True
14	Max number of APs	32.0000	3.2000	32.0000	32.0000	True
15	ISI coefficient of variation	0.2700	0.0270	0.2687	0.2687	True
16	ISI Fano factor	0.0014	0.0001	0.0014	0.0014	True
17	AP coefficient of variation	0.1700	0.0170	0.1669	0.1669	True
18	AP Fano factor	1.3200	0.1320	1.3197	1.3197	True
19	Sag ratio	1.1700	0.1170	1.1682	1.1682	True
20	Sag time (s)	0.2600	0.0260	0.2603	0.2603	True
21	Sag area (mV*s)	31.1600	3.1160	31.1597	31.1597	True
22	Rebound (mV)	17.8400	1.7840	17.8408	17.8408	True

	Max number of APs	Sag ratio	Rebound (mV)
family
CT	16.0	1.095	1.894
ET	17.0	1.252	3.089
IT	14.0	1.063	0.572
Lamp5	24.0	1.085	1.937
Pvalb	71.0	1.063	1.544
Sst	27.0	1.154	5.701
Vip	24.0	1.048	3.085

	Max number of APs	Sag ratio	Rebound (mV)
family
Lamp5	24.0	1.085	1.937
Pvalb	71.0	1.063	1.544
Sst	27.0	1.154	5.701
Vip	24.0	1.048	3.085
excitatory	15.0	1.079	1.075

Scala, Kobak et al. (2021): electrophysiology and metadata of dandiset 000008¶

Loading¶

Data check¶

Link claims¶

link-metadata-table¶

link-ephys-features¶

Cohort and recording protocol¶

cohort-cells-mice¶

cohort-ages¶

cohort-genotypes¶

stimulus-protocol¶

soma-depth-count¶

Electrophysiological features¶

ed4-example-cell¶

ephys-feature-count¶

fig2-family-features¶

fig3d-families-ephys¶

ed6-it-et-room-features¶

Transcriptomic coverage¶

fig1a-ttype-coverage¶

fig1b-depth-layers¶

l45-it1-boundary¶

Family separability with morphology¶

fig3d-morph-combined¶

Variability of individual t-types¶

fig5a-ttypes-between-families¶

fig5b-ttype-variability¶

ed8-leiden-entropy¶

Recording procedure¶

hold-time-yield¶

Blocked claims¶

Summary¶

	field	non-missing	distinct	range	example
0	RNA_family	1328	10		low quality
1	cell_id	1328	1328		20180420_sample_2
2	cortical_thickness_4x	1283	1128	16.5 to 675.1	108.5
3	cortical_thickness_um	1283	860	414.8 to 706.5	608.5
4	exon_reads	1328	1328	2 to 1.817e+07	393455
5	inferred_layer	1328	4	1 to 6	2/3
6	intergenic_reads	1328	1322	1 to 2.327e+06	222014
7	intron_reads	1328	1327	4 to 1.327e+07	1275287
8	number_of_genes_detected	1328	1250	11 to 2.483e+04	5789
9	sequencing_batch	1328	12	1 to 12	2
10	slice_id	1328	1260		20180420_slice_2
11	targeted_layer	1328	4	1 to 6	2/3
12	subject.age	1328	102		P148D
13	subject.cre	1328	18		WT
14	subject.date_of_birth	1328	104		2017-11-23T00:00:00-05:00
15	subject.genotype	1328	11		wt/wt
16	subject.sex	1328	2		F
17	subject.subject_id	1328	266		mouse_AAYYT
18	experimenter	1328	2		['Matteo Bernabucci']
19	institution	1328	1		Baylor College of Medicine
20	lab	1328	1		Tolias
21	ALM_VISp_top3	1231	920		L4 IT VISp Rspo1 (82.0%), L5 IT ALM Pld5
22	RNA_type	1231	76		L4/5 IT_1
23	RNA_type_confidence	1231	84	0.14 to 1	0.5
24	RNA_type_top3	1231	1171		L4/5 IT_1 (49.6%), L2/3 IT_3 (31.1%), L5
25	subject.species	1326	1		Mus musculus
26	soma_depth_um	0	0

	median	fraction > 20	fraction < -20
onset	157.875	1.000	0.000
onset_shifted_50ms	0.093	0.006	0.003
offset	-142.167	0.000	1.000
offset_shifted_50ms	-0.298	0.000	0.002

	file field	table column	equal	equal in the paper's spelling	of
0	subject.subject_id	Mouse	1328	1328	1328
1	slice_id	Slice	1328	1328	1328
2	subject.genotype	Mouse genotype	1328	1328	1328
3	subject.sex	Mouse gender	1328	1328	1328
4	subject.cre	Cre	1328	1328	1328
5	targeted_layer	Targeted layer	1328	1328	1328
6	inferred_layer	Inferred layer	1328	1328	1328
7	exon_reads	Exon reads	1328	1328	1328
8	intron_reads	Intron reads	1328	1328	1328
9	intergenic_reads	Intergenic reads	1328	1328	1328
10	sequencing_batch	Sequencing batch	1328	1328	1328
11	number_of_genes_detected	Number of genes detected	1328	1328	1328
12	RNA_type_confidence	RNA type confidence	1328	1328	1328
13	ALM_VISp_top3	ALM/VISp top-3	1328	1328	1328
14	RNA_family	RNA family	1279	1328	1328
15	RNA_type	RNA type	1279	1328	1328
16	cortical_thickness_4x	Soma depth (4x)	1328	1328	1328
17	cortical_thickness_um	Cortical thickness (4x)	1328	1328	1328
18	cortical_thickness_4x	Soma depth (µm)	43	43	1328
19	cortical_thickness_um	Cortical thickness (µm)	0	0	1328

	RNA family	RNA type	Exclusion reasons
Cell
20180315_sample_6	IT	L2/3 IT_3	mismatch (FS firing but excitatory RNA)
20180321_sample_4	IT	L2/3 IT_3	mismatch (MC firing+morphology but excitatory ...
20181116_sample_10	low quality	NaN	mismatch (FS firing but excitatory RNA)
20190325_sample_2	low quality	NaN	mismatch (FS firing but excitatory RNA)
20190611_sample_1	IT	L4/5 IT_1	mismatch (L1 firing+morphology but excitatory ...
20190611_sample_11	IT	L6 IT_1	mismatch (FS firing+morphology but excitatory ...
20190704_sample_2	Sst	Sst Hpse	outside of M1

	feature	equal	of	largest difference
0	Resting membrane potential (mV)	1325	1328	11.4026
1	Input resistance (MOhm)	1325	1328	33.7153
2	Membrane time constant (ms)	1324	1328	3.3352
3	AP threshold (mV)	1325	1328	11.1932
4	AP amplitude (mV)	1325	1328	17.5440
5	AP width (ms)	1325	1328	0.0800
6	Upstroke-to-downstroke ratio	1325	1328	1.0914
7	Afterhyperpolarization (mV)	1325	1328	6.1918
8	Afterdepolarization (mV)	1328	1328	0.0000
9	ISI adaptation index	1325	1328	0.1593
10	Max number of APs	1324	1328	75.0000
11	Rheobase (pA)	1287	1328	40.0000
12	Sag ratio	1325	1328	0.0802
13	Latency (ms)	1325	1328	53.5600
14	Latency @ +20pA current (ms)	1325	1328	11.7200
15	Spike frequency adaptation	1323	1328	0.1250
16	ISI Fano factor	1321	1328	0.0026
17	ISI coefficient of variation	1321	1328	0.0720
18	ISI average adaptation index	1325	1328	0.0441
19	Rebound (mV)	1325	1328	2.7692
20	Sag time (s)	1326	1328	0.4645
21	Sag area (mV*s)	1326	1328	40.4043
22	AP amplitude adaptation index	1325	1328	0.0589
23	AP amplitude average adaptation index	1325	1328	0.0073
24	AP Fano factor	1321	1328	7.3368
25	AP coefficient of variation	1321	1328	0.3210
26	Burstiness	1327	1328	0.0881
27	Wildness	1325	1328	7.0000
28	Rebound number of APs	1328	1328	0.0000

	mouse	ttype	Max number of APs	Input resistance (MOhm)	Max number of APs (published)	Input resistance (MOhm) (published)	n_sweeps
20190610_sample_3	mouse_CCRZS	Lamp5 Pdlim5_2	31.0	199.976754	31	199.976754	80
20190821_sample_2	mouse_CUEES	L4/5 IT_1	22.0	195.145801	22	195.145801	80
20190906_sample_6	mouse_GSHJT	Lamp5 Slc35d3	30.0	239.453514	30	239.453514	80
20190513_sample_1	mouse_HGZZB	Lamp5 Egln3_1	25.0	242.381887	25	242.381887	65
20190128_sample_3	mouse_IGDZU		153.0	170.931986	78	178.614682	41
20180215_sample_1	mouse_LBJDN	Pvalb Reln	51.0	217.842707	50	184.127387	80
20190823_sample_5	mouse_NZDYM		29.0	199.239622	28	199.239622	80
20190222_sample_5	mouse_ZGNZE	L6 CT Cpa6	13.0	209.529841	11	240.644673	80

	mice	cre_labels
genotype
Gnb4-IRES2-CreER2-D knock-in/wt; Ai9/wt	1	[GNB4-]
Npy-Cre/wt; Ai9/wt	2	[NPY+]
Pvalb-Cre/wt; Ai9/wt	76	[PV+, PV-]
Scnn1a-Tg3-Cre/wt; Ai9/wt	1	[SCNN1]
Slc17a8-IRES2-Cre-D/wt; Ai9/wt	6	[SLC17a8+, SLC17a8-]
Sst-IRES-Cre/wt; Ai9/wt	75	[SST+, SST-]
Vglut3-Cre/wt; Ai9/wt	2	[SLC17a8-iCre+, SLC17a8-iCre-]
Viaat-Cre/wt; Ai9/wt	24	[VIAAT+, VIAAT-]
Vip-IRES-Cre/wt; Ai9/wt	46	[VIP+, VIP-]
Vipr2-IRES2-Cre-D/wt; Ai9/wt	7	[VIPR2+, VIPR2-]
wt/wt	26	[WT]

	line (paper)	genotype in the files	paper	dandiset	match
0	wild type	wt/wt	27	26	False
1	Viaat-Cre	Viaat-Cre/wt; Ai9/wt	24	24	True
2	Sst-Cre	Sst-IRES-Cre/wt; Ai9/wt	75	75	True
3	Vip-Cre	Vip-IRES-Cre/wt; Ai9/wt	46	46	True
4	Pvalb-Cre	Pvalb-Cre/wt; Ai9/wt	76	76	True
5	Npy-Cre	Npy-Cre/wt; Ai9/wt	2	2	True
6	Vipr2-Cre	Vipr2-IRES2-Cre-D/wt; Ai9/wt	7	7	True
7	Slc17a8-Cre	Slc17a8-IRES2-Cre-D/wt; Ai9/wt	6	6	True
8	Gnb4-Cre	Gnb4-IRES2-CreER2-D knock-in/wt; Ai9/wt	1	1	True
9	Slc17a8-iCre	Vglut3-Cre/wt; Ai9/wt	2	2	True

	Max number of APs	Sag ratio	Rebound (mV)
family
CT	15.5	1.094	1.889
ET	17.5	1.258	3.283
IT	14.0	1.063	0.583
Lamp5	22.0	1.083	2.500
Pvalb	71.0	1.064	1.535
Sst	26.5	1.158	5.892
Vip	24.5	1.045	3.045

	Lamp5	Vip	Sst	Pvalb	IT	ET	CT	n
Lamp5	65.6	10.0	14.4	0.0	2.2	0.0	7.8	90
Vip	0.7	85.6	11.1	0.7	0.7	0.0	1.3	153
Sst	0.7	1.8	91.5	5.2	0.0	0.0	0.7	271
Pvalb	0.0	0.3	2.1	97.6	0.0	0.0	0.0	287
IT	0.4	0.4	0.4	0.0	91.1	1.2	6.5	246
ET	0.0	0.0	0.0	0.0	49.0	49.0	2.0	49
CT	2.0	1.0	4.9	0.0	6.9	0.0	85.3	102

	linear	lower	higher	nearest
Q1	64.0	64	64	64
median	75.0	75	75	75
Q3	100.0	100	100	100

	feature	U	p	median IT	median ET	n IT	n ET
0	Sag ratio	615.0	7.871150e-14	1.084591	1.279730	111	46
1	Input resistance (MOhm)	4338.0	5.886350e-12	150.192214	79.514414	111	46
2	AP width (ms)	4246.5	6.436852e-11	1.560000	1.280000	111	46
3	Upstroke-to-downstroke ratio	4146.0	8.155200e-10	4.106889	3.514868	111	46
4	Sag time (s)	1492.0	1.523547e-07	0.000000	0.376120	111	46
5	Sag area (mV*s)	1596.0	2.191959e-06	0.000000	33.493019	111	46
6	Rebound (mV)	1441.5	1.827660e-05	1.632227	3.399744	111	46
7	Afterhyperpolarization (mV)	1535.0	8.700223e-05	-3.406979	-0.435464	111	46
8	Latency (ms)	3342.0	2.357557e-03	123.520000	93.400000	111	46
9	AP threshold (mV)	3312.0	3.440733e-03	-35.425793	-39.194165	111	46
10	Rheobase (pA)	1828.5	5.052147e-03	80.000000	100.000000	111	46
11	AP amplitude (mV)	1850.0	6.741152e-03	89.303866	92.343735	111	46

	cells	inhibitory mean diagonal	seven-family mean diagonal
ephys	1317.0	0.851	0.808
morphology	635.0	0.737	0.753
combined	627.0	0.905	0.872

	cells	inhibitory mean diagonal	seven-family mean diagonal
ephys	1317.0	0.851	0.808
morphology	635.0	0.737	0.753
combined	627.0	0.905	0.874

	t-type	n with depth	median	layer	name allows	agrees	median (µm)	agrees (µm)
0	L2/3 IT_3	83	0.169	2/3	2/3	True	0.171	True
1	L4/5 IT_1	22	0.279	2/3	2/3/5	True	0.277	True
2	L4/5 IT_2	33	0.360	5	2/3/5	True	0.365	True
3	L5 ET_1	32	0.546	5	5	True	0.540	True
4	L5 ET_2	4	0.386	5	5	True	0.381	True
5	L5 ET_3	4	0.467	5	5	True	0.471	True
6	L5 ET_4	9	0.420	5	5	True	0.459	True
7	L5 IT_1	58	0.446	5	5	True	0.447	True
8	L5 IT_2	7	0.527	5	5	True	0.523	True
9	L5 IT_3	7	0.545	5	5	True	0.535	True
10	L5 IT_4	4	0.757	6	5	False	0.767	False
11	L6 CT Cpa6	65	0.820	6	6	True	0.814	True
12	L6 CT Pou3f2	16	0.801	6	6	True	0.796	True
13	L6 IT_1	9	0.750	6	6	True	0.757	True
14	L6 IT_2	25	0.822	6	6	True	0.809	True
15	L6b Col6a1	8	0.895	6	6	True	0.893	True
16	L6b Ror1	5	0.911	6	6	True	0.938	True
17	L6b Shisa6_1	4	0.907	6	6	True	0.881	True
18	L6b Shisa6_2	3	0.909	6	6	True	0.900	True

	cells	inhibitory mean diagonal	seven-family mean diagonal
ephys	1320.0	0.847	0.806
morphology	635.0	0.737	0.753
combined	628.0	0.901	0.869

	family	cells	counted	Lamp5	Vip	Sst	Pvalb	IT	ET	CT	wrong
t-type
L2/3 IT_3	IT	79	79	0.00	0.00	0.00	0.00	0.97	0.00	0.03	0.03
L4/5 IT_1	IT	22	22	0.00	0.00	0.00	0.00	0.91	0.00	0.09	0.09
L4/5 IT_2	IT	33	33	0.03	0.00	0.00	0.00	0.82	0.00	0.15	0.18
L5 ET_1	ET	32	32	0.00	0.00	0.00	0.00	0.59	0.38	0.03	0.62
L5 IT_1	IT	59	59	0.00	0.02	0.02	0.00	0.93	0.02	0.02	0.07
L6 CT Cpa6	CT	62	62	0.02	0.00	0.03	0.00	0.02	0.00	0.94	0.06
L6 CT Pou3f2	CT	17	17	0.06	0.00	0.00	0.00	0.24	0.00	0.71	0.29
L6 IT_2	IT	25	25	0.00	0.00	0.00	0.00	0.80	0.00	0.20	0.20
Lamp5 Egln3_1	Lamp5	10	10	0.10	0.40	0.50	0.00	0.00	0.00	0.00	0.90
Lamp5 Egln3_2	Lamp5	10	10	0.60	0.00	0.20	0.00	0.00	0.00	0.20	0.40
Lamp5 Pdlim5_2	Lamp5	21	21	0.81	0.05	0.05	0.00	0.05	0.00	0.05	0.19
Lamp5 Slc35d3	Lamp5	38	38	0.79	0.05	0.05	0.00	0.03	0.00	0.08	0.21
Pvalb Calb1_1	Pvalb	16	16	0.00	0.00	0.06	0.94	0.00	0.00	0.00	0.06
Pvalb Calb1_2	Pvalb	26	26	0.00	0.00	0.00	1.00	0.00	0.00	0.00	0.00
Pvalb Egfem1	Pvalb	11	11	0.00	0.00	0.18	0.82	0.00	0.00	0.00	0.18
Pvalb Gpr149	Pvalb	31	31	0.00	0.00	0.00	1.00	0.00	0.00	0.00	0.00
Pvalb Il1rapl2	Pvalb	84	84	0.00	0.00	0.00	1.00	0.00	0.00	0.00	0.00
Pvalb Kank4	Pvalb	20	20	0.00	0.00	0.00	1.00	0.00	0.00	0.00	0.00
Pvalb Reln	Pvalb	76	76	0.00	0.00	0.00	1.00	0.00	0.00	0.00	0.00
Pvalb Vipr2_2	Pvalb	17	17	0.00	0.06	0.06	0.88	0.00	0.00	0.00	0.12
Sst Calb2	Sst	61	61	0.00	0.00	1.00	0.00	0.00	0.00	0.00	0.00
Sst Chodl	Sst	14	14	0.07	0.14	0.71	0.00	0.00	0.00	0.07	0.29
Sst Crhr2_1	Sst	16	16	0.00	0.00	1.00	0.00	0.00	0.00	0.00	0.00
Sst Crhr2_2	Sst	11	11	0.00	0.00	1.00	0.00	0.00	0.00	0.00	0.00
Sst Etv1	Sst	10	10	0.00	0.00	1.00	0.00	0.00	0.00	0.00	0.00
Sst Hpse	Sst	47	47	0.00	0.00	0.98	0.00	0.00	0.00	0.02	0.02
Sst Myh8_1	Sst	16	16	0.00	0.00	1.00	0.00	0.00	0.00	0.00	0.00
Sst Pvalb Calb2	Sst	26	26	0.04	0.00	0.58	0.38	0.00	0.00	0.00	0.42
Sst Th_3	Sst	18	18	0.00	0.06	0.83	0.11	0.00	0.00	0.00	0.17
Vip Chat_1	Vip	16	16	0.00	0.94	0.00	0.00	0.06	0.00	0.00	0.06
Vip Gpc3	Vip	24	24	0.00	0.92	0.04	0.00	0.00	0.00	0.04	0.08
Vip Mybpc1_2	Vip	34	34	0.00	0.76	0.24	0.00	0.00	0.00	0.00	0.24
Vip Mybpc1_3	Vip	16	16	0.00	0.50	0.38	0.06	0.00	0.00	0.06	0.50
Vip Serpinf1_1	Vip	10	10	0.00	0.90	0.10	0.00	0.00	0.00	0.00	0.10
Vip Sncg	Vip	23	23	0.00	1.00	0.00	0.00	0.00	0.00	0.00	0.00

Claim	Outcome
link-metadata-table	reproduced
link-ephys-features	reproduced
cohort-cells-mice	reproduced
cohort-ages	reproduced
cohort-genotypes	partially reproduced
stimulus-protocol	partially reproduced
soma-depth-count	reproduced
ed4-example-cell	reproduced
ephys-feature-count	reproduced
fig2-family-features	reproduced
fig3d-families-ephys	reproduced
ed6-it-et-room-features	reproduced
ed6-physiological-temperature	blocked (methods underspecified)
fig1a-ttype-coverage	partially reproduced
fig1b-depth-layers	partially reproduced
l45-it1-boundary	reproduced
fig4e-pvalb-axon-shape	blocked (methods underspecified)
fig3d-morph-combined	reproduced (with the paper's morphology tables)
fig5a-ttypes-between-families	reproduced
fig5b-ttype-variability	reproduced
ed8-leiden-entropy	not reproduced
morph-reconstructions	blocked (data unreadable)
hold-time-yield	reproduced (with the paper's table)