septum-mec/actions/lfp_speed/data/10_lfp_speed.ipynb

420 lines
20 KiB
Plaintext

{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [],
"source": [
"%load_ext autoreload\n",
"%autoreload 2"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"14:47:12 [I] klustakwik KlustaKwik2 version 0.2.6\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject\n",
" return f(*args, **kwds)\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject\n",
" return f(*args, **kwds)\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/importlib/_bootstrap.py:219: RuntimeWarning: numpy.ufunc size changed, may indicate binary incompatibility. Expected 192 from C header, got 216 from PyObject\n",
" return f(*args, **kwds)\n"
]
}
],
"source": [
"import os\n",
"import pathlib\n",
"import numpy as np\n",
"import matplotlib.pyplot as plt\n",
"from matplotlib import colors\n",
"import seaborn as sns\n",
"import re\n",
"import shutil\n",
"import pandas as pd\n",
"import scipy.stats\n",
"\n",
"import exdir\n",
"import expipe\n",
"from distutils.dir_util import copy_tree\n",
"import septum_mec\n",
"import spatial_maps as sp\n",
"import head_direction.head as head\n",
"import septum_mec.analysis.data_processing as dp\n",
"import septum_mec.analysis.registration\n",
"from septum_mec.analysis.plotting import violinplot, despine\n",
"from spatial_maps.fields import (\n",
" find_peaks, calculate_field_centers, separate_fields_by_laplace, \n",
" map_pass_to_unit_circle, calculate_field_centers, distance_to_edge_function, \n",
" which_field, compute_crossings)\n",
"from phase_precession import cl_corr\n",
"from spike_statistics.core import permutation_resampling\n",
"import matplotlib.mlab as mlab\n",
"import scipy.signal as ss\n",
"from scipy.interpolate import interp1d\n",
"from septum_mec.analysis.plotting import regplot\n",
"from skimage import measure\n",
"from tqdm.notebook import tqdm_notebook as tqdm\n",
"tqdm.pandas()\n",
"\n",
"import pycwt"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"max_speed = 1 # m/s only used for speed score\n",
"min_speed = 0.02 # m/s only used for speed score\n",
"position_sampling_rate = 1000 # for interpolation\n",
"position_low_pass_frequency = 6 # for low pass filtering of position\n",
"\n",
"box_size = [1.0, 1.0]\n",
"bin_size = 0.02\n",
"\n",
"speed_binsize = 0.02\n",
"\n",
"stim_mask = True\n",
"baseline_duration = 600"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [],
"source": [
"data_loader = dp.Data(\n",
" position_sampling_rate=position_sampling_rate, \n",
" position_low_pass_frequency=position_low_pass_frequency,\n",
" box_size=box_size, bin_size=bin_size, \n",
" stim_mask=stim_mask, baseline_duration=baseline_duration\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"project_path = dp.project_path()\n",
"project = expipe.get_project(project_path)\n",
"actions = project.actions\n",
"\n",
"output_path = pathlib.Path(\"output\") / \"lfp-speed\"\n",
"(output_path / \"statistics\").mkdir(exist_ok=True, parents=True)\n",
"(output_path / \"figures\").mkdir(exist_ok=True, parents=True)"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [],
"source": [
"identify_neurons = actions['identify-neurons']\n",
"sessions = pd.read_csv(identify_neurons.data_path('sessions'))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [],
"source": [
"channel_groups = []\n",
"for i, row in sessions.iterrows():\n",
" for ch in range(8):\n",
" row['channel_group'] = ch\n",
" channel_groups.append(row.to_dict())"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [],
"source": [
"sessions = pd.DataFrame(channel_groups)"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [],
"source": [
"def signaltonoise(a, axis=0, ddof=0):\n",
" a = np.asanyarray(a)\n",
" m = a.mean(axis)\n",
" sd = a.std(axis=axis, ddof=ddof)\n",
" return np.where(sd == 0, 0, m / sd)\n",
"\n",
"\n",
"def remove_artifacts(anas, spikes=None, width=500, threshold=2, sampling_rate=None, fillval=0):\n",
" sampling_rate = sampling_rate or anas.sampling_rate.magnitude\n",
" times = np.arange(anas.shape[0]) / sampling_rate\n",
" anas = np.array(anas)\n",
" if anas.ndim == 1:\n",
" anas = np.reshape(anas, (anas.size, 1))\n",
" assert len(times) == anas.shape[0]\n",
" nchan = anas.shape[1]\n",
" if spikes is not None:\n",
" spikes = np.array(spikes)\n",
" for ch in range(nchan):\n",
" idxs, = np.where(abs(anas[:, ch]) > threshold)\n",
" for idx in idxs:\n",
" if spikes is not None:\n",
" t0 = times[idx-width]\n",
" stop = idx+width\n",
" if stop > len(times) - 1:\n",
" stop = len(times) - 1 \n",
" t1 = times[stop]\n",
" mask = (spikes > t0) & (spikes < t1)\n",
" spikes = spikes[~mask]\n",
" anas[idx-width:idx+width, ch] = fillval\n",
" if spikes is not None:\n",
" spikes = spikes[spikes <= times[-1]]\n",
" return anas, times, spikes\n",
" else:\n",
" return anas, times\n",
" \n",
"def find_theta_peak(p, f, f1, f2):\n",
" if np.all(np.isnan(p)):\n",
" return np.nan, np.nan\n",
" mask = (f > f1) & (f < f2)\n",
" p_m = p[mask]\n",
" f_m = f[mask]\n",
" peaks = find_peaks(p_m)\n",
" idx = np.argmax(p_m[peaks])\n",
" return f_m[peaks[idx]], p_m[peaks[idx]]"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [],
"source": [
"def zscore(a):\n",
" return (a - a.mean()) / a.std()\n",
"# return a"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"output = exdir.File(output_path / 'results')\n",
"\n",
"mother = pycwt.Morlet(80)\n",
"NFFT = 2056\n",
"\n",
"def process(row, flim=[6,10]):\n",
" name = row['action'] + '-' + str(row['channel_group'])\n",
" \n",
" lfp = data_loader.lfp(row.action, row.channel_group)\n",
" sample_rate = lfp.sampling_rate.magnitude\n",
" sampling_period = 1 / sample_rate\n",
" x, y, t, speed = map(data_loader.tracking(row.action).get, ['x', 'y', 't', 'v'])\n",
" cleaned_lfp, times = remove_artifacts(lfp)\n",
" speed = interp1d(t, speed, bounds_error=False, fill_value='extrapolate')(times)\n",
" peak_amp = {}\n",
" for i, ch in enumerate(cleaned_lfp.T):\n",
" pxx, freqs = mlab.psd(ch, Fs=lfp.sampling_rate.magnitude, NFFT=4000)\n",
" f, p = find_theta_peak(pxx, freqs, *flim)\n",
" peak_amp[i] = p\n",
"\n",
" theta_channel = max(peak_amp, key=lambda x: peak_amp[x])\n",
" signal = zscore(cleaned_lfp[:,theta_channel])\n",
" \n",
" assert np.array\n",
" \n",
" if name in output:\n",
" return\n",
" \n",
" results = output.require_group(name)\n",
" freqs = np.arange(*flim, .1)\n",
" wave, scales, freqs, coi, fft, fftfreqs = pycwt.cwt(\n",
" signal, sampling_period, freqs=freqs, wavelet=mother)\n",
" \n",
" power = (np.abs(wave)) ** 2\n",
" power /= scales[:, None] #rectify the power spectrum according to the suggestions proposed by Liu et al. (2007)\n",
" \n",
" theta_freq = np.array([freqs[i] for i in np.argmax(power, axis=0)])\n",
" theta_power = np.mean(power, axis=0)\n",
"\n",
" speed_bins = np.arange(min_speed, max_speed + speed_binsize, speed_binsize)\n",
" ia = np.digitize(speed, bins=speed_bins, right=True)\n",
" mean_freq = np.zeros_like(speed_bins)\n",
" mean_power = np.zeros_like(speed_bins)\n",
" for i in range(len(speed_bins)):\n",
" mean_freq[i] = np.mean(theta_freq[ia==i])\n",
" mean_power[i] = np.mean(theta_power[ia==i])\n",
" \n",
" freq_score = np.corrcoef(speed, theta_freq)[0,1]\n",
" power_score = np.corrcoef(speed, theta_power)[0,1]\n",
" \n",
" results.attrs = {\n",
" 'freq_score': float(freq_score),\n",
" 'sample_rate': float(sample_rate),\n",
" 'power_score': float(power_score),\n",
" 'action': row['action'],\n",
" 'channel_group': int(row['channel_group']),\n",
" 'max_speed': max_speed,\n",
" 'min_speed': min_speed,\n",
" 'position_low_pass_frequency': position_low_pass_frequency\n",
" }\n",
" \n",
" results.create_dataset('wavelet_power', data=power)\n",
" results.create_dataset('wavelet_freqs', data=freqs)\n",
" results.create_dataset('theta_freq', data=theta_freq)\n",
" results.create_dataset('theta_power', data=theta_power)\n",
" results.create_dataset('speed', data=speed)\n",
" results.create_dataset('mean_freq', data=mean_freq)\n",
" results.create_dataset('mean_power', data=mean_power)\n",
" results.create_dataset('speed_bins', data=speed_bins)\n"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {
"scrolled": false
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "fab7e86d753341aa9a74ffcd56a128d6",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
"HBox(children=(IntProgress(value=0, max=696), HTML(value='')))"
]
},
"metadata": {},
"output_type": "display_data"
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/fromnumeric.py:3257: RuntimeWarning: Mean of empty slice.\n",
" out=out, **kwargs)\n",
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/numpy/core/_methods.py:161: RuntimeWarning: invalid value encountered in double_scalars\n",
" ret = ret.dtype.type(ret / rcount)\n"
]
},
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n"
]
}
],
"source": [
"sessions.progress_apply(process, axis=1);"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Store results in Expipe action"
]
},
{
"cell_type": "code",
"execution_count": 20,
"metadata": {},
"outputs": [],
"source": [
"action = project.require_action(\"lfp_speed\")"
]
},
{
"cell_type": "code",
"execution_count": 22,
"metadata": {},
"outputs": [
{
"ename": "DistutilsFileError",
"evalue": "could not create '/media/storage/expipe/septum-mec/actions/lfp_speed/data/results.exdir/1849-110319-3-5/attributes.yaml': No such file or directory",
"output_type": "error",
"traceback": [
"\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
"\u001b[0;31mFileNotFoundError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m/usr/lib/python3.6/distutils/file_util.py\u001b[0m in \u001b[0;36m_copy_file_contents\u001b[0;34m(src, dst, buffer_size)\u001b[0m\n\u001b[1;32m 40\u001b[0m \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 41\u001b[0;31m \u001b[0mfdst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdst\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'wb'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 42\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mOSError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mFileNotFoundError\u001b[0m: [Errno 2] No such file or directory: '/media/storage/expipe/septum-mec/actions/lfp_speed/data/results.exdir/1849-110319-3-5/attributes.yaml'",
"\nDuring handling of the above exception, another exception occurred:\n",
"\u001b[0;31mDistutilsFileError\u001b[0m Traceback (most recent call last)",
"\u001b[0;32m<ipython-input-22-150bddb90459>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[1;32m 1\u001b[0m \u001b[0mpower_spectrum_density_action\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m\"results\"\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m\"results.exdir\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 2\u001b[0;31m \u001b[0mcopy_tree\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0moutput_path\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mstr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maction\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdata_path\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
"\u001b[0;32m/usr/lib/python3.6/distutils/dir_util.py\u001b[0m in \u001b[0;36mcopy_tree\u001b[0;34m(src, dst, preserve_mode, preserve_times, preserve_symlinks, update, verbose, dry_run)\u001b[0m\n\u001b[1;32m 170\u001b[0m copy_tree(src_name, dst_name, preserve_mode,\n\u001b[1;32m 171\u001b[0m \u001b[0mpreserve_times\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpreserve_symlinks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupdate\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 172\u001b[0;31m verbose=verbose, dry_run=dry_run))\n\u001b[0m\u001b[1;32m 173\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 174\u001b[0m copy_file(src_name, dst_name, preserve_mode,\n",
"\u001b[0;32m/usr/lib/python3.6/distutils/dir_util.py\u001b[0m in \u001b[0;36mcopy_tree\u001b[0;34m(src, dst, preserve_mode, preserve_times, preserve_symlinks, update, verbose, dry_run)\u001b[0m\n\u001b[1;32m 170\u001b[0m copy_tree(src_name, dst_name, preserve_mode,\n\u001b[1;32m 171\u001b[0m \u001b[0mpreserve_times\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpreserve_symlinks\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupdate\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 172\u001b[0;31m verbose=verbose, dry_run=dry_run))\n\u001b[0m\u001b[1;32m 173\u001b[0m \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 174\u001b[0m copy_file(src_name, dst_name, preserve_mode,\n",
"\u001b[0;32m/usr/lib/python3.6/distutils/dir_util.py\u001b[0m in \u001b[0;36mcopy_tree\u001b[0;34m(src, dst, preserve_mode, preserve_times, preserve_symlinks, update, verbose, dry_run)\u001b[0m\n\u001b[1;32m 174\u001b[0m copy_file(src_name, dst_name, preserve_mode,\n\u001b[1;32m 175\u001b[0m \u001b[0mpreserve_times\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mupdate\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mverbose\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mverbose\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 176\u001b[0;31m dry_run=dry_run)\n\u001b[0m\u001b[1;32m 177\u001b[0m \u001b[0moutputs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdst_name\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 178\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/lib/python3.6/distutils/file_util.py\u001b[0m in \u001b[0;36mcopy_file\u001b[0;34m(src, dst, preserve_mode, preserve_times, update, link, verbose, dry_run)\u001b[0m\n\u001b[1;32m 149\u001b[0m \u001b[0;31m# Otherwise (non-Mac, not linking), copy the file contents and\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 150\u001b[0m \u001b[0;31m# (optionally) copy the times and mode.\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 151\u001b[0;31m \u001b[0m_copy_file_contents\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdst\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m 152\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0mpreserve_mode\u001b[0m \u001b[0;32mor\u001b[0m \u001b[0mpreserve_times\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 153\u001b[0m \u001b[0mst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mos\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msrc\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;32m/usr/lib/python3.6/distutils/file_util.py\u001b[0m in \u001b[0;36m_copy_file_contents\u001b[0;34m(src, dst, buffer_size)\u001b[0m\n\u001b[1;32m 42\u001b[0m \u001b[0;32mexcept\u001b[0m \u001b[0mOSError\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0me\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 43\u001b[0m raise DistutilsFileError(\n\u001b[0;32m---> 44\u001b[0;31m \"could not create '%s': %s\" % (dst, e.strerror))\n\u001b[0m\u001b[1;32m 45\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m 46\u001b[0m \u001b[0;32mwhile\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
"\u001b[0;31mDistutilsFileError\u001b[0m: could not create '/media/storage/expipe/septum-mec/actions/lfp_speed/data/results.exdir/1849-110319-3-5/attributes.yaml': No such file or directory"
]
}
],
"source": [
"power_spectrum_density_action.data[\"results\"] = \"results.exdir\"\n",
"copy_tree(output_path, str(action.data_path()))"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"septum_mec.analysis.registration.store_notebook(action, \"10_lfp_speed.ipynb\")"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"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.6.8"
}
},
"nbformat": 4,
"nbformat_minor": 2
}