Reproducing Azcorra, Gaertner et al. (2023) from dandiset 001038¶
Paper. Azcorra M, Gaertner Z, Davidson C, et al., Awatramani R, Dombeck DA (2023). Unique functional responses differentially map onto genetic subtypes of dopamine neurons. Nature Neuroscience 26, 1762-1774. doi:10.1038/s41593-023-01401-9. This attempt used the full text, which is open access (CC BY 4.0), from PubMed Central (PMC10545540).
Dandiset. DANDI:001038, version 0.241111.1719. The paper cites no DANDI version (its data statement points to Zenodo records that the dandiset was converted from). This is the only published version; the draft holds the same 1238 files with the same blobs and differs only in metadata (it fixes a malformed DOI link to the paper and adds affiliations), none of which the analysis uses.
What the dandiset holds. Fiber photometry from 102 mice of five genotypes (Vglut2, Calb1, Anxa1, Aldh1a1 and DAT lines), one file per recording depth, with normalized dF/F, its 405 nm control, treadmill velocity and acceleration, and for sessions with stimuli the reward, air puff, light and licking events. These are the data of Figures 2 to 7 and Extended Data Figures 1 and 6 to 10. The snRNA-seq of Figure 1 is in GEO (GSE222558) and the histology images on Zenodo; two snRNA-seq claims are tested here on the GEO matrices, counted apart from the dandiset.
Code read. The paper's analysis code, https://github.com/DombeckLab/Azcorra2023 at commit
fad12b25313f21611dd70c0056fb8bee48beeda1 (MATLAB; preprocessing, exclusion criteria and one script
per figure group), and the conversion code, https://github.com/catalystneuro/dombeck-lab-to-nwb at
commit 5794df1e13b1a44ab2f21ad7e5b38ef8b85940f7. The helper azcorra.py ports the code's
selection and measures; its docstring says where each parameter comes from. Where the Methods
and the code disagree (response windows, timing), the Methods give the main result and the code
a variant.
| Claim | Data | Status |
|---|---|---|
| recording-counts-locomotion | dandiset | tested |
| recording-counts-stimuli | dandiset | tested |
| fig2f-xcorr-direction | dandiset | tested |
| fig2i-trough-lag | dandiset | tested |
| fig2jk-pca-variance | dandiset | tested |
| fig2l-pc-angles | dandiset | tested |
| fig2g-event-counts | dandiset | tested |
| ed7f-prob-transient-after-event | dandiset | tested |
| ed6f-dec-peak-timing | dandiset | tested |
| fig5ab-transient-scaling | dandiset | tested |
| fig4g-reward-response | dandiset | tested |
| fig4g-anxa1-no-reward | dandiset | tested |
| fig4g-airpuff-response | dandiset | tested |
| fig4h-reward-vs-airpuff | dandiset | tested |
| fig4i-reward-size | dandiset | tested |
| ed8f-rest-reward-same | dandiset | tested |
| fig6c-kmeans | dandiset | tested |
| fig7f-soma-axon-correlation | dandiset | tested |
| ed9c-snc-stimulus-responses | dandiset | tested |
| ed9k-snc-variance-explained | dandiset | tested |
| ed1h-aldh1a1-depth | dandiset | tested |
| ed1i-aldh1a1-reward-locomotion | dandiset | tested |
| fig3a-overlap-locomotion | dandiset | tested |
| fig4l-overlap-stimuli | dandiset | tested |
| snrna-nuclei-qc | elsewhere (GEO GSE222558) | tested |
| snrna-anxa1-within-aldh1a1 | elsewhere (GEO GSE222558) | tested |
| snrna-clusters | elsewhere (GEO GSE222558) | blocked: not_runnable (Seurat pipeline in R) |
| scrna-meta-integration | elsewhere (four public datasets) | blocked: not_runnable (Seurat and LIGER in R) |
| histology-anxa1-cre | elsewhere (Zenodo images) | blocked: other (manual cell counting on images) |
The notebook reads all 1238 files in parallel processes and loads the series of the 653 files with a fiber above the signal-to-noise threshold; the run time is printed at the end.
import os
for v in ("OMP_NUM_THREADS", "OPENBLAS_NUM_THREADS", "MKL_NUM_THREADS"):
os.environ.setdefault(v, "1")
import time
import warnings
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy import stats
from sklearn.cluster import KMeans
import azcorra as az
warnings.filterwarnings("ignore", category=RuntimeWarning)
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
ORDER = ["Vglut2", "Calb1", "Anxa1", "DAT", "Aldh1a1"]
COLOR = dict(zip(ORDER, ["#2a78d6", "#eb6834", "#1baf7a", "#7a7a72", "#e87ba4"]))
def p_sr(x, factor=1):
"""Two-sided Wilcoxon signed-rank p-value times a Bonferroni factor, capped at 1."""
x = np.asarray(x, float)
x = x[~np.isnan(x)]
if len(x) < 2:
return np.nan
return min(1.0, stats.wilcoxon(x).pvalue * factor)
def p_rs(a, b, factor=1):
"""Two-sided rank-sum (Mann-Whitney U) p-value times a factor, capped at 1."""
return min(1.0, stats.mannwhitneyu(a, b, alternative="two-sided").pvalue * factor)
def mice(d):
return d.mouse.nunique()
started = time.time()
assets = az.list_assets()
recs = az.load_all(assets, workers=16)
print(f"read {len(recs)} files of dandiset {az.DANDISET_ID} version {az.DANDISET_VERSION} "
f"in {time.time() - started:.0f} s; series loaded for {sum(r['loaded'] for r in recs)} files "
f"(those with a fiber whose signal-to-noise ratio exceeds 10)")
read 1238 files of dandiset 001038 version 0.241111.1719 in 86 s; series loaded for 653 files (those with a fiber whose signal-to-noise ratio exceeds 10)
Data check¶
Each check of the metadata against the values is stated under its output.
df = az.recordings_table(recs)
files = pd.DataFrame([{k: r[k] for k in ("path", "subject", "sex", "genotype", "age", "prefix", "session_id", "has_events")}
for r in recs])
files["subtype"] = files.prefix.map(az.SUBTYPE)
files["mouse"] = files.subject.str.rstrip("L")
print("Files:", len(files), "; subject ids:", files.subject.nunique(), "; mice (trailing L removed):", files.mouse.nunique())
print("Archive summary: 105 subjects, species Mus musculus\n")
print(files.groupby("subtype").agg(files=("path", "size"), mice=("mouse", "nunique"),
with_stimuli=("has_events", "sum")).loc[ORDER].to_string())
print("\nGenotype per subtype:")
print(files.groupby("subtype").genotype.unique().to_string())
print("\nSex of files:", files.sex.value_counts().to_dict(), "; age values:", files.age.unique().tolist())
print("Subject ids ending in L:", sorted(files.subject[files.subject.str.endswith("L")].unique()))
print("Files per mouse: median", files.groupby("mouse").size().median(), ", range",
files.groupby("mouse").size().min(), "to", files.groupby("mouse").size().max())
print("\nFibers (one per column of the dF/F series):", len(df), "; files with two fibers:",
int((df.groupby("path").size() == 2).sum()))
print(pd.crosstab(df.location, df.target).to_string())
Files: 1238 ; subject ids: 105 ; mice (trailing L removed): 102
Archive summary: 105 subjects, species Mus musculus
files mice with_stimuli
subtype
Vglut2 265 24 97
Calb1 271 25 89
Anxa1 205 14 93
DAT 208 18 105
Aldh1a1 289 21 90
Genotype per subtype:
subtype
Aldh1a1 [Aldh1a1-iCre/Th-Flpo]
Anxa1 [Anxa1-iCre]
Calb1 [Calb1-IRES2-Cre/DAT-tTA/Ai93D]
DAT [DAT-IRES-Cre (RRID, IMSR_JAX:027178)]
Vglut2 [Vglut2-IRES-Cre/Th-Flpo (RRID, IMSR_JAX:016963)]
Sex of files: {'F': 612, 'M': 500, 'U': 126} ; age values: ['P8W/P16W']
Subject ids ending in L: ['Anxa-I765L', 'Anxa-I767L', 'Calb-2246L']
Files per mouse: median 12.0 , range 1 to 32
Fibers (one per column of the dF/F series): 2137 ; files with two fibers: 899
target axons in striatum cell bodies in SNc
location
DLS 112 0
DMS 112 0
DS 823 0
SNc 0 790
TS 300 0
The files hold 105 subject ids for 102 mice: three ids end in "L" (Anxa-I765L, Anxa-I767L,
Calb-2246L), which the paper's code merges with the same id without the letter (a mouse recorded
once from the left hemisphere); mice are counted that way here. The age is a range (P8W/P16W) for
every file. Each file is one recording at one fiber depth; a session (date) holds several.
loc = df.groupby("location")[["ml", "ap", "dv"]].agg(["mean", "min", "max"]).round(2)
print("Fiber coordinates by location label, columns in stored order (the file says AP, ML, DV):")
print(loc.to_string())
print("\nInjection coordinates of the indicator device, stored as AP, ML, DV: -3.25, 1.55, -3.8 (Methods: -3.25 mm caudal, +1.55 mm lateral)")
print("Distinct values of the second coordinate:", sorted(df.ap.dropna().round(2).unique().tolist()))
print("Fibers without coordinates:", int(df.ml.isna().sum()), "of", len(df),
"; SNc fibers without coordinates:", int(df[df.snc].ml.isna().sum()), "of", int(df.snc.sum()))
Fiber coordinates by location label, columns in stored order (the file says AP, ML, DV):
ml ap dv
mean min max mean min max mean min max
location
DLS 2.73 2.04 3.07 0.49 0.14 0.86 3.59 2.93 4.58
DMS 1.53 0.82 1.97 0.49 0.14 0.86 3.24 2.52 4.27
DS 1.76 1.15 2.98 0.57 -0.94 0.86 3.17 1.67 4.88
SNc 2.52 2.50 2.55 0.50 0.50 0.50 4.98 4.56 5.46
TS 2.73 2.15 3.16 -0.92 -1.34 0.50 3.47 2.16 5.10
Injection coordinates of the indicator device, stored as AP, ML, DV: -3.25, 1.55, -3.8 (Methods: -3.25 mm caudal, +1.55 mm lateral)
Distinct values of the second coordinate: [-1.34, -0.94, -0.58, 0.14, 0.5, 0.86]
Fibers without coordinates: 1281 of 2137 ; SNc fibers without coordinates: 785 of 790
The coordinates column is described as "(AP, ML, DV) mm relative to Bregma". The values do not
fit that order: the second value takes only the seven reference-slice positions of the Methods
(+0.86 to -1.34 mm) and is negative for the tail of the striatum (TS), as an AP coordinate must
be, and the first value is about 1.5 mm for DMS and 2.7 mm for DLS, as a lateral coordinate is.
The paper's code builds this vector as [ML, AP, DV] (FiberLocations.m). The data show the order is
(ML, AP, DV), and it is read that way here. Most SNc fibers have no coordinates.
L = [r for r in recs if r["loaded"]]
print("Files with series loaded:", len(L))
print("Samples per file at 100 Hz (minutes):", pd.Series([r["n_samples"] / 6000 for r in L]).describe().round(1).to_dict())
print("Velocity range over files (m/s):", round(min(r["vel_range"][0] for r in L), 3), "to", round(max(r["vel_range"][1] for r in L), 3))
print("Acceleration range over files (m/s2):", round(min(r["acc_range"][0] for r in L), 2), "to", round(max(r["acc_range"][1] for r in L), 2))
dmax = pd.Series([m for r in L for m in r["dff_max"]])
print("Maximum of each dF/F column:", dmax.min(), "to", dmax.max(), "; fraction of NaN samples, max:",
max(m for r in L for m in r["dff_nan"]))
acc = pd.Series([r["acc_check_maxabs"] for r in L])
print("\nAcceleration recomputed as diff(smooth(velocity, 6)) x 100 (data_processing.m), largest |difference| per file (m/s2):")
print(acc.describe().round(6).to_dict(), "; files above 1e-6:", int((acc > 1e-6).sum()))
Files with series loaded: 653
Samples per file at 100 Hz (minutes): {'count': 653.0, 'mean': 16.0, 'std': 4.1, 'min': 0.3, '25%': 14.2, '50%': 14.9, '75%': 17.8, 'max': 37.1}
Velocity range over files (m/s): -3.335 to 3.277
Acceleration range over files (m/s2): -65.02 to 66.61
Maximum of each dF/F column: 1.0 to 1.0 ; fraction of NaN samples, max: 4.999625028122891e-05
Acceleration recomputed as diff(smooth(velocity, 6)) x 100 (data_processing.m), largest |difference| per file (m/s2):
{'count': 653.0, 'mean': 0.007343, 'std': 0.079154, 'min': 0.0, '25%': 0.0, '50%': 0.0, '75%': 0.0, 'max': 1.748186} ; files above 1e-6: 13
The dF/F series reach exactly 1 in every column: they are the paper's normalized dF/F (0 to 1; the
paper reports the same quantity on a 0 to 100 scale). The files call them "DF/F calculated from
green fluorescence" with unit "n.a." and do not say they are normalized; the normalization factor
is in the column normalized_fluorescence, whose name suggests a normalized value. Velocity is in
m/s and the acceleration is the code's derivative of the smoothed velocity in all but 13 files,
which differ by up to 1.75 m/s2. The sampling rate is 100 Hz in every series.
mo = pd.Series([r["movonoff_agree"] for r in L if "movonoff_agree" in r])
print("Movement bouts (MovOnOff) recomputed from velocity against the stored Events rows, fraction of samples that agree:")
print(mo.describe().round(5).to_dict())
w = pd.DataFrame([r["wheel_check"] for r in L])
print("\nAcceleration/deceleration onsets: stored", int(w.acc_file.sum()), "/", int(w.dec_file.sum()),
"; recomputed", int(w.acc_re.sum()), "/", int(w.dec_re.sum()), "; in both", int(w.acc_both.sum()), "/", int(w.dec_both.sum()))
print("Files where the recomputed onsets include ones not stored:", int(((w.acc_re > w.acc_both) | (w.dec_re > w.dec_both)).sum()))
a405 = df.dropna(subset=["acc405_re"])
print("\n405 nm peak cross-correlation with acceleration, recomputed over running against the stored column:")
print("|difference|", (a405.acc405_re - a405.acc405_file).abs().describe().round(4).to_dict())
g = df[df.passes_snr]
print("Stored value present (the code computes it only with >100 s running) against running time recomputed > 100 s:")
print(pd.crosstab(g.time_run > 100, g.acc405_file.notna(), rownames=["run > 100 s"], colnames=["stored"]).to_string())
Movement bouts (MovOnOff) recomputed from velocity against the stored Events rows, fraction of samples that agree:
{'count': 279.0, 'mean': 0.99995, 'std': 0.00014, 'min': 0.99919, '25%': 0.99998, '50%': 0.99999, '75%': 1.0, 'max': 1.0}
Acceleration/deceleration onsets: stored 62649 / 61595 ; recomputed 63560 / 62362 ; in both 62649 / 61595
Files where the recomputed onsets include ones not stored: 146
405 nm peak cross-correlation with acceleration, recomputed over running against the stored column:
|difference| {'count': 552.0, 'mean': 0.0029, 'std': 0.0062, 'min': 0.0, '25%': 0.0, '50%': 0.0003, '75%': 0.0036, 'max': 0.0611}
Stored value present (the code computes it only with >100 s running) against running time recomputed > 100 s:
stored False True
run > 100 s
False 275 5
True 0 552
The recomputed movement bouts agree with the stored ones on more than 99.9% of samples, and the
recomputed acceleration and deceleration onsets contain every stored one. The extra recomputed
onsets fall in the 5 s after rewards: the stored onsets were also excluded after rewards the mouse
did not consume, which the Events table does not list (it keeps consumed rewards only, tagged long
or short). The same exclusion applies to the running periods, so running time here can be slightly
longer than the code's after unconsumed rewards. The 405 nm cross-correlation recomputed from the
recomputed running periods matches the stored value (median difference near 0), and the stored
value is present exactly where the recomputed running time exceeds 100 s, apart from the few cases
in the table. The column cross_correlation_with_acceleration is described as "the
cross-correlation between dF/F and acceleration"; it holds the peak cross-correlation of the
405 nm control with acceleration, the code's Acc405.
c = pd.DataFrame([r["clock"] for r in L if "clock" in r])
print("Around consumed rewards: licking (fraction of samples) and acceleration (m/s2), 1 s after minus 1 s before,")
print("at the stored times and at the same times shifted by 5 to 15 s; mean over", len(c), "files with rewards")
print(c.agg(["mean", "std"]).round(3).to_string())
print("Files where licking rises after the reward:", int((c.lick_true > 0).sum()), "of", len(c),
"; with shifted times:", int((c.lick_shift > 0).sum()))
Around consumed rewards: licking (fraction of samples) and acceleration (m/s2), 1 s after minus 1 s before,
at the stored times and at the same times shifted by 5 to 15 s; mean over 255 files with rewards
lick_true lick_shift acc_true acc_shift
mean 0.427 -0.006 -0.166 -0.003
std 0.163 0.041 0.159 0.084
Files where licking rises after the reward: 254 of 255 ; with shifted times: 95
The reward times share the clock of the licking and treadmill series: licking rises after
consumed rewards (by about 0.43 of samples) and mice decelerate, while at shifted times both
contrasts are near 0. The dF/F series share the 100 Hz clock of velocity by construction (one
starting_time and rate). That the photometry is aligned to the stimuli is consistent with the
stimulus responses tested below, not shown independently of them.
hashes = pd.Series([h for r in L for h in r["dff_hash"]])
vh = pd.Series([r["vel_hash"] for r in L])
print("dF/F columns with an identical copy elsewhere:", int(hashes.duplicated(keep=False).sum()), "of", len(hashes))
print("Files with an identical velocity trace elsewhere:", int(vh.duplicated(keep=False).sum()), "of", len(vh))
print("Files never loaded (no fiber with signal-to-noise ratio above 10):", len(recs) - len(L))
print("\nSignal-to-noise ratio column: zero in", int((df.snr == 0).sum()), "of", len(df), "fibers; above 10 in", int(df.passes_snr.sum()))
ev = pd.DataFrame([r["n_events"] for r in L if r["has_events"]])
print("Events per file with stimuli (median):", ev.median().to_dict())
nn = df.assign(n_norm=df.norm.map(len)).groupby("path").n_norm.max()
print("\nValues per fiber in normalized_fluorescence:", nn.value_counts().sort_index().to_dict(), "(files by the largest count)")
for v in (az.DANDISET_VERSION, "draft"):
md = az.version_metadata(v)
print(f"\nDandiset metadata, version {v}:")
for r in md.get("relatedResource", []):
print(f" related resource {r.get('relation')}: url {r.get('url')}")
s = md["assetsSummary"]
print(" assetsSummary measurementTechnique:", [x.get("name") for x in s.get("measurementTechnique", [])],
"; variableMeasured:", s.get("variableMeasured"))
dF/F columns with an identical copy elsewhere: 0 of 1116
Files with an identical velocity trace elsewhere: 0 of 653
Files never loaded (no fiber with signal-to-noise ratio above 10): 585
Signal-to-noise ratio column: zero in 994 of 2137 fibers; above 10 in 832
Events per file with stimuli (median): {'reward': 19.0, 'long': 8.0, 'short': 9.0, 'rest': 3.0, 'airpuff': 11.0, 'light': 11.0, 'lick': 585.0}
Values per fiber in normalized_fluorescence: {1: 1115, 2: 123} (files by the largest count)
Dandiset metadata, version 0.241111.1719:
related resource dcite:IsDescribedBy: url https://doi.org/https://doi.org/10.1038/s41593-023-01401-9
related resource dcite:IsSupplementedBy: url https://github.com/DombeckLab/Azcorra2023/
related resource dcite:IsDerivedFrom: url https://zenodo.org/records/7871634
related resource dcite:IsDerivedFrom: url https://zenodo.org/records/7871982
related resource dcite:IsCompiledBy: url https://github.com/catalystneuro/dombeck-lab-to-nwb
assetsSummary measurementTechnique: ['analytical technique'] ; variableMeasured: ['ProcessingModule']
Dandiset metadata, version draft: related resource dcite:IsDescribedBy: url https://doi.org/10.1038/s41593-023-01401-9 related resource dcite:IsSupplementedBy: url https://github.com/DombeckLab/Azcorra2023/ related resource dcite:IsDerivedFrom: url https://zenodo.org/records/7871634 related resource dcite:IsDerivedFrom: url https://zenodo.org/records/7871982 related resource dcite:IsCompiledBy: url https://github.com/catalystneuro/dombeck-lab-to-nwb assetsSummary measurementTechnique: ['analytical technique'] ; variableMeasured: ['ProcessingModule']
No recording is deposited twice: no dF/F column or velocity trace is repeated among the 653 files
loaded (the others are not analysed). The signal-to-noise ratio is 0 for 46% of fibers, the code's
value for traces with too few transients. Rewards appear twice in the Events table, once tagged
long or short and once more, 0.01 s later, tagged rest when delivered at rest; the code's consumed
reward times are the long and short rows plus one bin. The column normalized_fluorescence
holds one value per fiber in 1115 files and two in 123. In the published version the related
resource that describes the dandiset, the paper, has the malformed URL
https://doi.org/https://doi.org/..., which the draft fixes; in both, the asset summary lists
only "analytical technique" as measurement technique and "ProcessingModule" as variable measured.
recording-counts-locomotion¶
The striatal recordings that pass the locomotion criteria number 42 from 12 Vglut2 mice, 22 from 6 Calb1 mice, 47 from 10 Anxa1 mice, 74 from 14 DAT mice and 75 from 14 Aldh1a1 mice.
paper_rec = {"Vglut2": 42, "Calb1": 22, "Anxa1": 47, "DAT": 74, "Aldh1a1": 75}
paper_mice_legend = {"Vglut2": 12, "Calb1": 6, "Anxa1": 10, "DAT": 14, "Aldh1a1": 14}
paper_mice_methods = {"Vglut2": 12, "Calb1": 6, "Anxa1": 9, "DAT": 14, "Aldh1a1": 14}
strf = df[df.striatum & df.passes_snr]
steps = pd.DataFrame({
"snr>10": strf.groupby("subtype").size(),
"+run>100s": strf[strf.time_run > 100].groupby("subtype").size(),
"+405<0.1": strf[(strf.time_run > 100) & (strf.acc405_file < 0.1)].groupby("subtype").size(),
"+artifacts<=5%": strf[strf.loco].groupby("subtype").size()}).loc[ORDER]
print("Striatal fibers left after each criterion:")
print(steps.to_string())
loco = df[df.striatum & df.loco].copy()
t = loco.groupby("subtype").agg(recordings=("path", "size"), mice=("mouse", "nunique")).loc[ORDER]
t["paper_recordings"] = pd.Series(paper_rec)
t["within_10pct"] = (t.recordings - t.paper_recordings).abs() <= 0.1 * t.paper_recordings
t["paper_mice_legend"] = pd.Series(paper_mice_legend)
t["paper_mice_methods"] = pd.Series(paper_mice_methods)
print("\n", t.to_string())
print("\nRecording counts within 10%:", int(t.within_10pct.sum()), "of 5")
print("\nVariant thresholds on the 405 nm cross-correlation (recordings):")
for thr in (0.05, 0.1):
print(f" < {thr}:", strf[(strf.time_run > 100) & (strf.acc405_file < thr) & (strf.bad_all == False)].groupby("subtype").size().loc[ORDER].to_dict())
print("Per mouse (Anxa1):", loco[loco.subtype == "Anxa1"].groupby("mouse").size().to_dict())
# breakdown of the selection: by location, fibers per file and stimuli, per mouse, and variants
strf = strf.assign(n_fibers=strf.path.map(df.groupby("path").size()), date=strf.session_id.str[:8])
stepmask = {"snr>10": strf.index == strf.index, "+run>100s": strf.time_run > 100,
"+405<0.1": (strf.time_run > 100) & (strf.acc405_file < 0.1), "+artifacts<=5%": strf.loco}
for col in ("location", "n_fibers", "has_events"):
print(f"\nStriatal fibers left after each criterion, by {col}:")
tab = pd.concat({k: strf[m].groupby([col, "subtype"]).size() for k, m in stepmask.items()}, axis=1).fillna(0).astype(int)
tab = tab.swaplevel().sort_index()
print(tab.loc[[s for s in ORDER if s in tab.index.get_level_values(0)]].to_string())
print("\nSelected recordings per mouse:")
for s in ORDER:
print(f" {s}:", loco[loco.subtype == s].groupby("mouse").size().to_dict())
sel = strf[strf.loco]
print("\nVariant counts (recordings):")
print(" one per mouse, date and depth:", sel.drop_duplicates(["mouse", "date", "depth", "location"]).groupby("subtype").size().reindex(ORDER).to_dict())
print(" one per mouse and depth: ", sel.drop_duplicates(["mouse", "depth", "location"]).groupby("subtype").size().reindex(ORDER).to_dict())
print(" recorded before submission (2022-09-23):", sel[sel.date < "20220923"].groupby("subtype").size().reindex(ORDER).to_dict())
st_ = df[df.striatum & df.stim & (df.n_AirPuff >= az.N_EVENTS_MIN)]
print(" air puff recordings (which match the paper) made after submission:",
st_[st_.session_id.str[:8] >= "20220923"].groupby("subtype").size().to_dict())
print(" paper:", paper_rec)
Striatal fibers left after each criterion:
snr>10 +run>100s +405<0.1 +artifacts<=5%
subtype
Vglut2 71 51 50 50
Calb1 54 29 28 28
Anxa1 132 84 82 82
DAT 146 103 93 89
Aldh1a1 144 87 85 85
recordings mice paper_recordings within_10pct paper_mice_legend paper_mice_methods
subtype
Vglut2 50 14 42 False 12 12
Calb1 28 7 22 False 6 6
Anxa1 82 9 47 False 10 9
DAT 89 14 74 False 14 14
Aldh1a1 85 15 75 False 14 14
Recording counts within 10%: 0 of 5
Variant thresholds on the 405 nm cross-correlation (recordings):
< 0.05: {'Vglut2': 39, 'Calb1': 24, 'Anxa1': 73, 'DAT': 68, 'Aldh1a1': 73}
< 0.1: {'Vglut2': 50, 'Calb1': 28, 'Anxa1': 82, 'DAT': 89, 'Aldh1a1': 85}
Per mouse (Anxa1): {'Anxa-2559': 7, 'Anxa-A003': 1, 'Anxa-I765': 6, 'Anxa-I767': 19, 'Anxa-J074': 14, 'Anxa-J075': 7, 'Anxa-J077': 9, 'Anxa-J078': 11, 'Anxa-J079': 8}
Striatal fibers left after each criterion, by location:
snr>10 +run>100s +405<0.1 +artifacts<=5%
subtype location
Vglut2 DS 6 4 4 4
TS 65 47 46 46
Calb1 DLS 5 0 0 0
DMS 10 0 0 0
DS 39 29 28 28
Anxa1 DLS 6 2 2 2
DMS 8 5 5 5
DS 118 77 75 75
DAT DS 121 82 75 74
TS 25 21 18 15
Aldh1a1 DLS 12 1 1 1
DMS 16 6 6 6
DS 116 80 78 78
Striatal fibers left after each criterion, by n_fibers:
snr>10 +run>100s +405<0.1 +artifacts<=5%
subtype n_fibers
Vglut2 1 26 23 22 22
2 45 28 28 28
Calb1 1 28 18 17 17
2 26 11 11 11
Anxa1 2 132 84 82 82
DAT 1 30 18 13 12
2 116 85 80 77
Aldh1a1 1 61 44 42 42
2 83 43 43 43
Striatal fibers left after each criterion, by has_events:
snr>10 +run>100s +405<0.1 +artifacts<=5%
subtype has_events
Vglut2 False 42 31 30 30
True 29 20 20 20
Calb1 False 36 14 13 13
True 18 15 15 15
Anxa1 False 71 45 43 43
True 61 39 39 39
DAT False 71 50 41 37
True 75 53 52 52
Aldh1a1 False 114 62 61 61
True 30 25 24 24
Selected recordings per mouse:
Vglut2: {'VGlut-2070': 2, 'VGlut-2292': 3, 'VGlut-2293': 4, 'VGlut-2295': 2, 'VGlut-2298': 1, 'VGlut-2299': 1, 'VGlut-2443': 4, 'VGlut-A997': 3, 'VGlut-B001': 6, 'VGlut-B411': 5, 'VGlut-B413': 6, 'VGlut-B612': 3, 'VGlut-B615': 8, 'VGlut-B617': 2}
Calb1: {'Calb-1789': 7, 'Calb-1790': 3, 'Calb-1791': 3, 'Calb-1792': 8, 'Calb-1877': 2, 'Calb-2290': 4, 'Calb-S448': 1}
Anxa1: {'Anxa-2559': 7, 'Anxa-A003': 1, 'Anxa-I765': 6, 'Anxa-I767': 19, 'Anxa-J074': 14, 'Anxa-J075': 7, 'Anxa-J077': 9, 'Anxa-J078': 11, 'Anxa-J079': 8}
DAT: {'Dat-2111': 10, 'Dat-2114': 8, 'Dat-2115': 1, 'Dat-D001': 8, 'Dat-D002': 4, 'Dat-D087': 8, 'Dat-D088': 4, 'Dat-D092': 8, 'Dat-D093': 4, 'Dat-D095': 7, 'Dat-D170': 3, 'Dat-D171': 5, 'Dat-F554': 5, 'Dat-G557': 14}
Aldh1a1: {'Raldh-1572': 14, 'Raldh-1573': 10, 'Raldh-1697': 7, 'Raldh-1698': 5, 'Raldh-1846': 7, 'Raldh-1847': 2, 'Raldh-1848': 3, 'Raldh-1849': 5, 'Raldh-2076': 7, 'Raldh-2220': 6, 'Raldh-2221': 1, 'Raldh-2276': 12, 'Raldh-2281': 2, 'Raldh-2306': 1, 'Raldh-2403': 3}
Variant counts (recordings):
one per mouse, date and depth: {'Vglut2': 47, 'Calb1': 28, 'Anxa1': 57, 'DAT': 89, 'Aldh1a1': 78}
one per mouse and depth: {'Vglut2': 33, 'Calb1': 21, 'Anxa1': 44, 'DAT': 85, 'Aldh1a1': 59}
recorded before submission (2022-09-23): {'Vglut2': 50, 'Calb1': 27, 'Anxa1': 75, 'DAT': 89, 'Aldh1a1': 85}
air puff recordings (which match the paper) made after submission: {'Anxa1': 14, 'Calb1': 2}
paper: {'Vglut2': 42, 'Calb1': 22, 'Anxa1': 47, 'DAT': 74, 'Aldh1a1': 75}
Result: not reproduced. With every criterion of the Methods and the code applied, more striatal recordings pass than the paper reports, for all five groups: 50 against 42 (Vglut2), 28 against 22 (Calb1), 82 against 47 (Anxa1), 89 against 74 (DAT) and 85 against 75 (Aldh1a1); none is within 10%. The mouse counts are close (14, 7, 9, 14 and 15 against 12, 6, 10 or 9, 14 and 14). A stricter 405 nm threshold of 0.05 does not give the paper's counts either. The stimulus counts below match exactly with the same files, so the dandiset holds the paper's recordings; which of them entered Figure 2 cannot be told from the files. The paper gives no counts by location, fiber setup or session type, so where the excess lies cannot be measured; the tables only show that our selection includes recordings of every kind, and for some groups one kind holds all of them (every Anxa1 recording is from a dual-fiber file, every Calb1 recording is labelled DS). Counting one recording per mouse, date and depth gives 47, 28, 57, 89 and 78, and dropping recordings made after the paper's submission (16 air puff recordings that match the paper's counts were made after it) gives 50, 27, 75, 89 and 85; neither gives the paper's counts. The locomotion claims below use the larger selection.
recording-counts-stimuli¶
The striatal recordings in the reward analysis number 28 (Vglut2), 17 (Calb1), 51 (Anxa1) and 63 (DAT), and in the air puff analysis 29, 17, 57 and 69.
stim = df[df.striatum & df.stim].copy()
stim["has_rew"] = stim.n_Reward >= az.N_EVENTS_MIN
stim["has_puff"] = stim.n_AirPuff >= az.N_EVENTS_MIN
t = pd.DataFrame({
"reward": stim[stim.has_rew].groupby("subtype").size(),
"reward_mice": stim[stim.has_rew].groupby("subtype").mouse.nunique(),
"paper_reward": pd.Series({"Vglut2": 28, "Calb1": 17, "Anxa1": 51, "DAT": 63, "Aldh1a1": 30}),
"paper_reward_mice": pd.Series({"Vglut2": 11, "Calb1": 8, "Anxa1": 8, "DAT": 11, "Aldh1a1": 8}),
"airpuff": stim[stim.has_puff].groupby("subtype").size(),
"airpuff_mice": stim[stim.has_puff].groupby("subtype").mouse.nunique(),
"paper_airpuff": pd.Series({"Vglut2": 29, "Calb1": 17, "Anxa1": 57, "DAT": 69, "Aldh1a1": 30}),
"paper_airpuff_mice": pd.Series({"Vglut2": 12, "Calb1": 8, "Anxa1": 8, "DAT": 11, "Aldh1a1": 8}),
}).loc[ORDER]
print(t.to_string())
main4 = t.loc[["Vglut2", "Calb1", "Anxa1", "DAT"]]
ok = ((main4.reward - main4.paper_reward).abs() <= 0.1 * main4.paper_reward).sum() + \
((main4.airpuff - main4.paper_airpuff).abs() <= 0.1 * main4.paper_airpuff).sum()
print("\nCounts within 10% (of 8):", int(ok))
print("Stimulus recordings without enough consumed rewards (air puffs kept), paper 1/29, 0/17, 6/57:",
{s: f"{int((d.has_puff & ~d.has_rew).sum())}/{int(d.has_puff.sum())}" for s, d in stim.groupby("subtype")})
reward reward_mice paper_reward paper_reward_mice airpuff airpuff_mice paper_airpuff paper_airpuff_mice
Vglut2 28 11 28 11 29 12 29 12
Calb1 17 8 17 8 17 8 17 8
Anxa1 51 8 51 8 57 8 57 8
DAT 63 11 63 11 69 11 69 11
Aldh1a1 30 8 30 8 30 8 30 8
Counts within 10% (of 8): 8
Stimulus recordings without enough consumed rewards (air puffs kept), paper 1/29, 0/17, 6/57: {'Aldh1a1': '0/30', 'Anxa1': '6/57', 'Calb1': '0/17', 'DAT': '10/69', 'Vglut2': '1/29'}
Result: reproduced. All eight counts equal the paper's: reward 28, 17, 51 and 63 recordings and air puff 29, 17, 57 and 69 for Vglut2, Calb1, Anxa1 and DAT, with the paper's mouse counts, and the recordings without consumed rewards are the paper's 1/29, 0/17 and 6/57.
fig2f-xcorr-direction¶
Cross-correlations between dF/F and acceleration show a trough at positive lags for Calb1 and Vglut2 axons (transients after decelerations) and a peak at positive lags for Anxa1 axons.
loco = loco[loco.cc470.map(lambda c: isinstance(c, np.ndarray))].copy()
lags = np.arange(-az.LAG, az.LAG + 1) / az.HZ
loco["ccint"] = loco.cc470.map(lambda c: np.trapezoid(c[az.LAG:], dx=1 / az.HZ))
loco["ccint405"] = loco.cc405.map(lambda c: np.trapezoid(c[az.LAG:], dx=1 / az.HZ))
rows = []
for s in ORDER:
d = loco[loco.subtype == s]
pm = d.groupby("mouse").ccint.mean()
rows.append({"subtype": s, "n": len(d), "mice": mice(d), "mean_int_0_1s": d.ccint.mean(),
"median": d.ccint.median(), "frac_negative": (d.ccint < 0).mean(),
"p_x3": p_sr(d.ccint, 3) if s in ("Vglut2", "Calb1", "Anxa1") else p_sr(d.ccint),
"405_mean": d.ccint405.mean(), "405_p": p_sr(d.ccint405),
"mice_mean_negative": f"{int((pm < 0).sum())}/{len(pm)}"})
print(pd.DataFrame(rows).round(4).to_string(index=False))
print("\nMean cross-correlation at lags 0.1 to 0.5 s:", {s: round(float(np.mean(np.stack(loco[loco.subtype == s].cc470.values)[:, 110:151])), 4) for s in ORDER})
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
for s in ["Vglut2", "Calb1", "Anxa1", "DAT"]:
m = np.stack(loco[loco.subtype == s].cc470.values)
m4 = np.stack(loco[loco.subtype == s].cc405.values)
ax[0].plot(lags, m.mean(0), color=COLOR[s], lw=2, label=f"{s} (n={len(m)})")
ax[1].plot(lags, m4.mean(0), color=COLOR[s], lw=2, label=s)
for a, t_ in zip(ax, ["470 nm dF/F vs acceleration", "405 nm control vs acceleration"]):
a.axhline(0, color="0.7", lw=0.8); a.axvline(0, color="0.7", lw=0.8); a.set_title(t_); a.set_xlabel("lag (s); dF/F after acceleration at positive lags")
ax[0].set_ylabel("mean cross-correlation"); ax[0].legend(frameon=False, fontsize=8)
plt.tight_layout(); plt.show()
allm = np.concatenate([np.stack(loco[loco.subtype == s].cc470.values).mean(0) for s in ORDER])
print("Plotted means range from", round(allm.min(), 3), "to", round(allm.max(), 3), "; all points within the axes")
subtype n mice mean_int_0_1s median frac_negative p_x3 405_mean 405_p mice_mean_negative
Vglut2 50 14 -0.0334 -0.0295 0.9400 0.0000 0.0050 0.0000 14/14
Calb1 28 7 -0.0362 -0.0358 0.9286 0.0000 0.0047 0.0156 6/7
Anxa1 82 9 0.0142 0.0134 0.2561 0.0000 0.0016 0.2232 2/9
DAT 89 14 -0.0133 -0.0151 0.6404 0.0004 0.0054 0.0000 10/14
Aldh1a1 85 15 0.0167 0.0160 0.2824 0.0000 0.0003 0.7342 1/15
Mean cross-correlation at lags 0.1 to 0.5 s: {'Vglut2': -0.0483, 'Calb1': -0.0596, 'Anxa1': 0.0242, 'DAT': -0.012, 'Aldh1a1': 0.0408}
Plotted means range from -0.117 to 0.061 ; all points within the axes
Result: reproduced. The integral of the cross-correlation over lags 0 to 1 s is negative for Vglut2 (mean -0.033, 94% of recordings negative) and Calb1 (-0.036, 93%) and positive for Anxa1 (0.014, 26% negative), each with p < 0.0001 after multiplying by 3. The 405 nm control integrals are about a tenth as large (0.002 to 0.005). The direction holds in 14 of 14 Vglut2 mice and 6 of 7 Calb1 mice, and Anxa1 is positive in 7 of 9 mice. The analysis uses the larger selection (recording-counts-locomotion).
fig2i-trough-lag¶
Calb1 transients follow decelerations with a shorter lag than Vglut2 transients; the trough of the cross-correlation is at 0.42 s for Vglut2 and 0.17 s for Calb1 on average (P = 1e-6).
def trough_methods(c):
return np.argmin(c[az.LAG:]) / az.HZ # lag of the minimum, 0 to 1 s
def trough_code(c):
i = np.argmin(c[az.LAG - 1:]) + 1 # code: sig(neg0:end), 1-based index / hz
return min(i, 100) / az.HZ
loco["lag_m"] = loco.cc470.map(trough_methods)
loco["lag_c"] = loco.cc470.map(trough_code)
for col, name in (("lag_m", "Methods (main)"), ("lag_c", "code (variant)")):
v, c_ = loco[loco.subtype == "Vglut2"][col], loco[loco.subtype == "Calb1"][col]
p = p_rs(v, c_)
print(f"{name}: Vglut2 mean {v.mean():.3f} s (sem {v.sem():.3f}, n={len(v)}), Calb1 {c_.mean():.3f} s (sem {c_.sem():.3f}, n={len(c_)}); "
f"rank-sum p = {p:.2g}, x3 = {min(1, 3 * p):.2g}; difference of means {v.mean() - c_.mean():.3f} s")
print("Anxa1 (reference, its trough is not the measure):", round(loco[loco.subtype == "Anxa1"].lag_m.mean(), 3))
print("\nPer mouse (Methods):")
print(loco[loco.subtype.isin(["Vglut2", "Calb1"])].groupby(["subtype", "mouse"]).lag_m.agg(["mean", "size"]).round(3).to_string())
pmv = loco[loco.subtype == "Vglut2"].groupby("mouse").lag_m.mean(); pmc = loco[loco.subtype == "Calb1"].groupby("mouse").lag_m.mean()
print(f"Mouse means: Vglut2 {pmv.mean():.3f} s ({len(pmv)} mice), Calb1 {pmc.mean():.3f} s ({len(pmc)} mice), rank-sum p = {p_rs(pmv, pmc):.2g}")
Methods (main): Vglut2 mean 0.416 s (sem 0.031, n=50), Calb1 0.140 s (sem 0.023, n=28); rank-sum p = 6.9e-09, x3 = 2.1e-08; difference of means 0.276 s
code (variant): Vglut2 mean 0.436 s (sem 0.031, n=50), Calb1 0.160 s (sem 0.023, n=28); rank-sum p = 6.9e-09, x3 = 2.1e-08; difference of means 0.276 s
Anxa1 (reference, its trough is not the measure): 0.222
Per mouse (Methods):
mean size
subtype mouse
Calb1 Calb-1789 0.119 7
Calb-1790 0.190 3
Calb-1791 0.097 3
Calb-1792 0.121 8
Calb-1877 0.105 2
Calb-2290 0.112 4
Calb-S448 0.610 1
Vglut2 VGlut-2070 0.215 2
VGlut-2292 0.377 3
VGlut-2293 0.548 4
VGlut-2295 0.370 2
VGlut-2298 0.420 1
VGlut-2299 0.250 1
VGlut-2443 0.348 4
VGlut-A997 0.240 3
VGlut-B001 0.667 6
VGlut-B411 0.422 5
VGlut-B413 0.320 6
VGlut-B612 0.247 3
VGlut-B615 0.464 8
VGlut-B617 0.530 2
Mouse means: Vglut2 0.387 s (14 mice), Calb1 0.193 s (7 mice), rank-sum p = 0.0056
Result: partially reproduced. By the Methods' timing the trough is at 0.416 s for Vglut2, within tolerance, and 0.140 s for Calb1, below the tolerance of 0.153 to 0.187 s; the difference is significant (rank-sum p = 2.1e-08 after multiplying by 3). With the code's timing, which reports 0.02 s later, both means match the paper (0.436 and 0.160 s). Over mouse means the difference holds (0.387 against 0.193 s, p = 0.0056).
fig2jk-pca-variance¶
In a PCA of the cross-correlations of all Vglut2, Calb1 and Anxa1 striatal recordings, PC1 explains 64.2% and PC2 20.1% of the variance.
pure = loco[loco.subtype.isin(["Vglut2", "Calb1", "Anxa1"])]
X = np.stack(pure.cc470.values)
coeff, explained = az.pca_uncentered(X)
print(f"Uncentered PCA of {X.shape[0]} recordings x {X.shape[1]} lags: PC1 {explained[0]:.2f}%, PC2 {explained[1]:.2f}%, together {explained[:2].sum():.2f}%")
Xc = X - X.mean(0)
_, sc, _ = np.linalg.svd(Xc, full_matrices=False)
print("Variant, centered PCA:", np.round(sc[:2] ** 2 / np.sum(sc ** 2) * 100, 1))
rng = np.random.default_rng(20231)
boot = []
for _ in range(200):
i = rng.integers(0, len(X), len(X))
boot.append(az.pca_uncentered(X[i])[1][:2])
boot = np.array(boot)
print("Bootstrap over recordings (200 resamples), 2.5 to 97.5 percentiles: PC1", np.round(np.percentile(boot[:, 0], [2.5, 97.5]), 1),
"PC2", np.round(np.percentile(boot[:, 1], [2.5, 97.5]), 1))
loco["pc1"] = np.stack(loco.cc470.values) @ coeff[:, 0]
loco["pc2"] = np.stack(loco.cc470.values) @ coeff[:, 1]
Uncentered PCA of 160 recordings x 201 lags: PC1 58.32%, PC2 22.09%, together 80.41% Variant, centered PCA: [58. 22.]
Bootstrap over recordings (200 resamples), 2.5 to 97.5 percentiles: PC1 [53.9 63. ] PC2 [18.5 26.3]
Result: reproduced. PC1 explains 58.32% and PC2 22.09% of the variance (80.41% together), against 64.2% and 20.1%; both are inside their tolerance, PC1 by 0.5 points and PC2 by 0.02 points, so the match is marginal. The bootstrap over recordings gives 53.9 to 63.0% for PC1 and 18.5 to 26.3% for PC2. The selection holds 160 recordings, more than the paper's 111.
fig2l-pc-angles¶
The PC1/PC2 angles of the cross-correlations separate the subtypes, with mean angles of 141 degrees for Vglut2, 208 for Calb1 and 239 for Anxa1.
sd = loco[["pc1", "pc2"]].std(ddof=1) # over all selected recordings of the five groups, as the code does
loco["angle"] = np.mod(np.degrees(np.arctan2(loco.pc2 / sd.pc2, loco.pc1 / sd.pc1)), 360)
opened = np.where(loco.angle < 45, loco.angle + 360, loco.angle) - 45
loco["angle_open"] = opened
print("Mean angle (degrees):", {s: round(float(loco[loco.subtype == s].angle.mean()), 1) for s in ORDER})
print("Circular mean (degrees), for reference:", {s: round(float(np.mod(np.degrees(stats.circmean(np.radians(loco[loco.subtype == s].angle))), 360)), 1) for s in ORDER})
pairs = [("Vglut2", "Calb1"), ("Vglut2", "Anxa1"), ("Calb1", "Anxa1")]
for a, b in pairs:
x, y = loco[loco.subtype == a].angle_open, loco[loco.subtype == b].angle_open
print(f"{a} vs {b}: rank-sum p x3 = {p_rs(x, y, 3):.2g}; median difference {y.median() - x.median():.1f} degrees")
print("Per mouse mean angles:")
print(loco[loco.subtype.isin(["Vglut2", "Calb1", "Anxa1"])].groupby(["subtype", "mouse"]).angle.mean().round(0).unstack(0).to_string())
fig, ax = plt.subplots(figsize=(5, 4))
for s in ["Vglut2", "Calb1", "Anxa1", "DAT"]:
d = loco[loco.subtype == s]
ax.scatter(d.pc1, d.pc2, s=14, color=COLOR[s], label=s, alpha=0.8)
ax.axhline(0, color="0.7", lw=0.8); ax.axvline(0, color="0.7", lw=0.8)
ax.set_xlabel("PC1 score"); ax.set_ylabel("PC2 score"); ax.legend(frameon=False, fontsize=8)
plt.show()
print("PC1 range", loco.pc1.min().round(2), "to", loco.pc1.max().round(2), "; PC2 range", loco.pc2.min().round(2), "to", loco.pc2.max().round(2), "; all within the axes")
Mean angle (degrees): {'Vglut2': 140.1, 'Calb1': 208.8, 'Anxa1': 233.3, 'DAT': 229.3, 'Aldh1a1': 196.7}
Circular mean (degrees), for reference: {'Vglut2': 136.4, 'Calb1': 223.0, 'Anxa1': 301.3, 'DAT': 244.4, 'Aldh1a1': 333.6}
Vglut2 vs Calb1: rank-sum p x3 = 7.9e-09; median difference 92.4 degrees
Vglut2 vs Anxa1: rank-sum p x3 = 7.7e-13; median difference 156.9 degrees
Calb1 vs Anxa1: rank-sum p x3 = 1.3e-05; median difference 64.5 degrees
Per mouse mean angles:
subtype Anxa1 Calb1 Vglut2
mouse
Anxa-2559 223.0 NaN NaN
Anxa-A003 120.0 NaN NaN
Anxa-I765 56.0 NaN NaN
Anxa-I767 265.0 NaN NaN
Anxa-J074 264.0 NaN NaN
Anxa-J075 229.0 NaN NaN
Anxa-J077 152.0 NaN NaN
Anxa-J078 293.0 NaN NaN
Anxa-J079 275.0 NaN NaN
Calb-1789 NaN 230.0 NaN
Calb-1790 NaN 169.0 NaN
Calb-1791 NaN 241.0 NaN
Calb-1792 NaN 207.0 NaN
Calb-1877 NaN 221.0 NaN
Calb-2290 NaN 219.0 NaN
Calb-S448 NaN 34.0 NaN
VGlut-2070 NaN NaN 186.0
VGlut-2292 NaN NaN 119.0
VGlut-2293 NaN NaN 101.0
VGlut-2295 NaN NaN 147.0
VGlut-2298 NaN NaN 144.0
VGlut-2299 NaN NaN 146.0
VGlut-2443 NaN NaN 163.0
VGlut-A997 NaN NaN 177.0
VGlut-B001 NaN NaN 93.0
VGlut-B411 NaN NaN 128.0
VGlut-B413 NaN NaN 185.0
VGlut-B612 NaN NaN 182.0
VGlut-B615 NaN NaN 125.0
VGlut-B617 NaN NaN 124.0
PC1 range -1.32 to 1.39 ; PC2 range -0.79 to 0.83 ; all within the axes
Result: reproduced. The mean angles are 140.1 degrees (Vglut2), 208.8 (Calb1) and 233.3 (Anxa1), against 141, 208 and 239, and all three pairs differ (rank-sum p x3 = 7.9e-09, 7.7e-13 and 1.3e-05). Per mouse, the Vglut2 angles lie between 93 and 186 degrees and most Calb1 mice between 169 and 241 (one, with a single recording, at 34); three of nine Anxa1 mice have mean angles below 160 degrees (56, 120 and 152). The arithmetic mean on 0 to 360 degrees follows the code; the circular means are 136.4, 223.0 and 301.3, so the Anxa1 mean depends on where the circle is cut. The tests use angles opened at 45 degrees, as the Methods state.
fig2g-event-counts¶
Triggered averages on large accelerations and decelerations average 135 +/- 77 accelerations and 138 +/- 101 decelerations per recording.
for name in ("AccOn", "DecOn"):
n = loco[f"n_{name}"]
n = n[n >= az.N_EVENTS_MIN]
print(f"{name}: mean {n.mean():.1f} +/- {n.std():.1f} (s.d.) over {len(n)} recordings; by subtype",
loco[loco[f'n_{name}'] >= 3].groupby('subtype')[f'n_{name}'].mean().round(0).to_dict())
AccOn: mean 131.4 +/- 75.1 (s.d.) over 334 recordings; by subtype {'Aldh1a1': 125.0, 'Anxa1': 124.0, 'Calb1': 165.0, 'DAT': 123.0, 'Vglut2': 153.0}
DecOn: mean 129.2 +/- 95.0 (s.d.) over 334 recordings; by subtype {'Aldh1a1': 131.0, 'Anxa1': 86.0, 'Calb1': 192.0, 'DAT': 137.0, 'Vglut2': 148.0}
Result: reproduced. The recordings average 131.4 +/- 75.1 accelerations and 129.2 +/- 95.0 decelerations (334 recordings), against 135 +/- 77 and 138 +/- 101.
ed7f-prob-transient-after-event¶
The probability that a dF/F increase follows an acceleration is 57.5% for Anxa1, and that one follows a deceleration is 62.4% for Vglut2 and 62.3% for Calb1.
for s, ev in (("Vglut2", "DecOn"), ("Calb1", "DecOn"), ("Anxa1", "AccOn"), ("Vglut2", "AccOn"), ("Calb1", "AccOn"), ("Anxa1", "DecOn")):
v = loco[loco.subtype == s][f"pos_{ev}"].dropna()
pm = loco[loco.subtype == s].groupby("mouse")[f"pos_{ev}"].mean()
print(f"{s} {ev}: mean {v.mean():.1f}% over {len(v)} recordings (range {v.min():.0f} to {v.max():.0f}); "
f"mouse means {pm.min():.0f} to {pm.max():.0f}%; signed-rank against 50%: p = {p_sr(v - 50):.2g}")
print("The last three rows are the opposite event type, as a reference.")
Vglut2 DecOn: mean 61.8% over 50 recordings (range 34 to 86); mouse means 50 to 81%; signed-rank against 50%: p = 2.5e-07 Calb1 DecOn: mean 63.9% over 28 recordings (range 36 to 84); mouse means 36 to 75%; signed-rank against 50%: p = 2.3e-06 Anxa1 AccOn: mean 55.8% over 82 recordings (range 35 to 76); mouse means 49 to 67%; signed-rank against 50%: p = 4e-07 Vglut2 AccOn: mean 38.8% over 50 recordings (range 24 to 59); mouse means 24 to 47%; signed-rank against 50%: p = 4e-09 Calb1 AccOn: mean 43.6% over 28 recordings (range 25 to 66); mouse means 36 to 53%; signed-rank against 50%: p = 0.00045 Anxa1 DecOn: mean 38.7% over 82 recordings (range 16 to 70); mouse means 29 to 57%; signed-rank against 50%: p = 3.6e-11 The last three rows are the opposite event type, as a reference.
Result: reproduced. A dF/F increase follows 61.8% of decelerations for Vglut2, 63.9% for Calb1, and 55.8% of accelerations for Anxa1, against 62.4, 62.3 and 57.5%; each is above 50% (p < 3e-06). For the opposite event type the values are below 50% (38.7 to 43.6%).
ed6f-dec-peak-timing¶
In averages triggered on decelerations, the dF/F peak comes later for Vglut2 (0.35 s) than for Calb1 (0.23 s), P = 0.01.
def peak_methods(a):
t0 = az.NEG - 1
return np.argmax(a[t0:t0 + 101]) / az.HZ
def peak_code(a):
s201 = a[5:206] # the code keeps columns 6 to 206 of 210
i = np.argmax(s201[99:]) + 1 # sig = -sig; min(sig(neg0:end))
return min(i, 100) / az.HZ
dd = loco[loco.subtype.isin(["Vglut2", "Calb1"]) & loco.avg_DecOn.map(lambda a: isinstance(a, np.ndarray))].copy()
dd["pk_m"] = dd.avg_DecOn.map(peak_methods)
dd["pk_c"] = dd.avg_DecOn.map(peak_code)
for col, name in (("pk_m", "Methods (main)"), ("pk_c", "code (variant)")):
v, c_ = dd[dd.subtype == "Vglut2"][col], dd[dd.subtype == "Calb1"][col]
print(f"{name}: Vglut2 {v.mean():.3f} s (n={len(v)}), Calb1 {c_.mean():.3f} s (n={len(c_)}), rank-sum p = {p_rs(v, c_):.2g}; difference {v.mean() - c_.mean():.3f} s")
print(dd.groupby(["subtype", "mouse"]).pk_m.mean().round(2).unstack(0).to_string())
print("\nPeak of the mean deceleration-triggered average over recordings (Methods timing), and spread of per-recording peaks:")
paper_pk = {"Vglut2": 0.35, "Calb1": 0.23}
tt = (np.arange(az.NEG + az.POS) - (az.NEG - 1)) / az.HZ
fig, ax = plt.subplots(figsize=(6, 3.2))
for s in ("Vglut2", "Calb1"):
x = dd[dd.subtype == s]
m = np.stack(x.avg_DecOn.values).mean(0)
pk = peak_methods(m)
q = x.pk_m.quantile([0.25, 0.5, 0.75]).round(2).tolist()
lo, hi = paper_pk[s] * 0.9, paper_pk[s] * 1.1
print(f" {s}: peak of the mean {pk:.2f} s; per recording mean {x.pk_m.mean():.3f} s, median {q[1]} s (quartiles {q[0]} to {q[2]}), "
f"skewness {stats.skew(x.pk_m):.2f}; paper tolerance {lo:.3f} to {hi:.3f} s, median inside: {lo <= q[1] <= hi}; "
f"recordings above {hi:.3f} s: {int((x.pk_m > hi).sum())}, below {lo:.3f} s: {int((x.pk_m < lo).sum())}")
print(f" "
f"at 0 s {int((x.pk_m == 0).sum())}, at 1 s {int((x.pk_m >= 1).sum())} of {len(x)}; paper mean {paper_pk[s]} s")
ax.plot(tt, m, color=COLOR[s], lw=2, label=f"{s} (n={len(x)})")
ax.axvline(paper_pk[s], color=COLOR[s], ls="--", lw=1)
ax.axvline(x.pk_m.mean(), color=COLOR[s], ls=":", lw=1.5)
ax.axvline(0, color="0.7", lw=0.8)
ax.set_xlabel("time from deceleration onset (s)"); ax.set_ylabel("normalized dF/F")
ax.set_title("dashed: paper's mean peak time; dotted: mean here")
ax.legend(frameon=False, fontsize=8); plt.tight_layout(); plt.show()
allm = np.concatenate([np.stack(dd[dd.subtype == s].avg_DecOn.values).mean(0) for s in ("Vglut2", "Calb1")])
print("Plotted means range from", round(allm.min(), 3), "to", round(allm.max(), 3), "; all points within the axes")
Methods (main): Vglut2 0.476 s (n=50), Calb1 0.320 s (n=28), rank-sum p = 0.001; difference 0.156 s
code (variant): Vglut2 0.486 s (n=50), Calb1 0.330 s (n=28), rank-sum p = 0.001; difference 0.156 s
subtype Calb1 Vglut2
mouse
Calb-1789 0.35 NaN
Calb-1790 0.14 NaN
Calb-1791 0.37 NaN
Calb-1792 0.43 NaN
Calb-1877 0.23 NaN
Calb-2290 0.28 NaN
Calb-S448 0.00 NaN
VGlut-2070 NaN 0.29
VGlut-2292 NaN 0.39
VGlut-2293 NaN 0.60
VGlut-2295 NaN 0.50
VGlut-2298 NaN 0.65
VGlut-2299 NaN 0.33
VGlut-2443 NaN 0.29
VGlut-A997 NaN 0.39
VGlut-B001 NaN 0.59
VGlut-B411 NaN 0.69
VGlut-B413 NaN 0.44
VGlut-B612 NaN 0.35
VGlut-B615 NaN 0.40
VGlut-B617 NaN 0.76
Peak of the mean deceleration-triggered average over recordings (Methods timing), and spread of per-recording peaks:
Vglut2: peak of the mean 0.44 s; per recording mean 0.476 s, median 0.5 s (quartiles 0.32 to 0.65), skewness -0.07; paper tolerance 0.315 to 0.385 s, median inside: False; recordings above 0.385 s: 31, below 0.315 s: 12
at 0 s 1, at 1 s 1 of 50; paper mean 0.35 s
Calb1: peak of the mean 0.27 s; per recording mean 0.320 s, median 0.25 s (quartiles 0.22 to 0.32), skewness 1.48; paper tolerance 0.207 to 0.253 s, median inside: True; recordings above 0.253 s: 13, below 0.207 s: 4
at 0 s 2, at 1 s 1 of 28; paper mean 0.23 s
Plotted means range from 0.033 to 0.095 ; all points within the axes
Result: partially reproduced. The dF/F peak after decelerations comes later for Vglut2 than for Calb1 (0.476 against 0.320 s, p = 0.001), as stated, but both means are about 0.1 s later than the paper's 0.35 and 0.23 s and outside tolerance. The code's timing adds 0.01 s and does not change this. The two subtypes differ in how. For Vglut2 the whole distribution is later: the median is 0.5 s, the peak of the mean triggered average 0.44 s, and 31 of 50 recordings peak after 0.385 s. For Calb1 the per-recording peaks are skewed to the right (skewness 1.48, mean 0.320 s at the upper quartile of 0.32 s): the median of 0.25 s and the peak of the mean (0.27 s) are close to the paper's 0.23 s, the median within tolerance, and the gap in the mean comes from a tail of late peaks (13 of 28 above 0.253 s). Recordings at the window edges are few (2 of 50 and 3 of 28 at 0 or 1 s). The cause of the Vglut2 difference is not found.
fig5ab-transient-scaling¶
Transient amplitude grows with the size of the deceleration (Vglut2 214%, Calb1 243%) or acceleration (Anxa1 206%) from the smallest to the largest quintile (P = 0.01, 0.002, 0.002).
def scaling_table(col):
rows = []
for s, kind in (("Vglut2", "dec"), ("Calb1", "dec"), ("Anxa1", "acc"), ("Vglut2", "acc"), ("Calb1", "acc"), ("Anxa1", "dec")):
amps = np.array([a[kind] for a in loco[loco.subtype == s][col] if a.get(kind) is not None])
q1, q5 = amps[:, 0], amps[:, 4]
ok = ~(np.isnan(q1) | np.isnan(q5))
rows.append({"subtype": s, "events": kind, "n_rec": len(amps), "n_paired": int(ok.sum()),
**{f"q{i + 1}": np.nanmean(amps[:, i]) * 100 for i in range(5)},
"ratio_pct": np.nanmean(q5) / np.nanmean(q1) * 100,
"ratio_paired_pct": q5[ok].mean() / q1[ok].mean() * 100,
"p_x3": min(1, stats.wilcoxon(q1[ok], q5[ok]).pvalue * 3)})
return pd.DataFrame(rows)
print("Mean transient amplitude per quintile of event size (normalized dF/F, 0 to 100 scale).")
print("Main: the recomputed onsets that are also stored, which carry the code's exclusion after every reward:")
print(scaling_table("size_amps").round(3).to_string(index=False))
print("\nVariant: every recomputed onset (5 s excluded after consumed rewards only):")
print(scaling_table("size_amps_all").round(3).to_string(index=False))
print("The first three rows of each table are the claim; the others use the opposite event type.")
Mean transient amplitude per quintile of event size (normalized dF/F, 0 to 100 scale). Main: the recomputed onsets that are also stored, which carry the code's exclusion after every reward: subtype events n_rec n_paired q1 q2 q3 q4 q5 ratio_pct ratio_paired_pct p_x3 Vglut2 dec 50 50 3.238 4.314 4.918 6.540 7.812 241.269 241.269 0.002 Calb1 dec 28 28 3.247 6.268 6.836 7.053 7.995 246.233 246.233 0.000 Anxa1 acc 82 81 2.726 3.351 5.852 4.546 6.322 231.883 231.883 0.001 Vglut2 acc 50 50 -1.977 -4.256 -4.376 -4.142 -4.993 252.514 252.514 0.142 Calb1 acc 28 28 -2.527 -2.995 -2.720 -3.572 -1.508 59.646 59.646 1.000 Anxa1 dec 82 81 -2.671 -4.746 -4.369 -6.184 -5.419 202.883 202.883 0.051 Variant: every recomputed onset (5 s excluded after consumed rewards only):
subtype events n_rec n_paired q1 q2 q3 q4 q5 ratio_pct ratio_paired_pct p_x3 Vglut2 dec 50 50 3.258 4.302 4.981 6.512 7.797 239.326 239.326 0.002 Calb1 dec 28 28 3.035 6.249 6.795 7.245 8.058 265.551 265.551 0.000 Anxa1 acc 82 81 2.731 3.303 5.835 4.671 6.070 222.262 222.262 0.001 Vglut2 acc 50 50 -1.954 -4.438 -4.356 -4.122 -4.968 254.268 254.268 0.135 Calb1 acc 28 28 -2.117 -2.813 -2.482 -3.154 -1.656 78.192 78.192 1.000 Anxa1 dec 82 81 -2.636 -4.773 -4.475 -6.228 -5.387 204.337 204.337 0.068 The first three rows of each table are the claim; the others use the opposite event type.
Result: partially reproduced. With the recomputed onsets that are also stored, which carry the code's exclusion after every reward, transient amplitude grows from the smallest to the largest quintile in all three subtypes, by 241% (Vglut2 decelerations), 246% (Calb1) and 232% (Anxa1 accelerations), with p x3 = 0.002, below 0.001 and 0.001. Calb1 is within tolerance of 243%; Vglut2 is above 192.6 to 235.4% and Anxa1 above 185.4 to 226.6%. With every recomputed onset the ratios are 239, 266 and 222% (Vglut2 outside tolerance). The ratio over the recordings with both quintiles equals the ratio over all recordings to the precision printed.
fig4g-reward-response¶
Vglut2 and Calb1 axons respond to unexpected rewards (means 7.9 and 12.4 normalized dF/F s), as do DAT recordings (5.9).
def resp_table(d, name, groups=("Vglut2", "Calb1", "Anxa1", "DAT")):
rows = []
for s in ORDER:
x = d[d.subtype == s]
x = x[x[f"n_{name}"] >= az.N_EVENTS_MIN]
m, c_ = x[f"resp_methods_{name}"] * 100, x[f"resp_code_{name}"] * 100
pm = (x.groupby("mouse")[f"resp_methods_{name}"].mean() * 100)
f = 4 if s in groups else 1
rows.append({"subtype": s, "n": len(x), "mice": mice(x), "mean": m.mean(), "sem": m.sem(),
"p_x4": p_sr(m, f), "frac_pos": (m > 0).mean(), "code_mean": c_.mean(), "code_p_x4": p_sr(c_, f),
"mouse_means_pos": f"{int((pm > 0).sum())}/{len(pm)}"})
return pd.DataFrame(rows).set_index("subtype")
rew = resp_table(stim, "Reward")
print("Reward responses (normalized dF/F s, 0 to 100 scale); Aldh1a1 not in the Bonferroni family:")
print(rew.round(4).to_string())
print("\nPaper: Vglut2 7.9, Calb1 12.4, Anxa1 -0.5, DAT 5.9")
print("Paper P for rewards: Figure 4g legend Vglut2 2e-5, Calb1 0.001; Results text 'Calb1+ and Vglut2+ ... P = 2e-5 and P = 0.001, respectively'")
print("Here (p x4):", rew.loc[["Vglut2", "Calb1"], "p_x4"].map(lambda v: f"{v:.2g}").to_dict())
Reward responses (normalized dF/F s, 0 to 100 scale); Aldh1a1 not in the Bonferroni family:
n mice mean sem p_x4 frac_pos code_mean code_p_x4 mouse_means_pos
subtype
Vglut2 28 11 8.7738 0.7023 0.0000 1.0000 7.8999 0.0000 11/11
Calb1 17 8 13.7288 1.7691 0.0002 0.9412 12.4012 0.0001 8/8
Anxa1 51 8 -0.4179 0.5887 0.1215 0.3725 -0.3927 0.1430 3/8
DAT 63 11 6.6895 0.9566 0.0000 0.7302 5.8776 0.0000 9/11
Aldh1a1 30 8 3.8301 1.2725 0.0106 0.6667 3.3996 0.0175 7/8
Paper: Vglut2 7.9, Calb1 12.4, Anxa1 -0.5, DAT 5.9
Paper P for rewards: Figure 4g legend Vglut2 2e-5, Calb1 0.001; Results text 'Calb1+ and Vglut2+ ... P = 2e-5 and P = 0.001, respectively'
Here (p x4): {'Vglut2': '3e-08', 'Calb1': '0.00018'}
Result: partially reproduced. Vglut2, Calb1 and DAT recordings respond to rewards (means 8.77, 13.73 and 6.69 with the Methods' windows, p x4 < 0.001), as stated, but each mean is just above its tolerance (upper bounds 8.69, 13.64 and 6.49). With the code's windows the means are 7.90, 12.40 and 5.88, equal to the paper's 7.9, 12.4 and 5.9 to the decimal, which suggests the paper's reward means were computed with the code's windows. For air puffs this holds for Vglut2, Calb1 and Anxa1 but not for DAT (fig4g-airpuff-response). Every Vglut2 and Calb1 mouse has a positive mean response.
fig4g-anxa1-no-reward¶
Reward responses are almost entirely absent from Anxa1 axons (mean -0.5, P = 0.1, not significant).
a = stim[(stim.subtype == "Anxa1") & (stim.n_Reward >= az.N_EVENTS_MIN)]
m = a.resp_methods_Reward * 100
print(f"Anxa1 reward responses: n = {len(m)} from {mice(a)} mice, mean {m.mean():.2f} (sem {m.sem():.2f}), median {m.median():.2f}, "
f"signed-rank p x4 = {p_sr(m, 4):.3g}; code windows mean {(a.resp_code_Reward * 100).mean():.2f}, p x4 = {p_sr(a.resp_code_Reward * 100, 4):.3g}")
print(f"Recordings with a response above 0: {(m > 0).mean():.2f}; above 5 (about the DAT mean): {int((m > 5).sum())}")
print("Per mouse:", (a.groupby("mouse").resp_methods_Reward.mean() * 100).round(2).to_dict())
Anxa1 reward responses: n = 51 from 8 mice, mean -0.42 (sem 0.59), median -1.18, signed-rank p x4 = 0.121; code windows mean -0.39, p x4 = 0.143
Recordings with a response above 0: 0.37; above 5 (about the DAT mean): 4
Per mouse: {'Anxa-2558': -2.09, 'Anxa-2559': 0.86, 'Anxa-I767': 0.03, 'Anxa-J074': 0.67, 'Anxa-J075': -1.3, 'Anxa-J077': -0.78, 'Anxa-J078': -1.35, 'Anxa-J079': -2.53}
Result: reproduced. Anxa1 reward responses do not differ from 0 by the paper's test (mean -0.42, p x4 = 0.12; code windows -0.39, p x4 = 0.14), against the paper's -0.5 and P = 0.1. 37% of recordings are above 0 and 4 of 51 above 5.
fig4g-airpuff-response¶
Air puffs evoke increases in Vglut2 (15.8), Calb1 (5.3) and DAT (5.3) axons and a decrease in Anxa1 axons (-3.7).
puff = resp_table(stim, "AirPuff")
print(puff.round(4).to_string())
print("\nPaper: Vglut2 15.8, Calb1 5.3, Anxa1 -3.7, DAT 5.3")
fig, ax = plt.subplots(1, 2, figsize=(9, 3.2), sharey=True)
for a_, tab, title in ((ax[0], rew, "Reward"), (ax[1], puff, "Air puff")):
for i, s in enumerate(ORDER):
a_.bar(i, tab.loc[s, "mean"], yerr=tab.loc[s, "sem"], color=COLOR[s], width=0.6)
a_.axhline(0, color="0.5", lw=0.8); a_.set_xticks(range(5), ORDER, fontsize=8); a_.set_title(title)
ax[0].set_ylabel("response, normalized dF/F s (mean +/- sem)")
plt.tight_layout(); plt.show()
print("Bar heights range", round(min(rew['mean'].min(), puff['mean'].min()), 2), "to", round(max(rew['mean'].max(), puff['mean'].max()), 2), "; all within the axes")
dp = stim[(stim.subtype == "DAT") & (stim.n_AirPuff >= az.N_EVENTS_MIN)].assign(
methods=lambda x: x.resp_methods_AirPuff * 100, code=lambda x: x.resp_code_AirPuff * 100)
print("\nDAT air puff responses per mouse (normalized dF/F s):")
print(dp.groupby("mouse").agg(n=("path", "size"), location=("location", lambda v: ",".join(sorted(set(v)))),
methods=("methods", "mean"), code=("code", "mean")).round(2).to_string())
print("\nDAT air puff responses by location label:")
print(dp.groupby("location").agg(n=("path", "size"), mice=("mouse", "nunique"), methods=("methods", "mean"), code=("code", "mean")).round(2).to_string())
pm = dp.groupby("mouse").code.mean()
print(f"Mean over mouse means (code windows) {pm.mean():.2f}; median over recordings {dp.code.median():.2f}")
n mice mean sem p_x4 frac_pos code_mean code_p_x4 mouse_means_pos subtype Vglut2 29 12 16.7749 1.0583 0.0000 1.0000 15.7849 0.0000 12/12 Calb1 17 8 5.4098 0.9405 0.0026 0.9412 5.3066 0.0026 7/8 Anxa1 57 8 -3.6555 0.4538 0.0000 0.1404 -3.6576 0.0000 0/8 DAT 69 11 2.6933 0.7334 0.0091 0.6087 2.4585 0.0134 8/11 Aldh1a1 30 8 -1.2085 0.9702 0.2286 0.3000 -1.2405 0.2286 3/8 Paper: Vglut2 15.8, Calb1 5.3, Anxa1 -3.7, DAT 5.3
Bar heights range -3.66 to 16.77 ; all within the axes
DAT air puff responses per mouse (normalized dF/F s):
n location methods code
mouse
Dat-2111 12 DS 2.48 2.35
Dat-2114 11 DS 6.33 6.03
Dat-2115 8 DS 5.83 5.57
Dat-D001 4 DS -3.37 -3.83
Dat-D002 5 DS 4.67 3.32
Dat-D087 3 DS,TS -0.26 -0.26
Dat-D093 1 DS 3.42 3.37
Dat-D095 1 DS 6.48 5.83
Dat-D171 2 DS,TS 4.93 4.82
Dat-F554 7 DS 2.74 2.89
Dat-G557 15 DS -0.55 -0.60
DAT air puff responses by location label:
n mice methods code
location
DS 66 11 2.76 2.52
TS 3 2 1.12 1.13
Mean over mouse means (code windows) 2.68; median over recordings 0.69
Result: partially reproduced. Air puffs increase Vglut2 (16.77), Calb1 (5.41) and DAT (2.69) signals and decrease Anxa1 signals (-3.66), all with p x4 < 0.01. Vglut2, Calb1 and Anxa1 are within tolerance of 15.8, 5.3 and -3.7 (and with the code's windows 15.78, 5.31, -3.66); the DAT mean is half the paper's 5.3 under either window (2.69 and 2.46). Per mouse, the DAT means (code windows) run from -3.83 to 6.03, three of 11 mice are at or below 0, among them Dat-G557 with 15 recordings (-0.60); the mean over mice is 2.68 and the median over recordings 0.69. The cause of the difference from the paper is not found.
fig4h-reward-vs-airpuff¶
Vglut2 axons respond more to air puffs than to rewards, and Calb1 axons more to rewards than to air puffs.
both = stim[(stim.n_Reward >= 3) & (stim.n_AirPuff >= 3)]
for s in ORDER:
x = both[both.subtype == s]
d = (x.resp_methods_AirPuff - x.resp_methods_Reward) * 100
dc = (x.resp_code_AirPuff - x.resp_code_Reward) * 100
print(f"{s}: n = {len(x)}, mean air puff minus reward {d.mean():.2f}, fraction puff > reward {(d > 0).mean():.2f}, "
f"p x4 = {p_sr(d, 4):.2g}; code windows {dc.mean():.2f}, p x4 = {p_sr(dc, 4):.2g}")
Vglut2: n = 28, mean air puff minus reward 8.03, fraction puff > reward 0.86, p x4 = 9.1e-06; code windows 7.90, p x4 = 1.6e-06 Calb1: n = 17, mean air puff minus reward -8.32, fraction puff > reward 0.18, p x4 = 0.0034; code windows -7.09, p x4 = 0.0054 Anxa1: n = 51, mean air puff minus reward -3.11, fraction puff > reward 0.22, p x4 = 0.00029; code windows -3.18, p x4 = 3.9e-05 DAT: n = 59, mean air puff minus reward -3.81, fraction puff > reward 0.29, p x4 = 0.0014; code windows -3.16, p x4 = 0.0025 Aldh1a1: n = 30, mean air puff minus reward -5.04, fraction puff > reward 0.13, p x4 = 2.1e-05; code windows -4.64, p x4 = 5.7e-06
Result: reproduced. Vglut2 recordings respond more to air puffs than to rewards (mean difference 8.03, 86% of recordings, p x4 = 9.1e-06) and Calb1 recordings more to rewards (-8.32, 18% favour air puffs, p x4 = 0.0034).
fig4i-reward-size¶
Calb1 axons respond more to large than to small rewards (mean difference 3.9, P = 9e-3), as do DAT recordings (1.9, P = 8e-5), while Vglut2 (0.9, P = 0.6) and Anxa1 (0.04, P = 1) do not.
sz = stim[(stim.n_RewardShort >= 3) & (stim.n_RewardLong >= 3)]
for s in ORDER:
x = sz[sz.subtype == s]
d = (x.resp_methods_RewardLong - x.resp_methods_RewardShort) * 100
dc = (x.resp_code_RewardLong - x.resp_code_RewardShort) * 100
pm = x.assign(d=d).groupby("mouse").d.mean()
print(f"{s}: n = {len(x)} from {mice(x)} mice, large minus small {d.mean():.2f} (sem {d.sem():.2f}), p x4 = {p_sr(d, 4):.2g}; "
f"code windows {dc.mean():.2f}, p x4 = {p_sr(dc, 4):.2g}; mice with mean > 0: {int((pm > 0).sum())}/{len(pm)}")
print("Paper n: Vglut2 25, Calb1 14, Anxa1 42, DAT 55")
print("\nVariant: minimum number of rewards of each size (Methods windows):")
for k in (1, 2, 3):
x = stim[(stim.n_RewardShort >= k) & (stim.n_RewardLong >= k)]
out = {}
for s in ("Vglut2", "Calb1", "Anxa1", "DAT"):
y = x[x.subtype == s]
d = (y.resp_methods_RewardLong - y.resp_methods_RewardShort) * 100
out[s] = f"n={len(y)} diff={d.mean():.2f} p_x4={p_sr(d, 4):.2g}"
print(f" at least {k}:", out)
Vglut2: n = 25 from 11 mice, large minus small 1.23 (sem 0.56), p x4 = 0.22; code windows 0.95, p x4 = 0.63; mice with mean > 0: 7/11 Calb1: n = 14 from 6 mice, large minus small 4.43 (sem 0.93), p x4 = 0.0049; code windows 3.93, p x4 = 0.0093; mice with mean > 0: 5/6
Anxa1: n = 37 from 8 mice, large minus small 0.22 (sem 0.61), p x4 = 1; code windows 0.04, p x4 = 1; mice with mean > 0: 4/8
DAT: n = 50 from 10 mice, large minus small 2.33 (sem 0.41), p x4 = 3.7e-06; code windows 1.83, p x4 = 2.8e-05; mice with mean > 0: 9/10
Aldh1a1: n = 28 from 8 mice, large minus small 1.94 (sem 0.76), p x4 = 0.087; code windows 1.44, p x4 = 0.23; mice with mean > 0: 5/8
Paper n: Vglut2 25, Calb1 14, Anxa1 42, DAT 55
Variant: minimum number of rewards of each size (Methods windows):
at least 1: {'Vglut2': 'n=26 diff=1.23 p_x4=0.22', 'Calb1': 'n=16 diff=4.43 p_x4=0.0049', 'Anxa1': 'n=47 diff=0.22 p_x4=1', 'DAT': 'n=58 diff=2.33 p_x4=3.7e-06'}
at least 2: {'Vglut2': 'n=26 diff=1.23 p_x4=0.22', 'Calb1': 'n=15 diff=4.43 p_x4=0.0049', 'Anxa1': 'n=40 diff=0.22 p_x4=1', 'DAT': 'n=53 diff=2.33 p_x4=3.7e-06'}
at least 3: {'Vglut2': 'n=25 diff=1.23 p_x4=0.22', 'Calb1': 'n=14 diff=4.43 p_x4=0.0049', 'Anxa1': 'n=37 diff=0.22 p_x4=1', 'DAT': 'n=50 diff=2.33 p_x4=3.7e-06'}
Result: partially reproduced. The significance pattern holds: Calb1 (4.43, p x4 = 0.0049) and DAT (2.33, p x4 = 3.7e-06) respond more to large rewards, Vglut2 (1.23, p x4 = 0.22) and Anxa1 (0.22, p x4 = 1) do not. The Calb1 and DAT differences are above tolerance with the Methods' windows; with the code's windows they are 3.93 and 1.83, within it. Fewer Anxa1 and DAT recordings have 3 or more rewards of each size (37 and 50) than the paper uses (42 and 55); requiring 1 or 2 of each gives 47 and 58, or 40 and 53, so no minimum gives the paper's counts, and the differences and p-values are the same under all three.
ed8f-rest-reward-same¶
Reward responses of Vglut2 and Calb1 axons are the same for rewards delivered at rest as for the others (P = 1 for both).
rr = stim[(stim.n_RewardRest >= 3) & (stim.n_RewardNoRest >= 3)]
for s in ORDER:
x = rr[rr.subtype == s]
d = (x.resp_methods_RewardRest - x.resp_methods_RewardNoRest) * 100
dc = (x.resp_code_RewardRest - x.resp_code_RewardNoRest) * 100
print(f"{s}: n = {len(x)} from {mice(x)} mice, at rest minus not at rest {d.mean():.2f} (sem {d.sem():.2f}; not-at-rest mean "
f"{(x.resp_methods_RewardNoRest * 100).mean():.2f}), p x4 = {p_sr(d, 4):.2g}; code windows {dc.mean():.2f}, p x4 = {p_sr(dc, 4):.2g}")
print("Paper n (rest): Vglut2 8 (6 mice), Calb1 10 (6), Anxa1 42 (8), DAT 42 (10)")
rc = [r["rest_check"] for r in L if "rest_check" in r]
off = np.concatenate([c["start_offset_s"] for c in rc])
dr = np.concatenate([c["dur_rest_s"] for c in rc]); dw = np.concatenate([c["dur_reward_s"] for c in rc])
print(f"\nRest rows in {len(rc)} files: start minus matched reward start, s: min {off.min():.3f}, max {off.max():.3f}; "
f"duration rest minus matched reward, |max| {np.abs(dr - dw).max():.3f} s")
print("Durations (s), rest rows:", pd.Series(dr).describe().round(3).to_dict())
print("Durations (s), matched reward rows:", pd.Series(dw).describe().round(3).to_dict())
mism = [c for c in rc if c["norest"] != c["rewards"] - c["rest"]]
print(f"RewardNoRest onsets {sum(c['norest'] for c in rc)} against rewards minus rest {sum(c['rewards'] - c['rest'] for c in rc)}; "
f"files where they differ: {len(mism)}; rest samples outside a reward: {sum(c['rest_outside_reward'] for c in rc)}")
Vglut2: n = 8 from 3 mice, at rest minus not at rest 0.45 (sem 1.25; not-at-rest mean 9.38), p x4 = 1; code windows 0.36, p x4 = 1 Calb1: n = 8 from 4 mice, at rest minus not at rest -0.15 (sem 0.71; not-at-rest mean 15.35), p x4 = 1; code windows 0.38, p x4 = 1
Anxa1: n = 36 from 8 mice, at rest minus not at rest -1.39 (sem 0.74; not-at-rest mean 0.09), p x4 = 0.43; code windows -1.40, p x4 = 0.45
DAT: n = 33 from 6 mice, at rest minus not at rest -1.13 (sem 0.94; not-at-rest mean 6.93), p x4 = 0.4; code windows -1.00, p x4 = 0.37
Aldh1a1: n = 15 from 6 mice, at rest minus not at rest -1.64 (sem 1.66; not-at-rest mean 4.14), p x4 = 1; code windows -2.04, p x4 = 1
Paper n (rest): Vglut2 8 (6 mice), Calb1 10 (6), Anxa1 42 (8), DAT 42 (10)
Rest rows in 207 files: start minus matched reward start, s: min 0.010, max 0.010; duration rest minus matched reward, |max| 0.000 s
Durations (s), rest rows: {'count': 1393.0, 'mean': 0.035, 'std': 0.031, 'min': 0.01, '25%': 0.01, '50%': 0.02, '75%': 0.07, 'max': 0.34}
Durations (s), matched reward rows: {'count': 1393.0, 'mean': 0.035, 'std': 0.031, 'min': 0.01, '25%': 0.01, '50%': 0.02, '75%': 0.07, 'max': 0.34}
RewardNoRest onsets 2647 against rewards minus rest 2647; files where they differ: 0; rest samples outside a reward: 0
Result: reproduced. Responses to rewards at rest do not differ from the others by the paper's test for Vglut2 (difference 0.45, p x4 = 1, n = 8) or Calb1 (-0.15, p x4 = 1, n = 8), both small next to the not-at-rest means of 9.38 and 15.35. Anxa1 and DAT give p x4 = 0.43 and 0.4, as the paper's 0.4 and 0.2 are not significant. The counts differ from the paper's for three groups: Calb1 8 against 10, Anxa1 36 against 42 and DAT 33 against 42, and the Vglut2 recordings come from 3 mice against 6 (Calb1 4 against 6, DAT 6 against 10); the cause is not found. The rest rows sit exactly 0.01 s after the reward rows they tag, with the same durations, so removing them from the rewards leaves 2647 onsets, equal to the rewards minus the rest rewards in every file.
fig6c-kmeans¶
Unsupervised k-means clustering on locomotion (PC1, PC2 scores), reward and air puff responses assigns recordings to clusters that match the genetic subtypes with 91% accuracy.
kd = loco[loco.has_events & loco.subtype.isin(["Vglut2", "Calb1", "Anxa1"]) &
(loco.n_Reward >= 3) & (loco.n_AirPuff >= 3)].copy()
print("Recordings with all four measures:", kd.subtype.value_counts().reindex(["Vglut2", "Calb1", "Anxa1"]).to_dict(),
"(paper 16, 11, 16); mice:", kd.groupby("subtype").mouse.nunique().to_dict())
lab = kd.subtype.map({"Vglut2": 1, "Calb1": 2, "Anxa1": 4}).values
def kmeans_accuracy(F, seed, n_init=1):
Z = (F - F.mean(0)) / F.std(0, ddof=1)
k = KMeans(3, init="k-means++", n_init=n_init, random_state=seed).fit_predict(Z) + 1
k2 = k.copy()
for s in (1, 2, 4): # the code: median cluster label of each subtype
m = np.median(k[lab == s])
k2[k == m] = s
per = {s: float(np.mean(k2[lab == s] == s)) for s in (1, 2, 4)}
return float(np.mean(k2 == lab)), per
for resp, name in (("resp_methods", "Methods windows (main)"), ("resp_code", "code windows")):
F = np.column_stack([kd.pc1, kd.pc2, kd[f"{resp}_Reward"], kd[f"{resp}_AirPuff"]])
acc1 = [kmeans_accuracy(F, seed)[0] for seed in range(100)]
acc10 = [kmeans_accuracy(F, seed, 10)[0] for seed in range(20)]
print(f"{name}: median accuracy over 100 seeds {np.median(acc1) * 100:.1f}% (range {min(acc1) * 100:.0f} to {max(acc1) * 100:.0f}%); "
f"with 10 starts, median over 20 seeds {np.median(acc10) * 100:.1f}%")
F = np.column_stack([kd.pc1, kd.pc2, kd.resp_methods_Reward, kd.resp_methods_AirPuff])
seeds = [kmeans_accuracy(F, seed) for seed in range(100)]
per = pd.DataFrame([s[1] for s in seeds]).rename(columns={1: "Vglut2", 2: "Calb1", 4: "Anxa1"})
print("Per-subtype accuracy, median over seeds:", (per.median() * 100).round(1).to_dict(), "; chance 33%")
from scipy.optimize import linear_sum_assignment
def kmeans_accuracy_matched(F, seed):
Z = (F - F.mean(0)) / F.std(0, ddof=1)
k = KMeans(3, init="k-means++", n_init=1, random_state=seed).fit_predict(Z)
conf = np.array([[np.sum((k == c) & (lab == s)) for s in (1, 2, 4)] for c in range(3)])
r, c = linear_sum_assignment(-conf)
return conf[r, c].sum() / len(lab)
acc_m = [kmeans_accuracy(F, seed)[0] for seed in range(100)]
acc_h = [kmeans_accuracy_matched(F, seed) for seed in range(100)]
dup = sum(len(set(np.median(KMeans(3, init="k-means++", n_init=1, random_state=seed).fit_predict((F - F.mean(0)) / F.std(0, ddof=1))[lab == s])
for s in (1, 2, 4))) < 3 for seed in range(100))
print(f"Variant, one-to-one mapping of clusters to subtypes: median accuracy {np.median(acc_h) * 100:.1f}% "
f"(range {min(acc_h) * 100:.0f} to {max(acc_h) * 100:.0f}%); code mapping {np.median(acc_m) * 100:.1f}%. "
f"Seeds where two subtypes share a median cluster: {dup} of 100")
Recordings with all four measures: {'Vglut2': 19, 'Calb1': 14, 'Anxa1': 35} (paper 16, 11, 16); mice: {'Anxa1': 6, 'Calb1': 6, 'Vglut2': 9}
Methods windows (main): median accuracy over 100 seeds 92.6% (range 41 to 93%); with 10 starts, median over 20 seeds 92.6% code windows: median accuracy over 100 seeds 92.6% (range 44 to 93%); with 10 starts, median over 20 seeds 92.6%
Per-subtype accuracy, median over seeds: {'Vglut2': 89.5, 'Calb1': 92.9, 'Anxa1': 94.3} ; chance 33%
Variant, one-to-one mapping of clusters to subtypes: median accuracy 92.6% (range 54 to 93%); code mapping 92.6%. Seeds where two subtypes share a median cluster: 5 of 100
Result: reproduced. The median k-means accuracy over 100 seeds is 92.6% (Vglut2 89.5%, Calb1 92.9%, Anxa1 94.3%), against 91% and chance at 33%; single starts range from 41% to 93%, and 10 starts give 92.6%. The code's windows give the same. The clustering uses 19, 14 and 35 recordings, more than the paper's 16, 11 and 16, as the locomotion selection is larger. Mapping clusters to subtypes one to one, in place of the code's median label, gives the same median, 92.6%; the code's mapping sends two subtypes to one cluster in 5 of 100 seeds.
fig7f-soma-axon-correlation¶
Within genetic subtypes, simultaneously recorded SNc somas and striatal axons are highly correlated (mean peak cross-correlation 0.65 Vglut2, 0.67 Calb1, 0.58 Anxa1), more than in DAT recordings (0.37).
du = pd.DataFrame([{"path": r["path"], "subtype": az.SUBTYPE[r["prefix"]], "mouse": r["subject"].rstrip("L"), **r["dual"]}
for r in recs if "dual" in r])
print("Dual SNc + striatum files with both fibers above SNR 10:", du.subtype.value_counts().to_dict())
du["ok"] = (du.rg405 < az.RG405_MAX) & ~du.bad_any
du = du[du.ok]
du["peak405"] = du.cc405.map(np.max)
t = du.groupby("subtype").agg(n=("peak", "size"), mice=("mouse", "nunique"), mean_peak=("peak", "mean"), sem=("peak", "sem"),
mean_peak405=("peak405", "mean")).reindex(ORDER)
dat = du[du.subtype == "DAT"].peak
t["p_vs_DAT"] = [p_rs(du[du.subtype == s].peak, dat) if s != "DAT" else np.nan for s in t.index]
t["p_vs_DAT_x3"] = np.minimum(1, t.p_vs_DAT * 3)
t["paper_mean"] = pd.Series({"Vglut2": 0.65, "Calb1": 0.67, "Anxa1": 0.58, "DAT": 0.37})
t["paper_n"] = pd.Series({"Vglut2": 11, "Calb1": 5, "Anxa1": 43, "DAT": 35})
print(t.round(4).to_string())
print("\nPer mouse mean peak:")
print(du.groupby(["subtype", "mouse"]).peak.mean().round(2).to_string())
Dual SNc + striatum files with both fibers above SNR 10: {'Anxa1': 52, 'DAT': 37, 'Aldh1a1': 32, 'Vglut2': 12, 'Calb1': 5}
n mice mean_peak sem mean_peak405 p_vs_DAT p_vs_DAT_x3 paper_mean paper_n
subtype
Vglut2 11 4 0.6504 0.0338 -0.0110 0.0001 0.0003 0.65 11.0
Calb1 5 2 0.6717 0.0175 0.0171 0.0001 0.0003 0.67 5.0
Anxa1 40 7 0.5295 0.0277 0.0064 0.0006 0.0017 0.58 43.0
DAT 33 5 0.3771 0.0327 -0.0086 NaN NaN 0.37 35.0
Aldh1a1 28 9 0.4913 0.0342 0.0129 0.0273 0.0819 NaN NaN
Per mouse mean peak:
subtype mouse
Aldh1a1 Raldh-1572 0.33
Raldh-1573 0.22
Raldh-1847 0.44
Raldh-1849 0.57
Raldh-2071 0.56
Raldh-2076 0.54
Raldh-2276 0.49
Raldh-2281 0.64
Raldh-2403 0.05
Anxa1 Anxa-I765 0.30
Anxa-I767 0.45
Anxa-J074 0.59
Anxa-J075 0.66
Anxa-J077 0.40
Anxa-J078 0.55
Anxa-J079 0.61
Calb1 Calb-1789 0.65
Calb-1792 0.68
DAT Dat-2111 0.31
Dat-2114 0.32
Dat-2115 0.38
Dat-F554 0.21
Dat-G557 0.47
Vglut2 VGlut-2292 0.78
VGlut-2293 0.67
VGlut-2295 0.55
VGlut-2298 0.50
Result: reproduced. The mean peak soma-axon cross-correlation is 0.650 (Vglut2), 0.672 (Calb1), 0.530 (Anxa1) and 0.377 (DAT), all within tolerance of 0.65, 0.67, 0.58 and 0.37, and each subtype exceeds DAT (rank-sum p x3 = 0.0003, 0.0003 and 0.0017). The 405 nm control peaks are near 0 (-0.011 to 0.017). The counts are 11, 5, 40 and 33 pairs against 11, 5, 43 and 35.
ed9c-snc-stimulus-responses¶
Like their axons, Vglut2 and Calb1 somas in SNc respond to rewards and air puffs, while reward responses are largely undetectable in Anxa1 somas.
snc_stim = df[df.snc & df.stim].copy()
print("SNc reward responses:")
print(resp_table(snc_stim, "Reward").round(4).to_string())
print("\nSNc air puff responses:")
print(resp_table(snc_stim, "AirPuff").round(4).to_string())
print("\nPaper n: reward 25, 10, 23, 39; air puff 25, 10, 25, 47 (Vglut2, Calb1, Anxa1, DAT)")
SNc reward responses:
n mice mean sem p_x4 frac_pos code_mean code_p_x4 mouse_means_pos
subtype
Vglut2 25 9 10.1963 0.8714 0.0000 1.0000 8.9365 0.0000 9/9
Calb1 10 5 12.3392 1.7247 0.0078 1.0000 11.0999 0.0078 5/5
Anxa1 22 5 0.5167 0.6528 1.0000 0.6364 0.5426 1.0000 2/5
DAT 39 7 8.1553 0.7889 0.0000 1.0000 6.9573 0.0000 7/7
Aldh1a1 40 11 9.3690 0.9088 0.0000 0.9250 8.4438 0.0000 11/11
SNc air puff responses:
n mice mean sem p_x4 frac_pos code_mean code_p_x4 mouse_means_pos
subtype
Vglut2 25 9 14.7174 1.3089 0.0000 1.0000 13.5306 0.0000 9/9
Calb1 10 5 4.6760 0.9423 0.0078 1.0000 4.4292 0.0078 5/5
Anxa1 25 5 -0.9483 0.6380 1.0000 0.4000 -0.8557 1.0000 2/5
DAT 47 8 3.4957 0.7323 0.0000 0.7872 3.2767 0.0000 8/8
Aldh1a1 41 11 5.7725 1.2837 0.0001 0.7561 5.4357 0.0001 10/11
Paper n: reward 25, 10, 23, 39; air puff 25, 10, 25, 47 (Vglut2, Calb1, Anxa1, DAT)
Result: reproduced. In SNc, Vglut2 somas respond to rewards (10.20) and air puffs (14.72) and Calb1 somas too (12.34 and 4.68), all with p x4 < 0.01, while the Anxa1 reward response is not significant (0.52, p x4 = 1). The counts match the legend except Anxa1 rewards (22 against 23).
ed9k-snc-variance-explained¶
The striatal PCs explain most of the variance of the SNc cross-correlations, 77.6% (PC1 52.4%, PC2 25.2%).
snc_loco = df[df.snc & df.loco & df.cc470.map(lambda c: isinstance(c, np.ndarray))]
print("SNc recordings in the locomotion selection:", snc_loco.groupby("subtype").agg(n=("path", "size"), mice=("mouse", "nunique")).reindex(ORDER).to_dict("index"))
print("Paper: Vglut2 28 (11 mice), Calb1 6 (3), Anxa1 34 (8), Aldh1a1 42 (13), DAT 31 (8)")
print("Paper: Extended Data Fig. 9k legend 77.6% (PC1 52.4%, PC2 25.2%); Methods PC1 53.2%, PC2 24.3%")
Xs = np.stack(snc_loco[snc_loco.subtype.isin(["Vglut2", "Calb1", "Anxa1"])].cc470.values)
var = np.diag(coeff.T @ (Xs.T @ Xs) @ coeff) / (len(Xs) - 1)
ex = var / var.sum() * 100
print(f"Striatal PCs on {len(Xs)} SNc recordings: PC1 {ex[0]:.2f}%, PC2 {ex[1]:.2f}%, together {ex[:2].sum():.2f}%")
own = az.pca_uncentered(Xs)[1]
print(f"Reference, the SNc recordings' own PCA: PC1 {own[0]:.1f}%, PC2 {own[1]:.1f}%")
tot = np.trace(Xs.T @ Xs) / (len(Xs) - 1)
print(f"Striatal PCs used: {coeff.shape[1]} (one per striatal recording); they capture {var.sum() / tot * 100:.2f}% of trace(Xs'Xs)")
print(f"With the trace as denominator: PC1 {var[0] / tot * 100:.2f}%, PC2 {var[1] / tot * 100:.2f}%, together {var[:2].sum() / tot * 100:.2f}%")
SNc recordings in the locomotion selection: {'Vglut2': {'n': 26, 'mice': 11}, 'Calb1': {'n': 6, 'mice': 3}, 'Anxa1': {'n': 30, 'mice': 7}, 'DAT': {'n': 28, 'mice': 8}, 'Aldh1a1': {'n': 40, 'mice': 13}}
Paper: Vglut2 28 (11 mice), Calb1 6 (3), Anxa1 34 (8), Aldh1a1 42 (13), DAT 31 (8)
Paper: Extended Data Fig. 9k legend 77.6% (PC1 52.4%, PC2 25.2%); Methods PC1 53.2%, PC2 24.3%
Striatal PCs on 62 SNc recordings: PC1 48.09%, PC2 26.96%, together 75.05%
Reference, the SNc recordings' own PCA: PC1 58.3%, PC2 22.8%
Striatal PCs used: 160 (one per striatal recording); they capture 100.00% of trace(Xs'Xs)
With the trace as denominator: PC1 48.09%, PC2 26.96%, together 75.05%
Result: reproduced. The striatal PCs explain 75.05% of the variance of 62 SNc cross-correlations, PC1 48.09% and PC2 26.96%. These are within 10% of the legend's reading (77.6, 52.4, 25.2); under the Methods' reading (53.2 and 24.3) PC2 is outside. The SNc selection has fewer recordings than the paper for Vglut2 and Anxa1 (26 and 30 against 28 and 34). The 160 striatal components capture 100.00% of the SNc data's total variance, so the code's denominator (the sum over the components) equals the total and the percentages are not inflated.
ed1h-aldh1a1-depth¶
Aldh1a1 axons are heterogeneous; dorsal recordings are acceleration-correlated (Type 1) and more ventral ones deceleration-correlated (Type 2).
al = loco[loco.subtype == "Aldh1a1"]
r_, p_ = stats.spearmanr(al.depth, al.ccint)
print(f"Aldh1a1: n = {len(al)} from {mice(al)} mice; Spearman rho(depth, integral 0-1 s) = {r_:.3f}, p = {p_:.2g}")
for lo, hi in ((0, 2.5), (2.5, 3.0), (3.0, 9)):
x = al[(al.depth >= lo) & (al.depth < hi)]
print(f" depth {lo} to {hi} mm: n = {len(x)}, mean integral {x.ccint.mean():.4f}, fraction positive {(x.ccint > 0).mean():.2f}")
for s in ("DAT", "Anxa1", "Calb1"):
x = loco[loco.subtype == s]
print(f"Reference {s}: rho = {stats.spearmanr(x.depth, x.ccint)[0]:.3f}, p = {stats.spearmanr(x.depth, x.ccint)[1]:.2g}")
print("Within-mouse (Aldh1a1 mice with 4 or more recordings), rho per mouse:",
{m: round(float(stats.spearmanr(x.depth, x.ccint)[0]), 2) for m, x in al.groupby("mouse") if len(x) >= 4})
print("\nAldh1a1 fiber_depth_in_mm:", al.depth.describe().round(2).to_dict())
cd = df.dropna(subset=["dv"])
print(f"All fibers with coordinates: n = {len(cd)}; depth minus DV coordinate: median {(cd.depth - cd.dv).median():.2f} mm, "
f"range {(cd.depth - cd.dv).min():.2f} to {(cd.depth - cd.dv).max():.2f}; Spearman rho(depth, DV) = {stats.spearmanr(cd.depth, cd.dv)[0]:.2f}")
print("Aldh1a1 with coordinates:", al.dv.notna().sum(), "of", len(al), "; DV range", al.dv.min(), "to", al.dv.max())
print("Variant bands (deep band only):")
for lo in (2.75, 2.9, 3.0, 3.2):
x = al[al.depth >= lo]
print(f" depth >= {lo} mm: n = {len(x)}, mean integral {x.ccint.mean():.4f}, fraction negative {(x.ccint < 0).mean():.2f}")
Aldh1a1: n = 85 from 15 mice; Spearman rho(depth, integral 0-1 s) = -0.569, p = 1.3e-08
depth 0 to 2.5 mm: n = 60, mean integral 0.0242, fraction positive 0.85
depth 2.5 to 3.0 mm: n = 20, mean integral -0.0014, fraction positive 0.35
depth 3.0 to 9 mm: n = 5, mean integral 0.0001, fraction positive 0.60
Reference DAT: rho = -0.331, p = 0.0015
Reference Anxa1: rho = -0.285, p = 0.0094
Reference Calb1: rho = -0.602, p = 0.00071
Within-mouse (Aldh1a1 mice with 4 or more recordings), rho per mouse: {'Raldh-1572': -0.78, 'Raldh-1573': -0.95, 'Raldh-1697': -0.93, 'Raldh-1698': -0.95, 'Raldh-1846': -0.7, 'Raldh-1849': 0.1, 'Raldh-2076': -0.85, 'Raldh-2220': -0.52, 'Raldh-2276': -0.64}
Aldh1a1 fiber_depth_in_mm: {'count': 85.0, 'mean': 2.26, 'std': 0.54, 'min': 1.6, '25%': 1.9, '50%': 2.15, '75%': 2.5, 'max': 3.9}
All fibers with coordinates: n = 856; depth minus DV coordinate: median -0.83 mm, range -1.45 to 1.03; Spearman rho(depth, DV) = 0.92
Aldh1a1 with coordinates: 80 of 85 ; DV range 2.072096551941404 to 4.800964255147742
Variant bands (deep band only):
depth >= 2.75 mm: n = 19, mean integral -0.0043, fraction negative 0.63
depth >= 2.9 mm: n = 18, mean integral -0.0036, fraction negative 0.61
depth >= 3.0 mm: n = 5, mean integral 0.0001, fraction negative 0.40
depth >= 3.2 mm: n = 4, mean integral -0.0005, fraction negative 0.50
Result: partially reproduced. Aldh1a1 cross-correlation integrals fall with depth (Spearman rho = -0.569, p = 1.3e-08) and are positive above 2.5 mm (mean 0.024, 85% positive), as stated; but the 5 recordings below 3 mm average 0.0001, not below 0, so the second part of the criterion fails. The depth trend holds within 8 of 9 Aldh1a1 mice with 4 or more recordings. Aldh1a1 depths run from 1.6 to 3.9 mm; across fibers with coordinates the depth follows the DV coordinate (Spearman rho = 0.92) and lies a median 0.83 mm above it, as a depth from the brain surface would. With the deep band starting at 2.75 or 2.9 mm (19 and 18 recordings) the mean integral is negative (-0.0043 and -0.0036), so the failing part rests on the 5 recordings deeper than 3 mm.
ed1i-aldh1a1-reward-locomotion¶
In Aldh1a1 recordings, larger reward responses go with deceleration-correlated signaling, and small or negative reward responses with acceleration-correlated signaling.
ai = loco[(loco.subtype == "Aldh1a1") & loco.has_events & (loco.n_Reward >= 3)]
r_, p_ = stats.spearmanr(ai.resp_methods_Reward, ai.ccint)
rc, pc = stats.spearmanr(ai.resp_code_Reward, ai.ccint)
print(f"Aldh1a1 recordings with both: n = {len(ai)} from {mice(ai)} mice; Spearman rho(reward response, integral) = {r_:.3f}, p = {p_:.2g}; code windows rho = {rc:.3f}, p = {pc:.2g}")
print("Reward response by sign of the integral:", (ai.groupby(ai.ccint > 0).resp_methods_Reward.mean() * 100).round(2).to_dict())
Aldh1a1 recordings with both: n = 24 from 8 mice; Spearman rho(reward response, integral) = -0.509, p = 0.011; code windows rho = -0.520, p = 0.0092
Reward response by sign of the integral: {False: 9.59, True: 1.9}
Result: reproduced. In 24 Aldh1a1 recordings, larger reward responses go with more negative (deceleration-correlated) cross-correlation integrals (Spearman rho = -0.509, p = 0.011); the mean reward response is 9.59 where the integral is negative and 1.90 where it is positive.
fig3a-overlap-locomotion¶
Where Calb1 and Anxa1 axons overlap in striatum, Calb1 recordings remain deceleration-correlated and Anxa1 recordings acceleration-correlated.
SCALE, SEC = 35.6667, 204.4
C1 = ((91 - 37) / SCALE, 101 / SCALE + 0.67, 0.86) # FiberLocations.m offsets, slice 1
C2 = ((306 - SEC - 36) / SCALE, 101 / SCALE + 0.68, 0.50) # slice 2
print("Circle centres (ML, DV, AP):", np.round(C1, 3), np.round(C2, 3), "; radius 0.5 mm")
def in_overlap(d):
a = np.isclose(d.ap, C1[2]) & ((d.ml - C1[0]) ** 2 + (d.dv - C1[1]) ** 2 <= 0.25)
b = np.isclose(d.ap, C2[2]) & ((d.ml - C2[0]) ** 2 + (d.dv - C2[1]) ** 2 <= 0.25)
return a | b
ov = loco[in_overlap(loco) & loco.subtype.isin(["Calb1", "Anxa1"])]
for s in ("Calb1", "Anxa1"):
x = ov[ov.subtype == s]
print(f"{s}: n = {len(x)} from {mice(x)} mice, mean integral {x.ccint.mean():.4f}, values {np.round(np.sort(x.ccint.values), 3).tolist()}")
c_, a_ = ov[ov.subtype == "Calb1"].ccint, ov[ov.subtype == "Anxa1"].ccint
print(f"Rank-sum Calb1 vs Anxa1: p = {p_rs(c_, a_):.2g}; paper n Calb1 5 (4 mice), Anxa1 16 (5 mice)")
Circle centres (ML, DV, AP): [1.514 3.502 0.86 ] [1.839 3.512 0.5 ] ; radius 0.5 mm Calb1: n = 17 from 4 mice, mean integral -0.0269, values [-0.061, -0.055, -0.047, -0.039, -0.037, -0.034, -0.03, -0.03, -0.029, -0.027, -0.022, -0.016, -0.012, -0.01, -0.006, -0.003, 0.003] Anxa1: n = 20 from 6 mice, mean integral 0.0092, values [-0.032, -0.024, -0.024, -0.013, -0.013, -0.01, -0.007, -0.003, -0.002, 0.007, 0.011, 0.015, 0.018, 0.018, 0.019, 0.023, 0.041, 0.049, 0.05, 0.063] Rank-sum Calb1 vs Anxa1: p = 0.00012; paper n Calb1 5 (4 mice), Anxa1 16 (5 mice)
Result: reproduced. In the overlap region the Calb1 integrals average -0.027 (16 of 17 negative) and the Anxa1 integrals 0.009, and the two differ (rank-sum p = 0.00012). The region holds 17 Calb1 recordings from 4 mice and 20 Anxa1 from 6, against the paper's 5 from 4 mice and 16 from 5 mice: as many Calb1 mice and one more Anxa1 mouse, and more recordings, as in recording-counts-locomotion. Whether the mice are the paper's cannot be told from the counts.
fig4l-overlap-stimuli¶
Where Calb1 and Anxa1 axons overlap, Calb1 recordings respond to rewards and Anxa1 recordings do not, and Anxa1 recordings decrease to air puffs.
os_ = stim[in_overlap(stim) & stim.subtype.isin(["Calb1", "Anxa1"])]
for s in ("Calb1", "Anxa1"):
x = os_[os_.subtype == s]
xr, xp = x[x.n_Reward >= 3], x[x.n_AirPuff >= 3]
print(f"{s}: reward n = {len(xr)} ({mice(xr)} mice) mean {(xr.resp_methods_Reward * 100).mean():.2f}; "
f"air puff n = {len(xp)} ({mice(xp)} mice) mean {(xp.resp_methods_AirPuff * 100).mean():.2f}")
cr = os_[(os_.subtype == "Calb1") & (os_.n_Reward >= 3)].resp_methods_Reward
ar = os_[(os_.subtype == "Anxa1") & (os_.n_Reward >= 3)].resp_methods_Reward
print(f"Rank-sum reward Calb1 vs Anxa1: p = {p_rs(cr, ar):.2g}; paper n Calb1 9 (4 mice), Anxa1 13 rewards and 17 air puffs (5 mice)")
Calb1: reward n = 9 (4 mice) mean 14.35; air puff n = 9 (4 mice) mean 6.58 Anxa1: reward n = 13 (5 mice) mean 0.94; air puff n = 17 (5 mice) mean -2.06 Rank-sum reward Calb1 vs Anxa1: p = 0.00052; paper n Calb1 9 (4 mice), Anxa1 13 rewards and 17 air puffs (5 mice)
Result: reproduced. In the overlap region Calb1 recordings respond to rewards (mean 14.35) and Anxa1 recordings barely (0.94), the two differ (rank-sum p = 0.00052), and Anxa1 recordings decrease to air puffs (-2.06). The counts equal the paper's: Calb1 9 from 4 mice, Anxa1 13 for rewards and 17 for air puffs from 5 mice, which also confirms the region's conversion to mm.
snrna-nuclei-qc¶
Single-nucleus RNA-seq of dopamine neurons from five DAT-Cre, Sun1-GFP mice yields 12,065 nuclei after quality control, with a mean of 3,435 UMIs and 1,683 genes per nucleus.
These data are not in the dandiset; they are the paper's GEO deposit GSE222558, read here from its two 10x count matrices.
t0 = time.time()
nuc = az.load_geo_counts()
print(f"GEO GSE222558: {len(nuc)} barcodes in two samples, read in {time.time() - t0:.0f} s; by sample", nuc["sample"].value_counts().to_dict())
qc = nuc[nuc.n_features >= 500]
print(f"Nuclei with 500 or more genes: {len(qc)} (by sample {qc['sample'].value_counts().to_dict()}); "
f"mean UMI {qc.n_umi.mean():.0f}, mean genes {qc.n_features.mean():.0f}")
for thr in (200, 400, 600, 800):
q = nuc[nuc.n_features >= thr]
print(f" threshold {thr}: {len(q)} nuclei, mean UMI {q.n_umi.mean():.0f}, mean genes {q.n_features.mean():.0f}")
GEO GSE222558: 12151 barcodes in two samples, read in 4 s; by sample {'M2': 6594, 'F2': 5557}
Nuclei with 500 or more genes: 12065 (by sample {'M2': 6556, 'F2': 5509}); mean UMI 3436, mean genes 1684
threshold 200: 12151 nuclei, mean UMI 3416, mean genes 1675
threshold 400: 12150 nuclei, mean UMI 3416, mean genes 1675
threshold 600: 11967 nuclei, mean UMI 3458, mean genes 1693
threshold 800: 11724 nuclei, mean UMI 3509, mean genes 1714
Result: reproduced. Nuclei with 500 or more genes number 12,065, with a mean of 3,436 UMIs and 1,684 genes, against 12,065, 3,435 and 1,683.
snrna-anxa1-within-aldh1a1¶
Anxa1 expression is limited to a subset of Aldh1a1-expressing dopamine neurons.
an, al_ = qc.Anxa1 > 0, qc.Aldh1a1 > 0
tab = pd.crosstab(an, al_, rownames=["Anxa1 detected"], colnames=["Aldh1a1 detected"])
print(tab.to_string())
odds, p = stats.fisher_exact(tab.values)
print(f"Fisher exact: odds ratio {odds:.2f}, p = {p:.2g}")
print(f"Aldh1a1 detected in {al_[an].mean():.3f} of Anxa1-detected nuclei (all nuclei {al_.mean():.3f})")
print(f"Anxa1 detected in {an[al_].mean():.3f} of Aldh1a1-detected nuclei (all nuclei {an.mean():.3f})")
print("Dopaminergic markers detected (Th, Slc6a3):", round((qc.Th > 0).mean(), 3), round((qc.Slc6a3 > 0).mean(), 3))
for thr in (2, 3):
a2, b2 = qc.Anxa1 >= thr, qc.Aldh1a1 > 0
print(f"Anxa1 with {thr}+ UMIs: n = {int(a2.sum())}, Aldh1a1 detected in {b2[a2].mean():.3f} of them")
Aldh1a1 detected False True Anxa1 detected False 8215 3334 True 119 397 Fisher exact: odds ratio 8.22, p = 1.6e-107 Aldh1a1 detected in 0.769 of Anxa1-detected nuclei (all nuclei 0.309) Anxa1 detected in 0.106 of Aldh1a1-detected nuclei (all nuclei 0.043) Dopaminergic markers detected (Th, Slc6a3): 0.654 0.735 Anxa1 with 2+ UMIs: n = 140, Aldh1a1 detected in 0.771 of them Anxa1 with 3+ UMIs: n = 53, Aldh1a1 detected in 0.811 of them
Result: reproduced. Aldh1a1 is detected in 76.9% of Anxa1-detected nuclei (30.9% of all), the two are associated (odds ratio 8.22, p = 1.6e-107), and Anxa1 is detected in 10.6% of Aldh1a1-detected nuclei. With 2 or 3 Anxa1 UMIs required, 77.1% and 81.1% have Aldh1a1.
Summary¶
| Claim | Outcome |
|---|---|
| recording-counts-locomotion | not reproduced |
| recording-counts-stimuli | reproduced |
| fig2f-xcorr-direction | reproduced |
| fig2i-trough-lag | partially reproduced |
| fig2jk-pca-variance | reproduced |
| fig2l-pc-angles | reproduced |
| fig2g-event-counts | reproduced |
| ed7f-prob-transient-after-event | reproduced |
| ed6f-dec-peak-timing | partially reproduced |
| fig5ab-transient-scaling | partially reproduced |
| fig4g-reward-response | partially reproduced |
| fig4g-anxa1-no-reward | reproduced |
| fig4g-airpuff-response | partially reproduced |
| fig4h-reward-vs-airpuff | reproduced |
| fig4i-reward-size | partially reproduced |
| ed8f-rest-reward-same | reproduced |
| fig6c-kmeans | reproduced |
| fig7f-soma-axon-correlation | reproduced |
| ed9c-snc-stimulus-responses | reproduced |
| ed9k-snc-variance-explained | reproduced |
| ed1h-aldh1a1-depth | partially reproduced |
| ed1i-aldh1a1-reward-locomotion | reproduced |
| fig3a-overlap-locomotion | reproduced |
| fig4l-overlap-stimuli | reproduced |
| snrna-nuclei-qc (GEO) | reproduced |
| snrna-anxa1-within-aldh1a1 (GEO) | reproduced |
| snrna-clusters (GEO) | blocked: not_runnable |
| scrna-meta-integration (other datasets) | blocked: not_runnable |
| histology-anxa1-cre (Zenodo images) | blocked: other |
A failure to reproduce here may be a failure of this attempt, not of the paper: see the reasoning and the differences from the paper's method.
How the analysis differs from the paper's. It is a Python port of the paper's MATLAB code. Response windows and peak times follow the Methods, which differ from the code by 0.05 s and by one or two bins; the code's choices are shown as variants. With them the reward means equal the paper's to the decimal and the Calb1 trough lag (0.160 against 0.17 s) is within tolerance. The locomotion selection, applied as the Methods and code state, keeps more recordings than the paper reports, so the locomotion and clustering claims rest on more recordings than the paper's. Running periods after rewards the mouse did not consume cannot be excluded, as those rewards are not in the Events table; for the event-size analysis the stored onsets, which carry that exclusion, are used. P-values come from scipy's Wilcoxon and Mann-Whitney tests, whose exact and approximate methods may differ slightly from MATLAB's. k-means uses scikit-learn and is summarized over 100 seeds.
What was found about the dandiset. It holds the paper's photometry: the stimulus analyses
recover the paper's recording and mouse counts exactly, and with the code's windows its reward
means to the decimal and its air puff means within tolerance except DAT (2.46 against 5.3);
the DAT reward size difference is 1.83 against 1.9. The selection flags of the paper's table (running time, artifacts) are not
in the files but can be recomputed, and the recomputed values match the stored 405 nm
cross-correlations. Which recordings entered the locomotion figures cannot be recovered: the stated
criteria select more. The coordinates column is stored as (ML, AP, DV), not the (AP, ML, DV) its
description says. cross_correlation_with_acceleration holds the 405 nm control's peak
cross-correlation, not that of dF/F. The dF/F series are normalized to their largest transient,
which the files do not state. Consumed rewards only are in the Events table. Three subject ids
ending in L are mice that also have ids without it, so the archive's 105 subjects are 102 mice.
No recording is deposited twice. normalized_fluorescence holds two values per fiber in 123
files. The published version's link to the paper is malformed and its asset summary names no
photometry technique or variable. The paper itself gives the Calb1 and Vglut2 reward P-values the
other way round in the Results and the Figure 4g legend, 10 Anxa1 mice in the Figure 2f legend and
9 in the Methods, and different PC percentages for Extended Data Figure 9k in the legend and the
Methods (printed in the claim sections).
print(f"Notebook ran in {(time.time() - started) / 60:.1f} min")
Notebook ran in 1.5 min