diff --git a/code/ProblemSet-DataDetective.ipynb b/code/ProblemSet-DataDetective.ipynb new file mode 100644 index 0000000..a2af7bd --- /dev/null +++ b/code/ProblemSet-DataDetective.ipynb @@ -0,0 +1,1902 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "835d3761", + "metadata": {}, + "source": [ + "

SWDB Problem Set: Becoming a Data Detective

\n", + "

From someone else's figure to your own analysis

\n", + "

Works with any SWDB dataset — bring the one your chosen figure came from.

" + ] + }, + { + "cell_type": "markdown", + "id": "b76af5ad", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

How this problem set works

\n", + "\n", + "This morning you explored a dataset and made figures. Those figures are now posted on Slack.\n", + "\n", + "**Your starting point is one of your classmates' figures.** Pick any figure from the channel, along\n", + "with the dataset it came from — ideally one you did *not* work on this morning.\n", + "\n", + "| Part | Task |\n", + "| --- | --- |\n", + "| 1 | Load their dataset and find the pieces the figure needs |\n", + "| 2 | Reproduce the figure, and interrogate what it shows |\n", + "| 3 | Align activity to event onsets: raster and PSTH |\n", + "| 4 | Signal and noise correlations, and whether to trust them |\n", + "\n", + "You already have the data-access skills for Part 1 from this morning's tutorial. This problem set is\n", + "about what comes after loading: **shaping data, and checking whether the result means anything.**\n", + "\n", + "**Deliverable:** a short README naming the figure and dataset you chose, the decisions you made at\n", + "each step, and an honest assessment of what your numbers do and do not support.\n", + "\n", + "Every dataset is different, and the notebook does not know which one you picked. The code\n", + "cells are prompts, not templates — you write what goes in them, using the access patterns from\n", + "this morning. Only a few things are given: the imports, and two helper functions from the tutorial.\n", + "\n", + "The differences you will run into are not cosmetic. Across the datasets in this workshop:\n", + "\n", + "- **Recording modality** — a continuous calcium signal in some, discrete spike times in\n", + " others. Spikes need binning before anything here applies.\n", + "- **Sampling rate** — from a few Hz to tens of kHz, which sets what timing you can resolve.\n", + "- **Number of neurons** — tens to thousands, which changes what is tractable in one pass.\n", + "- **Stimulus structure** — many conditions with few repeats, few conditions with many, or no\n", + " sensory stimulus at all.\n", + "- **What was recorded alongside** — running, licking, pupil, reward; some datasets have all of\n", + " it, some none.\n", + "- **Where things live in the file** — container and column names differ, and so does which\n", + " container holds the trial table.\n", + "\n", + "None of that is written on the outside of the file. **You have to look.** Part of each prompt is\n", + "deciding whether the analysis it asks for even applies to your dataset — and saying so when it\n", + "does not.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "a9965fde", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Taking it slow: Analysis step by step

\n", + "\n", + "You can now generate an analysis faster than you can check one. Ask an LLM for a correlation matrix\n", + "and you will have one in thirty seconds, beautifully formatted, with a colorbar.\n", + "\n", + "The problem is that a result computed on four trials can look exactly like a result computed on four\n", + "hundred. A bug can look exactly like a finding. A correlation computed in a window where nothing\n", + "happened can look exactly like a real effect.\n", + "\n", + "So the questions to keep asking are:\n", + "\n", + "- **What is actually in this file?** Not what you assume — what is there.\n", + "- **Does this dataset support the question I am asking?**\n", + "- **How is the data being transformed?** Plot the data after each step.\n", + "- **What would make this result wrong?** Name it before you see the answer.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "0dfe350b", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "b9da37e5", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Part 1: Load the dataset and find the pieces you need

\n", + "\n", + "Same access pattern as this morning: find your dataset's mount under /data, locate a\n", + "session's NWB file, then dot and bracket notation into the containers.\n", + "\n", + "**Your classmate's figure tells you what to look for.** Before you open anything, list the pieces the\n", + "figure needs — neural activity, plus whatever else it plots: a behavioral trace, epoch\n", + "boundaries, trial times, stimulus identity.\n", + "\n", + "Then find each one, and note the ones that turn out not to exist. **A piece being absent is a\n", + "finding about the dataset, not a failure.** Some datasets have no running wheel, no pupil camera, no\n", + "visual stimulus at all. You will build the figure from what is there.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50d8a203", + "metadata": {}, + "outputs": [], + "source": [ + "import os\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "import matplotlib.pyplot as plt\n", + "import pynwb\n", + "from scipy import stats\n", + "\n", + "pd.set_option('display.width', 200)\n", + "pd.set_option('display.max_columns', 30)\n", + "\n", + "data_dir = '/data'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a5cd672", + "metadata": {}, + "outputs": [], + "source": [ + "# List the datasets mounted under /data." + ] + }, + { + "cell_type": "markdown", + "id": "34f16e31", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Start from the metadata table, not the file tree. Each dataset has a metadata CSV in\n", + "/code/metadata/ — one row per session, with subject, session type, date and the\n", + "asset name. Read that first and choose a session from it, because the filename alone will not tell you\n", + "which imaging stage or task condition you are looking at.\n", + "\n", + "Then build the path: the NWB lives inside that dataset's mount under /data/.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "76df3e44", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: read your dataset's metadata CSV from /code/metadata and look at what it\n", + "# offers: how many subjects, how many session types, how many sessions each." + ] + }, + { + "cell_type": "markdown", + "id": "cd09d433", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Which session does your classmate's figure come from? Use the table to find it\n", + "— subject, session type, date — and say what you filtered on.\n", + "\n", + "Look at what the table offers before you filter. How many subjects, how many session types, how many\n", + "sessions each? That inventory is the first thing you know about the dataset.\n", + "\n", + "**Then ask what kind of neurons you are recording from.** This is not a detail — it decides\n", + "what your population average means. Check the transgenic line, the virus, and any other metadata\n", + "describing what was labeled (`nwb.subject.genotype`, the imaging plane's `indicator`, the session\n", + "metadata table).\n", + "\n", + "- **Imaging.** You see only the cells expressing the calcium indicator. A pan-excitatory driver\n", + " gives you a very different population from an interneuron-specific one, and \"population activity\"\n", + " in each case means something different.\n", + "- **Electrophysiology.** A probe records whatever is near it, so the recording is not cell-type\n", + " specific by default. But a line or virus may still be present for **optotagging** — light\n", + " activation used to identify a targeted cell type among the recorded units. If so, there may be a\n", + " column marking which units were tagged.\n", + "\n", + "Write down what is labeled in your session, and say what population your averages are actually\n", + "averaging over.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3995769b", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: filter the table to the session behind your figure, take one row as\n", + "# `session`, and say what you filtered on." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "79a72f42", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: examine the column values of the session you selected\n", + "# what is the session type, genotype, the targeted structure, etc. " + ] + }, + { + "cell_type": "markdown", + "id": "15085858", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Now build the path. The metadata table's name column is normally the session's\n", + "folder name inside the mount, so you can go straight there rather than searching. Inside that folder\n", + "sits one NWB store — either a single .nwb file (HDF5) or a directory\n", + "(zarr). List the folder, keep the entry with nwb in its name, and check you got exactly\n", + "one before continuing.\n", + "\n", + "
\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d81da35f", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: set `dataset_dir` to the mount holding your dataset (one of the names\n", + "# printed above), then join it with your session's folder name to get `session_dir`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3cb59b2", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: list `session_dir`, keep entries with 'nwb' in the name, and assert you\n", + "# found exactly one. Watch for sidecar files that also contain 'nwb'. Join the\n", + "# match onto `session_dir` to get `nwb_path`.\n" + ] + }, + { + "cell_type": "markdown", + "id": "9922f8c6", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Did you get exactly one match? More than one usually means several processing\n", + "generations of the same session are attached — check which you picked. Zero means the session\n", + "in the table is not mounted in this capsule, which is worth knowing before you debug anything else.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4fd1343a", + "metadata": {}, + "outputs": [], + "source": [ + "# Open the NWB file. Name it `nwb`." + ] + }, + { + "cell_type": "markdown", + "id": "32bd96d7", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Find the data the figure needs

\n", + "\n", + "A handful of containers hold almost everything. Which one holds what **varies by dataset**, so list\n", + "them all before you index into any of them.\n", + "\n", + "| container | commonly holds |\n", + "| --- | --- |\n", + "| `processing` | processed neural activity — in some datasets also behavior |\n", + "| `intervals` | epoch tables, trial tables, stimulus presentation tables |\n", + "| `stimulus` | stimulus templates — but in some datasets, the trial tables too |\n", + "| `acquisition` | raw acquired signals |\n", + "| `events` | discrete behavioral and stimulus events, in some datasets |\n", + "\n", + "Row three is not hypothetical: some datasets put their trial tables in `stimulus` and leave\n", + "`intervals` holding only epochs. If you look in one container, find nothing, and conclude the data\n", + "is missing, you will be wrong. **Print them all.**\n", + "\n", + "The `events` row needs its own warning. It is optional — plenty of files do not have one, and\n", + "`nwb.processing` will not reveal it either way, because it is reached by its own accessor\n", + "(`nwb.events`, or `nwb.get_all_events()` for a single table across all event types). When it *is*\n", + "present it holds **behavioral and stimulus events** — licks, rewards, stimulus changes —\n", + "each a timestamped row with an `event_type` column. It does **not** hold neural events. Where a file\n", + "has no events table, the same information is usually in a `processing` behavior module or implicit\n", + "in columns of the trials table.\n", + "\n", + "

“Events” means two different things

\n", + "\n", + "The word is overloaded in NWB, and the two meanings live in different places.\n", + "\n", + "1. Neural events — inside a `processing` plane. A plane usually holds several\n", + "representations of the same neurons: raw fluorescence, neuropil-corrected, dF/F, and often events.\n", + "Events are the output of running deconvolution on dF/F — an attempt to recover the\n", + "discrete firing that produced the slow calcium signal. Stored as an array with the same shape and\n", + "same timestamps as dF/F, but mostly zeros: nonzero only where an event was detected, the\n", + "value carrying its inferred magnitude. Treat the nonzero samples as spike-like events, not as a\n", + "continuous trace. The name is not standardised — one dataset calls it events,\n", + "another event_timeseries, and some have none at all and give you only dF/F.\n", + "\n", + "2. Behavioral / task events — a separate table. Discrete, timestamped occurrences during\n", + "the session: licks, rewards, stimulus changes. These may sit in an events table reached through\n", + "`nwb.events` or `nwb.get_all_events()`, in a `processing` behavior module, or be implicit in columns\n", + "of the trials table. Unlike neural events, these are measured, not inferred.\n", + "\n", + "A container is not always visible from the top level, so print the interfaces inside each processing\n", + "module too — and remember `nwb.processing` will not show you an events table reached by its own\n", + "accessor.\n", + "\n", + "
\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "765a84d8", + "metadata": {}, + "outputs": [], + "source": [ + "# What is in this file? Print the containers before you index into any of them:\n", + "# processing, intervals, acquisition, stimulus. Then look INSIDE each processing\n", + "# module -- the listing above only gives you the module names. Check for an events\n", + "# table too; it has its own accessor and will not appear in any of the four.\n" + ] + }, + { + "cell_type": "markdown", + "id": "41eb1f54", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Is your dataset continuous or spiking? This is the first fork in the road, and it\n", + "changes what \"activity\" even means.\n", + "\n", + "Continuous (calcium imaging, LFP): a `(n_timepoints, n_cells)` array already exists in the\n", + "file. Find it and you are done.\n", + "\n", + "Spiking (Neuropixels, sorted electrophysiology): there is no such array. Each unit carries its\n", + "own list of spike times, usually in a `units` table, and you must bin them yourself —\n", + "choose a bin width, count spikes per bin, divide by the width to get a rate in spikes/s. Everything\n", + "downstream then works the same way.\n", + "\n", + "Two decisions come with spiking data, and neither has a default:\n", + "\n", + "- Which units. Spike sorting produces more units than you should analyze. There will be\n", + " quality-control columns (`is_qc_pass`, `firing_rate`, `presence_ratio`, `snr`) and often an\n", + " anatomical label. Select on them explicitly and say what you selected — a session can drop\n", + " from thousands of units to dozens, and the ones you drop change your answer.\n", + "- Bin width. Too wide blurs the response; too narrow leaves mostly-empty bins and noisy\n", + " single-trial estimates. Try a few and see how much your answer moves.\n", + "\n", + "
\n",
+    "bin_width = 0.010                                    # seconds -- your decision\n",
+    "edges = np.arange(0, t_end + bin_width, bin_width)\n",
+    "counts, _ = np.histogram(one_unit_spike_times, bins=edges)\n",
+    "rate = counts / bin_width                            # spikes/s\n",
+    "bin_centres = edges[:-1] + bin_width / 2\n",
+    "
\n", + "\n", + "Sparse binned spikes behave like a deconvolved calcium trace: sharper in time, and noisy per\n", + "trial.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "1317510b", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Two things to check as you pull out the activity trace.\n", + "\n", + "Timestamps. Some datasets store an explicit `timestamps` array; others store a sampling\n", + "`rate` and a `starting_time`, and you reconstruct the times yourself. Everything downstream needs\n", + "real times in seconds, so check which you have — `series.timestamps` is `None` when the file\n", + "uses a rate.\n", + "\n", + "Lazy loading. NWB data objects do not load until you index them. That is what lets you open a\n", + "50 GB file instantly, but it means `data.std()` may fail where `np.std(data)` works. Convert\n", + "with `np.asarray()` once you know the array is small enough to hold, or slice first.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce7e3b99", + "metadata": {}, + "outputs": [], + "source": [ + "# Find the neural activity in this file and pull out two things:\n", + "# `dff` -- the (n_timepoints, n_cells) trace array\n", + "# `ts` -- the matching times in seconds\n", + "# Not every dataset stores a timestamps array; see the note above.\n", + "# Print the shapes, the frame rate, and the session duration as a sanity check." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ace8196", + "metadata": {}, + "outputs": [], + "source": [ + "# Pull out the other pieces your chosen figure needs -- stimulus/trial tables,\n", + "# and any behavioral traces the dataset has. Tables become DataFrames with\n", + "# .to_dataframe(); timeseries have .data and .timestamps.\n", + "# Print the shape of each, and note anything that turns out not to exist." + ] + }, + { + "cell_type": "markdown", + "id": "583526cb", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Quality control: which cells or units belong in the analysis?

\n", + "\n", + "Segmentation and spike sorting are automated, and both over-produce. An ophys plane contains ROIs the\n", + "classifier thinks are not cell bodies; a sorted probe contains units that drift, that are barely\n", + "above noise, or that are two neurons merged. The activity matrix you just loaded usually contains\n", + "all of them.\n", + "\n", + "Pipelines record their own verdicts. For imaging they live on the ROI table beside the masks; for\n", + "electrophysiology, on the units table. The columns differ by pipeline and by dataset — boolean\n", + "flags, continuous probabilities, morphology metrics, contamination estimates — so there is no\n", + "list to memorise. Print the columns and see what your dataset offers.\n", + "\n", + "Filtering is not automatically the right move, and the criteria are yours to justify. But\n", + "inheriting the unfiltered set by default is a decision you made without noticing, and it is the\n", + "kind that never appears in a methods section.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71bdf349", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: find the per-cell quality table for your dataset -- a plane segmentation\n", + "# for imaging, nwb.units for electrophysiology -- and print its scalar columns:\n", + "# which are flags, which are continuous scores, and what each would exclude." + ] + }, + { + "cell_type": "markdown", + "id": "ad6b13e7", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Does your dataset carry per-cell or per-unit quality metrics? Report what the\n", + "columns are, how many entries each flag would exclude, and whether the activity matrix is already\n", + "filtered or contains everything.\n", + "\n", + "Then decide. Whatever you choose, **state the criterion and the count you dropped** — that\n", + "sentence belongs in your methods.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5ddd96e", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: apply your QC criterion. Check the table length matches the activity\n", + "# matrix first, apply the SAME mask to every per-cell array you loaded, and print\n", + "# how many cells you dropped." + ] + }, + { + "cell_type": "markdown", + "id": "33843692", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Plotting a long recording. A whole session at a fine sampling rate can be hundreds of\n", + "thousands of points — slow to draw and impossible to read. Plot a slice instead, but choose the\n", + "slice from the data rather than picking a round number: an arbitrary window can easily contain no\n", + "activity at all, and an empty panel looks identical to a broken one.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0ce932b5", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot one cell's trace against time, as a sanity check on what you loaded.\n", + "# Choose the cell deliberately rather than taking index 0, and say how you chose.\n", + "# Watch out for cells that are entirely NaN.\n", + "# If the recording is long, plot a slice -- and pick the slice from the data." + ] + }, + { + "cell_type": "markdown", + "id": "057cd94e", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Look at that trace for a few seconds before moving on. Is anything about it\n", + "surprising? Would you have noticed if you had skipped straight to the analysis?\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "093087fc", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Set up the main variables for this dataset

\n", + "\n", + "Point these names at the equivalent pieces of your own NWB file. Later sections reference them,\n", + "so getting them right here saves repeating yourself — but edit anything you like as you go.\n", + "This is your notebook now.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f1fbee26", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the main variables for this dataset. Later sections use these names,\n", + "# so getting them right here saves repeating yourself -- but edit anything you\n", + "# like as you go.\n", + "activity = ... # (n_timepoints, n_cells)\n", + "timestamps = ... # (n_timepoints,) in seconds\n", + "events = ... # one row per event / trial / presentation\n", + "\n", + "# A second activity representation, if your dataset has one (deconvolved\n", + "# events, spike estimates). Set to None if it does not.\n", + "activity_events = ...\n" + ] + }, + { + "cell_type": "markdown", + "id": "b8d7adc3", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Print the columns of your stimulus table. Which describe *what was\n", + "presented*, which describe *what the animal did*, and which are bookkeeping?\n", + "\n", + "Note any column whose meaning you cannot guess — that is a databook lookup for your README.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f50525e", + "metadata": {}, + "outputs": [], + "source": [ + "# Print the columns of your event table, then look at the first few rows." + ] + }, + { + "cell_type": "markdown", + "id": "b6e8635f", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "15dbba5d", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Part 2: Reproduce the figure, and interrogate what it shows

\n", + "\n", + "You have your classmate's figure. You do not have their code, and you may not have a caption either.\n", + "\n", + "Before you write anything, write down what you think the figure shows. One or two sentences,\n", + "in your notebook, as a claim someone could disagree with: \"activity is higher during X than during\n", + "Y\", \"the response is larger on this trial type\", \"these two signals rise together.\"\n", + "\n", + "Two reasons this comes first. It commits you to an interpretation before the data can talk you into\n", + "one — and it converts a picture into something you can actually test. A figure cannot be right\n", + "or wrong. A claim can.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "10a7584d", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Write your claim about the figure you picked, in the cell below, before you\n", + "write any code.\n", + "\n", + "Be specific enough to be wrong. \"There is neural activity\" is not a claim; \"population activity is\n", + "higher in the second half of the session\" is.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "855780a7", + "metadata": {}, + "source": [ + "_Your claim:_\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "afa5dc32", + "metadata": {}, + "outputs": [], + "source": [ + "# For each timeseries your figure uses: the typical sampling interval, how many\n", + "# gaps are much larger than it, and how much of the session has no data." + ] + }, + { + "cell_type": "markdown", + "id": "672b7f36", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Now rebuild it

\n", + "\n", + "Get the pieces the figure needs and plot them. You will not match it exactly — different\n", + "smoothing, different colors, a different subset of cells — and that is fine. What matters is\n", + "that the structure you see is the same structure they saw.\n", + "\n", + "If you cannot rebuild some element because the dataset does not contain it, note that and rebuild\n", + "what you can.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2418ea7", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute whatever your chosen figure shows.\n", + "# For a population_rate average: the mean across cells at each timepoint.\n", + "# Name it `population_rate`, and check its shape before you plot." + ] + }, + { + "cell_type": "markdown", + "id": "9ee42d5b", + "metadata": {}, + "source": [ + "
\n", + "\n", + "To shade the epochs we need their start and stop times. Where epochs live varies by dataset:\n", + "sometimes an `epoch_name` column on the stimulus table, sometimes a separate epochs table.\n", + "\n", + "**Check that the column you group by actually varies.** If it takes one value, you will get a single\n", + "block spanning the session — a figure that looks fine and is wrong.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c3698fbc", + "metadata": {}, + "outputs": [], + "source": [ + "# Does your event table carry an epoch column? If so, how many distinct\n", + "# values does it actually take?" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed0fe478", + "metadata": {}, + "outputs": [], + "source": [ + "# Build a table of epoch boundaries with one row per epoch, indexed by name\n", + "# and sorted in time. Name it `epochs`, with `start_time` and `stop_time`\n", + "# columns -- the shading helper below expects those." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47b46b60", + "metadata": {}, + "outputs": [], + "source": [ + "# Shade each epoch a different color -- same helper as the tutorial\n", + "colors = dict(zip(epochs.index, plt.cm.Pastel1.colors))\n", + "\n", + "\n", + "def shade_epoch_blocks(ax):\n", + " \"\"\"Shade each epoch on `ax`, one colour per epoch label.\n", + "\n", + " Epochs are the coarse structure of the session -- which stimulus block or\n", + " task phase was running. Shading them behind a trace shows at a glance\n", + " whether a change in activity lines up with a change in what was happening.\n", + " \"\"\"\n", + " for label, row in epochs.iterrows():\n", + " # zorder=0 keeps the shading BEHIND the data; alpha so the trace on top\n", + " # stays readable. label= puts each epoch in the legend once.\n", + " ax.axvspan(row.start_time, row.stop_time, color=colors[label],\n", + " alpha=0.5, zorder=0, label=label)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2a9253b", + "metadata": {}, + "outputs": [], + "source": [ + "# Build your version of the figure: the traces stacked on a shared time axis,\n", + "# with the epochs shaded (use `shade_epoch_blocks`). Include only the streams your\n", + "# dataset actually has." + ] + }, + { + "cell_type": "markdown", + "id": "509a4b43", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Now test the claim you wrote above — do not eyeball it.\n", + "\n", + "Turn your sentence into a number you can check. If it compares epochs, compute the mean in each one,\n", + "alongside how long each epoch lasted, when in the session it happened, and what the animal was doing.\n", + "If it compares something else, compute the equivalent.\n", + "\n", + "Before you look: **what would make this comparison unfair?** Write your answer down first, then see\n", + "whether the table bears it out.\n", + "\n", + "Then go back and mark your claim as supported, contradicted, or untestable with this data. All three\n", + "are legitimate outcomes and all three belong in your README.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "326af115", + "metadata": {}, + "outputs": [], + "source": [ + "# For each epoch compute the mean activity, and alongside it the things that\n", + "# could confound the comparison: how long the epoch lasted, how many samples\n", + "# that is, when in the session it happened, and what the animal was doing.\n", + "# Build a DataFrame with one row per epoch." + ] + }, + { + "cell_type": "markdown", + "id": "72c9c230", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "3ac7d615", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Part 3: Align activity to event onsets

\n", + "\n", + "The session overview shows everything at once, which means it shows very little. To see a response\n", + "you need to **align** activity to the times when something happened, and look across repeats.\n", + "\n", + "\"Something happened\" need not be a visual stimulus. It might be a sound, an optogenetic pulse, a\n", + "reward, a lick, or the start of a trial. Anything with a repeatable onset time works the same way\n", + "— and the rest of this notebook says \"event\" rather than \"stimulus\" for that reason.\n", + "\n", + "This morning's tutorial averaged across presentations. Here we look at what the average hides.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7fa6db01", + "metadata": {}, + "outputs": [], + "source": [ + "# How many times was each chosen_condition presented?" + ] + }, + { + "cell_type": "markdown", + "id": "0f119c7c", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Are all of these events the same kind of event?

\n", + "\n", + "An event table usually contains rows that are **not equivalent trials**. Depending on the dataset\n", + "that might be first versus repeated presentations, rewarded versus unrewarded trials, different\n", + "stimulus families, trials the animal responded to versus ignored, blocks recorded before and after\n", + "a manipulation, or blank and omitted entries that are not events at all.\n", + "\n", + "This matters before you align anything, for two reasons:\n", + "\n", + "- **Response magnitude can differ several-fold between trial types.** Averaging them together dilutes\n", + " the response toward whichever type is most numerous — which is often the weakest one.\n", + "- **Trial types differ in what else is happening.** Reward, licking, and arousal ride along with some\n", + " trial types and not others, so a difference you attribute to the stimulus may not be about the\n", + " stimulus.\n", + "\n", + "Find the columns in your table that distinguish trial types, and count them.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f02accfa", + "metadata": {}, + "outputs": [], + "source": [ + "# Which columns in your event table distinguish different KINDS of trial?\n", + "# Find them and count the rows of each kind." + ] + }, + { + "cell_type": "markdown", + "id": "f17c032d", + "metadata": {}, + "source": [ + "
\n", + "\n", + "To compare them we need to cut a window of data around each onset. Same\n", + "`align_to_event_times` helper as this morning's tutorial.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3e418107", + "metadata": {}, + "outputs": [], + "source": [ + "def align_to_event_times(data, timestamps, event_times, pre=0.5, post=1.5):\n", + " \"\"\"Cut a window of data around each event time.\n", + "\n", + " data : array with time along the first axis\n", + " timestamps : time of each row of data, in seconds\n", + " event_times : times to align to, in seconds\n", + " pre, post : seconds before and after each event\n", + "\n", + " Returns (aligned_windows, window_time_axis) where the time axis is in\n", + " seconds relative to the event, and there is one window per usable_cells event.\n", + " \"\"\"\n", + " # Sampling interval. Median, not mean: one gap in the recording would\n", + " # inflate a mean and silently shrink every window.\n", + " dt = np.median(np.diff(timestamps))\n", + "\n", + " # Convert the requested seconds into a number of samples. int() truncates,\n", + " # so a window that is not a whole number of samples comes out slightly\n", + " # short -- check this if you need exact window edges.\n", + " n_pre, n_post = int(pre / dt), int(post / dt)\n", + "\n", + " aligned_windows = []\n", + " for event_time in event_times:\n", + " # Index of the first sample AT OR AFTER the event. side='left' returns\n", + " # the insertion point, so timestamps[i] >= event_time always.\n", + " #\n", + " # Do NOT round to the nearest sample: that pulls roughly half the\n", + " # trials one sample EARLIER than the event, which smears the onset and\n", + " # can make a real response look like it starts before the stimulus.\n", + " # Landing just after is honest -- the bias is one-directional and at\n", + " # most one sample.\n", + " i = np.searchsorted(timestamps, event_time, side='left')\n", + "\n", + " # Skip events too close to either end of the recording to fill a whole\n", + " # window. This drops trials SILENTLY, so compare\n", + " # aligned_windows.shape[0] against len(event_times) afterwards.\n", + " if i - n_pre >= 0 and i + n_post <= len(timestamps):\n", + " # Slice is n_pre + n_post samples long. Index n_pre within the\n", + " # window is the first sample at/after the event, i.e. t = 0.\n", + " aligned_windows.append(data[i - n_pre:i + n_post])\n", + "\n", + " # Time axis in seconds relative to the event. Starts at -n_pre*dt, which\n", + " # can be slightly later than -pre because of the truncation above.\n", + " window_time_axis = np.arange(-n_pre, n_post) * dt\n", + " return np.array(aligned_windows), window_time_axis\n" + ] + }, + { + "cell_type": "markdown", + "id": "518b20dd", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Pick two trial types from your table and align the population average to\n", + "each separately, then plot both on the same axes.\n", + "\n", + "Write down your prediction first: do you expect a difference, and how large?\n", + "\n", + "Then choose which type to carry forward, and one condition within it. Name the things below, because\n", + "the rest of Part 3 refers to them:\n", + "\n", + "| name | what it holds |\n", + "| --- | --- |\n", + "| `stimulus_onset_times` | onset times of ALL trials of your chosen type |\n", + "| `onset_times` | onset times of the one condition you picked |\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b428d085", + "metadata": {}, + "outputs": [], + "source": [ + "# Pick two kinds of trial and compare them: align the population_rate average to\n", + "# each set of onset_times separately and plot both on the same axes.\n", + "# Subtract each window's own pre-onset baseline before averaging.\n", + "# Name the two onset arrays `first_type_onset_times` and `second_type_onset_times`.\n", + "\n", + "# Note on the baseline: exclude the sample immediately before onset. With\n", + "# binned or sampled data that sample can straddle the event, so including\n", + "# it puts part of the response into the baseline and shrinks what you\n", + "# measure. `window_time_axis < -bin_width` rather than `< 0`.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "929d6d88", + "metadata": {}, + "outputs": [], + "source": [ + "# Choose the trial type to carry forward, and one chosen_condition within it.\n", + "# Name them:\n", + "# `stimulus_onset_times` -- onset_times of ALL trials of that type\n", + "# `onset_times` -- onset_times of the one chosen_condition you picked" + ] + }, + { + "cell_type": "markdown", + "id": "5bfe8592", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Which cell or unit to look at?

\n", + "\n", + "Whatever your dataset calls them — ROIs in an imaging plane, sorted units on a probe —\n", + "taking the first one in the table is an arbitrary choice you did not disclose. Ranking by how strongly\n", + "they respond is a *different* undisclosed choice unless you say so. Pick deliberately and write down\n", + "how you picked.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "1b842877", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Which signal do you align? Most datasets ship more than one representation of the\n", + "same activity, and the choice is yours — but it is a choice, and it changes what the figures\n", + "show.\n", + "\n", + "\n", + "\n", + "\n", + "\n", + "
ΔF/F (imaging)Continuous fluorescence. Carries the indicator's rise and\n", + "decay, so a brief response is smeared forward by hundreds of milliseconds, and slow drift shared\n", + "across the field of view inflates correlations between any two cells. Every timepoint has a\n", + "value.
Deconvolved events (imaging)An estimate of when the cell actually fired, with\n", + "the indicator kinetics removed. Temporally tighter, and mostly exact zeros — so single-trial\n", + "estimates are much noisier even though the trial average looks cleaner.
Spike times (electrophysiology)Discrete times, no continuous trace at all. You\n", + "choose a bin width to get a matrix, and that width is a real analysis decision: too fine and every\n", + "bin is empty, too coarse and you lose the timing you came for.
\n", + "\n", + "None of these is the correct one. A question about response latency or duration is badly served\n", + "by ΔF/F; a question needing a reliable per-trial number is badly served by a sparse signal. Pick\n", + "one, say why, and if you have time run the analysis twice and compare — that comparison is\n", + "usually more informative than either result alone.\n", + "\n", + "Set the choice in one place so switching it is a one-line edit rather than a rewrite.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a17dafe", + "metadata": {}, + "outputs": [], + "source": [ + "# EDIT: which signal to align for the figures below -- the deconvolved events if\n", + "# your dataset has them, otherwise the continuous trace. Set `aligned_signal` and\n", + "# a `signal_label` string for the axis labels, and say why you chose it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e796db6", + "metadata": {}, + "outputs": [], + "source": [ + "# Choose an example cell to look at, and say how you chose it.\n", + "# Set `pre` and `post` (seconds before/after onset) and name the cell `example_roi`.\n", + "# Re-derive the index against the CURRENT activity matrix -- an index from an\n", + "# earlier cell may not refer to the same neuron." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e2b89cb", + "metadata": {}, + "outputs": [], + "source": [ + "# Cut a window around every onset for your example cell, using\n", + "# `align_to_event_times`. Name the results `aligned_windows` and `t`, and print the shape." + ] + }, + { + "cell_type": "markdown", + "id": "b71b2ed3", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Compare the number of windows you got back against the number of onsets you\n", + "asked for. Are they the same?\n", + "\n", + "If not, read the helper again and work out where the missing trials went — then decide whether\n", + "losing them matters for your analysis.\n", + "\n", + "This is worth doing every time you call something that returns one row per trial. A function that\n", + "quietly returns fewer rows than you gave it will not raise an error; it will just make your\n", + "n smaller than you think it is.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f0d53ee", + "metadata": {}, + "outputs": [], + "source": [ + "# Compare the number of onset_times you asked for against the number of aligned_windows\n", + "# that came back." + ] + }, + { + "cell_type": "markdown", + "id": "2306d820", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Raster and PSTH

\n", + "\n", + "The raster shows every trial; the PSTH is their average. Plot them together so you can see what the\n", + "average discards.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b19071f", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the raster and the PSTH side by side: every trial as a heatmap, and\n", + "# the trial average with a measure of spread." + ] + }, + { + "cell_type": "markdown", + "id": "37b0809f", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Now do it for every cell and plot the result as a heatmap, sorted by\n", + "response magnitude. How many cells respond at all?\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "242adec6", + "metadata": {}, + "outputs": [], + "source": [ + "# Do it for every cell: build a (n_cells, n_timepoints) array of\n", + "# trial-averaged responses. Name it `responses`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2ac363c7", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot `responses` as a heatmap, sorted by response magnitude.\n", + "# Then plot it again with each cell's own pre-onset baseline subtracted,\n", + "# and compare the two panels." + ] + }, + { + "cell_type": "markdown", + "id": "1430b5f2", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

The same analysis on a different signal

\n", + "\n", + "Skip this section if your dataset has only one representation of activity. A probe recording\n", + "gives you spike times and nothing else — there is no second signal to compare against, and\n", + "saying so in your write-up is the correct answer here, not a gap.\n", + "\n", + "If you do have two — a continuous trace and a deconvolved estimate, most commonly — they\n", + "are not interchangeable, and running the same analysis on both is the cheapest way to find out how\n", + "much your conclusion depends on that choice.\n", + "\n", + "Check what your dataset has before assuming. List the interfaces in the processing container\n", + "and see whether a second per-cell timeseries is there at all.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f50f4424", + "metadata": {}, + "outputs": [], + "source": [ + "# Clean tiny float noise to exact zeros so \"fraction exactly zero\" means what\n", + "# it says, then compare the two representations: shapes, sparsity, shared clock." + ] + }, + { + "cell_type": "markdown", + "id": "248dc786", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** If your dataset has two activity representations, align both to the same\n", + "onsets and plot the trial-averaged population response side by side. What differs — the\n", + "duration, the shape, the size relative to baseline?\n", + "\n", + "If it has only one, note that in your README and move on.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "798195b0", + "metadata": {}, + "outputs": [], + "source": [ + "# Build a list of the activity representations you have, then plot the aligned\n", + "# population_rate response for each. If you only have one, say so and move on." + ] + }, + { + "cell_type": "markdown", + "id": "36c28085", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "0fc19301", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Part 4: Signal and noise correlations

\n", + "\n", + "

First, the math

\n", + "\n", + "The Pearson correlation between two variables $x$ and $y$ is\n", + "\n", + "$$ r = \\frac{\\sum_i (x_i - \\bar{x})(y_i - \\bar{y})}\n", + " {\\sqrt{\\sum_i (x_i - \\bar{x})^2}\\;\\sqrt{\\sum_i (y_i - \\bar{y})^2}} $$\n", + "\n", + "In words:\n", + "\n", + "1. **Center** each variable by subtracting its mean.\n", + "2. **Multiply** the centered values pointwise and sum — large and positive when they vary\n", + " together, negative when oppositely, near zero when unrelated.\n", + "3. **Normalize** by each variable's spread, forcing the result between -1 and +1.\n", + "\n", + "Compute it once by hand before running it thousands of times.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e58de86e", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute the correlation between two cells' traces BY HAND, in three steps:\n", + "# 1. centre each variable (subtract its mean)\n", + "# 2. multiply the centred values pointwise and sum\n", + "# 3. normalise by the spread of each\n", + "# Then check your answer against np.corrcoef." + ] + }, + { + "cell_type": "markdown", + "id": "eb57af3d", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Two consequences that matter for everything below:\n", + "\n", + "- $r$ says nothing about response **size**, only whether two things move together.\n", + "- $r$ is computed over a set of paired observations, and **how many observations you have determines\n", + " how noisy $r$ is** — but the value itself gives you no clue how many there were.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "34039ac5", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** What does a given value of $r$ look like? Simulate pairs with known\n", + "correlations and plot them.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5aaed716", + "metadata": {}, + "outputs": [], + "source": [ + "# Simulate pairs of variables with known correlations (try 0, 0.2, 0.5, 0.9)\n", + "# and plot each as a scatter, titled with its measured r." + ] + }, + { + "cell_type": "markdown", + "id": "99087112", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Two reasons neurons are correlated

\n", + "\n", + "- **Signal correlation.** Do they respond similarly *across conditions*? Correlate the two neurons'\n", + " tuning curves — their average response to each condition.\n", + "- **Noise correlation.** When the *same* condition repeats, do they fluctuate together around their\n", + " own averages? Subtract each condition's mean and correlate the residuals.\n", + "\n", + "A \"condition\" is whatever your event table repeats: an image, a grating direction, a tone, a\n", + "photostimulation target, a task context. All that matters is that it recurs enough times to average\n", + "over.\n", + "\n", + "Same data, different thing averaged over:\n", + "\n", + "| | what is correlated | one observation is |\n", + "| --- | --- | --- |\n", + "| signal | condition means | one condition |\n", + "| noise | within-condition residuals | one trial |\n", + "\n", + "That last column matters more than anything else in this notebook.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "a7bf3465", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Step 1: choose which events to use

\n", + "\n", + "Not every event is comparable to every other. Decide which subset is a fair comparison and write down\n", + "why.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef61874b", + "metadata": {}, + "outputs": [], + "source": [ + "# Restrict to comparable events, and write down WHY -- this choice belongs in\n", + "# your methods. Name the results:\n", + "# `all_onset_times` -- the onset times you keep\n", + "# `labels` -- the chosen_condition label for each of those onset_times" + ] + }, + { + "cell_type": "markdown", + "id": "0bf2b4f2", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Step 2: one number per trial per neuron

\n", + "\n", + "We need a `(n_trials, n_cells)` matrix. Average each aligned window over a response window, and\n", + "subtract a **baseline** from just before onset — otherwise each trial's \"response\" includes\n", + "wherever the cell happened to be sitting beforehand, and those levels drift together across the\n", + "population from bleaching, arousal, and movement.\n", + "\n", + "Choosing the two windows is dataset-specific. The response window should cover the response\n", + "your Part 3 plot showed — look at it rather than copying a number from here, since a calcium\n", + "signal and a spike rate need very different windows. The baseline window should sit in the gap\n", + "before onset, and must **exclude any stimulation artifact**: with optogenetics or electrical\n", + "stimulation the frames around the pulse can be unusable, so leave a margin on both sides.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "636791ff", + "metadata": {}, + "outputs": [], + "source": [ + "# Build the (n_trials, n_cells) response matrix `trial_response_matrix`: for each\n", + "# trial, mean activity in a response window minus a pre-onset baseline.\n", + "#\n", + "# BUILD IT IN STEPS, in separate cells, checking as you go. The solutions\n", + "# notebooks do it this way for a reason -- a shape printed at the end of a loop\n", + "# tells you almost nothing about whether the arithmetic inside was right.\n", + "# a) pick the two windows, and print how many SAMPLES each one holds\n", + "# b) do one trial, one cell, by hand -- print the actual values you average\n", + "# c) one trial, all cells -- check the row length equals n_cells\n", + "# d) loop over trials -- count how many rows came out incomplete\n", + "# e) drop those rows from the matrix AND from the labels with the SAME mask,\n", + "# then assert the two lengths match\n", + "#\n", + "# Note on the baseline: exclude the sample immediately before onset. With\n", + "# binned or sampled data that sample can straddle the event, so including it\n", + "# puts part of the response into the baseline and shrinks what you measure --\n", + "# `window_time_axis < -bin_width` rather than `< 0`. Then check how many\n", + "# samples are actually left: a \"baseline\" of one sample is not an average.\n" + ] + }, + { + "cell_type": "markdown", + "id": "84ab91ca", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Step 3: tuning curves — look before correlating

\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8992eaed", + "metadata": {}, + "outputs": [], + "source": [ + "# Build the condition_mean_response curves: the per-chosen_condition mean response of each cell.\n", + "# Name the chosen_condition list `conditions` and the array `condition_mean_response`,\n", + "# shaped (n_conditions, n_cells)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e71ccdb", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the condition_mean_response curves before correlating anything: a few cells as lines,\n", + "# and all cells as a heatmap." + ] + }, + { + "cell_type": "markdown", + "id": "6dfbf7ad", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** How many numbers make up one neuron's tuning curve?\n", + "\n", + "That is how many paired observations each signal correlation gets. Write it down.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ff2650b", + "metadata": {}, + "outputs": [], + "source": [ + "# How many numbers make up one condition_mean_response curve, and how many trials are available\n", + "# for the noise correlations? Print both." + ] + }, + { + "cell_type": "markdown", + "id": "d3fe7dda", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Step 4: residuals — look before correlating

\n", + "\n", + "Subtract **each condition's own mean**, not the grand mean. Subtracting the grand mean would leave\n", + "the differences between conditions in the residuals, making your \"noise\" correlation partly a signal\n", + "correlation.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "415fe1c4", + "metadata": {}, + "outputs": [], + "source": [ + "# Build the residuals: subtract each chosen_condition's OWN mean from its trials.\n", + "# Name the array `residuals`, and check that its overall mean is ~0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dde4cee7", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the raw responses and the residuals for one cell, side by side, so you\n", + "# can see what subtracting the chosen_condition means removed." + ] + }, + { + "cell_type": "markdown", + "id": "fc28dc8b", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Step 5: correlate

\n", + "\n", + "`np.corrcoef` correlates **rows**, so transpose to get cells rather than trials. Getting this\n", + "backwards produces a plausible matrix of entirely the wrong thing — check the output shape\n", + "against the number of cells.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "75007877", + "metadata": {}, + "outputs": [], + "source": [ + "# Compute the two correlation matrices, `signal_corr_matrix` and `noise_corr_matrix`.\n", + "# Watch the orientation: np.corrcoef correlates ROWS.\n", + "# Take each pair once with np.triu_indices -- name the index `pairs`, and the\n", + "# extracted values `signal_values` and `noise_values`.\n", + "# Print the mean of each, with how many observations went into it." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f0c7503f", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the two matrices side by side, plus signal against noise correlation\n", + "# for every pair. Scale each matrix to its own range so neither saturates." + ] + }, + { + "cell_type": "markdown", + "id": "b426fa0f", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Is this result trustworthy?

\n", + "\n", + "Every number so far is a point estimate with no error bar. The single most useful check: **would you\n", + "get the same answer with half the data?**\n", + "\n", + "Split trials in half at random, compute the correlations on each half separately, and correlate the\n", + "two halves' answers. Split **within each condition** so both halves see every condition.\n", + "\n", + "Three outcomes, and all three are informative:\n", + "\n", + "- **One high, one low** — trust the high one, and say why the other is not trustworthy.\n", + "- **Both high** — you have enough data for both; proceed.\n", + "- **Both near zero** — report that. It usually means the condition variable you chose does not\n", + " organise these neurons' responses, however well-balanced it looked in the inventory. That is a\n", + " real result about your dataset, and it is a better README than a matrix you cannot defend.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "650615b8", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Before running it — which do you expect to be more reliable, signal or\n", + "noise correlations? Look back at the observation counts you wrote down.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b183a126", + "metadata": {}, + "outputs": [], + "source": [ + "def signal_and_noise_correlations(responses, labels):\n", + " \"\"\"Signal and noise correlation matrices from a set of trials.\n", + "\n", + " responses : (n_trials, n_cells) one response value per trial per cell\n", + " labels : (n_trials,) which condition each trial belongs to\n", + "\n", + " Signal correlation = do two cells prefer the same conditions?\n", + " Noise correlation = do two cells co-vary trial to trial WITHIN a\n", + " condition, once the condition mean is removed?\n", + " \"\"\"\n", + " conditions = np.unique(labels)\n", + "\n", + " # TUNING: one row per condition, holding that condition's mean response for\n", + " # every cell. Averaging over trials is what removes trial-to-trial noise\n", + " # and leaves the stimulus preference -- the \"signal\".\n", + " condition_means = np.vstack([responses[labels == c].mean(axis=0)\n", + " for c in conditions])\n", + "\n", + " # RESIDUALS: each trial minus its own condition's mean. What remains is\n", + " # everything the condition does NOT explain -- the \"noise\". Subtracting the\n", + " # condition mean is essential: skip it and the condition structure leaks\n", + " # into the noise matrix and inflates it.\n", + " residuals = responses.astype(float).copy()\n", + " for c in conditions:\n", + " in_condition = labels == c\n", + " residuals[in_condition] -= responses[in_condition].mean(axis=0)\n", + "\n", + " # .T because np.corrcoef correlates ROWS: we want cell-by-cell matrices,\n", + " # and cells are the columns of both arrays.\n", + " #\n", + " # Note the very different sample sizes feeding these two matrices: signal\n", + " # is estimated from len(conditions) numbers per cell, noise from\n", + " # len(labels) trials. That asymmetry is why they differ so much in\n", + " # reliability even though both render as equally convincing heatmaps.\n", + " return np.corrcoef(condition_means.T), np.corrcoef(residuals.T)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d9e41a9f", + "metadata": {}, + "outputs": [], + "source": [ + "# Split-half split_half_reliability. Split the trials in half WITHIN each chosen_condition,\n", + "# compute the correlation matrices on each half with `signal_and_noise_correlations`, and\n", + "# correlate the two halves' answers (scipy.stats.spearmanr on `pairs`).\n", + "# Repeat ~10 times; report the mean and spread for signal and for noise." + ] + }, + { + "cell_type": "markdown", + "id": "b060a99a", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Does the signal you chose change the answer?

\n", + "\n", + "Everything so far used one representation of activity. If your dataset provides a second one, repeat\n", + "the whole chain on it and compare the numbers that matter. If it provides only one, note that and\n", + "move on.\n", + "\n", + "To repeat the chain you need the response-matrix construction as a reusable function rather than a\n", + "one-off block — so wrap it, the same way you wrapped the correlations.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7b3b1433", + "metadata": {}, + "outputs": [], + "source": [ + "def trial_by_cell_responses(A):\n", + " \"\"\"Build the (n_trials, n_cells) baseline-subtracted response matrix.\n", + "\n", + " One number per trial per cell: mean activity in the response window minus\n", + " mean activity in the baseline window. Every correlation below is computed\n", + " from this matrix, so both window choices propagate into every later result.\n", + " \"\"\"\n", + " response_rows = []\n", + " for t0 in all_onset_times:\n", + " # Boolean masks selecting the samples in each window for this trial.\n", + " # >= start and < end so the two windows never share a sample.\n", + " in_response = (timestamps >= t0 + response_window[0]) & (timestamps < t0 + response_window[1])\n", + " in_baseline = (timestamps >= t0 + baseline_window[0]) & (timestamps < t0 + baseline_window[1])\n", + "\n", + " # nanmean, not mean: a single all-NaN cell would otherwise propagate\n", + " # NaN across the whole row and silently cost you every trial.\n", + " # A trial at the very start of the recording can have an empty\n", + " # baseline window -- fill it with NaN and drop it below.\n", + " response_rows.append(np.nanmean(A[in_response], axis=0) - np.nanmean(A[in_baseline], axis=0)\n", + " if in_response.sum() and in_baseline.sum()\n", + " else np.full(A.shape[1], np.nan))\n", + "\n", + " responses = np.array(response_rows)\n", + "\n", + " # Drop trials with any missing cell. Report the count if it is not zero:\n", + " # trials vanishing here is exactly the kind of silent loss to check for.\n", + " return responses[~np.isnan(responses).any(axis=1)]\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c096bf9e", + "metadata": {}, + "outputs": [], + "source": [ + "# Run the whole chain on each activity representation your dataset has, and\n", + "# compare: mean signal and noise correlation, their split-half reliabilities,\n", + "# and what fraction of the single-trial responses are exactly zero." + ] + }, + { + "cell_type": "markdown", + "id": "4331f937", + "metadata": {}, + "source": [ + "
\n", + "\n", + "One column may not be the whole condition.\n", + "\n", + "A column can look like a clean condition variable — many levels, perfectly balanced —\n", + "while the stimulus varied in some other way at the same time. Two trials sharing that column's\n", + "value are then not repeats of the same thing, and averaging them together destroys the tuning you\n", + "were trying to measure.\n", + "\n", + "Receptive-field mapping is the classic case: orientation is balanced, but the stimulus also moves\n", + "around the screen, so \"144 repeats of 45°\" is really a handful of repeats at each of many\n", + "positions. The same trap appears whenever a design crosses two factors and you only notice one.\n", + "\n", + "Check for it by asking what else varies across the trials you just called identical. Group by your\n", + "condition column, look at the other columns within a group, and see whether they are constant. If\n", + "they are not, either restrict to one level of the other factor, or make the condition the\n", + "combination of both.\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "4b539b97", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Exercise:** Signal correlations need a condition that repeats. Does your dataset have\n", + "one?\n", + "\n", + "Inventory the candidate columns: how many distinct values, how many repeats, how balanced.\n", + "\n", + "Then answer **two separate questions**, because they can disagree:\n", + "\n", + "1. **Is the analysis possible?** Does some column have enough conditions with enough repeats?\n", + "2. **Is it meaningful?** Does that column label something you would expect neurons to be tuned\n", + " *to*, in a way that a correlation across condition means would capture?\n", + "\n", + "A column can pass the first test and fail the second. State a verdict on both, and check it against\n", + "your reliability numbers.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c66a5bf3", + "metadata": {}, + "outputs": [], + "source": [ + "# Inventory the candidate chosen_condition columns in your event table: for each,\n", + "# how many distinct values, the repeats of the least and most common, and how\n", + "# balanced. Then state your verdict on both questions above." + ] + }, + { + "cell_type": "markdown", + "id": "e3e56e03", + "metadata": {}, + "source": [ + "---" + ] + }, + { + "cell_type": "markdown", + "id": "c633f716", + "metadata": {}, + "source": [ + "
\n", + "\n", + "

Summary

\n", + "\n", + "

The process

\n", + "\n", + "1. **Find out what is in the file** before analyzing it — and check that the dataset supports\n", + " your question. Sometimes the answer is no.\n", + "2. **Plot the data after each transformation.** Single trials before averages; tuning curves before\n", + " correlations.\n", + "3. **Name every decision.** Event subset, condition column, response window, baseline. Each is a\n", + " fork, and each belongs in your methods.\n", + "4. **Try to break your own result.** Split the data in half and see if the answer survives.\n", + "5. **Let the dataset answer back.** If the check says your result is noise, or the dataset has no\n", + " variable that supports your question, that is the finding. Report it rather than reaching for the\n", + " analysis you planned to run.\n", + "\n", + "

Traps this notebook demonstrated

\n", + "\n", + "| trap | how you catch it |\n", + "| --- | --- |\n", + "| A result from few observations looks like one from many | split-half reliability |\n", + "| A well-balanced condition variable that means nothing | reliability, not the inventory |\n", + "| A condition column that hides a second varying factor | group by it, check what else moves |\n", + "| Analyzing units that should have been dropped | select on quality columns, and say so |\n", + "| A helper function silently drops data | compare output shape to input |\n", + "| A column exists but carries no information | check that it actually varies |\n", + "| One bad trial turns every cell's score into NaN | count your NaNs; use `nanmean` |\n", + "| Epoch comparisons confounded with time and behavior | check durations, order, behavior |\n", + "| An example cell chosen to look good | state your selection rule |\n", + "| Data looks absent but is stored elsewhere | look in every container first |\n", + "| An index from an earlier cell after reshaping the data | re-derive indices, never carry them |\n", + "\n", + "

Why this matters

\n", + "\n", + "You can generate an analysis faster than you can validate one. The only defense is to know your data\n", + "well enough that a wrong answer looks wrong to **you** — because it will not look wrong to the\n", + "code, and it will not look wrong on the plot.\n", + "\n", + "
" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "name": "python", + "version": "3.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/code/metadata/visual_learning_metadata.ipynb b/code/metadata/visual_learning_metadata.ipynb new file mode 100644 index 0000000..566cc02 --- /dev/null +++ b/code/metadata/visual_learning_metadata.ipynb @@ -0,0 +1,1694 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "d561dd1b", + "metadata": {}, + "source": [ + "# Visual Learning session metadata\n", + "\n", + "Builds `visual_learning_session_metadata.csv` — one row per session for the six\n", + "Visual Learning mice." + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "0a6da650", + "metadata": {}, + "outputs": [], + "source": [ + "import re\n", + "import time\n", + "from datetime import datetime\n", + "\n", + "import numpy as np\n", + "import pandas as pd\n", + "\n", + "pd.set_option('display.width', 220)\n", + "pd.set_option('display.max_columns', 40)" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "fa59053b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "https://api.allenneuraldynamics.org/v2/metadata_index/data_assets\n" + ] + } + ], + "source": [ + "from aind_data_access_api.document_db import MetadataDbClient\n", + "\n", + "API_GATEWAY_HOST = \"api.allenneuraldynamics.org\"\n", + "OUTPUT_DIR = '/data/metadata'\n", + "DATABASE = 'metadata_index'\n", + "COLLECTION = 'data_assets'\n", + "\n", + "docdb_api_client = MetadataDbClient(\n", + " host=API_GATEWAY_HOST,\n", + " version=\"v2\",\n", + " database=DATABASE,\n", + " collection=COLLECTION,\n", + ")\n", + "print(docdb_api_client._base_url)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "713baaa7", + "metadata": {}, + "outputs": [], + "source": [ + "# The cohort is defined by subject: five different project_name values are\n", + "# interleaved across the same six mice.\n", + "VISUAL_LEARNING_MICE = ['782149', '790322', '788406', '800792', '800995', '804363']\n", + "\n", + "# Processed asset names end in _processed__