{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Storing estimated HRFs in snirf files\n", "\n", "This notebook estimates the HRF in a finger-tapping experiment by blockaveraging and then stores the result in a snirf file." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:21.893747Z", "iopub.status.busy": "2024-09-02T16:04:21.893595Z", "iopub.status.idle": "2024-09-02T16:04:23.463127Z", "shell.execute_reply": "2024-09-02T16:04:23.462605Z" } }, "outputs": [], "source": [ "from pathlib import Path\n", "import tempfile\n", "\n", "import matplotlib.pyplot as p\n", "import numpy as np\n", "import xarray as xr\n", "from snirf import Snirf\n", "\n", "import cedalion\n", "import cedalion.datasets\n", "import cedalion.io\n", "import cedalion.nirs\n", "\n", "xr.set_options(display_max_rows=3, display_values_threshold=50, display_expand_data=False)\n", "np.set_printoptions(precision=4)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Load a finger-tapping dataset \n", "\n", "For this demo we load an example finger-tapping recording through `cedalion.datasets.get_fingertapping`. The file contains a single NIRS element with one block of raw amplitude data. " ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:23.465682Z", "iopub.status.busy": "2024-09-02T16:04:23.465324Z", "iopub.status.idle": "2024-09-02T16:04:23.572857Z", "shell.execute_reply": "2024-09-02T16:04:23.572332Z" } }, "outputs": [], "source": [ "rec = cedalion.datasets.get_fingertapping()" ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:23.575149Z", "iopub.status.busy": "2024-09-02T16:04:23.574976Z", "iopub.status.idle": "2024-09-02T16:04:23.578530Z", "shell.execute_reply": "2024-09-02T16:04:23.578081Z" } }, "outputs": [], "source": [ "# Rename events\n", "rec.stim.cd.rename_events( {\n", " \"1.0\" : \"control\",\n", " \"2.0\" : \"Tapping/Left\",\n", " \"3.0\" : \"Tapping/Right\"\n", "})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Calculate concentrations" ] }, { "cell_type": "code", "execution_count": 4, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:23.580489Z", "iopub.status.busy": "2024-09-02T16:04:23.580204Z", "iopub.status.idle": "2024-09-02T16:04:23.636684Z", "shell.execute_reply": "2024-09-02T16:04:23.636157Z" } }, "outputs": [], "source": [ "dpf = xr.DataArray([6, 6], dims=\"wavelength\", coords={\"wavelength\" : rec[\"amp\"].wavelength})\n", "rec[\"od\"] = - np.log( rec[\"amp\"] / rec[\"amp\"].mean(\"time\") )\n", "rec[\"conc\"] = cedalion.nirs.beer_lambert(rec[\"amp\"], rec.geo3d, dpf)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Frequency filtering and splitting into epochs" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:23.639128Z", "iopub.status.busy": "2024-09-02T16:04:23.638951Z", "iopub.status.idle": "2024-09-02T16:04:23.718534Z", "shell.execute_reply": "2024-09-02T16:04:23.717970Z" } }, "outputs": [], "source": [ "rec[\"conc_freqfilt\"] = rec[\"conc\"].cd.freq_filter(fmin=0.02, fmax=0.5, butter_order=4)\n", "\n", "\n", "cf_epochs = rec[\"conc_freqfilt\"].cd.to_epochs(\n", " rec.stim, # stimulus dataframe\n", " [\"Tapping/Left\", \"Tapping/Right\"], # select events\n", " before=5, # seconds before stimulus\n", " after=20, # seconds after stimulus\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Blockaveraging to estimate the HRFs" ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:23.720740Z", "iopub.status.busy": "2024-09-02T16:04:23.720563Z", "iopub.status.idle": "2024-09-02T16:04:23.739441Z", "shell.execute_reply": "2024-09-02T16:04:23.738982Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray 'concentration' (trial_type: 2, chromo: 2, channel: 28,\n",
       "                                   reltime: 196)> Size: 176kB\n",
       "[µM] -0.05399 -0.05393 -0.05398 -0.05408 ... -0.005362 -0.005991 -0.007054\n",
       "Coordinates: (3/6)\n",
       "  * chromo      (chromo) <U3 24B 'HbO' 'HbR'\n",
       "  * channel     (channel) object 224B 'S1D1' 'S1D2' 'S1D3' ... 'S8D8' 'S8D16'\n",
       "    ...          ...\n",
       "  * trial_type  (trial_type) object 16B 'Tapping/Left' 'Tapping/Right'
" ], "text/plain": [ " Size: 176kB\n", "[µM] -0.05399 -0.05393 -0.05398 -0.05408 ... -0.005362 -0.005991 -0.007054\n", "Coordinates: (3/6)\n", " * chromo (chromo) >" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/plain": [ ">" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(s.nirs)\n", "display(s.nirs[0].data)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### Stim DataFrame" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:24.914711Z", "iopub.status.busy": "2024-09-02T16:04:24.914539Z", "iopub.status.idle": "2024-09-02T16:04:24.917797Z", "shell.execute_reply": "2024-09-02T16:04:24.917385Z" } }, "outputs": [ { "data": { "text/plain": [ ">" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "display(s.nirs[0].stim)" ] }, { "cell_type": "code", "execution_count": 12, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:24.919395Z", "iopub.status.busy": "2024-09-02T16:04:24.919254Z", "iopub.status.idle": "2024-09-02T16:04:24.928881Z", "shell.execute_reply": "2024-09-02T16:04:24.928476Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
onsetdurationvaluetrial_type
061.8245.01.0control
187.2965.01.0control
2181.5045.01.0control
3275.7125.01.0control
4544.6405.01.0control
...............
872234.1125.01.0Tapping/Right
882334.8485.01.0Tapping/Right
892362.2405.01.0Tapping/Right
902677.5045.01.0Tapping/Right
912908.0325.01.0Tapping/Right
\n", "

92 rows × 4 columns

\n", "
" ], "text/plain": [ " onset duration value trial_type\n", "0 61.824 5.0 1.0 control\n", "1 87.296 5.0 1.0 control\n", "2 181.504 5.0 1.0 control\n", "3 275.712 5.0 1.0 control\n", "4 544.640 5.0 1.0 control\n", ".. ... ... ... ...\n", "87 2234.112 5.0 1.0 Tapping/Right\n", "88 2334.848 5.0 1.0 Tapping/Right\n", "89 2362.240 5.0 1.0 Tapping/Right\n", "90 2677.504 5.0 1.0 Tapping/Right\n", "91 2908.032 5.0 1.0 Tapping/Right\n", "\n", "[92 rows x 4 columns]" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "cedalion.io.snirf.stim_to_dataframe(s.nirs[0].stim)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "#### MeasurementList DataFrame" ] }, { "cell_type": "code", "execution_count": 13, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:24.930632Z", "iopub.status.busy": "2024-09-02T16:04:24.930453Z", "iopub.status.idle": "2024-09-02T16:04:24.970631Z", "shell.execute_reply": "2024-09-02T16:04:24.970103Z" } }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexwavelengthIndexdataTypedataUnit
01111volt
11121volt
21911volt
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex wavelengthIndex dataType dataUnit\n", "0 1 1 1 1 volt\n", "1 1 1 2 1 volt\n", "2 1 9 1 1 volt" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexwavelengthIndexdataTypedataUnit
5381521volt
548811volt
558821volt
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex wavelengthIndex dataType dataUnit\n", "53 8 15 2 1 volt\n", "54 8 8 1 1 volt\n", "55 8 8 2 1 volt" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexwavelengthIndexdataTypedataUnitdataTypeLabel
011199999dimensionlessdOD
111299999dimensionlessdOD
219199999dimensionlessdOD
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex wavelengthIndex dataType dataUnit \\\n", "0 1 1 1 99999 dimensionless \n", "1 1 1 2 99999 dimensionless \n", "2 1 9 1 99999 dimensionless \n", "\n", " dataTypeLabel \n", "0 dOD \n", "1 dOD \n", "2 dOD " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexwavelengthIndexdataTypedataUnitdataTypeLabel
53815299999dimensionlessdOD
5488199999dimensionlessdOD
5588299999dimensionlessdOD
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex wavelengthIndex dataType dataUnit \\\n", "53 8 15 2 99999 dimensionless \n", "54 8 8 1 99999 dimensionless \n", "55 8 8 2 99999 dimensionless \n", "\n", " dataTypeLabel \n", "53 dOD \n", "54 dOD \n", "55 dOD " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexdataTypedataUnitdataTypeLabel
01199999micromolarHbO
11199999micromolarHbR
21999999micromolarHbO
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex dataType dataUnit dataTypeLabel\n", "0 1 1 99999 micromolar HbO\n", "1 1 1 99999 micromolar HbR\n", "2 1 9 99999 micromolar HbO" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexdataTypedataUnitdataTypeLabel
5381599999micromolarHbR
548899999micromolarHbO
558899999micromolarHbR
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex dataType dataUnit dataTypeLabel\n", "53 8 15 99999 micromolar HbR\n", "54 8 8 99999 micromolar HbO\n", "55 8 8 99999 micromolar HbR" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexdataTypedataUnitdataTypeLabel
01199999micromolarHbO
11199999micromolarHbR
21999999micromolarHbO
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex dataType dataUnit dataTypeLabel\n", "0 1 1 99999 micromolar HbO\n", "1 1 1 99999 micromolar HbR\n", "2 1 9 99999 micromolar HbO" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexdataTypedataUnitdataTypeLabel
5381599999micromolarHbR
548899999micromolarHbO
558899999micromolarHbR
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex dataType dataUnit dataTypeLabel\n", "53 8 15 99999 micromolar HbR\n", "54 8 8 99999 micromolar HbO\n", "55 8 8 99999 micromolar HbR" ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexdataTypedataUnitdataTypeLabeldataTypeIndex
01199999micromolarHRF HbO3
11199999micromolarHRF HbR3
21999999micromolarHRF HbO3
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex dataType dataUnit dataTypeLabel \\\n", "0 1 1 99999 micromolar HRF HbO \n", "1 1 1 99999 micromolar HRF HbR \n", "2 1 9 99999 micromolar HRF HbO \n", "\n", " dataTypeIndex \n", "0 3 \n", "1 3 \n", "2 3 " ] }, "metadata": {}, "output_type": "display_data" }, { "data": { "text/html": [ "
\n", "\n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "
sourceIndexdetectorIndexdataTypedataUnitdataTypeLabeldataTypeIndex
10981599999micromolarHRF HbR4
1108899999micromolarHRF HbO4
1118899999micromolarHRF HbR4
\n", "
" ], "text/plain": [ " sourceIndex detectorIndex dataType dataUnit dataTypeLabel \\\n", "109 8 15 99999 micromolar HRF HbR \n", "110 8 8 99999 micromolar HRF HbO \n", "111 8 8 99999 micromolar HRF HbR \n", "\n", " dataTypeIndex \n", "109 4 \n", "110 4 \n", "111 4 " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "for i_data in range(len(s.nirs[0].data)):\n", " df_ml = cedalion.io.snirf.measurement_list_to_dataframe(\n", " s.nirs[0].data[i_data].measurementList, drop_none=True\n", " )\n", " display(df_ml.head(3))\n", " display(df_ml.tail(3))" ] }, { "cell_type": "code", "execution_count": 14, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:24.972477Z", "iopub.status.busy": "2024-09-02T16:04:24.972177Z", "iopub.status.idle": "2024-09-02T16:04:24.978907Z", "shell.execute_reply": "2024-09-02T16:04:24.978430Z" } }, "outputs": [], "source": [ "s.close()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Read HRFs from snirf file\n", "\n", "Note: read_snirf names the time dimension `time` whereas in `blockaverage` it was called `reltime`. Need to agree on a convention." ] }, { "cell_type": "code", "execution_count": 15, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:24.980774Z", "iopub.status.busy": "2024-09-02T16:04:24.980481Z", "iopub.status.idle": "2024-09-02T16:04:25.444628Z", "shell.execute_reply": "2024-09-02T16:04:25.444132Z" } }, "outputs": [ { "data": { "text/plain": [ "[]" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "hrf_recs = cedalion.io.read_snirf(snirf_fname)\n", "display(hrf_recs)" ] }, { "cell_type": "code", "execution_count": 16, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:25.446461Z", "iopub.status.busy": "2024-09-02T16:04:25.446291Z", "iopub.status.idle": "2024-09-02T16:04:25.457585Z", "shell.execute_reply": "2024-09-02T16:04:25.457175Z" }, "scrolled": false }, "outputs": [ { "data": { "text/html": [ "
\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "\n", "
<xarray.DataArray (channel: 28, chromo: 2, trial_type: 2, time: 196)> Size: 176kB\n",
       "[µM] -0.05399 -0.05393 -0.05398 -0.05408 ... -0.005362 -0.005991 -0.007054\n",
       "Coordinates: (3/7)\n",
       "  * time        (time) float64 2kB -4.992 -4.864 -4.736 ... 19.71 19.84 19.97\n",
       "    samples     (time) int64 2kB 0 1 2 3 4 5 6 7 ... 189 190 191 192 193 194 195\n",
       "    ...          ...\n",
       "  * trial_type  (trial_type) <U13 104B 'Tapping/Left' 'Tapping/Right'\n",
       "Attributes:\n",
       "    data_type_group:  processed HRF concentrations
" ], "text/plain": [ " Size: 176kB\n", "[µM] -0.05399 -0.05393 -0.05398 -0.05408 ... -0.005362 -0.005991 -0.007054\n", "Coordinates: (3/7)\n", " * time (time) float64 2kB -4.992 -4.864 -4.736 ... 19.71 19.84 19.97\n", " samples (time) int64 2kB 0 1 2 3 4 5 6 7 ... 189 190 191 192 193 194 195\n", " ... ...\n", " * trial_type (trial_type) " ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ba = read_blockaverage\n", "\n", "f,ax = p.subplots(5,6, figsize=(12,10))\n", "ax = ax.flatten()\n", "for i_ch, ch in enumerate(ba.channel.values):\n", " for ls, trial_type in zip([\"-\", \"--\"], ba.trial_type): \n", " ax[i_ch].plot(ba.time, ba.sel(chromo=\"HbO\", trial_type=trial_type, channel=ch), \"r\", lw=2, ls=ls)\n", " ax[i_ch].plot(ba.time, ba.sel(chromo=\"HbR\", trial_type=trial_type, channel=ch), \"b\", lw=2, ls=ls)\n", " ax[i_ch].grid(1)\n", " ax[i_ch].set_title(ch)\n", " ax[i_ch].set_ylim(-.3, .6)\n", "\n", "p.tight_layout()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Tidy up" ] }, { "cell_type": "code", "execution_count": 19, "metadata": { "execution": { "iopub.execute_input": "2024-09-02T16:04:27.804738Z", "iopub.status.busy": "2024-09-02T16:04:27.804430Z", "iopub.status.idle": "2024-09-02T16:04:27.815518Z", "shell.execute_reply": "2024-09-02T16:04:27.815072Z" } }, "outputs": [], "source": [ "del tmpdir" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.11.8" } }, "nbformat": 4, "nbformat_minor": 2 }