From b4aece8a3ece5d009a6f2bf431d4b39097b097c3 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 12:15:02 +0100 Subject: [PATCH 1/9] Fixing truncated test pulse If the test pulse is truncated, skip only the first change in the difference. --- ipfx/stim_features.py | 9 ++++++++- tests/test_stim_features.py | 5 +++++ 2 files changed, 13 insertions(+), 1 deletion(-) diff --git a/ipfx/stim_features.py b/ipfx/stim_features.py index 0158442a..0f15718b 100644 --- a/ipfx/stim_features.py +++ b/ipfx/stim_features.py @@ -9,7 +9,14 @@ def get_stim_characteristics(i, t, test_pulse=True): di = np.diff(i) di_idx = np.flatnonzero(di) # != 0 - start_idx_idx = 2 if test_pulse else 0 # skip the first up/down (test pulse) if present + + if test_pulse: + if len(di_idx) % 2 == 1: # skip the truncated test pulse + start_idx_idx = 1 + else: + start_idx_idx = 2 + else: + start_idx_idx = 0 if len(di_idx[start_idx_idx:]) == 0: # if no stimulus is found return None, None, 0.0, None, None diff --git a/tests/test_stim_features.py b/tests/test_stim_features.py index ed7557c2..510f38f5 100644 --- a/tests/test_stim_features.py +++ b/tests/test_stim_features.py @@ -43,6 +43,11 @@ [0, 1, 1, 0], False, (1, 2, 1, 1, 2) + ), + ( + [1, 1, 0, 0, 0, 2, 2, 2, 2, 0, 0], + True, + (5, 4, 2, 5, 8) ) ] From 13d111fab2c6ba30c7dec699cb835012410e3234 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 12:30:12 +0100 Subject: [PATCH 2/9] Fixing nan parameter for xlim Skip nan values to avoid error caused when it is passed to xlim as a parameter. --- ipfx/plot_qc_figures.py | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/ipfx/plot_qc_figures.py b/ipfx/plot_qc_figures.py index 949eecfe..c564b0b9 100644 --- a/ipfx/plot_qc_figures.py +++ b/ipfx/plot_qc_figures.py @@ -65,8 +65,8 @@ def plot_single_ap_values(data_set, sweep_numbers, lims_features, sweep_features else: rheo_sn = cell_features["long_squares"]["rheobase_sweep"]["sweep_number"] rheo_spike = get_spikes(sweep_features, rheo_sn)[0] - voltages = [ rheo_spike[f] for f in voltage_features] - times = [ rheo_spike[f] for f in time_features] + voltages = [rheo_spike[f] for f in voltage_features if not np.isnan(rheo_spike[f])] + times = [rheo_spike[f] for f in time_features if not np.isnan(rheo_spike[f])] plt.figure(figs[0].number) plt.scatter(range(len(voltages)), voltages, color='gray') @@ -123,13 +123,13 @@ def plot_single_ap_values(data_set, sweep_numbers, lims_features, sweep_features if nspikes: if type_name != "long_square" and nspikes: - voltages = [spikes[0][f] for f in voltage_features] - times = [spikes[0][f] for f in time_features] + voltages = [spikes[0][f] for f in voltage_features if not np.isnan(spikes[0][f])] + times = [spikes[0][f] for f in time_features if not np.isnan(spikes[0][f])] else: rheo_sn = cell_features["long_squares"]["rheobase_sweep"]["sweep_number"] rheo_spike = get_spikes(sweep_features, rheo_sn)[0] - voltages = [rheo_spike[f] for f in voltage_features] - times = [rheo_spike[f] for f in time_features] + voltages = [rheo_spike[f] for f in voltage_features if not np.isnan(rheo_spike[f])] + times = [rheo_spike[f] for f in time_features if not np.isnan(rheo_spike[f])] plt.scatter(times, voltages, color='red', zorder=20) @@ -144,7 +144,7 @@ def plot_single_ap_values(data_set, sweep_numbers, lims_features, sweep_features if type_name == "ramp": if nspikes: - plt.xlim(spikes[0]["threshold_t"] - 0.002, spikes[0]["fast_trough_t"] + 0.01) + plt.xlim(times[0] - 0.002, times[-2] + 0.01) elif type_name == "short_square": plt.xlim(stim_start_shifted - 0.002, stim_start_shifted + stim_dur + 0.01) elif type_name == "long_square": From 44658a745bb270d64e18480b845e9fe600d7da12 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 14:41:56 +0100 Subject: [PATCH 3/9] Fix negative stimulus issue Handling above baseline stimulus with negative value --- ipfx/subthresh_features.py | 2 +- tests/test_subthresh_features.py | 34 ++++++++++++++++++++++++++++++++ 2 files changed, 35 insertions(+), 1 deletion(-) diff --git a/ipfx/subthresh_features.py b/ipfx/subthresh_features.py index 7010e47b..62d554a8 100644 --- a/ipfx/subthresh_features.py +++ b/ipfx/subthresh_features.py @@ -54,7 +54,7 @@ def voltage_deflection(t, v, i, start, end, deflect_type=None): if deflect_type is None: if i is not None: halfway_index = tsu.find_time_index(t, (end - start) / 2. + start) - if i[halfway_index] >= 0: + if (i[halfway_index] - i[start_index - 1]) >= 0: deflect_type = "max" else: deflect_type = "min" diff --git a/tests/test_subthresh_features.py b/tests/test_subthresh_features.py index 791d2248..bbd47560 100644 --- a/tests/test_subthresh_features.py +++ b/tests/test_subthresh_features.py @@ -1,5 +1,6 @@ import ipfx.subthresh_features as subf import numpy as np +import pytest def test_input_resistance(): @@ -76,3 +77,36 @@ def test_time_constant_noise_acceptance(): tau = subf.time_constant(t, v, i, start=start, end=end) assert np.isclose(actual_tau, tau, rtol=1e-3) + + +test_params = [ + ( + [-5, -5, -5, -2, 0, 1, 1, 1, -2 -5, -5], + [0, 0, 0, 1, 1, 1, 1, 1, 1, 0, 0], + 3, + 8, + (1, 5) + ), + ( + [-5, -5, -5, -7, -10, -12, -12, -12, -7 - 5, -5], + [0, 0, 0, -1, -1, -1, -1, -1, -1, 0, 0], + 3, + 8, + (-12, 5) + ), + ( + [-70, -70, -70, -50, -20, -30, -30, -30, -50, -70, -70], + [-70, -70, -70, -20, -20, -20, -20, -20, -20, -70, -70], + 3, + 8, + (-20, 4) + ), +] + + +@pytest.mark.parametrize('v, i, start, end, deflection_result', test_params) +def test_voltage_deflection(v, i, start, end, deflection_result): + t = np.arange(0, 10) + deflection_v, deflection_idx = subf.voltage_deflection(t, v, i, start, end) + assert deflection_v == deflection_result[0] + assert deflection_idx == deflection_result[1] From 983159d12cdfa9de022d75a927c25378cf9dbadc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 14:54:07 +0100 Subject: [PATCH 4/9] Modification of input resistance calculation calculating the V and I values as an average instead of the min value --- ipfx/subthresh_features.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/ipfx/subthresh_features.py b/ipfx/subthresh_features.py index 62d554a8..91b7bf14 100644 --- a/ipfx/subthresh_features.py +++ b/ipfx/subthresh_features.py @@ -156,9 +156,10 @@ def input_resistance(t_set, i_set, v_set, start, end, baseline_interval=0.1): v_vals = [] i_vals = [] for t, i, v, in zip(t_set, i_set, v_set): - v_peak, min_index = voltage_deflection(t, v, i, start, end, 'min') - v_vals.append(v_peak) - i_vals.append(i[min_index]) + v_stim_avg = tsu.average_voltage(v, t, start, end) + i_stim_avg = tsu.average_voltage(i, t, start, end) + v_vals.append(v_stim_avg) + i_vals.append(i_stim_avg) v = np.array(v_vals) i = np.array(i_vals) From e6e762b38f613ec25e10fa1e5834c52969f3a799 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 15:55:18 +0100 Subject: [PATCH 5/9] Fixing comparison of np.nan warning If any of the trough_indexes is nan, it will be ignored at the width level correction. --- ipfx/spike_features.py | 4 +++- tests/test_spike_features.py | 15 ++++++++++++++- 2 files changed, 17 insertions(+), 2 deletions(-) diff --git a/ipfx/spike_features.py b/ipfx/spike_features.py index 728d4be8..a70cfa7d 100644 --- a/ipfx/spike_features.py +++ b/ipfx/spike_features.py @@ -84,7 +84,9 @@ def find_widths(v, t, spike_indexes, peak_indexes, trough_indexes, clipped=None) # Some spikes in burst may have deep trough but short height, so can't use same # definition for width - width_levels[width_levels < v[spike_indexes]] = thresh_to_peak_levels[width_levels < v[spike_indexes]] + used_indexes = np.argwhere(use_indexes) + change_indexes = used_indexes[width_levels[use_indexes] < v[spike_indexes[use_indexes]]] + width_levels[change_indexes] = thresh_to_peak_levels[change_indexes] width_starts = np.zeros_like(trough_indexes) * np.nan width_starts[use_indexes] = np.array([pk - np.flatnonzero(v[pk:spk:-1] <= wl)[0] if diff --git a/tests/test_spike_features.py b/tests/test_spike_features.py index 6d5d3121..737b6b96 100644 --- a/tests/test_spike_features.py +++ b/tests/test_spike_features.py @@ -49,7 +49,7 @@ def test_width_calculation(spike_test_pair): expected_widths = np.array([0.000545, 0.000585]) assert np.allclose( - spkf.find_widths(v, t, spikes, peaks, troughs), + spkf.find_widths(v, t, spikes, peaks, troughs), expected_widths ) @@ -64,3 +64,16 @@ def test_width_calculation_too_many_troughs(spike_test_pair): with pytest.raises(er.FeatureError): spkf.find_widths(v, t, spikes, peaks, troughs) + + +def test_width_calculation_missing_troughs(spike_test_pair): + data = spike_test_pair + t = data[:, 0] + v = data[:, 1] + spikes = np.array([725, 3382]) + peaks = np.array([812, 3478]) + troughs = np.array([1089, np.nan]) + + found_width = spkf.find_widths(v, t, spikes, peaks, troughs) + assert np.isclose(found_width[0], 0.000545) + assert np.isnan(found_width[1]) From 26c897c7c5679d93f5e676e16ce2691d9513aa2b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 17:11:00 +0100 Subject: [PATCH 6/9] Fixing spike too close to the sweep end If the spike start is closer to the end than 2*max_interval, we drop it. --- ipfx/spike_detector.py | 4 ++++ tests/test_spike_detector.py | 13 +++++++++++++ 2 files changed, 17 insertions(+) diff --git a/ipfx/spike_detector.py b/ipfx/spike_detector.py index 34f0b2fa..5cfa8ed5 100644 --- a/ipfx/spike_detector.py +++ b/ipfx/spike_detector.py @@ -274,6 +274,10 @@ def check_thresholds_and_peaks(v, t, spike_indexes, peak_indexes, upstroke_index # Find the peak in a window twice the size of our allowed window spike = spike_indexes[i] + if t[spike] + 2 * max_interval > t[-1]: + drop_spikes.append(i) + continue + t_0 = tsu.find_time_index(t, t[spike] + 2 * max_interval) new_peak = np.argmax(v[spike:t_0]) + spike diff --git a/tests/test_spike_detector.py b/tests/test_spike_detector.py index 5ce0fa90..a95f1fd6 100644 --- a/tests/test_spike_detector.py +++ b/tests/test_spike_detector.py @@ -174,6 +174,19 @@ def test_check_spikes_and_peaks(): assert np.allclose(new_peaks, peaks[1:]) +def test_check_spikes_and_peaks_close_to_end(): + t = np.arange(0, 2000) * 5e-6 + v = np.zeros_like(t) + spikes = np.array([500]) + peaks = np.array([1510]) + upstrokes = np.array([700]) + dvdt = np.ones_like(t) + + new_spikes, new_peaks, new_upstrokes, clipped = spkd.check_thresholds_and_peaks(v, t, spikes, peaks, upstrokes, dvdt=dvdt) + assert np.allclose(new_spikes, spikes[:-1]) + assert np.allclose(new_peaks, peaks[1:]) + + def test_thresholds(spike_test_pair): data = spike_test_pair t = data[:, 0] From 03338addf9c51a31014ea65766357de54eee7820 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 17:24:25 +0100 Subject: [PATCH 7/9] Fix truncated test pulse at test epoch detection If the sweep has only a single test pulse, returns None instead of an exception. --- ipfx/epochs.py | 3 ++- tests/test_epochs.py | 6 ++++++ 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/ipfx/epochs.py b/ipfx/epochs.py index f15e3b93..51ee5091 100644 --- a/ipfx/epochs.py +++ b/ipfx/epochs.py @@ -149,7 +149,8 @@ def get_test_epoch(i,hz): return None if len(di_idx) == 1: - raise Exception("Cannot detect and end to the test pulse") + print("Cannot detect and end to the test pulse") + return None start_pulse_idx = di_idx[0] + 1 # shift by one to compensate for diff() end_pulse_idx = di_idx[1] diff --git a/tests/test_epochs.py b/tests/test_epochs.py index 6364f0dc..dc0913cb 100644 --- a/tests/test_epochs.py +++ b/tests/test_epochs.py @@ -72,6 +72,12 @@ def test_get_experiment_epoch(i, sampling_rate, expt_epoch): None ), + # truncated test pulse + ( + [1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0], + 4, + None + ), ] ) def test_get_test_epoch(i, sampling_rate, test_epoch): From 4b85ec60ee5248bad90a0f7d35a067a013dd430b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 17:50:26 +0100 Subject: [PATCH 8/9] Fix dvdt calculation of a sweep with low sample rate Calculate dvdt without filtering if the sampling rate is low. --- ipfx/time_series_utils.py | 25 ++++++++++++++----------- tests/test_time_series_utils.py | 8 ++++++++ 2 files changed, 22 insertions(+), 11 deletions(-) diff --git a/ipfx/time_series_utils.py b/ipfx/time_series_utils.py index 02491f08..f074cf1e 100644 --- a/ipfx/time_series_utils.py +++ b/ipfx/time_series_utils.py @@ -37,22 +37,25 @@ def calculate_dvdt(v, t, filter=None): dvdt : numpy array of time-derivative of voltage (V/s = mV/ms) """ - if has_fixed_dt(t) and filter: - delta_t = t[1] - t[0] - sample_freq = 1. / delta_t - filt_coeff = (filter * 1e3) / (sample_freq / 2.) # filter kHz -> Hz, then get fraction of Nyquist frequency - if filt_coeff < 0 or filt_coeff >= 1: - raise ValueError("bessel coeff ({:f}) is outside of valid range [0,1); cannot filter sampling frequency {:.1f} kHz with cutoff frequency {:.1f} kHz.".format(filt_coeff, sample_freq / 1e3, filter)) - b, a = signal.bessel(4, filt_coeff, "low") - v_filt = signal.filtfilt(b, a, v, axis=0) - dv = np.diff(v_filt) - else: + try: + if has_fixed_dt(t) and filter: + delta_t = t[1] - t[0] + sample_freq = 1. / delta_t + filt_coeff = (filter * 1e3) / (sample_freq / 2.) # filter kHz -> Hz, then get fraction of Nyquist frequency + if filt_coeff < 0 or filt_coeff >= 1: + raise ValueError("bessel coeff ({:f}) is outside of valid range [0,1); cannot filter sampling frequency {:.1f} kHz with cutoff frequency {:.1f} kHz.".format(filt_coeff, sample_freq / 1e3, filter)) + b, a = signal.bessel(4, filt_coeff, "low") + v_filt = signal.filtfilt(b, a, v, axis=0) + dv = np.diff(v_filt) + else: + dv = np.diff(v) + except ValueError: dv = np.diff(v) dt = np.diff(t) dvdt = 1e-3 * dv / dt # in V/s = mV/ms - # some data sources, such as neuron, occasionally report + # some data sources, such as neuron, occasionally report # duplicate timestamps, so we require that dt is not 0 return dvdt[np.fabs(dt) > sys.float_info.epsilon] diff --git a/tests/test_time_series_utils.py b/tests/test_time_series_utils.py index 1b2d2b76..28ba8f9d 100644 --- a/tests/test_time_series_utils.py +++ b/tests/test_time_series_utils.py @@ -19,6 +19,13 @@ def test_dvdt_no_filter(): assert np.allclose(tsu.calculate_dvdt(v, t), np.diff(v) / np.diff(t)) +def test_dvdt_with_filter_small_samplingrate(): + t = np.array([0, 1, 2, 3]) / 10000.0 + v = np.array([1, 1, 1, 1]) + + assert np.allclose(tsu.calculate_dvdt(v, t, 10), np.diff(v) / np.diff(t)) + + def test_fixed_dt(): t = [0, 1, 2, 3] assert tsu.has_fixed_dt(t) @@ -27,6 +34,7 @@ def test_fixed_dt(): t[0] -= 3. assert not tsu.has_fixed_dt(t) + def test_flatnotnan(): a = [1, 10, 12, 17, 13, 4, 8, np.nan, np.nan] From 2a270668dabc06e03ba415012b163c5b12c1330d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Bord=C3=A9=20S=C3=A1ndor?= Date: Wed, 18 Nov 2020 18:15:17 +0100 Subject: [PATCH 9/9] Fixed invalid troughs Some troughs fell onto the preceding peak, which caused an error at slicing. Before using it, drop those invalid putative spike features. --- ipfx/feature_extractor.py | 14 ++++++++++++++ 1 file changed, 14 insertions(+) diff --git a/ipfx/feature_extractor.py b/ipfx/feature_extractor.py index 08d4061c..c53aaf34 100644 --- a/ipfx/feature_extractor.py +++ b/ipfx/feature_extractor.py @@ -115,6 +115,20 @@ def process(self, t, v, i): # Spike list and thresholds have been refined - now find other features upstrokes = spkd.find_upstroke_indexes(v, t, thresholds, peaks, self.filter, dvdt) troughs = spkd.find_trough_indexes(v, t, thresholds, peaks, clipped, self.end) + not_nan = np.logical_not(np.logical_or(np.isnan(troughs), np.isnan(peaks))) + not_nan_idx = np.argwhere(not_nan) + valid_peak_tr_pair = not_nan_idx[peaks[not_nan_idx] < troughs[not_nan_idx]] + + troughs = troughs[valid_peak_tr_pair] + peaks = peaks[valid_peak_tr_pair] + upstrokes = upstrokes[valid_peak_tr_pair] + clipped = clipped[valid_peak_tr_pair] + thresholds = thresholds[valid_peak_tr_pair] + + if not thresholds.size: + # Save time if no spikes detected + return DataFrame() + downstrokes = spkd.find_downstroke_indexes(v, t, peaks, troughs, clipped, dvdt=dvdt) trough_details, clipped = spkf.analyze_trough_details(v, t, thresholds, peaks, clipped, self.end, dvdt=dvdt)