Skip to content

Commit

Permalink
Add SigmaPeak timestamp in SW detection when coupling=True
Browse files Browse the repository at this point in the history
  • Loading branch information
@raphaelvallat
raphaelvallat committed Jun 30, 2020
1 parent aa22ca5 commit 34d3d46
Show file tree
Hide file tree
Showing 5 changed files with 158 additions and 163 deletions.
124 changes: 35 additions & 89 deletions notebooks/05_sw_detection.ipynb
View file

Large diffs are not rendered by default.

76 changes: 45 additions & 31 deletions notebooks/12_spindles-SO_coupling.ipynb
View file
Original file line number Diff line number Diff line change
Expand Up @@ -63,7 +63,7 @@
"name": "stderr",
"output_type": "stream",
"text": [
"29-Jun-20 17:35:22 | WARNING | Hypnogram is SHORTER than data by 10.58 seconds. Padding hypnogram with last value to match data.size.\n"
"30-Jun-20 14:53:02 | WARNING | Hypnogram is SHORTER than data by 10.58 seconds. Padding hypnogram with last value to match data.size.\n"
]
},
{
Expand Down Expand Up @@ -135,6 +135,7 @@
" <th>PTP</th>\n",
" <th>Slope</th>\n",
" <th>Frequency</th>\n",
" <th>SigmaPeak</th>\n",
" <th>PhaseAtSigmaPeak</th>\n",
" <th>ndPAC</th>\n",
" <th>Stage</th>\n",
Expand All @@ -156,6 +157,7 @@
" <td>139.147</td>\n",
" <td>632.488</td>\n",
" <td>0.541</td>\n",
" <td>573.46</td>\n",
" <td>0.255</td>\n",
" <td>0.202</td>\n",
" <td>2</td>\n",
Expand All @@ -175,6 +177,7 @@
" <td>168.178</td>\n",
" <td>800.849</td>\n",
" <td>1.250</td>\n",
" <td>592.75</td>\n",
" <td>-1.174</td>\n",
" <td>0.467</td>\n",
" <td>2</td>\n",
Expand All @@ -194,6 +197,7 @@
" <td>89.251</td>\n",
" <td>153.880</td>\n",
" <td>0.781</td>\n",
" <td>631.77</td>\n",
" <td>1.340</td>\n",
" <td>0.422</td>\n",
" <td>2</td>\n",
Expand All @@ -213,6 +217,7 @@
" <td>126.089</td>\n",
" <td>700.497</td>\n",
" <td>1.111</td>\n",
" <td>644.13</td>\n",
" <td>1.175</td>\n",
" <td>0.083</td>\n",
" <td>2</td>\n",
Expand All @@ -232,6 +237,7 @@
" <td>91.355</td>\n",
" <td>397.195</td>\n",
" <td>1.075</td>\n",
" <td>644.13</td>\n",
" <td>1.175</td>\n",
" <td>0.154</td>\n",
" <td>2</td>\n",
Expand All @@ -256,6 +262,7 @@
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" <td>...</td>\n",
" </tr>\n",
" <tr>\n",
" <th>2850</th>\n",
Expand All @@ -270,6 +277,7 @@
" <td>92.155</td>\n",
" <td>460.774</td>\n",
" <td>1.587</td>\n",
" <td>20512.79</td>\n",
" <td>0.995</td>\n",
" <td>0.038</td>\n",
" <td>2</td>\n",
Expand All @@ -289,6 +297,7 @@
" <td>99.374</td>\n",
" <td>129.057</td>\n",
" <td>0.483</td>\n",
" <td>20520.65</td>\n",
" <td>-1.707</td>\n",
" <td>0.394</td>\n",
" <td>2</td>\n",
Expand All @@ -308,6 +317,7 @@
" <td>244.881</td>\n",
" <td>1020.336</td>\n",
" <td>0.725</td>\n",
" <td>20525.08</td>\n",
" <td>0.815</td>\n",
" <td>0.384</td>\n",
" <td>2</td>\n",
Expand All @@ -327,6 +337,7 @@
" <td>127.652</td>\n",
" <td>531.884</td>\n",
" <td>0.694</td>\n",
" <td>20532.43</td>\n",
" <td>-0.872</td>\n",
" <td>0.158</td>\n",
" <td>2</td>\n",
Expand All @@ -346,6 +357,7 @@
" <td>137.480</td>\n",
" <td>808.708</td>\n",
" <td>0.488</td>\n",
" <td>20540.00</td>\n",
" <td>2.743</td>\n",
" <td>0.200</td>\n",
" <td>2</td>\n",
Expand All @@ -354,7 +366,7 @@
" </tr>\n",
" </tbody>\n",
"</table>\n",
"<p>2855 rows × 16 columns</p>\n",
"<p>2855 rows × 17 columns</p>\n",
"</div>"
],
"text/plain": [
Expand All @@ -371,33 +383,33 @@
"2853 20531.47 20531.70 20531.94 20532.58 20532.91 1.44 \n",
"2854 20537.49 20538.62 20538.79 20539.02 20539.54 2.05 \n",
"\n",
" ValNegPeak ValPosPeak PTP Slope Frequency PhaseAtSigmaPeak \\\n",
"0 -73.967 65.180 139.147 632.488 0.541 0.255 \n",
"1 -112.435 55.743 168.178 800.849 1.250 -1.174 \n",
"2 -45.324 43.926 89.251 153.880 0.781 1.340 \n",
"3 -60.872 65.218 126.089 700.497 1.111 1.175 \n",
"4 -55.222 36.133 91.355 397.195 1.075 1.175 \n",
"... ... ... ... ... ... ... \n",
"2850 -66.765 25.390 92.155 460.774 1.587 0.995 \n",
"2851 -44.161 55.213 99.374 129.057 0.483 -1.707 \n",
"2852 -163.348 81.533 244.881 1020.336 0.725 0.815 \n",
"2853 -50.424 77.228 127.652 531.884 0.694 -0.872 \n",
"2854 -54.047 83.434 137.480 808.708 0.488 2.743 \n",
" ValNegPeak ValPosPeak PTP Slope Frequency SigmaPeak \\\n",
"0 -73.967 65.180 139.147 632.488 0.541 573.46 \n",
"1 -112.435 55.743 168.178 800.849 1.250 592.75 \n",
"2 -45.324 43.926 89.251 153.880 0.781 631.77 \n",
"3 -60.872 65.218 126.089 700.497 1.111 644.13 \n",
"4 -55.222 36.133 91.355 397.195 1.075 644.13 \n",
"... ... ... ... ... ... ... \n",
"2850 -66.765 25.390 92.155 460.774 1.587 20512.79 \n",
"2851 -44.161 55.213 99.374 129.057 0.483 20520.65 \n",
"2852 -163.348 81.533 244.881 1020.336 0.725 20525.08 \n",
"2853 -50.424 77.228 127.652 531.884 0.694 20532.43 \n",
"2854 -54.047 83.434 137.480 808.708 0.488 20540.00 \n",
"\n",
" ndPAC Stage Channel IdxChannel \n",
"0 0.202 2 Cz 0 \n",
"1 0.467 2 Cz 0 \n",
"2 0.422 2 Cz 0 \n",
"3 0.083 2 Cz 0 \n",
"4 0.154 2 Cz 0 \n",
"... ... ... ... ... \n",
"2850 0.038 2 Cz 0 \n",
"2851 0.394 2 Cz 0 \n",
"2852 0.384 2 Cz 0 \n",
"2853 0.158 2 Cz 0 \n",
"2854 0.200 2 Cz 0 \n",
" PhaseAtSigmaPeak ndPAC Stage Channel IdxChannel \n",
"0 0.255 0.202 2 Cz 0 \n",
"1 -1.174 0.467 2 Cz 0 \n",
"2 1.340 0.422 2 Cz 0 \n",
"3 1.175 0.083 2 Cz 0 \n",
"4 1.175 0.154 2 Cz 0 \n",
"... ... ... ... ... ... \n",
"2850 0.995 0.038 2 Cz 0 \n",
"2851 -1.707 0.394 2 Cz 0 \n",
"2852 0.815 0.384 2 Cz 0 \n",
"2853 -0.872 0.158 2 Cz 0 \n",
"2854 2.743 0.200 2 Cz 0 \n",
"\n",
"[2855 rows x 16 columns]"
"[2855 rows x 17 columns]"
]
},
"execution_count": 3,
Expand All @@ -415,9 +427,11 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"As explained in the documentation of the [yasa.sw_detect](https://raphaelvallat.com/yasa/build/html/generated/yasa.sw_detect.html), we get two coupling metrics:\n",
"As explained in the documentation of the [yasa.sw_detect](https://raphaelvallat.com/yasa/build/html/generated/yasa.sw_detect.html), we get some additional columns in the output dataframe:\n",
"\n",
"1) The ``SigmaPeak`` column contains the timestamp (in seconds from the beginning of the recording) where the sigma-filtered (12 - 16 Hz, see ``freq_sp``) amplitude is at its maximal. This is calculated separately for each detected slow-wave, using a **4-seconds epoch centered around the negative peak of the slow-wave**.\n",
"\n",
"1) The ``PhaseAtSigmaPeak`` column contains the phase (in radians) of the slow-wave filtered signal (0.3 - 2 Hz, see ``freq_sw``) where the sigma-filtered (12 - 16 Hz, see ``freq_sp``) amplitude is at its maximal. This is calculated separately for each detected slow-wave, using a **4-seconds epoch centered around the negative peak of the slow-wave**. Using the [Pingouin package](https://pingouin-stats.org/), we can then easily extract and visualize the direction and strength of coupling across all channels:"
"2) The ``PhaseAtSigmaPeak`` column contains the phase (in radians) of the slow-wave filtered signal (0.3 - 2 Hz, see ``freq_sw``) at ``SigmaPeak``. Using the [Pingouin package](https://pingouin-stats.org/), we can then easily extract and visualize the direction and strength of coupling across all channels:"
]
},
{
Expand Down Expand Up @@ -455,7 +469,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
"2) The ``ndPAC`` columns contains the normalized mean vector length (also called the normalized direct PAC, see [Ozkurt 2012](https://doi.org/10.1109/TBME.2012.2194783)). Note that ``ndPAC`` should be highly correlated the vector length of the ``PhaseAtSigmaPeak``. It may be more accurate though since it is calculated on the entire 4-seconds epoch.\n",
"3) The ``ndPAC`` columns contains the normalized mean vector length (also called the normalized direct PAC, see [Ozkurt 2012](https://doi.org/10.1109/TBME.2012.2194783)). Note that ``ndPAC`` should be highly correlated the vector length of the ``PhaseAtSigmaPeak``. It may be more accurate though since it is calculated on the entire 4-seconds epoch.\n",
"\n",
"For more details, please refer to:\n",
"\n",
Expand Down Expand Up @@ -491,7 +505,7 @@
{
"data": {
"text/plain": [
"0.21858686994116763"
"0.2185868699411677"
]
},
"execution_count": 6,
Expand Down
118 changes: 76 additions & 42 deletions yasa/main.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -212,7 +212,7 @@ def get_sync_events(self, center, time_before, time_after,
assert time_after >= 0
bef = int(self._sf * time_before)
aft = int(self._sf * time_after)
# Step size is determined by sf: 0.01 sec at 100 Hz, 0.002 sec at
# TODO: Step size is determined by sf: 0.01 sec at 100 Hz, 0.002 sec at
# 500 Hz, 0.00390625 sec at 256 Hz. Should we add a step_size=0.01
# option?
time = np.arange(-bef, aft + 1, dtype='int') / self._sf
Expand Down Expand Up @@ -949,12 +949,15 @@ def sw_detect(data, sf=None, ch_names=None, hypno=None, include=(2, 3),
coupling : boolean
If True, YASA will also calculate the phase-amplitude coupling between
the slow-waves phase and the spindles-related sigma band
amplitude. Specifically, two PAC metrics will be calculated:
amplitude. Specifically, the following columns will be added to the
output dataframe:
1. ``PhaseAtSigmaPeak``: the phase of the bandpas-filtered slow-wave
signal (in radians) at the maximum sigma peak amplitude within an
4-seconds epoch centered around the negative peak (through) of the
current slow-wave.
1. ``'SigmaPeak'``: The location (in seconds) of the maximum sigma peak
amplitude within a 4-seconds epoch centered around the negative peak
(through) of the current slow-wave.
2. ``PhaseAtSigmaPeak``: the phase of the bandpas-filtered slow-wave
signal (in radians) at ``'SigmaPeak'``.
Importantly, since ``PhaseAtSigmaPeak`` is expressed in radians,
one should use circular statistics to calculate the mean direction
Expand All @@ -966,7 +969,7 @@ def sw_detect(data, sf=None, ch_names=None, hypno=None, include=(2, 3),
mean_direction = pg.circ_mean(sw['PhaseAtSigmaPeak'])
vector_length = pg.circ_r(sw['PhaseAtSigmaPeak'])
2. ``ndPAC``: the normalized Mean Vector Length
3. ``ndPAC``: the normalized Mean Vector Length
(also called the normalized direct PAC, or ndPAC) within a 4-sec
epoch centered around the negative peak of the slow-wave.
Expand Down Expand Up @@ -1039,9 +1042,12 @@ def sw_detect(data, sf=None, ch_names=None, hypno=None, include=(2, 3),
* ``'Slope'``: Slope between ``NegPeak`` and ``MidCrossing`` (in uV/sec,
calculated on the ``freq_sw`` bandpass-filtered signal)
* ``'Frequency'``: Frequency of the slow-wave (= 1 / ``Duration``)
* ``'SigmaPeak'``: Location of the sigma peak amplitude within a 4-sec
epoch centered around the negative peak of the slow-wave. This is only
calculated when ``coupling=True``.
* ``'PhaseAtSigmaPeak'``: SW phase at max sigma amplitude within a
4-sec epoch centered the negative peak of the slow-wave. This is only
calculated when ``coupling=True``
4-sec epoch centered around the negative peak of the slow-wave. This is
only calculated when ``coupling=True``
* ``'ndPAC'``: Normalized direct PAC within a 4-sec epoch centered
the negative peak of the slow-wave. This is only calculated when
``coupling=True``
Expand Down Expand Up @@ -1224,52 +1230,80 @@ def sw_detect(data, sf=None, ch_names=None, hypno=None, include=(2, 3),
logger.warning('No SW were found in channel %s.', ch_names[i])
continue

# Create a dictionnary and then a dataframe
sw_params = OrderedDict({'Start': sw_start,
'NegPeak': sw_idx_neg,
'MidCrossing': sw_midcrossing,
'PosPeak': sw_idx_pos,
'End': sw_end,
'Duration': sw_dur,
'ValNegPeak': data_filt[i, idx_neg_peaks],
'ValPosPeak': data_filt[i, idx_pos_peaks],
'PTP': sw_ptp,
'Slope': sw_slope,
'Frequency': 1 / sw_dur,
'Stage': sw_sta,
})
# Filter good events
idx_neg_peaks = idx_neg_peaks[good_sw]
idx_pos_peaks = idx_pos_peaks[good_sw]
sw_start = sw_start[good_sw]
sw_idx_neg = sw_idx_neg[good_sw]
sw_midcrossing = sw_midcrossing[good_sw]
sw_idx_pos = sw_idx_pos[good_sw]
sw_end = sw_end[good_sw]
sw_dur = sw_dur[good_sw]
sw_ptp = sw_ptp[good_sw]
sw_slope = sw_slope[good_sw]
sw_sta = sw_sta[good_sw]

# Create a dictionnary
sw_params = OrderedDict({
'Start': sw_start,
'NegPeak': sw_idx_neg,
'MidCrossing': sw_midcrossing,
'PosPeak': sw_idx_pos,
'End': sw_end,
'Duration': sw_dur,
'ValNegPeak': data_filt[i, idx_neg_peaks],
'ValPosPeak': data_filt[i, idx_pos_peaks],
'PTP': sw_ptp,
'Slope': sw_slope,
'Frequency': 1 / sw_dur,
'Stage': sw_sta,
})

# Add phase (in radians) of slow-oscillation signal at maximum
# spindles-related sigma amplitude within a 4-seconds centered epochs.
if coupling:
# Get Phase and Amplitude for each centered epoch
# Get phase and amplitude for each centered epoch
# TODO: allow user-specified window size.
time_before = time_after = 2
bef = int(sf * time_before)
aft = int(sf * time_after)
# Center of each epoch is defined as the negative peak of the SW
n_peaks = idx_neg_peaks.shape[0]
idx, idx_nomask = get_centered_indices(data[i, :], idx_neg_peaks,
bef, aft)
# idx.shape = (len(idx_valid), bef + aft + 1)
idx, idx_valid = get_centered_indices(data[i, :], idx_neg_peaks,
bef, aft)
sw_pha_ev = sw_pha[i, idx]
sp_amp_ev = sp_amp[i, idx]
# Find SW phase at max sigma amplitude in epoch
idx_max_amp = sp_amp_ev.argmax(axis=1)[..., None]
pha_at_max = np.take_along_axis(sw_pha_ev, idx_max_amp, axis=1)
pha_at_max = np.squeeze(pha_at_max)
# 1) Find location of max sigma amplitude in epoch
idx_max_amp = sp_amp_ev.argmax(axis=1)
# Now we need to append it back to the original unmasked shape
pha_at_max_full = np.ones(n_peaks) * np.nan
pha_at_max_full[idx_nomask] = pha_at_max
sw_params['PhaseAtSigmaPeak'] = pha_at_max_full
# Normalized Direct PAC, without thresholding
ndp = tpm.ndpac(sw_pha_ev[None, ...], sp_amp_ev[None, ...], p=1)
ndp = np.squeeze(ndp)
ndp_full = np.ones(n_peaks) * np.nan
ndp_full[idx_nomask] = ndp
sw_params['ndPAC'] = ndp_full
# to avoid error when idx.shape[0] != idx_valid.shape, i.e.
# some epochs were out of data bounds.
sw_params['SigmaPeak'] = np.ones(n_peaks) * np.nan
# Timestamp at sigma peak, expressed in seconds from negative peak
# e.g. -0.39, 0.5, 1, 2 -- limits are [time_before, time_after]
time_sigpk = (idx_max_amp - bef) / sf
# convert to absolute time from beginning of the recording
# time_sigpk only includes valid epoch
time_sigpk_abs = sw_idx_neg[idx_valid] + time_sigpk
sw_params['SigmaPeak'][idx_valid] = time_sigpk_abs
# 2) PhaseAtSigmaPeak
# Find SW phase at max sigma amplitude in epoch
pha_at_max = np.squeeze(np.take_along_axis(sw_pha_ev,
idx_max_amp[..., None],
axis=1))
sw_params['PhaseAtSigmaPeak'] = np.ones(n_peaks) * np.nan
sw_params['PhaseAtSigmaPeak'][idx_valid] = pha_at_max
# 3) Normalized Direct PAC, without thresholding
ndp = np.squeeze(tpm.ndpac(sw_pha_ev[None, ...],
sp_amp_ev[None, ...], p=1))
sw_params['ndPAC'] = np.ones(n_peaks) * np.nan
sw_params['ndPAC'][idx_valid] = ndp
# Make sure that Stage is the last column of the dataframe
sw_params.move_to_end('Stage')

# Convert to dataframe, keeping only good events
df_chan = pd.DataFrame(sw_params)[good_sw]
df_chan = pd.DataFrame(sw_params)

# Remove all duplicates
df_chan = df_chan.drop_duplicates(subset=['Start'], keep=False)
Expand Down Expand Up @@ -1796,9 +1830,9 @@ def get_sync_events(self, center='Peak', time_before=0.4, time_after=0.4,
# Get location of peaks in data
peaks = (self._events[center] * self._sf).astype(int).to_numpy()
# Get centered indices (here we could use second channel as well).
idx, idx_nomask = get_centered_indices(data[0, :], peaks, bef, aft)
idx, idx_valid = get_centered_indices(data[0, :], peaks, bef, aft)
# If no good epochs are returned raise a warning
assert len(idx_nomask), (
assert len(idx_valid), (
'Time before and/or time after exceed data bounds, please '
'lower the temporal window around center.')

Expand Down
Loading

0 comments on commit 34d3d46

Please sign in to comment.