Gillon, Lecoq et al. (2023) and dandisets 000336 and 000037¶
Paper: Gillon CJ, Lecoq JA, et al. "Responses of pyramidal cell somata and apical dendrites in mouse
visual cortex over multiple days." Scientific Data 10, 287 (2023), doi:10.1038/s41597-023-02214-y.
The claims file (claims/000336/10.1038-s41597-023-02214-y.yaml) was written from the version of
record as served by Europe PMC (PMC10192216), with its tables and figure legends; Supplementary Table 1
was not available.
Target dataset: DANDI:000336, version 0.260325.1802, the only published version (the other is the
draft). It is the OpenScope "Dendritic Coupling" project. It names the paper among its related
resources, but the paper's data are on DANDI:000037 (OpenScope "Credit Assignment"), the dandiset the
paper's Data Records cite. The first claim tests whether 000336 holds any of the paper's sessions; all
other claims are marked data: elsewhere in the claims file and are tested on DANDI:000037.
Data used: DANDI:000037, version 0.230426.0054 (doi:10.48324/DANDI.000037/0.230426.0054). The paper
cites the dandiset without a version. The claims file names this version for every claim, and it is the
version published shortly before the paper (April 2023, paper published May 2023). The later version
0.240209.1623 adds only a CHANGES file (shown below), so the choice does not change the data read.
Only the basic files (_behavior+ophys.nwb) are read, streamed from the archive.
Outside input: the session table of the paper's code, mouse_df.csv of
github.com/colleenjg/OpenScope_CA_Analysis at commit 046e7fef5f452cdc0775fe54040b0f13aeae1e57, for the
session numbers (1 to 6 per mouse) and the QC status, which the files do not carry. The paper's Tables 1
and 2, copied from its text into the helper module, are used only as reference values.
| Claim | Tested here |
|---|---|
| dataset-holds-paper-sessions | yes, on 000336 |
| table1-subjects | yes |
| table2-sessions | yes |
| session-schedule | yes |
| imaging-planes | yes |
| roi-counts-per-session | yes |
| tracked-rois-per-plane | yes |
| fig8-example-matches | yes |
| fig8-ordering-variation | yes, in roi_matching_orderings.ipynb |
| stimulus-durations | yes |
| inconsistent-event-fraction | yes |
| fig5-stimulus-responses | yes |
| fig6-behaviour | yes |
| fig7a-snr-high-stable | yes |
| fig7b-signal-stable | yes |
| fig7c-correlations-similar | yes |
| high-correlation-pairs | yes |
| fig9-onset-offset-transients | yes |
In [1]:
import os
os.environ.setdefault("OMP_NUM_THREADS", "4")
import time
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.stats
import h5py, remfile, pynwb
import credassign as ca
pd.set_option("display.width", 200)
pd.set_option("display.max_columns", 40)
pd.set_option("display.max_rows", 200)
T0 = time.time()
The target, DANDI:000336¶
The archive's metadata, the file list, and one file opened and listed in full, including every column of every table.
In [2]:
m336 = ca.version_metadata(ca.TARGET_DANDISET, ca.TARGET_VERSION)
print(m336["name"])
print("published versions:", [v["version"] for v in ca.get_json(f"{ca.API}/dandisets/{ca.TARGET_DANDISET}/versions/")["results"]])
print("subjects (archive summary):", m336["assetsSummary"]["numberOfSubjects"], " files:", m336["assetsSummary"]["numberOfFiles"],
" bytes:", m336["assetsSummary"]["numberOfBytes"])
for r in m336["relatedResource"]:
print("related:", r["relation"], r.get("identifier", r["url"]), "|", r["name"])
a336 = ca.list_assets(ca.TARGET_DANDISET, ca.TARGET_VERSION)
a336["subject"] = a336.path.str.extract(r"sub-(\d+)/")[0]
a336["session"] = a336.path.str.extract(r"_ses-(\d+)")[0]
a336["kind"] = a336.path.str.extract(r"_([a-z+]+)\.nwb$")[0]
print(len(a336), "assets;", a336.subject.nunique(), "subjects;", a336.session.nunique(), "session IDs")
print(a336.groupby("kind").size().to_string())
print("subjects:", sorted(a336.subject.unique()))
Allen Institute Openscope - Dendritic Coupling project
published versions: ['draft', '0.260325.1802'] subjects (archive summary): 16 files: 2356 bytes: 42229353097330 related: dcite:IsDocumentedBy https://alleninstitute.github.io/openscope_databook/projects/dendritic_coupling.html | OpenScope Databook - Dendritic Coupling project documentation related: dcite:IsDocumentedBy https://github.com/AllenInstitute/openscope_dendritic_coupling | Public code repository for Dendritic Coupling project related: dcite:Describes 10.1016/j.tins.2012.11.006 | A cellular mechanism for cortical associations: an organizing principle for the cerebral cortex. related: dcite:Describes 1810.11393 | Dendritic cortical microcircuits approximate the backpropagation algorithm related: dcite:Describes 10.1016/j.neuron.2018.10.003 | Predictive Processing: A Canonical Cortical Computation related: dcite:Describes 10.1038/s41597-023-02214-y | Responses of pyramidal cell somata and apical dendrites in mouse visual cortex over multiple days
2356 assets; 16 subjects; 104 session IDs kind image 2 image+ophys 362 ophys 1992 subjects: ['621602', '621603', '631188', '634402', '635489', '641331', '644972', '647715', '655042', '655580', '656228', '660986', '661749', '662470', '690226', '698808']
In [3]:
# One file of 000336, the smallest processed ophys file, listed in full.
small = a336[a336.kind == "ophys"].sort_values("size").iloc[0]
print(small.path, small["size"])
with pynwb.NWBHDF5IO(file=ca.open_h5(small.asset_id), load_namespaces=True) as io:
nwb = io.read()
print("session_start_time:", nwb.session_start_time, "| session_id:", nwb.session_id, "| identifier:", nwb.identifier)
s = nwb.subject
print("subject:", s.subject_id, s.genotype, s.sex, s.age, "|", s.description)
for k, p in nwb.imaging_planes.items():
print("imaging plane", k, "|", p.description, "|", p.location, "| rate", p.imaging_rate)
print("acquisition:", {k: type(v).__name__ for k, v in nwb.acquisition.items()})
for name, mod in nwb.processing.items():
for dn, d in mod.data_interfaces.items():
print("processing", name, dn, type(d).__name__)
if hasattr(d, "colnames"):
print(" columns:", list(d.colnames))
for attr in ("roi_response_series", "time_series", "plane_segmentations"):
for sn, ser in (getattr(d, attr, None) or {}).items():
print(" ", sn, getattr(getattr(ser, "data", None), "shape", None), list(getattr(ser, "colnames", [])))
for k, t in nwb.intervals.items():
print("intervals", k, len(t), "rows; columns:", list(t.colnames))
sub-662470/sub-662470_ses-1253666600-acq-1253834948_ophys.nwb 89216248
session_start_time: 2023-03-13 15:58:48.095000+00:00 | session_id: None | identifier: 1253834948
subject: 662470 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt F P84.0D | external: 662470 donor_id: 1244572630 specimen_id: [(1244572640,)]
imaging plane imaging_plane | Failed Cell Segmentation | Failed Cell Segmentation | rate 10.71
imaging plane imaging_plane_1 | (512, 512) field of view in VISl at depth 16 um | VISl | rate 10.0
acquisition: {'EyeTracking': 'EllipseEyeTracking', 'v_in': 'TimeSeries', 'v_sig': 'TimeSeries'}
processing ophys images Images
processing ophys ophys_motion_correction_x TimeSeries
processing ophys ophys_motion_correction_y TimeSeries
processing running dx TimeSeries
processing running speed TimeSeries
processing running speed_unfiltered TimeSeries
processing stimulus timestamps TimeSeries
processing stimulus_ophys timestamps TimeSeries
intervals fixed_gabors_presentations 3000 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'stimulus_block', 'units', 'stimulus_index', 'OriSurp', 'PosSizesAll', 'colors', 'contrs', 'depths', 'elementMask', 'elementTex', 'fieldDepth', 'fieldPos', 'fieldShape', 'fieldSize', 'nElements', 'opacities', 'oris', 'phases', 'rgbs', 'sfs', 'sizes', 'xys', 'tags', 'timeseries']
intervals gratings_presentations 624 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'stimulus_block', 'color', 'contrast', 'mask', 'opacity', 'ori', 'size', 'units', 'stimulus_index', 'spatial_frequency', 'phase', 'tags', 'timeseries']
intervals movie_flower_fwd_presentations 16200 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'stimulus_block', 'frame', 'color', 'contrast', 'opacity', 'ori', 'size', 'units', 'stimulus_index', 'tags', 'timeseries']
intervals movie_touch_of_evil_fwd_presentations 16200 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'stimulus_block', 'frame', 'color', 'contrast', 'opacity', 'ori', 'size', 'units', 'stimulus_index', 'tags', 'timeseries']
intervals movie_worms_fwd_presentations 16200 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'stimulus_block', 'frame', 'color', 'contrast', 'opacity', 'ori', 'size', 'units', 'stimulus_index', 'tags', 'timeseries']
intervals rotate_gabors_presentations 3000 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'stimulus_block', 'units', 'stimulus_index', 'OriSurp', 'PosSizesAll', 'colors', 'contrs', 'depths', 'elementMask', 'elementTex', 'fieldDepth', 'fieldPos', 'fieldShape', 'fieldSize', 'nElements', 'opacities', 'oris', 'phases', 'rgbs', 'sfs', 'sizes', 'xys', 'tags', 'timeseries']
intervals spontaneous_presentations 183 rows; columns: ['start_time', 'stop_time', 'stimulus_name', 'tags', 'timeseries']
The paper's data, DANDI:000037¶
The files list, the difference between its published versions, and one basic file listed in full.
In [4]:
v37 = [v["version"] for v in ca.get_json(f"{ca.API}/dandisets/{ca.DANDISET}/versions/")["results"]]
print("versions of 000037:", v37)
a37 = ca.list_assets()
later = ca.list_assets(ca.DANDISET, "0.240209.1623")
both = a37.merge(later, on="path", how="outer", suffixes=("", "_later"), indicator=True)
print("only in the pinned version:", both.loc[both._merge == "left_only", "path"].tolist())
print("only in 0.240209.1623:", both.loc[both._merge == "right_only", "path"].tolist())
print("paths in both with a different size:", int((both.loc[both._merge == "both", "size"] != both.loc[both._merge == "both", "size_later"]).sum()))
a37["subject"] = a37.path.str.extract(r"sub-(\d+)/")[0]
a37["session"] = a37.path.str.extract(r"_ses-(\d+)_")[0]
a37["kind"] = a37.path.str.extract(r"_ses-\d+_(.*)\.nwb$")[0]
print(a37.groupby("kind").agg(n=("path", "size"), min_gb=("size", lambda x: x.min() / 1e9), max_gb=("size", lambda x: x.max() / 1e9)))
basic = a37[a37.kind == "behavior+ophys"].sort_values(["subject", "session"]).reset_index(drop=True)
print(len(basic), "basic files,", basic.session.nunique(), "sessions,", basic.subject.nunique(), "subjects")
print("each session has all three versions:", (a37.groupby("session").size() == 3).all())
versions of 000037: ['draft', '0.230426.0054', '0.240209.1623']
only in the pinned version: []
only in 0.240209.1623: ['CHANGES']
paths in both with a different size: 0
n min_gb max_gb
kind
behavior+image+ophys 50 1.589581 3.188968
behavior+ophys 50 0.136659 1.703330
obj-raw_behavior+image+ophys 50 42.864406 49.504858
50 basic files, 50 sessions, 13 subjects
each session has all three versions: True
In [5]:
# One basic file of 000037 listed in full (subject 418779, session 777496949).
row = basic[basic.session == "777496949"].iloc[0]
with pynwb.NWBHDF5IO(file=ca.open_h5(row.asset_id), load_namespaces=True) as io:
nwb = io.read()
print("session_start_time:", nwb.session_start_time, "| session_id:", nwb.session_id)
print("subject:", nwb.subject.subject_id, nwb.subject.genotype, nwb.subject.sex, nwb.subject.age)
for k, p in nwb.imaging_planes.items():
print("imaging plane", k, "| location", p.location, "| rate", p.imaging_rate, "| indicator", p.indicator)
for name, mod in nwb.processing.items():
for dn, d in mod.data_interfaces.items():
for attr in ("roi_response_series", "time_series", "plane_segmentations", "images"):
for sn, ser in (getattr(d, attr, None) or {}).items():
print(name, dn, sn, getattr(ser, "data", None).shape if hasattr(ser, "data") else "",
getattr(ser, "unit", ""), list(getattr(ser, "colnames", [])), getattr(ser, "description", ""))
for col in nwb.trials.columns:
print("trials column:", col.name, "|", col.description)
tr = nwb.trials.to_dataframe()
print(tr.groupby(["stimulus_type", "stimulus_template_name"]).size())
print("acquisition:", list(nwb.acquisition), "| stimulus:", list(nwb.stimulus))
session_start_time: 2018-11-13 01:51:48.328000-08:00 | session_id: 777496949
subject: 418779 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 F P84D
imaging plane ImagingPlane | location {'area': 'VISp', 'depth': '375'} | rate 30.0 | indicator GCaMP6f
behavior BehavioralTimeSeries running_velocity (251999,) cm/s [] Velocity of the mouse on a rotating disc.
behavior PupilTracking pupil_diameter (126727,) mm [] Diameter of the mouse pupil (right) facing the stimulus presentation screen.
behavior PupilTracking pupil_position_x (126727,) mm [] Position in x of the center of the mouse pupil (right) facing the stimulus presentation screen.
behavior PupilTracking pupil_position_y (126727,) mm [] Position in y of the center of the mouse pupil (right) facing the stimulus presentation screen.
ophys DfOverF RoiResponseSeries (126727, 15) Normalized fluorescence (A.U.) [] ROI traces
ophys ImageSegmentation PlaneSegmentation ['image_mask', 'tracking_id'] Segmentation for imaging plane (de Vries et al., 2019, Nat Neurosci)
ophys PlaneImages max_projection (512, 512) [] Normalized maximum projection of the motion corrected z-stack, averaged over 8 frames (hei x wid).
trials column: start_time | Start time of epoch, in seconds
trials column: stop_time | Stop time of epoch, in seconds
trials column: stimulus_type | Type of stimulus presented during the trial.
trials column: stimulus_template_name | Name of the stimulus template for the trial.
trials column: unexpected | Whether the stimulus features are expected (0) or unexpected (1).
trials column: gabor_frame | Gabor frame or set, defined by Gabor patch locations and sizes (A, B, C, D or U, G).
trials column: gabor_kappa | Kappa (dispersion) of the von Mises distribution used to sample Gabor patch orientations (degrees).
trials column: gabor_mean_orientation | Mean orientation of the von Mises distribution used to sample Gabor patch orientations (degrees, [0, 360].
trials column: gabor_number | Number of Gabor patches on-screen (before warping mask is applied).
trials column: gabor_locations_x_index | Index for VectorData 'gabor_locations_x'
trials column: gabor_locations_x | Locations along the x dimension (width) of the Gabor patch centers (in pixels), centered on the middle of the screen. All locations are for the unwarped stimulus.
trials column: gabor_locations_y_index | Index for VectorData 'gabor_locations_y'
trials column: gabor_locations_y | Locations along the y dimension (height) of Gabor patch centers (in pixels), centered on the middle of the screen. All locations are for the unwarped stimulus.
trials column: gabor_sizes_index | Index for VectorData 'gabor_sizes'
trials column: gabor_sizes | Gabor patch sizes (in pixels, with patch width defined, as in psychopy, as width across 6 standard deviations).
trials column: gabor_orientations_index | Index for VectorData 'gabor_orientations'
trials column: gabor_orientations | Orientations of the Gabor patches (degrees, [0, 360]).
trials column: main_flow_direction | Main direction of flow for the squares (all moving at 50 visual degrees/second). Flow directions are either rightward (nasal-to-temporal) or leftward (temporal-to-nasal).
trials column: num_frames_stim | Stop frame number for the stimulus frames (excluded).
trials column: num_frames_twop | Stop frame number for the two-photon frames (excluded).
trials column: square_size | Size of the squares (in pixels).
trials column: square_number | Total number of squares present, both in the off and on-screen areas (before warping mask is applied).
trials column: square_proportion_flipped | Proportion of squares moving in the direction opposite to the main flow (all moving at 50 visual degrees/second).
trials column: square_locations_x_index | Index for VectorData 'square_locations_x'
trials column: square_locations_x | Locations along the x dimension (width) of square centers for each frame of the trial (in pixels), (square x stimulus frame), centered on the middle of the screen. Note that squares disappear a bit off-screen, and are (re)initialized a bit off-screen. All locations are for the unwarped stimulus.
trials column: square_locations_y_index | Index for VectorData 'square_locations_y'
trials column: square_locations_y | Locations along the y dimension (height) of square centers for each frame of the trial (in pixels), (square x stimulus frame), centered on the middle of the screen. Note that squares are (re)initialized a bit off-screen area. All locations are for the unwarped stimulus.
trials column: start_frame_stim_template | Start frame number for the trial in the corresponding template.
trials column: start_frame_stim | Start frame number for the stimulus frames.
trials column: start_frame_twop | Start frame number for the two-photon frames.
trials column: stop_frame_stim | Stop frame number for the stimulus frames (excluded).
trials column: stop_frame_twop | Stop frame number for the two-photon frames (excluded).
stimulus_type stimulus_template_name
gabors gabors 6800
grayscreen grayscreen 4
visflow visflow_left 1020
visflow_right 1020
dtype: int64
acquisition: [] | stimulus: []
Data check¶
All 50 basic files are opened for their metadata (subject, session, plane, stream lengths and the trial
table). The ROI traces of the 33 analysed sessions (sessions 1 to 3 that passed QC, by the paper's
code table) are read once, and every measure the claims need is computed as they load (helper,
roi_measures). This takes about eight minutes.
In [6]:
mouse_df = ca.read_mouse_df()
mouse_df["sessid"] = mouse_df.sessid.astype(str)
mouse_df["mouseid"] = mouse_df.mouseid.astype(str)
table1, table2 = ca.paper_table1(), ca.paper_table2()
analysed = mouse_df[(mouse_df.pass_fail == "P") & mouse_df.sess_n.between(1, 3)]
ANALYSED = set(analysed.sessid)
SESS1 = set(analysed.loc[analysed.sess_n == 1, "sessid"])
print(len(ANALYSED), "analysed sessions from", analysed.mouseid.nunique(), "mice;", len(SESS1), "session-1 files")
meta, trials, roi, beh = {}, {}, {}, {}
for r in basic.itertuples():
f = ca.open_h5(r.asset_id)
meta[r.session] = ca.session_metadata(f)
trials[r.session] = ca.trials_table(f)
if r.session in ANALYSED:
roi[r.session] = ca.roi_measures(f, trials[r.session], fig5=r.session in SESS1)
beh[r.session] = ca.behaviour(f)
f.close()
print(f"loaded in {(time.time() - T0) / 60:.1f} min")
33
analysed sessions from 11 mice; 11 session-1 files
loaded in 12.0 min
In [7]:
md_ = pd.DataFrame([{k: v for k, v in m.items() if k not in ("tracking_ids", "twop_timestamps")} for m in meta.values()])
md_["session"] = list(meta)
md_["depth"] = md_.location.str.extract(r"'depth': '(\d+)'")[0].astype(int)
md_["area"] = md_.location.str.extract(r"'area': '(\w+)'")[0]
md_["line"] = np.where(md_.genotype.str.startswith("Cux2"), "L2/3", np.where(md_.genotype.str.startswith("Rbp4"), "L5", "?"))
md_["plane"] = np.where(md_.depth.isin([175, 375]), "soma", "dend")
md_ = md_.merge(mouse_df[["sessid", "sess_n", "pass_fail", "full_timestamp", "date", "plane", "line", "nrois", "nrois_tracked"]]
.rename(columns={"sessid": "session", "plane": "plane_code", "line": "line_code"}), on="session", how="left")
md_["analysed"] = md_.session.isin(ANALYSED)
print("files:", len(md_), "| subjects:", md_.subject_id.nunique(), "| sessions:", md_.session_id.nunique())
print("session_id equals the file name's session:", (md_.session_id == md_.session).all())
print("subject_id equals the file name's subject:", (md_.subject_id.values == basic.subject.values).all())
print(md_.groupby(["line", "plane"]).agg(mice=("subject_id", "nunique"), sessions=("session", "size"), analysed=("analysed", "sum")))
print("plane from depth agrees with the code table:", (md_.plane == md_.plane_code).all(),
"| line from genotype agrees:", ((md_.line == "L2/3") == md_.line_code.str.startswith("L23")).all())
print("areas:", md_.area.unique(), "| imaging rates (metadata):", md_.imaging_rate.unique())
print("2p frame interval, s: median over files", md_.twop_median_dt.median(), "range", md_.twop_median_dt.min(), md_.twop_median_dt.max())
print("2p frames per session:", md_.twop_n.min(), "to", md_.twop_n.max(), "| dF/F shape frames equal timestamps:",
all(m["dff_shape"][0] == m["twop_n"] for m in meta.values()), "| dF/F unit:", md_.dff_unit.unique())
print("running unit:", md_.run_unit.unique(), "| running sample interval (s):", md_.run_median_dt.min(), "to", md_.run_median_dt.max())
print("pupil unit:", md_.pupil_unit.unique(), "| pupil samples equal 2p frames:", (md_.pupil_n == md_.twop_n).all())
print("max projection shapes:", md_.max_projection_shape.astype(str).unique())
print("files whose pupil samples differ from their 2p frames:", md_.loc[md_.pupil_n != md_.twop_n, ["session", "pupil_n", "twop_n"]].to_dict("records"))
print("files without a tracking_id column:", int((~md_.has_tracking_column).sum()), "| their tracked count in Table 2:",
sorted(set(table2.set_index("session_id").loc[md_.loc[~md_.has_tracking_column, "session"], "n_tracked"])))
files: 50 | subjects: 13 | sessions: 50
session_id equals the file name's session: True
subject_id equals the file name's subject: True
mice sessions analysed
line plane
L2/3 dend 3 12 6
soma 3 12 9
L5 dend 4 16 9
soma 3 10 9
plane from depth agrees with the code table: True | line from genotype agrees: True
areas: <StringArray>
['VISp']
Length: 1, dtype: str | imaging rates (metadata): [30.]
2p frame interval, s: median over files 0.033249999999952706 range 0.03316999999969994 0.033249999999952706
2p frames per session: 126721 to 127588 | dF/F shape frames equal timestamps: True | dF/F unit: <StringArray>
['Normalized fluorescence (A.U.)']
Length: 1, dtype: str
running unit: <StringArray>
['cm/s']
Length: 1, dtype: str | running sample interval (s): 0.016679999999723805 to 0.01667999999995118
pupil unit: <StringArray>
['mm']
Length: 1, dtype: str | pupil samples equal 2p frames: False
max projection shapes: <StringArray>
['(512, 512)']
Length: 1, dtype: str
files whose pupil samples differ from their 2p frames: [{'session': '759872185', 'pupil_n': 61134, 'twop_n': 127588}]
files without a tracking_id column: 11 | their tracked count in Table 2: [0]
In [8]:
# The paper's code table (mouse_df.csv), the columns used, for the 50 sessions, beside Table 2.
cols = ["mouseid", "sessid", "sess_n", "date", "pass_fail", "depth", "plane", "line", "nrois", "nrois_tracked"]
mdf = mouse_df[mouse_df.sessid.isin(table2.session_id)][cols].sort_values(["mouseid", "sess_n"])
print(mdf.to_string(index=False))
t2_ = table2.set_index("session_id").loc[mdf.sessid]
print("code table rows for the 50 sessions of Table 2:", len(mdf))
print("agrees with Table 2 on: QC", ((mdf.pass_fail.values == "P") == (t2_.qc.values == "passed")).sum(),
"| date", (mdf.date.astype(str).values == t2_.date.values).sum(), "| depth", (mdf.depth.values == t2_.depth.values).sum(),
"| ROIs", (mdf.nrois.values == t2_.n_rois.values).sum(), "| tracked ROIs", (mdf.nrois_tracked.clip(lower=0).values == t2_.n_tracked.values).sum(), "of 50",
"(the code table gives -1 where Table 2 gives 0)")
# Session numbers: within each mouse, sess_n runs 1, 2, ... in order of Table 2's dates.
order_ok = all((g.sort_values("date").sess_n.values == np.arange(1, len(g) + 1)).all() for _, g in mdf.groupby("mouseid"))
print("sess_n numbers each mouse's sessions 1, 2, ... in date order:", order_ok)
mouseid sessid sess_n date pass_fail depth plane line nrois nrois_tracked 408021 758519303 1 20180926 P 175 soma L23-Cux2 96 59 408021 759189643 2 20180927 P 175 soma L23-Cux2 74 59 408021 759660390 3 20181001 P 175 soma L23-Cux2 107 59 411400 759666166 1 20181001 F 20 dend L5-Rbp4 942 -1 411400 759872185 2 20181002 F 20 dend L5-Rbp4 289 -1 411400 760269100 3 20181003 F 20 dend L5-Rbp4 524 -1 411400 761730740 4 20181009 P 20 dend L5-Rbp4 630 162 411400 762415169 5 20181011 P 20 dend L5-Rbp4 637 162 411400 763646681 6 20181015 P 20 dend L5-Rbp4 597 162 411424 761624763 1 20181009 P 175 soma L23-Cux2 87 55 411424 761944562 2 20181010 P 175 soma L23-Cux2 90 55 411424 762250376 3 20181011 P 175 soma L23-Cux2 80 55 411771 760260459 1 20181003 P 375 soma L5-Rbp4 90 47 411771 760659782 2 20181004 P 375 soma L5-Rbp4 70 47 411771 761269197 3 20181008 P 375 soma L5-Rbp4 79 47 412933 763949859 1 20181016 F 75 dend L23-Cux2 1041 -1 412933 764897534 2 20181017 F 75 dend L23-Cux2 948 -1 412933 765427689 3 20181018 F 75 dend L23-Cux2 836 -1 412933 766755831 4 20181022 P 50 dend L23-Cux2 344 98 412933 767254594 5 20181023 P 50 dend L23-Cux2 168 98 412933 768807532 6 20181025 P 50 dend L23-Cux2 250 98 413663 764704289 1 20181017 P 50 dend L23-Cux2 628 136 413663 765193831 2 20181018 P 50 dend L23-Cux2 348 136 413663 766502238 3 20181022 P 50 dend L23-Cux2 512 136 418779 777496949 1 20181112 P 375 soma L5-Rbp4 15 12 418779 778374308 2 20181113 P 375 soma L5-Rbp4 26 12 418779 779152062 3 20181114 P 375 soma L5-Rbp4 29 12 420011 777914830 1 20181113 P 20 dend L5-Rbp4 205 51 420011 778864809 2 20181114 P 20 dend L5-Rbp4 159 51 420011 779650018 3 20181115 P 20 dend L5-Rbp4 182 51 433414 826187862 1 20190220 P 75 dend L23-Cux2 727 118 433414 826773996 2 20190221 P 75 dend L23-Cux2 300 118 433414 827833392 3 20190222 P 75 dend L23-Cux2 333 118 433448 826338612 1 20190220 P 20 dend L5-Rbp4 1636 112 433448 826819032 2 20190221 P 20 dend L5-Rbp4 445 112 433448 828816509 3 20190225 P 20 dend L5-Rbp4 496 112 433448 829283315 4 20190226 P 20 dend L5-Rbp4 436 -1 433451 823453391 1 20190215 P 20 dend L5-Rbp4 966 353 433451 824434038 2 20190218 P 20 dend L5-Rbp4 1029 353 433451 825180479 3 20190219 P 20 dend L5-Rbp4 986 353 433458 826659257 1 20190221 P 375 soma L5-Rbp4 99 70 433458 827300090 2 20190222 P 375 soma L5-Rbp4 87 70 433458 828475005 3 20190225 P 375 soma L5-Rbp4 97 70 433458 829520904 4 20190227 P 375 soma L5-Rbp4 88 -1 440889 832883243 1 20190306 P 175 soma L23-Cux2 224 147 440889 833704570 2 20190307 P 175 soma L23-Cux2 224 147 440889 834403597 3 20190308 P 175 soma L23-Cux2 210 147 440889 836968429 4 20190314 P 175 soma L23-Cux2 205 -1 440889 837360280 5 20190315 F 175 soma L23-Cux2 217 -1 440889 838633305 6 20190318 F 175 soma L23-Cux2 227 -1 code table rows for the 50 sessions of Table 2: 50 agrees with Table 2 on: QC 50 | date 50 | depth 50 | ROIs 50 | tracked ROIs 50 of 50 (the code table gives -1 where Table 2 gives 0) sess_n numbers each mouse's sessions 1, 2, ... in date order: True
In [9]:
# Ranges of the variables used, analysed sessions.
r_ = pd.DataFrame({s: {"dff_min": v["dff_min"], "dff_max": v["dff_max"], "rois": v["n_rois"], "rois_nan_or_inf": v["n_bad"]}
for s, v in roi.items()}).T
print("dF/F over all ROIs and frames:", r_.dff_min.min(), "to", r_.dff_max.max(), "| sessions with NaN or Inf ROIs:", int((r_.rois_nan_or_inf > 0).sum()))
b_ = pd.DataFrame({s: {"run_min": np.min(v["run"]), "run_max": np.max(v["run"]), "run_dropped": v["run_dropped"],
"pup_min": np.min(v["pupil"]), "pup_max": np.max(v["pupil"]), "pup_nan": v["pupil_nan"],
"pup_dropped": v["pupil_dropped"]} for s, v in beh.items()}).T
print(b_.describe().loc[["min", "50%", "max"]])
# Duplicates: frame and ROI counts, and a fingerprint of the first 1000 frames of every ROI.
dup = md_.groupby(["twop_n", "n_rois"]).session.apply(list)
print("sessions sharing frame and ROI counts:", [v for v in dup if len(v) > 1])
fp = pd.Series({s: v["fingerprint"] for s, v in roi.items()})
print("analysed sessions sharing a trace fingerprint:", int(fp.duplicated(keep=False).sum()))
dF/F over all ROIs and frames: -10.346665467865007 to 30.209230369583594 | sessions with NaN or Inf ROIs: 0
run_min run_max run_dropped pup_min pup_max pup_nan pup_dropped
min -28.410246 11.064477 0.0 0.033900 0.736593 79.0 1.0
50% -10.081872 51.221538 0.0 0.218147 0.877274 504.0 9.0
max -1.899914 83.848043 12.0 0.530078 1.021494 2906.0 23.0
sessions sharing frame and ROI counts: []
analysed sessions sharing a trace fingerprint: 0
Clocks¶
The trial table gives each stimulus row its 2p frames (start_frame_twop), and the 2p, pupil and
running series carry timestamps. The check below compares the table's frame numbers with the 2p
timestamps, and then uses a contrast: the population response at the stimulus onsets of the table,
against the same frames shifted by 5 to 20 s, where no onset is.
In [10]:
lag = []
for s in sorted(ANALYSED):
t, ts = trials[s], meta[s]["twop_timestamps"]
fr = t.start_frame_twop.values.astype(int)
d = ts[fr] - t.start_time.values # 2p timestamp of each row's start frame minus the row's start time
lag.append({"session": s, "median_s": np.median(d), "min_s": d.min(), "max_s": d.max()})
lag = pd.DataFrame(lag)
print("2p timestamp of each row's start frame - row start time, s: medians", lag.median_s.min().round(4), "to", lag.median_s.max().round(4),
"| extremes", lag.min_s.min().round(4), "to", lag.max_s.max().round(4))
print("running timestamps span the stimuli in every session:",
all((meta[s]["run_first"] <= trials[s].start_time.min() + 1) and (meta[s]["run_last"] >= trials[s].stop_time.max() - 1) for s in ANALYSED))
print("pupil timestamps equal the 2p timestamps in every analysed session:",
all(meta[s]["pupil_n"] == meta[s]["twop_n"] and meta[s]["pupil_first"] == meta[s]["twop_first"] for s in ANALYSED))
2p timestamp of each row's start frame - row start time, s: medians -0.0174 to 0.0147 | extremes -0.0434 to 0.0363 running timestamps span the stimuli in every session: True pupil timestamps equal the 2p timestamps in every analysed session: True
In [11]:
# Contrast: per session, the mean over ROIs of (mean dF/F 0-2 s after - 2 s before) at the table's
# onset frames, and at the same frames shifted by -20, -10, -5, +5, +10 and +20 s, in units of its SEM.
z = pd.DataFrame({s: {sh: v["onset_shift"][sh][0] / v["onset_shift"][sh][1] for sh in v["onset_shift"]} for s, v in roi.items()}).T
z.columns = [f"{c / ca.FPS:+.0f} s" for c in z.columns]
print(z.describe().loc[["min", "50%", "max"]].round(2).to_string())
true_col = "+0 s"
best = z[true_col].abs() >= z.drop(columns=true_col).abs().max(axis=1)
print("sessions where the table's onsets give the largest |mean/SEM| of all shifts:", int(best.sum()), "of", len(z))
print("sessions where they do not:")
print(z[~best].round(2).to_string())
-20 s -10 s -5 s +0 s +5 s +10 s +20 s
min -5.15 -5.27 -5.74 1.63 -8.78 -8.33 -8.96
50% -0.74 -0.71 -0.08 9.28 -0.70 0.33 -0.11
max 4.82 6.62 2.91 31.28 5.31 4.93 4.77
sessions where the table's onsets give the largest |mean/SEM| of all shifts: 32 of 33
sessions where they do not:
-20 s -10 s -5 s +0 s +5 s +10 s +20 s
758519303 -2.36 -0.31 2.15 1.63 2.52 1.2 -0.42
The trial table's frame numbers agree with the 2p timestamps to within about one frame (printed above), the pupil series shares the 2p timestamps, and the running series spans the stimuli. In 32 of 33 sessions the population response at the table's stimulus onsets is larger, in units of its SEM, than at any shifted time; in the other, 758519303 (408021, session 1), the responses are weak at every shift (1.6 SEM at the onsets, 2.5 at +5 s), so the check is not informative there. The running series could not be checked by a contrast, as the claims use no running response.
dataset-holds-paper-sessions¶
The paper's data are the 50 sessions of Table 2, deposited on DANDI as dandiset 000037, with files named by subject and session ID. This claim checks whether dandiset 000336 holds any of those sessions.
Sessions are matched by subject and session ID from the file names of 000336 and the subject metadata of its files, against the IDs of Table 2.
In [12]:
paper_subjects, paper_sessions = set(table2.subject_id), set(table2.session_id)
print("Table 2:", len(paper_subjects), "subjects,", len(paper_sessions), "sessions")
print("000336 subjects in Table 2:", sorted(set(a336.subject) & paper_subjects))
print("000336 session IDs in Table 2:", sorted(set(a336.session.dropna()) & paper_sessions))
# subject metadata of one file per 000336 subject: the smallest processed (ophys) file, or for the
# subjects without one, the smallest raw (image+ophys) file
subj = []
for s, g in a336.groupby("subject"):
kind = "ophys" if (g.kind == "ophys").any() else "image+ophys"
r = g[g.kind == kind].sort_values("size").iloc[0]
with h5py.File(remfile.File(ca.s3_url(r.asset_id)), "r") as f:
subj.append({"subject": s, "file_kind": kind, "subject_id": f["general/subject/subject_id"][()].decode(),
"genotype": f["general/subject/genotype"][()].decode(),
"start": f["session_start_time"][()].decode()[:10]})
subj = pd.DataFrame(subj)
print(subj.to_string(index=False))
print("subjects opened:", len(subj), "of", a336.subject.nunique())
print("subject_id metadata in Table 2:", sorted(set(subj.subject_id) & paper_subjects))
print("000336 session start dates:", subj.start.min(), "to", subj.start.max(), "| Table 2 dates:", table2.date.min(), "to", table2.date.max())
print("000037 holds the 50 sessions of Table 2:", set(basic.session) == paper_sessions)
Table 2: 13 subjects, 50 sessions 000336 subjects in Table 2: [] 000336 session IDs in Table 2: []
subject file_kind subject_id genotype start 621602 ophys 621602 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2022-07-25 621603 ophys 621603 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2022-07-27 631188 ophys 631188 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2022-09-15 634402 ophys 634402 Cux2-CreERT2/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2022-10-04 635489 ophys 635489 Cux2-CreERT2/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2022-10-18 641331 ophys 641331 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2022-12-13 644972 ophys 644972 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-01-04 647715 ophys 647715 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-01-24 655042 ophys 655042 Cux2-CreERT2/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-01-26 655580 ophys 655580 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-02-08 656228 ophys 656228 Cux2-CreERT2/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-02-10 660986 ophys 660986 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-02-23 661749 ophys 661749 Cux2-CreERT2/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-03-07 662470 ophys 662470 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-03-13 690226 image+ophys 690226 Cux2-CreERT2/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-11-10 698808 image+ophys 698808 Rbp4-Cre_KL100/wt;Camk2a-tTA/wt;Ai93(TITL-GCaMP6f)/wt 2023-11-15 subjects opened: 16 of 16 subject_id metadata in Table 2: [] 000336 session start dates: 2022-07-25 to 2023-11-15 | Table 2 dates: 20180926 to 20190318 000037 holds the 50 sessions of Table 2: True
Result: not reproduced. None of the 16 subjects or 104 session IDs in the file names of DANDI:000336 is among the 13 subjects and 50 sessions of Table 2, nor is any subject ID in the metadata of the 16 subjects (one processed file each for 14, one raw file each for the 2 without processed files), and their sessions start between 2022-07-25 and 2023-11-15, against 2018 and 2019 for the paper's. 000336 holds a different project. All 50 sessions of Table 2 are in DANDI:000037, where the other claims are tested.
table1-subjects¶
13 mice were imaged, 3 in L2/3-D, 3 in L2/3-S, 4 in L5-D and 3 in L5-S, 7 male and 6 female, each with at least three imaging sessions.
Line from the genotype (Cux2 is L2/3, Rbp4 is L5), plane from the imaging plane depth, sex from the subject metadata, per subject, checked against Table 1. The date of birth is checked as the session date minus the stated age.
In [13]:
subj37 = md_.groupby("subject_id").agg(sex=("sex", "unique"), line=("line", "unique"), plane=("plane", "unique"),
sessions=("session", "size"), genotype=("genotype", "first"))
for c in ["sex", "line", "plane"]:
assert subj37[c].map(len).eq(1).all(), c
subj37[c] = subj37[c].str[0]
# date of birth as start date (local, see session-schedule) minus age in days
md_["start_utc"] = pd.to_datetime(md_.session_start_time.str[:19])
md_["local_date"] = md_.start_utc.dt.tz_localize("UTC").dt.tz_convert("America/Los_Angeles").dt.strftime("%Y%m%d")
md_["age_days"] = md_.age.str.extract(r"P(\d+)D")[0].astype(int)
dob = (pd.to_datetime(md_.local_date) - pd.to_timedelta(md_.age_days, unit="D")).dt.strftime("%Y%m%d")
subj37["dob_from_age"] = dob.groupby(md_.subject_id).agg(lambda x: sorted(set(x)))
t1 = table1.set_index("subject_id")
t1["plane"] = t1.plane.map({"somata": "soma", "dendrites": "dend"})
cmp1 = subj37.merge(t1, left_index=True, right_index=True, suffixes=("", "_table1"))
print(cmp1[["genotype", "sex", "sex_table1", "line", "line_table1", "plane", "plane_table1", "sessions", "dob_from_age", "dob"]].to_string())
print("sex, line and plane agree with Table 1 for every subject:",
all((cmp1[c] == cmp1[c + "_table1"]).all() for c in ["sex", "line", "plane"]))
grp = subj37.groupby(["line", "plane"]).size()
print(grp.to_string())
print("mice:", len(subj37), "| male:", (subj37.sex == "M").sum(), "| female:", (subj37.sex == "F").sum(),
"| fewest sessions per mouse:", subj37.sessions.min())
genotype sex sex_table1 line line_table1 plane plane_table1 sessions dob_from_age dob
subject_id
408021 Cux2-CreERT2;Camk2a-tTA;Ai93 M M L2/3 L2/3 soma soma 3 [20180623] 20180623
411400 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 F F L5 L5 dend dend 6 [20180711] 20180711
411424 Cux2-CreERT2;Camk2a-tTA;Ai93 F F L2/3 L2/3 soma soma 3 [20180711] 20180711
411771 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 M M L5 L5 soma soma 3 [20180713] 20180713
412933 Cux2-CreERT2;Camk2a-tTA;Ai93 M M L2/3 L2/3 dend dend 6 [20180718] 20180718
413663 Cux2-CreERT2;Camk2a-tTA;Ai93 M M L2/3 L2/3 dend dend 3 [20180721] 20180721
418779 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 F F L5 L5 soma soma 3 [20180820] 20180820
420011 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 F F L5 L5 dend dend 3 [20180826] 20180826
433414 Cux2-CreERT2;Camk2a-tTA;Ai93 F F L2/3 L2/3 dend dend 3 [20181106] 20181106
433448 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 M M L5 L5 dend dend 4 [20181106] 20181106
433451 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 M M L5 L5 dend dend 3 [20181106] 20181106
433458 Rbp4-Cre_KL100;Camk2a-tTA;Ai93 M M L5 L5 soma soma 4 [20181106] 20181106
440889 Cux2-CreERT2;Camk2a-tTA;Ai93 F F L2/3 L2/3 soma soma 6 [20181212] 20181212
sex, line and plane agree with Table 1 for every subject: True
line plane
L2/3 dend 3
soma 3
L5 dend 4
soma 3
mice: 13 | male: 7 | female: 6 | fewest sessions per mouse: 3
Result: reproduced. The files hold 13 mice: 3 L2/3-D, 3 L2/3-S, 4 L5-D and 3 L5-S, 7 male and 6 female, each with at least 3 sessions. Sex, line and plane agree with Table 1 for every subject, and the date of birth implied by each file's age agrees with Table 1.
table2-sessions¶
50 imaging sessions were recorded from the 13 mice. Two mice have sessions that failed QC; the sessions 1 to 3 that passed QC are 33 files.
The files carry no QC status; it comes from the paper's code table (mouse_df.csv), checked here
against Table 2. Readings for the mice with failed sessions: (a) mice whose sessions 1 to 3 include a
failed session, (b) mice with any failed session.
In [14]:
t2 = table2.set_index("session_id")
qc = mouse_df.set_index("sessid")
print("sessions in the files:", md_.session.nunique(), "| in Table 2:", len(t2), "| same set:", set(md_.session) == set(t2.index))
print("code table QC agrees with Table 2:", ((qc.loc[t2.index, "pass_fail"] == "P") == (t2.qc == "passed")).all())
failed = t2[t2.qc == "failed"]
print("failed sessions:", failed.groupby("subject_id").size().to_dict())
s123 = mouse_df[mouse_df.sess_n.between(1, 3)]
read_a = sorted(s123.loc[s123.pass_fail == "F", "mouseid"].unique())
read_b = sorted(failed.subject_id.unique())
print("(a) mice with a failed session among sessions 1 to 3:", len(read_a), read_a)
print("(b) mice with any failed session:", len(read_b), read_b)
print("sessions 1 to 3 that passed QC:", len(ANALYSED), "files, from", analysed.mouseid.nunique(), "mice")
print("sessions per mouse in the files:", md_.groupby("subject_id").size().to_dict())
sessions in the files: 50 | in Table 2: 50 | same set: True
code table QC agrees with Table 2: True
failed sessions: {'411400': 3, '412933': 3, '440889': 2}
(a) mice with a failed session among sessions 1 to 3: 2 ['411400', '412933']
(b) mice with any failed session: 3 ['411400', '412933', '440889']
sessions 1 to 3 that passed QC: 33 files, from 11 mice
sessions per mouse in the files: {'408021': 3, '411400': 6, '411424': 3, '411771': 3, '412933': 6, '413663': 3, '418779': 3, '420011': 3, '433414': 3, '433448': 4, '433451': 3, '433458': 4, '440889': 6}
Result: reproduced. The files hold 50 sessions, the set of Table 2. Under reading (a), 2 mice (411400 and 412933) have a failed session among sessions 1 to 3; under reading (b) the count is 3, with 440889, whose sessions 5 and 6 failed. 33 sessions 1 to 3 passed QC. The QC status is not in the files; it comes from the paper's code table, which agrees with Table 2.
session-schedule¶
After habituation each mouse had three 70-minute imaging sessions within three to six days, at most one per day.
The files' session_start_time is stored with a -08:00 offset, but its clock time equals the UTC start
time of the code table (full_timestamp), so the stored offset is wrong: the start times are UTC. Local
dates are taken by converting them from UTC to Pacific time, and compared with Table 2's dates. Readings:
which sessions, (a) sessions 1 to 3 as the code table numbers them, (b) the three tracked sessions (a
nonzero tracked count in the files); span, (i) days from first to third, (ii) calendar days inclusive;
duration, (x) first to last 2p frame, (y) first stimulus row start to last stop, grey screen included.
In [15]:
md_["start_clock"] = md_.session_start_time.str[:19].str.replace("-", "").str.replace(":", "")
print("stored offsets:", md_.session_start_time.str[-6:].unique())
print("start clock time equals the code table's UTC timestamp in every file:", (md_.start_clock == md_.full_timestamp).all())
md_["utc_date"] = md_.start_utc.dt.strftime("%Y%m%d")
print("local dates equal Table 2's dates:", (md_.set_index("session").local_date == t2.date.reindex(md_.session).values).all(),
"| UTC dates equal them:", (md_.set_index("session").utc_date == t2.date.reindex(md_.session).values).sum(), "of", len(md_))
md_["dur_x"] = (md_.twop_last - md_.twop_first) / 60
md_["dur_y"] = [(trials[s].stop_time.max() - trials[s].start_time.min()) / 60 for s in md_.session]
print("duration (x), min:", md_.dur_x.min().round(2), "to", md_.dur_x.max().round(2),
"| (y):", md_.dur_y.min().round(2), "to", md_.dur_y.max().round(2))
def spans(sel, datecol):
out = {}
for mouse, g in sel.groupby("subject_id"):
d = pd.to_datetime(g[datecol]).sort_values()
out[mouse] = {"n": len(d), "distinct": d.nunique() == len(d), "days_i": (d.iloc[-1] - d.iloc[0]).days,
"days_ii": (d.iloc[-1] - d.iloc[0]).days + 1, "dates": " ".join(d.dt.strftime("%m%d"))}
return pd.DataFrame(out).T
sel_a = md_[md_.sess_n.between(1, 3)]
sel_b = md_[md_.n_tracked > 0]
for name, sel in (("a", sel_a), ("b", sel_b)):
for datecol in ("local_date", "utc_date"):
sp = spans(sel, datecol)
ok_i = sp.distinct & sp.days_i.between(3, 6) & (sp.n == 3)
ok_ii = sp.distinct & sp.days_ii.between(3, 6) & (sp.n == 3)
print(f"\nreading ({name}), {datecol}: mice {len(sp)}, three sessions each: {(sp.n == 3).all()}; "
f"span (i) in 3-6 for {ok_i.sum()}, span (ii) in 3-6 for {ok_ii.sum()}")
if datecol == "local_date":
print(sp.to_string())
for name, sel in (("a", sel_a), ("b", sel_b)):
print(f"durations under ({name}): (x) {sel.dur_x.min():.2f} to {sel.dur_x.max():.2f} min (median {sel.dur_x.median():.2f}), "
f"(y) {sel.dur_y.min():.2f} to {sel.dur_y.max():.2f} min (median {sel.dur_y.median():.2f})")
stored offsets: <StringArray>
['-07:00', '-08:00']
Length: 2, dtype: str
start clock time equals the code table's UTC timestamp in every file: True
local dates equal Table 2's dates: True | UTC dates equal them: 46 of 50
duration (x), min: 70.23 to 70.53 | (y): 70.06 to 70.37
reading (a), local_date: mice 13, three sessions each: True; span (i) in 3-6 for 6, span (ii) in 3-6 for 13
n distinct days_i days_ii dates
408021 3 True 5 6 0926 0927 1001
411400 3 True 2 3 1001 1002 1003
411424 3 True 2 3 1009 1010 1011
411771 3 True 5 6 1003 1004 1008
412933 3 True 2 3 1016 1017 1018
413663 3 True 5 6 1017 1018 1022
418779 3 True 2 3 1112 1113 1114
420011 3 True 2 3 1113 1114 1115
433414 3 True 2 3 0220 0221 0222
433448 3 True 5 6 0220 0221 0225
433451 3 True 4 5 0215 0218 0219
433458 3 True 4 5 0221 0222 0225
440889 3 True 2 3 0306 0307 0308
reading (a), utc_date: mice 13, three sessions each: True; span (i) in 3-6 for 7, span (ii) in 3-6 for 13
reading (b), local_date: mice 13, three sessions each: True; span (i) in 3-6 for 8, span (ii) in 3-6 for 12
n distinct days_i days_ii dates
408021 3 True 5 6 0926 0927 1001
411400 3 True 6 7 1009 1011 1015
411424 3 True 2 3 1009 1010 1011
411771 3 True 5 6 1003 1004 1008
412933 3 True 3 4 1022 1023 1025
413663 3 True 5 6 1017 1018 1022
418779 3 True 2 3 1112 1113 1114
420011 3 True 2 3 1113 1114 1115
433414 3 True 2 3 0220 0221 0222
433448 3 True 5 6 0220 0221 0225
433451 3 True 4 5 0215 0218 0219
433458 3 True 4 5 0221 0222 0225
440889 3 True 2 3 0306 0307 0308
reading (b), utc_date: mice 13, three sessions each: True; span (i) in 3-6 for 8, span (ii) in 3-6 for 12
durations under (a): (x) 70.23 to 70.53 min (median 70.23), (y) 70.06 to 70.37 min (median 70.06)
durations under (b): (x) 70.23 to 70.32 min (median 70.23), (y) 70.06 to 70.16 min (median 70.06)
Result: reproduced. Under readings (a) and (ii), sessions 1 to 3 of each of the 13 mice fall on three distinct local dates within 3 to 6 calendar days, and every session lasts 70.2 to 70.5 minutes under duration reading (x) and 70.1 to 70.4 under (y). Under span (i) (days from the first session to the third) only 6 of 13 mice are within 3 to 6 days, as 7 had their sessions on three consecutive days (2 days). With the tracked sessions (b), 411400 spans 7 calendar days. The local dates match Table 2 for all 50 sessions once the stored start times are read as UTC; the UTC dates match for 46.
imaging-planes¶
Somata were imaged at 175 um (L2/3) and 375 um (L5) and distal apical dendrites at 50 to 75 um (L2/3) and 20 um (L5), at 30 Hz with 512 x 512 pixels over a 400 um field of view.
In [16]:
d = md_.set_index("session")
print("depth equals Table 2 in", (d.depth == t2.depth.reindex(d.index)).sum(), "of", len(d), "sessions")
mis = d[d.depth != t2.depth.reindex(d.index)]
print("sessions whose depth differs from Table 2:")
print(pd.DataFrame({"subject": mis.subject_id, "sess_n": mis.sess_n, "file_depth": mis.depth, "table2_depth": t2.depth.reindex(mis.index),
"code_table_depth": mouse_df.set_index("sessid").depth.reindex(mis.index)}).to_string())
print("raw imaging plane location strings of 433414:")
for s, loc in md_.loc[md_.subject_id == "433414", ["session", "location"]].values:
print(" ", s, loc)
code_depth = mouse_df.set_index("sessid").depth.reindex(d.index)
print("code table depth equals the file depth in", (code_depth == d.depth).sum(), "of", len(d),
"sessions, and equals Table 2 in", (code_depth == t2.depth.reindex(d.index)).sum())
print(md_.groupby(["line", "plane"]).depth.unique().to_string())
print("imaging rate metadata:", md_.imaging_rate.unique(), "Hz | from timestamps:",
(1 / md_.twop_median_dt).min().round(3), "to", (1 / md_.twop_median_dt).max().round(3), "Hz")
print("max projection image sizes:", md_.max_projection_shape.astype(str).value_counts().to_dict())
print("field of view in um: not stated in the files (imaging plane grid spacing not set)")
depth equals Table 2 in 47 of 50 sessions
sessions whose depth differs from Table 2:
subject sess_n file_depth table2_depth code_table_depth
session
826187862 433414 1 60 75 75
826773996 433414 2 60 75 75
827833392 433414 3 60 75 75
raw imaging plane location strings of 433414:
826187862 {'area': 'VISp', 'depth': '60'}
826773996 {'area': 'VISp', 'depth': '60'}
827833392 {'area': 'VISp', 'depth': '60'}
code table depth equals the file depth in 47 of 50 sessions, and equals Table 2 in 50
line plane
L2/3 dend [75, 50, 60]
soma [175]
L5 dend [20]
soma [375]
imaging rate metadata: [30.] Hz | from timestamps: 30.075 to 30.148 Hz
max projection image sizes: {'(512, 512)': 50}
field of view in um: not stated in the files (imaging plane grid spacing not set)
Result: partially reproduced. The imaging rate is 30 Hz in the metadata and 30.08 to 30.15 Hz from
the timestamps, and every plane image is 512 x 512 pixels. The depth equals Table 2 in 47 of 50
sessions: the three sessions of 433414 (L2/3-D) are at 60 um in the files and at 75 um in Table 2 and
in the paper's code table. The 400 um field of view is not in the files. The raw location string of each of the three 433414 files reads
'depth': '60', so 60 um is the deposited value, not a parsing artefact; the code table gives 75 um
there and equals Table 2 in all 50 sessions. The code table (printed in the data check) also agrees with
Table 2 on QC, date, ROI count and tracked ROI count in all 50 sessions, and its session numbers run 1,
2, ... in date order within each mouse.
roi-counts-per-session¶
15 to 224 ROIs per session were retained in the somatic planes and 159 to 1636 in the dendritic planes.
ROI count per session from the files, against Table 2. Readings: (a) the 33 analysed sessions, (b) every session that passed QC, (c) all 50 sessions. Variant: ROIs left after removing those with NaN or Inf values (traces read for the 33 analysed sessions only).
In [17]:
d = md_.set_index("session")
print("ROI count equals Table 2 in", (d.n_rois == t2.n_rois.reindex(d.index)).sum(), "of", len(d), "sessions")
passed = set(t2.index[t2.qc == "passed"])
for name, sel in (("a", d[d.index.isin(ANALYSED)]), ("b", d[d.index.isin(passed)]), ("c", d)):
g = sel.groupby("plane").n_rois.agg(["min", "max", "size"])
print(f"({name})", g.to_dict("index"))
clean = pd.Series({s: v["n_rois"] - v["n_bad"] for s, v in roi.items()})
print("variant, after NaN/Inf removal (a):", clean.groupby(d.loc[clean.index, "plane"]).agg(["min", "max"]).to_dict("index"),
"| ROIs removed:", int(sum(v["n_bad"] for v in roi.values())))
ROI count equals Table 2 in 50 of 50 sessions
(a) {'dend': {'min': 159, 'max': 1636, 'size': 15}, 'soma': {'min': 15, 'max': 224, 'size': 18}}
(b) {'dend': {'min': 159, 'max': 1636, 'size': 22}, 'soma': {'min': 15, 'max': 224, 'size': 20}}
(c) {'dend': {'min': 159, 'max': 1636, 'size': 28}, 'soma': {'min': 15, 'max': 227, 'size': 22}}
variant, after NaN/Inf removal (a): {'dend': {'min': 159, 'max': 1636}, 'soma': {'min': 15, 'max': 224}} | ROIs removed: 0
Result: reproduced. Every session's ROI count equals Table 2. Under readings (a) and (b) the somatic counts run from 15 to 224 and the dendritic from 159 to 1636, the paper's ranges exactly; under reading (c) the largest somatic count is 227, inside the widened range. No ROI of the 33 analysed sessions has NaN or Inf values, so the variant gives the same counts.
tracked-rois-per-plane¶
ROI matching across the three sessions gave 254 tracked ROIs in L2/3-D, 261 in L2/3-S, 516 in L5-D and 129 in L5-S, about 32 to 75% of the theoretical maximum.
Tracked ROIs are those with a tracking index (tracking_id not NaN). Per mouse, the tracked count must
be the same in its three tracked sessions, and the indices must run over the same set. Theoretical
maximum readings: (a) per plane, the sum over mice of the smallest ROI count among the three sessions;
(b) per plane, the sum of the mean count; (c) per mouse, as (a).
In [18]:
trk = []
for s, m in meta.items():
ids = m["tracking_ids"]
trk.append({"session": s, "subject_id": m["subject_id"], "n_tracked": m["n_tracked"], "n_rois": m["n_rois"],
"ids": tuple(sorted(ids[~np.isnan(ids)].astype(int)))})
trk = pd.DataFrame(trk).merge(md_[["session", "line", "plane"]], on="session")
trk["group"] = trk.line + "-" + trk.plane.map({"soma": "S", "dend": "D"})
print("tracked count equals Table 2 in", (trk.set_index("session").n_tracked == t2.n_tracked.reindex(trk.session).values).sum(), "of", len(trk))
tracked_sess = trk[trk.n_tracked > 0]
per_mouse = tracked_sess.groupby("subject_id").agg(group=("group", "first"), sessions=("session", "size"),
counts=("n_tracked", lambda x: sorted(set(x))),
same_ids=("ids", lambda x: len(set(x)) == 1),
min_rois=("n_rois", "min"), mean_rois=("n_rois", "mean"))
per_mouse["n_tracked"] = per_mouse.counts.str[0]
per_mouse["analysed"] = per_mouse.index.isin(analysed.mouseid)
per_mouse["tracked_in_sessions_1_3"] = [set(tracked_sess.loc[tracked_sess.subject_id == m, "session"]) <= set(analysed.sessid)
for m in per_mouse.index]
print(per_mouse.to_string())
pm = per_mouse[per_mouse.analysed]
by_plane = pm.groupby("group").agg(tracked=("n_tracked", "sum"), max_a=("min_rois", "sum"), max_b=("mean_rois", "sum"))
by_plane["frac_a_%"] = 100 * by_plane.tracked / by_plane.max_a
by_plane["frac_b_%"] = 100 * by_plane.tracked / by_plane.max_b
print(by_plane.round(1).to_string())
frac_c = 100 * pm.n_tracked / pm.min_rois
print("(c) per mouse:", frac_c.round(1).to_dict())
for name, x in (("a", by_plane["frac_a_%"]), ("b", by_plane["frac_b_%"]), ("c", frac_c)):
print(f"reading ({name}): {x.min():.1f} to {x.max():.1f} %, inside 28.8 to 82.5: {x.min() >= 28.8 and x.max() <= 82.5}")
tracked count equals Table 2 in 50 of 50
group sessions counts same_ids min_rois mean_rois n_tracked analysed tracked_in_sessions_1_3
subject_id
408021 L2/3-S 3 [59] True 74 92.333333 59 True True
411400 L5-D 3 [162] True 597 621.333333 162 False False
411424 L2/3-S 3 [55] True 80 85.666667 55 True True
411771 L5-S 3 [47] True 70 79.666667 47 True True
412933 L2/3-D 3 [98] True 168 254.000000 98 False False
413663 L2/3-D 3 [136] True 348 496.000000 136 True True
418779 L5-S 3 [12] True 15 23.333333 12 True True
420011 L5-D 3 [51] True 159 182.000000 51 True True
433414 L2/3-D 3 [118] True 300 453.333333 118 True True
433448 L5-D 3 [112] True 445 859.000000 112 True True
433451 L5-D 3 [353] True 966 993.666667 353 True True
433458 L5-S 3 [70] True 87 94.333333 70 True True
440889 L2/3-S 3 [147] True 210 219.333333 147 True True
tracked max_a max_b frac_a_% frac_b_%
group
L2/3-D 254 648 949.3 39.2 26.8
L2/3-S 261 364 397.3 71.7 65.7
L5-D 516 1570 2034.7 32.9 25.4
L5-S 129 172 197.3 75.0 65.4
(c) per mouse: {'408021': 79.7, '411424': 68.8, '411771': 67.1, '413663': 39.1, '418779': 80.0, '420011': 32.1, '433414': 39.3, '433448': 25.2, '433451': 36.5, '433458': 80.5, '440889': 70.0}
reading (a): 32.9 to 75.0 %, inside 28.8 to 82.5: True
reading (b): 25.4 to 65.7 %, inside 28.8 to 82.5: False
reading (c): 25.2 to 80.5 %, inside 28.8 to 82.5: False
Result: reproduced. Summed over the 11 analysed mice, the files hold 254 (L2/3-D), 261 (L2/3-S), 516 (L5-D) and 129 (L5-S) tracked ROIs, the paper's counts, and each mouse has the same tracking indices in its three tracked sessions. Under reading (a) of the theoretical maximum the fractions run from 32.9 to 75.0%, inside the widened range; under (b) they run from 25.4 to 65.7% and under (c) from 25.2 to 80.5%, both outside it.
fig8-example-matches¶
In the example L2/3-D mouse, taking the union of matches over all session orderings and removing four conflicting matches gave 136 final matches; in the example L5-S mouse, removing two conflicts gave 47.
The examples are mice 413663 (L2/3-D) and 411771 (L5-S), as the claims file gives them from the paper's code. The conflict counts cannot be tested from the files.
In [19]:
for mouse in ("413663", "411771"):
r = per_mouse.loc[mouse]
print(mouse, r.group, "tracked ROIs:", r.n_tracked, "| sessions:", r.sessions, "| same indices:", r.same_ids)
413663 L2/3-D tracked ROIs: 136 | sessions: 3 | same indices: True 411771 L5-S tracked ROIs: 47 | sessions: 3 | same indices: True
Result: reproduced. The files give 136 tracked ROIs for 413663 and 47 for 411771, the paper's final match counts. The conflicts removed (4 and 2) are not in the files and are not tested.
fig8-ordering-variation¶
The number of matched ROIs depends on the order in which sessions are aligned, and this variation is far smaller for somata than for dendrites.
Tested in the second notebook, roi_matching_orderings.ipynb, which reruns the matching with the public
n-way matching package over the six session orderings of each tracked mouse.
stimulus-durations¶
Each imaging session lasts 70 minutes, with 34 minutes of Gabor sequences and 17 minutes of visual flow in each direction, each stimulus period flanked by 30 s of grey screen.
From the trial table of each analysed session: each stimulus period is its first row start to its last row stop, per Gabor and per flow direction; grey periods are the grey-screen rows.
In [20]:
dur = []
for s in sorted(ANALYSED):
t = trials[s]
row = {"session": s}
for name, sel in (("gabors", t.stimulus_type == "gabors"),
("flow_left", t.stimulus_template_name == "visflow_left"),
("flow_right", t.stimulus_template_name == "visflow_right")):
row[name + "_min"] = (t.stop_time[sel].max() - t.start_time[sel].min()) / 60
g = t[t.stimulus_type == "grayscreen"]
row["grey_n"] = len(g)
row["grey_min_s"], row["grey_max_s"] = (g.stop_time - g.start_time).min(), (g.stop_time - g.start_time).max()
order = t.stimulus_template_name.values
row["order"] = " ".join(pd.Series(order)[pd.Series(order) != pd.Series(order).shift()].str.replace("visflow_", "").str.replace("grayscreen", "grey"))
dur.append(row)
dur = pd.DataFrame(dur).set_index("session")
print(dur.drop(columns="order").describe().loc[["min", "50%", "max"]].round(3))
print(dur.order.value_counts().to_string())
ok = (dur.gabors_min.sub(34).abs() <= 3.4) & dur.flow_left_min.sub(17).abs().le(1.7) & dur.flow_right_min.sub(17).abs().le(1.7) \
& dur.grey_min_s.sub(30).abs().le(3) & dur.grey_max_s.sub(30).abs().le(3)
print("sessions within tolerance:", ok.sum(), "of", len(ok))
gabors_min flow_left_min flow_right_min grey_n grey_min_s grey_max_s min 34.027 17.013 17.013 4.0 29.983 30.029 50% 34.028 17.014 17.014 4.0 29.983 30.031 max 34.127 17.043 17.016 4.0 29.983 30.039 order grey gabors grey right grey left grey 11 grey gabors grey left grey right grey 9 grey left grey right grey gabors grey 7 grey right grey left grey gabors grey 6 sessions within tolerance: 33 of 33
Result: reproduced. In all 33 sessions the Gabor period lasts 34.03 to 34.13 minutes, each flow direction 17.01 to 17.04 minutes, and each of the four grey periods 29.98 to 30.04 s, one before and after each stimulus period.
inconsistent-event-fraction¶
Inconsistent sequences are about 7% of Gabor sequences and 5% of visual flow time, in about 30 blocks per session, each 30 to 90 s of consistent stimulus followed by 3 to 6 s (Gabor) or 2 to 4 s (visual flow) of inconsistent stimulus.
Gabor fraction: sequences with a U image over all sequences (D or U). Flow fraction: time of rows with reversed squares over all flow time. A block is a run of consistent rows followed by a run of inconsistent rows, within the Gabor period and within each flow direction (the stimulus code samples them so). Durations are from row start and stop times. Readings: (a) mean of per-session fractions, (b) pooled; blocks (i) per stimulus, (ii) per session.
In [21]:
def runs(t):
# consistent and inconsistent run durations, in order, for one stimulus period (rows in time order)
u = t.unexpected.values == 1
edges = np.flatnonzero(np.diff(u.astype(int))) + 1
bounds = np.r_[0, edges, len(u)]
out = []
for a, b in zip(bounds[:-1], bounds[1:]):
out.append((bool(u[a]), t.stop_time.values[b - 1] - t.start_time.values[a]))
return out
inc = []
cons_d, inc_g, inc_f, block_runs = [], [], [], []
for s in sorted(ANALYSED):
t = trials[s]
gab = t[t.stimulus_type == "gabors"]
n_d, n_u = (gab.gabor_frame == "D").sum(), (gab.gabor_frame == "U").sum()
vf = t[t.stimulus_type == "visflow"]
vf_dur = vf.stop_time - vf.start_time
blocks = {}
for name, sel in (("gabors", gab), ("visflow_left", vf[vf.stimulus_template_name == "visflow_left"]),
("visflow_right", vf[vf.stimulus_template_name == "visflow_right"])):
r = runs(sel)
# a block: a consistent run followed by an inconsistent one; a trailing consistent run is not a block
nb = 0
for (u1, d1), (u2, d2) in zip(r[:-1], r[1:]):
if not u1 and u2:
nb += 1
cons_d.append((name, d1))
(inc_g if name == "gabors" else inc_f).append(d2)
block_runs += [(s, name, "consistent", d1), (s, name, "inconsistent", d2)]
blocks[name] = nb
blocks[name + "_trailing_s"] = r[-1][1] if not r[-1][0] else 0.0
blocks[name + "_first_inconsistent"] = r[0][0]
inc.append({"session": s, "gabor_U": n_u, "gabor_seq": n_d + n_u, "gabor_frac_%": 100 * n_u / (n_d + n_u),
"flow_incons_s": vf_dur[vf.unexpected == 1].sum(), "flow_s": vf_dur.sum(),
"flow_frac_%": 100 * vf_dur[vf.unexpected == 1].sum() / vf_dur.sum(),
"blocks_gabor": blocks["gabors"], "blocks_flow": blocks["visflow_left"] + blocks["visflow_right"],
"trailing_gabor_s": blocks["gabors_trailing_s"],
"starts_inconsistent": any(blocks[k + "_first_inconsistent"] for k in ("gabors", "visflow_left", "visflow_right"))})
inc = pd.DataFrame(inc).set_index("session")
inc["blocks_session"] = inc.blocks_gabor + inc.blocks_flow
print(inc.describe().loc[["mean", "min", "max"]].round(2).to_string())
print("(a) mean of per-session fractions: Gabor %.2f%%, flow %.2f%%" % (inc["gabor_frac_%"].mean(), inc["flow_frac_%"].mean()))
print("(b) pooled: Gabor %.2f%%, flow %.2f%%" % (100 * inc.gabor_U.sum() / inc.gabor_seq.sum(), 100 * inc.flow_incons_s.sum() / inc.flow_s.sum()))
print("(i) per stimulus, mean blocks: Gabor %.1f, flow %.1f; (ii) per session %.1f" % (inc.blocks_gabor.mean(), inc.blocks_flow.mean(), inc.blocks_session.mean()))
print("blocks per session (reading i): Gabor", inc.blocks_gabor.min(), "to", inc.blocks_gabor.max(),
"| flow", inc.blocks_flow.min(), "to", inc.blocks_flow.max())
print("sessions outside 27 to 33 blocks: Gabor", int((~inc.blocks_gabor.between(27, 33)).sum()),
"| flow", int((~inc.blocks_flow.between(27, 33)).sum()), "of", len(inc))
print(inc[["blocks_gabor", "blocks_flow"]].T.to_string())
cd = pd.DataFrame(cons_d, columns=["stim", "s"])
print("consistent run before each inconsistent run, s:", cd.groupby(cd.stim.str[:7]).s.agg(["min", "max", "size"]).round(2).to_dict("index"))
print("inconsistent runs, s: Gabor %.2f to %.2f (n %d); flow %.2f to %.2f (n %d)" % (min(inc_g), max(inc_g), len(inc_g), min(inc_f), max(inc_f), len(inc_f)))
print("Gabor inconsistent run lengths (s, rounded to 0.1):", pd.Series(np.round(inc_g, 1)).value_counts().sort_index().to_dict())
print("flow inconsistent run lengths (s, rounded to 0.1):", pd.Series(np.round(inc_f, 1)).value_counts().sort_index().to_dict())
out_c = cd[(cd.s < 30) | (cd.s > 90)]
print("consistent runs outside 30-90 s:", len(out_c), "of", len(cd), "| furthest outside by",
round(max((30 - out_c.s).max(), (out_c.s - 90).max()), 3) if len(out_c) else 0, "s")
out_g = [x for x in inc_g if not 3 <= x <= 6]
out_f = [x for x in inc_f if not 2 <= x <= 4]
print("inconsistent runs outside the stated range: Gabor", len(out_g), "of", len(inc_g), "(largest", round(max(out_g, default=0), 3), "s)",
"| flow", len(out_f), "of", len(inc_f), "(range", round(min(out_f, default=0), 3), "to", round(max(out_f, default=0), 3), "s)")
# The Methods sample durations on a grid (multiples of 1.5 s for Gabor, 1 s for flow); on that grid:
grid = np.where(cd.stim == "gabors", 1.5, 1.0)
cd_grid = np.round(cd.s / grid) * grid
print("on the sampling grid: consistent runs outside 30-90 s:", int(((cd_grid < 30) | (cd_grid > 90)).sum()),
"| Gabor inconsistent outside 3-6 s:", int(sum(not 3 <= np.round(x / 1.5) * 1.5 <= 6 for x in inc_g)),
"| flow inconsistent outside 2-4 s:", int(sum(not 2 <= np.round(x) <= 4 for x in inc_f)))
far = pd.DataFrame(block_runs, columns=["session", "stim", "run", "s"])
step = np.where(far.stim == "gabors", 1.5, 1.0)
far["grid"] = np.round(far.s / step) * step
far["off_s"] = far.s - far.grid
far["off_frames_60Hz"] = far.off_s * 60
print("runs furthest from their grid value:")
print(far.reindex(far.off_s.abs().sort_values(ascending=False).index).head(5).round(3).to_string(index=False))
print("runs more than 0.25 s from their grid value:", int((far.off_s.abs() > 0.25).sum()), "of", len(far),
"| by session:", far[far.off_s.abs() > 0.25].session.value_counts().to_dict())
# Measured over grid duration, per session (consistent runs): a session whose display ran slow
# stretches every run in proportion to its length.
far["stretch_per_mille"] = 1000 * far.off_s / far.grid
st_ = far[far.run == "consistent"].groupby("session").stretch_per_mille.agg(["median", "min", "max"])
print("consistent runs, measured over grid duration minus 1, per mille: median over sessions",
round(st_["median"].median(), 2), "| sessions with the largest median:")
print(st_.sort_values("median", ascending=False).head(4).round(2).to_string())
print("largest distance of a measured run from its grid value, s:",
round(max(np.abs(cd.s - cd_grid).max(), max(abs(x - np.round(x / 1.5) * 1.5) for x in inc_g), max(abs(x - np.round(x)) for x in inc_f)), 3))
print("per session: Gabor fraction", inc["gabor_frac_%"].min().round(2), "to", inc["gabor_frac_%"].max().round(2),
"%; flow fraction", inc["flow_frac_%"].min().round(2), "to", inc["flow_frac_%"].max().round(2), "%")
gabor_U gabor_seq gabor_frac_% flow_incons_s flow_s flow_frac_% blocks_gabor blocks_flow trailing_gabor_s blocks_session
mean 94.21 1360.0 6.93 94.98 2041.71 4.65 31.48 32.18 0.0 63.67
min 82.00 1360.0 6.03 80.07 2041.61 3.92 29.00 29.00 0.0 59.00
max 108.00 1360.0 7.94 112.09 2043.38 5.49 34.00 35.00 0.0 68.00
(a) mean of per-session fractions: Gabor 6.93%, flow 4.65%
(b) pooled: Gabor 6.93%, flow 4.65%
(i) per stimulus, mean blocks: Gabor 31.5, flow 32.2; (ii) per session 63.7
blocks per session (reading i): Gabor 29 to 34 | flow 29 to 35
sessions outside 27 to 33 blocks: Gabor 2 | flow 11 of 33
session 758519303 759189643 759660390 760260459 760659782 761269197 761624763 761944562 762250376 764704289 765193831 766502238 777496949 777914830 778374308 778864809 779152062 779650018 823453391 824434038 825180479 826187862 826338612 826659257 826773996 826819032 827300090 827833392 828475005 828816509 832883243 833704570 834403597
blocks_gabor 32 29 34 31 30 30 33 33 33 31 31 32 34 30 33 30 33 32 33 31 32 32 32 30 30 33 31 32 31 30 31 30 30
blocks_flow 31 34 33 32 29 34 34 30 35 33 34 29 30 32 31 29 30 32 30 34 34 34 35 31 33 31 33 31 33 34 30 35 32
consistent run before each inconsistent run, s: {'gabors': {'min': 30.02, 'max': 90.21, 'size': 1039}, 'visflow': {'min': 30.02, 'max': 90.07, 'size': 1062}}
inconsistent runs, s: Gabor 2.99 to 6.06 (n 1039); flow 1.99 to 4.01 (n 1062)
Gabor inconsistent run lengths (s, rounded to 0.1): {3.0: 349, 4.5: 348, 4.6: 1, 6.0: 339, 6.1: 2}
flow inconsistent run lengths (s, rounded to 0.1): {2.0: 387, 3.0: 342, 4.0: 333}
consistent runs outside 30-90 s: 49 of 2101 | furthest outside by 0.215 s
inconsistent runs outside the stated range: Gabor 377 of 1039 (largest 6.057 s) | flow 347 of 1062 (range 1.994 to 4.014 s)
on the sampling grid: consistent runs outside 30-90 s: 0 | Gabor inconsistent outside 3-6 s: 0 | flow inconsistent outside 2-4 s: 0
runs furthest from their grid value:
session stim run s grid off_s off_frames_60Hz
832883243 gabors consistent 88.887 88.5 0.387 23.225
833704570 gabors consistent 73.852 73.5 0.352 21.091
833704570 gabors consistent 79.836 79.5 0.336 20.151
833704570 gabors consistent 81.331 81.0 0.331 19.878
833704570 gabors consistent 72.316 72.0 0.316 18.941
runs more than 0.25 s from their grid value: 22 of 4202 | by session: {'833704570': 18, '832883243': 3, '826187862': 1}
consistent runs, measured over grid duration minus 1, per mille: median over sessions 0.81 | sessions with the largest median:
median min max
session
826187862 1.24 0.78 3.23
832883243 1.23 0.68 5.21
833704570 0.82 0.77 5.51
761944562 0.81 -2.18 1.10
largest distance of a measured run from its grid value, s: 0.387
per session: Gabor fraction 6.03 to 7.94 %; flow fraction 3.92 to 5.49 %
Result: reproduced. Inconsistent sequences are 6.93% of Gabor sequences and 4.65% of flow time (readings a and b agree; per session 6.03 to 7.94% and 3.92 to 5.49%). There are on average 31.5 Gabor blocks and 32.2 flow blocks per session (reading i, within 30 +- 3; reading ii, 63.7 per session, is not). Measured from row times, 49 of 2101 consistent runs exceed 90 s by up to 0.22 s, and 377 of 1039 Gabor and 347 of 1062 flow inconsistent runs fall outside 3 to 6 s and 2 to 4 s by at most 0.06 s; on the sampling grid of the Methods (1.5 s and 1 s steps) every run is within the stated ranges. The block count is judged on the mean per session: per session it runs from 29 to 34 for Gabor and 29 to 35 for flow, so 2 of 33 sessions (Gabor) and 11 of 33 (flow) fall outside 27 to 33. The run furthest from its grid value, a consistent Gabor run of 88.887 s in session 832883243 (440889), lies 0.387 s from 88.5 s; consistent runs are longer than their grid values by a median of 0.81 per mille, and the 22 runs more than 0.25 s off all lie in the three sessions with the largest stretch (1.24, 1.23 and 0.82 per mille medians, up to 5.5 per mille), consistent with the display running slightly slow in those sessions, not with a different design. The outcome rests on two choices the criterion does not fix: durations are judged on the sampling grid of the Methods, and the block count on its mean over sessions. On row times, or with the block count required in every session, the claim would not hold.
fig5-stimulus-responses¶
Somatic and dendritic ROIs respond to both the Gabor and the visual flow stimuli, many with increased fluorescence at stimulus onset, and their responses differ between consistent and inconsistent stimuli.
Session 1 of the 11 analysed mice, ROIs pooled per plane. Gabor: per ROI, mean dF/F over 0 to 0.6 s after D against after U. Flow, each direction: mean over 0 to 2 s after against 2 s before inconsistent flow onset. ROI selection (a) all ROIs, (b) the Fig. 5 legend's criteria. Tests (i) two-sided Wilcoxon signed-rank across ROIs, p < 0.05, (ii) |mean difference| > 2 SEM. Unscaled dF/F as the claims file says; the code's scaled traces are shown as a variant.
In [22]:
def two_tests(x, y):
dlt = np.asarray(y) - np.asarray(x)
n = len(dlt)
if n < 2:
return n, np.nan, np.nan, np.nan, False, False
p = scipy.stats.wilcoxon(dlt).pvalue if np.any(dlt != 0) else 1.0
mean, sem = dlt.mean(), dlt.std(ddof=1) / np.sqrt(n)
return n, mean, sem, p, p < 0.05, abs(mean) > 2 * sem
rows5 = []
for key in ("fig5", "fig5_scaled"):
for s in sorted(SESS1):
m = meta[s]
grp = md_.set_index("session").loc[s]
for stim, v in roi[s][key].items():
for sel_name in ("a", "b", "b_u") if stim == "gabors" else ("a", "b"):
keep = {"a": np.ones(len(v["consistent"]), bool), "b": v["selected"], "b_u": v.get("selected_u")}[sel_name]
rows5.append({"traces": key, "plane": grp.line + "-" + ("S" if grp.plane == "soma" else "D"),
"mouse": m["subject_id"], "stim": stim, "selection": sel_name,
"consistent": v["consistent"][keep], "inconsistent": v["inconsistent"][keep]})
r5 = pd.DataFrame(rows5)
res5 = []
for (traces, plane, stim, sel), g in r5.groupby(["traces", "plane", "stim", "selection"]):
x, y = np.concatenate(g.consistent.values), np.concatenate(g.inconsistent.values)
n, mean, sem, p, wil, sem2 = two_tests(x, y)
res5.append({"traces": traces, "plane": plane, "stim": stim, "selection": sel, "n_rois": n,
"mean_incons_minus_cons": mean, "sem": sem, "wilcoxon_p": p, "i_wilcoxon": wil, "ii_2sem": sem2})
res5 = pd.DataFrame(res5)
with pd.option_context("display.float_format", "{:.3g}".format):
print(res5[res5.traces == "fig5"].drop(columns="traces").to_string(index=False))
summary5 = res5[res5.selection != "b_u"].groupby(["traces", "selection", "plane"]).agg(any_i=("i_wilcoxon", "any"), any_ii=("ii_2sem", "any"))
print(summary5.to_string())
print("variant b_u (Gabor selection on inconsistent sequences, Fig. 5b legend), Gabor only:")
with pd.option_context("display.float_format", "{:.3g}".format):
print(res5[res5.selection == "b_u"].drop(columns=["selection"]).to_string(index=False))
# With b_u in place of b for Gabor: does every plane still have a difference under both tests?
alt = res5[(res5.selection == "b_u") | ((res5.selection == "b") & (res5.stim != "gabors"))]
print("selection b with b_u for Gabor, every plane has a difference:",
alt.groupby(["traces", "plane"]).agg(any_i=("i_wilcoxon", "any"), any_ii=("ii_2sem", "any")).all().to_dict())
for traces in ("fig5", "fig5_scaled"):
s5 = summary5.loc[traces]
print(traces, "| every reading pair has a difference in all four planes:", bool(s5.all().all()))
plane stim selection n_rois mean_incons_minus_cons sem wilcoxon_p i_wilcoxon ii_2sem
L2/3-D gabors a 1355 -0.00245 0.000376 3.06e-14 True True
L2/3-D gabors b 54 -0.0224 0.00454 3.37e-07 True True
L2/3-D gabors b_u 35 -0.0121 0.0078 0.368 False False
L2/3-D visflow left (nasal) a 1355 0.0081 0.000449 2.9e-74 True True
L2/3-D visflow left (nasal) b 56 0.0374 0.00392 7.55e-11 True True
L2/3-D visflow right (temp) a 1355 0.00743 0.000489 1.49e-60 True True
L2/3-D visflow right (temp) b 66 0.0354 0.00268 1.64e-12 True True
L2/3-S gabors a 407 0.00118 0.00185 0.717 False False
L2/3-S gabors b 16 0.0054 0.00752 0.159 False False
L2/3-S gabors b_u 20 0.00442 0.0146 0.784 False False
L2/3-S visflow left (nasal) a 407 0.00823 0.00226 1.06e-14 True True
L2/3-S visflow left (nasal) b 23 0.0808 0.0148 2.38e-07 True True
L2/3-S visflow right (temp) a 407 0.00773 0.00187 3.31e-17 True True
L2/3-S visflow right (temp) b 23 0.0686 0.0117 2.38e-07 True True
L5-D gabors a 2807 -0.000229 0.000338 0.883 False False
L5-D gabors b 124 -0.0104 0.00379 0.0255 True True
L5-D gabors b_u 95 0.0042 0.00419 0.543 False False
L5-D visflow left (nasal) a 2807 0.0014 0.000395 2.71e-08 True True
L5-D visflow left (nasal) b 110 0.0371 0.00314 8.75e-20 True True
L5-D visflow right (temp) a 2807 -7.22e-05 0.00039 0.464 False False
L5-D visflow right (temp) b 109 0.0273 0.00181 1.28e-19 True True
L5-S gabors a 204 0.0173 0.00589 0.000394 True True
L5-S gabors b 8 0.0691 0.0518 0.547 False False
L5-S gabors b_u 13 0.141 0.0749 0.0105 True False
L5-S visflow left (nasal) a 204 0.00696 0.00281 4.12e-08 True True
L5-S visflow left (nasal) b 13 0.0786 0.0146 0.000244 True True
L5-S visflow right (temp) a 204 -0.00173 0.00168 0.751 False False
L5-S visflow right (temp) b 10 0.0155 0.0019 0.00195 True True
any_i any_ii
traces selection plane
fig5 a L2/3-D True True
L2/3-S True True
L5-D True True
L5-S True True
b L2/3-D True True
L2/3-S True True
L5-D True True
L5-S True True
fig5_scaled a L2/3-D True True
L2/3-S True True
L5-D True True
L5-S True True
b L2/3-D True True
L2/3-S True True
L5-D True True
L5-S True True
variant b_u (Gabor selection on inconsistent sequences, Fig. 5b legend), Gabor only:
traces plane stim n_rois mean_incons_minus_cons sem wilcoxon_p i_wilcoxon ii_2sem
fig5 L2/3-D gabors 35 -0.0121 0.0078 0.368 False False
fig5 L2/3-S gabors 20 0.00442 0.0146 0.784 False False
fig5 L5-D gabors 95 0.0042 0.00419 0.543 False False
fig5 L5-S gabors 13 0.141 0.0749 0.0105 True False
fig5_scaled L2/3-D gabors 62 -0.00975 0.00641 0.53 False False
fig5_scaled L2/3-S gabors 21 0.0496 0.0534 0.946 False False
fig5_scaled L5-D gabors 116 -0.00355 0.00486 0.572 False False
fig5_scaled L5-S gabors 13 0.428 0.261 0.00464 True False
selection b with b_u for Gabor, every plane has a difference: {'any_i': True, 'any_ii': True}
fig5 | every reading pair has a difference in all four planes: True
fig5_scaled | every reading pair has a difference in all four planes: True
In [23]:
# Selection (b) per session: ROIs selected out of all, and the session thresholds (75th percentiles;
# SNR median). Unscaled traces.
selrows = []
for s in sorted(SESS1):
grp = md_.set_index("session").loc[s]
for stim, v in roi[s]["fig5"].items():
for name, sel_key, thr_key in (("b", "selected", "thresholds"), ("b_u", "selected_u", "thresholds_u")):
if sel_key not in v:
continue
selrows.append({"session": s, "mouse": grp.subject_id, "plane": grp.line + "-" + ("S" if grp.plane == "soma" else "D"),
"stim": stim, "selection": name, "selected": int(v[sel_key].sum()), "rois": len(v[sel_key]),
**{f"thr_{k}": t for k, t in v[thr_key].items()}})
selrows = pd.DataFrame(selrows)
with pd.option_context("display.float_format", "{:.3g}".format):
print(selrows.to_string(index=False))
print("sessions x stimuli with no ROI selected:", selrows.groupby("selection").apply(lambda g: int((g.selected == 0).sum())).to_dict())
session mouse plane stim selection selected rois thr_snr thr_corr thr_sd thr_skew thr_diff
758519303 408021 L2/3-S gabors b 2 96 3.85 0.0567 0.0227 1 NaN
758519303 408021 L2/3-S gabors b_u 2 96 3.85 0.0941 0.0396 0.649 NaN
758519303 408021 L2/3-S visflow left (nasal) b 4 96 3.85 0.0107 NaN NaN 0.00688
758519303 408021 L2/3-S visflow right (temp) b 8 96 3.85 0.0137 NaN NaN 0.0203
760260459 411771 L5-S gabors b 4 90 3.65 0.0336 0.0165 0.731 NaN
760260459 411771 L5-S gabors b_u 5 90 3.65 0.073 0.0259 0.841 NaN
760260459 411771 L5-S visflow left (nasal) b 3 90 3.65 0.00696 NaN NaN 0.0158
760260459 411771 L5-S visflow right (temp) b 2 90 3.65 0.00804 NaN NaN 0.00557
761624763 411424 L2/3-S gabors b 3 87 3.01 0.0604 0.0239 0.824 NaN
761624763 411424 L2/3-S gabors b_u 6 87 3.01 0.0689 0.0263 0.791 NaN
761624763 411424 L2/3-S visflow left (nasal) b 8 87 3.01 0.0305 NaN NaN 0.0183
761624763 411424 L2/3-S visflow right (temp) b 5 87 3.01 0.0226 NaN NaN 0.0201
764704289 413663 L2/3-D gabors b 25 628 2.39 0.00631 0.0108 0.498 NaN
764704289 413663 L2/3-D gabors b_u 17 628 2.39 0.00541 0.0252 0.313 NaN
764704289 413663 L2/3-D visflow left (nasal) b 27 628 2.39 0.0159 NaN NaN 0.0165
764704289 413663 L2/3-D visflow right (temp) b 31 628 2.39 0.0109 NaN NaN 0.013
777496949 418779 L5-S gabors b 0 15 3.55 0.0102 0.0108 0.564 NaN
777496949 418779 L5-S gabors b_u 1 15 3.55 0.00873 0.0114 0.371 NaN
777496949 418779 L5-S visflow left (nasal) b 3 15 3.55 0.0262 NaN NaN 0.000536
777496949 418779 L5-S visflow right (temp) b 1 15 3.55 0.00646 NaN NaN 0.0041
777914830 420011 L5-D gabors b 13 205 2.54 0.0019 0.0146 0.406 NaN
777914830 420011 L5-D gabors b_u 5 205 2.54 0.00386 0.0423 0.335 NaN
777914830 420011 L5-D visflow left (nasal) b 5 205 2.54 0.0104 NaN NaN 0.00794
777914830 420011 L5-D visflow right (temp) b 7 205 2.54 0.0208 NaN NaN 0.0097
823453391 433451 L5-D gabors b 54 966 2.58 0.0125 0.018 0.696 NaN
823453391 433451 L5-D gabors b_u 49 966 2.58 0.00845 0.0293 0.408 NaN
823453391 433451 L5-D visflow left (nasal) b 33 966 2.58 0.0108 NaN NaN 0.0067
823453391 433451 L5-D visflow right (temp) b 30 966 2.58 0.0137 NaN NaN 0.0107
826187862 433414 L2/3-D gabors b 29 727 2.39 0.0127 0.0171 0.405 NaN
826187862 433414 L2/3-D gabors b_u 18 727 2.39 0.0119 0.033 0.33 NaN
826187862 433414 L2/3-D visflow left (nasal) b 29 727 2.39 0.0144 NaN NaN 0.0147
826187862 433414 L2/3-D visflow right (temp) b 35 727 2.39 0.0119 NaN NaN 0.0157
826338612 433448 L5-D gabors b 57 1636 2.54 0.00399 0.0143 0.271 NaN
826338612 433448 L5-D gabors b_u 41 1636 2.54 0.00556 0.0397 0.315 NaN
826338612 433448 L5-D visflow left (nasal) b 72 1636 2.54 0.0112 NaN NaN 0.0108
826338612 433448 L5-D visflow right (temp) b 72 1636 2.54 0.0116 NaN NaN 0.00686
826659257 433458 L5-S gabors b 4 99 3 0.0328 0.0162 0.838 NaN
826659257 433458 L5-S gabors b_u 7 99 3 0.101 0.0287 0.886 NaN
826659257 433458 L5-S visflow left (nasal) b 7 99 3 0.0396 NaN NaN 0.0223
826659257 433458 L5-S visflow right (temp) b 7 99 3 0.00621 NaN NaN 0.0107
832883243 440889 L2/3-S gabors b 11 224 2.89 0.045 0.0117 0.72 NaN
832883243 440889 L2/3-S gabors b_u 12 224 2.89 0.0454 0.0156 0.543 NaN
832883243 440889 L2/3-S visflow left (nasal) b 11 224 2.89 0.0221 NaN NaN 0.0115
832883243 440889 L2/3-S visflow right (temp) b 10 224 2.89 0.0176 NaN NaN 0.00955
sessions x stimuli with no ROI selected: {'b': 1, 'b_u': 0}
In [24]:
# By mouse: Gabor U minus D (all ROIs, unscaled), and the per-mouse Wilcoxon p.
bym = []
for _, g in r5[(r5.traces == "fig5") & (r5.selection == "a")].iterrows():
n, mean, sem, p, wil, sem2 = two_tests(g.consistent, g.inconsistent)
bym.append({"plane": g.plane, "mouse": g.mouse, "stim": g.stim, "n": n, "mean": mean, "sem": sem, "p": p})
with pd.option_context("display.float_format", "{:.3g}".format):
print(pd.DataFrame(bym).pivot_table(index=["plane", "mouse"], columns="stim", values=["mean", "p"]).to_string())
mean p
stim gabors visflow left (nasal) visflow right (temp) gabors visflow left (nasal) visflow right (temp)
plane mouse
L2/3-D 413663 -0.00664 0.00981 0.00686 5.53e-50 3.64e-50 3.38e-23
433414 0.00117 0.00662 0.00793 0.000455 3.04e-28 3.21e-39
L2/3-S 408021 0.00546 -0.00358 0.0117 0.00245 0.31 0.00197
411424 -0.00424 0.0173 0.0137 9.53e-07 5.06e-10 8.81e-08
440889 0.00145 0.00975 0.0037 0.967 1.17e-08 1.21e-09
L5-D 420011 -6.62e-05 0.00221 0.00231 0.524 0.0248 0.134
433448 0.00521 0.00242 -0.00178 1.01e-60 1.7e-08 5.5e-08
433451 -0.00948 -0.000488 0.00232 6.61e-78 0.287 2.73e-14
L5-S 411771 0.0221 0.00786 -0.00127 0.00204 5.51e-07 0.774
418779 -0.00249 -0.00943 -0.00814 0.151 0.0554 0.277
433458 0.016 0.00862 -0.00118 0.00841 0.000276 0.219
Result: reproduced. Under all four reading pairs (selection a or b, test i or ii), each of the four planes has at least one stimulus with a difference between consistent and inconsistent responses, and the same holds with the code's scaled traces. With all ROIs (a), the Gabor difference is significant in L2/3-D and L5-S but not in L2/3-S or L5-D, where only the flow differs. Under selection (b) the flow selection is made on the after-minus-before difference itself, so its test is not independent of the selection; that selection is the claims file's own. Selecting Gabor ROIs on their inconsistent (U) sequence responses, as the Fig. 5b legend can be read, gives no Gabor difference in L2/3-D, L2/3-S or L5-D under either test (L5-S, 13 ROIs, Wilcoxon p = 0.011, mean within 2 SEM), but every plane keeps a difference through flow, so the outcome is the same. Selection (b) picks 0 to 72 ROIs per session and stimulus; one session (418779, Gabor) selects none.
fig6-behaviour¶
Running velocity and pupil diameter were recorded in every analysed session, for the 11 mice whose data are shown aggregated in Fig. 6b.
Each analysed session must have a running velocity trace (cm/s) and a pupil diameter trace (mm) that span the session, here the stimulus period from the first to the last trial-table row. Running is median-filtered (5 frames) and single-frame changes of at least 50 cm/s are dropped; pupil frames with a single-frame change of at least 0.05 mm are dropped (Methods).
In [25]:
cov = []
for s in sorted(ANALYSED):
m, t, b = meta[s], trials[s], beh[s]
t0, t1 = t.start_time.min(), t.stop_time.max()
cov.append({"session": s, "mouse": m["subject_id"], "run_unit": m["run_unit"], "pupil_unit": m["pupil_unit"],
"run_spans": m["run_first"] <= t0 + 1 and m["run_last"] >= t1 - 1,
"pupil_spans": m["pupil_first"] <= t0 + 1 and m["pupil_last"] >= t1 - 1,
"run_n": m["run_n"], "run_kept_%": 100 * len(b["run"]) / m["run_n"],
"pupil_nan_%": 100 * b["pupil_nan"] / m["pupil_n"], "pupil_kept_%": 100 * len(b["pupil"]) / m["pupil_n"],
"run_median": np.median(b["run"]), "pupil_median": np.median(b["pupil"])})
cov = pd.DataFrame(cov)
print(cov.drop(columns=["session"]).describe().loc[["min", "50%", "max"]].round(3).to_string())
print("units:", cov.run_unit.unique(), cov.pupil_unit.unique())
ok6 = cov[cov.run_spans & cov.pupil_spans & (cov.run_unit == "cm/s") & (cov.pupil_unit == "mm")]
print("sessions with both traces spanning the session:", len(ok6), "of", len(cov), "| mice:", ok6.mouse.nunique())
print(cov.groupby("mouse")[["pupil_nan_%", "run_median", "pupil_median"]].median().round(3).to_string())
run_n run_kept_% pupil_nan_% pupil_kept_% run_median pupil_median
min 251999.0 99.995 0.062 97.699 0.006 0.314
50% 251999.0 100.000 0.396 99.601 0.130 0.487
max 251999.0 100.000 2.293 99.937 33.141 0.678
units: <StringArray>
['cm/s']
Length: 1, dtype: str <StringArray>
['mm']
Length: 1, dtype: str
sessions with both traces spanning the session: 33 of 33 | mice: 11
pupil_nan_% run_median pupil_median
mouse
408021 0.986 0.191 0.509
411424 0.381 3.543 0.643
411771 0.331 0.019 0.318
413663 0.247 0.078 0.487
418779 1.038 0.125 0.447
420011 0.660 0.159 0.464
433414 0.122 0.015 0.487
433448 0.133 19.983 0.629
433451 1.004 0.055 0.425
433458 0.499 4.198 0.612
440889 0.389 0.196 0.405
In [26]:
fig, ax = plt.subplots(2, 3, figsize=(11, 5), sharey="row")
for col, sn in enumerate((1, 2, 3)):
ss = analysed.loc[analysed.sess_n == sn, "sessid"]
run = [beh[s]["run"] for s in ss]
pup = [beh[s]["pupil"] for s in ss]
ax[0, col].hist(run, bins=np.linspace(-20, 80, 51), stacked=True, color=plt.cm.Greys(np.linspace(0.3, 0.9, len(run))))
ax[1, col].hist(pup, bins=np.linspace(0, 1.5, 51), stacked=True, color=plt.cm.Greys(np.linspace(0.3, 0.9, len(pup))))
ax[0, col].set_title(f"session {sn}")
out_r = sum(int(((x < -20) | (x > 80)).sum()) for x in run)
out_p = sum(int(((x < 0) | (x > 1.5)).sum()) for x in pup)
print(f"session {sn}: running {min(x.min() for x in run):.1f} to {max(x.max() for x in run):.1f} cm/s, {out_r} values off the axes; "
f"pupil {min(x.min() for x in pup):.3f} to {max(x.max() for x in pup):.3f} mm, {out_p} off the axes")
ax[0, 0].set_ylabel("running frames"); ax[1, 0].set_ylabel("pupil frames")
for a in ax[0]: a.set_xlabel("running velocity (cm/s)")
for a in ax[1]: a.set_xlabel("pupil diameter (mm)")
fig.tight_layout()
session 1: running -28.4 to 74.1 cm/s, 147 values off the axes; pupil 0.082 to 1.021 mm, 0 off the axes
session 2: running -22.1 to 69.1 cm/s, 44 values off the axes; pupil 0.034 to 0.962 mm, 0 off the axes
session 3: running -18.4 to 83.8 cm/s, 146 values off the axes; pupil 0.102 to 0.919 mm, 0 off the axes
Result: reproduced. All 33 analysed sessions, from 11 mice, have a running velocity trace in cm/s and a pupil diameter trace in mm that span the stimulus period; 0.06 to 2.3% of pupil frames are NaN.
fig7a-snr-high-stable¶
ROI SNR is relatively high, above 1, in every plane and session, and stable over the three days.
SNR per ROI as in the paper's code, ROIs pooled across mice within each line, plane and session number (12 groups). "Above 1": (a) group median > 1, (b) group 25th percentile > 1. "Stable": (i) within each plane, each session's median lies inside every other session's 25th to 75th percentile box; (ii) each session's median within 10% of the mean of the three medians.
In [27]:
def per_roi(key):
rows = []
for s, v in roi.items():
m = md_.set_index("session").loc[s]
rows.append(pd.DataFrame({"value": v[key], "session": s, "mouse": m.subject_id, "sess_n": m.sess_n,
"plane": m.line + "-" + ("S" if m.plane == "soma" else "D")}))
return pd.concat(rows, ignore_index=True)
def group_stats(df):
g = df.groupby(["plane", "sess_n"]).value
return pd.DataFrame({"n": g.size(), "q25": g.quantile(0.25), "median": g.median(), "q75": g.quantile(0.75),
"q05": g.quantile(0.05), "q95": g.quantile(0.95)})
def stability(st):
out = {}
for plane, g in st.groupby(level="plane"):
med, q25, q75 = g["median"].values, g.q25.values, g.q75.values
inside = all(q25[j] <= med[i] <= q75[j] for i in range(3) for j in range(3) if i != j)
within = bool(np.all(np.abs(med - med.mean()) <= 0.1 * med.mean()))
out[plane] = {"i_inside_boxes": inside, "ii_within_10%": within,
"max_dev_from_mean_%": 100 * np.max(np.abs(med - med.mean())) / med.mean()}
return pd.DataFrame(out).T
snr = per_roi("snr")
st = group_stats(snr)
print(st.round(3).to_string())
print("(a) every group median > 1:", bool((st["median"] > 1).all()), "| (b) every group 25th percentile > 1:", bool((st.q25 > 1).all()),
"| smallest 25th percentile:", st.q25.min().round(3))
print(stability(st).to_string())
print("ROIs with SNR <= 1:", int((snr.value <= 1).sum()), "of", len(snr))
print(snr.groupby(["plane", "mouse", "sess_n"]).value.median().unstack().round(2).to_string())
n q25 median q75 q05 q95
plane sess_n
L2/3-D 1 1355 2.337 2.388 2.470 2.284 2.925
2 648 2.335 2.382 2.460 2.292 3.075
3 845 2.318 2.372 2.442 2.271 2.885
L2/3-S 1 407 2.529 3.090 4.593 2.285 9.067
2 388 2.479 2.924 4.173 2.273 7.947
3 397 2.475 2.831 3.769 2.273 7.605
L5-D 1 2807 2.455 2.552 2.688 2.367 3.178
2 1633 2.438 2.540 2.707 2.356 3.341
3 1664 2.431 2.519 2.660 2.350 3.132
L5-S 1 204 2.793 3.430 4.632 2.411 8.130
2 183 2.673 3.226 4.161 2.359 6.885
3 205 2.606 3.188 4.138 2.308 7.467
(a) every group median > 1: True | (b) every group 25th percentile > 1: True | smallest 25th percentile: 2.318
i_inside_boxes ii_within_10% max_dev_from_mean_%
L2/3-D True True 0.365022
L2/3-S True True 4.799678
L5-D True True 0.70787
L5-S True True 4.533112
ROIs with SNR <= 1: 0 of 10736
sess_n 1 2 3
plane mouse
L2/3-D 413663 2.39 2.39 2.36
433414 2.39 2.38 2.38
L2/3-S 408021 3.85 3.80 3.01
411424 3.01 2.97 2.84
440889 2.89 2.75 2.73
L5-D 420011 2.54 2.52 2.51
433448 2.54 2.47 2.48
433451 2.58 2.58 2.54
L5-S 411771 3.65 3.38 3.41
418779 3.55 3.13 2.95
433458 3.00 3.10 2.99
In [28]:
# The SNR of pure Gaussian noise (B1 of the first verification): unit-SD noise traces of a session's
# length with mean 0, +0.5 and -0.5 SD. For Gaussian noise of mean m and SD sd the SNR is about
# m / sd + E[Z | Z > 1.645] = m / sd + 2.063, so the measure exceeds 1 for noise alone unless the
# trace's median lies below about -1.06 noise SD.
rng_snr = np.random.default_rng(71)
n_frames = int(md_[md_.analysed].twop_n.median())
floor = []
for m0 in (0.0, 0.5, -0.5):
x = rng_snr.normal(m0, 1.0, size=(50, n_frames))
v = np.array([ca.calculate_snr(r)[0] for r in x])
floor.append({"noise_mean_sd": m0, "snr_median": np.median(v), "snr_min": v.min(), "snr_max": v.max(),
"analytic": m0 + scipy.stats.norm.pdf(1.645) / 0.05})
print(pd.DataFrame(floor).round(3).to_string(index=False))
# Each ROI's own noise floor: m / sd + 2.063, with m its median and sd its noise SD (signal / SNR).
fl = []
for s, v in roi.items():
sd = v["signal"] / v["snr"]
m = md_.set_index("session").loc[s]
fl.append(pd.DataFrame({"plane": m.line + "-" + ("S" if m.plane == "soma" else "D"), "sess_n": m.sess_n,
"snr": v["snr"], "floor": v["median"] / sd + scipy.stats.norm.pdf(1.645) / 0.05}))
fl = pd.concat(fl, ignore_index=True)
fl["above_floor"] = fl.snr - fl.floor
print(fl.groupby(["plane", "sess_n"]).agg(n=("snr", "size"), snr_median=("snr", "median"), floor_median=("floor", "median"),
above_floor_median=("above_floor", "median"),
frac_above_floor_by_1=("above_floor", lambda x: (x > 1).mean())).round(3).to_string())
noise_mean_sd snr_median snr_min snr_max analytic
0.0 2.063 2.050 2.080 2.062
0.5 2.563 2.549 2.579 2.562
-0.5 1.562 1.547 1.574 1.562
n snr_median floor_median above_floor_median frac_above_floor_by_1
plane sess_n
L2/3-D 1 1355 2.388 2.062 0.326 0.041
2 648 2.382 2.062 0.319 0.051
3 845 2.372 2.062 0.309 0.031
L2/3-S 1 407 3.090 2.062 1.028 0.506
2 388 2.924 2.062 0.862 0.448
3 397 2.831 2.062 0.769 0.398
L5-D 1 2807 2.552 2.062 0.489 0.070
2 1633 2.540 2.062 0.477 0.088
3 1664 2.519 2.062 0.456 0.064
L5-S 1 204 3.430 2.062 1.368 0.618
2 183 3.226 2.062 1.164 0.568
3 205 3.188 2.062 1.126 0.546
In [29]:
# The helper's calculate_snr and rolling_mean against the functions of the paper's code
# (github.com/colleenjg/util, util/math_util.py, at ca.UTIL_COMMIT), and its pairwise correlations against
# the code's own correlation step (get_correlation, analysis/misc_analys.py at ca.CODE_COMMIT), on every
# ROI of one somatic and one dendritic session-1 file (418779 and 420011).
paper_fn = ca.paper_code_functions()
cmp_rows = []
for s in sorted(SESS1):
if meta[s]["subject_id"] not in ("418779", "420011"):
continue
with ca.open_h5(basic.loc[basic.session == s, "asset_id"].iloc[0]) as f:
dff = f["processing/ophys/DfOverF/RoiResponseSeries/data"][:].T
dff = dff[np.all(np.isfinite(dff), axis=1)]
snr_paper = np.array([paper_fn["calculate_snr"](x) for x in dff])
rm_paper, rm_here = paper_fn["rolling_mean"](dff, win=4), ca.rolling_mean(dff, 4)
c_paper = ca.paper_get_correlation(dff)
cmp_rows.append({"session": s, "mouse": meta[s]["subject_id"], "rois": len(dff),
"max_abs_diff_snr": np.abs(snr_paper - roi[s]["snr"]).max(),
"max_abs_diff_rolling_mean": np.abs(rm_paper - rm_here).max(),
"pairs": len(c_paper), "max_abs_diff_corr": np.abs(c_paper - roi[s]["corr"]).max()})
del dff, rm_paper, rm_here
print(pd.DataFrame(cmp_rows).to_string(index=False))
session mouse rois max_abs_diff_snr max_abs_diff_rolling_mean pairs max_abs_diff_corr 777496949 418779 15 0.0 8.881784e-16 105 3.445616e-09 777914830 420011 205 0.0 1.776357e-15 20910 2.898118e-08
Result: reproduced. Every group's median and 25th percentile of SNR are above 1 (smallest 25th percentile 2.32; no ROI has SNR at or below 1), and in all four planes each session's median lies within every other session's interquartile box and within 4.8% of the mean of the three medians, so every reading pair holds. The criterion's threshold of 1 is low for this measure: pure Gaussian noise gives an SNR of 2.06 (1.56 and 2.56 with the median 0.5 SD below or above the noise mean), matching the analytic m / sd + 2.063. Every group's median SNR lies above its ROIs' own Gaussian noise floor (median floor 2.062), by 0.31 to 0.33 in L2/3-D, 0.46 to 0.49 in L5-D and 0.77 to 1.37 in the somatic planes; 3 to 9% of dendritic and 40 to 62% of somatic ROIs exceed it by more than 1. The outcome follows the criterion, but SNR above 1 alone would hold for noise. The helper's calculate_snr and rolling mean equal those of the paper's code (util/math_util.py at commit 2b05f954), and its pairwise correlations equal those of the code's own correlation step (get_correlation, analysis/misc_analys.py at commit 046e7fe), on every ROI and pair of 418779 and 420011 session 1 (105 and 20910 pairs; largest differences 0, under 2e-15 and under 3e-8, the last from storing correlations as float32).
fig7b-signal-stable¶
The mean dF/F signal of the ROIs is stable over days.
Signal is the mean of an ROI's dF/F values above its noise threshold. Readings (i) and (ii) as for SNR.
In [30]:
sig = per_roi("signal")
st_b = group_stats(sig)
print(st_b.round(4).to_string())
print(stability(st_b).to_string())
print(sig.groupby(["plane", "mouse", "sess_n"]).value.median().unstack().round(3).to_string())
n q25 median q75 q05 q95
plane sess_n
L2/3-D 1 1355 0.3616 0.4674 0.5749 0.2551 0.7414
2 648 0.3524 0.4481 0.5357 0.2440 0.6903
3 845 0.3334 0.4323 0.5345 0.2418 0.7068
L2/3-S 1 407 0.1601 0.2015 0.2839 0.1328 0.6628
2 388 0.1573 0.1908 0.2653 0.1299 0.5760
3 397 0.1620 0.1871 0.2424 0.1297 0.5403
L5-D 1 2807 0.3749 0.5168 0.7086 0.2512 0.9884
2 1633 0.3515 0.4461 0.5815 0.2511 0.8850
3 1664 0.3537 0.4541 0.6041 0.2554 0.9197
L5-S 1 204 0.2137 0.2529 0.3551 0.1671 0.5970
2 183 0.1857 0.2216 0.2891 0.1311 0.4606
3 205 0.1954 0.2367 0.3090 0.1620 0.5594
i_inside_boxes ii_within_10% max_dev_from_mean_%
L2/3-D True True 4.02721
L2/3-S True True 4.350542
L5-D True True 9.413938
L5-S True True 6.694551
sess_n 1 2 3
plane mouse
L2/3-D 413663 0.405 0.417 0.392
433414 0.522 0.482 0.490
L2/3-S 408021 0.248 0.218 0.186
411424 0.213 0.212 0.200
440889 0.185 0.182 0.184
L5-D 420011 0.521 0.526 0.557
433448 0.592 0.443 0.444
433451 0.440 0.435 0.444
L5-S 411771 0.240 0.188 0.217
418779 0.226 0.221 0.262
433458 0.278 0.252 0.257
Result: reproduced. Both stability readings hold in all four planes: each session's median signal lies inside the other sessions' boxes and within 9.4% (L5-D, the largest) of the mean of the three medians.
fig7c-correlations-similar¶
The distributions of pairwise ROI correlations are very similar for apical dendrite ROIs and for somatic ROIs, so the dendrites are unlikely to be heavily over-segmented.
Pearson correlations over all ROI pairs of each session (4-sample rolling mean), pooled across the mice of each line, plane and session. Six comparisons (line x session), dendritic against somatic: (i) the dendritic median inside the somatic 25th to 75th percentile range; (ii) the dendritic 99th percentile at most 10% above the somatic 99th percentile.
In [31]:
corr_rows = []
pooled = {}
for (line, plane, sn), g in md_[md_.analysed].groupby(["line", "plane", "sess_n"]):
c = np.concatenate([roi[s]["corr"] for s in g.session])
pooled[(line, plane, sn)] = c
q = np.quantile(c, [0.25, 0.5, 0.75, 0.99])
corr_rows.append({"line": line, "plane": plane, "sess_n": sn, "pairs": c.size, "q25": q[0], "median": q[1], "q75": q[2],
"q99": q[3], "frac_above_0.8_%": 100 * np.mean(c > 0.8)})
cr = pd.DataFrame(corr_rows).set_index(["line", "sess_n", "plane"])
with pd.option_context("display.float_format", "{:.4g}".format):
print(cr.to_string())
comp = []
for (line, sn), g in cr.groupby(level=["line", "sess_n"]):
d, s = g.xs("dend", level="plane").iloc[0], g.xs("soma", level="plane").iloc[0]
comp.append({"line": line, "sess_n": sn, "dend_median": d["median"], "soma_q25": s.q25, "soma_q75": s.q75,
"i": s.q25 <= d["median"] <= s.q75, "dend_q99": d.q99, "soma_q99": s.q99, "ii": d.q99 <= 1.1 * s.q99,
"ks_stat": scipy.stats.ks_2samp(pooled[(line, "dend", sn)], pooled[(line, "soma", sn)]).statistic})
comp = pd.DataFrame(comp)
with pd.option_context("display.float_format", "{:.4g}".format):
print(comp.to_string(index=False))
print("(i) holds in", comp.i.sum(), "of 6; (ii) holds in", comp.ii.sum(), "of 6")
pairs q25 median q75 q99 frac_above_0.8_%
line sess_n plane
L2/3 1 dend 460779 0.01972 0.02773 0.036 0.08179 0.02105
2 dend 105228 0.01795 0.0245 0.03172 0.1007 0.02756
3 dend 186094 0.01918 0.02692 0.03533 0.0957 0.004299
1 soma 33277 0.01484 0.0337 0.05675 0.1746 0.01503
2 soma 31682 0.02195 0.04531 0.07191 0.1774 0.009469
3 soma 30776 0.02558 0.04791 0.07344 0.1867 0.009748
L5 1 dend 1824435 0.009622 0.01619 0.02392 0.08058 0.00581
2 dend 640257 0.01112 0.01881 0.0285 0.1185 0.0114
3 dend 624836 0.01111 0.01856 0.028 0.1138 0.006242
1 soma 8961 0.007197 0.02341 0.04338 0.1791 0
2 soma 6481 0.01538 0.03058 0.05013 0.1897 0
3 soma 8143 0.01898 0.04351 0.08251 0.2067 0
line sess_n dend_median soma_q25 soma_q75 i dend_q99 soma_q99 ii ks_stat L2/3 1 0.02773 0.01484 0.05675 True 0.08179 0.1746 True 0.2705 L2/3 2 0.0245 0.02195 0.07191 True 0.1007 0.1774 True 0.4548 L2/3 3 0.02692 0.02558 0.07344 True 0.0957 0.1867 True 0.4364 L5 1 0.01619 0.007197 0.04338 True 0.08058 0.1791 True 0.2641 L5 2 0.01881 0.01538 0.05013 True 0.1185 0.1897 True 0.2871 L5 3 0.01856 0.01898 0.08251 False 0.1138 0.2067 True 0.4325 (i) holds in 5 of 6; (ii) holds in 6 of 6
In [32]:
fig, ax = plt.subplots(1, 2, figsize=(10, 3.5), sharey=True)
bins = np.linspace(-0.3, 1, 66)
for k, line in enumerate(("L2/3", "L5")):
for plane, ls in (("dend", "-"), ("soma", "--")):
for sn in (1, 2, 3):
c = pooled[(line, plane, sn)]
ax[k].hist(c, bins=bins, histtype="step", density=True, ls=ls, label=f"{plane} s{sn}")
print(line, plane, sn, "range", c.min().round(3), "to", c.max().round(3), "| pairs off the axes:", int(((c < -0.3) | (c > 1)).sum()))
ax[k].set_yscale("log"); ax[k].set_title(line); ax[k].set_xlabel("pairwise correlation")
ax[0].legend(fontsize=7); ax[0].set_ylabel("density")
fig.tight_layout()
L2/3 dend 1 range -0.263 to 0.953 | pairs off the axes: 0 L2/3 dend 2 range -0.107 to 0.904 | pairs off the axes: 0 L2/3 dend 3 range -0.433 to 0.886 | pairs off the axes: 1 L2/3 soma 1 range -0.087 to 0.899 | pairs off the axes: 0 L2/3 soma 2 range -0.1 to 0.938 | pairs off the axes: 0 L2/3 soma 3 range -0.091 to 0.897 | pairs off the axes: 0 L5 dend 1 range -0.295 to 0.955 | pairs off the axes: 0 L5 dend 2 range -0.318 to 0.948 | pairs off the axes: 1 L5 dend 3 range -0.531 to 0.919 | pairs off the axes: 4 L5 soma 1 range -0.107 to 0.603 | pairs off the axes: 0 L5 soma 2 range -0.062 to 0.578 | pairs off the axes: 0 L5 soma 3 range -0.13 to 0.618 | pairs off the axes: 0
Result: blocked (methods underspecified). Reading (ii) holds in all six comparisons: the dendritic 99th percentile is below the somatic one in each (0.081 to 0.119 against 0.175 to 0.207). Reading (i) holds in 5 of 6: in L5 session 3 the dendritic median, 0.0186, is just below the somatic 25th percentile, 0.0190. The criterion asks both readings to hold in all six for a reproduction and blocks the claim when they disagree.
high-correlation-pairs¶
About 0.01% of possible pairs of L2/3 dendritic ROIs have correlations above 0.8.
Readings for the denominator: (a) pairs pooled over the L2/3-D sessions 1 to 3 of the analysed mice; (b) mean of per-session fractions; (c) session 1 only (pooled over its mice).
In [33]:
l23d = md_[md_.analysed & (md_.line == "L2/3") & (md_.plane == "dend")]
per_s = pd.DataFrame({s: {"mouse": meta[s]["subject_id"], "sess_n": l23d.set_index("session").loc[s, "sess_n"],
"pairs": roi[s]["corr"].size, "above": int(np.sum(roi[s]["corr"] > 0.8))} for s in l23d.session}).T
per_s["frac_%"] = 100 * per_s.above.astype(float) / per_s.pairs.astype(float)
print(per_s.to_string())
a = 100 * per_s.above.sum() / per_s.pairs.sum()
b = per_s["frac_%"].mean()
s1 = per_s[per_s.sess_n == 1]
c = 100 * s1.above.sum() / s1.pairs.sum()
bymouse = per_s.groupby("mouse")[["pairs", "above"]].sum()
bymouse["frac_%"] = 100 * bymouse.above.astype(float) / bymouse.pairs.astype(float)
print(bymouse.to_string())
for m_out in bymouse.index:
rest = per_s[per_s.mouse != m_out]
print(f"pooled without {m_out}: {100 * rest.above.sum() / rest.pairs.sum():.4f}%")
print(f"(a) pooled {a:.4f}% (b) mean of sessions {b:.4f}% (c) session 1 {c:.4f}% | tolerance 0.005 to 0.015%")
other = cr["frac_above_0.8_%"].groupby(level=["line", "plane"]).mean()
print("mean over sessions of the pooled fraction above 0.8 (%), other planes:", other.round(4).to_dict())
mouse sess_n pairs above frac_%
764704289 413663 1 196878 90 0.045714
765193831 413663 2 60378 26 0.043062
766502238 413663 3 130816 7 0.005351
826187862 433414 1 263901 7 0.002653
826773996 433414 2 44850 3 0.006689
827833392 433414 3 55278 1 0.001809
pairs above frac_%
mouse
413663 388072 123 0.031695
433414 364029 11 0.003022
pooled without 413663: 0.0030%
pooled without 433414: 0.0317%
(a) pooled 0.0178% (b) mean of sessions 0.0175% (c) session 1 0.0211% | tolerance 0.005 to 0.015%
mean over sessions of the pooled fraction above 0.8 (%), other planes: {('L2/3', 'dend'): 0.0176, ('L2/3', 'soma'): 0.0114, ('L5', 'dend'): 0.0078, ('L5', 'soma'): 0.0}
Result: not reproduced. 0.0178% of L2/3-D pairs have correlations above 0.8 when pooled over the six sessions (reading a), 0.0175% as a mean of sessions (b), and 0.0211% in session 1 (c), all above the 0.005 to 0.015% the claim allows. Per session the fraction ranges from 0.0018 to 0.046%. By mouse, 413663 gives 0.032% (123 of the 134 pairs above 0.8) and 433414 0.0030%; without 413663 the pooled fraction is 0.0030%, without 433414 0.032%, both outside the tolerance as well.
fig9-onset-offset-transients¶
Mean dF/F traces show distinct transients aligned with the stimulus onsets and offsets (grey screen to stimulus and back), in every plane and session.
Per ROI, the mean over events of (mean dF/F 0 to 2 s after - mean over 2 s before), for onsets (grey to stimulus) and offsets (stimulus to grey), both stimuli. Readings (i) |mean across ROIs| > 2 SEM, (ii) two-sided Wilcoxon signed-rank across ROIs, p < 0.05. Variants: Gabor onsets and offsets only, and the code's scaled traces.
In [34]:
res9 = []
for key in ("onset_diff", "offset_diff", "onset_diff_gabor", "offset_diff_gabor", "onset_diff_scaled", "offset_diff_scaled"):
for (line, plane, sn), g in md_[md_.analysed].groupby(["line", "plane", "sess_n"]):
dlt = np.concatenate([roi[s][key] for s in g.session])
n, mean, sem, p, wil, sem2 = two_tests(np.zeros_like(dlt), dlt)
res9.append({"measure": key, "group": f"{line}-{'S' if plane == 'soma' else 'D'}", "sess_n": sn, "n_rois": n,
"mean": mean, "sem": sem, "wilcoxon_p": p, "i_2sem": sem2, "ii_wilcoxon": wil})
res9 = pd.DataFrame(res9)
with pd.option_context("display.float_format", "{:.3g}".format):
print(res9[res9.measure.isin(["onset_diff", "offset_diff"])].to_string(index=False))
print(res9.groupby("measure")[["i_2sem", "ii_wilcoxon"]].sum().astype(int).astype(str).add(" of 12").to_string())
measure group sess_n n_rois mean sem wilcoxon_p i_2sem ii_wilcoxon
onset_diff L2/3-D 1 1355 0.051 0.00153 3.99e-193 True True
onset_diff L2/3-D 2 648 0.0517 0.0019 6.29e-93 True True
onset_diff L2/3-D 3 845 0.0647 0.00251 3.51e-131 True True
onset_diff L2/3-S 1 407 0.0487 0.00508 9.33e-46 True True
onset_diff L2/3-S 2 388 0.0601 0.00537 6.15e-55 True True
onset_diff L2/3-S 3 397 0.0862 0.00615 6.86e-64 True True
onset_diff L5-D 1 2807 0.103 0.00283 0 True True
onset_diff L5-D 2 1633 0.0615 0.00248 1.92e-213 True True
onset_diff L5-D 3 1664 0.0691 0.00292 5.03e-194 True True
onset_diff L5-S 1 204 0.0738 0.00754 4.85e-30 True True
onset_diff L5-S 2 183 0.0622 0.00678 8.1e-29 True True
onset_diff L5-S 3 205 0.0726 0.00746 6.49e-29 True True
offset_diff L2/3-D 1 1355 -0.0132 0.000904 3.87e-69 True True
offset_diff L2/3-D 2 648 -0.0248 0.00123 1.39e-69 True True
offset_diff L2/3-D 3 845 -0.0161 0.00121 4.34e-62 True True
offset_diff L2/3-S 1 407 -0.0278 0.00328 2.75e-46 True True
offset_diff L2/3-S 2 388 -0.0227 0.00413 1.02e-38 True True
offset_diff L2/3-S 3 397 -0.0289 0.00425 2.64e-45 True True
offset_diff L5-D 1 2807 -0.0163 0.000993 6.13e-109 True True
offset_diff L5-D 2 1633 -0.0243 0.00108 6.4e-159 True True
offset_diff L5-D 3 1664 -0.0144 0.00119 1.52e-88 True True
offset_diff L5-S 1 204 -0.0348 0.00368 1.17e-29 True True
offset_diff L5-S 2 183 -0.0193 0.00863 3.16e-25 True True
offset_diff L5-S 3 205 -0.0244 0.00603 3.73e-26 True True
i_2sem ii_wilcoxon
measure
offset_diff 12 of 12 12 of 12
offset_diff_gabor 7 of 12 9 of 12
offset_diff_scaled 12 of 12 12 of 12
onset_diff 12 of 12 12 of 12
onset_diff_gabor 12 of 12 12 of 12
onset_diff_scaled 12 of 12 12 of 12
In [35]:
# By mouse (unscaled, both stimuli): mean difference and Wilcoxon p.
bym9 = []
for key in ("onset_diff", "offset_diff"):
for s in sorted(ANALYSED):
m = md_.set_index("session").loc[s]
n, mean, sem, p, wil, sem2 = two_tests(np.zeros(len(roi[s][key])), roi[s][key])
bym9.append({"measure": key, "plane": m.line + "-" + m.plane, "mouse": m.subject_id, "sess_n": m.sess_n, "mean": mean, "p": p})
bym9 = pd.DataFrame(bym9)
with pd.option_context("display.float_format", "{:.2g}".format):
print(bym9.pivot_table(index=["plane", "mouse"], columns=["measure", "sess_n"], values="mean").to_string())
print("mouse-sessions with p < 0.05: onset", int((bym9[bym9.measure == "onset_diff"].p < 0.05).sum()),
"| offset", int((bym9[bym9.measure == "offset_diff"].p < 0.05).sum()), "of", len(ANALYSED))
measure offset_diff onset_diff
sess_n 1 2 3 1 2 3
plane mouse
L2/3-dend 413663 -0.013 -0.033 -0.019 0.064 0.05 0.067
433414 -0.013 -0.016 -0.012 0.04 0.054 0.06
L2/3-soma 408021 -0.029 -0.062 -0.046 0.026 0.069 0.062
411424 -0.041 -0.032 -0.041 0.062 0.033 0.074
440889 -0.022 -0.0057 -0.015 0.053 0.068 0.1
L5-dend 420011 -0.0059 -0.01 -0.0049 0.04 0.07 0.074
433448 -0.0097 -0.014 -0.013 0.14 0.11 0.12
433451 -0.03 -0.031 -0.017 0.058 0.039 0.041
L5-soma 411771 -0.036 -0.0046 -0.0032 0.042 0.073 0.044
418779 -0.045 -0.023 -0.017 0.033 0.055 0.075
433458 -0.032 -0.03 -0.044 0.11 0.056 0.095
mouse-sessions with p < 0.05: onset 33 | offset 33 of 33
In [36]:
# Mean traces (+-SEM across ROIs) around onsets and offsets, per group, unscaled, both stimuli.
tt = np.arange(-ca.PRE, ca.POST) / ca.FPS
fig, ax = plt.subplots(2, 4, figsize=(13, 5), sharex=True)
for k, (line, plane) in enumerate([("L2/3", "dend"), ("L2/3", "soma"), ("L5", "dend"), ("L5", "soma")]):
for r, key in enumerate(("onset_trace", "offset_trace")):
for sn in (1, 2, 3):
ss = md_[md_.analysed & (md_.line == line) & (md_.plane == plane) & (md_.sess_n == sn)].session
tr_ = np.concatenate([roi[s][key] for s in ss])
mu, se = tr_.mean(0), tr_.std(0, ddof=1) / np.sqrt(len(tr_))
ax[r, k].plot(tt, mu, label=f"session {sn}"); ax[r, k].fill_between(tt, mu - se, mu + se, alpha=0.3)
print(line, plane, key, sn, "mean trace range", mu.min().round(4), "to", mu.max().round(4))
ax[r, k].axvline(0, ls="--", c="k"); ax[r, k].set_title(f"{line} {plane} {key.split('_')[0]}", fontsize=9)
ax[0, 0].legend(fontsize=7)
for a in ax[1]: a.set_xlabel("time from event (s)")
for a in ax[:, 0]: a.set_ylabel("dF/F")
fig.tight_layout()
L2/3 dend onset_trace 1 mean trace range -0.026 to 0.1326 L2/3 dend onset_trace 2 mean trace range -0.0299 to 0.1137 L2/3 dend onset_trace 3 mean trace range -0.0455 to 0.1298 L2/3 dend offset_trace 1 mean trace range -0.0075 to 0.0346 L2/3 dend offset_trace 2 mean trace range -0.0159 to 0.0461 L2/3 dend offset_trace 3 mean trace range -0.0135 to 0.0551 L2/3 soma onset_trace 1 mean trace range -0.0244 to 0.107 L2/3 soma onset_trace 2 mean trace range -0.024 to 0.1036
L2/3 soma onset_trace 3 mean trace range -0.031 to 0.1431 L2/3 soma offset_trace 1 mean trace range -0.0138 to 0.041 L2/3 soma offset_trace 2 mean trace range -0.0013 to 0.0431 L2/3 soma offset_trace 3 mean trace range -0.0116 to 0.0538 L5 dend onset_trace 1 mean trace range -0.0321 to 0.219 L5 dend onset_trace 2 mean trace range -0.0313 to 0.1137 L5 dend onset_trace 3 mean trace range -0.0223 to 0.1311 L5 dend offset_trace 1 mean trace range 0.017 to 0.0545 L5 dend offset_trace 2 mean trace range 0.0003 to 0.0538 L5 dend offset_trace 3 mean trace range -0.0042 to 0.0578 L5 soma onset_trace 1 mean trace range -0.0153 to 0.1124 L5 soma onset_trace 2 mean trace range -0.0179 to 0.0944 L5 soma onset_trace 3 mean trace range -0.0237 to 0.1151 L5 soma offset_trace 1 mean trace range -0.0166 to 0.0493 L5 soma offset_trace 2 mean trace range -0.012 to 0.0426 L5 soma offset_trace 3 mean trace range -0.0204 to 0.0483
Result: reproduced. Both readings find a transient at onset and at offset in all 12 groups. The
mean dF/F rises after stimulus onsets (mean differences 0.049 to 0.103 dF/F) and is lower after offsets
than before them (-0.013 to -0.035). With Gabor events only, onsets give a difference in 12 of 12
groups and offsets in 7 (i) and 9 (ii) of 12. The Fig. 9 legends name the Gabor sequence onset and
offset only, while the text names both stimuli; whether the claim should be judged on Gabor events alone
is asked of a reader (questions/000336.yaml). On that variant both readings agree that some groups
lack an offset transient, which the criterion calls partially reproduced.
Summary¶
| Claim | Outcome |
|---|---|
| dataset-holds-paper-sessions | not reproduced: 000336 holds none of the paper's sessions |
| table1-subjects | reproduced (on 000037) |
| table2-sessions | reproduced (on 000037) |
| session-schedule | reproduced (on 000037), under readings (a) and (ii) |
| imaging-planes | partially reproduced (on 000037): 433414 is at 60 um in the files, 75 um in Table 2 |
| roi-counts-per-session | reproduced (on 000037) |
| tracked-rois-per-plane | reproduced (on 000037), fractions under reading (a) |
| fig8-example-matches | reproduced (on 000037) |
| fig8-ordering-variation | reproduced (on 000037), in roi_matching_orderings.ipynb, resting on assumptions |
| stimulus-durations | reproduced (on 000037) |
| inconsistent-event-fraction | reproduced (on 000037), blocks under reading (i), resting on assumptions (durations on the sampling grid, block count as a mean over sessions) |
| fig5-stimulus-responses | reproduced (on 000037) |
| fig6-behaviour | reproduced (on 000037) |
| fig7a-snr-high-stable | reproduced (on 000037); pure Gaussian noise has SNR about 2.06 on this measure |
| fig7b-signal-stable | reproduced (on 000037) |
| fig7c-correlations-similar | blocked, methods underspecified: the criterion's two readings disagree in one of six comparisons |
| high-correlation-pairs | not reproduced (on 000037): about 0.018% of pairs, against 0.01 +- 0.005% |
| fig9-onset-offset-transients | reproduced (on 000037) |
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: the target dandiset, 000336, does not hold the paper's data, so every finding was tested on DANDI:000037 (version 0.230426.0054). Session numbers and QC come from the paper's code table, as the files carry neither. Onset, offset and Fig. 5 measures use unscaled dF/F, as the claims file says, with the code's scaled traces as a variant. The Fig. 5 test and ROI selection, the stability and similarity readings and the theoretical maximum of tracked ROIs are the claims file's readings, since the paper gives no tests. The ordering analysis uses the public matching package, not the paper's modified one, with the overlap constraint read as intersection over union.
What was found about the datasets:
- DANDI:000336 names the paper as a related resource but holds a different OpenScope project (16 other mice, recorded 2022 to 2023); the paper's data are in DANDI:000037.
- In DANDI:000037 the stored
session_start_timecarries a Pacific offset (-07:00 or -08:00) on a clock time that is the UTC start time, so the stored times are 7 to 8 hours late; read as written, 4 of 50 sessions fall on a different date from Table 2. - The three sessions of 433414 give a depth of 60 um, where Table 2 and the paper's code table give 75 um.
- The files carry no QC status or session number, which Table 2 and the paper's code give.
- Session 759872185 (411400, failed QC) has 61134 pupil samples for 127588 imaging frames.
Checks added after the first verification: the subject metadata of all 16 subjects of 000336; the paper's
code table in full beside Table 2; the raw depth strings of 433414; per-session block counts and the runs
furthest from the sampling grid; the Fig. 5 selection counts and thresholds, and a Gabor selection on
inconsistent sequences; a pure-noise control for SNR and a comparison of the helper's SNR, rolling mean and
correlations with the paper's code; and the L2/3-D high-correlation fraction by mouse. In
roi_matching_orderings.ipynb: the package's own run on all orderings of two mice, and the matches lost
only for lacking a direct pairwise link. After the second verification: the correlations checked against
the paper's own correlation step (get_correlation), and the ordering analysis run on all 13 tracked mice
as well as the 11 analysed ones.
In [37]:
print(f"notebook ran in {(time.time() - T0) / 60:.1f} min")
notebook ran in 12.7 min