import os
os.environ.setdefault("OMP_NUM_THREADS", "1")
import time
import math
import itertools
import numpy as np
import pandas as pd
import scipy.stats as ss
from scipy.spatial.distance import pdist
import matplotlib
matplotlib.use("Agg")
from IPython.display import display
import scala

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

dandiset 000008 version 0.211014.0809

assets = scala.list_assets()
t = time.time()
meta, feats, stim, binned, errors = scala.load_all(assets, workers=8)
print(f"read and extracted {len(meta)} files in {(time.time() - t) / 60:.1f} min; errors: {len(errors)}")
tab = scala.paper_meta_table()
excluded = set(tab.index[tab["Exclusion reasons"].notna()])
dups = scala.duplicate_recordings(meta)
dropped = list(dups.cell_id[~dups.kept])
t = scala.cell_table(meta, feats, excluded, drop=dropped)       # each recording once
t_paper = scala.cell_table(meta, feats, excluded)                # all cells, as the paper
X17 = scala.features17(t)
S = scala.standardized17(X17)
S_paper = scala.standardized17(scala.features17(t_paper))
ex, it, gl = scala.patchseq_counts()
ref = scala.reference("10X_cells_v2_AIBS")
names, _ = scala.reference_names()
fams = scala.family_of_type(names)
rnaseqTools = scala.paper_module("rnaseqTools")

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

print("dandiset cells:", len(t), "| mice:", t.mouse.nunique(), "| with all 17 features:", len(S),
      "| with a t-type and not excluded:", int(((t.ttype != "") & ~t.excluded).sum()), "| excluded (table):", int(t.excluded.sum()))
print("feature extraction errors:", len(errors))
print("duplicated recordings (see scala_2021_ephys.ipynb); features left out for:", {c: t.loc[c, "ttype"] for c in dropped},
      "| t-types of both members:", sorted(set(t.loc[dups.cell_id, "ttype"])))
print("count tables: genes x cells", ex.shape, "| cells of the dandiset in them:", int(t.index.isin(ex.columns).sum()))
gd = pd.to_numeric(meta.set_index("cell_id").number_of_genes_detected)
det_ = ((ex[gd.index].values + it[gd.index].values) > 0).sum(axis=0)
print("'number_of_genes_detected' in the files equal to the genes detected in the count tables:", int((gd.values == det_).sum()), "of", len(gd))
for k in ["neurons", "pvsst", "viplamp", "exc"]:
    print(f"10x v2 reference '{k}': cells x genes {ref[k]['counts'].shape}, t-types present {len(np.unique(ref[k]['clusters']))}")
print("reference t-type names (updated as the code does):", len(names), "| dandiset t-types not among them:",
      sorted(set(t.ttype) - set(names) - {""}))
for k in ["viplamp", "pvsst", "exc"]:
    Zk = scala.tsne_atlas(k)
    print(f"t-SNE atlas '{k}': {Zk.shape}, rows match the reference cells: {Zk.shape[0] == ref[k]['counts'].shape[0]}, "
          f"x {Zk[:, 0].min():.0f} to {Zk[:, 0].max():.0f}, y {Zk[:, 1].min():.0f} to {Zk[:, 1].max():.0f}")
print("features (standardized 17):", list(S.columns))

dandiset cells: 1328 | mice: 266 | with all 17 features: 1317 | with a t-type and not excluded: 1226 | excluded (table): 7
feature extraction errors: 0
duplicated recordings (see scala_2021_ephys.ipynb); features left out for: {'20190731_sample_8': 'Pvalb Il1rapl2', '20181116_sample_4': 'Pvalb Gpr149', '20190829_sample_6': 'Lamp5 Slc35d3'} | t-types of both members: ['Lamp5 Slc35d3', 'Pvalb Gpr149', 'Pvalb Il1rapl2']
count tables: genes x cells (42466, 1329) | cells of the dandiset in them: 1328

'number_of_genes_detected' in the files equal to the genes detected in the count tables: 1328 of 1328
10x v2 reference 'neurons': cells x genes (121423, 1000), t-types present 87
10x v2 reference 'pvsst': cells x genes (12083, 500), t-types present 29
10x v2 reference 'viplamp': cells x genes (15511, 500), t-types present 27
10x v2 reference 'exc': cells x genes (93829, 500), t-types present 31
reference t-type names (updated as the code does): 88 | dandiset t-types not among them: []

t-SNE atlas 'viplamp': (15511, 2), rows match the reference cells: True, x -68 to 66, y -58 to 63
t-SNE atlas 'pvsst': (12083, 2), rows match the reference cells: True, x -69 to 68, y -59 to 70

t-SNE atlas 'exc': (93829, 2), rows match the reference cells: True, x -65 to 60, y -63 to 67
features (standardized 17): ['Resting membrane potential (mV)', 'Input resistance (MOhm)', 'Membrane time constant (ms)', 'AP threshold (mV)', 'AP amplitude (mV)', 'AP width (ms)', 'Upstroke-to-downstroke ratio', 'Afterhyperpolarization (mV)', 'ISI adaptation index', 'Max number of APs', 'Rheobase (pA)', 'Sag ratio', 'Latency (ms)', 'ISI coefficient of variation', 'Rebound (mV)', 'AP amplitude adaptation index', 'AP coefficient of variation']

sparseRRR = scala.paper_module("sparseRRR")
import scipy.sparse as sp

def standardize(A):
    return (A - A.mean(axis=0)) / A.std(axis=0)

def rrr_data(tt):
    order_cells = [c for c in tab.index if c in tt.index]   # the paper's cell order
    F = scala.features17(tt).reindex(order_cells)
    sel = F.notna().all(axis=1) & (tt.loc[order_cells, "ttype"] != "") & ~tt.loc[order_cells, "excluded"]
    cells_ = list(F.index[sel])
    exs_, its_ = ex[cells_].values.T.astype(float), it[cells_].values.T.astype(float)
    genes_ = rnaseqTools.geneSelection(sp.csr_matrix(exs_), n=1000, threshold=32, decay=1.5, plot=False)
    Xg_ = exs_ / gl.exon_bp.values + its_ / (gl.intron_bp.values + 0.001)
    Xg_ = Xg_ / Xg_.sum(axis=1, keepdims=True) * 1e6
    return cells_, exs_, its_, genes_, Xg_

cells_rrr, exs, its, genes_rrr, Xg = rrr_data(t)
print("cells:", len(cells_rrr), "from", t.loc[cells_rrr, "mouse"].nunique(), "mice")

def rrr_run(X, logs, lams, cells_=None):
    cells_ = cells_rrr if cells_ is None else cells_
    Y = scala.features17(t_paper.loc[cells_], log=logs).values
    Y = standardize(Y)
    out = []
    for lam in lams:
        w, v = sparseRRR.relaxed_elastic_rrr(X, Y, rank=5, lambdau=lam, alpha=1)
        r = sparseRRR.elastic_rrr_cv(X, Y, rank=5, reps=1, folds=10, alphas=np.array([1]), lambdas=np.array([lam]))
        corr = np.nanmean(r[4], axis=0).squeeze()
        out.append({"lambda": lam, "genes (full data)": int(np.sum(w[:, 0] != 0)), "genes (mean over folds)": float(np.nanmean(r[2])),
                    "R2": float(np.nanmean(r[1])), "corr 1": corr[0], "corr 2": corr[1], "corr 3": corr[2], "corr 4": corr[3], "corr 5": corr[4]})
    return pd.DataFrame(out)

X = standardize(np.log2(Xg[:, genes_rrr] + 1))
print("X:", X.shape)
rrr_main = rrr_run(X, scala.FEATURES_LOG, [0.5, 0.45])
rrr_code = rrr_run(X, scala.FEATURES_LOG[:3], [0.5, 0.45])
print("main reading (four features log-transformed, Methods):"); display(rrr_main.round(3))
# The folds as sparseRRR.elastic_rrr_cv builds them (reps=1): np.random.seed(42), one permutation of the cells,
# then 10 contiguous blocks of int(n / 10) permuted positions as test sets; the remainder is never a test cell.
n_cv = len(cells_rrr)
np.random.seed(42)
perm = np.random.permutation(n_cv)
fs = int(n_cv / 10)
print(f"folds: {n_cv} cells, test blocks of {fs}, cells never in a test block: {n_cv - 10 * fs} {perm[10 * fs:].tolist()}")
for k_ in range(10):
    print(f"  fold {k_}: test cells (indices in the paper's cell order) {perm[k_ * fs:(k_ + 1) * fs][:6].tolist()} ...")
print("variant (three log-transformed, as the code's RRR notebook):"); display(rrr_code.round(3))
w, v = sparseRRR.relaxed_elastic_rrr(X, standardize(scala.features17(t.loc[cells_rrr]).values), rank=5, lambdau=0.5, alpha=1)
print("genes selected at lambda 0.5:", list(ex.index[genes_rrr][np.where(w[:, 0] != 0)[0]]))
exp = np.array([0.38, 0.90, 0.74, 0.67]); tol = np.array([0.038, 0.09, 0.074, 0.067])
cp, _, _, gp, Xgp = rrr_data(t_paper)
rrr_paper = rrr_run(standardize(np.log2(Xgp[:, gp] + 1)), scala.FEATURES_LOG, [0.5, 0.45], cells_=cp)
print(f"all cells, as the paper ({len(cp)} cells), main reading:"); display(rrr_paper.round(3))
for label, df in [("main", rrr_main), ("variant", rrr_code), ("all cells, main", rrr_paper)]:
    for _, row in df.iterrows():
        vals = np.array([row.R2, row["corr 1"], row["corr 2"], row["corr 3"]])
        print(f"{label}, lambda {row['lambda']}: within tolerance {[bool(x) for x in np.abs(vals - exp) <= tol]}")

Chosen offset: 7.59

cells: 1216 from 264 mice
X: (1216, 1000)

1.

.

.

.

.

.

.

.

import io
ion = pd.read_csv(io.StringIO(scala._get(f"{scala.MINI_ATLAS}/data/ion-channel-genes-group-177.csv").text))["Approved symbol"].values
genes_ion = np.array([g.upper() in ion for g in ex.index]) & (((exs + its) > 0).sum(axis=0) >= 10)
print("ion channel genes in the list:", len(ion), "| detected in at least 10 of these cells:", int(genes_ion.sum()))
Xi = standardize(np.log2(Xg[:, genes_ion] + 1))
ion_main = rrr_run(Xi, scala.FEATURES_LOG, [0.303])
ion_code = rrr_run(Xi, scala.FEATURES_LOG[:3], [0.303])
print("main reading:"); display(ion_main.round(3))
print("variant (three log-transformed):"); display(ion_code.round(3))
exp = np.array([0.33, 0.86, 0.71, 0.56]); tol = np.array([0.033, 0.086, 0.071, 0.056])
for label, df in [("main", ion_main), ("variant", ion_code)]:
    row = df.iloc[0]
    vals = np.array([row.R2, row["corr 1"], row["corr 2"], row["corr 3"]])
    print(f"{label}: within tolerance {[bool(x) for x in np.abs(vals - exp) <= tol]}")

ion channel genes in the list: 328 | detected in at least 10 of these cells: 307

1.

.

.

.

.

.

.

.

.

def mean_rna(subset, min_cells=10):
    out = {}
    cl = ref[subset]["clusters"]
    for c in range(len(names)):
        ind = cl == c
        if ind.sum() >= min_cells:
            out[names[c]] = np.asarray(np.log1p(ref[subset]["counts"][ind].toarray()).mean(axis=0)).ravel()
    return out

rna = {k: mean_rna(k) for k in ["neurons", "pvsst", "viplamp", "exc"]}


def pooled_within(tt, SS, drop):
    base = (tt.ttype != "") & ~tt.excluded & tt.index.isin(SS.index)
    sel = base & ~tt.ttype.isin(["Sst Chodl", "Pvalb Vipr2_2"]) if drop else base
    cnt = tt[sel].ttype.value_counts()
    types = [n for n in names if cnt.get(n, 0) >= 5 and n in rna["neurons"]]
    P = np.array([SS.loc[tt[sel & (tt.ttype == n)].index].mean().values for n in types])
    R = np.array([rna["neurons"][n] for n in types])
    Dt, Dph = pdist(R, "correlation"), pdist(P)
    fam_codes = np.unique(fams[[list(names).index(n) for n in types]], return_inverse=True)[1]
    within = pdist(fam_codes[:, None]) == 0
    r = ss.pearsonr(Dt[within], Dph[within])
    print(f"t-types {len(types)}, within-family pairs {int(within.sum())}, r = {r.statistic:.3f} (p = {r.pvalue:.2g}); "
          f"all pairs r = {np.corrcoef(Dt, Dph)[0, 1]:.3f}; t-types with at least 5 cells but fewer than 10 reference cells: "
          f"{[n for n in cnt.index if cnt[n] >= 5 and n not in rna['neurons']]}; transcriptomic distances "
          f"{Dt[within].min():.3f} to {Dt[within].max():.3f}, electrophysiological {Dph[within].min():.2f} to {Dph[within].max():.2f}")
    return r.statistic


for tlabel, tt, SS in [("each recording once", t, S), ("all cells, as the paper", t_paper, S_paper)]:
    for reading, drop in [("(a) code: Sst Chodl and Pvalb Vipr2_2 left out", True), ("(b) both kept", False)]:
        print(f"{tlabel}, {reading}:", end=" ")
        pooled_within(tt, SS, drop)

each recording once, (a) code: Sst Chodl and Pvalb Vipr2_2 left out: t-types 50, within-family pairs 200, r = 0.597 (p = 9.9e-21); all pairs r = 0.633; t-types with at least 5 cells but fewer than 10 reference cells: []; transcriptomic distances 0.042 to 0.339, electrophysiological 0.73 to 4.77
each recording once, (b) both kept: t-types 52, within-family pairs 222, r = 0.563 (p = 6e-20); all pairs r = 0.622; t-types with at least 5 cells but fewer than 10 reference cells: []; transcriptomic distances 0.042 to 0.418, electrophysiological 0.73 to 4.77
all cells, as the paper, (a) code: Sst Chodl and Pvalb Vipr2_2 left out:

t-types 50, within-family pairs 200, r = 0.597 (p = 1.1e-20); all pairs r = 0.633; t-types with at least 5 cells but fewer than 10 reference cells: []; transcriptomic distances 0.042 to 0.339, electrophysiological 0.73 to 4.77
all cells, as the paper, (b) both kept: t-types 52, within-family pairs 222, r = 0.562 (p = 6.5e-20); all pairs r = 0.622; t-types with at least 5 cells but fewer than 10 reference cells: []; transcriptomic distances 0.042 to 0.418, electrophysiological 0.73 to 4.77

group = {"Pvalb": "pvsst", "Sst": "pvsst", "Vip": "viplamp", "Lamp5": "viplamp", "IT": "exc", "CT": "exc"}

def within_families(tt, SS, show=False):
    rows, feat_r = [], {}
    for fam_ in ["Pvalb", "Sst", "Vip", "Lamp5", "IT", "CT"]:
        sel = (tt.family == fam_) & ~tt.excluded & tt.index.isin(SS.index) & ~tt.ttype.isin(["Sst Chodl", "Pvalb Vipr2_2"])
        cnt = tt[sel].ttype.value_counts()
        types = [n for n in names if cnt.get(n, 0) >= 5 and n in rna[group[fam_]]]
        row = {"family": fam_, "t-types": len(types), "pairs": len(types) * (len(types) - 1) // 2}
        if len(types) >= 3:
            P = np.array([SS.loc[tt[sel & (tt.ttype == n)].index].mean().values for n in types])
            R = np.array([rna[group[fam_]][n] for n in types])
            Dt = pdist(R, "correlation")
            row["r"] = ss.pearsonr(Dt, pdist(P)).statistic
            feat_r[fam_] = pd.Series([ss.pearsonr(Dt, pdist(P[:, [i]])).statistic for i in range(17)], index=SS.columns)
        rows.append(row)
    fam_tab = pd.DataFrame(rows).set_index("family")
    FR = pd.DataFrame(feat_r)
    if show:
        display(fam_tab.round(3))
        display(FR.round(2))
    counted = fam_tab.dropna(subset=["r"])
    npos = int((counted.r > 0).sum())
    many = int(((FR > 0).sum(axis=1) >= len(FR.columns) / 2).sum())
    print(f"families counted (at least 3 t-types): {len(counted)}; with r > 0: {npos}; reading (a) at least two: {npos >= 2}; "
          f"reading (b) at least four: {npos >= 4}; features with r > 0 in at least half of the {len(FR.columns)} counted "
          f"families: {many} of 17 (many: {many >= 8.5})")


for tlabel, tt, SS in [("each recording once", t, S), ("all cells, as the paper", t_paper, S_paper)]:
    print(f"--- {tlabel}")
    within_families(tt, SS, show=(tlabel == "each recording once"))

--- each recording once

families counted (at least 3 t-types): 6; with r > 0: 6; reading (a) at least two: True; reading (b) at least four: True; features with r > 0 in at least half of the 6 counted families: 16 of 17 (many: True)
--- all cells, as the paper

families counted (at least 3 t-types): 6; with r > 0: 6; reading (a) at least two: True; reading (b) at least four: True; features with r > 0 in at least half of the 6 counted families: 16 of 17 (many: True)

um = (pd.to_numeric(tab["Soma depth (µm)"], errors="coerce") / pd.to_numeric(tab["Cortical thickness (µm)"], errors="coerce").replace(0, np.nan)).reindex(t.index)
for label, dn in [("files (4x)", t.depth_norm), ("paper's table (µm)", um)]:
    sel = (t.family == "IT") & ~t.excluded & dn.notna()
    cnt = t[sel].ttype.value_counts()
    types = [n for n in names if cnt.get(n, 0) >= 5]
    dep = np.array([dn[t[sel & (t.ttype == n)].index].mean() for n in types])
    R = np.array([rna["exc"][n] for n in types])
    r = ss.pearsonr(pdist(R, "correlation"), pdist(dep[:, None]))
    print(f"{label}: IT cells with depth {int(sel.sum())}, t-types {len(types)} {types}; r = {r.statistic:.3f} (p = {r.pvalue:.2g}); "
          f"mean depths {np.round(dep, 3).tolist()}")

files (4x): IT cells with depth 249, t-types 8 ['L2/3 IT_3', 'L4/5 IT_1', 'L4/5 IT_2', 'L5 IT_1', 'L5 IT_2', 'L5 IT_3', 'L6 IT_1', 'L6 IT_2']; r = 0.690 (p = 4.9e-05); mean depths [0.185, 0.282, 0.368, 0.453, 0.525, 0.548, 0.725, 0.827]
paper's table (µm): IT cells with depth 249, t-types 8 ['L2/3 IT_3', 'L4/5 IT_1', 'L4/5 IT_2', 'L5 IT_1', 'L5 IT_2', 'L5 IT_3', 'L6 IT_1', 'L6 IT_2']; r = 0.698 (p = 3.6e-05); mean depths [0.185, 0.281, 0.368, 0.451, 0.515, 0.54, 0.772, 0.822]

def position(cells_, subset):
    Zr = scala.tsne_atlas(subset)
    pos = rnaseqTools.map_to_tsne(ref[subset]["counts"], ref[subset]["genes"], ex[cells_].values.T.astype(float), ex.index.values,
                                  Zr, knn=10, newIntronCounts=it[cells_].values.T.astype(float),
                                  newExonLengths=gl.exon_bp.values, newIntronLengths=gl.intron_bp.values,
                                  normalizeNew=True, normalizeReference=False)
    return pos, Zr

from scipy.spatial.distance import correlation
sel = (t.family == "Vip") & ~t.excluded & t.index.isin(S.index)
vip = list(t[sel].index)
pos, Zr = position(vip, "viplamp")
refn = names[ref["viplamp"]["clusters"]]
sncg = np.array([n.startswith("Sncg") for n in refn])
gpc3 = refn == "Vip Gpc3"
c0, c1 = Zr[sncg].mean(axis=0), Zr[gpc3].mean(axis=0)
u = (c1 - c0) / np.linalg.norm(c1 - c0)
proj = (pos - c0) @ u
tau = t.loc[vip, "Membrane time constant (ms)"].values
ra = ss.spearmanr(tau, proj)
pa = ss.spearmanr(tau, proj, alternative="less").pvalue
print(f"Vip cells: {len(vip)} from {t.loc[vip, 'mouse'].nunique()} mice; tau {tau.min():.1f} to {tau.max():.1f} ms; "
      f"projection {proj.min():.1f} to {proj.max():.1f} (Sncg centroid 0, Vip Gpc3 centroid {np.linalg.norm(c1 - c0):.1f})")
print(f"(a) Spearman rho(tau, position towards Vip Gpc3) = {ra.statistic:.3f}, one-sided p (negative) = {pa:.2g}")
L = np.log1p(ref["viplamp"]["counts"].toarray())
sn_cells = L[sncg].mean(axis=0)
sn_types = np.mean([L[refn == n].mean(axis=0) for n in np.unique(refn[sncg])], axis=0)
cnt = t[sel].ttype.value_counts()
types = [n for n in names if cnt.get(n, 0) >= 5]
med = [np.median(t.loc[t[sel & (t.ttype == n)].index, "Membrane time constant (ms)"]) for n in types]
for label, target in [("mean over Sncg reference cells", sn_cells), ("mean of the Sncg t-type means", sn_types)]:
    dist = [correlation(L[refn == n].mean(axis=0), target) for n in types]
    rb = ss.spearmanr(med, dist)
    print(f"(b) {label}: {len(types)} Vip t-types; Spearman rho(median tau, distance to Sncg) = {rb.statistic:.3f}, "
          f"one-sided p (negative) = {ss.spearmanr(med, dist, alternative='less').pvalue:.3f}")
    print("    ", dict(zip(types, [f"tau {m:.1f}, distance {d_:.3f}" for m, d_ in zip(med, dist)])))

Using a common set of 490 genes.

Vip cells: 153 from 51 mice; tau 3.1 to 29.7 ms; projection 6.2 to 95.3 (Sncg centroid 0, Vip Gpc3 centroid 91.1)
(a) Spearman rho(tau, position towards Vip Gpc3) = -0.582, one-sided p (negative) = 1.6e-15
(b) mean over Sncg reference cells: 9 Vip t-types; Spearman rho(median tau, distance to Sncg) = -0.383, one-sided p (negative) = 0.154
     {'Vip Sncg': 'tau 19.0, distance 0.119', 'Vip Serpinf1_1': 'tau 8.8, distance 0.299', 'Vip Htr1f': 'tau 10.0, distance 0.360', 'Vip Gpc3': 'tau 5.6, distance 0.441', 'Vip C1ql1': 'tau 5.5, distance 0.403', 'Vip Mybpc1_1': 'tau 11.8, distance 0.438', 'Vip Mybpc1_2': 'tau 14.7, distance 0.362', 'Vip Mybpc1_3': 'tau 13.1, distance 0.471', 'Vip Chat_1': 'tau 8.4, distance 0.480'}
(b) mean of the Sncg t-type means: 9 Vip t-types; Spearman rho(median tau, distance to Sncg) = -0.367, one-sided p (negative) = 0.166
     {'Vip Sncg': 'tau 19.0, distance 0.148', 'Vip Serpinf1_1': 'tau 8.8, distance 0.331', 'Vip Htr1f': 'tau 10.0, distance 0.394', 'Vip Gpc3': 'tau 5.6, distance 0.462', 'Vip C1ql1': 'tau 5.5, distance 0.426', 'Vip Mybpc1_1': 'tau 11.8, distance 0.466', 'Vip Mybpc1_2': 'tau 14.7, distance 0.391', 'Vip Mybpc1_3': 'tau 13.1, distance 0.490', 'Vip Chat_1': 'tau 8.4, distance 0.503'}

# Layer as the claims file's opening comment defines it (normalized depth from the files at the Nissl boundaries,
# targeted layer without a depth; scala.cell_table `layer`); variants: the files' inferred_layer, and the same rule
# on the paper's µm depths.
um_nd = (pd.to_numeric(tab["Soma depth (µm)"], errors="coerce")
         / pd.to_numeric(tab["Cortical thickness (µm)"], errors="coerce").replace(0, np.nan)).reindex(t.index)
layer_defs = [("depth from the files", t.layer), ("files' inferred_layer", t.inferred_layer),
              ("depth from the paper's µm table", pd.Series(scala.layer_from_depth(um_nd.values, t.targeted_layer.values), index=t.index))]
for lname, lay in layer_defs:
    sel = (lay == "5") & (t.family == "Sst") & ~t.excluded & t.index.isin(S.index) & (t.ttype != "Sst Chodl")
    cnt = t[sel].ttype.value_counts()
    types = [n for n in names if cnt.get(n, 0) >= 5 and n in rna["pvsst"]]
    P = np.array([S.loc[t[sel & (t.ttype == n)].index].mean().values for n in types])
    R = np.array([rna["pvsst"][n] for n in types])
    Dt = pdist(R, "correlation")
    r_full = ss.pearsonr(Dt, pdist(P))
    rbd = np.array([t.loc[t[sel & (t.ttype == n)].index, "Rebound (mV)"].mean() for n in types])
    r_reb = ss.pearsonr(Dt, pdist(rbd[:, None]))
    print(f"layer {lname}: L5 Sst cells {int(sel.sum())} from {t[sel].mouse.nunique()} mice; t-types with at least 5: {len(types)}; "
          f"pairs: {len(Dt)}; all 17 features r = {r_full.statistic:.3f}, p = {r_full.pvalue:.3g}; "
          f"rebound alone r = {r_reb.statistic:.3f}, p = {r_reb.pvalue:.3g}")
    if lname == layer_defs[0][0]:
        print(f"    t-type mean rebound (mV): {dict(zip(types, np.round(rbd, 1)))}")

layer depth from the files: L5 Sst cells 143 from 63 mice; t-types with at least 5: 11; pairs: 55; all 17 features r = 0.339, p = 0.0113; rebound alone r = 0.275, p = 0.0422
    t-type mean rebound (mV): {'Sst Myh8_1': np.float64(13.3), 'Sst Myh8_3': np.float64(11.3), 'Sst Htr1a': np.float64(6.1), 'Sst Etv1': np.float64(11.9), 'Sst Pvalb Etv1': np.float64(19.0), 'Sst Crhr2_1': np.float64(7.3), 'Sst Crhr2_2': np.float64(8.6), 'Sst Hpse': np.float64(3.6), 'Sst Calb2': np.float64(5.2), 'Sst Pvalb Calb2': np.float64(3.3), 'Sst Th_1': np.float64(5.1)}
layer files' inferred_layer: L5 Sst cells 143 from 64 mice; t-types with at least 5: 11; pairs: 55; all 17 features r = 0.347, p = 0.00944; rebound alone r = 0.275, p = 0.042

layer depth from the paper's µm table: L5 Sst cells 143 from 64 mice; t-types with at least 5: 11; pairs: 55; all 17 features r = 0.347, p = 0.00944; rebound alone r = 0.275, p = 0.042

from scipy.spatial import cKDTree
sel = (t.ttype == "Vip Mybpc1_2") & ~t.excluded & t.index.isin(S.index)
myb = list(t[sel].index)
pos, Zr = position(myb, "viplamp")
vs = refn == "Vip Sncg"
dist = cKDTree(Zr[vs]).query(pos)[0]
rb = t.loc[myb, "Rebound (mV)"].values
ra = ss.spearmanr(rb, dist)
print(f"Vip Mybpc1_2 cells: {len(myb)} from {t.loc[myb, 'mouse'].nunique()} mice; rebound {rb.min():.1f} to {rb.max():.1f} mV; "
      f"distance to the nearest Vip Sncg reference cell {dist.min():.1f} to {dist.max():.1f}")
print(f"(a) Spearman rho(rebound, distance) = {ra.statistic:.3f}, one-sided p (positive) = {ss.spearmanr(rb, dist, alternative='greater').pvalue:.3f}")
gg = sorted(set(ref["viplamp"]["genes"]) & set(ex.index))
gi = [list(ref["viplamp"]["genes"]).index(g) for g in gg]
pi = [ex.index.get_loc(g) for g in gg]
T = np.log2(ref["viplamp"]["counts"][:, gi].toarray() + 1)
target = T[vs].mean(axis=0)
E = ex.loc[gg, myb].values.T / (gl.exon_bp.values[pi] / 1000) + it.loc[gg, myb].values.T / ((gl.intron_bp.values[pi] + .001) / 1000)
cc = np.array([np.corrcoef(np.log2(e + 1), target)[0, 1] for e in E])
rb_b = ss.spearmanr(rb, cc)
print(f"(b) genes {len(gg)}; correlation with the Vip Sncg mean of log2(UMI + 1) {cc.min():.3f} to {cc.max():.3f}; Spearman rho(rebound, correlation) = "
      f"{rb_b.statistic:.3f}, one-sided p (negative) = {ss.spearmanr(rb, cc, alternative='less').pvalue:.3f}")
target2 = np.log2(np.asarray(ref["viplamp"]["counts"][:, gi][vs].mean(axis=0)).ravel() + 1)
cc2 = np.array([np.corrcoef(np.log2(e + 1), target2)[0, 1] for e in E])
print(f"(b) variant, against log2 of the mean UMI count of Vip Sncg: Spearman rho = {ss.spearmanr(rb, cc2).statistic:.3f}, "
      f"one-sided p (negative) = {ss.spearmanr(rb, cc2, alternative='less').pvalue:.3f}")
# Variants with library-size normalization of the Patch-seq counts before log2(x + 1), over all genes:
# (i) the length-normalized counts above scaled to a million per cell; (ii) raw exon plus intron counts per million.
E_all = ex[myb].values.T / (gl.exon_bp.values / 1000) + it[myb].values.T / ((gl.intron_bp.values + .001) / 1000)
C_all = (ex[myb].values + it[myb].values).T.astype(float)
for vname, A in [("length-normalized, per million", E_all), ("counts per million, no length normalization", C_all)]:
    A = A / A.sum(axis=1, keepdims=True) * 1e6
    ccv = np.array([np.corrcoef(np.log2(a + 1), target)[0, 1] for a in A[:, pi]])
    print(f"(b) variant, {vname}: Spearman rho = {ss.spearmanr(rb, ccv).statistic:.3f}, "
          f"one-sided p (negative) = {ss.spearmanr(rb, ccv, alternative='less').pvalue:.3f}")
vsn = t[(t.ttype == "Vip Sncg") & ~t.excluded & t.index.isin(S.index)]
print(f"median rebound: Vip Sncg {vsn['Rebound (mV)'].median():.2f} mV (n = {len(vsn)}), Vip Mybpc1_2 {np.median(rb):.2f} mV")
print("variability: s.d. of input resistance, tau and rebound in Vip Mybpc1_2:",
      t.loc[myb, ["Input resistance (MOhm)", "Membrane time constant (ms)", "Rebound (mV)"]].std().round(2).to_dict())

Using a common set of 490 genes.

Vip Mybpc1_2 cells: 34 from 22 mice; rebound 0.3 to 19.8 mV; distance to the nearest Vip Sncg reference cell 0.6 to 30.8
(a) Spearman rho(rebound, distance) = 0.356, one-sided p (positive) = 0.019
(b) genes 490; correlation with the Vip Sncg mean of log2(UMI + 1) 0.245 to 0.605; Spearman rho(rebound, correlation) = 0.067, one-sided p (negative) = 0.648
(b) variant, against log2 of the mean UMI count of Vip Sncg: Spearman rho = 0.066, one-sided p (negative) = 0.644
(b) variant, length-normalized, per million: Spearman rho = 0.068, one-sided p (negative) = 0.648
(b) variant, counts per million, no length normalization: Spearman rho = 0.016, one-sided p (negative) = 0.536
median rebound: Vip Sncg 2.26 mV (n = 23), Vip Mybpc1_2 11.31 mV
variability: s.d. of input resistance, tau and rebound in Vip Mybpc1_2: {'Input resistance (MOhm)': 168.57, 'Membrane time constant (ms)': 6.49, 'Rebound (mV)': 5.56}

sel = (t.ttype == "Sst Pvalb Calb2") & ~t.excluded & t.index.isin(S.index)
spc = list(t[sel].index)
pos, Zr = position(spc, "pvsst")
rn = names[ref["pvsst"]["clusters"]]
c0 = Zr[np.array([n.startswith("Sst") for n in rn])].mean(axis=0)
c1 = Zr[np.array([n.startswith("Pvalb") for n in rn])].mean(axis=0)
u = (c1 - c0) / np.linalg.norm(c1 - c0)
proj = (pos - c0) @ u
ap = t.loc[spc, "Max number of APs"].values
ra = ss.spearmanr(ap, proj)
print(f"Sst Pvalb Calb2 cells: {len(spc)} from {t.loc[spc, 'mouse'].nunique()} mice; max APs {ap.min():.0f} to {ap.max():.0f}; "
      f"projection {proj.min():.1f} to {proj.max():.1f} (Sst centroid 0, Pvalb centroid {np.linalg.norm(c1 - c0):.1f})")
print(f"(a) Spearman rho(max APs, position towards Pvalb) = {ra.statistic:.3f}, one-sided p (positive) = "
      f"{ss.spearmanr(ap, proj, alternative='greater').pvalue:.3f}")
hi = ap > np.median(ap)
k_hi = int(hi.sum())
print(f"max APs, sorted: {sorted(ap.astype(int).tolist())}; median {np.median(ap):.1f}; cells equal to the median: "
      f"{int((ap == np.median(ap)).sum())}; above: {k_hi}; below: {int((ap < np.median(ap)).sum())}")
print("cells above the median:", sorted(np.array(spc)[hi].tolist()))
obs = pdist(pos[hi]).mean()
# Reading (b) as the criterion fixes it (answer of 2026-10-21): exact over all subsets of the same size when
# there are at most 1,000,000, otherwise 100,000 random subsets with p = (k + 1) / (N + 1), and (b) holds
# only when p < 0.05 and the upper bound of the 95% Clopper-Pearson interval is also below 0.05.
n_sub = math.comb(len(ap), k_hi)
print(f"(b) cells above the median: {k_hi} of {len(ap)}; mean pairwise distance {obs:.3f}; subsets of that size: {n_sub:,}")
if n_sub <= 1_000_000:
    null = np.array([pdist(pos[list(s)]).mean() for s in itertools.combinations(range(len(ap)), k_hi)])
    k = int(np.sum(null <= obs))
    p_b, ci_b = k / n_sub, (k / n_sub, k / n_sub)
    print(f"    exact over all {n_sub:,} subsets: p = {p_b:.4f}")
else:
    N = 100_000
    rng = np.random.default_rng(55)
    null = np.array([pdist(pos[rng.choice(len(ap), k_hi, replace=False)]).mean() for _ in range(N)])
    k = int(np.sum(null <= obs))
    p_b = (k + 1) / (N + 1)
    ci_b = ss.binomtest(k, N).proportion_ci(confidence_level=0.95, method="exact")
    ci_b = (ci_b.low, ci_b.high)
    print(f"    {N:,} random subsets (seed 55): {k} at or below the observed distance; p = (k + 1) / (N + 1) = {p_b:.4f}; "
          f"95% Clopper-Pearson interval {ci_b[0]:.4f} to {ci_b[1]:.4f}")
print(f"    null mean pairwise distance {null.mean():.3f} (5th percentile {np.percentile(null, 5):.3f}); "
      f"observed minus null mean {obs - null.mean():.3f}")
holds_a = ra.statistic > 0
holds_b = (p_b < 0.05) and (ci_b[1] < 0.05)
print(f"(a) holds: {holds_a}; (b) holds (p < 0.05 and upper bound < 0.05): {holds_b}")
# Variants, not the criterion's: the earlier 1,000-permutation estimate over six seeds, to show its spread
for seed in [55, 1, 2, 3, 4, 5]:
    r_ = np.random.default_rng(seed)
    nl = np.array([pdist(pos[r_.choice(len(ap), k_hi, replace=False)]).mean() for _ in range(1000)])
    print(f"   variant, 1,000 permutations with seed {seed}: p = {(np.sum(nl <= obs) + 1) / 1001:.3f}")
# Variants, not the criterion's: the same 100,000-subset estimate with other seeds, and the exact p-value
# over all subsets, to show how far the criterion's estimate is from its own threshold
D = np.linalg.norm(pos[:, None, :] - pos[None, :, :], axis=-1)
def mean_dist(idx):
    return D[idx[:, :, None], idx[:, None, :]].sum(axis=(1, 2)) / (k_hi * (k_hi - 1))
for seed in [56, 57, 58, 59, 60]:
    r_ = np.random.default_rng(seed)
    idx = np.argsort(r_.random((100_000, len(ap))), axis=1)[:, :k_hi]
    kk = int(np.sum(mean_dist(idx) <= obs + 1e-12))
    ci = ss.binomtest(kk, 100_000).proportion_ci(confidence_level=0.95, method="exact")
    print(f"   variant, 100,000 subsets with seed {seed}: p = {(kk + 1) / 100_001:.4f}, interval {ci.low:.4f} to {ci.high:.4f}, "
          f"(b) holds {((kk + 1) / 100_001 < 0.05) and (ci.high < 0.05)}")
k_all, n_all = 0, 0
combos = itertools.combinations(range(len(ap)), k_hi)
while True:
    chunk = np.fromiter(itertools.islice(combos, 500_000), dtype=np.dtype((np.int64, k_hi)))
    if len(chunk) == 0:
        break
    k_all += int(np.sum(mean_dist(chunk) <= obs + 1e-12))
    n_all += len(chunk)
print(f"   variant, exact over all {n_all:,} subsets: {k_all:,} at or below the observed distance, p = {k_all / n_all:.4f}")

Using a common set of 493 genes.

Sst Pvalb Calb2 cells: 26 from 19 mice; max APs 15 to 87; projection 15.8 to 37.0 (Sst centroid 0, Pvalb centroid 57.4)
(a) Spearman rho(max APs, position towards Pvalb) = 0.363, one-sided p (positive) = 0.034
max APs, sorted: [15, 19, 20, 22, 23, 24, 25, 25, 26, 27, 27, 30, 31, 34, 37, 38, 41, 42, 43, 45, 46, 47, 53, 57, 76, 87]; median 32.5; cells equal to the median: 0; above: 13; below: 13
cells above the median: ['20180820_sample_7', '20181002_sample_9', '20181003_sample_10', '20181003_sample_9', '20181004_sample_5', '20181004_sample_8', '20181011_sample_3', '20190131_sample_5', '20190204_sample_17', '20190205_sample_6', '20190221_sample_4', '20190318_sample_4', '20191010_sample_2']
(b) cells above the median: 13 of 26; mean pairwise distance 7.054; subsets of that size: 10,400,600

    100,000 random subsets (seed 55): 4852 at or below the observed distance; p = (k + 1) / (N + 1) = 0.0485; 95% Clopper-Pearson interval 0.0472 to 0.0499
    null mean pairwise distance 8.493 (5th percentile 7.067); observed minus null mean -1.438
(a) holds: True; (b) holds (p < 0.05 and upper bound < 0.05): True
   variant, 1,000 permutations with seed 55: p = 0.052
   variant, 1,000 permutations with seed 1: p = 0.045
   variant, 1,000 permutations with seed 2: p = 0.057

   variant, 1,000 permutations with seed 3: p = 0.048
   variant, 1,000 permutations with seed 4: p = 0.055
   variant, 1,000 permutations with seed 5: p = 0.050

   variant, 100,000 subsets with seed 56: p = 0.0500, interval 0.0486 to 0.0513, (b) holds False
   variant, 100,000 subsets with seed 57: p = 0.0498, interval 0.0484 to 0.0511, (b) holds False

   variant, 100,000 subsets with seed 58: p = 0.0509, interval 0.0495 to 0.0522, (b) holds False
   variant, 100,000 subsets with seed 59: p = 0.0496, interval 0.0483 to 0.0510, (b) holds False

   variant, 100,000 subsets with seed 60: p = 0.0514, interval 0.0500 to 0.0528, (b) holds False

   variant, exact over all 10,400,600 subsets: 516,815 at or below the observed distance, p = 0.0497

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

run time: 17.6 min

Claim	Tested here
fig2-rrr	yes
ed5-rrr-ion-channels	yes
fig4c-distance-correlation	yes
ed7-within-families	yes
fig4d-it-depth	yes
fig4a-vip-tau	yes
fig4b-sst-rebound	yes
fig5c-vip-mybpc1-rebound	yes
fig5d-sst-pvalb-calb2-firing	yes

Claim	Outcome
fig2-rrr	reproduced
ed5-rrr-ion-channels	reproduced
fig4c-distance-correlation	reproduced
ed7-within-families	reproduced
fig4d-it-depth	reproduced
fig4a-vip-tau	blocked (methods underspecified)
fig4b-sst-rebound	reproduced
fig5c-vip-mybpc1-rebound	blocked (methods underspecified)
fig5d-sst-pvalb-calb2-firing	reproduced

Scala, Kobak et al. (2021): dandiset 000008 with the paper's gene counts and reference data¶

Loading¶

Reduced-rank regression¶

fig2-rrr¶

ed5-rrr-ion-channels¶

Within-family structure¶

fig4c-distance-correlation¶

ed7-within-families¶

fig4d-it-depth¶

fig4a-vip-tau¶

fig4b-sst-rebound¶

fig5c-vip-mybpc1-rebound¶

fig5d-sst-pvalb-calb2-firing¶

Summary¶

	lambda	genes (full data)	genes (mean over folds)	R2	corr 1	corr 2	corr 3	corr 4	corr 5
0	0.50	24	23.2	0.370	0.895	0.725	0.648	0.384	0.298
1	0.45	29	27.8	0.382	0.898	0.744	0.669	0.405	0.316

	lambda	genes (full data)	genes (mean over folds)	R2	corr 1	corr 2	corr 3	corr 4	corr 5
0	0.50	26	23.7	0.367	0.893	0.731	0.649	0.385	0.297
1	0.45	28	27.1	0.376	0.896	0.744	0.665	0.402	0.313

	lambda	genes (full data)	genes (mean over folds)	R2	corr 1	corr 2	corr 3	corr 4	corr 5
0	0.50	24	23.5	0.371	0.895	0.725	0.655	0.392	0.272
1	0.45	29	27.8	0.381	0.899	0.739	0.668	0.417	0.291

	t-types	pairs	r
family
Pvalb	8	28	0.747
Sst	14	91	0.494
Vip	9	36	0.403
Lamp5	5	10	0.781
IT	8	28	0.328
CT	4	6	0.967

	Pvalb	Sst	Vip	Lamp5	IT	CT
Resting membrane potential (mV)	0.34	0.14	0.05	0.75	0.29	0.35
Input resistance (MOhm)	0.67	0.13	0.61	0.57	0.05	0.90
Membrane time constant (ms)	0.61	0.15	0.40	0.94	-0.01	0.25
AP threshold (mV)	0.08	0.08	0.02	0.21	0.16	0.84
AP amplitude (mV)	0.78	0.19	-0.00	-0.17	-0.13	-0.32
AP width (ms)	0.58	-0.08	0.15	0.47	0.12	0.29
Upstroke-to-downstroke ratio	0.56	0.03	0.05	0.89	0.03	0.83
Afterhyperpolarization (mV)	0.35	0.27	0.25	0.77	0.09	0.44
ISI adaptation index	0.14	0.15	0.10	0.24	0.08	0.28
Max number of APs	0.09	0.12	0.10	0.16	0.01	-0.20
Rheobase (pA)	0.68	0.14	-0.13	0.20	0.17	0.92
Sag ratio	0.63	-0.03	0.19	0.46	0.12	0.92
Latency (ms)	0.08	-0.10	-0.02	0.79	0.03	0.64
ISI coefficient of variation	-0.10	0.41	0.30	0.18	0.25	-0.46
Rebound (mV)	0.64	0.28	0.03	0.66	0.20	0.87
AP amplitude adaptation index	-0.12	0.08	0.15	0.62	-0.03	-0.69
AP coefficient of variation	0.27	0.46	0.06	-0.06	0.22	-0.24