Reproducing Chartrand et al. (2023) from dandiset 000630¶
Paper. Chartrand T, Dalley R, Close J, et al. (2023). Morphoelectric and transcriptomic divergence of the layer 1 interneuron repertoire in human versus mouse neocortex. Science 382, eadf0805. doi:10.1126/science.adf0805. The claims file was written from the author manuscript in PubMed Central (PMC11864503) with its supplementary Materials and Methods, and from the paper's code, AllenInstitute/patchseq_human_L1 at commit 7777857c1b89b19735e45fd9c9bd7493cc1d805c, which this entry reads for the cell labels.
Dandiset. DANDI:000630, version
0.230915.2257, "Human L1 patch-seq electrophysiology": 210 NWB files, one human cell each. This is
the version the paper cites (claims file, opening comment) and the only published version; the
draft has the same 210 assets. The mouse comparison data (DANDI:000020) are used in
mouse_comparison.ipynb.
What is tested here. The files hold current clamp recordings (long and short squares, ramps,
in some cells chirps) and, in 38 files, nucleated-patch voltage clamp sweeps. The features are
computed here with IPFX from the long square sweeps (reading (a) of the claims file's open detail
O1); features that rest on code the paper did not release (irregularity ratio, depolarizing sag,
ADP, the ramp and short square features, and others) come from the paper's per-cell table
(reading (b)). A feature computed here is used in reading (a) when it agrees with the paper's
value within 1% for at least half of the cells, an assumption of this attempt (the claims file
does not say when a recomputed feature counts as computed); the five adaptation ratios fail it and
are the paper's values in both readings (a variant, "a, ratios", uses the ratios computed here where
they bear on a claim). Every claim on the human cells is run on both readings and on a variant
(a') of reading (a) with the features extracted without the Methods' 500 ms rms noise criterion,
which the paper's own values show its pipeline did not apply as measured here (see the cells after
dandiset-cell-link).
Cell types come from the paper's cell table, linked to the files by the claim dandiset-cell-link.
| Claim | Notebook | Status |
|---|---|---|
| dandiset-cell-link (added link claim) | this | tested |
| patchseq-counts | this | tested |
| fig2c-subclass-ephys-differences, fig2b-subclass-ephys-profiles, figs3a-sag-adaptation, fig2c-lamp5-gradient | this | tested |
| figs6-donor-condition-ephys, figs6c-region-ephys | this | tested |
| fig4b-mc4r-sag-irregularity, fig4b-adarb2-irregular-firing, fig5b-tnfaip8l3-ephys | this | tested |
| fig6a-ephys-dprime-subclass-structure, fig6d-gmm-subclass-f1, figs8a-classifier-cv | this | tested |
| figs8b-classifier-holdout | this | blocked: the Amsterdam and Szeged recordings are not in the dandiset (two later DANDI datasets are checked) |
| mouse-cell-link (added link claim), fig3d-ephys-species, fig3g-variability-dimensions, fig4d-mouse-no-rosehip-ephys, fig5d-mouse-pax6-doublets | mouse_comparison.ipynb | tested |
| fig3e-tea-validation | mouse_comparison.ipynb | blocked: no public mouse TEa recordings |
| fig3f-a-type-conductance | mouse_comparison.ipynb | blocked: the nucleated-patch sweeps are in the files (38 human, 25 mouse), but the patch surface areas the densities need are not public |
| fig6b-dprime-correlation | none | blocked: not attempted in this run; its inputs are public |
fig2d-subclass-morphology, fig4b-mc4r-axon, fig5b-tnfaip8l3-dendrites, figs6a-tumor-dendrite-branching (data: elsewhere) |
morphology_tables.ipynb | tested on the paper's deposited morphology features |
fig3b-l1-cytoarchitecture (data: elsewhere) |
none | blocked: the code rescales the mouse values with scanner resolutions that are not public |
fig1e-subclass-proportions (data: elsewhere) |
mouse_comparison.ipynb | blocked: the public mouse table holds only cells of L1 t-types, not the deeper t-types the test counts |
12 claims on transcriptomic references, mouse morphology and morphology-based models (data: elsewhere) |
none | blocked: not attempted in this run; their inputs are public |
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
from scipy import stats
import l1ephys as L
pd.set_option("display.width", 220)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
warnings.filterwarnings("ignore")
started = time.time()
meta = L.version_metadata()
assets = L.list_assets()
print(f"dandiset {L.DANDISET_ID} version {L.DANDISET_VERSION}: {meta['name']!r}, "
f"{len(assets)} assets, {assets['size'].sum() / 1e9:.2f} GB")
feats = L.all_cell_features(assets, workers=12, slow_noise_variant=True)
print(f"features extracted from {len(feats)} files in {time.time() - started:.0f} s")
labels = L.load_labels()
paper_raw = L.load_paper_features()
print(f"paper's tables (commit {L.LABELS_COMMIT[:7]}): {len(labels)} labelled cells, "
f"{len(paper_raw)} cells in aibs_features_E.csv")
dandiset 000630 version 0.230915.2257: 'Human L1 patch-seq electrophysiology', 210 assets, 8.07 GB
features extracted from 210 files in 445 s
paper's tables (commit 7777857): 419 labelled cells, 263 cells in aibs_features_E.csv
Data check¶
Each check of the metadata against the values is stated under its output.
s = meta["assetsSummary"]
print("Archive summary: files", s["numberOfFiles"], "| cells", s["numberOfCells"], "| subjects",
s["numberOfSubjects"], "| species", sorted({x["name"] for x in s["species"]}))
print("Files read:", len(feats), "| errors:", int(feats.get("error", pd.Series(dtype=str)).notna().sum()))
print("Distinct subject_id:", feats.subject_id.nunique(), "| distinct session_id:", feats.session_id.nunique(),
"| distinct cell names:", feats.cell_name.nunique())
print("Species in files:", feats.species.value_counts().to_dict())
print("Institution:", feats.institution.value_counts().to_dict())
print("Sex:", feats.sex.value_counts().to_dict())
age_y = feats.age.str.extract(r"P([\d.]+)D")[0].astype(float) / 365.25
print(f"Age: {age_y.min():.1f} to {age_y.max():.1f} years (stated as P<days>D)")
print("Files per subject: ", feats.subject_id.value_counts().describe()[["min", "50%", "max"]].to_dict())
Archive summary: files 210 | cells 210 | subjects 53 | species ['Homo sapiens - Human', 'Human']
Files read: 210 | errors: 0
Distinct subject_id: 53 | distinct session_id: 210 | distinct cell names: 210
Species in files: {'Homo Sapiens': 192, 'Homo sapiens': 18}
Institution: {'Allen Institute for Brain Science': 210}
Sex: {'M': 138, 'F': 71, 'unknown': 1}
Age: 0.0 to 82.9 years (stated as P<days>D)
Files per subject: {'min': 1.0, '50%': 3.0, 'max': 19.0}
The archive's summary counts 210 files, 210 cells and 53 subjects, all human, which the files
confirm. Each file holds one cell. subject_id is a LIMS donor id, and ages are stated in days.
The species is spelled two ways ("Homo Sapiens" and "Homo sapiens"). The two files of donor
922967568 state an age of P0.0D (0 years); the paper's table gives that donor's age as unknown, so
0 stands for a missing value. One file states sex "unknown" where the table gives M.
# Stimulus codes present, and the sweeps used for features
codes = feats.codes.str.split(";").explode().str.split(":").str[0].str.replace("_DA_0", "")
print("Files holding each stimulus code (top 30):")
print(codes.value_counts().head(30).to_string())
print("\nLong square sweeps per file:", feats.n_long_square.describe()[["min", "50%", "max"]].to_dict())
print("Sweeps failing the QC applied here, per file:", feats.n_sweeps_qc_failed.describe()[["min", "50%", "max"]].to_dict())
print("Sweeps used per file:", feats.n_sweeps_used.describe()[["min", "50%", "max"]].to_dict())
print("Long square sweeps stored with zeros before the end of the stimulus (left out):",
int(feats.n_sweeps_zero_padded.sum()), "in", int((feats.n_sweeps_zero_padded > 0).sum()), "files")
print("Stimulus onset (s):", feats.stim_start_s.value_counts().to_dict())
print("Stimulus duration (s):", feats.stim_dur_s.round(3).value_counts().to_dict())
Files holding each stimulus code (top 30):
codes
EXTPBREAKN180424 159
EXTPCllATT180424 159
EXTPEXPEND180424 159
X1PS_SubThresh 159
X2LP_Search 159
X3LP_Rheo 159
X5SP_Search 159
X6SP_Rheo 159
X7Ramp 159
EXTPINBATH180424 158
EXTPSMOKET180424 157
EXTPGGAEND180424 144
X4PS_SupraThresh 142
C1SQCAPCHK141203 100
C2CHIRP180503 65
C1LSFINEST150112 50
C1RP25PR1S141203 50
C1SSFINEST150112 50
EXTPBREAKN141203 50
EXTPCllATT141203 50
EXTPINBATH141203 50
EXTPSMOKET141203 50
EXTPEXPEND141203 49
C2SSTRIPLE171103 48
EXTPGGAEND141203 46
Blip_finder 40
Rheobase_finder 40
subthreshold_171 40
NucVCLS0 38
NucVCLS1Leak 35
Long square sweeps per file: {'min': 9.0, '50%': 20.0, 'max': 179.0}
Sweeps failing the QC applied here, per file: {'min': 0.0, '50%': 5.0, 'max': 133.0}
Sweeps used per file: {'min': 2.0, '50%': 13.0, 'max': 70.0}
Long square sweeps stored with zeros before the end of the stimulus (left out): 2586 in 145 files
Stimulus onset (s): {0.6: 160, 1.02: 50}
Stimulus duration (s): {1.0: 210}
The long squares are 1 s steps (Methods), starting 0.6 s or 1.02 s into the sweep depending on
the protocol version. Sweeps that stopped early hold zeros (0 mV) for the samples not acquired, where
IPFX expects missing values; IPFX's integrity check then passes them and a step to 0 mV is read as
a spike. They are left out here, as counted above. Files also hold short squares, ramps and in some cells chirps and
nucleated-patch voltage clamp protocols (NucVC...), which are not used here.
# Units. Voltages are stored with conversion to volts and currents to amperes; the values read
# here are in mV and pA. bias_current is stated in amperes in the NWB standard.
import h5py, io, requests
r = requests.get(assets.url[0], timeout=300)
f = h5py.File(io.BytesIO(r.content), "r")
k = [k for k in f["acquisition"] if "bias_current" in f["acquisition"][k]][0]
g = f["acquisition"][k]
print("data unit:", g["data"].attrs["unit"], "conversion:", g["data"].attrs["conversion"],
"| IGOR units:", list(g["data"].attrs.get("IGORWaveUnits", [])))
print("bias_current value:", float(np.ravel(g["bias_current"][()])[0]),
"| attrs:", {a: g["bias_current"].attrs[a] for a in g["bias_current"].attrs})
f.close()
print("\nVoltage ranges of features computed here: v_baseline", feats.v_baseline.min().round(1), "to",
feats.v_baseline.max().round(1), "mV; rheobase", feats.rheobase_i.min(), "to", feats.rheobase_i.max(), "pA")
data unit: volts conversion: 0.001 | IGOR units: ['mV', 'ms']
bias_current value: -2.8304920196533203 | attrs: {'unit': 'pA'}
Voltage ranges of features computed here: v_baseline -77.2 to -49.5 mV; rheobase 7.999999968033578 to 599.9999976025183 pA
The bias current is stored in pA, with a unit attribute of "pA", where the NWB schema fixes
the unit of bias_current as amperes; a reader that follows the schema gets values 10^12 times
too large. The QC here reads it as pA. Membrane potentials (-80 to -50 mV baselines) and stimulus amplitudes
are in the expected ranges.
# Clock: the hyperpolarizing response starts at the stimulus onset in the same series.
print("Median voltage change from 100 ms before onset to 100-200 ms after it, hyperpolarizing sweeps (mV):")
print(feats[["clock_dv_true_mv", "clock_dv_shifted_mv"]].describe().round(2).to_string())
print("Cells with a negative change at the true onset:", int((feats.clock_dv_true_mv < -1).sum()), "of",
int(feats.clock_dv_true_mv.notna().sum()))
Median voltage change from 100 ms before onset to 100-200 ms after it, hyperpolarizing sweeps (mV):
clock_dv_true_mv clock_dv_shifted_mv
count 210.00 210.00
mean -26.08 -0.01
std 8.36 0.18
min -56.96 -0.67
25% -31.11 -0.13
50% -24.88 -0.01
75% -19.62 0.10
max -10.38 0.76
Cells with a negative change at the true onset: 210 of 210
Stimulus and response are in the same file with the same sampling rate, and the response follows the stimulus onset: the hyperpolarizing response appears at the true onset and not with the onset taken 0.3 s early. This supports the shared clock of stimulus and response within each sweep.
# Recordings deposited twice: identical feature vectors or identical session start times.
cols = ["sag", "input_resistance", "tau", "v_baseline", "rheobase_i", "threshold_v_rheo", "upstroke_rheo"]
vec = feats[cols].round(6).fillna(-999).apply(tuple, axis=1)
print("Files sharing all of", cols, ":", int(vec.duplicated(keep=False).sum()))
print("Files sharing a session start time:", int(feats.session_start.duplicated(keep=False).sum()))
Files sharing all of ['sag', 'input_resistance', 'tau', 'v_baseline', 'rheobase_i', 'threshold_v_rheo', 'upstroke_rheo'] : 0 Files sharing a session start time: 0
No two files share their features or their start time; no recording is deposited twice.
The dataset's cells and the paper's cell table (link claim dandiset-cell-link)¶
The dandiset's files can be matched one to one to the cells of the paper's per-cell table (human_l1_dataset_2023_02_06.csv), and the matched records agree on what both describe.
The files carry no t-type or subclass. Every claim that compares cell types takes the labels from the paper's table, so they depend on this link.
feats["session_id_int"] = feats.session_id.astype("int64")
by_roi = labels.reset_index().dropna(subset=["ephys_roi_id"]).astype({"ephys_roi_id": "int64"}).set_index("ephys_roi_id")
linked = feats.join(by_roi, on="session_id_int", rsuffix="_table")
print("Files matched by session_id = Ephys_Roi_Result id of the table's nwb_path:",
int(linked.spec_id.notna().sum()), "of", len(linked))
same_name = linked.cell_name == linked.cell_name_table
print("Of these, cell name in the file equal to the table's cell_name:", int(same_name.sum()))
print(linked.loc[~same_name, ["path", "cell_name", "cell_name_table"]].to_string(index=False))
print("Distinct table cells matched:", linked.spec_id.nunique())
sex_ok = linked.sex.str[0] == linked.donor_sex.astype(str).str[0]
print("Sex agrees with the table's donor_sex:", int(sex_ok.sum()), "of", int(linked.donor_sex.notna().sum()))
print("Records whose sex differs:")
print(linked.loc[~sex_ok, ["cell_name", "subject_id", "sex", "donor_sex"]].to_string(index=False))
# Age: the file states days (P<days>D), the table whole years. The days are whole years times 365
# (checked below), so the file's age in years is days / 365, compared exactly.
age_days = linked.age.str.extract(r"P([\d.]+)D")[0].astype(float)
age_tab = linked.donor_age.str.extract(r"(\d+)")[0].astype(float)
print("Files whose age in days is a whole number of 365-day years:", int((age_days % 365 == 0).sum()), "of", int(age_days.notna().sum()))
age_file = age_days / 365
print("Age in years equal to the table's donor_age:", int((age_file == age_tab).sum()), "of",
int(age_tab.notna().sum()), "| largest difference:", float((age_file - age_tab).abs().max()))
age_off = age_tab.notna() & (age_file != age_tab)
print("Records whose age differs:")
print(linked.loc[age_off, ["cell_name", "subject_id", "donor_name", "age", "donor_age"]].assign(
years=age_file[age_off]).to_string(index=False))
print("Files with an age where the table has none:", int((age_days.notna() & age_tab.isna()).sum()),
linked.loc[age_days.notna() & age_tab.isna(), "age"].unique().tolist())
print("Donors per subject_id:", linked.groupby("subject_id").donor_name.nunique().max())
Files matched by session_id = Ephys_Roi_Result id of the table's nwb_path: 210 of 210
Of these, cell name in the file equal to the table's cell_name: 204
path cell_name cell_name_table
sub-1122692698/sub-1122692698_ses-1122718412_icephys.nwb 2021_08_17_112853 H21.03.315.11.11.01.01
sub-1122692698/sub-1122692698_ses-1122740196_icephys.nwb 2021_08_17_140401 H21.03.315.11.11.01.02
sub-1122692698/sub-1122692698_ses-1122745954_icephys.nwb 2021_08_17_142825 H21.03.315.11.11.01.03
sub-1122692698/sub-1122692698_ses-1122753392_icephys.nwb 2021_08_17_151852 H21.03.315.11.09.01.01
sub-1122692698/sub-1122692698_ses-1122759713_icephys.nwb 2021_08_17_155235 H21.03.315.11.09.01.02
sub-1154304730/sub-1154304730_ses-1154338526_icephys.nwb 2022_01_27_191454 H22.26.404.11.06.02.01
Distinct table cells matched: 210
Sex agrees with the table's donor_sex: 209 of 210
Records whose sex differs:
cell_name subject_id sex donor_sex
H19.03.321.12.04.01.04 977091139 unknown M
Files whose age in days is a whole number of 365-day years: 210 of 210
Age in years equal to the table's donor_age: 208 of 208 | largest difference: 0.0
Records whose age differs:
Empty DataFrame
Columns: [cell_name, subject_id, donor_name, age, donor_age, years]
Index: []
Files with an age where the table has none: 2 ['P0.0D']
Donors per subject_id: 1
# Features computed here against the paper's per-cell values (reading (a) against (b) of O1)
computed = L.with_adapt_order(linked.set_index("spec_id"))
fcols = L.feature_columns(paper_raw)
cmp_rows = []
for c in [c for c in fcols if c in computed.columns]:
a, b = computed[c].astype(float), paper_raw[c].reindex(computed.index).astype(float)
ok = a.notna() & b.notna()
close = (np.abs(a[ok] - b[ok]) <= 0.01 * np.abs(b[ok]) + 1e-9)
cmp_rows.append({"feature": c, "n_both": int(ok.sum()), "only_here": int((a.notna() & b.isna()).sum()),
"only_paper": int((a.isna() & b.notna()).sum()),
"within_1pct": round(close.mean(), 3) if ok.any() else np.nan,
"spearman": round(stats.spearmanr(a[ok], b[ok])[0], 3) if ok.sum() > 2 else np.nan})
cmp = pd.DataFrame(cmp_rows).set_index("feature")
print(cmp.sort_values("within_1pct").to_string())
# The order of the adaptation ratios (O1): which order matches the paper's table
for name in ["upstroke", "downstroke", "width", "isi"]:
b = paper_raw[f"{name}_adapt_ratio"].reindex(computed.index)
for order in ["3over1", "1over3"]:
a = computed[f"{name}_adapt_ratio_{order}"]
ok = a.notna() & b.notna()
print(f"{name}_adapt_ratio, order {order}: median |a/b - 1| = {np.median(np.abs(a[ok] / b[ok] - 1)):.3f}, "
f"spearman {stats.spearmanr(a[ok], b[ok])[0]:.3f} (n={ok.sum()})")
n_both only_here only_paper within_1pct spearman feature isi_adapt_ratio 206 1 1 0.039 0.724 downstroke_adapt_ratio 206 1 1 0.165 0.867 upstroke_adapt_ratio 206 1 1 0.189 0.922 width_adapt_ratio 206 1 1 0.199 0.906 threshold_v_adapt_ratio 206 1 1 0.374 0.874 peak_v_adapt_ratio 206 1 1 0.534 0.936 fi_fit_slope 208 0 2 0.707 0.917 input_resistance 209 0 1 0.837 0.984 tau 209 0 1 0.842 0.989 fast_trough_deltav_hero 209 0 0 0.852 0.991 latency_hero 210 0 0 0.852 0.927 avg_rate_hero 210 0 0 0.857 0.898 first_isi_inv_hero 206 1 1 0.859 0.902 isi_cv_hero 206 1 1 0.859 0.976 adapt_hero 205 1 0 0.859 0.936 downstroke_hero 209 0 0 0.861 0.968 trough_v_hero 209 0 0 0.861 0.979 mean_isi_hero 206 1 1 0.864 0.887 median_isi_hero 206 1 1 0.864 0.886 fast_trough_v_hero 209 0 0 0.871 0.990 peak_v_hero 210 0 0 0.871 0.984 width_hero 209 0 0 0.871 0.973 upstroke_downstroke_ratio_hero 209 0 0 0.880 0.989 upstroke_hero 210 0 0 0.881 0.974 peak_deltav_hero 210 0 0 0.895 0.974 threshold_v_hero 210 0 0 0.919 0.994 peak_v_rheo 210 0 0 0.933 0.987 sag 210 0 0 0.933 0.971 latency_rheo 210 0 0 0.933 0.967 mean_isi_rheo 127 1 4 0.937 0.976 first_isi_inv_rheo 127 1 4 0.937 0.969 avg_rate_rheo 210 0 0 0.938 0.953 fast_trough_deltav_rheo 208 0 1 0.938 0.996 peak_deltav_rheo 210 0 0 0.938 0.977 downstroke_rheo 208 0 1 0.942 0.993 width_rheo 208 0 1 0.942 0.990 trough_v_rheo 208 0 1 0.942 0.973 upstroke_rheo 210 0 0 0.943 0.992 threshold_v_rheo 210 0 0 0.943 0.973 fast_trough_v_rheo 208 0 1 0.947 0.989 vm_for_sag 210 0 0 0.952 0.961 upstroke_downstroke_ratio_rheo 208 0 1 0.957 0.988 v_baseline 210 0 0 0.971 0.997 rheobase_i 210 0 0 0.976 0.971 upstroke_adapt_ratio, order 3over1: median |a/b - 1| = 0.027, spearman 0.922 (n=206) upstroke_adapt_ratio, order 1over3: median |a/b - 1| = 0.411, spearman -0.922 (n=206) downstroke_adapt_ratio, order 3over1: median |a/b - 1| = 0.032, spearman 0.867 (n=206) downstroke_adapt_ratio, order 1over3: median |a/b - 1| = 0.425, spearman -0.867 (n=206) width_adapt_ratio, order 3over1: median |a/b - 1| = 0.029, spearman 0.906 (n=206) width_adapt_ratio, order 1over3: median |a/b - 1| = 0.279, spearman -0.906 (n=206) isi_adapt_ratio, order 3over1: median |a/b - 1| = 0.140, spearman 0.724 (n=206) isi_adapt_ratio, order 1over3: median |a/b - 1| = 0.549, spearman -0.724 (n=206)
# Reading (a): a feature computed here is used when it agrees with the paper's value for most
# cells (within 1% for at least half of them); otherwise, and for every feature not computed
# here, the paper's value is used (reading (b)). The 50% rule is an assumption of this attempt;
# the claims file's O1 says only that (a) is the recomputation and (b) the paper's values.
COMPUTED = [c for c in cmp.index if cmp.loc[c, "within_1pct"] >= 0.5]
NOT_USED = [c for c in cmp.index if c not in COMPUTED]
print(f"{len(fcols)} features in the paper's table; computed here and used: {len(COMPUTED)}; "
f"computed but not agreeing (paper's value used): {NOT_USED}")
print("Features taken from the paper's table:", len(fcols) - len(COMPUTED))
hyb = L.apply_code_filters(L.hybrid_table(computed, paper_raw, COMPUTED))
pap = L.apply_code_filters(paper_raw)
# Variants of reading (a): (a') the same features extracted without the Methods' 500 ms rms noise
# criterion (see the next cell), and (a, ratios) with the five adaptation ratios computed here in
# place of the paper's values.
computed_ns = L.with_adapt_order(linked.set_index("spec_id")[[c for c in linked.columns if c.startswith("noslow__")]]
.rename(columns=lambda c: c[len("noslow__"):]))
hyb_ns = L.apply_code_filters(L.hybrid_table(computed_ns, paper_raw, COMPUTED))
hyb_r = L.apply_code_filters(L.hybrid_table(computed, paper_raw, COMPUTED + NOT_USED))
print("Cells in the paper's table not in the dandiset:", len(paper_raw.index.difference(computed.index)))
94 features in the paper's table; computed here and used: 39; computed but not agreeing (paper's value used): ['isi_adapt_ratio', 'width_adapt_ratio', 'upstroke_adapt_ratio', 'downstroke_adapt_ratio', 'threshold_v_adapt_ratio'] Features taken from the paper's table: 55
Cells in the paper's table not in the dandiset: 53
# Where the features computed here disagree with the paper's values: for each feature, the cells
# outside 1% and how many of them have a different rheobase; for the initial firing rate at
# rheobase + 40 pA, each such cell with the sweep amplitudes on both sides. The paper's table
# gives no sweep amplitudes, so its hero sweep is found as the sweep of the file whose inverse
# first ISI equals the paper's value within 1%.
rheo_diff = (computed.rheobase_i - paper_raw.rheobase_i.reindex(computed.index)).abs() > 1
rows = []
for c in COMPUTED:
a, b = computed[c].astype(float), paper_raw[c].reindex(computed.index).astype(float)
ok = a.notna() & b.notna()
out = ok & ~(np.abs(a - b) <= 0.01 * np.abs(b) + 1e-9)
rows.append({"feature": c, "n_both": int(ok.sum()), "outside_1pct": int(out.sum()),
"of_which_rheobase_differs": int((out & rheo_diff).sum()),
"median_rel_diff_outside": round(float(np.median(np.abs(a[out] / b[out] - 1))), 3) if out.any() else np.nan})
print(pd.DataFrame(rows).set_index("feature").sort_values("outside_1pct", ascending=False).to_string())
def sweep_list(s):
return [tuple(float(v) for v in x.split(":")) for x in s.split(";")] if isinstance(s, str) and s else []
f = "first_isi_inv_hero"
a, b = computed[f].astype(float), paper_raw[f].reindex(computed.index).astype(float)
diff = (a.notna() != b.notna()) | (a.notna() & b.notna() & ~(np.abs(a - b) <= 0.01 * np.abs(b) + 1e-9))
rows = []
for sid in computed.index[diff]:
sw = sweep_list(computed.loc[sid, "spiking_sweeps"])
at = sorted({round(amp, 1) for amp, n, inv in sw if np.isfinite(inv) and np.isfinite(b[sid]) and abs(inv - b[sid]) <= 0.01 * abs(b[sid])})
hero = round(computed.loc[sid, "stim_amp_hero"], 1)
n_at_hero = sum(abs(amp - hero) < 0.5 for amp, n, inv in sw)
if not np.isfinite(b[sid]) or not np.isfinite(a[sid]):
why = "value on one side only"
elif any(abs(x - hero) < 0.5 for x in at):
why = "another repeat of the same amplitude"
elif at:
why = "another sweep amplitude"
else:
why = "no sweep here gives the paper's value"
rows.append({"specimen_id": sid, "here": round(a[sid], 1), "paper": round(b[sid], 1),
"rheobase_here": computed.loc[sid, "rheobase_i"], "rheobase_paper": round(paper_raw.loc[sid, "rheobase_i"], 1),
"hero_amp_here": hero, "sweeps_at_that_amp": n_at_hero,
"amp_giving_paper_value": ", ".join(f"{x:g}" for x in at), "explanation": why})
fx = pd.DataFrame(rows).set_index("specimen_id")
print(f"\n{f}: {int(diff.sum())} of {int((a.notna() | b.notna()).sum())} cells differ by more than 1% or have a value on one side only")
print(fx.to_string())
print("\nBy explanation:", fx.explanation.value_counts().to_dict())
in_fig2c = computed.in_strict_set & (computed.collaborator == "AIBS") & computed["t-type"].map(L.HOMOLOGY).isin(L.SUBCLASSES)
print("Of these, in the test set of fig2c (P, four core subclasses):", int(in_fig2c.reindex(fx.index).sum()), "of", int(in_fig2c.sum()),
"| by explanation:", fx[in_fig2c.reindex(fx.index).values].explanation.value_counts().to_dict())
n_both outside_1pct of_which_rheobase_differs median_rel_diff_outside
feature
peak_v_adapt_ratio 206 96 3 0.027
fi_fit_slope 208 61 2 0.143
input_resistance 209 34 2 0.033
tau 209 33 2 0.043
latency_hero 210 31 2 0.302
fast_trough_deltav_hero 209 31 2 0.086
avg_rate_hero 210 30 2 0.272
first_isi_inv_hero 206 29 2 0.315
isi_cv_hero 206 29 2 0.150
adapt_hero 205 29 2 0.385
downstroke_hero 209 29 2 0.054
trough_v_hero 209 29 2 0.035
mean_isi_hero 206 28 2 0.239
median_isi_hero 206 28 2 0.239
width_hero 209 27 2 0.032
peak_v_hero 210 27 2 0.066
fast_trough_v_hero 209 27 2 0.028
upstroke_downstroke_ratio_hero 209 25 2 0.036
upstroke_hero 210 25 2 0.040
peak_deltav_hero 210 22 2 0.032
threshold_v_hero 210 17 1 0.037
latency_rheo 210 14 5 0.284
sag 210 14 2 0.228
peak_v_rheo 210 14 5 0.165
avg_rate_rheo 210 13 4 0.600
peak_deltav_rheo 210 13 5 0.063
fast_trough_deltav_rheo 208 13 4 0.136
width_rheo 208 12 3 0.044
upstroke_rheo 210 12 4 0.086
downstroke_rheo 208 12 4 0.061
threshold_v_rheo 210 12 4 0.031
trough_v_rheo 208 12 4 0.037
fast_trough_v_rheo 208 11 3 0.044
vm_for_sag 210 10 1 0.068
upstroke_downstroke_ratio_rheo 208 9 3 0.056
mean_isi_rheo 127 8 3 0.313
first_isi_inv_rheo 127 8 3 0.351
v_baseline 210 6 2 0.025
rheobase_i 210 5 5 1.000
first_isi_inv_hero: 31 of 208 cells differ by more than 1% or have a value on one side only
here paper rheobase_here rheobase_paper hero_amp_here sweeps_at_that_amp amp_giving_paper_value explanation
specimen_id
1005131537 55.4 30.5 20.0 20.0 74.0 1 54 another sweep amplitude
1010273140 44.2 40.2 10.0 10.0 50.0 2 50 another repeat of the same amplitude
1154338548 50.3 73.4 14.0 14.0 64.0 3 64 another repeat of the same amplitude
569844356 56.9 70.8 20.0 20.0 50.0 1 no sweep here gives the paper's value
602196815 22.8 NaN 30.0 30.0 70.0 3 value on one side only
611901177 29.7 47.9 20.0 20.0 50.0 1 no sweep here gives the paper's value
650130090 38.8 13.8 20.0 20.0 60.0 1 30 another sweep amplitude
652081101 46.4 12.7 10.0 10.0 50.0 1 no sweep here gives the paper's value
653864596 60.4 34.4 110.0 30.0 110.0 2 no sweep here gives the paper's value
653907140 100.0 54.8 90.0 10.0 90.0 1 no sweep here gives the paper's value
672186592 34.3 45.7 10.0 10.0 50.0 1 no sweep here gives the paper's value
672646211 11.9 64.9 90.0 90.0 90.0 1 no sweep here gives the paper's value
672659580 17.0 3.3 30.0 30.0 70.0 1 30 another sweep amplitude
685814167 NaN 1.0 30.0 30.0 30.0 1 value on one side only
712885453 7.6 38.8 10.0 10.0 10.0 1 no sweep here gives the paper's value
737989339 27.4 17.3 10.0 10.0 60.0 1 40 another sweep amplitude
769331519 40.7 36.5 20.0 20.0 70.0 11 70 another repeat of the same amplitude
805729797 26.6 25.9 20.0 20.0 60.0 11 60 another repeat of the same amplitude
805732962 69.8 71.7 10.0 10.0 50.0 4 50 another repeat of the same amplitude
811804107 12.9 13.6 10.0 10.0 20.0 7 20 another repeat of the same amplitude
811830204 18.1 22.4 30.0 30.0 70.0 1 90 another sweep amplitude
832529199 22.4 17.7 40.0 40.0 80.0 2 80 another repeat of the same amplitude
832594616 41.9 60.2 10.0 10.0 50.0 1 70 another sweep amplitude
832614788 40.2 29.9 20.0 20.0 70.0 2 70 another repeat of the same amplitude
835611255 60.1 47.4 10.0 10.0 70.0 1 no sweep here gives the paper's value
835857228 57.4 37.3 10.0 10.0 70.0 1 no sweep here gives the paper's value
854543457 37.9 38.5 20.0 20.0 60.0 7 60 another repeat of the same amplitude
857650458 35.1 33.9 10.0 10.0 50.0 1 no sweep here gives the paper's value
868243321 12.2 48.1 10.0 10.0 20.0 15 no sweep here gives the paper's value
868257910 33.0 31.0 10.0 10.0 50.0 2 50 another repeat of the same amplitude
880478829 69.7 48.9 10.0 10.0 62.0 1 no sweep here gives the paper's value
By explanation: {"no sweep here gives the paper's value": 13, 'another repeat of the same amplitude': 10, 'another sweep amplitude': 6, 'value on one side only': 2}
Of these, in the test set of fig2c (P, four core subclasses): 28 of 179 | by explanation: {'another repeat of the same amplitude': 10, "no sweep here gives the paper's value": 10, 'another sweep amplitude': 6, 'value on one side only': 2}
# The sweep at rheobase + 40 pA. IPFX (find_hero_sweep) takes the lowest spiking sweep 39 to 61 pA
# above rheobase; where no sweep passing QC is in that range, it takes the spiking sweep closest to
# rheobase + 39 pA, which can be the rheobase sweep itself. Here each such cell is shown with its
# value as run and with the sweeps that fail only the Methods' 500 ms rms noise criterion (below
# 0.5 mV) kept (variant a').
off = computed.stim_amp_hero - computed.rheobase_i
off_ns = computed_ns.stim_amp_hero - computed_ns.rheobase_i
fallback = (off < 39) | (off > 61)
pv = paper_raw.first_isi_inv_hero.reindex(computed.index)
def agrees_paper(x):
return (np.abs(x - pv) <= 0.01 * np.abs(pv) + 1e-9) | (x.isna() & pv.isna())
fbt = pd.DataFrame({"rheobase": computed.rheobase_i, "hero_amp": computed.stim_amp_hero, "value": computed.first_isi_inv_hero,
"hero_amp_a'": computed_ns.stim_amp_hero, "value_a'": computed_ns.first_isi_inv_hero, "paper": pv,
"agrees": agrees_paper(computed.first_isi_inv_hero), "agrees_a'": agrees_paper(computed_ns.first_isi_inv_hero),
"in_fig2c": in_fig2c})[fallback].round(1)
print(f"Cells whose sweep at rheobase + 40 pA is IPFX's fallback (no sweep passing QC 39 to 61 pA above rheobase): "
f"{int(fallback.sum())} of {len(off)}; agreeing with the paper's value within 1%: {int(fbt.agrees.sum())}; "
f"in the test set of fig2c: {int(fbt.in_fig2c.sum())}")
print(fbt.to_string())
print("Of the fallback cells that disagree, agreeing in variant a' (the paper's value comes from a sweep that fails "
"only the 500 ms noise criterion here):", int((~fbt.agrees & fbt["agrees_a'"]).sum()), "of", int((~fbt.agrees).sum()))
print("Fallback cells in variant a':", int(((off_ns < 39) | (off_ns > 61)).sum()))
agree = pd.DataFrame({name: {c: float(((np.abs(X[c].astype(float) - paper_raw[c].reindex(X.index)) <= 0.01 * np.abs(paper_raw[c].reindex(X.index)) + 1e-9)
[X[c].notna() & paper_raw[c].reindex(X.index).notna()]).mean()) for c in COMPUTED}
for name, X in [("a", computed), ("a'", computed_ns)]})
ap, aa = agree["a'"], agree["a"]
print(f"\nShare of cells within 1% of the paper's value, over the {len(COMPUTED)} features of reading (a): median "
f"{aa.median():.3f} as run, {ap.median():.3f} in variant a'; higher in a' for {int((ap > aa).sum())} features, "
f"lower for {int((ap < aa).sum())}")
print(agree.round(3).sort_values("a").to_string())
Cells whose sweep at rheobase + 40 pA is IPFX's fallback (no sweep passing QC 39 to 61 pA above rheobase): 16 of 210; agreeing with the paper's value within 1%: 7; in the test set of fig2c: 11
rheobase hero_amp value hero_amp_a' value_a' paper agrees agrees_a' in_fig2c
spec_id
1122759850 30.0 58.0 30.3 58.0 30.3 30.3 True True True
569844356 20.0 50.0 56.9 70.0 70.8 70.8 False True False
596893815 20.0 50.0 22.5 50.0 22.5 22.5 True True True
611901177 20.0 50.0 29.7 70.0 47.9 47.9 False True False
653810044 50.0 80.0 32.2 80.0 32.2 32.2 True True True
653864596 110.0 110.0 60.4 70.0 34.4 34.4 False True True
653907140 90.0 90.0 100.0 50.0 54.8 54.8 False True True
672646211 90.0 90.0 11.9 130.0 64.9 64.9 False True False
685814167 30.0 30.0 NaN 90.0 4.1 1.0 False False True
712885453 10.0 10.0 7.6 50.0 38.8 38.8 False True True
732141075 10.0 30.0 35.1 30.0 35.1 35.1 True True False
805751552 10.0 10.0 6.9 10.0 6.9 6.9 True True True
811804107 10.0 20.0 12.9 20.0 12.9 13.6 False False True
835648767 10.0 20.0 12.1 20.0 12.1 12.1 True True True
857657437 600.0 600.0 NaN 600.0 NaN NaN True True False
868243321 10.0 20.0 12.2 50.0 48.1 48.1 False True True
Of the fallback cells that disagree, agreeing in variant a' (the paper's value comes from a sweep that fails only the 500 ms noise criterion here): 7 of 9
Fallback cells in variant a': 8
Share of cells within 1% of the paper's value, over the 39 features of reading (a): median 0.895 as run, 0.933 in variant a'; higher in a' for 38 features, lower for 1
a a'
peak_v_adapt_ratio 0.534 0.536
fi_fit_slope 0.707 0.695
input_resistance 0.837 0.862
tau 0.842 0.862
fast_trough_deltav_hero 0.852 0.895
latency_hero 0.852 0.895
avg_rate_hero 0.857 0.900
isi_cv_hero 0.859 0.899
adapt_hero 0.859 0.902
first_isi_inv_hero 0.859 0.899
trough_v_hero 0.861 0.904
downstroke_hero 0.861 0.904
median_isi_hero 0.864 0.903
mean_isi_hero 0.864 0.903
peak_v_hero 0.871 0.905
fast_trough_v_hero 0.871 0.909
width_hero 0.871 0.900
upstroke_downstroke_ratio_hero 0.880 0.914
upstroke_hero 0.881 0.914
peak_deltav_hero 0.895 0.933
threshold_v_hero 0.919 0.943
sag 0.933 0.952
peak_v_rheo 0.933 0.962
latency_rheo 0.933 0.957
mean_isi_rheo 0.937 0.977
first_isi_inv_rheo 0.937 0.977
fast_trough_deltav_rheo 0.938 0.966
peak_deltav_rheo 0.938 0.962
avg_rate_rheo 0.938 0.971
trough_v_rheo 0.942 0.966
width_rheo 0.942 0.962
downstroke_rheo 0.942 0.966
upstroke_rheo 0.943 0.957
threshold_v_rheo 0.943 0.957
fast_trough_v_rheo 0.947 0.966
vm_for_sag 0.952 0.967
upstroke_downstroke_ratio_rheo 0.957 0.971
v_baseline 0.971 0.981
rheobase_i 0.976 0.990
The sweep at rheobase + 40 pA is IPFX's fallback in 16 cells: no sweep passing QC lies 39 to 61 pA above rheobase, so IPFX takes the spiking sweep closest to rheobase + 39 pA, at times the rheobase sweep itself. In 7 of these the value agrees with the paper's, so the paper's pipeline fell back in the same way. In 7 of the 9 that disagree, the paper's value comes from a sweep in range that fails here only the Methods' criterion of rms noise below 0.5 mV over 500 ms, as IPFX measures it on these files. Without that criterion (variant a') the features agree with the paper's values for more cells in 38 of the 39 features of reading (a) (median share 0.933 against 0.895). The paper's own values thus show that its pipeline kept sweeps that this criterion, as measured here, leaves out. Variant a' is run alongside readings (a) and (b) below.
Result: partially reproduced. All 210 files match a cell of the paper's table through the
LIMS identifier (session_id equals the Ephys_Roi_Result number in the table's nwb_path), one
to one. 204 files carry the table's cell name; 6 carry a time stamp in its place. Age agrees for
all 208 cells with a known age: the files state ages as whole years times 365 days, and the one
mismatch of the first run came from dividing by 365.25. Two files of donor 922967568 state an age of
P0.0D where the table has none, and one file (H19.03.321.12.04.01.04) states sex "unknown" where
the table gives M, so the records do not all agree, and the criterion's partial level applies (at
least 90% of files match and agree). Of the 44 features computed here, the 39 used in reading (a)
agree with the paper's values within 1% for 53 to 98% of cells (Spearman 0.886 to 0.997); the five
adaptation ratios not used agree for 4 to 37%, and the 50 features computed by code the paper did
not release are the paper's values. Of the 31 cells where the initial firing rate at rheobase + 40
pA differs, 10 take another repeat of the same amplitude here, 6 another amplitude, 13 have no
sweep here that gives the paper's value, and 2 have a value on one side only. The next cells show
the main cause of the last group: in 16 cells no sweep passing QC lies 39 to 61 pA above rheobase
and IPFX falls back to the spiking sweep closest to rheobase + 39 pA, at times the rheobase sweep.
In 7 of the 9 of these cells that disagree with the paper, the paper's value comes from a sweep in
range that fails here only the 500 ms noise criterion.
# Sets used below (claims file, "Sets used below"; l1_load.py)
cells = linked.set_index("spec_id")
cells["homology_type"] = pd.Categorical(cells["t-type"].map(L.HOMOLOGY),
categories=["LAMP5", "MC4R", "PAX6", "L1 VIP", "other"], ordered=True)
cells["l1_type"] = cells.homology_type.notna()
cells["core_l1_type"] = cells.homology_type.isin(L.SUBCLASSES)
cells["failed_fx_long_squares"] = paper_raw["failed_fx_long_squares"].reindex(cells.index)
cells["P"] = cells.in_strict_set & (cells.collaborator == "AIBS") & (cells.failed_fx_long_squares == False)
print("Dandiset cells:", len(cells), "| in the strict L1 set:", int(cells.in_strict_set.sum()),
"| AIBS:", int((cells.collaborator == "AIBS").sum()), "| in P:", int(cells.P.sum()))
lab_P = labels[labels.in_strict_set & (labels.collaborator == "AIBS")].join(paper_raw[["failed_fx_long_squares"]])
lab_P = lab_P[lab_P.failed_fx_long_squares == False]
print("P in the paper's table:", len(lab_P), "| of which in the dandiset:", int(lab_P.index.isin(cells.index).sum()))
print("Strict L1 set in the paper's table by site:", labels[labels.in_strict_set].collaborator.value_counts().to_dict())
print("Cells by t-type in P:")
print(cells[cells.P]["t-type"].value_counts().to_string())
print("Donors:", cells.donor_name.nunique(), "| cells per donor:", cells.donor_name.value_counts().describe()[["min", "50%", "max"]].to_dict())
def table(kind):
# Cells with their labels and the features of reading (a) (hybrid) or (b) (paper)
feat = {"a": hyb, "a'": hyb_ns, "a, ratios": hyb_r, "b": pap}[kind]
return cells.drop(columns=[c for c in fcols if c in cells.columns]).join(feat[fcols])
TA, TB = table("a"), table("b")
TA2, TAR = table("a'"), table("a, ratios")
# Features of the Amsterdam (Huib) and Szeged (Gabor) cells from the paper's code. They are not
# in the dandiset and are used only where a claim's method needs the other sites (stated there).
tamas = L.apply_code_filters(L.read_table("tamas_features_E.csv", index_col=0))
mansv = L.apply_code_filters(L.read_table("mansvelder_features_E.csv", index_col=0))
sites = pd.concat([
labels[labels.collaborator == "Gabor"].join(tamas, on="sample_id"),
labels[labels.collaborator == "Huib"].join(mansv, on="sample_id"),
labels[labels.collaborator == "AIBS"].join(hyb)], sort=False)
sites_ns = pd.concat([sites[sites.collaborator != "AIBS"], sites[sites.collaborator == "AIBS"].drop(columns=fcols).join(hyb_ns[fcols])], sort=False)
sites["homology_type"] = pd.Categorical(sites["t-type"].map(L.HOMOLOGY), categories=L.SUBCLASSES + ["other"])
sites_ns["homology_type"] = pd.Categorical(sites_ns["t-type"].map(L.HOMOLOGY), categories=L.SUBCLASSES + ["other"])
Dandiset cells: 210 | in the strict L1 set: 210 | AIBS: 210 | in P: 210
P in the paper's table: 210 | of which in the dandiset: 210
Strict L1 set in the paper's table by site: {'AIBS': 222, 'Gabor': 28, 'Huib': 22}
Cells by t-type in P:
t-type
LAMP5 NMBR 71
PAX6 CDH12 31
ADARB2 MC4R 26
VIP TSPAN12 15
LAMP5 LCP2 13
LAMP5 DBP 12
SST CHRNA4 9
VIP PCDH20 9
VIP SYT6 7
SST BAGE2 6
PAX6 TNFAIP8L3 2
VIP LBH 2
LAMP5 CA1 2
VIP ADAMTSL1 2
VIP PENK 2
VIP CHRNA6 1
Donors: 53 | cells per donor: {'min': 1.0, '50%': 3.0, 'max': 19.0}
The dandiset holds exactly the paper's primary set P: the 210 Seattle (AIBS) cells of the strict L1 set whose long square features were extracted. The other 62 cells of the strict set are not in it: the 50 Amsterdam (22) and Szeged (28) cells, and 12 Seattle cells that failed feature extraction.
patchseq-counts¶
The human L1 patch-seq dataset has 250 cells passing transcriptomic quality, 194 with electrophysiology passing quality control and 71 reconstructed; by subclass, 96 LAMP5, 35 MC4R, 33 PAX6, 15 L1 VIP and 15 other cells have electrophysiology.
P = cells[cells.P]
counts = P[P.l1_type].homology_type.value_counts().reindex(["LAMP5", "MC4R", "PAX6", "L1 VIP", "other"])
print("Cells of P in the ten L1 t-types:", int(P.l1_type.sum()), "(expected 194 +/- 19.4)")
print(pd.DataFrame({"count": counts, "expected": [96, 35, 33, 15, 15]}).to_string())
print("Dandiset files not in P:", len(cells) - len(P), "| cells of P without an L1 t-type:", int((~P.l1_type).sum()),
P[~P.l1_type]["t-type"].value_counts().to_dict())
print("By donor: cells of P in L1 t-types per donor, median",
P[P.l1_type].donor_name.value_counts().median(), "range",
P[P.l1_type].donor_name.value_counts().min(), "to", P[P.l1_type].donor_name.value_counts().max())
Cells of P in the ten L1 t-types: 194 (expected 194 +/- 19.4)
count expected
homology_type
LAMP5 96 96
MC4R 35 35
PAX6 33 33
L1 VIP 15 15
other 15 15
Dandiset files not in P: 0 | cells of P without an L1 t-type: 16 {'VIP SYT6': 7, 'VIP LBH': 2, 'LAMP5 CA1': 2, 'VIP ADAMTSL1': 2, 'VIP PENK': 2, 'VIP CHRNA6': 1}
By donor: cells of P in L1 t-types per donor, median 2.0 range 1 to 19
Result: reproduced. The dandiset holds 194 cells of P in the ten L1 t-types, with 96 LAMP5, 35 MC4R, 33 PAX6, 15 L1 VIP and 15 other, the counts of Table S2. The other 16 files of P are cells of non-L1 t-types. The 250 and 71 counts need transcriptomes and reconstructions and are not tested.
fig2c-subclass-ephys-differences¶
Human L1 subclasses differ in electrophysiology, notably in sag and initial firing rate (Kruskal-Wallis, FDR-corrected p < 1e-7).
As l1_subclass_distinctness.ipynb cells 8 and 10: features with more than 40% non-missing values
among l1_type cells of P; Kruskal-Wallis across the four core subclasses; Benjamini-Hochberg;
Dunn's post hoc (code's Mann-Whitney shown as a variant). Run on reading (a) (features computed
here, the paper's values for the others) and (b) (the paper's values). The initial firing rate
(first_isi_inv_hero) is computed here in reading (a); the cells where it differs from the
paper's value are listed under dandiset-cell-link. Variant (a') is reading (a) with the
features extracted without the Methods' 500 ms rms noise criterion.
fig2c = {}
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
data = T[T.P & T.l1_type]
feats40 = [c for c in fcols if data[c].notna().mean() > 0.4]
sub = T[T.P & T.core_l1_type].copy()
sub["homology_type"] = sub.homology_type.cat.remove_unused_categories()
res, pairs = L.run_kw_dunn(sub, feats40, "homology_type", posthoc="dunn")
res_mw, pairs_mw = L.run_kw_dunn(sub, feats40, "homology_type", posthoc="mw")
fig2c[kind] = (res, pairs, pairs_mw, feats40, sub)
print(f"Reading ({kind}): {len(feats40)} features, {len(sub)} cells", sub.homology_type.value_counts().to_dict())
print(f" features with FDR p < 0.05: {(res.pval_fdr < 0.05).sum()}, < 1e-7: {(res.pval_fdr < 1e-7).sum()}")
print(res.loc[["sag", "first_isi_inv_hero"], ["KW_H", "pval", "pval_fdr", "epsilon2", "n"]].to_string())
print(" ten smallest p:"); print(res.head(10)[["pval_fdr", "epsilon2"]].to_string())
Reading (a): 92 features, 179 cells {'LAMP5': 96, 'MC4R': 35, 'PAX6': 33, 'L1 VIP': 15}
features with FDR p < 0.05: 56, < 1e-7: 15
KW_H pval pval_fdr epsilon2 n
feature
sag 101.574795 7.126244e-22 6.556144e-20 0.570645 179
first_isi_inv_hero 35.696649 8.679941e-08 4.697380e-07 0.201676 178
ten smallest p:
pval_fdr epsilon2
feature
sag 6.556144e-20 0.570645
upstroke_adapt_ratio 1.011960e-15 0.458376
sag_area 2.200193e-14 0.413571
sag_peak_t 4.975098e-13 0.374761
threshold_v_adapt_ratio 3.812129e-12 0.352951
input_resistance 4.215510e-09 0.269765
width_adapt_ratio 9.026421e-09 0.258640
downstroke_adapt_ratio 9.026421e-09 0.257592
fast_trough_deltav_hero 9.143029e-09 0.252446
upstroke_downstroke_ratio_rheo 9.143029e-09 0.253412
Reading (a'): 92 features, 179 cells {'LAMP5': 96, 'MC4R': 35, 'PAX6': 33, 'L1 VIP': 15}
features with FDR p < 0.05: 56, < 1e-7: 16
KW_H pval pval_fdr epsilon2 n
feature
sag 102.974980 3.562303e-22 3.277319e-20 0.578511 179
first_isi_inv_hero 46.764899 3.899752e-10 3.900865e-09 0.262724 179
ten smallest p:
pval_fdr epsilon2
feature
sag 3.277319e-20 0.578511
upstroke_adapt_ratio 1.011960e-15 0.458376
sag_area 2.200193e-14 0.413571
sag_peak_t 4.975098e-13 0.374761
threshold_v_adapt_ratio 3.812129e-12 0.352951
tau 3.900865e-09 0.265620
latency_hero 3.900865e-09 0.263901
input_resistance 3.900865e-09 0.264371
first_isi_inv_hero 3.900865e-09 0.262724
fast_trough_deltav_hero 3.900865e-09 0.261764
Reading (b): 92 features, 179 cells {'LAMP5': 96, 'MC4R': 35, 'PAX6': 33, 'L1 VIP': 15}
features with FDR p < 0.05: 56, < 1e-7: 17
KW_H pval pval_fdr epsilon2 n
feature
sag 102.487897 4.534084e-22 4.171357e-20 0.575775 179
first_isi_inv_hero 43.426098 1.998179e-09 1.671204e-08 0.245345 178
ten smallest p:
pval_fdr epsilon2
feature
sag 4.171357e-20 0.575775
upstroke_adapt_ratio 1.011960e-15 0.458376
sag_area 2.200193e-14 0.413571
sag_peak_t 4.975098e-13 0.374761
threshold_v_adapt_ratio 3.812129e-12 0.352951
fast_trough_deltav_hero 3.075715e-09 0.270348
tau 6.163789e-09 0.260608
width_adapt_ratio 7.221137e-09 0.258640
input_resistance 7.221137e-09 0.256499
downstroke_adapt_ratio 7.221137e-09 0.257592
print("By donor: number of donors per subclass in the test set (cells come from",
fig2c["a"][4].donor_name.nunique(), "donors)")
print(fig2c["a"][4].groupby("homology_type", observed=True).donor_name.nunique().to_string())
# Leave-one-donor-out check of the two features: largest FDR p over runs without each donor
sub = fig2c["a"][4]
worst = {}
for f in ["sag", "first_isi_inv_hero"]:
ps = []
for d in sub.donor_name.unique():
s = sub[sub.donor_name != d].dropna(subset=[f])
ps.append(stats.kruskal(*[g[f].values for _, g in s.groupby("homology_type", observed=True)]).pvalue)
worst[f] = max(ps)
print("Largest uncorrected KW p when any one donor is left out:", {k: f"{v:.2g}" for k, v in worst.items()})
By donor: number of donors per subclass in the test set (cells come from 49 donors) homology_type LAMP5 35 MC4R 23 PAX6 20 L1 VIP 13
Largest uncorrected KW p when any one donor is left out: {'sag': '3.1e-19', 'first_isi_inv_hero': '3.5e-05'}
Result: reproduced, resting on an assumption. Sag differs across the
four subclasses with FDR p = 6.6e-20 on the features computed here. The initial firing rate at
rheobase + 40 pA also differs, with FDR p = 4.7e-7 on the features computed here with the Methods'
sweep QC, above the 1e-7 of the claim; with the paper's own per-cell values it is 1.7e-8, below
it. The computed initial firing rate differs from the paper's value by more than 1%, or has a
value on one side only, in 28 of the 179 cells of this test: 10 take another repeat of the same
amplitude, 6 another amplitude, 10 have no sweep here that gives the paper's value and 2 have a
value on one side only. In 11 cells of this test no sweep passing QC lies 39 to 61 pA above
rheobase and IPFX falls back to a lower sweep; 6 of them disagree with the paper, and in 4 of these
the paper's value comes from a sweep in range that fails here only the Methods' 500 ms noise
criterion (table after dandiset-cell-link). Without that criterion (variant a') the features agree with the
paper's values for more cells in 38 of the 39 features, and the initial firing rate has FDR p =
3.9e-9, below 1e-7, so the claim holds in full. The two recomputations give different outcomes. The
paper's own per-cell values show that its pipeline kept these sweeps, so variant (a') is the better
supported recomputation, and the outcome rests on that assumption, which departs from the Methods'
QC as IPFX measures it; with the criterion as written the claim is partially reproduced. The third
round of verification agreed with the change to reproduced. The authors are asked how the 500 ms
noise was measured. With any one donor left out (reading a), the largest uncorrected p is 3.1e-19
for sag and 3.5e-5 for the initial firing rate; this check was not run on (a').
fig2b-subclass-ephys-profiles¶
LAMP5 cells show small sag and adaptation, PAX6 cells minimal sag and a high initial firing rate, while MC4R and L1 VIP cells show strong sag.
for kind in ["a", "a'", "b"]:
res, pairs, pairs_mw, feats40, sub = fig2c[kind]
med = sub.groupby("homology_type", observed=True)[["sag", "first_isi_inv_hero"]].median()
n = sub.groupby("homology_type", observed=True)[["sag", "first_isi_inv_hero"]].count()
print(f"Reading ({kind}): medians (n)"); print(med.round(3).astype(str).add(" (").add(n.astype(str)).add(")").to_string())
i = med.sag.loc[["MC4R", "L1 VIP"]].min() > med.sag.loc[["LAMP5", "PAX6"]].max()
ii = med.sag.idxmin() == "PAX6"
# (iii) as amended by the reader's answer (questions/000630.yaml): PAX6 above LAMP5 and MC4R
iii = med.first_isi_inv_hero.loc["PAX6"] > med.first_isi_inv_hero.loc[["LAMP5", "MC4R"]].max()
iii_before = med.first_isi_inv_hero.idxmax() == "PAX6"
print(f" (i) MC4R and L1 VIP sag above LAMP5 and PAX6: {i}; (ii) PAX6 lowest sag: {ii}; "
f"(iii) PAX6 initial rate above LAMP5 and MC4R: {iii}")
print(f" criterion before the answer, (iii) PAX6 highest of the four: {iii_before}")
for f in ["sag", "first_isi_inv_hero"]:
print(f" Dunn FDR p, {f}:", {k: f"{v:.2g}" for k, v in pairs[f].items()})
print(f" range of sag {sub.sag.min():.3f} to {sub.sag.max():.3f}; first_isi_inv_hero "
f"{sub.first_isi_inv_hero.min():.1f} to {sub.first_isi_inv_hero.max():.1f} Hz")
Reading (a): medians (n)
sag first_isi_inv_hero
homology_type
LAMP5 0.116 (96) 32.073 (96)
MC4R 0.336 (35) 30.508 (34)
PAX6 0.051 (33) 49.261 (33)
L1 VIP 0.277 (15) 54.645 (15)
(i) MC4R and L1 VIP sag above LAMP5 and PAX6: True; (ii) PAX6 lowest sag: True; (iii) PAX6 initial rate above LAMP5 and MC4R: True
criterion before the answer, (iii) PAX6 highest of the four: False
Dunn FDR p, sag: {'LAMP5 - MC4R': '2.1e-13', 'LAMP5 - PAX6': '0.00063', 'LAMP5 - L1 VIP': '3.1e-05', 'MC4R - PAX6': '1.7e-18', 'MC4R - L1 VIP': '0.34', 'PAX6 - L1 VIP': '3e-09'}
Dunn FDR p, first_isi_inv_hero: {'LAMP5 - MC4R': '0.78', 'LAMP5 - PAX6': '6.8e-06', 'LAMP5 - L1 VIP': '0.0008', 'MC4R - PAX6': '3.4e-05', 'MC4R - L1 VIP': '0.0008', 'PAX6 - L1 VIP': '1'}
range of sag 0.003 to 0.546; first_isi_inv_hero 6.9 to 176.1 Hz
Reading (a'): medians (n)
sag first_isi_inv_hero
homology_type
LAMP5 0.117 (96) 32.073 (96)
MC4R 0.336 (35) 29.922 (35)
PAX6 0.051 (33) 49.261 (33)
L1 VIP 0.277 (15) 54.645 (15)
(i) MC4R and L1 VIP sag above LAMP5 and PAX6: True; (ii) PAX6 lowest sag: True; (iii) PAX6 initial rate above LAMP5 and MC4R: True
criterion before the answer, (iii) PAX6 highest of the four: False
Dunn FDR p, sag: {'LAMP5 - MC4R': '2.8e-13', 'LAMP5 - PAX6': '0.00032', 'LAMP5 - L1 VIP': '3.7e-05', 'MC4R - PAX6': '5.7e-19', 'MC4R - L1 VIP': '0.33', 'PAX6 - L1 VIP': '1.8e-09'}
Dunn FDR p, first_isi_inv_hero: {'LAMP5 - MC4R': '0.68', 'LAMP5 - PAX6': '6.3e-08', 'LAMP5 - L1 VIP': '0.00029', 'MC4R - PAX6': '5.2e-07', 'MC4R - L1 VIP': '0.00029', 'PAX6 - L1 VIP': '0.7'}
range of sag 0.003 to 0.546; first_isi_inv_hero 4.1 to 176.1 Hz
Reading (b): medians (n)
sag first_isi_inv_hero
homology_type
LAMP5 0.116 (96) 31.817 (96)
MC4R 0.336 (35) 28.758 (34)
PAX6 0.051 (33) 49.261 (33)
L1 VIP 0.277 (15) 54.645 (15)
(i) MC4R and L1 VIP sag above LAMP5 and PAX6: True; (ii) PAX6 lowest sag: True; (iii) PAX6 initial rate above LAMP5 and MC4R: True
criterion before the answer, (iii) PAX6 highest of the four: False
Dunn FDR p, sag: {'LAMP5 - MC4R': '9.7e-14', 'LAMP5 - PAX6': '0.00083', 'LAMP5 - L1 VIP': '2.9e-05', 'MC4R - PAX6': '1.4e-18', 'MC4R - L1 VIP': '0.31', 'PAX6 - L1 VIP': '3.8e-09'}
Dunn FDR p, first_isi_inv_hero: {'LAMP5 - MC4R': '0.39', 'LAMP5 - PAX6': '1.2e-07', 'LAMP5 - L1 VIP': '0.0028', 'MC4R - PAX6': '2.8e-07', 'MC4R - L1 VIP': '0.0017', 'PAX6 - L1 VIP': '0.39'}
range of sag 0.003 to 0.546; first_isi_inv_hero 1.0 to 176.1 Hz
res, pairs, pairs_mw, feats40, sub = fig2c["a"]
fig, axs = plt.subplots(1, 2, figsize=(9, 3.5))
for ax, f, lab in zip(axs, ["sag", "first_isi_inv_hero"], ["sag ratio", "initial firing rate, rheobase + 40 pA (Hz)"]):
groups = [g[f].dropna().values for _, g in sub.groupby("homology_type", observed=True)]
ax.boxplot(groups, tick_labels=list(sub.homology_type.cat.categories), showfliers=True)
ax.set_ylabel(lab)
plt.tight_layout(); plt.show()
print("All points are within the axes (matplotlib autoscale).")
All points are within the axes (matplotlib autoscale).
Result: reproduced. Part (iii) of the criterion was amended by a reader (questions/000630.yaml) from "PAX6 has the highest median initial firing rate" to "PAX6 above LAMP5 and MC4R". Parts (i) and (ii) hold under both readings: the MC4R (0.336) and L1 VIP (0.277) medians of sag exceed those of LAMP5 (0.116) and PAX6 (0.051), and PAX6 has the lowest. Part (iii) as amended holds under both readings and in variant (a'): PAX6 has a median initial firing rate of 49.3 Hz, above LAMP5 (31.8 to 32.1) and MC4R (28.8 to 30.5), with Dunn FDR p of 3.4e-5 or less. L1 VIP is higher than PAX6 (54.6 Hz under both readings; Dunn FDR p = 0.39 to 1), so the criterion before the answer, PAX6 highest of the four, fails.
figs3a-sag-adaptation¶
Spike adaptation properties are strongly and inversely related to sag across human L1 cells (R2 = 0.74).
As cells 8, 16 and 17: cells of P with l1_type; OLS of upstroke_adapt_ratio on sag with HC3 errors.
rows = []
# Reading (a) uses the paper's adaptation ratio (it is not among the features computed and used),
# with sag computed here; the ratio computed here is shown in both orders.
for kind, T, col in [("a: sag here, paper's ratio (third over first)", TA, "upstroke_adapt_ratio"),
("a': sag here without the 500 ms noise QC, paper's ratio (third over first)", TA2, "upstroke_adapt_ratio"),
("sag and ratio here, third over first", TAR, "upstroke_adapt_ratio"),
("sag and ratio here, first over third", TA.assign(upstroke_adapt_ratio=computed["upstroke_adapt_ratio_1over3"].reindex(TA.index)), "upstroke_adapt_ratio"),
("b: paper's values (third over first)", TB, "upstroke_adapt_ratio")]:
data = T[T.P & T.l1_type].dropna(subset=["sag", col])
fit = smf.ols(f"{col} ~ sag", data=data).fit(cov_type="HC3")
rows.append({"reading": kind, "n": int(fit.nobs), "R2": round(fit.rsquared, 3), "slope": round(fit.params["sag"], 3),
"p (HC3)": f"{fit.f_pvalue:.2g}", "ratio range": f"{data[col].min():.2f} to {data[col].max():.2f}"})
print(pd.DataFrame(rows).to_string(index=False))
# Direction (criterion: cells with more sag show less change of upstroke from the first to the
# third spike). With the ratio third over first, a ratio of 1 is no change and a ratio below 1 a
# smaller upstroke on the third spike; the ratios here are all below 1, so a lower ratio is more change.
for r in rows:
more_change = (r["slope"] > 0) if "first over third" in r["reading"] else (r["slope"] < 0)
print(f" {r['reading']}: slope {r['slope']}, so cells with more sag show "
f"{'more' if more_change else 'less'} change of upstroke from the first to the third spike")
data = TA[TA.P & TA.l1_type].dropna(subset=["sag", "upstroke_adapt_ratio"])
print("Per subclass (reading a), R2 within each:")
for k, g in data.groupby("homology_type", observed=True):
if len(g) > 3:
print(f" {k}: n={len(g)}, R2={smf.ols('upstroke_adapt_ratio ~ sag', data=g).fit().rsquared:.2f}")
for kind, T in [("a", TA), ("b", TB)]:
d = T[T.P & T.l1_type].groupby("t-type")[["sag", "upstroke_adapt_ratio"]].median()
fit = smf.ols("upstroke_adapt_ratio ~ sag", data=d).fit()
print(f"Variant, t-type medians ({kind}): {len(d)} t-types, R2 {fit.rsquared:.2f}, slope {fit.params['sag']:.2f}, p {fit.f_pvalue:.2g}")
plt.figure(figsize=(4, 3.5)); plt.scatter(data.sag, data.upstroke_adapt_ratio, s=6)
plt.xlabel("sag ratio"); plt.ylabel("upstroke adaptation ratio (3rd/1st)"); plt.show()
reading n R2 slope p (HC3) ratio range
a: sag here, paper's ratio (third over first) 192 0.334 -0.608 1.3e-12 0.38 to 0.97
a': sag here without the 500 ms noise QC, paper's ratio (third over first) 192 0.346 -0.622 1.4e-13 0.38 to 0.97
sag and ratio here, third over first 192 0.204 -0.440 1.6e-07 0.36 to 0.98
sag and ratio here, first over third 192 0.149 0.771 4.4e-05 1.02 to 2.76
b: paper's values (third over first) 192 0.359 -0.637 3.3e-14 0.38 to 0.97
a: sag here, paper's ratio (third over first): slope -0.608, so cells with more sag show more change of upstroke from the first to the third spike
a': sag here without the 500 ms noise QC, paper's ratio (third over first): slope -0.622, so cells with more sag show more change of upstroke from the first to the third spike
sag and ratio here, third over first: slope -0.44, so cells with more sag show more change of upstroke from the first to the third spike
sag and ratio here, first over third: slope 0.771, so cells with more sag show more change of upstroke from the first to the third spike
b: paper's values (third over first): slope -0.637, so cells with more sag show more change of upstroke from the first to the third spike
Per subclass (reading a), R2 within each:
LAMP5: n=96, R2=0.22
MC4R: n=33, R2=0.20 PAX6: n=33, R2=0.11 L1 VIP: n=15, R2=0.00 other: n=15, R2=0.00 Variant, t-type medians (a): 10 t-types, R2 0.28, slope -0.50, p 0.11 Variant, t-type medians (b): 10 t-types, R2 0.31, slope -0.52, p 0.097
Result: not reproduced. The regression of the upstroke adaptation ratio (third over first spike, the order that matches the paper's table) on sag gives R2 = 0.33 (p = 1.3e-12) with sag computed here and the paper's ratio (reading (a), as the ratio computed here is not used there), and R2 = 0.36 (p = 3.3e-14) on the paper's own per-cell values, against 0.74 +/- 0.074. With both sag and the ratio computed here R2 is 0.20 (third over first) or 0.15 (first over third); variant (a') gives 0.35. All the paper's ratios are below 1 (0.38 to 0.97), so the upstroke falls from the first to the third spike, and the slope is negative (-0.61 in reading (a), -0.64 on the paper's values): cells with more sag have a lower ratio, a larger fall, so their upstroke changes more from the first to the third spike. The criterion asks for less change, so its direction also fails. The paper says only "inverse relationship" (Fig. S3A), which the negative slope of the third-over-first ratio matches; whether the criterion's direction misreads the paper is asked of a reader (questions/000630.yaml). It cannot change the outcome, which fails on R2. On the medians of the ten t-types R2 is 0.28 to 0.31.
fig2c-lamp5-gradient¶
Within the LAMP5 subclass, sag and initial firing rate decrease from the NMBR to the DBP to the LCP2 t-type.
As cell 59 (which calls run_correlations with filter_dropouts=False, as here; the
function's default is shown as a variant): LAMP5 cells of the strict L1 set, t-type coded LCP2 = 0, DBP = 1, NMBR = 2; Spearman
correlation of every feature whose name has no "3"; Benjamini-Hochberg. The code also includes
the Amsterdam and Szeged LAMP5 cells, which the dandiset does not hold.
order = {"LAMP5 LCP2": 0, "LAMP5 DBP": 1, "LAMP5 NMBR": 2}
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
data = T[T.in_strict_set & (T.homology_type == "LAMP5")].assign(cluster_i=lambda d: d["t-type"].map(order))
feats_no3 = [c for c in fcols if "3" not in c]
res = L.run_correlations(data, feats_no3, "cluster_i", filter_dropouts=False)
print(f"Reading ({kind}): {len(data)} LAMP5 cells", data["t-type"].value_counts().to_dict(),
f"| {len(res)} features, {int((res.pval_fdr < 0.001).sum())} with FDR p < 0.001")
print(res.loc[["sag", "first_isi_inv_hero"]].to_string())
print(" medians by t-type:", data.groupby("t-type")[["sag", "first_isi_inv_hero"]].median().round(3).to_dict())
print(" donors by t-type:", data.groupby("t-type").donor_name.nunique().to_dict())
d3 = sites[sites.in_strict_set & (sites.homology_type == "LAMP5")].assign(cluster_i=lambda d: d["t-type"].map(order))
res3 = L.run_correlations(d3, [c for c in fcols if "3" not in c and c in d3], "cluster_i", filter_dropouts=False)
print(f"Variant with the Amsterdam and Szeged LAMP5 cells as the code has them (their features from the paper's tables): "
f"{len(d3)} cells", d3.collaborator.value_counts().to_dict())
print(res3.loc[["sag", "first_isi_inv_hero"]].to_string())
# The code's cell 59 calls run_correlations with filter_dropouts=False, as here. With the function's
# default (True) a feature with fewer than 5 values above 1 is left out first; sag is a ratio below 1.
data = TA[TA.in_strict_set & (TA.homology_type == "LAMP5")].assign(cluster_i=lambda d: d["t-type"].map(order))
resT = L.run_correlations(data, [c for c in fcols if "3" not in c], "cluster_i", filter_dropouts=True)
print(f"Variant with filter_dropouts=True (reading a): {len(resT)} features kept; sag kept: {'sag' in resT.index}; "
f"first_isi_inv_hero FDR p {resT.loc['first_isi_inv_hero', 'pval_fdr']:.2g}; sag values above 1: {int((data.sag > 1).sum())}")
Reading (a): 96 LAMP5 cells {'LAMP5 NMBR': 71, 'LAMP5 LCP2': 13, 'LAMP5 DBP': 12} | 93 features, 14 with FDR p < 0.001
r pval n pval_fdr
feature
sag 0.276556 0.006380 96 0.021975
first_isi_inv_hero 0.456559 0.000003 96 0.000026
medians by t-type: {'sag': {'LAMP5 DBP': 0.114, 'LAMP5 LCP2': 0.077, 'LAMP5 NMBR': 0.125}, 'first_isi_inv_hero': {'LAMP5 DBP': 18.571, 'LAMP5 LCP2': 22.452, 'LAMP5 NMBR': 33.996}}
donors by t-type: {'LAMP5 DBP': 10, 'LAMP5 LCP2': 12, 'LAMP5 NMBR': 27}
Reading (a'): 96 LAMP5 cells {'LAMP5 NMBR': 71, 'LAMP5 LCP2': 13, 'LAMP5 DBP': 12} | 93 features, 14 with FDR p < 0.001
r pval n pval_fdr
feature
sag 0.291429 0.003966 96 0.014755
first_isi_inv_hero 0.453410 0.000003 96 0.000030
medians by t-type: {'sag': {'LAMP5 DBP': 0.114, 'LAMP5 LCP2': 0.077, 'LAMP5 NMBR': 0.125}, 'first_isi_inv_hero': {'LAMP5 DBP': 18.571, 'LAMP5 LCP2': 22.452, 'LAMP5 NMBR': 33.996}}
donors by t-type: {'LAMP5 DBP': 10, 'LAMP5 LCP2': 12, 'LAMP5 NMBR': 27}
Reading (b): 96 LAMP5 cells {'LAMP5 NMBR': 71, 'LAMP5 LCP2': 13, 'LAMP5 DBP': 12} | 93 features, 14 with FDR p < 0.001
r pval n pval_fdr
feature
sag 0.359737 0.000318 96 0.001974
first_isi_inv_hero 0.458262 0.000003 96 0.000025
medians by t-type: {'sag': {'LAMP5 DBP': 0.092, 'LAMP5 LCP2': 0.077, 'LAMP5 NMBR': 0.125}, 'first_isi_inv_hero': {'LAMP5 DBP': 19.554, 'LAMP5 LCP2': 21.626, 'LAMP5 NMBR': 33.996}}
donors by t-type: {'LAMP5 DBP': 10, 'LAMP5 LCP2': 12, 'LAMP5 NMBR': 27}
Variant with the Amsterdam and Szeged LAMP5 cells as the code has them (their features from the paper's tables): 113 cells {'AIBS': 99, 'Gabor': 10, 'Huib': 4}
r pval n pval_fdr
feature
sag 0.314007 0.000935 108 0.003906
first_isi_inv_hero 0.385549 0.000038 108 0.000351
Variant with filter_dropouts=True (reading a): 31 features kept; sag kept: False; first_isi_inv_hero FDR p 3.2e-05; sag values above 1: 0
Result: not reproduced. Both features have a positive Spearman rho with the code LCP2 = 0, DBP = 1, NMBR = 2 (sag 0.28 here, 0.29 in variant (a') and 0.36 on the paper's values, initial firing rate 0.45 to 0.46), the direction the criterion asks for; the medians are highest in NMBR, but the initial firing rate is not monotone (DBP 18.6 Hz, LCP2 22.5 Hz). The initial firing rate meets the claim's threshold (FDR p = 2.6e-5 here, 3.0e-5 in (a'), 2.5e-5 on the paper's values); sag does not (FDR p = 0.022 here, 0.015 in (a'), 0.0020 on the paper's values, against 0.001). With the Amsterdam and Szeged LAMP5 cells added from the paper's tables, as the code includes them, sag has FDR p = 0.0039, still above 0.001. With the function's default dropout filter, sag (all values below 1) is left out of the test. The criterion names no partial level, and sag fails under every reading, so the claim is not reproduced, by the rule applied to fig4b-adarb2-irregular-firing (changed from partially reproduced after the third round of verification). Even the paper's own values give sag FDR p above 0.001 under this test, which is asked of a reader.
figs6-donor-condition-ephys¶
No electrophysiological feature of human L1 cells differs with donor medical condition, sex or age after FDR correction.
A claim of absence, judged on the paper's tests (l1_metadata_effects.ipynb cells 50, 51, 99, 100,
108 to 110): cells of P with l1_type; features more than 40% non-missing; epileptic or not
(code's epileptic()) and sex by two-sided exact Mann-Whitney; age by Spearman with the dropout
filter; Benjamini-Hochberg within each. The table's medical_conditions, donor_sex and
donor_age stand in for the LIMS records the code reads.
def epileptic(x):
if pd.isna(x):
return x
return "epileptic" if "epilep" in x.lower() else "non-epileptic"
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
data = T[T.P & T.l1_type].copy()
data["epileptic"] = data.medical_conditions.map(epileptic)
data["sex_t"] = data.donor_sex
data["age_years"] = data.donor_age.str.extract(r"(\d+)")[0].astype(float)
feats40 = [c for c in fcols if data[c].notna().mean() > 0.4]
print(f"Reading ({kind}): {len(data)} cells, {len(feats40)} features")
print(" epileptic:", data.epileptic.value_counts(dropna=False).to_dict(), "| donors with a condition:",
data.dropna(subset=["epileptic"]).donor_name.nunique(), "of", data.donor_name.nunique())
print(" sex:", data.sex_t.value_counts(dropna=False).to_dict(), "| age known:", int(data.age_years.notna().sum()))
for name, res in [("medical condition", L.run_twosamp(data, feats40, "epileptic")),
("sex", L.run_twosamp(data, feats40, "sex_t"))]:
print(f" {name}: smallest p {res.pval_mw.min():.3g} ({res.pval_mw.idxmin()}, r={res.loc[res.pval_mw.idxmin(), 'mw_r']:.2f}), "
f"smallest FDR p {res.pval_mw_fdr.min():.3g}, features with FDR p < 0.05: {int((res.pval_mw_fdr < 0.05).sum())}")
res = L.run_correlations(data, feats40, "age_years", filter_dropouts=True)
print(f" age: {len(res)} features after the dropout filter; smallest p {res.pval.min():.3g} ({res.pval.idxmin()}, "
f"rho={res.loc[res.pval.idxmin(), 'r']:.2f}), smallest FDR p {res.pval_fdr.min():.3g}, "
f"features with FDR p < 0.05: {int((res.pval_fdr < 0.05).sum())}")
dropped = [c for c in feats40 if c not in res.index]
print(f" age, features left out by the dropout filter: {len(dropped)}, including sag: {'sag' in dropped}, "
f"adaptation ratios: {sum('adapt_ratio' in c for c in dropped)}")
res_nf = L.run_correlations(data, feats40, "age_years", filter_dropouts=False)
print(f" age, variant without the dropout filter: {len(res_nf)} features; smallest p {res_nf.pval.min():.3g} "
f"({res_nf.pval.idxmin()}), smallest FDR p {res_nf.pval_fdr.min():.3g}, features with FDR p < 0.05: "
f"{int((res_nf.pval_fdr < 0.05).sum())}; sag p {res_nf.loc['sag', 'pval']:.2g}")
Reading (a): 194 cells, 92 features
epileptic: {nan: 108, 'epileptic': 65, 'non-epileptic': 21} | donors with a condition: 36 of 52
sex: {'M': 129, 'F': 65} | age known: 192
medical condition: smallest p 0.00267 (sag, r=0.43), smallest FDR p 0.194, features with FDR p < 0.05: 0 sex: smallest p 0.0113 (rheobase_i, r=0.22), smallest FDR p 0.982, features with FDR p < 0.05: 0 age: 36 features after the dropout filter; smallest p 0.0388 (peak_deltav_hero, rho=-0.15), smallest FDR p 0.426, features with FDR p < 0.05: 0 age, features left out by the dropout filter: 56, including sag: True, adaptation ratios: 3
age, variant without the dropout filter: 92 features; smallest p 0.0388 (peak_deltav_hero), smallest FDR p 0.726, features with FDR p < 0.05: 0; sag p 0.11
Reading (a'): 194 cells, 92 features
epileptic: {nan: 108, 'epileptic': 65, 'non-epileptic': 21} | donors with a condition: 36 of 52
sex: {'M': 129, 'F': 65} | age known: 192
medical condition: smallest p 0.00728 (sag, r=0.39), smallest FDR p 0.194, features with FDR p < 0.05: 0 sex: smallest p 0.0203 (rheobase_i, r=0.20), smallest FDR p 0.99, features with FDR p < 0.05: 0 age: 36 features after the dropout filter; smallest p 0.02 (peak_deltav_hero, rho=-0.17), smallest FDR p 0.426, features with FDR p < 0.05: 0 age, features left out by the dropout filter: 56, including sag: True, adaptation ratios: 3
age, variant without the dropout filter: 92 features; smallest p 0.02 (peak_deltav_hero), smallest FDR p 0.805, features with FDR p < 0.05: 0; sag p 0.18
Reading (b): 194 cells, 92 features
epileptic: {nan: 108, 'epileptic': 65, 'non-epileptic': 21} | donors with a condition: 36 of 52
sex: {'M': 129, 'F': 65} | age known: 192
medical condition: smallest p 0.00582 (sag, r=0.40), smallest FDR p 0.241, features with FDR p < 0.05: 0 sex: smallest p 0.0175 (rheobase_i, r=0.21), smallest FDR p 0.996, features with FDR p < 0.05: 0 age: 36 features after the dropout filter; smallest p 0.0158 (peak_deltav_hero, rho=-0.17), smallest FDR p 0.451, features with FDR p < 0.05: 0 age, features left out by the dropout filter: 56, including sag: True, adaptation ratios: 3
age, variant without the dropout filter: 92 features; smallest p 0.0158 (peak_deltav_hero), smallest FDR p 0.69, features with FDR p < 0.05: 0; sag p 0.1
Result: reproduced. No feature differs with medical condition, sex or age after FDR correction under either reading (smallest FDR p 0.19, 0.98 and 0.43 on the features computed here) or in variant (a'). Two limits weaken this absence of an effect. The medical condition is known for only 86 of the 194 cells in the public table (65 epileptic, 21 not), so the condition test has less power than the paper's, which read the conditions from the Allen Institute's donor database. The dropout filter the criterion takes from the code leaves 36 of 92 features in the age test; sag and three adaptation ratios are among the 56 left out. Without the filter no feature reaches FDR p < 0.05 for age either (smallest FDR p 0.73; sag p = 0.11).
figs6c-region-ephys¶
Some electrophysiological features of human L1 cells, including input resistance, vary with brain region; sparsely sampled regions either match MTG or have larger input resistance.
As l1_metadata_effects.ipynb cells 15, 29 and 30. input_resistance_ss is not computed here, so
it comes from the paper's table under both readings.
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
df = T[T.l1_type & (T.collaborator == "AIBS") & (T.failed_fx_long_squares == False)].copy()
df.loc[df.donor_name == "H18.03.003", "structure"] = "CeG"
df.loc[df.donor_name == "H22.26.404", "structure"] = "SFG"
vc = df.structure.value_counts()
data = df[df.structure.isin(vc[vc >= 5].index) & (df.structure != "TemL")].copy()
data["structure"] = pd.Categorical(data.structure, categories=["ITG", "MTG", "MFG", "SFG"], ordered=True)
feats50 = [c for c in fcols if data[c].notna().mean() > 0.5]
res, pairs = L.run_kw_dunn(data, feats50, "structure")
print(f"Reading ({kind}): regions {data.structure.value_counts().to_dict()} (all regions {vc.to_dict()})")
print(f" {len(feats50)} features, {int((res.pval_fdr < 0.05).sum())} with FDR p < 0.05:", list(res.index[res.pval_fdr < 0.05]))
for f in ["input_resistance_ss", "sag_peak_t", "downstroke_rheo", "avg_rate_max"]:
if f in res.index:
print(f" {f}: FDR p {res.loc[f, 'pval_fdr']:.3g}; medians", data.groupby("structure", observed=True)[f].median().round(3).to_dict(),
"; Dunn FDR p", {k: f"{v:.2g}" for k, v in pairs.get(f, {}).items()})
print(" donors per region:", data.groupby("structure", observed=True).donor_name.nunique().to_dict())
Reading (a): regions {'MTG': 115, 'ITG': 23, 'SFG': 14, 'MFG': 5} (all regions {'MTG': 115, 'TemL': 26, 'ITG': 23, 'SFG': 14, 'MFG': 5, 'STG': 4, 'CeG': 2, 'PRG': 2, 'IFGtr': 1, 'LOrG': 1, 'OccL': 1})
83 features, 6 with FDR p < 0.05: ['sag_peak_t', 'downstroke_ramp', 'input_resistance_ss', 'downstroke_rheo', 'avg_rate_max', 'threshold_i_ramp']
input_resistance_ss: FDR p 0.0293; medians {'ITG': 291.047, 'MTG': 234.133, 'MFG': 367.196, 'SFG': 221.768} ; Dunn FDR p {'ITG - MTG': '0.075', 'ITG - MFG': '0.075', 'ITG - SFG': '0.024', 'MTG - MFG': '0.0073', 'MTG - SFG': '0.13', 'MFG - SFG': '0.003'}
sag_peak_t: FDR p 0.0156; medians {'ITG': 0.059, 'MTG': 0.05, 'MFG': 0.169, 'SFG': 0.044} ; Dunn FDR p {'ITG - MTG': '0.14', 'ITG - MFG': '0.075', 'ITG - SFG': '0.002', 'MTG - MFG': '0.0089', 'MTG - SFG': '0.0081', 'MFG - SFG': '0.00042'}
downstroke_rheo: FDR p 0.0293; medians {'ITG': -106.782, 'MTG': -93.36, 'MFG': -81.271, 'SFG': -113.841} ; Dunn FDR p {'ITG - MTG': '0.021', 'ITG - MFG': '0.014', 'ITG - SFG': '0.67', 'MTG - MFG': '0.081', 'MTG - SFG': '0.021', 'MFG - SFG': '0.014'}
avg_rate_max: FDR p 0.0293; medians {'ITG': 34.0, 'MTG': 27.0, 'MFG': 41.0, 'SFG': 27.0} ; Dunn FDR p {'ITG - MTG': '0.022', 'ITG - MFG': '0.31', 'ITG - SFG': '0.043', 'MTG - MFG': '0.022', 'MTG - SFG': '0.78', 'MFG - SFG': '0.025'}
donors per region: {'ITG': 4, 'MTG': 27, 'MFG': 3, 'SFG': 5}
Reading (a'): regions {'MTG': 115, 'ITG': 23, 'SFG': 14, 'MFG': 5} (all regions {'MTG': 115, 'TemL': 26, 'ITG': 23, 'SFG': 14, 'MFG': 5, 'STG': 4, 'CeG': 2, 'PRG': 2, 'IFGtr': 1, 'LOrG': 1, 'OccL': 1})
83 features, 13 with FDR p < 0.05: ['sag_peak_t', 'downstroke_ramp', 'input_resistance_ss', 'downstroke_rheo', 'avg_rate_max', 'threshold_i_ramp', 'downstroke_hero', 'width_suprathresh_ramp', 'sag', 'upstroke_short_square', 'downstroke_short_square', 'width_suprathresh_rheo', 'fast_trough_deltav_ramp']
input_resistance_ss: FDR p 0.0293; medians {'ITG': 291.047, 'MTG': 234.133, 'MFG': 367.196, 'SFG': 221.768} ; Dunn FDR p {'ITG - MTG': '0.075', 'ITG - MFG': '0.075', 'ITG - SFG': '0.024', 'MTG - MFG': '0.0073', 'MTG - SFG': '0.13', 'MFG - SFG': '0.003'}
sag_peak_t: FDR p 0.0156; medians {'ITG': 0.059, 'MTG': 0.05, 'MFG': 0.169, 'SFG': 0.044} ; Dunn FDR p {'ITG - MTG': '0.14', 'ITG - MFG': '0.075', 'ITG - SFG': '0.002', 'MTG - MFG': '0.0089', 'MTG - SFG': '0.0081', 'MFG - SFG': '0.00042'}
downstroke_rheo: FDR p 0.0293; medians {'ITG': -106.782, 'MTG': -93.36, 'MFG': -81.271, 'SFG': -113.841} ; Dunn FDR p {'ITG - MTG': '0.021', 'ITG - MFG': '0.014', 'ITG - SFG': '0.68', 'MTG - MFG': '0.079', 'MTG - SFG': '0.021', 'MFG - SFG': '0.014'}
avg_rate_max: FDR p 0.0293; medians {'ITG': 34.0, 'MTG': 27.0, 'MFG': 41.0, 'SFG': 27.0} ; Dunn FDR p {'ITG - MTG': '0.022', 'ITG - MFG': '0.31', 'ITG - SFG': '0.043', 'MTG - MFG': '0.022', 'MTG - SFG': '0.78', 'MFG - SFG': '0.025'}
donors per region: {'ITG': 4, 'MTG': 27, 'MFG': 3, 'SFG': 5}
Reading (b): regions {'MTG': 115, 'ITG': 23, 'SFG': 14, 'MFG': 5} (all regions {'MTG': 115, 'TemL': 26, 'ITG': 23, 'SFG': 14, 'MFG': 5, 'STG': 4, 'CeG': 2, 'PRG': 2, 'IFGtr': 1, 'LOrG': 1, 'OccL': 1})
83 features, 6 with FDR p < 0.05: ['sag_peak_t', 'downstroke_rheo', 'downstroke_ramp', 'input_resistance_ss', 'avg_rate_max', 'threshold_i_ramp']
input_resistance_ss: FDR p 0.0289; medians {'ITG': 291.047, 'MTG': 234.133, 'MFG': 367.196, 'SFG': 221.768} ; Dunn FDR p {'ITG - MTG': '0.075', 'ITG - MFG': '0.075', 'ITG - SFG': '0.024', 'MTG - MFG': '0.0073', 'MTG - SFG': '0.13', 'MFG - SFG': '0.003'}
sag_peak_t: FDR p 0.0156; medians {'ITG': 0.059, 'MTG': 0.05, 'MFG': 0.169, 'SFG': 0.044} ; Dunn FDR p {'ITG - MTG': '0.14', 'ITG - MFG': '0.075', 'ITG - SFG': '0.002', 'MTG - MFG': '0.0089', 'MTG - SFG': '0.0081', 'MFG - SFG': '0.00042'}
downstroke_rheo: FDR p 0.0289; medians {'ITG': -106.782, 'MTG': -93.418, 'MFG': -81.271, 'SFG': -113.841} ; Dunn FDR p {'ITG - MTG': '0.017', 'ITG - MFG': '0.013', 'ITG - SFG': '0.69', 'MTG - MFG': '0.084', 'MTG - SFG': '0.017', 'MFG - SFG': '0.013'}
avg_rate_max: FDR p 0.0293; medians {'ITG': 34.0, 'MTG': 27.0, 'MFG': 41.0, 'SFG': 27.0} ; Dunn FDR p {'ITG - MTG': '0.022', 'ITG - MFG': '0.31', 'ITG - SFG': '0.043', 'MTG - MFG': '0.022', 'MTG - SFG': '0.78', 'MFG - SFG': '0.025'}
donors per region: {'ITG': 4, 'MTG': 27, 'MFG': 3, 'SFG': 5}
Result: reproduced. Over ITG, MTG, MFG and SFG, input resistance (input_resistance_ss,
the paper's values under both readings) varies with region (FDR p = 0.029), as do sag_peak_t,
downstroke_rheo and avg_rate_max (FDR p 0.016 to 0.029). No region has a median input resistance
below MTG's at Dunn FDR p < 0.05: SFG is lower (222 against 234 MOhm) with p = 0.13, and MFG is
higher (367 MOhm, p = 0.0073). MFG holds 5 cells from 3 donors. Variant (a') gives the
same input resistance result; 13 features have FDR p < 0.05 there.
fig4b-mc4r-sag-irregularity¶
MC4R (rosehip) cells are distinguished from the rest of human L1 by strong sag and by irregular firing (Mann-Whitney, FDR p < 1e-4).
As l1_rosehip.ipynb cell 6: AIBS cells of the strict L1 set, group against all others; exact
Mann-Whitney per feature over all features; Benjamini-Hochberg. The irregularity ratio
(ir_late_mean) is custom code that is not public, so it comes from the paper's table (O1).
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
base = T[T.in_strict_set & (T.collaborator == "AIBS")]
for grp_name, mask in [("MC4R subclass", base.homology_type == "MC4R"),
("SST CHRNA4", base["t-type"] == "SST CHRNA4"),
("ADARB2 MC4R", base["t-type"] == "ADARB2 MC4R")]:
df = base.assign(in_type=mask.map({True: "in", False: "rest"}))
res = L.run_twosamp(df, fcols, "in_type")
med = df.groupby("in_type")[["sag", "ir_late_mean"]].median()
out = {f: f"FDR p {res.loc[f, 'pval_mw_fdr']:.2g}, median in {med.loc['in', f]:.3f} vs rest {med.loc['rest', f]:.3f}"
for f in ["sag", "ir_late_mean"]}
print(f"({kind}) {grp_name} (n={int(mask.sum())} of {len(df)}):", out)
(a) MC4R subclass (n=35 of 210): {'sag': 'FDR p 1.6e-18, median in 0.336 vs rest 0.118', 'ir_late_mean': 'FDR p 1.2e-05, median in 0.232 vs rest 0.135'}
(a) SST CHRNA4 (n=9 of 210): {'sag': 'FDR p 0.084, median in 0.234 vs rest 0.131', 'ir_late_mean': 'FDR p 0.6, median in 0.153 vs rest 0.139'}
(a) ADARB2 MC4R (n=26 of 210): {'sag': 'FDR p 3.1e-17, median in 0.374 vs rest 0.121', 'ir_late_mean': 'FDR p 4.2e-06, median in 0.421 vs rest 0.136'}
(a') MC4R subclass (n=35 of 210): {'sag': 'FDR p 8.6e-19, median in 0.336 vs rest 0.118', 'ir_late_mean': 'FDR p 1.2e-05, median in 0.232 vs rest 0.135'}
(a') SST CHRNA4 (n=9 of 210): {'sag': 'FDR p 0.081, median in 0.234 vs rest 0.131', 'ir_late_mean': 'FDR p 0.6, median in 0.153 vs rest 0.139'}
(a') ADARB2 MC4R (n=26 of 210): {'sag': 'FDR p 2.5e-17, median in 0.374 vs rest 0.123', 'ir_late_mean': 'FDR p 4.2e-06, median in 0.421 vs rest 0.136'}
(b) MC4R subclass (n=35 of 210): {'sag': 'FDR p 4.1e-19, median in 0.336 vs rest 0.118', 'ir_late_mean': 'FDR p 1.2e-05, median in 0.232 vs rest 0.135'}
(b) SST CHRNA4 (n=9 of 210): {'sag': 'FDR p 0.084, median in 0.234 vs rest 0.129', 'ir_late_mean': 'FDR p 0.65, median in 0.153 vs rest 0.139'}
(b) ADARB2 MC4R (n=26 of 210): {'sag': 'FDR p 1.6e-17, median in 0.374 vs rest 0.123', 'ir_late_mean': 'FDR p 4.2e-06, median in 0.421 vs rest 0.136'}
Result: reproduced. Under reading (a) of the group, the MC4R subclass against the rest of the Seattle strict L1 cells, sag (median 0.336 against 0.118, FDR p = 1.6e-18) and the irregularity ratio (0.232 against 0.135, FDR p = 1.2e-5) are higher, both below 1e-4. Under reading (b) the ADARB2 MC4R t-type alone meets it (FDR p 3.1e-17 and 4.2e-6); SST CHRNA4 alone does not. Variant (a') gives the same (FDR p 8.6e-19 and 1.2e-5). The irregularity ratio is the paper's value under both feature readings.
fig4b-adarb2-irregular-firing¶
Both rosehip t-types have strong sag, but only ADARB2 MC4R fires irregularly, with stronger and faster sag than SST CHRNA4.
As l1_rosehip.ipynb cell 9: cells of P grouped as LAMP5, SST CHRNA4, ADARB2 MC4R, PAX6, L1 VIP;
Dunn's tests on the seven pairs, Benjamini-Hochberg across them. ir_late_mean, sag_peak_t and
sag_tau come from the paper's table under both readings (not computed here).
pairs7 = [(x, y) for x in ["LAMP5", "PAX6", "L1 VIP"] for y in ["SST CHRNA4", "ADARB2 MC4R"]] + [("SST CHRNA4", "ADARB2 MC4R")]
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
df = T[T.P & T.l1_type & (T.homology_type != "other")].copy()
df["grp"] = np.where(df.homology_type == "MC4R", df["t-type"], df.homology_type.astype(str))
df["grp"] = pd.Categorical(df.grp, categories=["LAMP5", "SST CHRNA4", "ADARB2 MC4R", "PAX6", "L1 VIP"])
med = df.groupby("grp", observed=True)[["ir_late_mean", "sag", "sag_peak_t", "sag_tau"]].median()
print(f"Reading ({kind}): n", df.grp.value_counts(sort=False).to_dict()); print(med.round(4).to_string())
P7 = {f: L.dunn_pairs(df, f, "grp", pairs7) for f in ["ir_late_mean", "sag", "sag_peak_t", "sag_tau"]}
print(pd.DataFrame(P7).map(lambda v: f"{v:.2g}").to_string())
kw = {f: stats.kruskal(*[g[f].dropna() for _, g in df.groupby("grp", observed=True)]).pvalue for f in P7}
print(" KW p:", {k: f"{v:.2g}" for k, v in kw.items()})
i = (med.loc["ADARB2 MC4R", "ir_late_mean"] > med.loc["SST CHRNA4", "ir_late_mean"]) and P7["ir_late_mean"]["SST CHRNA4 - ADARB2 MC4R"] < 0.05
ii = not any((med.loc["SST CHRNA4", "ir_late_mean"] > med.loc[x, "ir_late_mean"]) and P7["ir_late_mean"][f"{x} - SST CHRNA4"] < 0.05
for x in ["LAMP5", "PAX6", "L1 VIP"])
iii = all(med.loc[t, "sag"] > med.loc[["LAMP5", "PAX6"], "sag"].max() for t in ["SST CHRNA4", "ADARB2 MC4R"])
iv_a = med.loc["ADARB2 MC4R", "sag"] > med.loc["SST CHRNA4", "sag"] and med.loc["ADARB2 MC4R", "sag_peak_t"] < med.loc["SST CHRNA4", "sag_peak_t"]
iv_b = med.loc["ADARB2 MC4R", "sag"] > med.loc["SST CHRNA4", "sag"] and med.loc["ADARB2 MC4R", "sag_tau"] < med.loc["SST CHRNA4", "sag_tau"]
print(f" (i) {i} (ii) {ii} (iii) {iii} (iv) with sag_peak_t {iv_a}, with sag_tau {iv_b}")
Reading (a): n {'LAMP5': 96, 'SST CHRNA4': 9, 'ADARB2 MC4R': 26, 'PAX6': 33, 'L1 VIP': 15}
ir_late_mean sag sag_peak_t sag_tau
grp
LAMP5 0.1139 0.1158 0.0484 0.1313
SST CHRNA4 0.1525 0.2338 0.0518 0.0896
ADARB2 MC4R 0.4211 0.3737 0.0447 0.0501
PAX6 0.1467 0.0512 0.1551 0.1102
L1 VIP 0.1842 0.2765 0.0485 0.0918
ir_late_mean sag sag_peak_t sag_tau
LAMP5 - SST CHRNA4 0.046 0.0026 0.64 0.12
LAMP5 - ADARB2 MC4R 1.4e-08 1.3e-12 0.11 9.7e-10
PAX6 - SST CHRNA4 0.74 3.7e-06 0.0023 0.74
PAX6 - ADARB2 MC4R 0.017 9.8e-18 1.6e-12 0.0062
L1 VIP - SST CHRNA4 0.88 0.86 0.64 0.62
L1 VIP - ADARB2 MC4R 0.18 0.23 0.24 0.12
SST CHRNA4 - ADARB2 MC4R 0.18 0.23 0.17 0.063
KW p: {'ir_late_mean': '4.3e-09', 'sag': '2e-21', 'sag_peak_t': '3e-14', 'sag_tau': '5.3e-09'}
(i) False (ii) False (iii) True (iv) with sag_peak_t True, with sag_tau True
Reading (a'): n {'LAMP5': 96, 'SST CHRNA4': 9, 'ADARB2 MC4R': 26, 'PAX6': 33, 'L1 VIP': 15}
ir_late_mean sag sag_peak_t sag_tau
grp
LAMP5 0.1139 0.1169 0.0484 0.1313
SST CHRNA4 0.1525 0.2338 0.0518 0.0896
ADARB2 MC4R 0.4211 0.3737 0.0447 0.0501
PAX6 0.1467 0.0506 0.1551 0.1102
L1 VIP 0.1842 0.2765 0.0485 0.0918
ir_late_mean sag sag_peak_t sag_tau
LAMP5 - SST CHRNA4 0.046 0.0027 0.64 0.12
LAMP5 - ADARB2 MC4R 1.4e-08 1.7e-12 0.11 9.7e-10
PAX6 - SST CHRNA4 0.74 2.3e-06 0.0023 0.74
PAX6 - ADARB2 MC4R 0.017 3.9e-18 1.6e-12 0.0062
L1 VIP - SST CHRNA4 0.88 0.87 0.64 0.62
L1 VIP - ADARB2 MC4R 0.18 0.23 0.24 0.12
SST CHRNA4 - ADARB2 MC4R 0.18 0.23 0.17 0.063
KW p: {'ir_late_mean': '4.3e-09', 'sag': '1e-21', 'sag_peak_t': '3e-14', 'sag_tau': '5.3e-09'}
(i) False (ii) False (iii) True (iv) with sag_peak_t True, with sag_tau True
Reading (b): n {'LAMP5': 96, 'SST CHRNA4': 9, 'ADARB2 MC4R': 26, 'PAX6': 33, 'L1 VIP': 15}
ir_late_mean sag sag_peak_t sag_tau
grp
LAMP5 0.1139 0.1158 0.0484 0.1313
SST CHRNA4 0.1525 0.2338 0.0518 0.0896
ADARB2 MC4R 0.4211 0.3737 0.0447 0.0501
PAX6 0.1467 0.0512 0.1551 0.1102
L1 VIP 0.1842 0.2765 0.0485 0.0918
ir_late_mean sag sag_peak_t sag_tau
LAMP5 - SST CHRNA4 0.046 0.0018 0.64 0.12
LAMP5 - ADARB2 MC4R 1.4e-08 8e-13 0.11 9.7e-10
PAX6 - SST CHRNA4 0.74 2.7e-06 0.0023 0.74
PAX6 - ADARB2 MC4R 0.017 1e-17 1.6e-12 0.0062
L1 VIP - SST CHRNA4 0.88 0.92 0.64 0.62
L1 VIP - ADARB2 MC4R 0.18 0.25 0.24 0.12
SST CHRNA4 - ADARB2 MC4R 0.18 0.26 0.17 0.063
KW p: {'ir_late_mean': '4.3e-09', 'sag': '1.4e-21', 'sag_peak_t': '3e-14', 'sag_tau': '5.3e-09'}
(i) False (ii) False (iii) True (iv) with sag_peak_t True, with sag_tau True
Result: not reproduced. Both MC4R t-types have a median sag above LAMP5 and PAX6 (iii), and ADARB2 MC4R has a stronger and faster sag than SST CHRNA4 under both readings of "faster" (iv). ADARB2 MC4R has a much higher median irregularity ratio than SST CHRNA4 (0.421 against 0.153), but with 9 SST CHRNA4 cells the Dunn FDR p is 0.18, so (i) fails; and SST CHRNA4 is above LAMP5 (0.153 against 0.114) with Dunn FDR p = 0.046, so (ii) fails. The criterion names no partial level, and the two parts that fail are the claim's point that only ADARB2 MC4R fires irregularly. Variant (a') gives the same. The irregularity ratio and the sag time measures are the paper's values under both readings.
fig5b-tnfaip8l3-ephys¶
PAX6 TNFAIP8L3 cells are distinguished from all other human L1 cells by a high initial firing rate at rheobase and an after-depolarization following the last spike.
As l1_tnfaip.ipynb cells 10, 12 and 14: cells of P, PAX6 TNFAIP8L3 against all others; features
more than 50% non-missing among TNFAIP8L3 cells, plus adp_deltav; exact Mann-Whitney;
Benjamini-Hochberg. postap_slope_rheo ("ADP slope") and adp_v_last_rheo come from the paper's
table (not computed here).
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
df = T[T.P].copy()
df["adp_deltav"] = df.adp_v_last_rheo - df.threshold_v_rheo
df["in_type"] = np.where(df["t-type"] == "PAX6 TNFAIP8L3", "in", "rest")
inn = df[df.in_type == "in"]
fl = [c for c in fcols if inn[c].notna().mean() > 0.5] + ["adp_deltav"]
res = L.run_twosamp(df, fl, "in_type")
print(f"Reading ({kind}): {len(inn)} TNFAIP8L3 cells against {int((df.in_type == 'rest').sum())}; {len(fl)} features")
for f in ["first_isi_inv_rheo", "postap_slope_rheo", "adp_deltav"]:
if f in res.index:
print(f" {f}: TNFAIP8L3 values {inn[f].round(2).tolist()}, rest median {df.loc[df.in_type == 'rest', f].median():.2f} "
f"(n={df.loc[df.in_type == 'rest', f].notna().sum()}), p {res.loc[f, 'pval_mw']:.2g}, FDR p {res.loc[f, 'pval_mw_fdr']:.2g}")
else:
print(f" {f}: not tested (fewer than half the TNFAIP8L3 cells have it)")
print(" smallest FDR p over features:", f"{res.pval_mw_fdr.min():.2g}")
Reading (a): 2 TNFAIP8L3 cells against 208; 82 features first_isi_inv_rheo: TNFAIP8L3 values [136.24, 122.85], rest median 6.47 (n=126), p 0.00049, FDR p 0.0081 postap_slope_rheo: TNFAIP8L3 values [2.4, 2.85], rest median -0.43 (n=207), p 9.2e-05, FDR p 0.0051 adp_deltav: TNFAIP8L3 values [-0.56, -0.09], rest median -11.02 (n=50), p 0.014, FDR p 0.07 smallest FDR p over features: 0.0051
Reading (a'): 2 TNFAIP8L3 cells against 208; 82 features first_isi_inv_rheo: TNFAIP8L3 values [136.24, 122.85], rest median 6.37 (n=127), p 0.00048, FDR p 0.0079 postap_slope_rheo: TNFAIP8L3 values [2.4, 2.85], rest median -0.43 (n=207), p 9.2e-05, FDR p 0.0051 adp_deltav: TNFAIP8L3 values [-0.56, -0.09], rest median -11.02 (n=50), p 0.014, FDR p 0.07 smallest FDR p over features: 0.0051
Reading (b): 2 TNFAIP8L3 cells against 208; 82 features first_isi_inv_rheo: TNFAIP8L3 values [136.24, 122.85], rest median 6.37 (n=129), p 0.00047, FDR p 0.0077 postap_slope_rheo: TNFAIP8L3 values [2.4, 2.85], rest median -0.43 (n=207), p 9.2e-05, FDR p 0.0051 adp_deltav: TNFAIP8L3 values [-0.56, -0.09], rest median -11.02 (n=50), p 0.014, FDR p 0.074 smallest FDR p over features: 0.0051
Result: reproduced. The two PAX6 TNFAIP8L3 cells have initial firing rates at rheobase of 136 and 123 Hz (computed here, equal to the paper's values) against a median of 6.5 Hz in the other cells (FDR p = 0.0081), and ADP slopes of 2.4 and 2.85 V/s (the paper's values) against -0.43 (FDR p = 0.0051); both are below 0.01 and higher in TNFAIP8L3. The ADP voltage variant (adp_deltav) has FDR p = 0.07. Variant (a') gives the same (FDR p 0.0079 and 0.0051).
fig6a-ephys-dprime-subclass-structure¶
Pairwise electrophysiological distinctness (d') of human L1 t-types is smaller within subclasses than between them; the only within-subclass pair with high d' is LAMP5 NMBR against LAMP5 LCP2.
As l1_dprime.ipynb cells 9 to 13 and cluster_pairs.py: cells of P; t-types with at least 5 cells (code) or 4 (legend); features more than 80% non-missing; preprocessing and random forest d' as in the helper. Cells of P outside the L1 t-types take part in the preprocessing, as in the code. Seeds 0, 1 and 2 (a standard choice; the code sets none).
within = {"LAMP5": ["LAMP5 NMBR", "LAMP5 DBP", "LAMP5 LCP2"], "MC4R": ["SST CHRNA4", "ADARB2 MC4R"],
"PAX6": ["PAX6 CDH12", "PAX6 TNFAIP8L3"]}
sub_of = {t: s for s, ts in within.items() for t in ts}
dp_results = {}
tt_all = TA[TA.P]["t-type"].where(TA[TA.P]["t-type"].isin(L.HOMOLOGY)).value_counts()
print("L1 t-types in P:", tt_all.to_dict())
print("t-types kept with at least 5 cells:", list(tt_all[tt_all >= 5].index), "| with at least 4:", list(tt_all[tt_all >= 4].index))
same_sets = set(tt_all[tt_all >= 5].index) == set(tt_all[tt_all >= 4].index)
print("The two rules keep the same t-types:", same_sets, "(so the legend's rule is run only when they differ)")
for kind, T in [("a", TA), ("a'", TA2), ("b", TB)]:
for min_n in ([5] if same_sets else [5, 4]):
d = T[T.P].copy()
d["tt"] = d["t-type"].where(d["t-type"].isin(L.HOMOLOGY))
vc = d.tt.value_counts()
d = d[d.tt.isin(vc[vc >= min_n].index) | d.tt.isna()]
if d.tt.isna().sum() < 5:
d = d[d.tt.notna()]
f80 = [c for c in fcols if d[c].notna().mean() > 0.8]
for seed in [0, 1, 2]:
D = L.pairwise_dprime(d, f80, "tt", seed)
rows = []
for i, a in enumerate(D.index):
for b in D.index[:i]:
same = sub_of.get(a) is not None and sub_of.get(a) == sub_of.get(b)
rows.append({"a": a, "b": b, "within": same, "dprime": D.loc[a, b]})
pr = pd.DataFrame(rows)
dp_results[(kind, min_n, seed)] = (D, pr, len(f80), d.tt.value_counts().to_dict())
for key, (D, pr, nf, n) in dp_results.items():
med_between = pr[~pr.within].dprime.median()
hi_a = pr[pr.within & (pr.dprime >= 2.2)]
hi_b = pr[pr.within & (pr.dprime > med_between)]
print(f"reading {key[0]}, min {key[1]} cells, seed {key[2]}: {nf} features; mean within {pr[pr.within].dprime.mean():.2f}, "
f"between {pr[~pr.within].dprime.mean():.2f} (median {med_between:.2f}); high within (>=2.2): "
f"{[a + ' / ' + b for a, b in zip(hi_a.a, hi_a.b)]}; above median between: {[a + ' / ' + b for a, b in zip(hi_b.a, hi_b.b)]}")
print("\nt-types (min 5):", dp_results[("a", 5, 0)][3])
print("\nWithin-subclass d', reading a, min 5, seed 0:")
print(dp_results[("a", 5, 0)][1].query("within").round(2).to_string(index=False))
print("\nd' matrix, reading a, min 5, seed 0:"); print(dp_results[("a", 5, 0)][0].round(2).to_string())
L1 t-types in P: {'LAMP5 NMBR': 71, 'PAX6 CDH12': 31, 'ADARB2 MC4R': 26, 'VIP TSPAN12': 15, 'LAMP5 LCP2': 13, 'LAMP5 DBP': 12, 'SST CHRNA4': 9, 'VIP PCDH20': 9, 'SST BAGE2': 6, 'PAX6 TNFAIP8L3': 2}
t-types kept with at least 5 cells: ['LAMP5 NMBR', 'PAX6 CDH12', 'ADARB2 MC4R', 'VIP TSPAN12', 'LAMP5 LCP2', 'LAMP5 DBP', 'SST CHRNA4', 'VIP PCDH20', 'SST BAGE2'] | with at least 4: ['LAMP5 NMBR', 'PAX6 CDH12', 'ADARB2 MC4R', 'VIP TSPAN12', 'LAMP5 LCP2', 'LAMP5 DBP', 'SST CHRNA4', 'VIP PCDH20', 'SST BAGE2']
The two rules keep the same t-types: True (so the legend's rule is run only when they differ)
reading a, min 5 cells, seed 0: 81 features; mean within 1.27, between 2.54 (median 2.62); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading a, min 5 cells, seed 1: 81 features; mean within 1.31, between 2.53 (median 2.60); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading a, min 5 cells, seed 2: 81 features; mean within 1.24, between 2.48 (median 2.40); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading a', min 5 cells, seed 0: 81 features; mean within 1.28, between 2.53 (median 2.55); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading a', min 5 cells, seed 1: 81 features; mean within 1.28, between 2.56 (median 2.62); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading a', min 5 cells, seed 2: 81 features; mean within 1.30, between 2.52 (median 2.44); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading b, min 5 cells, seed 0: 81 features; mean within 1.21, between 2.49 (median 2.47); high within (>=2.2): ['LAMP5 NMBR / LAMP5 LCP2']; above median between: []
reading b, min 5 cells, seed 1: 81 features; mean within 1.23, between 2.43 (median 2.44); high within (>=2.2): []; above median between: []
reading b, min 5 cells, seed 2: 81 features; mean within 1.12, between 2.39 (median 2.42); high within (>=2.2): []; above median between: []
t-types (min 5): {'LAMP5 NMBR': 71, 'PAX6 CDH12': 31, 'ADARB2 MC4R': 26, 'VIP TSPAN12': 15, 'LAMP5 LCP2': 13, 'LAMP5 DBP': 12, 'SST CHRNA4': 9, 'VIP PCDH20': 9, 'SST BAGE2': 6}
Within-subclass d', reading a, min 5, seed 0:
a b within dprime
LAMP5 LCP2 LAMP5 DBP True 0.32
LAMP5 NMBR LAMP5 DBP True 1.22
LAMP5 NMBR LAMP5 LCP2 True 2.39
SST CHRNA4 ADARB2 MC4R True 1.14
d' matrix, reading a, min 5, seed 0:
ADARB2 MC4R LAMP5 DBP LAMP5 LCP2 LAMP5 NMBR PAX6 CDH12 SST BAGE2 SST CHRNA4 VIP PCDH20 VIP TSPAN12
ADARB2 MC4R NaN 3.31 2.76 4.16 2.98 0.87 1.14 2.69 1.78
LAMP5 DBP 3.31 NaN 0.32 1.22 1.92 2.25 2.95 2.36 2.59
LAMP5 LCP2 2.76 0.32 NaN 2.39 2.83 3.52 2.39 2.39 3.96
LAMP5 NMBR 4.16 1.22 2.39 NaN 2.72 2.65 2.54 3.87 2.73
PAX6 CDH12 2.98 1.92 2.83 2.72 NaN 1.82 3.35 3.09 2.29
SST BAGE2 0.87 2.25 3.52 2.65 1.82 NaN 1.48 1.60 1.31
SST CHRNA4 1.14 2.95 2.39 2.54 3.35 1.48 NaN 2.78 0.95
VIP PCDH20 2.69 2.36 2.39 3.87 3.09 1.60 2.78 NaN 2.45
VIP TSPAN12 1.78 2.59 3.96 2.73 2.29 1.31 0.95 2.45 NaN
Result: not reproduced. Mean d' within subclasses (1.24 to 1.31) is about half that between them (2.48 to 2.54) on the features computed here, for seeds 0 to 2, and on the paper's values and in variant (a'); this part holds in every run. The criterion asks that, under both of its readings of "high", NMBR against LCP2 be the only high within-subclass pair. With "high" as d' of at least 2.2 it is so on the features computed here and in variant (a') for every seed (2.39 for seed 0), and on the paper's values for seed 0 only (no within-subclass pair reaches 2.2 for seeds 1 and 2). With "high" as above the median between-subclass d' (2.40 to 2.62), no within-subclass pair is high, NMBR against LCP2 included. That part fails under the second reading. Both t-type rules (at least 5 or 4 cells) keep the same nine t-types. The criterion names no partial level, so the claim is not reproduced, by the rule applied to fig4b-adarb2-irregular-firing (changed from partially reproduced after the third round of verification).
fig6d-gmm-subclass-f1¶
Two-cluster Gaussian mixtures on the best pair of electrophysiological features recover each subclass with F1 scores of 0.81 (LAMP5), 0.69 (MC4R), 0.89 (PAX6) and 0.5 (all others).
As clustering_classification.ipynb cells 16 to 21 and clustering.py: AIBS cells with l1_type and
soma in L1 (layer == '1'); features more than 60% non-missing; for each group every single
feature and every pair, a 2-component Gaussian mixture (k-means init, 5 inits, seed 0 before each
fit as the code sets), the component holding most of the group's cells taken as the group; the
best F1. Both feature readings of O1, and as variants (a') and reading (a) with the adaptation
ratios computed here ("a, ratios"), with the code's preprocessing. Soma in L1 is the depth table's layer, as
clustering_classification.ipynb cell 16 has it (no fill from the LIMS layer); the filled layer is a
variant. Preprocessing O3: (b) code, Yeo-Johnson and KNN imputation; (a) Methods, cells missing
more than 60% of features left out, then RobustScaler, KNN, LocalOutlierFactor < -20 set missing,
KNN.
from itertools import combinations
from joblib import Parallel, delayed
from sklearn.impute import KNNImputer
from sklearn.metrics import f1_score
from sklearn.mixture import GaussianMixture
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PowerTransformer
def gmm_f1(X, y):
np.random.seed(0)
pred = GaussianMixture(n_components=2, init_params="kmeans", n_init=5).fit_predict(X)
i_true = int(np.mean(pred[y]) > 0.5)
return f1_score(y, pred == i_true)
def best_subspace(Xall, names, y):
sets = [(i,) for i in range(len(names))] + list(combinations(range(len(names)), 2))
scores = Parallel(n_jobs=-1)(delayed(gmm_f1)(Xall[:, list(s)], y) for s in sets)
k = int(np.argmax(scores))
return scores[k], ", ".join(names[i] for i in sets[k]), len(sets)
# Soma in L1 is the depth table's layer == '1', as clustering_classification.ipynb (cell 16) has
# it, without filling missing layers from the LIMS layer (the code fills them only in
# l1_cross_species.ipynb); the filled rule is a variant.
def with_filled_layer(T):
return T.assign(layer=T.layer.fillna(T.layer_lims.astype("string").str.replace(".0", "", regex=False)))
gmm_out = {}
t0 = time.time()
for feat_kind, prep, layer_rule in [("a", "b", "code"), ("a", "a", "code"), ("b", "b", "code"), ("b", "a", "code"),
("a", "b", "filled"), ("a'", "b", "code"), ("a, ratios", "b", "code")]:
T = {"a": TA, "a'": TA2, "a, ratios": TAR, "b": TB}[feat_kind]
if layer_rule == "filled":
T = with_filled_layer(T)
d = T[(T.collaborator == "AIBS") & T.l1_type & (T.failed_fx_long_squares == False) & (T.layer == "1")].copy()
f60 = [c for c in fcols if d[c].notna().mean() > 0.6]
if prep == "a":
d = d[d[f60].notna().mean(axis=1) >= 0.4]
X = L.preproc_pipeline().fit_transform(d[f60].values)
else:
X = Pipeline([("norm", PowerTransformer()), ("impute", KNNImputer())]).fit_transform(d[f60])
grp = d.homology_type.astype(str).replace({"other": "all others", "L1 VIP": "all others"})
for g in ["LAMP5", "MC4R", "PAX6", "all others"]:
y = (grp == g).values
gmm_out[(feat_kind, prep, layer_rule, g)] = best_subspace(X, f60, y) + (int(y.sum()), len(d))
print(f"GMM search took {time.time() - t0:.0f} s")
exp = {"LAMP5": 0.81, "MC4R": 0.69, "PAX6": 0.89, "all others": 0.5}
for (feat_kind, prep, layer_rule, g), (f1, fs, nsets, ng, n) in gmm_out.items():
print(f"features ({feat_kind}), preprocessing ({prep}), layer {layer_rule}, {g}: n {ng} of {n}, best F1 {f1:.3f} (expected {exp[g]} +/- {exp[g] / 10:.3f}) "
f"on [{fs}] over {nsets} subspaces; within tolerance: {abs(f1 - exp[g]) <= exp[g] / 10 + 1e-9}")
GMM search took 1272 s features (a), preprocessing (b), layer code, LAMP5: n 86 of 177, best F1 0.806 (expected 0.81 +/- 0.081) on [upstroke_downstroke_ratio_ramp, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer code, MC4R: n 31 of 177, best F1 0.703 (expected 0.69 +/- 0.069) on [sag_tau, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer code, PAX6: n 32 of 177, best F1 0.892 (expected 0.89 +/- 0.089) on [fast_trough_deltav_ramp, sag_peak_t] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer code, all others: n 28 of 177, best F1 0.526 (expected 0.5 +/- 0.050) on [upstroke_downstroke_ratio_rheo, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (a), layer code, LAMP5: n 86 of 177, best F1 0.785 (expected 0.81 +/- 0.081) on [width_hero, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (a), layer code, MC4R: n 31 of 177, best F1 0.709 (expected 0.69 +/- 0.069) on [sag, downstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (a), layer code, PAX6: n 32 of 177, best F1 0.896 (expected 0.89 +/- 0.089) on [width_short_square, sag_peak_t] over 3486 subspaces; within tolerance: True features (a), preprocessing (a), layer code, all others: n 28 of 177, best F1 0.490 (expected 0.5 +/- 0.050) on [threshold_v_short_square, sag_area] over 3486 subspaces; within tolerance: True features (b), preprocessing (b), layer code, LAMP5: n 86 of 177, best F1 0.806 (expected 0.81 +/- 0.081) on [upstroke_downstroke_ratio_ramp, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (b), preprocessing (b), layer code, MC4R: n 31 of 177, best F1 0.693 (expected 0.69 +/- 0.069) on [sag_tau, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (b), preprocessing (b), layer code, PAX6: n 32 of 177, best F1 0.892 (expected 0.89 +/- 0.089) on [fast_trough_deltav_ramp, sag_peak_t] over 3486 subspaces; within tolerance: True features (b), preprocessing (b), layer code, all others: n 28 of 177, best F1 0.506 (expected 0.5 +/- 0.050) on [upstroke_downstroke_ratio_ramp, fast_trough_v_rheo] over 3486 subspaces; within tolerance: True features (b), preprocessing (a), layer code, LAMP5: n 86 of 177, best F1 0.802 (expected 0.81 +/- 0.081) on [sag, first_isi_inv_hero] over 3486 subspaces; within tolerance: True features (b), preprocessing (a), layer code, MC4R: n 31 of 177, best F1 0.712 (expected 0.69 +/- 0.069) on [sag, width_adapt_ratio] over 3486 subspaces; within tolerance: True features (b), preprocessing (a), layer code, PAX6: n 32 of 177, best F1 0.896 (expected 0.89 +/- 0.089) on [width_short_square, sag_peak_t] over 3486 subspaces; within tolerance: True features (b), preprocessing (a), layer code, all others: n 28 of 177, best F1 0.495 (expected 0.5 +/- 0.050) on [sag_area, upstroke_downstroke_ratio_hero] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer filled, LAMP5: n 95 of 189, best F1 0.840 (expected 0.81 +/- 0.081) on [rheobase_i, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer filled, MC4R: n 34 of 189, best F1 0.707 (expected 0.69 +/- 0.069) on [sag_tau, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer filled, PAX6: n 32 of 189, best F1 0.879 (expected 0.89 +/- 0.089) on [fast_trough_deltav_ramp, sag_peak_t] over 3486 subspaces; within tolerance: True features (a), preprocessing (b), layer filled, all others: n 28 of 189, best F1 0.494 (expected 0.5 +/- 0.050) on [upstroke_downstroke_ratio_hero, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a'), preprocessing (b), layer code, LAMP5: n 86 of 177, best F1 0.806 (expected 0.81 +/- 0.081) on [upstroke_downstroke_ratio_ramp, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a'), preprocessing (b), layer code, MC4R: n 31 of 177, best F1 0.703 (expected 0.69 +/- 0.069) on [sag_tau, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a'), preprocessing (b), layer code, PAX6: n 32 of 177, best F1 0.892 (expected 0.89 +/- 0.089) on [fast_trough_deltav_ramp, sag_peak_t] over 3486 subspaces; within tolerance: True features (a'), preprocessing (b), layer code, all others: n 28 of 177, best F1 0.519 (expected 0.5 +/- 0.050) on [upstroke_downstroke_ratio_rheo, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a, ratios), preprocessing (b), layer code, LAMP5: n 86 of 177, best F1 0.825 (expected 0.81 +/- 0.081) on [threshold_i_ramp, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a, ratios), preprocessing (b), layer code, MC4R: n 31 of 177, best F1 0.694 (expected 0.69 +/- 0.069) on [sag_tau, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True features (a, ratios), preprocessing (b), layer code, PAX6: n 32 of 177, best F1 0.892 (expected 0.89 +/- 0.089) on [fast_trough_deltav_ramp, sag_peak_t] over 3486 subspaces; within tolerance: True features (a, ratios), preprocessing (b), layer code, all others: n 28 of 177, best F1 0.511 (expected 0.5 +/- 0.050) on [upstroke_downstroke_ratio_rheo, upstroke_adapt_ratio] over 3486 subspaces; within tolerance: True
Result: reproduced. Under both preprocessing readings and both feature readings the best one- or two-feature mixture gives F1 within tolerance for every group. On the features computed here, the Methods' preprocessing gives LAMP5 0.785, MC4R 0.709, PAX6 0.896 and all others 0.490, and the code's preprocessing 0.806, 0.703, 0.892 and 0.526. On the paper's values alone the two give LAMP5 0.802 to 0.806, MC4R 0.693 to 0.712, PAX6 0.892 to 0.896 and all others 0.495 to 0.506. The closest to a bound is all others at 0.490 (floor 0.45). With the layer filled from the LIMS layer (189 cells) the values are 0.840, 0.707, 0.879 and 0.494; in variant (a') 0.806, 0.703, 0.892 and 0.519; with the adaptation ratios computed here 0.825, 0.694, 0.892 and 0.511; all within tolerance. Several of the best features (sag_tau, sag_peak_t, ramp and short square features) are the paper's values, as they are not computed here.
figs8a-classifier-cv¶
A multi-class classifier assigns human L1 subclasses from electrophysiology with 82% cross-validated accuracy, balanced across subclasses.
As clustering_classification.ipynb cells 26 to 35 and 44 to 45: cells of the strict L1 set in the four core subclasses with soma in L1, AIBS; logistic regression with balanced class weights after RobustScaler and KNN imputation (O3 (a) adds the LocalOutlierFactor step); repeated stratified 5-fold cross-validation, 10 repeats (seed 0). Features (a) more than 50% non-missing at every site, (b) the site-aligned subset. Choosing features by site needs the Amsterdam and Szeged feature tables of the paper's code (tamas_features_E.csv, mansvelder_features_E.csv); they are used only for this choice. Soma in L1 as cell 26 has it (depth table layer, not filled); the filled layer and feature reading (b) are variants.
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import confusion_matrix
from sklearn.model_selection import RepeatedStratifiedKFold
from sklearn.preprocessing import RobustScaler
sites_b = pd.concat([sites[sites.collaborator != "AIBS"], sites[sites.collaborator == "AIBS"].drop(columns=fcols).join(pap[fcols])], sort=False)
sites_b["homology_type"] = pd.Categorical(sites_b["t-type"].map(L.HOMOLOGY), categories=L.SUBCLASSES + ["other"])
def core_l1(S):
return S[S.in_strict_set & S.homology_type.isin(L.SUBCLASSES) & (S.failed_fx_long_squares == False) & (S.layer == "1")]
inh = core_l1(sites)
print("Cells by site (strict set, core subclasses, soma in L1):", inh.collaborator.value_counts().to_dict())
common = [c for c in fcols if inh.groupby("collaborator")[c].apply(lambda x: x.notna().mean()).min() > 0.5]
site_fit = []
for c in common:
fit = smf.ols(f"{c} ~ collaborator", data=inh).fit(cov_type="HC3")
site_fit.append({"feature": c, "rsquared": fit.rsquared, "p": fit.f_pvalue})
site_fit = pd.DataFrame(site_fit).set_index("feature")
aligned = list(site_fit[(site_fit.rsquared < 0.05) & (site_fit.p > 0.05)].index)
print(f"features (a) common to every site: {len(common)}; (b) site-aligned: {len(aligned)}")
def cv_scores(X, y, outliers):
steps = [("norm", RobustScaler()), ("impute", KNNImputer())]
if outliers:
steps = [("pre", L.preproc_pipeline())]
pipe = Pipeline(steps + [("class", LogisticRegression(class_weight="balanced"))])
cv = RepeatedStratifiedKFold(n_splits=5, n_repeats=10, random_state=0)
acc, cms = [], []
for tr, te in cv.split(X, y):
pipe.fit(X[tr], y[tr]); pred = pipe.predict(X[te])
acc.append(np.mean(pred == y[te])); cms.append(confusion_matrix(y[te], pred, labels=range(4)))
cm = np.mean(cms, axis=0); cm = cm / cm.sum(axis=1, keepdims=True)
return np.mean(acc), np.std(acc), np.mean(np.diag(cm)), cm
# Soma in L1 as clustering_classification.ipynb cell 26 has it (depth table layer, not filled);
# the layer filled from the LIMS layer is a variant. Feature reading (b) uses the paper's values
# for the Seattle cells.
for feat_kind, S, layer_rule in [("a", sites, "code"), ("a'", sites_ns, "code"), ("b", sites_b, "code"), ("a", sites.assign(
layer=sites.layer.fillna(sites.layer_lims.astype("string").str.replace(".0", "", regex=False))), "filled")]:
inh_k = core_l1(S)
data = inh_k[inh_k.collaborator == "AIBS"]
y = data.homology_type.cat.codes.values
print(f"feature reading ({feat_kind}), layer {layer_rule}: AIBS cells {len(data)}", data.homology_type.value_counts(sort=False).to_dict())
for fname, fl in [("a common", common), ("b aligned", aligned)]:
for prep, outl in [("code", False), ("Methods + LOF", True)]:
acc, sd, bal, cm = cv_scores(data[fl].values, y, outl)
print(f" features {fname}, preprocessing {prep}: accuracy {100 * acc:.1f}% (sd {100 * sd:.1f}), balanced {100 * bal:.1f}%; "
f"per class recall {np.round(np.diag(cm), 2).tolist()}")
print(" chance (largest class):", round(100 * pd.Series(y).value_counts(normalize=True).max(), 1), "%; balanced chance 25%")
Cells by site (strict set, core subclasses, soma in L1): {'AIBS': 164, 'Gabor': 17, 'Huib': 9}
features (a) common to every site: 72; (b) site-aligned: 35
feature reading (a), layer code: AIBS cells 164 {'LAMP5': 86, 'MC4R': 31, 'PAX6': 32, 'L1 VIP': 15, 'other': 0}
features a common, preprocessing code: accuracy 84.6% (sd 6.0), balanced 78.8%; per class recall [0.92, 0.76, 0.84, 0.63]
features a common, preprocessing Methods + LOF: accuracy 83.5% (sd 6.2), balanced 77.1%; per class recall [0.92, 0.72, 0.83, 0.61]
features b aligned, preprocessing code: accuracy 78.8% (sd 7.0), balanced 74.5%; per class recall [0.84, 0.78, 0.72, 0.63]
features b aligned, preprocessing Methods + LOF: accuracy 79.7% (sd 6.6), balanced 74.6%; per class recall [0.86, 0.76, 0.76, 0.61]
chance (largest class): 52.4 %; balanced chance 25%
feature reading (a'), layer code: AIBS cells 164 {'LAMP5': 86, 'MC4R': 31, 'PAX6': 32, 'L1 VIP': 15, 'other': 0}
features a common, preprocessing code: accuracy 85.0% (sd 5.6), balanced 79.0%; per class recall [0.92, 0.75, 0.86, 0.63]
features a common, preprocessing Methods + LOF: accuracy 84.1% (sd 5.7), balanced 77.4%; per class recall [0.92, 0.72, 0.85, 0.6]
features b aligned, preprocessing code: accuracy 79.4% (sd 6.7), balanced 75.5%; per class recall [0.84, 0.78, 0.76, 0.64]
features b aligned, preprocessing Methods + LOF: accuracy 80.1% (sd 6.5), balanced 75.4%; per class recall [0.85, 0.76, 0.78, 0.62]
chance (largest class): 52.4 %; balanced chance 25%
feature reading (b), layer code: AIBS cells 164 {'LAMP5': 86, 'MC4R': 31, 'PAX6': 32, 'L1 VIP': 15, 'other': 0}
features a common, preprocessing code: accuracy 84.0% (sd 6.0), balanced 77.2%; per class recall [0.92, 0.73, 0.85, 0.59]
features a common, preprocessing Methods + LOF: accuracy 83.5% (sd 6.0), balanced 76.8%; per class recall [0.91, 0.75, 0.82, 0.58]
features b aligned, preprocessing code: accuracy 79.6% (sd 7.0), balanced 75.5%; per class recall [0.84, 0.77, 0.78, 0.63]
features b aligned, preprocessing Methods + LOF: accuracy 80.0% (sd 7.3), balanced 75.9%; per class recall [0.84, 0.79, 0.79, 0.61]
chance (largest class): 52.4 %; balanced chance 25%
feature reading (a), layer filled: AIBS cells 176 {'LAMP5': 95, 'MC4R': 34, 'PAX6': 32, 'L1 VIP': 15, 'other': 0}
features a common, preprocessing code: accuracy 84.2% (sd 5.9), balanced 77.7%; per class recall [0.92, 0.73, 0.84, 0.62]
features a common, preprocessing Methods + LOF: accuracy 83.1% (sd 5.5), balanced 76.2%; per class recall [0.91, 0.7, 0.84, 0.59]
features b aligned, preprocessing code: accuracy 80.8% (sd 7.4), balanced 76.1%; per class recall [0.86, 0.77, 0.77, 0.64]
features b aligned, preprocessing Methods + LOF: accuracy 80.1% (sd 7.4), balanced 74.9%; per class recall [0.86, 0.73, 0.77, 0.63] chance (largest class): 54.0 %; balanced chance 25%
Result: reproduced. Cross-validated accuracy is 84.6% with the features common to all sites and the code's preprocessing, within 82 +/- 8.2; every other combination of features, preprocessing, feature reading (variant a' included) and layer rule gives 78.8 to 85.0%, also within it. Balanced accuracy is 74.5 to 79.0%. The largest class is 52.4% of cells (54.0% with the filled layer), and balanced chance is 25%; the paper states no chance level.
figs8b-classifier-holdout¶
Trained on the primary dataset with site-dependent features excluded, the classifier reaches 81% accuracy on held-out cells recorded at other sites.
Blocked: the held-out Amsterdam and Szeged cells are not in this dandiset. The cell below checks two other DANDI datasets with Szeged (Tamas lab) and Amsterdam (Mansvelder lab) human recordings, read from the archive's metadata only; neither is named by the paper or the claims file.
# DANDI 000228 (Szeged files, subjects named by donor) and 000337 (Amsterdam files, subjects and
# sessions named by LIMS ids), at their latest published versions; metadata and asset paths only.
held = core_l1(sites)
held = held[held.collaborator != "AIBS"]
print("Held-out cells (strict set, core subclasses, soma in L1, features extracted) by site:", held.collaborator.value_counts().to_dict())
for did, site in [("000228", "Gabor"), ("000337", "Huib")]:
ver = requests.get(f"{L.API}/dandisets/{did}/", timeout=60).json()["most_recent_published_version"]["version"]
m = L.version_metadata(did, ver)
a = L.list_assets(did, ver)
subj = a.path.str.extract(r"sub-([^/_]+)")[0]
ses = pd.to_numeric(a.path.str.extract(r"ses-(\d+)")[0], errors="coerce")
h = held[held.collaborator == site]
by_donor = h.donor_name.isin(set(subj.str.replace("-", ".")))
by_session = h.ephys_roi_id.notna() & h.ephys_roi_id.isin(set(ses.dropna()))
print(f"DANDI {did} version {ver} ({m['name']!r}; contributors {[c.get('name') for c in m.get('contributor', []) if c.get('schemaKey') == 'Person']}):")
print(f" {len(a)} files, {subj.nunique()} subjects, {int(ses.notna().sum())} files with a numeric session id; "
f"session part of the paths: {a.path.str.extract(r'ses-([^/_]+)')[0].fillna('none').value_counts().head(3).to_dict()}")
print(f" held-out {site} cells: {len(h)}; with their donor among the subjects: {int(by_donor.sum())}; "
f"with their LIMS session id among the files: {int(by_session.sum())}; with a LIMS id in the table: {int(h.ephys_roi_id.notna().sum())}")
Held-out cells (strict set, core subclasses, soma in L1, features extracted) by site: {'Gabor': 17, 'Huib': 9}
DANDI 000228 version 0.260829.2040 ('20220330_AIBS_Patchseq_human'; contributors ['Chartrand, Thomas', 'Kalmbach, Brian', 'Molnar, Gabor', 'Tamas, Gabor', 'Lein, Ed']):
91 files, 32 subjects, 0 files with a numeric session id; session part of the paths: {'none': 91}
held-out Gabor cells: 17; with their donor among the subjects: 15; with their LIMS session id among the files: 0; with a LIMS id in the table: 17
DANDI 000337 version 0.260917.2218 ('20220917_AIBS_Patchseq_human'; contributors ['Chartrand, Thomas', 'Kalmbach, Brian', 'Mansvelder, Huib', 'Lein, Ed']):
21 files, 12 subjects, 13 files with a numeric session id; session part of the paths: {'PLACEHOLDER': 8, '1053917650': 1, '1053917889': 1}
held-out Huib cells: 9; with their donor among the subjects: 0; with their LIMS session id among the files: 0; with a LIMS id in the table: 5
Result: blocked, data missing. The 26 held-out cells (17 Szeged, 9 Amsterdam) are not in this dandiset. DANDI 000228 (published 2026) holds Szeged files from the donors of 15 of the 17 Szeged cells, but its file paths carry no session id and the cells cannot be matched without matching features. DANDI 000337 holds no file with the LIMS session id of an Amsterdam held-out cell (5 of the 9 have one in the table); 8 of its 21 file paths carry the session id PLACEHOLDER. Neither dataset is named by the paper or the claims file; whether the claims file should name 000228 is asked of a reader.
Summary¶
| Claim | Outcome |
|---|---|
| dandiset-cell-link | partially reproduced |
| patchseq-counts | reproduced |
| fig2c-subclass-ephys-differences | reproduced, resting on an assumption |
| fig2b-subclass-ephys-profiles | reproduced |
| figs3a-sag-adaptation | not reproduced |
| fig2c-lamp5-gradient | not reproduced |
| figs6-donor-condition-ephys | reproduced |
| figs6c-region-ephys | reproduced |
| fig4b-mc4r-sag-irregularity | reproduced |
| fig4b-adarb2-irregular-firing | not reproduced |
| fig5b-tnfaip8l3-ephys | reproduced |
| fig6a-ephys-dprime-subclass-structure | not reproduced |
| fig6d-gmm-subclass-f1 | reproduced |
| figs8a-classifier-cv | reproduced |
| figs8b-classifier-holdout | blocked, data missing |
The mouse claims are in mouse_comparison.ipynb and four morphology claims in
morphology_tables.ipynb; the record of every claim is in attempts/000630-chartrand-2023-01.yaml.
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.
Changed after the third round of verification: fig2c-subclass-ephys-differences to reproduced, resting on the assumption that the paper's QC did not apply the 500 ms noise criterion as measured here (the paper's own values show its pipeline kept those sweeps; with the criterion as written the claim is partially reproduced), and fig2c-lamp5-gradient and fig6a-ephys-dprime-subclass-structure to not reproduced (their criteria name no partial level).
How the analysis differs from the paper's: the features are computed here with IPFX 2.1.2 from the long square sweeps, with the Methods' sweep QC as far as the files allow; the 39 features used agree with the paper's table within 1% for 53 to 98% of cells. The custom features the paper did not release, and the adaptation ratios, which agree for 4 to 37% of cells, are the paper's values. Where a feature at rheobase + 40 pA differs, the sweep chosen here is often another repeat or another amplitude than the paper's, or a lower sweep because the one in range fails the 500 ms noise criterion here. Cell labels come from the paper's table.
What was found about the dataset: it holds exactly the paper's primary set of 210 Seattle cells,
one per file, linked to the paper's table only through the LIMS number in session_id; the files
carry no t-type, subclass, layer or region, 6 files carry a time stamp in place of the cell name,
one file states sex unknown and two state an age of 0 days where the table differs, and the 50
Amsterdam and Szeged cells of the strict L1 set are not in it. Unacquired samples of sweeps that stopped
early are stored as zeros, and bias_current is stored in pA under a schema that fixes amperes.
print(f"notebook run time {(time.time() - started) / 60:.1f} min")
notebook run time 36.0 min