From 74ce76d307d08a0bb5a60e40f0bd9e7f0b1aafd5 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Wed, 19 May 2021 10:30:54 -0700 Subject: [PATCH 1/3] Trim trailing NaNs from sweeps This follows the pattern previously used for trailing zeros, trimming as data is loaded. Also reduces the logging of sweep integrity tags from warning to info, as these sweeps are now common and expected. --- ipfx/dataset/ephys_data_set.py | 7 ++++--- ipfx/qc_feature_extractor.py | 2 +- ipfx/stim_features.py | 1 + 3 files changed, 6 insertions(+), 4 deletions(-) diff --git a/ipfx/dataset/ephys_data_set.py b/ipfx/dataset/ephys_data_set.py index e49fcf79..0449b697 100644 --- a/ipfx/dataset/ephys_data_set.py +++ b/ipfx/dataset/ephys_data_set.py @@ -316,11 +316,12 @@ def get_sweep_data(self, sweep_number: int) -> Dict: response = sweep_data['response'] - nonzero = np.flatnonzero(response) - if len(nonzero) == 0: + valid_data = (response != 0) & ~np.isnan(response) + valid_i = np.flatnonzero(valid_data) + if len(valid_i) == 0: recording_end_idx = 0 else: - recording_end_idx = nonzero[-1] + 1 + recording_end_idx = valid_i[-1] + 1 sweep_data["response"] = response[:recording_end_idx] sweep_data["stimulus"] = sweep_data["stimulus"][:recording_end_idx] diff --git a/ipfx/qc_feature_extractor.py b/ipfx/qc_feature_extractor.py index 05868a34..46fde2fd 100644 --- a/ipfx/qc_feature_extractor.py +++ b/ipfx/qc_feature_extractor.py @@ -287,7 +287,7 @@ def sweep_qc_features(data_set): qc_features = current_clamp_sweep_qc_features(sweep, is_ramp) sweep_features.update(qc_features) else: - logging.warning("sweep {}: {}".format(sweep_num, tags)) + logging.info("sweep {}: {}".format(sweep_num, tags)) sweeps_features.append(sweep_features) diff --git a/ipfx/stim_features.py b/ipfx/stim_features.py index 01ff04c4..4045223d 100644 --- a/ipfx/stim_features.py +++ b/ipfx/stim_features.py @@ -32,6 +32,7 @@ def get_stim_characteristics(i, t, test_pulse=True): else: amplitude = float(peak_low) + # this is very unlikely, but possible for zero stim with mid-sweep NaNs if np.isnan(amplitude): return None, None, 0.0, None, None From 8d1b8bd75a9fef5386c885f8e7d00abcab76188d Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Wed, 19 May 2021 13:38:03 -0700 Subject: [PATCH 2/3] Remove recording epoch --- .../gallery/analysis_examples/lsq_analysis.py | 4 +--- .../analysis_examples/ramp_analysis.py | 9 +------- .../short_square_analysis.py | 4 +--- ipfx/data_set_features.py | 4 ---- ipfx/epochs.py | 21 ------------------- ipfx/qc_feature_extractor.py | 6 +++--- ipfx/script_utils.py | 4 ---- ipfx/sweep.py | 1 - tests/test_epochs.py | 19 ----------------- tests/test_sweep.py | 10 ++++----- 10 files changed, 11 insertions(+), 71 deletions(-) diff --git a/docs/gallery/analysis_examples/lsq_analysis.py b/docs/gallery/analysis_examples/lsq_analysis.py index 63f8b5cf..b1518cb1 100644 --- a/docs/gallery/analysis_examples/lsq_analysis.py +++ b/docs/gallery/analysis_examples/lsq_analysis.py @@ -36,9 +36,7 @@ ) long_square_sweeps = data_set.sweep_set(long_square_table.sweep_number) -# Select epoch corresponding to the actual recording from the sweeps -# and align sweeps so that the experiment would start at the same time -long_square_sweeps.select_epoch("recording") +# Align sweeps so that the experiment would start at the same time long_square_sweeps.align_to_start_of_epoch("experiment") # find the start and end time of the stimulus diff --git a/docs/gallery/analysis_examples/ramp_analysis.py b/docs/gallery/analysis_examples/ramp_analysis.py index af30480b..09ed755e 100644 --- a/docs/gallery/analysis_examples/ramp_analysis.py +++ b/docs/gallery/analysis_examples/ramp_analysis.py @@ -37,14 +37,7 @@ ) ramp_sweeps = data_set.sweep_set(ramp_table.sweep_number) -# Select epoch corresponding to the actual recording from the sweeps -# and align sweeps so that the experiment would start at the same time -ramp_sweeps.select_epoch("recording") -ramp_sweeps.align_to_start_of_epoch("experiment") - -# Select epoch corresponding to the actual recording from the sweeps -# and align sweeps so that the experiment would start at the same time -ramp_sweeps.select_epoch("recording") +# Align sweeps so that the experiment would start at the same time ramp_sweeps.align_to_start_of_epoch("experiment") # find the start and end time of the stimulus diff --git a/docs/gallery/analysis_examples/short_square_analysis.py b/docs/gallery/analysis_examples/short_square_analysis.py index 48a0208e..dd29de9a 100644 --- a/docs/gallery/analysis_examples/short_square_analysis.py +++ b/docs/gallery/analysis_examples/short_square_analysis.py @@ -35,9 +35,7 @@ ) short_square_sweeps = data_set.sweep_set(short_square_table.sweep_number) -# Select epoch corresponding to the actual recording from the sweeps -# and align sweeps so that the experiment would start at the same time -short_square_sweeps.select_epoch("recording") +# Align sweeps so that the experiment would start at the same time short_square_sweeps.align_to_start_of_epoch("experiment") # find the start and end time of the stimulus diff --git a/ipfx/data_set_features.py b/ipfx/data_set_features.py index 2d73aa0c..b49fdc6b 100644 --- a/ipfx/data_set_features.py +++ b/ipfx/data_set_features.py @@ -160,7 +160,6 @@ def extract_sweep_features(data_set, sweep_table): sweep_numbers = sorted(sweep_numbers) sweep_set = data_set.sweep_set(sweep_numbers) - sweep_set.select_epoch("recording") sweep_set.align_to_start_of_epoch("experiment") dp = detection_parameters(stimulus_name).copy() @@ -238,7 +237,6 @@ def extract_cell_long_square_features(data_set, subthresh_min_amp=None): logging.info("Assigned subthreshold minimum amplitude of %f.", subthresh_min_amp) lsq_sweeps = data_set.sweep_set(long_square_sweep_numbers) - lsq_sweeps.select_epoch("recording") lsq_sweeps.align_to_start_of_epoch("experiment") lsq_start, lsq_dur, _, _, _ = stf.get_stim_characteristics( @@ -278,7 +276,6 @@ def extract_cell_short_square_features(data_set): raise er.FeatureError("No short square sweeps available for feature extraction") ssq_sweeps = data_set.sweep_set(short_square_sweep_numbers) - ssq_sweeps.select_epoch("recording") ssq_sweeps.align_to_start_of_epoch("experiment") ssq_start, ssq_dur, _, _, _ = stf.get_stim_characteristics( @@ -312,7 +309,6 @@ def extract_cell_ramp_features(data_set): raise er.FeatureError("No ramp sweeps available for feature extraction") ramp_sweeps = data_set.sweep_set(ramp_sweep_numbers) - ramp_sweeps.select_epoch("recording") ramp_sweeps.align_to_start_of_epoch("experiment") ramp_start, ramp_dur, _, _, _ = stf.get_stim_characteristics(ramp_sweeps.sweeps[0].i, ramp_sweeps.sweeps[0].t) diff --git a/ipfx/epochs.py b/ipfx/epochs.py index f15e3b93..e4b6243f 100644 --- a/ipfx/epochs.py +++ b/ipfx/epochs.py @@ -51,27 +51,6 @@ def get_last_noise_epoch(idx1, hz): return idx1-int(NOISE_EPOCH * hz), idx1 -def get_recording_epoch(response): - """ - Detect response epoch defined as interval from start to the last non-nan value of the response - - Parameters - ---------- - response: float np.array - - Returns - ------- - start,end: int - indices of the epoch - """ - - if len(tsu.flatnotnan(response)) == 0: - end_idx = 0 - else: - end_idx = tsu.flatnotnan(response)[-1] - return 0, end_idx - - def get_sweep_epoch(response): """ Defined as interval including entire sweep diff --git a/ipfx/qc_feature_extractor.py b/ipfx/qc_feature_extractor.py index 46fde2fd..a333d9a2 100644 --- a/ipfx/qc_feature_extractor.py +++ b/ipfx/qc_feature_extractor.py @@ -303,8 +303,8 @@ def check_sweep_integrity(sweep, is_ramp): tags.append(F"{k} epoch is missing") if not is_ramp: - if sweep.epochs["recording"] and sweep.epochs["experiment"]: - if sweep.epochs["recording"][1] < sweep.epochs["experiment"][1]: + if sweep.epochs["sweep"] and sweep.epochs["experiment"]: + if sweep.epochs["sweep"][1] < sweep.epochs["experiment"][1]: tags.append("Recording stopped before completing the experiment epoch") if np.isnan(sweep.i).any(): @@ -353,7 +353,7 @@ def current_clamp_sweep_qc_features(sweep, is_ramp): # measure mean and rms of Vm at end of recording # do not check for ramps, because they do not have enough time to recover - _, rec_end_idx = ep.get_recording_epoch(voltage) + _, rec_end_idx = ep.get_sweep_epoch(voltage) if not is_ramp: idx0, idx1 = ep.get_last_stability_epoch(rec_end_idx, hz) diff --git a/ipfx/script_utils.py b/ipfx/script_utils.py index 8e03864f..a93cfced 100644 --- a/ipfx/script_utils.py +++ b/ipfx/script_utils.py @@ -179,7 +179,6 @@ def categorize_iclamp_sweeps(data_set, stimuli_names, sweep_qc_option="none", sp def validate_sweeps(data_set, sweep_numbers, extra_dur=0.2): check_sweeps = data_set.sweep_set(sweep_numbers) - check_sweeps.select_epoch("recording") valid_sweep_stim = [] start = None dur = None @@ -216,7 +215,6 @@ def preprocess_long_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, subthr lsq_sweeps, lsq_start, lsq_end = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur) if len(lsq_sweeps.sweeps) == 0: raise er.FeatureError("No long square sweeps were long enough or did not end early") - lsq_sweeps.select_epoch("recording") lsq_spx, lsq_spfx = dsf.extractors_for_sweeps( lsq_sweeps, @@ -239,7 +237,6 @@ def preprocess_short_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, spike ssq_sweeps, ssq_start, ssq_end = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur) if len(ssq_sweeps.sweeps) == 0: raise er.FeatureError("No short square sweeps were long enough or did not end early") - ssq_sweeps.select_epoch("recording") ssq_spx, ssq_spfx = dsf.extractors_for_sweeps(ssq_sweeps, est_window = [ssq_start, ssq_start + 0.001], @@ -258,7 +255,6 @@ def preprocess_ramp_sweeps(data_set, sweep_numbers): raise er.FeatureError("No ramp sweeps available for feature extraction") ramp_sweeps = data_set.sweep_set(sweep_numbers) - ramp_sweeps.select_epoch("recording") ramp_start, ramp_dur, _, _, _ = stf.get_stim_characteristics(ramp_sweeps.sweeps[0].i, ramp_sweeps.sweeps[0].t) ramp_spx, ramp_spfx = dsf.extractors_for_sweeps(ramp_sweeps, diff --git a/ipfx/sweep.py b/ipfx/sweep.py index 888d3b1f..f518251e 100644 --- a/ipfx/sweep.py +++ b/ipfx/sweep.py @@ -62,7 +62,6 @@ def detect_epochs(self): epoch_detectors = { "sweep": ep.get_sweep_epoch(self.response), - "recording": ep.get_recording_epoch(self.response), "experiment": ep.get_experiment_epoch(self._i, self.sampling_rate,test_pulse), "stim": ep.get_stim_epoch(self.stimulus, test_pulse), } diff --git a/tests/test_epochs.py b/tests/test_epochs.py index 6364f0dc..b2a5d120 100644 --- a/tests/test_epochs.py +++ b/tests/test_epochs.py @@ -116,25 +116,6 @@ def test_get_stim_epoch(i, stim_epoch): assert stim_epoch == ep.get_stim_epoch(i) -@pytest.mark.parametrize('response,' - 'recording_epoch', - [ - ( - [0, 0, 1, 1.5, 0, 0, 2, 3, 4, 1, np.nan, np.nan], - (0, 9) - ), - - # zero array - ( - [np.nan, np.nan, np.nan, np.nan, np.nan, np.nan], - (0, 0) - ), - - ] - ) -def test_get_recording_epoch(response, recording_epoch): - assert recording_epoch == ep.get_recording_epoch(response) - @pytest.mark.parametrize('response,' 'sweep_epoch', diff --git a/tests/test_sweep.py b/tests/test_sweep.py index b5a5ea58..ef20436c 100644 --- a/tests/test_sweep.py +++ b/tests/test_sweep.py @@ -6,8 +6,8 @@ def sweep(): i = [0,0,1,1,0,0,0,2,2,2,2,2,0,0,0,0] - v = [0,0,1,2,1,0,0,1,2,3,1,np.nan,np.nan,np.nan,np.nan,np.nan] - sampling_rate = 2 + v = [0,0,1,2,1,0,0,1,2,3,1,1,np.nan,np.nan,np.nan,np.nan] + sampling_rate = 10 dt = 1./sampling_rate t = np.arange(0,len(v))*dt @@ -19,9 +19,9 @@ def test_select_epoch(sweep): i_sweep = sweep.i v_sweep = sweep.v - sweep.select_epoch("recording") - assert np.all(sweep.i == [0,0,1,1,0,0,0,2,2,2,2]) - assert np.all(sweep.v == [0,0,1,2,1,0,0,1,2,3,1]) + sweep.select_epoch("stim") + assert np.all(sweep.i == [2,2,2,2,2]) + assert np.all(sweep.v == [1,2,3,1,1]) sweep.select_epoch("sweep") assert np.all(sweep.i == i_sweep) From 68168b18723579dcb55cbdd35e0623d55390ca01 Mon Sep 17 00:00:00 2001 From: Tom Chartrand Date: Wed, 19 May 2021 14:51:57 -0700 Subject: [PATCH 3/3] improved plot_data_set --- ipfx/bin/plot_ephys_nwb.py | 87 +++++++++++++++++++++++++++----------- 1 file changed, 62 insertions(+), 25 deletions(-) diff --git a/ipfx/bin/plot_ephys_nwb.py b/ipfx/bin/plot_ephys_nwb.py index 10a9c332..0a816c36 100755 --- a/ipfx/bin/plot_ephys_nwb.py +++ b/ipfx/bin/plot_ephys_nwb.py @@ -1,34 +1,68 @@ import sys +import warnings import matplotlib.pyplot as plt import numpy as np from ipfx.dataset.create import create_ephys_data_set +from ipfx.qc_feature_extractor import sweep_qc_features +from ipfx.utilities import drop_failed_sweeps from ipfx.stimulus import StimulusOntology import allensdk.core.json_utilities as ju - - -def plot_data_set(data_set, sweep_table, nwb_file_name): - +from typing import ( + Optional, List, Dict, Tuple, Collection, Sequence, Union +) + + +def plot_data_set(data_set, + clamp_mode: Optional[str] = None, + stimuli: Optional[Collection[str]] = None, + stimuli_exclude: Optional[Collection[str]] = None, + show_amps: Optional[bool] = True, + qc_sweeps: Optional[bool] = True, + figsize=(15, 7), + ): + nwb_file_name = str(data_set._data.nwb_file) + if qc_sweeps: + drop_failed_sweeps(data_set) + elif show_amps: + data_set.sweep_info = sweep_qc_features(data_set) + + sweep_table = data_set.filtered_sweep_table(clamp_mode=clamp_mode, stimuli=stimuli, stimuli_exclude=stimuli_exclude) + + if len(sweep_table)==0: + warnings.warn("No sweeps to plot") + return + height_ratios, width_ratios = axes_ratios(sweep_table) fig, ax = plt.subplots(len(height_ratios), 3, - figsize=(15, 7), + figsize=figsize, gridspec_kw={'height_ratios': height_ratios, 'width_ratios': width_ratios} ) + if len(height_ratios)==1: + # ensure 2d array + ax = ax[None, :] for fig_row, (stimulus_code, sweep_set_table) in enumerate(sweep_table.groupby("stimulus_code")): - sweep_numbers = sweep_set_table["sweep_number"].sort_values().values + sweep_set_table = sweep_set_table.copy().sort_values("sweep_number", ascending=False) + sweep_numbers = sweep_set_table["sweep_number"] ss = data_set.sweep_set(sweep_numbers) + if qc_sweeps: + ss.select_epoch('experiment') + annot = sweep_numbers.astype(str) + if show_amps: + annot += sweep_set_table['stimulus_amplitude'].apply(": {:.3g} pA".format) + ax_a = ax[fig_row,0] ax_i = ax[fig_row,1] ax_v = ax[fig_row,2] - plot_waveforms(ax_i, ss.i, ss.sampling_rate, ss.sweep_number) - plot_waveforms(ax_v, ss.v, ss.sampling_rate, ss.sweep_number) + plot_waveforms(ax_i, ss.i, ss.sampling_rate, annot) + plot_waveforms(ax_v, ss.v, ss.sampling_rate) ax_v.get_shared_x_axes().join(ax_i, ax_v) clamp_mode = sweep_set_table["clamp_mode"].values[0] - ax_a.text(0, 0.0, "%s, %s " % (stimulus_code, clamp_mode)) + ax_a.text(0, 0.0, "%s \n%s " % (stimulus_code, clamp_mode)) ax_a.axis('off') ax[0,0].set_title("Description") @@ -40,14 +74,15 @@ def plot_data_set(data_set, sweep_table, nwb_file_name): fig.suptitle("file: " + nwb_file_name, fontsize=12) mng = plt.get_current_fig_manager() - mng.resize(*mng.window.maxsize()) + if hasattr(mng, 'window'): + mng.resize(*mng.window.maxsize()) plt.subplots_adjust(left=0.01, right=0.98, bottom=0.02,top=0.92) def axes_ratios(sweep_table): height_ratios = [] - width_ratios = [1,3,3] + width_ratios = [1,4,4] for _, sweep_set_table in sweep_table.groupby("stimulus_code"): height_ratios.append(len(sweep_set_table.index)) @@ -55,26 +90,32 @@ def axes_ratios(sweep_table): return height_ratios, width_ratios -def plot_waveforms(ax, ys, rs, sn): +def plot_waveforms(ax, ys, rs, annotations=None): offset = 0 - - for y, r, sn in zip(ys, rs, sn): - - offset += get_vertical_offset(y) + dy = 0 + for i, (y, r) in enumerate(zip(ys, rs)): + if len(y)==0: + continue + y -= y[0] + dy = max(dy, get_vertical_offset(y)) + offset += dy y += offset x = np.arange(0, len(y)) / r ax.plot(x, y) - ax.text(x[0] - 0.05 * (x[-1] - x[0]), y[0], str(sn), fontsize=8) - ax.set_xlim(x[0], x[-1]) + if annotations is not None: + ax.text(x[0] - 0.01 * (x[-1] - x[0]), y[0], annotations.iloc[i], fontsize=8, ha='right') + # need to set limits to show all sweeps + # mode would be best but mean will do fine + ax.set_xlim(0, np.max([len(y) for y in ys])/r) customize_axis(ax) def customize_axis(ax): - ax.tick_params(axis="x", direction="in", pad=-10,labelsize=8) + ax.tick_params(axis="x", direction="in", pad=3, labelsize=8) ax.get_yaxis().set_visible(False) ax.spines['top'].set_visible(False) ax.spines['right'].set_visible(False) @@ -102,12 +143,8 @@ def main(): ont = StimulusOntology(ju.read(stimulus_ontology_file)) data_set = create_ephys_data_set(nwb_file=nwb_file) - vclamp_sweep_table = data_set.filtered_sweep_table(clamp_mode=data_set.VOLTAGE_CLAMP) - plot_data_set(data_set, vclamp_sweep_table, nwb_file) - - data_set = create_ephys_data_set(nwb_file=nwb_file) - iclamp_sweep_table = data_set.filtered_sweep_table(clamp_mode=data_set.CURRENT_CLAMP) - plot_data_set(data_set, iclamp_sweep_table, nwb_file) + plot_data_set(data_set, clamp_mode=data_set.VOLTAGE_CLAMP) + plot_data_set(data_set, clamp_mode=data_set.CURRENT_CLAMP) plt.show()