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/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() 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/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/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 05868a34..a333d9a2 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) @@ -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/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 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)