spike lfp relations
This commit is contained in:
Mikkel Elle Lepperød
parent
c1e2713401
commit
8ac5e2d59c
|
|
@ -0,0 +1,5 @@
|
|||
registered: '2019-10-15T17:17:48'
|
||||
data:
|
||||
results: results.csv
|
||||
notebook: 10-calculate-stimulus-spike-lfp-response.ipynb
|
||||
html: 10-calculate-stimulus-spike-lfp-response.html
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,467 @@
|
|||
{
|
||||
"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": [
|
||||
"16:56:09 [I] klustakwik KlustaKwik2 version 0.2.6\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"import matplotlib.pyplot as plt\n",
|
||||
"%matplotlib inline\n",
|
||||
"import spatial_maps as sp\n",
|
||||
"import septum_mec.analysis.data_processing as dp\n",
|
||||
"import septum_mec.analysis.registration\n",
|
||||
"import expipe\n",
|
||||
"import os\n",
|
||||
"import pathlib\n",
|
||||
"import numpy as np\n",
|
||||
"import exdir\n",
|
||||
"import pandas as pd\n",
|
||||
"import optogenetics as og\n",
|
||||
"import quantities as pq\n",
|
||||
"import shutil\n",
|
||||
"from distutils.dir_util import copy_tree\n",
|
||||
"import elephant as el\n",
|
||||
"import neo\n",
|
||||
"from scipy.signal import find_peaks\n",
|
||||
"from scipy.interpolate import interp1d\n",
|
||||
"from matplotlib import mlab\n",
|
||||
"\n",
|
||||
"from tqdm import tqdm_notebook as tqdm\n",
|
||||
"from tqdm._tqdm_notebook import tqdm_notebook\n",
|
||||
"tqdm_notebook.pandas()"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 3,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"data_loader = dp.Data()\n",
|
||||
"actions = data_loader.actions\n",
|
||||
"project = data_loader.project"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 4,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"output = pathlib.Path('output/stimulus-spike-lfp-response')\n",
|
||||
"(output / 'figures').mkdir(parents=True, exist_ok=True)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 5,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"identify_neurons = actions['identify-neurons']\n",
|
||||
"# sessions = pd.read_csv(identify_neurons.data_path('sessions'))\n",
|
||||
"units = pd.read_csv(identify_neurons.data_path('units'))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 6,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def get_lim(action_id):\n",
|
||||
" stim_times = data_loader.stim_times(action_id)\n",
|
||||
" if stim_times is None:\n",
|
||||
" return [0, np.inf]\n",
|
||||
" stim_times = np.array(stim_times)\n",
|
||||
" return [stim_times.min(), stim_times.max()]\n",
|
||||
"\n",
|
||||
"def get_mask(lfp, lim):\n",
|
||||
" return (lfp.times >= lim[0]) & (lfp.times <= lim[1])\n",
|
||||
"\n",
|
||||
"def zscore(a):\n",
|
||||
" return (a - a.mean(0)) / a.std(0)\n",
|
||||
"\n",
|
||||
"def compute_stim_freq(action_id):\n",
|
||||
" stim_times = data_loader.stim_times(action_id)\n",
|
||||
" if stim_times is None:\n",
|
||||
" return np.nan\n",
|
||||
" stim_times = np.array(stim_times)\n",
|
||||
" return 1 / np.mean(np.diff(stim_times))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 7,
|
||||
"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 clean_lfp(anas, width=500, threshold=2):\n",
|
||||
" anas = np.array(anas)\n",
|
||||
"\n",
|
||||
" for ch in range(anas.shape[1]):\n",
|
||||
" idxs, = np.where(abs(anas[:, ch]) > threshold)\n",
|
||||
" for idx in idxs:\n",
|
||||
" anas[idx-width:idx+width, ch] = 0 # TODO AR model prediction\n",
|
||||
" return anas"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": 8,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def compute_energy(p, f, f1, f2):\n",
|
||||
" if np.isnan(f1) or np.all(np.isnan(p)):\n",
|
||||
" return np.nan\n",
|
||||
" mask = (f > f1) & (f < f2)\n",
|
||||
" df = f[1] - f[0]\n",
|
||||
" return np.sum(p[mask]) * df"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"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": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def compute_half_width(p, f, m_p, m_f):\n",
|
||||
" if np.isnan(m_p):\n",
|
||||
" return np.nan, np.nan\n",
|
||||
" m_p_half = m_p / 2\n",
|
||||
" half_p = p - m_p_half\n",
|
||||
" idx_f = np.where(f <= m_f)[0].max()\n",
|
||||
" idxs_p1, = np.where(np.diff(half_p[:idx_f + 1] > 0) == 1)\n",
|
||||
" if len(idxs_p1) == 0:\n",
|
||||
" return np.nan, np.nan\n",
|
||||
" m1 = idxs_p1.max()\n",
|
||||
" idxs_p2, = np.where(np.diff(half_p[idx_f:] > 0) == 1)\n",
|
||||
" m2 = idxs_p2.min() + idx_f\n",
|
||||
" assert p[m1] < m_p_half < p[m1+1], (p[m1], m_p_half, p[m1+1])\n",
|
||||
" assert p[m2] > m_p_half > p[m2+1], (p[m2], m_p_half, p[m2+1])\n",
|
||||
" \n",
|
||||
" f1 = interp1d([half_p[m1], half_p[m1 + 1]], [f[m1], f[m1 + 1]])(0)\n",
|
||||
" f2 = interp1d([half_p[m2], half_p[m2 + 1]], [f[m2], f[m2 + 1]])(0)\n",
|
||||
" return f1, f2"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def compute_stim_peak(p, f, s_f):\n",
|
||||
" if np.isnan(s_f):\n",
|
||||
" return np.nan\n",
|
||||
" return interp1d(f, p)(s_f)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def compute_spike_lfp(anas, sptr, t_start, t_stop, NFFT):\n",
|
||||
"\n",
|
||||
" t_start = t_start * pq.s if t_start is not None else 0 * pq.s\n",
|
||||
" sampling_rate = anas.sampling_rate\n",
|
||||
" units = anas.units\n",
|
||||
" if t_start is not None and t_stop is not None:\n",
|
||||
" t_stop = t_stop * pq.s\n",
|
||||
" mask = (anas.times > t_start) & (anas.times < t_stop)\n",
|
||||
" anas = np.array(anas)[mask,:]\n",
|
||||
"\n",
|
||||
" cleaned_anas = zscore(clean_lfp(anas))\n",
|
||||
" \n",
|
||||
" cleaned_anas = neo.AnalogSignal(\n",
|
||||
" signal=cleaned_anas * units, sampling_rate=sampling_rate, t_start=t_start\n",
|
||||
" )\n",
|
||||
"\n",
|
||||
" sptr = neo.SpikeTrain(\n",
|
||||
" sptr.times[(sptr.times > t_start) & (sptr.times < cleaned_anas.times[-1])],\n",
|
||||
" t_start=t_start, t_stop=cleaned_anas.times[-1]\n",
|
||||
" )\n",
|
||||
"\n",
|
||||
" sigs, freqs = el.sta.spike_field_coherence(cleaned_anas, sptr, **{'nperseg': NFFT})\n",
|
||||
" return sigs, freqs, cleaned_anas"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"# action_id_0 = '1833-200619-1'\n",
|
||||
"# action_id_1 = '1833-200619-2'\n",
|
||||
"# lfp_0 = data_loader.lfp(action_id_0, 6)\n",
|
||||
"# lfp_1 = data_loader.lfp(action_id_1, 6)\n",
|
||||
"\n",
|
||||
"# lim_0 = get_lim(action_id_0)\n",
|
||||
"# lim_1 = get_lim(action_id_1)\n",
|
||||
"\n",
|
||||
"# sptrs_0 = data_loader.spike_trains(action_id_0, 6)\n",
|
||||
"# sptrs_1 = data_loader.spike_trains(action_id_1, 6)\n",
|
||||
"\n",
|
||||
"# coher_0, freq_0, clean_lfp_0 = compute_spike_lfp(lfp_0, sptrs_0[163], *lim_0, 4096)\n",
|
||||
"# coher_1, freq_1, clean_lfp_1 = compute_spike_lfp(lfp_1, sptrs_1[28], *lim_1, 4096)\n",
|
||||
"\n",
|
||||
"# best_channel_0 = np.argmax(signaltonoise(clean_lfp_0))\n",
|
||||
"# best_channel_1 = np.argmax(signaltonoise(clean_lfp_1))\n",
|
||||
"\n",
|
||||
"# plt.plot(freq_0, coher_0[:,best_channel_0])\n",
|
||||
"# plt.plot(freq_1, coher_1[:,best_channel_1])\n",
|
||||
"# plt.xlim(0,20)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"NFFT = 4096*2\n",
|
||||
"theta_band_f1, theta_band_f2 = 6, 10 "
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"def process(row):\n",
|
||||
" action_id = row['action']\n",
|
||||
" channel_group = row['channel_group']\n",
|
||||
" unit_name = row['unit_name']\n",
|
||||
" \n",
|
||||
" lfp = data_loader.lfp(action_id, channel_group)\n",
|
||||
" \n",
|
||||
" sptr = data_loader.spike_train(action_id, channel_group, unit_name)\n",
|
||||
" \n",
|
||||
" lim = get_lim(action_id)\n",
|
||||
"\n",
|
||||
" p_xys, freq, clean_lfp = compute_spike_lfp(lfp, sptr, *lim, NFFT=NFFT)\n",
|
||||
" \n",
|
||||
" snls = signaltonoise(clean_lfp)\n",
|
||||
" best_channel = np.argmax(snls)\n",
|
||||
" snl = snls[best_channel]\n",
|
||||
" p_xy = p_xys[:,best_channel]\n",
|
||||
" \n",
|
||||
" theta_f, theta_p_max = find_theta_peak(p_xy, freq, theta_band_f1, theta_band_f2)\n",
|
||||
" \n",
|
||||
" theta_energy = compute_energy(p_xy, freq, theta_band_f1, theta_band_f2) # theta band 6 - 10 Hz\n",
|
||||
" \n",
|
||||
" theta_half_f1, theta_half_f2 = compute_half_width(p_xy, freq, theta_p_max, theta_f)\n",
|
||||
" \n",
|
||||
" theta_half_width = theta_half_f2 - theta_half_f1\n",
|
||||
" \n",
|
||||
" theta_half_energy = compute_energy(p_xy, freq, theta_half_f1, theta_half_f2) # theta band 6 - 10 Hz\n",
|
||||
" \n",
|
||||
" # stim\n",
|
||||
" \n",
|
||||
" stim_freq = compute_stim_freq(action_id)\n",
|
||||
" \n",
|
||||
" stim_p_max = compute_stim_peak(p_xy, freq, stim_freq)\n",
|
||||
" \n",
|
||||
" stim_half_f1, stim_half_f2 = compute_half_width(p_xy, freq, stim_p_max, stim_freq)\n",
|
||||
" stim_half_width = stim_half_f2 - stim_half_f1\n",
|
||||
" \n",
|
||||
" stim_energy = compute_energy(p_xy, freq, stim_half_f1, stim_half_f2)\n",
|
||||
" \n",
|
||||
" result = pd.Series({\n",
|
||||
" 'signal_to_noise': snl,\n",
|
||||
" 'best_channel': best_channel,\n",
|
||||
" 'theta_freq': theta_f,\n",
|
||||
" 'theta_peak': theta_p_max,\n",
|
||||
" 'theta_energy': theta_energy,\n",
|
||||
" 'theta_half_f1': theta_half_f1, \n",
|
||||
" 'theta_half_f2': theta_half_f2,\n",
|
||||
" 'theta_half_width': theta_half_width,\n",
|
||||
" 'theta_half_energy': theta_half_energy,\n",
|
||||
" 'stim_freq': stim_freq,\n",
|
||||
" 'stim_p_max': stim_p_max,\n",
|
||||
" 'stim_half_f1': stim_half_f1, \n",
|
||||
" 'stim_half_f2': stim_half_f2,\n",
|
||||
" 'stim_half_width': stim_half_width,\n",
|
||||
" 'stim_energy': stim_energy\n",
|
||||
" })\n",
|
||||
" return result"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [
|
||||
{
|
||||
"data": {
|
||||
"application/vnd.jupyter.widget-view+json": {
|
||||
"model_id": "e28818ca79be42abaf9c2b6106580c75",
|
||||
"version_major": 2,
|
||||
"version_minor": 0
|
||||
},
|
||||
"text/plain": [
|
||||
"HBox(children=(IntProgress(value=0, max=1298), HTML(value='')))"
|
||||
]
|
||||
},
|
||||
"metadata": {},
|
||||
"output_type": "display_data"
|
||||
},
|
||||
{
|
||||
"name": "stderr",
|
||||
"output_type": "stream",
|
||||
"text": [
|
||||
"/home/mikkel/.virtualenvs/expipe/lib/python3.6/site-packages/scipy/signal/spectral.py:1578: RuntimeWarning: invalid value encountered in true_divide\n",
|
||||
" Cxy = np.abs(Pxy)**2 / Pxx / Pyy\n"
|
||||
]
|
||||
}
|
||||
],
|
||||
"source": [
|
||||
"results = units.merge(\n",
|
||||
" units.progress_apply(process, axis=1), \n",
|
||||
" left_index=True, right_index=True)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"%debug"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "markdown",
|
||||
"metadata": {},
|
||||
"source": [
|
||||
"# Save to expipe"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"action = project.require_action(\"stimulus-spike-lfp-response\")"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"action.modules['parameters'] = {\n",
|
||||
" 'NFFT': NFFT,\n",
|
||||
" 'theta_band_f1': theta_band_f1,\n",
|
||||
" 'theta_band_f2': theta_band_f2\n",
|
||||
"}"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"action.data['results'] = 'results.csv'\n",
|
||||
"results.to_csv(action.data_path('results'), index=False)"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"copy_tree(output, str(action.data_path()))"
|
||||
]
|
||||
},
|
||||
{
|
||||
"cell_type": "code",
|
||||
"execution_count": null,
|
||||
"metadata": {},
|
||||
"outputs": [],
|
||||
"source": [
|
||||
"septum_mec.analysis.registration.store_notebook(action, \"10-calculate-stimulus-spike-lfp-response.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
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
|
|
@ -0,0 +1,3 @@
|
|||
NFFT: 8192
|
||||
theta_band_f1: 6
|
||||
theta_band_f2: 10
|
||||
Loading…
Reference in New Issue