diff --git a/ipfx/bin/run_feature_collection.py b/ipfx/bin/run_feature_collection.py index ec265e0c..9101a023 100755 --- a/ipfx/bin/run_feature_collection.py +++ b/ipfx/bin/run_feature_collection.py @@ -24,29 +24,44 @@ class CollectFeatureParameters(ags.ArgSchema): default="lims", validate=lambda x: x in ["lims", "filesystem"] ) - - -def data_for_specimen_id(specimen_id, passed_only, data_source, ontology, file_list=None): + sweep_qc_option = ags.fields.String( + description=("Sweep-level QC option - " + "'none': use all sweeps; " + "'passed-only': only use passed sweeps; " + "'passed-except-delta-vm': also use sweeps whose only failure is delta Vm; " + "'passed-except-delta-vm-and-rms': as above but also re-check RMS"), + default='none' + ) + + +def data_for_specimen_id(specimen_id, sweep_qc_option, data_source, ontology, + sweep_qc_record, file_list=None): data_set = su.dataset_for_specimen_id(specimen_id, data_source, ontology, file_list) if type(data_set) is dict and "error" in data_set: logging.warning("Problem getting AibsDataSet for specimen {:d} from LIMS".format(specimen_id)) return {} try: - lsq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, ontology.long_square_names) - ssq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, ontology.short_square_names) - ramp_sweep_numbers = su.categorize_iclamp_sweeps(data_set, ontology.ramp_names) + lsq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, + ontology.long_square_names, sweep_qc_option=sweep_qc_option, + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) + ssq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, + ontology.short_square_names, sweep_qc_option=sweep_qc_option, + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) + ramp_sweep_numbers = su.categorize_iclamp_sweeps(data_set, + ontology.ramp_names, sweep_qc_option=sweep_qc_option, + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) except Exception as detail: - logging.warn("Exception when processing specimen {:d}".format(specimen_id)) - logging.warn(detail) + logging.warning("Exception when processing specimen {:d}".format(specimen_id)) + logging.warning(detail) # return {"error": {"type": "sweep_table", "details": traceback.format_exc(limit=1)}} return {} try: result = extract_features(data_set, ramp_sweep_numbers, ssq_sweep_numbers, lsq_sweep_numbers) except Exception as detail: - logging.warn("Exception when processing specimen {:d}".format(specimen_id)) - logging.warn(detail) + logging.warning("Exception when extracting features for specimen {:d}".format(specimen_id)) + logging.warning(detail) # return {"error": {"type": "processing", "details": traceback.format_exc(limit=1)}} return {} @@ -77,8 +92,7 @@ def extract_features(data_set, ramp_sweep_numbers, ssq_sweep_numbers, lsq_sweep_ (lsq_sweeps, basic_lsq_features, lsq_an, - lsq_start, - lsq_end) = su.preprocess_long_square_sweeps(data_set, lsq_sweep_numbers) + lsq_stim_timing) = su.preprocess_long_square_sweeps(data_set, lsq_sweep_numbers) features.update({ "input_resistance": basic_lsq_features["input_resistance"], @@ -234,7 +248,8 @@ def lin_sqrt_fit(x, y): def run_feature_collection(ids=None, project="T301", include_failed_sweeps=True, include_failed_cells=False, - output_file="", run_parallel=True, data_source="lims", file_list=None, **kwargs): + output_file="", run_parallel=True, data_source="lims", file_list=None, + sweep_qc_option="none", **kwargs): if ids is not None: specimen_ids = ids else: @@ -242,11 +257,25 @@ def run_feature_collection(ids=None, project="T301", include_failed_sweeps=True, logging.info("Number of specimens to process: {:d}".format(len(specimen_ids))) + # Build the sweep QC record used by sweep categorization. When sweep-level QC + # is being applied against LIMS, query the record; otherwise an empty record + # (with the expected columns) is sufficient. + sweep_qc_record_df = kwargs.get("sweep_qc_record_df", None) + if sweep_qc_record_df is None: + if sweep_qc_option != "none" and data_source == "lims": + sweep_qc_record = lq.get_sweep_states_and_tags_for_specimens(specimen_ids) + sweep_qc_record_df = pd.DataFrame(sweep_qc_record) + sweep_qc_record_df["tag_name"] = sweep_qc_record_df["tag_name"].fillna("None") + else: + sweep_qc_record_df = pd.DataFrame( + columns=["specimen_id", "sweep_number", "workflow_state", "tag_name"]) + ontology = StimulusOntology(ju.read(StimulusOntology.DEFAULT_STIMULUS_ONTOLOGY_FILE)) get_data_partial = partial(data_for_specimen_id, - passed_only=not include_failed_sweeps, + sweep_qc_option=sweep_qc_option, data_source=data_source, ontology=ontology, + sweep_qc_record=sweep_qc_record_df, file_list=file_list) if run_parallel: diff --git a/ipfx/bin/run_feature_vector_extraction.py b/ipfx/bin/run_feature_vector_extraction.py index af967650..d1b62d44 100755 --- a/ipfx/bin/run_feature_vector_extraction.py +++ b/ipfx/bin/run_feature_vector_extraction.py @@ -1,4 +1,5 @@ import numpy as np +import pandas as pd import argschema as ags import logging import traceback @@ -73,6 +74,7 @@ def data_for_specimen_id( sweep_qc_option, data_source, ontology, + sweep_qc_record, ap_window_length=0.005, target_sampling_rate=50000, file_list=None, @@ -113,12 +115,15 @@ def data_for_specimen_id( try: lsq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, ontology.long_square_names, sweep_qc_option=sweep_qc_option, - specimen_id=specimen_id) + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) (lsq_sweeps, lsq_features, _, - lsq_start, - lsq_end) = su.preprocess_long_square_sweeps(data_set, lsq_sweep_numbers) + lsq_stim_timing) = su.preprocess_long_square_sweeps(data_set, lsq_sweep_numbers) + + # Create stimulus timing dictionary keyed on sweep number + lsq_stim_timing_dict = {lsq_sweeps.sweeps[i].sweep_number: lsq_stim_timing[i] + for i in range(len(lsq_stim_timing))} except Exception as detail: logging.warning("Exception when preprocessing long square sweeps from specimen {:d}".format(specimen_id)) @@ -129,7 +134,7 @@ def data_for_specimen_id( try: ssq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, ontology.short_square_names, sweep_qc_option=sweep_qc_option, - specimen_id=specimen_id) + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) ssq_sweeps, ssq_features, _ = su.preprocess_short_square_sweeps(data_set, ssq_sweep_numbers) except Exception as detail: @@ -141,7 +146,7 @@ def data_for_specimen_id( try: ramp_sweep_numbers = su.categorize_iclamp_sweeps(data_set, ontology.ramp_names, sweep_qc_option=sweep_qc_option, - specimen_id=specimen_id) + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) ramp_sweeps, ramp_features, _ = su.preprocess_ramp_sweeps(data_set, ramp_sweep_numbers) except Exception as detail: @@ -162,20 +167,19 @@ def data_for_specimen_id( target_amps_for_step_subthresh = [-90, -70, -50, -30, -10] result["step_subthresh"] = fv.step_subthreshold( subthresh_hyperpol_dict, target_amps_for_step_subthresh, - lsq_start, lsq_end, amp_tolerance=5) + lsq_stim_timing_dict, amp_tolerance=5) result["subthresh_norm"] = fv.subthresh_norm(subthresh_hyperpol_dict, hyperpol_deflect_dict, - lsq_start, lsq_end) + lsq_stim_timing_dict) (subthresh_depol_dict, depol_deflect_dict) = fv.identify_subthreshold_depol_with_amplitudes(lsq_features, lsq_sweeps) result["subthresh_depol_norm"] = fv.subthresh_depol_norm( subthresh_depol_dict, depol_deflect_dict, - np.round(lsq_start, decimals=3), - np.round(lsq_end, decimals=3)) + lsq_stim_timing_dict) isi_sweep, isi_sweep_spike_info = fv.identify_sweep_for_isi_shape( - lsq_sweeps, lsq_features, lsq_end - lsq_start) - result["isi_shape"] = fv.isi_shape(isi_sweep, isi_sweep_spike_info, lsq_end) + lsq_sweeps, lsq_features, lsq_stim_timing_dict) + result["isi_shape"] = fv.isi_shape(isi_sweep, isi_sweep_spike_info, lsq_stim_timing_dict) # Calculate waveforms from each type of sweep spiking_ssq_sweep_list = [ssq_sweeps.sweeps[swp_ind] @@ -210,10 +214,15 @@ def data_for_specimen_id( result["first_ap_dv"] = np.hstack([ssq_ap_dv, lsq_ap_dv, ramp_ap_dv]) target_amplitudes = np.arange(0, 120, 20) - supra_info_list = fv.identify_suprathreshold_spike_info( - lsq_features, target_amplitudes, shift=10) - result["psth"] = fv.psth_vector(supra_info_list, lsq_start, lsq_end) - result["inst_freq"] = fv.inst_freq_vector(supra_info_list, lsq_start, lsq_end) + supra_info_list, supra_sweep_numbers = fv.identify_suprathreshold_spike_info( + lsq_features, target_amplitudes, + sweep_numbers=[swp.sweep_number for swp in lsq_sweeps.sweeps], + shift=10) + supra_lsq_stim_timing_list = [ + lsq_stim_timing_dict[sn] if sn is not None else None + for sn in supra_sweep_numbers] + result["psth"] = fv.psth_vector(supra_info_list, supra_lsq_stim_timing_list) + result["inst_freq"] = fv.inst_freq_vector(supra_info_list, supra_lsq_stim_timing_list) spike_feature_list = [ "upstroke_downstroke_ratio", @@ -224,7 +233,7 @@ def data_for_specimen_id( ] for feature in spike_feature_list: result["spiking_" + feature] = fv.spike_feature_vector(feature, - supra_info_list, lsq_start, lsq_end) + supra_info_list, supra_lsq_stim_timing_list) except Exception as detail: logging.warning("Exception when processing specimen {:d}".format(specimen_id)) logging.warning(detail) @@ -294,11 +303,25 @@ def run_feature_vector_extraction( ontology = StimulusOntology(ju.read(StimulusOntology.DEFAULT_STIMULUS_ONTOLOGY_FILE)) + # Build the sweep QC record used by sweep categorization. When sweep-level QC + # is being applied against LIMS, query the record; otherwise an empty record + # (with the expected columns) is sufficient. + sweep_qc_record_df = kwargs.get("sweep_qc_record_df", None) + if sweep_qc_record_df is None: + if sweep_qc_option != "none" and data_source == "lims": + sweep_qc_record = lq.get_sweep_states_and_tags_for_specimens(specimen_ids) + sweep_qc_record_df = pd.DataFrame(sweep_qc_record) + sweep_qc_record_df["tag_name"] = sweep_qc_record_df["tag_name"].fillna("None") + else: + sweep_qc_record_df = pd.DataFrame( + columns=["specimen_id", "sweep_number", "workflow_state", "tag_name"]) + logging.info("Number of specimens to process: {:d}".format(len(specimen_ids))) get_data_partial = partial(data_for_specimen_id, sweep_qc_option=sweep_qc_option, data_source=data_source, ontology=ontology, + sweep_qc_record=sweep_qc_record_df, ap_window_length=ap_window_length, file_list=file_list) diff --git a/ipfx/bin/run_feature_vector_extraction_flex.py b/ipfx/bin/run_feature_vector_extraction_flex.py new file mode 100644 index 00000000..3b395ca1 --- /dev/null +++ b/ipfx/bin/run_feature_vector_extraction_flex.py @@ -0,0 +1,489 @@ +import logging +import os +import json +import traceback + +import argschema as ags +import numpy as np +import pandas as pd + +import ipfx.lims_queries as lq +import ipfx.json_utilities as ju +import ipfx.script_utils as su +import ipfx.feature_vectors as fv + +from concurrent.futures import ProcessPoolExecutor +from functools import partial + +from ipfx.dataset.create import create_ephys_data_set + + +class StartEndDurationSchema(ags.schemas.DefaultSchema): + before = ags.fields.Float( + description="duration to extend before stimulus", + default=0.2, + ) + after = ags.fields.Float( + description="duration to extend after stimulus", + default=0.2, + ) + + +class ExtendDurationSchema(ags.schemas.DefaultSchema): + step_subthresh = ags.fields.Nested(StartEndDurationSchema, + required=True, + default={"before": 0.2, "after": 0.2}, + description="parameters for extending duration around step subthreshold analysis", + ) + subthresh_norm = ags.fields.Nested(StartEndDurationSchema, + required=True, + default={"before": 0.2, "after": 0.2}, + description="parameters for extending duration around normalized subthreshold analysis", + ) + + +class ApWaveformSchema(ags.schemas.DefaultSchema): + use = ags.fields.Boolean( + default=True, + description="whether to use AP from stimulus type", + ) + duration = ags.fields.Float( + default=0.003, + description="Duration after threshold for AP shape (s)", + ) + + +class ApWaveformForStimuliSchema(ags.schemas.DefaultSchema): + ssq = ags.fields.Nested(ApWaveformSchema, + default={}, + description="analysis parameters for short square AP waveform", + ) + lsq = ags.fields.Nested(ApWaveformSchema, + default={}, + description="analysis parameters for long square AP waveform", + ) + ramp = ags.fields.Nested(ApWaveformSchema, + default={}, + description="analysis parameters for ramp AP waveform", + ) + + +class CollectFeatureVectorParameters(ags.ArgSchema): + output_dir = ags.fields.OutputDir( + description="Destination directory for output files", + default="." + ) + output_code = ags.fields.String( + description="Code used for naming of output files", + default="test" + ) + specimen_id_file = ags.fields.InputFile( + description=("Input file of specimen IDs (one per line)"), + ) + nwb_path_file = ags.fields.InputFile( + description=("JSON file with paths to each specimen's NWB file - " + "if not supplied, LIMS will be queried for them"), + default=None, + allow_none=True, + ) + sweep_qc_record_file = ags.fields.InputFile( + description=("File with sweep QC status and tags - " + "if not supplied, LIMS will be queried for them"), + default=None, + allow_none=True, + ) + manual_fail_sweep_file = ags.fields.InputFile( + description=("File with manual sweep failure information"), + default=None, + allow_none=True, + ) + sweep_qc_option = ags.fields.String( + description=("Sweep-level QC option - " + "'none': use all sweeps; " + "'passed-only': check passed status with LIMS and " + "only used passed sweeps " + "'passed-except-delta-vm': check status with LIMS and " + "use passed sweeps and sweeps where only failure criterion is delta_vm" + "'passed-except-delta-vm-and-rms': check status with LIMS and " + "use passed sweeps and sweeps where only failure criterion is delta_vm," + "but also re-calculate RMS values with current code" + ), + default='none' + ) + extract_from_ramp = ags.fields.Boolean( + description="whether to run analysis on ramp sweep", + default=True, + ) + amp_tolerance = ags.fields.Float( + description="how much deviation from expected stimulus amplitudes is acceptable (in pA)", + default=4. + ) + additional_fvs = ags.fields.List(ags.fields.String, + allow_none=True, + default=[], + cli_as_single_argument=True, + ) + extend_durations = ags.fields.Nested(ExtendDurationSchema, + description="parameters for extending time windows for analyses", + default={}, + ) + ap_waveforms = ags.fields.Nested(ApWaveformForStimuliSchema, + default={}, + description="parameters for AP waveform analysis", + ) + needed_amplitudes = ags.fields.List( + ags.fields.Integer, + allow_none=True, + default=None, + cli_as_single_argument=True + ) + run_parallel = ags.fields.Boolean( + description="boolean - use multiprocessing", + default=True + ) + + + +def data_for_specimen_id( + specimen_id, + sweep_qc_option, + sweep_qc_record, + file_list, + ap_waveforms, + extend_durations, + extract_from_ramp, + additional_fvs, + target_sampling_rate=50000, + needed_amplitudes=None, + amp_tolerance=0., + manual_fail_sweeps=None, +): + """ + Extract feature vector from given cell identified by the specimen_id + Parameters + ---------- + specimen_id : int + cell identified + sweep_qc_option : str + see CollectFeatureVectorParameters input schema for details + sweep_qc_record: DataFrame + sweep status and error tag dataframe + data_source: str + see CollectFeatureVectorParameters input schema for details + target_sampling_rate : float + sampling rate + file_list : list of str + nwbfile names + Returns + ------- + dict : + features for a given cell specimen_id + + """ + logging.debug(f"Starting to process specimen id: {specimen_id}") + + try: + data_set = create_ephys_data_set(nwb_file=file_list[specimen_id]) + except Exception as detail: + logging.warning("Exception when creating data set for specimen {:d}".format(specimen_id)) + logging.warning(detail) + return {"error": {"type": "data_set", "details": traceback.format_exc(limit=None)}, "specimen_id": specimen_id} + + # Identify and preprocess long square sweeps + try: + lsq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, + data_set.ontology.long_square_names, sweep_qc_option=sweep_qc_option, + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) + + if manual_fail_sweeps is not None and specimen_id in manual_fail_sweeps: + lsq_sweep_numbers = np.array([sn for sn in lsq_sweep_numbers if sn not in manual_fail_sweeps[specimen_id]]) + + (lsq_sweeps, + lsq_features, + _, + lsq_stim_timing) = su.preprocess_long_square_sweeps(data_set, lsq_sweep_numbers) + + # Create stimulus timing dictionary keyed on sweep number + lsq_stim_timing_dict = {lsq_sweeps.sweeps[i].sweep_number: lsq_stim_timing[i] + for i in range(len(lsq_stim_timing))} + except Exception as detail: + logging.warning("Exception when preprocessing long square sweeps from specimen {:d}".format(specimen_id)) + logging.warning(detail) + return {"error": {"type": "sweep_table", "details": traceback.format_exc(limit=None)}, "specimen_id": specimen_id} + + + # Identify and preprocess short square sweeps + try: + ssq_sweep_numbers = su.categorize_iclamp_sweeps(data_set, + data_set.ontology.short_square_names, sweep_qc_option=sweep_qc_option, + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) + + if manual_fail_sweeps is not None and specimen_id in manual_fail_sweeps: + ssq_sweep_numbers = np.array([sn for sn in ssq_sweep_numbers if sn not in manual_fail_sweeps[specimen_id]]) + + ssq_sweeps, ssq_features, _ = su.preprocess_short_square_sweeps(data_set, + ssq_sweep_numbers) + except Exception as detail: + logging.warning("Exception when preprocessing short square sweeps from specimen {:d}".format(specimen_id)) + logging.warning(detail) + return {"error": {"type": "sweep_table", "details": traceback.format_exc(limit=None)}, "specimen_id": specimen_id} + + # Identify and preprocess ramp sweeps + if extract_from_ramp: + logging.debug("Identifying and processing ramp sweeps") + try: + ramp_sweep_numbers = su.categorize_iclamp_sweeps(data_set, + data_set.ontology.ramp_names, sweep_qc_option=sweep_qc_option, + specimen_id=specimen_id, sweep_qc_record=sweep_qc_record) + if manual_fail_sweeps is not None and specimen_id in manual_fail_sweeps: + ramp_sweep_numbers = np.array([sn for sn in ramp_sweep_numbers if sn not in manual_fail_sweeps[specimen_id]]) + ramp_sweeps, ramp_features, _ = su.preprocess_ramp_sweeps(data_set, + ramp_sweep_numbers) + except Exception as detail: + logging.warning("Exception when preprocessing ramp sweeps from specimen {:d}".format(specimen_id)) + logging.warning(detail) + return {"error": {"type": "sweep_table", "details": traceback.format_exc(limit=None), "specimen_id": specimen_id}} + + # Calculate desired feature vectors + result = {"id": [specimen_id]} # list because of how the results are accumulated for output + + try: + (subthresh_hyperpol_dict, + hyperpol_deflect_dict) = fv.identify_subthreshold_hyperpol_with_amplitudes(lsq_features, + lsq_sweeps) + target_amps_for_step_subthresh = [-90, -70, -50, -30, -10] + result["step_subthresh"] = fv.step_subthreshold( + subthresh_hyperpol_dict, target_amps_for_step_subthresh, + lsq_stim_timing_dict, amp_tolerance=amp_tolerance, + extend_duration_before=extend_durations["step_subthresh"]["before"], + extend_duration_after=extend_durations["step_subthresh"]["after"], + ) + result["subthresh_norm"] = fv.subthresh_norm(subthresh_hyperpol_dict, hyperpol_deflect_dict, + lsq_stim_timing_dict, + extend_duration_before=extend_durations["subthresh_norm"]["before"], + extend_duration_after=extend_durations["subthresh_norm"]["after"], + ) + if "subthresh_rebound" in additional_fvs: + result["subthresh_rebound"] = fv.subthresh_rebound( + subthresh_hyperpol_dict, + lsq_stim_timing_dict, dur=0.3, + ) + + (subthresh_depol_dict, + depol_deflect_dict) = fv.identify_subthreshold_depol_with_amplitudes(lsq_features, + lsq_sweeps) + result["subthresh_depol_norm"] = fv.subthresh_depol_norm(subthresh_depol_dict, + depol_deflect_dict, lsq_stim_timing_dict) + isi_sweep, isi_sweep_spike_info = fv.identify_sweep_for_isi_shape( + lsq_sweeps, lsq_features, lsq_stim_timing_dict) + result["isi_shape"] = fv.isi_shape(isi_sweep, isi_sweep_spike_info, lsq_stim_timing_dict) + + if result["isi_shape"] is None: + # Failed to calculate a shape for the first value; try other sweeps + exclude_sweeps_for_isi = [] + while result["isi_shape"] is None: + exclude_sweeps_for_isi.append(isi_sweep.sweep_number) + isi_sweep, isi_sweep_spike_info = fv.identify_sweep_for_isi_shape( + lsq_sweeps, lsq_features, lsq_stim_timing_dict, exclude_sweep_numbers=exclude_sweeps_for_isi) + result["isi_shape"] = fv.isi_shape(isi_sweep, isi_sweep_spike_info, lsq_stim_timing_dict) + + + # Calculate waveforms from each type of sweep - if multiple sweeps, use the earliest + ap_v_list = [] + ap_dv_list = [] + + if ap_waveforms["ssq"]["use"]: + spiking_ssq_sweep_list = [ssq_sweeps.sweeps[swp_ind] + for swp_ind in ssq_features["common_amp_sweeps"].index] + spiking_ssq_info_list = [ssq_features["spikes_set"][swp_ind] + for swp_ind in ssq_features["common_amp_sweeps"].index] + ssq_ap_v, ssq_ap_dv = fv.first_ap_vectors(spiking_ssq_sweep_list[:1], + spiking_ssq_info_list[:1], + target_sampling_rate=target_sampling_rate, + window_length=ap_waveforms["ssq"]["duration"], + skip_clipped=True) + ap_v_list.append(ssq_ap_v) + ap_dv_list.append(ssq_ap_dv) + + if ap_waveforms["lsq"]["use"]: + rheo_ind = lsq_features["rheobase_sweep"].name + sweep = lsq_sweeps.sweeps[rheo_ind] + lsq_ap_v, lsq_ap_dv = fv.first_ap_vectors([sweep], + [lsq_features["spikes_set"][rheo_ind]], + target_sampling_rate=target_sampling_rate, + window_length=ap_waveforms["lsq"]["duration"]) + ap_v_list.append(lsq_ap_v) + ap_dv_list.append(lsq_ap_dv) + + if extract_from_ramp and ap_waveforms["ramp"]["use"]: + spiking_ramp_sweep_list = [ramp_sweeps.sweeps[swp_ind] + for swp_ind in ramp_features["spiking_sweeps"].index] + spiking_ramp_info_list = [ramp_features["spikes_set"][swp_ind] + for swp_ind in ramp_features["spiking_sweeps"].index] + ramp_ap_v, ramp_ap_dv = fv.first_ap_vectors(spiking_ramp_sweep_list[:1], + spiking_ramp_info_list[:1], + target_sampling_rate=target_sampling_rate, + window_length=ap_waveforms["ramp"]["duration"], + skip_clipped=True) + ap_v_list.append(ramp_ap_v) + ap_dv_list.append(ramp_ap_dv) + + # Combine so that differences can be assessed by analyses like sPCA + result["first_ap_v"] = np.hstack(ap_v_list) + result["first_ap_dv"] = np.hstack(ap_dv_list) + + target_amplitudes = np.arange(0, 100, 10) + supra_info_list, supra_sweep_numbers = fv.identify_suprathreshold_spike_info( + lsq_features, + target_amplitudes, + sweep_numbers=[swp.sweep_number for swp in lsq_sweeps.sweeps], + shift=None, + amp_tolerance=amp_tolerance, + needed_amplitudes=needed_amplitudes + ) + + supra_lsq_stim_timing_list = [lsq_stim_timing_dict[sn] if sn is not None else None for sn in supra_sweep_numbers] + + actual_amps = [int(a) for a, si in zip(target_amplitudes, supra_info_list) if si is not None] + actual_rheobase_i = int(lsq_features["rheobase_i"]) + + result["long_squares_data_info"] = {"rheobase_i": actual_rheobase_i, "amplitudes_with_data": actual_amps} + + result["psth"] = fv.psth_vector(supra_info_list, supra_lsq_stim_timing_list) + result["inst_freq"] = fv.inst_freq_vector(supra_info_list, supra_lsq_stim_timing_list) + + spike_feature_list = [ + "upstroke_downstroke_ratio", + "peak_v", + "fast_trough_v", + "threshold_v", + "width", + ] + for feature in spike_feature_list: + result["spiking_" + feature] = fv.spike_feature_vector(feature, + supra_info_list, supra_lsq_stim_timing_list) + except Exception as detail: + logging.warning("Exception when processing specimen {:d}".format(specimen_id)) + logging.warning(detail) + return {"error": {"type": "processing", "details": traceback.format_exc(limit=None)}, "specimen_id": specimen_id} + + logging.info(f"Successfully processed {specimen_id}") + + # Flush the LRU cache for the data_set object + if hasattr(data_set, "_data") and hasattr(data_set._data, "_get_series"): + data_set._data._get_series.cache_clear() + + return result + + +def run_feature_vector_extraction( + specimen_ids, + output_dir, + output_code, + sweep_qc_option, + file_list, + sweep_qc_record_df, + manual_fail_sweep_dict, + extract_from_ramp, + amp_tolerance, + additional_fvs, + extend_durations, + ap_waveforms, + needed_amplitudes, + run_parallel=True, + ): + """ + Extract feature vectors from a list of cells and save results + """ + + get_data_partial = partial(data_for_specimen_id, + sweep_qc_option=sweep_qc_option, + needed_amplitudes=needed_amplitudes, + amp_tolerance=amp_tolerance, + ap_waveforms=ap_waveforms, + extend_durations=extend_durations, + extract_from_ramp=extract_from_ramp, + additional_fvs=additional_fvs, + file_list=file_list, + sweep_qc_record=sweep_qc_record_df, + manual_fail_sweeps=manual_fail_sweep_dict) + + logging.info("Number of specimens to process: {:d}".format(len(specimen_ids))) + if run_parallel: + with ProcessPoolExecutor(max_workers=os.cpu_count() - 1) as executor: + results = executor.map(get_data_partial, specimen_ids) + else: + results = map(get_data_partial, specimen_ids) + + used_ids, results, error_set = su.filter_results(specimen_ids, results) + + logging.info("Finished with {:d} processed specimens".format(len(used_ids))) + + results_dict = su.organize_results(used_ids, results, skip_keys=["long_squares_data_info"]) + + su.save_results_to_h5(used_ids, results_dict, output_dir, output_code) + + su.save_errors_to_json(error_set, output_dir, output_code) + + amp_file_name = os.path.join(output_dir, "fv_amplitudes_with_data_{}.json".format(output_code)) + + amp_info = {spec_id: r["long_squares_data_info"] for spec_id, r in zip(used_ids, results)} + + with open(amp_file_name, "w") as f: + json.dump(amp_info, f) + + logging.info("Finished saving") + + +def main(args): + ids = np.genfromtxt(args["specimen_id_file"], dtype=int).tolist() + + nwb_path_file = args["nwb_path_file"] + if nwb_path_file is None: + file_list = lq.get_nwb_file_paths_for_specimen_ids(ids) + else: + with open(nwb_path_file, "r") as f: + file_list = json.load(f) + + # convert string keys to integers + file_list = {int(k): v for k, v in file_list.items()} + + sweep_qc_record_file = args["sweep_qc_record_file"] + if sweep_qc_record_file is None: + sweep_qc_record = lq.get_sweep_states_and_tags_for_specimens(ids) + sweep_qc_record_df = pd.DataFrame(sweep_qc_record) + sweep_qc_record_df["tag_name"] = sweep_qc_record_df["tag_name"].fillna("None") + else: + sweep_qc_record_df = pd.read_csv(sweep_qc_record_file, index_col=0) + sweep_qc_record_df["tag_name"] = sweep_qc_record_df["tag_name"].fillna("None") + + manual_fail_sweep_file = args["manual_fail_sweep_file"] + if manual_fail_sweep_file is not None: + manual_fail_df = pd.read_csv(manual_fail_sweep_file) + manual_fail_sweep_dict = {} + for specimen_id in manual_fail_df.specimen_id.unique(): + sweeps_for_specimen = manual_fail_df.loc[manual_fail_df.specimen_id == specimen_id, "sweep_number"].tolist() + manual_fail_sweep_dict[specimen_id] = sweeps_for_specimen + else: + manual_fail_sweep_dict = None + + run_feature_vector_extraction( + specimen_ids=ids, + output_dir=args["output_dir"], + output_code=args["output_code"], + sweep_qc_option=args["sweep_qc_option"], + file_list=file_list, + sweep_qc_record_df=sweep_qc_record_df, + manual_fail_sweep_dict=manual_fail_sweep_dict, + extract_from_ramp=args["extract_from_ramp"], + amp_tolerance=args["amp_tolerance"], + additional_fvs=args["additional_fvs"], + extend_durations=args["extend_durations"], + ap_waveforms=args["ap_waveforms"], + needed_amplitudes=args["needed_amplitudes"], + run_parallel=args["run_parallel"], + ) + +if __name__ == "__main__": + module = ags.ArgSchemaParser(schema_type=CollectFeatureVectorParameters) + main(module.args) diff --git a/ipfx/defaults/stimulus_ontology.json b/ipfx/defaults/stimulus_ontology.json index 7f5145ae..d45ccb33 100644 --- a/ipfx/defaults/stimulus_ontology.json +++ b/ipfx/defaults/stimulus_ontology.json @@ -1252,7 +1252,8 @@ [ [ "code", - "X1PS_SubThresh" + "X1PS_SubThresh", + "Y1PS_SubThresh" ], [ "name", @@ -1282,7 +1283,8 @@ [ [ "code", - "X4PS_SupraThresh" + "X4PS_SupraThresh", + "Y4PS_SupraThresh_DA_1" ], [ "name", @@ -1292,7 +1294,8 @@ [ [ "code", - "X6SP_Rheo" + "X6SP_Rheo", + "Y6SP_Rheo" ], [ "name", @@ -1306,7 +1309,10 @@ ], [ "code", - "X7Ramp" + "X7Ramp", + "X7RAMP", + "X0NRMP", + "Y7RAMP" ] ], [ @@ -1372,7 +1378,8 @@ [ [ "code", - "X5SP_Search" + "X5SP_Search", + "Y5SP_Search" ], [ "name", @@ -2058,5 +2065,20 @@ "name", "Chirp A Threshold" ] + ], + [ + [ + "code", + "C2CHIRPB", + "C2CHIRPB180503" + ], + [ + "core", + "Core 2" + ], + [ + "name", + "Chirp B Threshold" + ] ] -] \ No newline at end of file +] diff --git a/ipfx/epochs.py b/ipfx/epochs.py index 98bfa8e3..3eaffee2 100644 --- a/ipfx/epochs.py +++ b/ipfx/epochs.py @@ -41,12 +41,12 @@ def get_last_stability_epoch(idx1, hz): return idx0, idx1 -def get_first_noise_epoch(idx, hz): +def get_noise_epoch_from_start(idx, hz): return idx, idx + int(NOISE_EPOCH * hz) -def get_last_noise_epoch(idx1, hz): +def get_noise_epoch_from_end(idx1, hz): return idx1-int(NOISE_EPOCH * hz), idx1 @@ -113,7 +113,7 @@ def get_stim_epoch(i, test_pulse=True): if test_pulse: di_idx = di_idx[2:] # drop the first up/down (test pulse) if present - if len(di_idx) < 2: # if no stimulus is found + if len(di_idx) == 0: # if no stimulus is found return None start_idx = di_idx[0] + 1 # shift by one to compensate for diff() diff --git a/ipfx/feature_extractor.py b/ipfx/feature_extractor.py index 08d4061c..e8398ef3 100644 --- a/ipfx/feature_extractor.py +++ b/ipfx/feature_extractor.py @@ -84,14 +84,30 @@ def __init__(self, start=None, end=None, filter=10., self.thresh_frac = thresh_frac self.reject_at_stim_start_interval = reject_at_stim_start_interval - def process(self, t, v, i): + def process(self, t, v, i, sweep_index=None): dvdt = tsu.calculate_dvdt(v, t, self.filter) + if type(self.start) is list: + if sweep_index is not None: + start = self.start[sweep_index] + else: + start = self.start[0] + else: + start = self.start + + if type(self.end) is list: + if sweep_index is not None: + end = self.end[sweep_index] + else: + end = self.end[0] + else: + end = self.end + # Basic features of spikes - putative_spikes = spkd.detect_putative_spikes(v, t, self.start, self.end, + putative_spikes = spkd.detect_putative_spikes(v, t, start, end, dv_cutoff=self.dv_cutoff, dvdt=dvdt) - peaks = spkd.find_peak_indexes(v, t, putative_spikes, self.end) + peaks = spkd.find_peak_indexes(v, t, putative_spikes, end) putative_spikes, peaks = spkd.filter_putative_spikes(v, t, putative_spikes, peaks, self.min_height, self.min_peak, dvdt=dvdt) @@ -105,7 +121,7 @@ def process(self, t, v, i): dvdt=dvdt) thresholds, peaks, upstrokes, clipped = spkd.check_thresholds_and_peaks(v, t, thresholds, peaks, - upstrokes, self.start, self.end, self.max_interval, + upstrokes, start, end, self.max_interval, dvdt=dvdt, reject_at_stim_start_interval=self.reject_at_stim_start_interval) if not thresholds.size: @@ -114,9 +130,9 @@ 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) + troughs = spkd.find_trough_indexes(v, t, thresholds, peaks, clipped, end) 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, + trough_details, clipped = spkf.analyze_trough_details(v, t, thresholds, peaks, clipped, end, dvdt=dvdt) widths = spkf.find_widths(v, t, thresholds, peaks, trough_details[1], clipped) @@ -261,34 +277,52 @@ def __init__(self, start, end, self.sag_baseline_interval = sag_baseline_interval self.peak_width = peak_width - def process(self, t, v, i, spikes_df, extra_features=None, exclude_clipped=False): - features = strf.basic_spike_train_features(t, spikes_df, self.start, self.end, exclude_clipped=exclude_clipped) + def process(self, t, v, i, spikes_df, extra_features=None, exclude_clipped=False, sweep_index=None): if self.start is None: self.start = 0.0 + if type(self.start) is list: + if sweep_index is not None: + start = self.start[sweep_index] + else: + start = self.start[0] + else: + start = self.start + + if type(self.end) is list: + if sweep_index is not None: + end = self.end[sweep_index] + else: + end = self.end[0] + else: + end = self.end + + features = strf.basic_spike_train_features(t, spikes_df, start, end, exclude_clipped=exclude_clipped) + if extra_features is None: extra_features = [] if 'peak_deflect' in extra_features: - features['peak_deflect'] = subf.voltage_deflection(t, v, i, self.start, self.end, self.deflect_type) + features['peak_deflect'] = subf.voltage_deflection(t, v, i, start, end, self.deflect_type) if 'stim_amp' in extra_features: - features['stim_amp'] = self.stim_amp_fn(t, i, self.start) if self.stim_amp_fn else None + features['stim_amp'] = self.stim_amp_fn(t, i, start) if self.stim_amp_fn else None if 'v_baseline' in extra_features: - features['v_baseline'] = subf.baseline_voltage(t, v, self.start, self.baseline_interval, self.filter_frequency) + features['v_baseline'] = subf.baseline_voltage( + t, v, start, self.baseline_interval, self.filter_frequency) if 'sag' in extra_features: - features['sag'] = subf.sag(t, v, i, self.start, self.end, self.peak_width, self.sag_baseline_interval) + features['sag'] = subf.sag(t, v, i, start, end, self.peak_width, self.sag_baseline_interval) if features["avg_rate"] > 0: if 'pause' in extra_features: - features['pause'] = strf.pause(t, spikes_df, self.start, self.end, self.pause_cost_weight) + features['pause'] = strf.pause(t, spikes_df, start, end, self.pause_cost_weight) if 'burst' in extra_features: features['burst'] = strf.burst(t, spikes_df, self.burst_tol, self.pause_cost) if 'delay' in extra_features: - features['delay'] = strf.delay(t, v, spikes_df, self.start, self.end) + features['delay'] = strf.delay(t, v, spikes_df, start, end) return features diff --git a/ipfx/feature_vectors.py b/ipfx/feature_vectors.py index 126a4efe..d870ca53 100644 --- a/ipfx/feature_vectors.py +++ b/ipfx/feature_vectors.py @@ -1,10 +1,12 @@ import numpy as np import logging from scipy import stats +from scipy.signal import savgol_filter from . import data_set_features as dsf from . import stimulus_protocol_analysis as spa from . import time_series_utils as tsu from . import error as er +from ipfx.script_utils import StimulusTiming def identify_subthreshold_hyperpol_with_amplitudes(features, sweeps): @@ -93,9 +95,9 @@ def identify_subthreshold_depol_with_amplitudes(features, sweeps): return amp_sweep_dict, deflect_dict -def step_subthreshold(amp_sweep_dict, target_amps, start, end, - extend_duration=0.2, subsample_interval=0.01, - amp_tolerance=0.): +def step_subthreshold(amp_sweep_dict, target_amps, stim_timing_dict, + extend_duration_before=0.2, extend_duration_after=0.2, subsample_interval=0.01, + amp_tolerance=0., remove_transients=True): """ Subsample set of subthreshold step responses including regions before and after step Parameters @@ -124,14 +126,26 @@ def step_subthreshold(amp_sweep_dict, target_amps, start, end, subsampled_dict = {} for amp in amp_sweep_dict: swp = amp_sweep_dict[amp] - start_index = tsu.find_time_index(swp.t, start - extend_duration) + stim_timing = stim_timing_dict[swp.sweep_number] + + if remove_transients: + v = swp.v.copy() + start_index = tsu.find_time_index(swp.t, stim_timing.start) + end_index = tsu.find_time_index(swp.t, stim_timing.end) + v_clean = _remove_transients(v[start_index:end_index], swp.t[start_index:end_index]) + v[start_index:end_index] = v_clean + else: + v = swp.v + + start_index = tsu.find_time_index(swp.t, stim_timing.start - extend_duration_before) + end_index = tsu.find_time_index(swp.t, stim_timing.end + extend_duration_after) delta_t = swp.t[1] - swp.t[0] subsample_width = int(np.round(subsample_interval / delta_t)) - end_index = tsu.find_time_index(swp.t, end + extend_duration) - subsampled_v = _subsample_average(swp.v[start_index:end_index], subsample_width) + subsampled_v = _subsample_average(v[start_index:end_index], subsample_width) subsampled_dict[amp] = subsampled_v - extend_length = int(np.round(extend_duration / subsample_interval)) + extend_length_before = int(np.round(extend_duration_before / subsample_interval)) + extend_length_after = int(np.round(extend_duration_after / subsample_interval)) available_amps = np.array(list(subsampled_dict.keys())) output_list = [] for amp in target_amps: @@ -159,25 +173,68 @@ def step_subthreshold(amp_sweep_dict, target_amps, start, end, logging.debug("interpolating for amp {} with lower {} and upper {}".format(amp, lower_amp, upper_amp)) avg = (subsampled_dict[lower_amp] + subsampled_dict[upper_amp]) / 2. scale = amp / ((lower_amp + upper_amp) / 2.) - base_v = avg[:extend_length].mean() - avg[extend_length:-extend_length] = (avg[extend_length:-extend_length] - base_v) * scale + base_v + base_v = avg[:extend_length_before].mean() + avg[extend_length_before:-extend_length_after] = (avg[extend_length_before:-extend_length_after] - base_v) * scale + base_v elif lower_amp != 0: logging.debug("interpolating for amp {} from lower {}".format(amp, lower_amp)) avg = subsampled_dict[lower_amp].copy() scale = amp / lower_amp - base_v = avg[:extend_length].mean() - avg[extend_length:] = (avg[extend_length:] - base_v) * scale + base_v + base_v = avg[:extend_length_before].mean() + avg[extend_length_before:] = (avg[extend_length_before:] - base_v) * scale + base_v elif upper_amp != 0: logging.debug("interpolating for amp {} from upper {}".format(amp, upper_amp)) avg = subsampled_dict[upper_amp].copy() scale = amp / upper_amp - base_v = avg[:extend_length].mean() - avg[extend_length:] = (avg[extend_length:] - base_v) * scale + base_v + base_v = avg[:extend_length_before].mean() + avg[extend_length_before:] = (avg[extend_length_before:] - base_v) * scale + base_v output_list.append(avg) return np.hstack(output_list) +def _remove_transients(v, t, dvdt_thresh=2.0, window_width=400): + dvdt = savgol_filter(v, 50, 2, deriv=1, delta=1e3 * (t[1] - t[0])) # mV/ms, smoothed + v_clean = v.copy() + + # use an adaptive threshold to avoid treating the start of the step as a transient + t_envelope = 1e3 * (t - t[0]) + dvdt_thresh_envelope = 13 * np.exp(-t_envelope / 5) + dvdt_thresh + + while np.any(np.abs(dvdt) > dvdt_thresh_envelope): + exceed_inds = np.flatnonzero(np.abs(dvdt) > dvdt_thresh_envelope) + exceed_peak_ind = np.argmax(np.abs(dvdt)[exceed_inds]) + dvdt_peak_ind = exceed_inds[exceed_peak_ind] + peak_ind = np.nanargmax(np.abs(v[dvdt_peak_ind - window_width:dvdt_peak_ind + window_width])) + peak_ind += dvdt_peak_ind - window_width + + search_start = min(peak_ind, dvdt_peak_ind) + transient_start_index = np.flatnonzero(np.abs(dvdt[search_start::-1]) < dvdt_thresh / 10)[0] + transient_start_index = search_start - transient_start_index + + transient_base_avg = np.mean(v[transient_start_index - window_width:transient_start_index]) + transient_base_range = 3 * np.std(v[transient_start_index - window_width:transient_start_index]) + + search_start = max(peak_ind, dvdt_peak_ind) + + baseline_return = np.flatnonzero(np.abs(v[search_start:] - transient_base_avg) < transient_base_range) + if len(baseline_return) > 0: + transient_end_index = baseline_return[0] + transient_end_index += search_start + else: + transient_end_index = len(t) - 1 + + # blank out the transient + v_clean[transient_start_index:transient_end_index + 1] = np.nan + dvdt[transient_start_index:transient_end_index + 1] = np.nan + + nan_mask = np.isnan(v_clean) + nan_ind = np.nonzero(nan_mask) + not_nan_ind = np.nonzero(~nan_mask) + + v_clean[nan_ind] = np.interp(t[nan_ind], t[not_nan_ind], v_clean[not_nan_ind]) + return v_clean + + def _subsample_average(x, width): """Downsamples x by averaging `width` points""" @@ -185,8 +242,9 @@ def _subsample_average(x, width): return avg -def subthresh_norm(amp_sweep_dict, deflect_dict, start, end, target_amp=-101., - extend_duration=0.2, subsample_interval=0.01): +def subthresh_norm(amp_sweep_dict, deflect_dict, stim_timing_dict, target_amp=-101., + extend_duration_before=0.2, extend_duration_after=0.2, + subsample_interval=0.01, remove_transients=True): """ Subthreshold step response closest to target amplitude normalized to baseline and peak deflection Parameters @@ -219,18 +277,81 @@ def subthresh_norm(amp_sweep_dict, deflect_dict, start, end, target_amp=-101., base, deflect_v = deflect_dict[matching_amp] delta = base - deflect_v - start_index = tsu.find_time_index(swp.t, start - extend_duration) + stim_timing = stim_timing_dict[swp.sweep_number] + + if remove_transients: + v = swp.v.copy() + start_index = tsu.find_time_index(swp.t, stim_timing.start) + end_index = tsu.find_time_index(swp.t, stim_timing.end) + v_clean = _remove_transients(v[start_index:end_index], swp.t[start_index:end_index]) + v[start_index:end_index] = v_clean + else: + v = swp.v + + + start_index = tsu.find_time_index(swp.t, stim_timing.start - extend_duration_before) + end_index = tsu.find_time_index(swp.t, stim_timing.end + extend_duration_after) delta_t = swp.t[1] - swp.t[0] subsample_width = int(np.round(subsample_interval / delta_t)) - end_index = tsu.find_time_index(swp.t, end + extend_duration) - subsampled_v = _subsample_average(swp.v[start_index:end_index], subsample_width) + subsampled_v = _subsample_average(v[start_index:end_index], subsample_width) subsampled_v -= base subsampled_v /= delta return subsampled_v -def subthresh_depol_norm(amp_sweep_dict, deflect_dict, start, end, +def subthresh_rebound(amp_sweep_dict, stim_timing_dict, dur=0.2, target_amp=-101., + subsample_interval=0.01, psth_bin_width=20): + """ Subthreshold step rebound response closest to target amplitude + + Parameters + ---------- + amp_sweep_dict: dict + Amplitude-sweep pairs + start: float + start stimulus interval (seconds) + dur: float + duration of rebound interval (seconds) + target_amp: float (optional, default=-101) + Search target for amplitude (pA) + extend_duration_[before, after]: float (optional, default 0.05) + Durations to extend sweep on either side of interval (seconds) + subsample_interval: float (optional, default 0.01) + Size of subsampled bins (seconds) + + Returns + ------- + subsampled_v: array + Subsampled, normalized voltage trace + """ + available_amps = np.array(list(amp_sweep_dict.keys())) + + sweep_ind = np.argmin(np.abs(available_amps - target_amp)) + matching_amp = available_amps[sweep_ind] + swp = amp_sweep_dict[matching_amp] + + swp_start = stim_timing_dict[swp.sweep_number].start + + # Find spikes in rebound interval + spx, spfx = dsf.extractors_for_sweeps( + SweepSet(sweeps=[swp]), + start=swp_start, + end=swp_start + dur, + min_peak=-25 + ) + spike_data = spx.process(swp.t, swp.v, swp.i, sweep_index=0) + + if spike_data.shape[0] == 0: + # No spikes + spike_data = pd.DataFrame(columns=["threshold_t"]) + + rebound_psth = psth_vector( + [spike_data], [StimulusTiming(start=swp_start, end=swp_start + dur, dur=dur)], width=psth_bin_width, duration=dur) + + return rebound_psth + + +def subthresh_depol_norm(amp_sweep_dict, deflect_dict, stim_timing_dict, extend_duration=0.2, subsample_interval=0.01, steady_state_interval=0.1): """ Largest positive-going subthreshold step response that does not evoke spikes, normalized to baseline and steady-state at end of step @@ -258,33 +379,39 @@ def subthresh_depol_norm(amp_sweep_dict, deflect_dict, start, end, Subsampled, normalized voltage trace """ - if (end - start) < steady_state_interval: - raise ValueError("steady state interval cannot exceed stimulus interval") if len(amp_sweep_dict) == 0: logging.debug("No subthreshold depolarizing sweeps found - returning all-nan response") + # Weird because not really using this sweep, but will just take the first sweep + # in the dict + stim_timing = list(stim_timing_dict.values())[0] + # create all-nan response of appropriate length - total_interval = extend_duration * 2 + (end - start) + total_interval = extend_duration * 2 + stim_timing.dur length = int(total_interval / subsample_interval) return np.ones(length) * np.nan available_amps = list(amp_sweep_dict.keys()) max_amp = np.max(available_amps) swp = amp_sweep_dict[max_amp] + stim_timing = stim_timing_dict[swp.sweep_number] + + if stim_timing.dur < steady_state_interval: + raise ValueError("steady state interval cannot exceed stimulus interval") base, _ = deflect_dict[max_amp] - interval_start_index = tsu.find_time_index(swp.t, end - steady_state_interval) - interval_end_index = tsu.find_time_index(swp.t, end) + interval_start_index = tsu.find_time_index(swp.t, stim_timing.end - steady_state_interval) + interval_end_index = tsu.find_time_index(swp.t, stim_timing.end) steady_state_v = swp.v[interval_start_index:interval_end_index].mean() delta = steady_state_v - base - start_index = tsu.find_time_index(swp.t, start - extend_duration) + start_index = tsu.find_time_index(swp.t, stim_timing.start - extend_duration) + end_index = tsu.find_time_index(swp.t, stim_timing.end + extend_duration) delta_t = swp.t[1] - swp.t[0] subsample_width = int(np.round(subsample_interval / delta_t)) - end_index = tsu.find_time_index(swp.t, end + extend_duration) subsampled_v = _subsample_average(swp.v[start_index:end_index], subsample_width) subsampled_v -= base subsampled_v /= delta @@ -292,7 +419,8 @@ def subthresh_depol_norm(amp_sweep_dict, deflect_dict, start, end, return subsampled_v -def identify_sweep_for_isi_shape(sweeps, features, duration, min_spike=5): +def identify_sweep_for_isi_shape(sweeps, features, stim_timing_dict, min_spike=5, + exclude_sweep_numbers=[]): """ Find lowest-amplitude spiking sweep that has at least min_spike or else sweep with most spikes @@ -315,9 +443,16 @@ def identify_sweep_for_isi_shape(sweeps, features, duration, min_spike=5): Spike info for selected sweep """ sweep_table = features["sweeps"] - mask_supra = (sweep_table["avg_rate"].values > 0) & (sweep_table["stim_amp"] > 0) - supra_table = sweep_table.loc[mask_supra, :] + sweep_numbers = np.array([s.sweep_number for s in sweeps.sweeps], dtype=int) + mask_exclude = np.array([s.sweep_number in exclude_sweep_numbers for s in sweeps.sweeps]) + + include_table = sweep_table.loc[~mask_exclude, :] + mask_supra = (include_table["avg_rate"].values > 0) & (include_table["stim_amp"] > 0) + supra_table = include_table.loc[mask_supra, :] + duration = np.array([stim_timing_dict[sn].dur for sn in sweep_numbers])[~mask_exclude][mask_supra] + amps = np.rint(supra_table["stim_amp"].values) + n_spikes = supra_table["avg_rate"].values * duration # Pick out the sweep to get the ISI shape @@ -337,13 +472,13 @@ def identify_sweep_for_isi_shape(sweeps, features, duration, min_spike=5): only_one_spike = True selection_index = np.argmin(amps) - selected_sweep = np.array(sweeps.sweeps)[mask_supra][selection_index] + selected_sweep = np.array(sweeps.sweeps)[~mask_exclude][mask_supra][selection_index] info_index = supra_table.index.tolist()[selection_index] selected_spike_info = features["spikes_set"][info_index] return selected_sweep, selected_spike_info -def isi_shape(sweep, spike_info, end, n_points=100, steady_state_interval=0.1, +def isi_shape(sweep, spike_info, stim_timing_dict, n_points=100, steady_state_interval=0.1, single_return_tolerance=1., single_max_duration=0.1): """ Average interspike voltage trajectory with normalized duration, aligned to threshold @@ -373,8 +508,13 @@ def isi_shape(sweep, spike_info, end, n_points=100, steady_state_interval=0.1, Averaged, threshold-aligned, duration-normalized voltage trace """ + spike_info = spike_info.copy() + + # only consider non-clipped spikes + spike_info = spike_info.loc[~spike_info["clipped"], :] n_spikes = spike_info.shape[0] + stim_timing = stim_timing_dict[sweep.sweep_number] if n_spikes > 1: threshold_indexes = spike_info["threshold_index"].values threshold_voltages = spike_info["threshold_v"].values @@ -394,11 +534,11 @@ def isi_shape(sweep, spike_info, end, n_points=100, steady_state_interval=0.1, threshold_v = spike_info["threshold_v"][0] fast_trough_index = spike_info["fast_trough_index"].astype(int)[0] fast_trough_t = spike_info["fast_trough_t"][0] - stim_end_index = tsu.find_time_index(sweep.t, end) - if fast_trough_t < end - steady_state_interval: + stim_end_index = tsu.find_time_index(sweep.t, stim_timing.end) + if fast_trough_t < stim_timing.end - steady_state_interval: max_end_index = tsu.find_time_index(sweep.t, sweep.t[fast_trough_index] + single_max_duration) - std_start_index = tsu.find_time_index(sweep.t, end - steady_state_interval) + std_start_index = tsu.find_time_index(sweep.t, stim_timing.end - steady_state_interval) steady_state_v = sweep.v[std_start_index:stim_end_index].mean() above_ss_ind = np.flatnonzero(sweep.v[fast_trough_index:] >= steady_state_v - single_return_tolerance) @@ -580,8 +720,8 @@ def first_ap_waveform(sweep, spikes, length_in_points): return sweep.v[start_index:end_index] -def identify_suprathreshold_spike_info(features, target_amplitudes, - shift=None, amp_tolerance=0): +def identify_suprathreshold_spike_info(features, target_amplitudes, sweep_numbers=None, + shift=None, amp_tolerance=0, needed_amplitudes=None): """ Find spike information for sweeps matching desired amplitudes relative to rheobase Parameters @@ -606,12 +746,17 @@ def identify_suprathreshold_spike_info(features, target_amplitudes, spike_data = features["spikes_set"] sweeps_to_use = _identify_suprathreshold_indices( - features, target_amplitudes, shift, amp_tolerance) - return [spike_data[ind] if ind is not None else None for ind in sweeps_to_use] + features, target_amplitudes, shift, amp_tolerance, needed_amplitudes=needed_amplitudes) + if sweep_numbers is not None: + used_sweep_numbers = [sweep_numbers[ind] if ind is not None else None for ind in sweeps_to_use] + else: + used_sweep_numbers = None + return [spike_data[ind] if ind is not None else None for ind in sweeps_to_use], used_sweep_numbers + def identify_suprathreshold_sweeps(sweeps, features, target_amplitudes, - shift=None, amp_tolerance=0): + shift=None, amp_tolerance=0, needed_amplitudes=None): """ Find spike information for sweeps matching desired amplitudes relative to rheobase Parameters @@ -628,7 +773,8 @@ def identify_suprathreshold_sweeps(sweeps, features, target_amplitudes, A value of None means that no shift is attempted. amp_tolerance: float (optional, default 0) Tolerance for matching amplitude (pA) - + needed_amplitudes: list-like (optional, default None) + Subset of `target_amplitudes` of which at least two are required Returns ------- sweeps: list @@ -637,12 +783,12 @@ def identify_suprathreshold_sweeps(sweeps, features, target_amplitudes, """ sweeps_to_use = _identify_suprathreshold_indices( - features, target_amplitudes, shift, amp_tolerance) + features, target_amplitudes, shift, amp_tolerance, needed_amplitudes=needed_amplitudes) return [sweeps.sweeps[ind] if ind is not None else None for ind in sweeps_to_use] def _identify_suprathreshold_indices(features, target_amplitudes, - shift=None, amp_tolerance=0): + shift=None, amp_tolerance=0, needed_amplitudes=None): """ Find indices for sweeps matching desired amplitudes relative to rheobase Parameters @@ -680,21 +826,68 @@ def _identify_suprathreshold_indices(features, target_amplitudes, sweeps_to_use = _spiking_sweeps_at_levels(amps, sweep_indexes, target_amplitudes, amp_tolerance) + orig_rheo_ind = sweeps_to_use[0] n_matches = np.sum([s is not None for s in sweeps_to_use]) + if needed_amplitudes is not None: + n_needed_matches = np.sum([s is not None for s, a in zip(sweeps_to_use, target_amplitudes) if a in needed_amplitudes]) + else: + n_needed_matches = n_matches - if len(target_amplitudes) > 1 and n_matches <= 1 and shift is not None: + if len(target_amplitudes) > 1 and n_needed_matches <= 1 and shift is not None: logging.debug("Found only one spiking sweep that matches expected amplitude levels; attempting to shift by {} pA".format(shift)) sweeps_to_use = _spiking_sweeps_at_levels(amps - shift, sweep_indexes, target_amplitudes, amp_tolerance) n_matches = np.sum([s is not None for s in sweeps_to_use]) + if needed_amplitudes is not None: + n_needed_matches = np.sum([s is not None for s, a in zip(sweeps_to_use, target_amplitudes) if a in needed_amplitudes]) + else: + n_needed_matches = n_matches + + # Compensate for earlier issue where shifting could lose the rheobase sweep + if sweeps_to_use[0] is None: + sweeps_to_use[0] = orig_rheo_ind + n_needed_matches += 1 + + if len(target_amplitudes) > 1 and n_needed_matches <= 1: + # last ditch - see if number of spikes on rheo and next highest available are same (or +20 is less) and try those + alt_target_amplitudes = sorted(np.unique(amps)) + alt_sweeps_to_use = _spiking_sweeps_at_levels(amps, sweep_indexes, + alt_target_amplitudes, amp_tolerance=0) + alt_rheo_ind = alt_sweeps_to_use[0] + logging.debug(alt_target_amplitudes) + logging.debug(sweeps_to_use) + start_range = 1 + keep_going = True + found_match = False + for i, try_ind in enumerate(alt_sweeps_to_use[1:]): + logging.debug(f"Trying sweep at index {try_ind}") + logging.debug(f"Sweep at {alt_target_amplitudes[i + 1]} had {np.round(sweep_table.at[try_ind, 'avg_rate'])} spikes/s vs rheo {np.round(sweep_table.at[alt_rheo_ind, 'avg_rate'])}") + if np.round(sweep_table.at[alt_rheo_ind, "avg_rate"]) >= np.round(sweep_table.at[try_ind, "avg_rate"]): + alt_shift = np.round(sweep_table.at[try_ind, "stim_amp"] - sweep_table.at[alt_rheo_ind, "stim_amp"]) + logging.debug(f"Trying shift of {alt_shift} pA") + sweeps_to_use = _spiking_sweeps_at_levels(amps - alt_shift, sweep_indexes, + target_amplitudes, amp_tolerance) + n_matches = np.sum([s is not None for s in sweeps_to_use]) + if needed_amplitudes is not None: + n_needed_matches = np.sum([s is not None for s, a in zip(sweeps_to_use, target_amplitudes) if a in needed_amplitudes]) + else: + n_needed_matches = n_matches + if len(target_amplitudes) > 1 and n_needed_matches <= 1: + logging.debug("No match yet") + else: + found_match = True + logging.info("Had to shift by {} pA to get more than one matching sweep".format(alt_shift)) + break + else: + break - if len(target_amplitudes) > 1 and n_matches <= 1: - raise er.FeatureError("Could not find at least two spiking sweeps matching requested amplitude levels") + if not found_match: + raise er.FeatureError(f"Could not find at least two spiking sweeps matching requested amplitude levels (available: {amps})") return sweeps_to_use -def psth_vector(spike_info_list, start, end, width=50): +def psth_vector(spike_info_list, stim_timing_list, width=50, duration=1.0): """ Create binned "PSTH"-like feature vector based on spike times, concatenated across sweeps @@ -702,9 +895,9 @@ def psth_vector(spike_info_list, start, end, width=50): ---------- spike_info_list: list Spike info DataFrames for each sweep - start: float + start: float or list Start of stimulus interval (seconds) - end: float + end: float or list End of stimulus interval (seconds) width: float (optional, default 50) Bin width in ms @@ -716,32 +909,40 @@ def psth_vector(spike_info_list, start, end, width=50): """ vector_list = [] - for si in spike_info_list: + if spike_info_list[0] is None: + logging.warning("Rheobase sweep appears to be missing") + + one_ms = 0.001 + n_bins = int(duration / one_ms) // width + for si, stim_timing in zip(spike_info_list, stim_timing_list): if si is None: vector_list.append(None) continue thresh_t = si["threshold_t"] spike_count = np.ones_like(thresh_t) - one_ms = 0.001 - - # round to nearest ms to deal with float approximations - duration = np.round(end, decimals=3) - np.round(start, decimals=3) + # only use actual duration to check against requested duration + if np.abs(stim_timing.dur - duration) > one_ms: + logging.warning(f"Actual duration ({stim_timing.dur}) does not match duration specified for analysis ({duration})") - n_bins = int(duration / one_ms) // width - bin_edges = np.linspace(start, end, n_bins + 1) # includes right edge, so adding one to desired bin number + bin_edges = np.linspace(stim_timing.start, stim_timing.start + duration, n_bins + 1) # includes right edge, so adding one to desired bin number bin_width = bin_edges[1] - bin_edges[0] - output = stats.binned_statistic(thresh_t, - spike_count, - statistic='sum', - bins=bin_edges)[0] - output[np.isnan(output)] = 0 - output /= bin_width # convert to spikes/s + + if len(thresh_t) == 0: + # No spikes - return all 0 vector_list + output = np.zeros(n_bins) + else: + output = stats.binned_statistic(thresh_t, + spike_count, + statistic='sum', + bins=bin_edges)[0] + output[np.isnan(output)] = 0 + output /= bin_width # convert to spikes/s vector_list.append(output) output_vector = _combine_and_interpolate(vector_list) return output_vector -def inst_freq_vector(spike_info_list, start, end, width=20): +def inst_freq_vector(spike_info_list, stim_timing_list, width=20, gap_factor=4, duration=1.0): """ Create binned instantaneous frequency feature vector, concatenated across sweeps @@ -749,13 +950,16 @@ def inst_freq_vector(spike_info_list, start, end, width=20): ---------- spike_info_list: list Spike info DataFrames for each sweep - start: float + start: float or list Start of stimulus interval (seconds) - end: float + end: float or list End of stimulus interval (seconds) width: float (optional, default 20) Bin width in ms - + gap_factor: int, default 4 + Factor to multiply average ISI by to determine if a gap should + be interpolated over or set to zero spikes/s (i.e., if the gap exceeds + avg_isi * gap_factor, it will be set to zero). Returns ------- @@ -763,29 +967,50 @@ def inst_freq_vector(spike_info_list, start, end, width=20): Concatenated vector of binned instantaneous firing rates (spikes/s) """ + if spike_info_list[0] is None: + logging.warning("Rheobase sweep appears to be missing") + vector_list = [] - for si in spike_info_list: + one_ms = 0.001 + n_bins = int(duration / one_ms) // width + for si, stim_timing in zip(spike_info_list, stim_timing_list): if si is None: vector_list.append(None) continue thresh_t = si["threshold_t"].values - inst_freq, inst_freq_times = _inst_freq_feature(thresh_t, start, end) - - one_ms = 0.001 + inst_freq, inst_freq_times = _inst_freq_feature(thresh_t, stim_timing.start, stim_timing.end) - # round to nearest ms to deal with float approximations - duration = np.round(end, decimals=3) - np.round(start, decimals=3) + # only use actual duration to check against requested duration + if np.abs(stim_timing.dur - duration) > one_ms: + logging.warning(f"Actual duration ({stim_timing.dur}) does not match duration specified for analysis ({duration})") - n_bins = int(duration / one_ms) // width - bin_edges = np.linspace(start, end, n_bins + 1) # includes right edge, so adding one to desired bin number + bin_edges = np.linspace(stim_timing.start, stim_timing.start + duration, n_bins + 1) # includes right edge, so adding one to desired bin number bin_width = bin_edges[1] - bin_edges[0] output = stats.binned_statistic(inst_freq_times, inst_freq, bins=bin_edges)[0] - nan_ind = np.isnan(output) + + # Check for long gaps without spikes + nan_ind = np.flatnonzero(np.isnan(output)) + consecutive_sections = np.split(nan_ind, np.where(np.diff(nan_ind) != 1)[0] + 1) + + avg_isi = 1 / np.mean(inst_freq) + for sec in consecutive_sections: + if len(sec) == 0: + continue + gap_length = (sec[-1] - sec[0] + 1) * bin_width + if gap_length > gap_factor * avg_isi: + # Since gap is long, set beginning and end of gap to 0 spikes/s + output[sec[0]] = 0 + output[sec[-1]] = 0 + + # Get mask for interpolation + nan_mask = np.isnan(output) + + # Interpolate missing values x = np.arange(len(output)) - output[nan_ind] = np.interp(x[nan_ind], x[~nan_ind], output[~nan_ind]) + output[nan_mask] = np.interp(x[nan_mask], x[~nan_mask], output[~nan_mask]) vector_list.append(output) output_vector = _combine_and_interpolate(vector_list) @@ -793,7 +1018,8 @@ def inst_freq_vector(spike_info_list, start, end, width=20): return output_vector -def spike_feature_vector(feature, spike_info_list, start, end, width=20): + +def spike_feature_vector(feature, spike_info_list, stim_timing_list, width=20, duration=1.0): """ Create binned feature vector for specified features, concatenated across sweeps @@ -803,9 +1029,9 @@ def spike_feature_vector(feature, spike_info_list, start, end, width=20): Name of feature found in members of spike_info_list spike_info_list: list Spike info DataFrames for each sweep - start: float + start: float or list Start of stimulus interval (seconds) - end: float + end: float or list End of stimulus interval (seconds) width: float (optional, default 20) Bin width in ms @@ -816,8 +1042,14 @@ def spike_feature_vector(feature, spike_info_list, start, end, width=20): Concatenated vector of binned spike features """ + if spike_info_list[0] is None: + logging.warning("Rheobase sweep appears to be missing") + + one_ms = 0.001 + n_bins = int(duration / one_ms) // width + vector_list = [] - for si in spike_info_list: + for si, stim_timing in zip(spike_info_list, stim_timing_list): if si is None: vector_list.append(None) continue @@ -828,16 +1060,20 @@ def spike_feature_vector(feature, spike_info_list, start, end, width=20): else: feature_values = si[feature].values mask = ~si["clipped"].values - thresh_t = thresh_t[mask] - feature_values = feature_values[mask] - - one_ms = 0.001 + if np.sum(mask) == 0: + if np.all(np.isnan(feature_values)): + logging.warning(f"All spikes were clipped (n={len(thresh_t)}) and had NaN values for feature {feature}; dropping sweep for this feature") + vector_list.append(None) + continue + else: + thresh_t = thresh_t[mask] + feature_values = feature_values[mask] - # round to nearest ms to deal with float approximations - duration = np.round(end, decimals=3) - np.round(start, decimals=3) + # only use actual duration to check against requested duration + if np.abs(stim_timing.dur - duration) > one_ms: + logging.warning(f"Actual duration ({actual_duration}) does not match duration specified for analysis ({duration})") - n_bins = int(duration / one_ms) // width - bin_edges = np.linspace(start, end, n_bins + 1) # includes right edge, so adding one to desired bin number + bin_edges = np.linspace(stim_timing.start, stim_timing.start + duration, n_bins + 1) # includes right edge, so adding one to desired bin number bin_width = bin_edges[1] - bin_edges[0] output = stats.binned_statistic(thresh_t, @@ -852,15 +1088,22 @@ def spike_feature_vector(feature, spike_info_list, start, end, width=20): return output_vector -def _spiking_sweeps_at_levels(amps, sweep_indexes, target_amplitudes, - amp_tolerance): + +def _spiking_sweeps_at_levels( + amps, sweep_indexes, target_amplitudes, amp_tolerance): """Search for sweep indexes that match target amplitudes""" sweeps_to_use = [] for target_amp in target_amplitudes: + # find exact 0 relative amplitude sweep for rheobase (which must exist); otherwise use amp_tolerance + if target_amp == 0: + used_amp_tolerance = 0 + else: + used_amp_tolerance = amp_tolerance + found_match = False for amp, swp_ind in zip(amps, sweep_indexes): - if (np.abs(amp - target_amp) <= amp_tolerance): + if (np.abs(amp - target_amp) <= used_amp_tolerance) and swp_ind not in sweeps_to_use: found_match = True sweeps_to_use.append(swp_ind) logging.debug("Using amplitude {} for target {}".format(amp, target_amp)) @@ -871,6 +1114,7 @@ def _spiking_sweeps_at_levels(amps, sweep_indexes, target_amplitudes, return sweeps_to_use + def _consolidated_long_square_indexes(sweep_table): """Identify a single sweep for each stimulus amplitude if an amplitude is repeated @@ -931,21 +1175,16 @@ def _inst_freq_feature(threshold_t, start, end): This function attempts to estimate a semi-continuous instantanteous firing rate from a set of interspike intervals (ISIs) and spike times. It makes - several assumptions: + several assumptions/methodological decisions: + - It only estimates the firing frequency when at the times of spikes. - It assumes that the instantaneous firing rate at the start of the stimulus interval is the inverse of the latency to the first spike. - It estimates the firing rate at each spike as the average of the ISIs on each side of the spike - - If the time between the end of the interval and the last spike is less - than the last true interspike interval, it sets the instantaneous rate of - that last spike and of the end of the interval to the inverse of the last ISI. - Therefore, the instantaneous rate would not "jump" just because the - stimulus interval ends. - - However, if the time between the end of the interval and the last spike is - longer than the final ISI, it assumes there may have been a spike just - after the end of the interval. Therefore, it essentially returns an upper - bound on the estimated rate. - + - It does not consider the interval between the last spike and the end of + the stimulus as an interspike interval (because there is no spike at the end). + Consequently, it only uses the last actual ISI for the estimated rate + of the last spike of the train. Parameters ---------- @@ -968,22 +1207,16 @@ def _inst_freq_feature(threshold_t, start, end): inst_inv_rate = [] time_points = [] isis = [(threshold_t[0] - start)] + np.diff(threshold_t).tolist() - isis = isis + [max(isis[-1], end - threshold_t[-1])] - # Estimate at start of stimulus interval - inst_inv_rate.append(isis[0]) - time_points.append(start) + # Add an ISI for the end so the average for the last spike will be its + # prior ISI + isis = isis + [isis[-1]] # Estimate for each spike time for t, pre_isi, post_isi in zip(threshold_t, isis[:-1], isis[1:]): inst_inv_rate.append((pre_isi + post_isi) / 2) time_points.append(t) - # Estimate for end of stimulus interval - inst_inv_rate.append(isis[-1]) - time_points.append(end) - - inst_firing_rate = 1 / np.array(inst_inv_rate) time_points = np.array(time_points) return inst_firing_rate, time_points diff --git a/ipfx/lims_queries.py b/ipfx/lims_queries.py index 625c2bae..2dcd4fef 100755 --- a/ipfx/lims_queries.py +++ b/ipfx/lims_queries.py @@ -1,6 +1,6 @@ import os import logging -import pg8000 +import pg8000.dbapi from ipfx.string_utils import to_str @@ -39,7 +39,7 @@ def _connect(timeout=TIMEOUT): credentials = dict((k, os.environ.get(env_var, LIMS_DB_CREDENTIAL_DEFAULTS[env_var])) for k, env_var in LIMS_DB_CREDENTIAL_MAP.items()) - conn = pg8000.connect( + conn = pg8000.dbapi.connect( user=credentials["user"], host=credentials["host"], database=credentials["dbname"], @@ -70,7 +70,7 @@ def _select(cursor, query, parameters=None): if parameters is None: cursor.execute(query) else: - pg8000.paramstyle = 'numeric' + pg8000.dbapi.paramstyle = 'numeric' cursor.execute(query, parameters) columns = [ to_str(d[0]) for d in cursor.description ] return [ dict(zip(columns, c)) for c in cursor.fetchall() ] @@ -179,9 +179,11 @@ def get_specimen_info_from_lims_by_id(specimen_id): def get_nwb_path_from_lims(ephys_roi_result): """ - Try to find NWBIgor file preferentially - If not found, look for a processed NWB file + Try to find EphysNWB2 file preferentially + If not found, find NWBIgor file + If also not found, look for a processed NWB file + well known file type ID for EphysNWB2 files is 1016154283 well known file type ID for NWB files is 475137571 well known file type ID for NWBIgor files is 570280085 @@ -198,9 +200,15 @@ def get_nwb_path_from_lims(ephys_roi_result): result = query(""" SELECT f.filename, f.storage_directory FROM well_known_files f - WHERE f.attachable_type = 'EphysRoiResult' AND f.attachable_id = %s AND f.well_known_file_type_id = 570280085 + WHERE f.attachable_type = 'EphysRoiResult' AND f.attachable_id = %s AND f.well_known_file_type_id = 1016154283 """ % (ephys_roi_result,)) + if len(result) == 0: + result = query(""" + SELECT f.filename, f.storage_directory FROM well_known_files f + WHERE f.attachable_type = 'EphysRoiResult' AND f.attachable_id = %s AND f.well_known_file_type_id = 570280085 + """ % (ephys_roi_result,)) + if len(result) == 0: logging.warning("Fall back to looking for NWB type") @@ -219,6 +227,8 @@ def get_nwb_path_from_lims(ephys_roi_result): return None + + def get_igorh5_path_from_lims(ephys_roi_result): sql = """ @@ -259,3 +269,47 @@ def project_specimen_ids(project, passed_only=True): results = query(SQL) sp_ids = [d["id"] for d in results] return sp_ids + + +def get_nwb_file_paths_for_specimen_ids(specimen_ids): + """ Get file path for each provided specimen ID + + Note: only returns NWB2 file paths + well known file type ID for EphysNWB2 files is 1016154283 + """ + + sql = """ + select specimens.id, f.filename, f.storage_directory + from specimens + join ephys_roi_results err on err.id = specimens.ephys_roi_result_id + join well_known_files f on f.attachable_id = err.id + where specimens.id = any(:1) + and f.attachable_type = 'EphysRoiResult' + and f.well_known_file_type_id = 1016154283 + """ + try: + result = query(sql, (set(specimen_ids), )) + except pg8000.dbapi.ProgrammingError as e: + print(f"Error Message: {e}") + # Inspect the raw response from the PostgreSQL server + if hasattr(e, 'args') and len(e.args) > 0: + print(f"Server Payload: {e.args}") + + file_list = {r["id"]: os.path.join(r["storage_directory"], r["filename"]) + for r in result} + return file_list + + +def get_sweep_states_and_tags_for_specimens(specimen_ids): + sql = """ + select swp.specimen_id, swp.sweep_number, swp.workflow_state, tag.name as tag_name + from ephys_sweeps swp + left join ephys_sweep_tags_ephys_sweeps estes on estes.ephys_sweep_id = swp.id + left join ephys_sweep_tags tag on tag.id = estes.ephys_sweep_tag_id + where swp.specimen_id = any(:1) + order by swp.specimen_id, swp.sweep_number + """ + result = query(sql, (set(specimen_ids), )) + + return result + diff --git a/ipfx/qc_feature_extractor.py b/ipfx/qc_feature_extractor.py index 05868a34..8b8cb105 100644 --- a/ipfx/qc_feature_extractor.py +++ b/ipfx/qc_feature_extractor.py @@ -95,17 +95,37 @@ def extract_clamp_seal(data_set, tags, manual_values=None): ontology = data_set.ontology + # Try to find the break-in sweeps, so that only sweeps pre-break-in will + # be considered + try: - seal_sweep_number = data_set.get_sweep_numbers(ontology.seal_names,"VoltageClamp")[-1] - seal_data = data_set.sweep(seal_sweep_number) + breakin_sweep_number = data_set.get_sweep_numbers(ontology.breakin_names, "VoltageClamp")[-1] + except: + breakin_sweep_number = None - seal_gohm = qcf.measure_seal(seal_data.v, - seal_data.i, - seal_data.t) + try: + seal_sweep_numbers = data_set.get_sweep_numbers(ontology.seal_names,"VoltageClamp") - if seal_gohm is None or not np.isfinite(seal_gohm): - raise er.FeatureError("Could not compute seal") + if breakin_sweep_number is not None: + seal_sweep_numbers = [s for s in seal_sweep_numbers if s < breakin_sweep_number] + # Find the maximum seal value encountered, in case break-in happened + # in later "cell-attached" sweeps + + seal_values = [] + for sn in seal_sweep_numbers: + seal_data = data_set.sweep(sn) + + seal_gohm = qcf.measure_seal(seal_data.v, + seal_data.i, + seal_data.t) + if seal_gohm is None or not np.isfinite(seal_gohm): + continue + seal_values.append(seal_gohm) + + if len(seal_values) == 0: + raise er.FeatureError("Could not compute seal") + seal_gohm = max(seal_values) except IndexError as e: # seal is not available, for whatever reason. log error tags.append("Seal is not available") @@ -345,11 +365,6 @@ def current_clamp_sweep_qc_features(sweep, is_ramp): current = sweep.i hz = sweep.sampling_rate - expt_start_idx, _ = ep.get_experiment_epoch(current, hz) - # measure noise before stimulus - idx0, idx1 = ep.get_first_noise_epoch(expt_start_idx, hz) # count from the beginning of the experiment - _, qc_features["pre_noise_rms_mv"] = qcf.measure_vm(voltage[idx0:idx1]) - # measure mean and rms of Vm at end of recording # do not check for ramps, because they do not have enough time to recover @@ -359,7 +374,7 @@ def current_clamp_sweep_qc_features(sweep, is_ramp): idx0, idx1 = ep.get_last_stability_epoch(rec_end_idx, hz) mean_last_stability_epoch, _ = qcf.measure_vm(voltage[idx0:idx1]) - idx0, idx1 = ep.get_last_noise_epoch(rec_end_idx, hz) + idx0, idx1 = ep.get_noise_epoch_from_end(rec_end_idx, hz) _, rms_last_noise_epoch = qcf.measure_vm(voltage[idx0:idx1]) else: rms_last_noise_epoch = None @@ -375,6 +390,10 @@ def current_clamp_sweep_qc_features(sweep, is_ramp): idx0, idx1 = ep.get_first_stability_epoch(stim_start_idx, hz) mean_first_stability_epoch, rms_first_stability_epoch = qcf.measure_vm(voltage[idx0:idx1]) + # measure noise before stimulus + idx0, idx1 = ep.get_noise_epoch_from_end(idx1, hz) + _, qc_features["pre_noise_rms_mv"] = qcf.measure_vm(voltage[idx0:idx1]) + qc_features["pre_vm_mv"] = mean_first_stability_epoch qc_features["slow_vm_mv"] = mean_first_stability_epoch qc_features["slow_noise_rms_mv"] = rms_first_stability_epoch @@ -383,4 +402,3 @@ def current_clamp_sweep_qc_features(sweep, is_ramp): return qc_features - diff --git a/ipfx/qc_features.py b/ipfx/qc_features.py index 7db47347..c8ca65d3 100644 --- a/ipfx/qc_features.py +++ b/ipfx/qc_features.py @@ -1,4 +1,5 @@ import numpy as np +from scipy.optimize import curve_fit def measure_blowout(v, idx0): @@ -16,18 +17,32 @@ def measure_electrode_0(curr, hz, t=0.005): else: return None -def measure_seal(v, curr, t): +def measure_seal(v, curr, t, post_transient_shift_ms=0.15): + avg_i, avg_v, slice_t, rel_up_ind, rel_down_ind = average_cell_attached_pulses( + v * 1e-3, curr * 1e-12, t) + return 1e-9 * get_r_from_stable_pulse_response_fit( + avg_v, avg_i, slice_t, + rel_up_ind, rel_down_ind, post_transient_shift_ms=post_transient_shift_ms) - return 1e-9 * get_r_from_stable_pulse_response(v*1e-3, curr*1e-12, t) - -def measure_input_resistance(v, curr, t): - - return 1e-6 * get_r_from_stable_pulse_response(v*1e-3, curr*1e-12, t) +def measure_input_resistance(v, curr, t, post_transient_shift_ms=1.0): + avg_i, avg_v, slice_t, rel_up_ind, rel_down_ind = average_whole_cell_pulses( + v * 1e-3, curr * 1e-12, t) + r_a = 1e-6 * get_r_from_peak_pulse_response( + avg_v, avg_i, slice_t, + rel_up_ind, rel_down_ind) + r_tot = 1e-6 * get_r_from_stable_pulse_response_fit( + avg_v, avg_i, slice_t, + rel_up_ind, rel_down_ind, post_transient_shift_ms=post_transient_shift_ms) + return r_tot - r_a def measure_initial_access_resistance(v, curr, t): - return 1e-6 * get_r_from_peak_pulse_response(v*1e-3, curr*1e-12, t) + avg_i, avg_v, slice_t, rel_up_ind, rel_down_ind = average_whole_cell_pulses( + v * 1e-3, curr * 1e-12, t) + return 1e-6 * get_r_from_peak_pulse_response( + avg_v, avg_i, slice_t, + rel_up_ind, rel_down_ind) def measure_vm(vals): @@ -47,70 +62,75 @@ def measure_vm_delta(mean_start, mean_end): return None -def get_r_from_stable_pulse_response(v, i, t): - """Compute input resistance from the stable pulse response - - Parameters - ---------- - v : float membrane voltage (V) - i : float input current (A) - t : time (s) +def get_r_from_stable_pulse_response_fit(avg_v, avg_i, t, + relative_up_ind, relative_down_ind, post_transient_shift_ms): + if avg_i is None: + # Not enough pulses to average + return np.nan - Returns - ------- - ir: float input resistance - """ - - up_idx, down_idx = get_square_pulse_idx(v) dt = t[1] - t[0] one_ms = int(0.001 / dt) - r = [] - for ii in range(len(up_idx)): - # take average v and i one ms before start - end = up_idx[ii] - 1 - start = end - one_ms + post_transient_shift = int(one_ms * post_transient_shift_ms) - avg_v_base = np.mean(v[start:end]) - avg_i_base = np.mean(i[start:end]) + # baseline - take average v and i one ms before start + end = relative_up_ind - 1 + start = end - one_ms + avg_v_base = np.mean(avg_v[start:end]) + avg_i_base = np.mean(avg_i[start:end]) - # take average v and i one ms before end - end = down_idx[ii]-1 - start = end - one_ms + t_window = (t[relative_up_ind + post_transient_shift:relative_down_ind] - + t[relative_up_ind + post_transient_shift]) + i_window = avg_i[relative_up_ind + post_transient_shift:relative_down_ind] * 1e12 + guess = ( + max(i_window[0] - i_window[-1], 100), + 1e3, + max(i_window[-1], avg_i_base * 1e12 + 5) + ) - avg_v_steady = np.mean(v[start:end]) - avg_i_steady = np.mean(i[start:end]) + popt, pcov = curve_fit( + _exp_curve, t_window, i_window, + p0=guess, bounds=([0, 0, avg_i_base * 1e12], [np.inf, np.inf, np.inf]) + ) + pred = _exp_curve(t_window, *popt) * 1e-12 - r_instance = (avg_v_steady-avg_v_base) / (avg_i_steady-avg_i_base) + # steady-state - take average v and i one ms before end + end = relative_down_ind - 1 + start = end - one_ms + avg_v_steady = np.mean(avg_v[start:end]) + avg_i_steady = np.mean(avg_i[start:end]) + avg_i_steady = np.mean(pred[-one_ms:]) - r.append(r_instance) + r = (avg_v_steady - avg_v_base) / (avg_i_steady - avg_i_base) - return np.mean(r) + return r -def get_r_from_peak_pulse_response(v, i, t): - - up_idx, down_idx = get_square_pulse_idx(v) +def get_r_from_peak_pulse_response(avg_v, avg_i, t, + relative_up_ind, relative_down_ind): + if avg_i is None: + # Not enough pulses to average + return np.nan dt = t[1] - t[0] one_ms = int(0.001 / dt) - r = [] - for ii in range(len(up_idx)): - # take average v and i one ms before - end = up_idx[ii] - 1 - start = end - one_ms - avg_v_base = np.mean(v[start:end]) - avg_i_base = np.mean(i[start:end]) - # take average v and i one ms before end - start = up_idx[ii] - end = down_idx[ii] - 1 - idx = start + np.argmax(i[start:end]) - avg_v_peak = v[idx] - avg_i_peak = i[idx] - r_instance = (avg_v_peak-avg_v_base) / (avg_i_peak-avg_i_base) - r.append(r_instance) - return np.mean(r) + # take average v and i one ms before start of pulse + end = relative_up_ind - 1 + start = end - one_ms + avg_v_base = np.mean(avg_v[start:end]) + avg_i_base = np.mean(avg_i[start:end]) + + # find peak i during the pulse + start = relative_up_ind + end = relative_down_ind - 1 + idx = start + np.argmax(avg_i[start:end]) + avg_v_peak = avg_v[idx] + avg_i_peak = avg_i[idx] + r = (avg_v_peak - avg_v_base) / (avg_i_peak - avg_i_base) + + return r + def get_square_pulse_idx(v): @@ -139,3 +159,160 @@ def get_square_pulse_idx(v): assert up_ix < down_ix, "Negative square pulse" return up_idx, down_idx + + +def average_pulses(v, i, t, cell_attached, access_cutoff_for_breakin=100, rmse_cutoff=50e-12): + up_idx, down_idx = get_square_pulse_idx(v) + if len(up_idx) == 0: + return None, None, None, None, None + + dt = t[1] - t[0] + one_ms = int(0.001 / dt) + + pulses_i_all = [] + pulses_i_no_spikes = [] + pulses_i_no_late_spikes = [] + pulses_v = [] + + for u, d in zip(up_idx, down_idx): + start_idx = u - one_ms * 2 + end_idx = d + one_ms * 2 + slice_t = dt * np.arange(end_idx - start_idx) + relative_up_ind = u - start_idx + relative_down_ind = d - start_idx + + end = relative_up_ind - 1 + start = end - one_ms + i_baseline = np.mean(i[start_idx:end_idx][start:end]) + v_baseline = np.mean(v[start_idx:end_idx][start:end]) + + # Check for noisy baseline + rmse = np.sqrt(((i_baseline - i[start_idx:end_idx][start:end]) ** 2).mean()) + if rmse > rmse_cutoff: + continue + + if cell_attached: + # Check if it's already broken in + peak_cap_ind = np.argmax(i[start_idx:end_idx][relative_up_ind:relative_up_ind + int(one_ms * 0.2)]) + relative_up_ind + i_peak = i[start_idx:end_idx][peak_cap_ind] + v_peak = v[start_idx:end_idx][peak_cap_ind] + + est_access_resistance = 1e6 * (v_peak - v_baseline) / (i_peak - i_baseline) + if est_access_resistance < access_cutoff_for_breakin: + break + + # check for contaminating spikes + if cell_attached: + spikes = detect_cell_attached_spikes( + i[start_idx:end_idx][relative_up_ind:relative_down_ind] - i_baseline, + slice_t[relative_up_ind:relative_down_ind]) + else: + spikes = detect_escaped_spikes( + i[start_idx:end_idx][relative_up_ind:relative_down_ind] - i_baseline, + slice_t[relative_up_ind:relative_down_ind]) + + pulses_v.append(v[start_idx:end_idx]) + pulses_i_all.append(i[start_idx:end_idx]) + if len(spikes) == 0: + pulses_i_no_spikes.append(i[start_idx:end_idx]) + else: + spike_inds = np.array([s[0] for s in spikes]) + if not np.any(spike_inds > relative_down_ind - relative_up_ind - 2 * one_ms): + pulses_i_no_late_spikes.append(i[start_idx:end_idx]) + + + if len(pulses_i_no_spikes) > 0: +# print(f"Using {len(pulses_i_no_spikes)} pulses with no spikes") + avg_i = np.vstack(pulses_i_no_spikes).mean(axis=0) + elif len(pulses_i_no_late_spikes) > 0: +# print(f"Using {len(pulses_i_no_late_spikes)} pulses with no late spikes") + avg_i = np.vstack(pulses_i_no_late_spikes).mean(axis=0) + elif len(pulses_i_all) > 0: +# print(f"Using {len(pulses_i_no_late_spikes)} pulses; all had late spikes") + avg_i = np.vstack(pulses_i_all).mean(axis=0) + else: + avg_i = None + + if len(pulses_v) > 0: + avg_v = np.vstack(pulses_v).mean(axis=0) + else: + avg_v = None + + return avg_i, avg_v, slice_t, relative_up_ind, relative_down_ind + + +def average_cell_attached_pulses(v, i, t, access_cutoff_for_breakin=100, rmse_cutoff=50e-12): + return average_pulses(v, i, t, + cell_attached=True, + access_cutoff_for_breakin=access_cutoff_for_breakin, + rmse_cutoff=rmse_cutoff + ) + + +def average_whole_cell_pulses(v, i, t, rmse_cutoff=50e-12): + return average_pulses(v, i, t, + cell_attached=False, + rmse_cutoff=rmse_cutoff + ) + + +def detect_cell_attached_spikes(i, t, min_spike_amp = -30e-12): + # i should be baselined + dt = t[1] - t[0] + one_ms = int(0.001 / dt) + + # start detection after peak of capacitance transient + cap_peak_ind = np.argmax(i[0:int(one_ms * 0.2)]) + + putative_spikes = np.flatnonzero(np.diff(np.less_equal(i[cap_peak_ind:], min_spike_amp).astype(int)) == 1) + cap_peak_ind + + last_spike_t = t[0] - 0.001 + spikes = [] + for spike_ind in putative_spikes: + if t[spike_ind] < last_spike_t + 0.001: + # too close to previous spike + continue + + # find peak in 1 ms window + peak_ind = np.argmin(i[spike_ind:spike_ind + one_ms]) + spike_ind + peak_amp = i[peak_ind] + + # check for biphasic - does it go at least 50% of min amplitude above baseline within a millisecond of negative-going peak? + pos_peak_ind = np.argmax(i[peak_ind:peak_ind + one_ms]) + peak_ind + pos_peak_amp = i[pos_peak_ind] + if i[pos_peak_ind] >= 0.5 * -min_spike_amp: + # count it as a spike + spikes.append((peak_ind, peak_amp)) + last_spike_t = t[peak_ind] + + return spikes + + + +def detect_escaped_spikes(i, t, min_spike_amp = -500e-12): + # i should be baselined + dt = t[1] - t[0] + one_ms = int(0.001 / dt) + + # start detection after peak of capacitance transient + cap_peak_ind = np.argmax(i[0:int(one_ms * 0.2)]) + + putative_spikes = np.flatnonzero(np.diff(np.less_equal(i[cap_peak_ind:], min_spike_amp).astype(int)) == 1) + cap_peak_ind + + last_spike_t = t[0] - 0.001 + spikes = [] + for spike_ind in putative_spikes: + if t[spike_ind] < last_spike_t + 0.001: + # too close to previous spike + continue + + # find peak in 1 ms window + peak_ind = np.argmin(i[spike_ind:spike_ind + one_ms]) + spike_ind + peak_amp = i[peak_ind] + spikes.append((peak_ind, peak_amp)) + + return spikes + + +def _exp_curve(x, a, inv_tau, y0): + return y0 + a * np.exp(-inv_tau * x) diff --git a/ipfx/script_utils.py b/ipfx/script_utils.py index 042433e7..5aac8616 100755 --- a/ipfx/script_utils.py +++ b/ipfx/script_utils.py @@ -7,16 +7,23 @@ import pandas as pd import h5py +from typing import NamedTuple import ipfx.lims_queries as lq import ipfx.stim_features as stf import ipfx.stimulus_protocol_analysis as spa import ipfx.data_set_features as dsf import ipfx.time_series_utils as tsu import ipfx.error as er +import ipfx.qc_feature_extractor as qc_fex +import ipfx.qc_feature_evaluator as qc_feval from ipfx.stimulus import StimulusType from ipfx.sweep import SweepSet from ipfx.dataset.create import create_ephys_data_set +class StimulusTiming(NamedTuple): + start: float + end: float + dur: float def lims_nwb_information(specimen_id): _, roi_id, _ = lq.get_specimen_info_from_lims_by_id(specimen_id) @@ -74,94 +81,203 @@ def dataset_for_specimen_id(specimen_id, data_source, ontology, file_list=None): return data_set -def categorize_iclamp_sweeps(data_set, stimuli_names, sweep_qc_option="none", specimen_id=None): - exist_sql = """ - select swp.sweep_number from ephys_sweeps swp - where swp.specimen_id = :1 - and swp.sweep_number = any(:2) - """ +def categorize_iclamp_sweeps(data_set, stimuli_names, sweep_qc_record, + sweep_qc_option="none", specimen_id=None): - passed_sql = """ - select swp.sweep_number from ephys_sweeps swp - where swp.specimen_id = :1 - and swp.sweep_number = any(:2) - and swp.workflow_state like '%%passed' - """ - - passed_except_delta_vm_sql = """ - select swp.sweep_number, tag.name - from ephys_sweeps swp - join ephys_sweep_tags_ephys_sweeps estes on estes.ephys_sweep_id = swp.id - join ephys_sweep_tags tag on tag.id = estes.ephys_sweep_tag_id - where swp.specimen_id = :1 - and swp.sweep_number = any(:2) - """ - - iclamp_st = data_set.filtered_sweep_table(clamp_mode=data_set.CURRENT_CLAMP, stimuli=stimuli_names) + my_sweep_qc_record = sweep_qc_record.loc[sweep_qc_record["specimen_id"] == specimen_id] + iclamp_st = data_set.filtered_sweep_table( + clamp_mode=data_set.CURRENT_CLAMP, stimuli=stimuli_names) if iclamp_st.shape[0] == 0: return np.array([]) + sweep_num_list = iclamp_st["sweep_number"].sort_values().unique().tolist() if sweep_qc_option == "none": - return iclamp_st["sweep_number"].sort_values().values - elif sweep_qc_option == "lims-passed-only": - # check that sweeps exist in LIMS - sweep_num_list = iclamp_st["sweep_number"].sort_values().tolist() - results = lq.query(exist_sql, (specimen_id, sweep_num_list)) - res_nums = pd.DataFrame(results, columns=["sweep_number"])["sweep_number"].tolist() + return np.array(sweep_num_list) + elif sweep_qc_option in ("passed-only", "passed-except-delta-vm", "passed-except-delta-vm-and-rms"): + # check that sweeps exist in sweep QC record not_checked_list = [] for swp_num in sweep_num_list: - if swp_num not in res_nums: - logging.debug("Could not find sweep {:d} from specimen {:d} in LIMS for QC check".format(swp_num, specimen_id)) + if swp_num not in my_sweep_qc_record["sweep_number"].unique(): not_checked_list.append(swp_num) + if len(not_checked_list) > 0: + sweep_num_list = [sn for sn in sweep_num_list if sn not in not_checked_list] + logging.warning(f"Could not find {len(not_checked_list)} sweeps from specimen {specimen_id} in QC record ({stimuli_names})") + # note: choosing not to include unchecked sweeps in returned list # Get passed sweeps - results = lq.query(passed_sql, (specimen_id, sweep_num_list)) - results_df = pd.DataFrame(results, columns=["sweep_number"]) - passed_sweep_nums = results_df["sweep_number"].values - return np.sort(np.hstack([passed_sweep_nums, np.array(not_checked_list)])) # deciding to keep non-checked sweeps for now - elif sweep_qc_option == "lims-passed-except-delta-vm": - # check that sweeps exist in LIMS - sweep_num_list = iclamp_st["sweep_number"].sort_values().tolist() - results = lq.query(exist_sql, (specimen_id, sweep_num_list)) - res_nums = pd.DataFrame(results, columns=["sweep_number"])["sweep_number"].tolist() - - not_checked_list = [] - for swp_num in sweep_num_list: - if swp_num not in res_nums: - logging.debug("Could not find sweep {:d} from specimen {:d} in LIMS for QC check".format(swp_num, specimen_id)) - not_checked_list.append(swp_num) + passed_record = my_sweep_qc_record.loc[ + my_sweep_qc_record["sweep_number"].isin(sweep_num_list) & + my_sweep_qc_record["workflow_state"].str.endswith("passed"), :] + passed_sweep_nums = passed_record["sweep_number"].unique() - # get straight-up passed sweeps - results = lq.query(passed_sql, (specimen_id, sweep_num_list)) - results_df = pd.DataFrame(results, columns=["sweep_number"]) - passed_sweep_nums = results_df["sweep_number"].values + if sweep_qc_option == "passed-only": + return np.sort(passed_sweep_nums).astype(int) # also get sweeps that only fail due to delta Vm failed_sweep_list = list(set(sweep_num_list) - set(passed_sweep_nums)) if len(failed_sweep_list) == 0: - return np.sort(passed_sweep_nums) - results = lq.query(passed_except_delta_vm_sql, (specimen_id, failed_sweep_list)) - results_df = pd.DataFrame(results, columns=["sweep_number", "name"]) - - # not all cells have tagged QC status - if there are no tags assume the - # fail call is correct and exclude those sweeps - tagged_mask = np.array([sn in results_df["sweep_number"].tolist() for sn in failed_sweep_list]) + return np.sort(passed_sweep_nums).astype(int) + + # check if only tag is "Vm delta" + also_passing_nums = [] + for sn in failed_sweep_list: + non_delta_vm_tag_record = my_sweep_qc_record.loc[ + (my_sweep_qc_record["sweep_number"] == sn) & + (~my_sweep_qc_record["tag_name"].str.startswith("Vm delta")) & + (my_sweep_qc_record["tag_name"] != "Blowout is not available"), :] # don't fail for blowout unavailable because we are considering patch-seq sweeps + if non_delta_vm_tag_record.shape[0] == 0: + also_passing_nums.append(sn) + + if sweep_qc_option == "passed-except-delta-vm": + return np.sort(np.hstack([ + passed_sweep_nums, + np.array(also_passing_nums), + ])).astype(int) + + # Don't use LIMS-calculated RMS fail/pass - recalculate here # otherwise, check for having an error tag that isn't 'Vm delta' - # and exclude those sweeps - has_non_delta_tags = np.array([np.any((results_df["sweep_number"].values == sn) & - (results_df["name"].values != "Vm delta")) for sn in failed_sweep_list]) - - also_passing_nums = np.array(failed_sweep_list)[tagged_mask & ~has_non_delta_tags] + # or one of the RMS tags and exclude those sweeps + rms_check_sweep_nums = [] + for sn in set(failed_sweep_list) - set(also_passing_nums): + non_delta_vm_or_rms_tag_record = my_sweep_qc_record.loc[ + (my_sweep_qc_record["sweep_number"] == sn) & + (~my_sweep_qc_record["tag_name"].str.startswith("Vm delta")) & + (~my_sweep_qc_record["tag_name"].str.startswith("slow noise")) & + (~my_sweep_qc_record["tag_name"].str.startswith("pre-noise")) & + (~my_sweep_qc_record["tag_name"].str.startswith("post-noise")) & + (my_sweep_qc_record["tag_name"] != "Blowout is not available"), # don't fail for blowout unavailable because we are considering patch-seq sweeps + :] + if non_delta_vm_or_rms_tag_record.shape[0] == 0: + rms_check_sweep_nums.append(sn) + + if len(rms_check_sweep_nums) == 0: + # if no sweeps need to be checked, skip the rest + return np.sort(np.hstack([ + passed_sweep_nums, + np.array(also_passing_nums), + ])).astype(int) + + # Now re-check each sweep's RMS + qc_criteria = qc_feval.load_default_qc_criteria() + + # Read the lab notebook for the RMS criteria used for the sweep + lnr = data_set._data.notebook + numeric_fields = [c.decode('utf-8') for c in lnr.colname_number[0]] + short_rms_fields = [f for f in numeric_fields if "S-RMS Threshold" in f] + long_rms_fields = [f for f in numeric_fields if "L-RMS Threshold" in f] + + pass_rms_nums = [] + for sn in rms_check_sweep_nums: + is_ramp = "Ramp" == iclamp_st.at[sn, "stimulus_name"] + + # Short RMS criterion + s_rms_threshold = None + for f in short_rms_fields: + if lnr.get_value(f, sn, None) is not None: + s_rms_threshold = lnr.get_value(f, sn, None) * 1e3 # from V to mV + break + if s_rms_threshold is None: + s_rms_threshold = qc_criteria["pre_noise_rms_mv_max"] + + # Long RMS criterion + l_rms_threshold = None + for f in long_rms_fields: + if lnr.get_value(f, sn, None) is not None: + l_rms_threshold = lnr.get_value(f, sn, None) * 1e3 # from V to mV + break + if l_rms_threshold is None: + l_rms_threshold = qc_criteria["slow_noise_rms_mv_max"] + + qc_features = qc_fex.current_clamp_sweep_qc_features( + data_set.sweep(sn), + is_ramp + ) + + if is_ramp: + if ((qc_features["pre_noise_rms_mv"] < s_rms_threshold) & + (qc_features["slow_noise_rms_mv"] < l_rms_threshold)): + pass_rms_nums.append(sn) + else: + if ((qc_features["pre_noise_rms_mv"] < s_rms_threshold) & + (qc_features["post_noise_rms_mv"] < s_rms_threshold) & + (qc_features["slow_noise_rms_mv"] < l_rms_threshold)): + pass_rms_nums.append(sn) + + if sweep_qc_option == "passed-except-delta-vm-and-rms": + return np.sort(np.hstack([ + passed_sweep_nums, + np.array(also_passing_nums), + np.array(pass_rms_nums), + ])).astype(int) - return np.sort(np.hstack([passed_sweep_nums, also_passing_nums, np.array(not_checked_list)])) else: raise ValueError("Invalid sweep-level QC option {}".format(sweep_qc_option)) 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 = [] + stim_timing = [] + for swp in check_sweeps.sweeps: + if len(swp.t) == 0: + valid_sweep_stim.append(False) + continue + + swp_start, swp_dur, _, _, _ = stf.get_stim_characteristics(swp.i, swp.t) + if swp_start is None: + valid_sweep_stim.append(False) + stim_timing.append(None) + else: + valid_sweep_stim.append(True) + stim_timing.append(StimulusTiming( + start=swp_start, + end=swp_start + swp_dur, + dur=swp_dur + )) + if len(stim_timing) == 0: + # Could not find any sweeps to define stimulus interval + return None, None + + + # Check that all sweeps are long enough and not ended early + good_sweeps = [] + good_stim_timing = [] + for s, v, swp_stim_timing in zip(check_sweeps.sweeps, valid_sweep_stim, stim_timing): + if not v: + logging.debug(f"Sweep {s.sweep_number} not valid stim") + continue + if s.t[-1] < swp_stim_timing.end + extra_dur: + logging.debug(f"Sweep {s.sweep_number} not long enough after end") + continue + if np.all(s.v[tsu.find_time_index(s.t, swp_stim_timing.end) - 100:tsu.find_time_index(s.t, swp_stim_timing.end)] == 0): + logging.debug(f"Sweep {s.sweep_number} end of stim interval was all zero") + continue + good_sweeps.append(s) + good_stim_timing.append(swp_stim_timing) + + if len(good_sweeps) == 0: + return None, None + + # Check for consistent stimulus intervals + + if not np.all(np.isclose([s.dur for s in good_stim_timing], good_stim_timing[0].dur)): + logging.warning("Sweeps in set do not all have the same duration") + + if not np.all(np.isclose([s.start for s in good_stim_timing], good_stim_timing[0].start)): + logging.debug("Stimulus start times are not identical across sweeps in set") + + if not np.all(np.isclose([s.end for s in good_stim_timing], good_stim_timing[0].end)): + logging.debug("Stimulus end times are not identical across sweeps in set") + + return SweepSet(sweeps=good_sweeps), good_stim_timing + + +def validate_ramp_sweeps(data_set, sweep_numbers, min_ramp_dur=0.1): + check_sweeps = data_set.sweep_set(sweep_numbers) + check_sweeps.select_epoch("recording") valid_sweep_stim = [] start = None dur = None @@ -171,7 +287,9 @@ def validate_sweeps(data_set, sweep_numbers, extra_dur=0.2): continue swp_start, swp_dur, _, _, _ = stf.get_stim_characteristics(swp.i, swp.t) - if swp_start is None: + if swp_start is None or swp_dur is None: + valid_sweep_stim.append(False) + elif swp_dur < min_ramp_dur: valid_sweep_stim.append(False) else: start = swp_start @@ -179,30 +297,27 @@ def validate_sweeps(data_set, sweep_numbers, extra_dur=0.2): valid_sweep_stim.append(True) if start is None: # Could not find any sweeps to define stimulus interval - return [], None, None - - end = start + dur + return None - # Check that all sweeps are long enough and not ended early + # Check that all sweeps are long enough and did not end early good_sweeps = [s for s, v in zip(check_sweeps.sweeps, valid_sweep_stim) - if s.t[-1] >= end + extra_dur - and v is True - and not np.all(s.v[tsu.find_time_index(s.t, end)-100:tsu.find_time_index(s.t, end)] == 0)] - return SweepSet(sweeps=good_sweeps), start, end + if v is True] + return SweepSet(sweeps=good_sweeps) def preprocess_long_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, subthresh_min_amp=-100.): if len(sweep_numbers) == 0: raise er.FeatureError("No long square sweeps available for feature extraction") - lsq_sweeps, lsq_start, lsq_end = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur) + lsq_sweeps, lsq_stim_timing = 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, - start=lsq_start, - end=lsq_end, + start=[s.start for s in lsq_stim_timing], + end=[s.end for s in lsq_stim_timing], min_peak=-25, **dsf.detection_parameters(StimulusType.LONG_SQUARE) ) @@ -210,21 +325,32 @@ def preprocess_long_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, subthr subthresh_min_amp=subthresh_min_amp) lsq_features = lsq_an.analyze(lsq_sweeps) - return lsq_sweeps, lsq_features, lsq_an, lsq_start, lsq_end + return lsq_sweeps, lsq_features, lsq_an, lsq_stim_timing def preprocess_short_square_sweeps(data_set, sweep_numbers, extra_dur=0.2, spike_window=0.05): if len(sweep_numbers) == 0: raise er.FeatureError("No short square sweeps available for feature extraction") - ssq_sweeps, ssq_start, ssq_end = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur) - if len(ssq_sweeps.sweeps) == 0: + ssq_sweeps, ssq_stim_timing = validate_sweeps(data_set, sweep_numbers, extra_dur=extra_dur) + if ssq_sweeps is None or 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") + est_window = [ + [s.start for s in ssq_stim_timing], + [s.start + 0.001 for s in ssq_stim_timing] + ] + extractor_end = [] + for idx, swp in enumerate(ssq_sweeps.sweeps): + if swp.t[-1] < ssq_stim_timing[idx].end + spike_window: + extractor_end.append(swp.t[-1] - 0.001) + else: + extractor_end.append(ssq_stim_timing[idx].end + spike_window) ssq_spx, ssq_spfx = dsf.extractors_for_sweeps(ssq_sweeps, - est_window = [ssq_start, ssq_start + 0.001], - start=ssq_start, - end=ssq_end + spike_window, + est_window=est_window, + start=[s.start for s in ssq_stim_timing], + end=extractor_end, reject_at_stim_start_interval=0.0002, **dsf.detection_parameters(StimulusType.SHORT_SQUARE)) ssq_an = spa.ShortSquareAnalysis(ssq_spx, ssq_spfx) @@ -237,9 +363,23 @@ def preprocess_ramp_sweeps(data_set, sweep_numbers): if len(sweep_numbers) == 0: raise er.FeatureError("No ramp sweeps available for feature extraction") - ramp_sweeps = data_set.sweep_set(sweep_numbers) - - ramp_start, ramp_dur, _, _, _ = stf.get_stim_characteristics(ramp_sweeps.sweeps[0].i, ramp_sweeps.sweeps[0].t) + ramp_sweeps = validate_ramp_sweeps(data_set, sweep_numbers) + if ramp_sweeps is None or len(ramp_sweeps.sweeps) == 0: + raise er.FeatureError("No ramp sweeps were long enough") + ramp_sweeps.select_epoch("recording") + + starts = [] + durs = [] + for swp in ramp_sweeps.sweeps: + ramp_start, ramp_dur, _, _, _ = stf.get_stim_characteristics( + swp.i, swp.t) + starts.append(ramp_start) + durs.append(ramp_dur) + + if np.all(np.isclose(starts, starts[0])): + ramp_start = starts[0] + else: + ramp_start = starts ramp_spx, ramp_spfx = dsf.extractors_for_sweeps(ramp_sweeps, start = ramp_start, **dsf.detection_parameters(StimulusType.RAMP)) @@ -251,16 +391,16 @@ def preprocess_ramp_sweeps(data_set, sweep_numbers): def filter_results(specimen_ids, results): filtered_set = [(i, r) for i, r in zip(specimen_ids, results) if not "error" in r.keys()] - error_set = [{"id": i, "error": d} for i, d in zip(specimen_ids, results) if "error" in d.keys()] + error_set = [d for d in results if "error" in d] if len(filtered_set) == 0: logging.info("No specimens had results") - return + return None, None, None used_ids, results = zip(*filtered_set) return used_ids, results, error_set -def organize_results(specimen_ids, results): +def organize_results(specimen_ids, results, skip_keys=[]): """Build dictionary of results, filling data from cells with appropriate-length nan arrays where needed""" result_sizes = {} @@ -268,10 +408,15 @@ def organize_results(specimen_ids, results): all_keys = np.unique(np.concatenate([list(r.keys()) for r in results])) for k in all_keys: + if k in skip_keys: + continue if k not in result_sizes: - for r in results: + for r, sp_id in zip(results, specimen_ids): if k in r and r[k] is not None: - result_sizes[k] = len(r[k]) + if k not in result_sizes: + result_sizes[k] = len(r[k]) + elif len(r[k]) != result_sizes[k]: + logging.warning(f"found result with length {len(r[k])} when expecting length {result_sizes[k]} for {k}; specimen ID {sp_id}") data = np.array([r[k] if k in r else np.nan * np.zeros(result_sizes[k]) for r in results]) output[k] = data diff --git a/ipfx/spike_features.py b/ipfx/spike_features.py index 728d4be8..fd4d673a 100644 --- a/ipfx/spike_features.py +++ b/ipfx/spike_features.py @@ -355,13 +355,20 @@ def estimate_adjusted_detection_parameters(v_set, t_set, interval_start, interva if len(v_set) == 0: raise er.FeatureError("t_set and v_set are empty") - start_index = tsu.find_time_index(t_set[0], interval_start) - end_index = tsu.find_time_index(t_set[0], interval_end) - maxes = [] ends = [] dv_set = [] - for v, t in zip(v_set, t_set): + for idx, (v, t) in enumerate(zip(v_set, t_set)): + if type(interval_start) is list: + start_t = interval_start[idx] + end_t = interval_end[idx] + else: + start_t = interval_start + end_t = interval_end + + start_index = tsu.find_time_index(t_set[idx], start_t) + end_index = tsu.find_time_index(t_set[idx], end_t) + dv = tsu.calculate_dvdt(v, t, filter) dv_set.append(dv) maxes.append(dv[start_index:end_index].max()) diff --git a/ipfx/stim_features.py b/ipfx/stim_features.py index 01ff04c4..a21073e8 100644 --- a/ipfx/stim_features.py +++ b/ipfx/stim_features.py @@ -24,8 +24,8 @@ def get_stim_characteristics(i, t, test_pulse=True): stim = i[start_idx:end_idx+1] - peak_high = max(stim) - peak_low = min(stim) + peak_high = np.max(stim) + peak_low = np.min(stim) if abs(peak_high) > abs(peak_low): amplitude = float(peak_high) diff --git a/ipfx/stimulus_protocol_analysis.py b/ipfx/stimulus_protocol_analysis.py index 6fb4978b..0ae58820 100644 --- a/ipfx/stimulus_protocol_analysis.py +++ b/ipfx/stimulus_protocol_analysis.py @@ -70,11 +70,12 @@ def mean_features_first_spike(self, spikes_set, features_list=None): def analyze_basic_features(self, sweep_set, extra_sweep_features=None, exclude_clipped=False): self._spikes_set = [] - for sweep in sweep_set.sweeps: - self._spikes_set.append(self.spx.process(sweep.t, sweep.v, sweep.i)) + for idx, sweep in enumerate(sweep_set.sweeps): + self._spikes_set.append(self.spx.process(sweep.t, sweep.v, sweep.i, sweep_index=idx)) - self._sweep_features = pd.DataFrame([ self.sptx.process(sweep.t, sweep.v, sweep.i, spikes, extra_sweep_features, exclude_clipped=exclude_clipped) - for sweep, spikes in zip(sweep_set.sweeps, self._spikes_set) ]) + self._sweep_features = pd.DataFrame([ + self.sptx.process(sweep.t, sweep.v, sweep.i, spikes, extra_sweep_features, exclude_clipped=exclude_clipped, sweep_index=idx) + for idx, (sweep, spikes) in enumerate(zip(sweep_set.sweeps, self._spikes_set))]) def reset_basic_features(self): self._spikes_set = None @@ -192,6 +193,9 @@ def analyze_subthreshold(self, sweep_set): calc_subthresh_features = subthreshold_sweep_features[ (subthreshold_sweep_features["stim_amp"] < self.SUBTHRESH_MAX_AMP) & \ (subthreshold_sweep_features["stim_amp"] > self.subthresh_min_amp) ].copy() + # only use the last-acquired sweeps + calc_subthresh_features.drop_duplicates("stim_amp", keep="last", inplace=True) + if len(calc_subthresh_features) == 0: error_string = F"No subthreshold long square sweeps with stim_amp " \ F"in range [{self.subthresh_min_amp,self.SUBTHRESH_MAX_AMP}] " \ @@ -203,9 +207,19 @@ def analyze_subthreshold(self, sweep_set): calc_subthresh_ss = SweepSet([sweep_set.sweeps[i] for i in calc_subthresh_features.index.values]) - median_peak_time = np.median([s.t[subf.voltage_deflection(s.t, s.v, s.i, self.spx.start, self.spx.end, "min")[1]] - for s in calc_subthresh_ss.sweeps]) - taus = [ subf.time_constant(s.t, s.v, s.i, self.spx.start, self.spx.end, median_peak_time, self.tau_frac, self.sptx.baseline_interval) for s in calc_subthresh_ss.sweeps ] + + # Handle non-identical start/end times + if type(self.spx.start) is list: + starts = [self.spx.start[i] for i in calc_subthresh_features.index.values] + ends = [self.spx.end[i] for i in calc_subthresh_features.index.values] + else: + starts = [self.spx.start] * len(calc_subthresh_ss.sweeps) + ends = [self.spx.end] * len(calc_subthresh_ss.sweeps) + + median_peak_time = np.median([s.t[subf.voltage_deflection(s.t, s.v, s.i, start, end, "min")[1]] + for s, start, end in zip(calc_subthresh_ss.sweeps, starts, ends)]) + taus = [ subf.time_constant(s.t, s.v, s.i, start, end, median_peak_time, self.tau_frac, self.sptx.baseline_interval) + for s, start, end in zip(calc_subthresh_ss.sweeps, starts, ends)] calc_subthresh_features['tau'] = taus @@ -213,7 +227,7 @@ def analyze_subthreshold(self, sweep_set): features["input_resistance"] = subf.input_resistance(calc_subthresh_ss.t, calc_subthresh_ss.i, calc_subthresh_ss.v, - self.spx.start, self.spx.end, + starts, ends, self.sptx.baseline_interval) features["tau"] = np.nanmean(calc_subthresh_features['tau']) diff --git a/ipfx/subthresh_features.py b/ipfx/subthresh_features.py index 7010e47b..7e463c1c 100644 --- a/ipfx/subthresh_features.py +++ b/ipfx/subthresh_features.py @@ -2,6 +2,7 @@ import logging from . import time_series_utils as tsu from scipy.optimize import curve_fit +from scipy.signal import savgol_filter from . import error as er def baseline_voltage(t, v, start, baseline_interval=0.1, baseline_detect_thresh=0.3, filter_frequency=1.0): @@ -27,7 +28,7 @@ def baseline_voltage(t, v, start, baseline_interval=0.1, baseline_detect_thresh= return np.nan -def voltage_deflection(t, v, i, start, end, deflect_type=None): +def voltage_deflection(t, v, i, start, end, deflect_type=None, reject_transients=False, smoothing=True): """Measure deflection (min or max, between start and end if specified). Parameters @@ -63,9 +64,85 @@ def voltage_deflection(t, v, i, start, end, deflect_type=None): deflect_func = deflect_dispatch[deflect_type] - v_window = v[start_index:end_index] + if smoothing: + v_smooth = savgol_filter(v, 40, 2) + v_window = v_smooth[start_index:end_index] + else: + v_window = v[start_index:end_index] + deflect_index = deflect_func(v_window) + start_index + # Try to automatically detect and reject transients if requested + # - look near the peak for high dv/dt values and block them out + if reject_transients: + nan_deflect_dispatch = { + "min": np.nanargmin, + "max": np.nanargmax, + } + nan_deflect_func = nan_deflect_dispatch[deflect_type] + + # use an adaptive threshold to avoid treating the start of the step as a transient + dvdt_thresh = 1.0 # mV/ms + t_envelope = 1e3 * (t[start_index:end_index] - t[start_index]) + dvdt_thresh_envelope = 5 * np.exp(-t_envelope / 2) + dvdt_thresh + + window_width = 400 + dvdt = savgol_filter(v, 50, 2, deriv=1, delta=1e3 * (t[1] - t[0])) # mV/ms, smoothed + if smoothing: + temp_v = v_smooth.copy() + else: + temp_v = v.copy() + + window_start = max(deflect_index - window_width, start_index) + window_end = min(deflect_index + window_width, end_index) + dvdt_window = dvdt[window_start:window_end] + peak_dvdt_ind = np.argmax(np.abs(dvdt_window)) + peak_dvdt = dvdt_window[peak_dvdt_ind] + peak_dvdt_ind += window_start + + iter_count = 0 + max_iter = 500 + while np.abs(peak_dvdt) > dvdt_thresh_envelope[peak_dvdt_ind - start_index]: + # found a peak to reject + iter_count += 1 + if iter_count > max_iter: + break + # determine extent of transient + + # find start + search_start = min(deflect_index, peak_dvdt_ind) + transient_start_index = np.flatnonzero(np.abs(dvdt[search_start:start_index - 1:-1]) < dvdt_thresh_envelope[search_start - start_index::-1] / 5)[0] + transient_start_index = search_start - transient_start_index + + transient_base_avg = np.mean(v[transient_start_index - window_width * 2:transient_start_index]) + transient_base_range = 3 * np.std(v[transient_start_index - window_width:transient_start_index]) + + # find end + search_start = max(deflect_index, peak_dvdt_ind) + baseline_return = np.flatnonzero(np.abs(v[search_start:] - transient_base_avg) < transient_base_range) + if len(baseline_return) > 0: + transient_end_index = baseline_return[0] + transient_end_index += search_start + else: + transient_end_index = len(t) - 1 + + # blank out the transient + + temp_v[transient_start_index:transient_end_index + 1] = np.nan + dvdt[transient_start_index:transient_end_index + 1] = np.nan + + # find a new peak + v_window = temp_v[start_index:end_index] + deflect_index = nan_deflect_func(v_window) + start_index + + # check the dv/dt around the new peak + window_start = max(deflect_index - window_width, start_index) + window_end = min(deflect_index + window_width, end_index) + dvdt_window = dvdt[window_start:window_end] + peak_dvdt_ind = np.nanargmax(np.abs(dvdt_window)) + peak_dvdt = dvdt_window[peak_dvdt_ind] + peak_dvdt_ind += window_start + return v[deflect_index], deflect_index @@ -94,8 +171,17 @@ def time_constant(t, v, i, start, end, max_fit_end=None, ------- tau : membrane time constant in seconds """ - # Assumes this is being done on a hyperpolarizing step - v_peak, peak_index = voltage_deflection(t, v, i, start, end, "min") + + # Check if this is being done on a hyperpolarizing step + check_index = tsu.find_time_index(t, end - baseline_interval) + stim_amp = i[check_index] + if stim_amp < 0: + reject_transients = True + else: + reject_transients = False + + v_peak, peak_index = voltage_deflection( + t, v, i, start, end, "min", reject_transients=reject_transients) if max_fit_end is not None: max_peak_index = tsu.find_time_index(t, max_fit_end) peak_index = min(max_peak_index, peak_index) @@ -139,7 +225,18 @@ def sag(t, v, i, start, end, peak_width=0.005, baseline_interval=0.03): ------- sag : fraction that membrane potential relaxes back to baseline """ - v_peak, peak_index = voltage_deflection(t, v, i, start, end, "min") + + # Check if actually hyperpolarizing (otherwise don't use reject_transients option + # for voltage deflection calculation, since there may be APs) + check_index = tsu.find_time_index(t, end - baseline_interval) + stim_amp = i[check_index] + if stim_amp < 0: + reject_transients = True + else: + reject_transients = False + + v_peak, peak_index = voltage_deflection(t, v, i, + start, end - baseline_interval, "min", reject_transients=reject_transients) v_peak_avg = tsu.average_voltage(v, t, start=t[peak_index] - peak_width / 2., end=t[peak_index] + peak_width / 2.) v_baseline = baseline_voltage(t, v, start, baseline_interval=baseline_interval) @@ -155,22 +252,31 @@ 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) + + if type(start) is list: + starts = start + ends = end + else: + starts = [start] * len(t_set) + ends = [end] * len(t_set) + + for t, i, v, start, end in zip(t_set, i_set, v_set, starts, ends): + v_peak, min_index = voltage_deflection( + t, v, i, start, end, 'min', reject_transients=True) + v_baseline = baseline_voltage(t, v, start, baseline_interval=baseline_interval) + v_vals.append(v_peak - v_baseline) i_vals.append(i[min_index]) v = np.array(v_vals) i = np.array(i_vals) if len(v) == 1: - # If there's just one sweep, we'll have to use its own baseline to estimate - # the input resistance - v = np.append(v, baseline_voltage(t_set[0], v_set[0], start, baseline_interval=baseline_interval)) + # If there's just one sweep, have to add zero deflection, zero stimulus point + v = np.append(v, 0.) i = np.append(i, 0.) A = np.vstack([i, np.ones_like(i)]).T - m, c = np.linalg.lstsq(A, v,rcond=None)[0] + m, c = np.linalg.lstsq(A, v, rcond=None)[0] return m * 1e3 diff --git a/ipfx/time_series_utils.py b/ipfx/time_series_utils.py index bb290ae1..0a23a6d5 100755 --- a/ipfx/time_series_utils.py +++ b/ipfx/time_series_utils.py @@ -1,4 +1,5 @@ import sys +from functools import lru_cache import numpy as np import scipy.signal as signal @@ -10,7 +11,7 @@ def find_time_index(t, t_0): Parameters ---------- - t : time array + t : time array (assumed monotonically non-decreasing) t_0 : time point to find an index Returns @@ -19,10 +20,26 @@ def find_time_index(t, t_0): """ assert t[0] <= t_0 <= t[-1], "Given time ({:f}) is outside of time range ({:f}, {:f})".format(t_0, t[0], t[-1]) - idx = np.argmin(abs(t - t_0)) + # t is sorted (guaranteed by the assert above), so a binary search finds + # the closest sample in O(log n) rather than scanning the whole array. + idx = int(np.searchsorted(t, t_0)) + if idx > 0 and (idx == len(t) or abs(t[idx - 1] - t_0) <= abs(t[idx] - t_0)): + # prefer the lower index on ties to match np.argmin's behavior + idx -= 1 return idx +@lru_cache(maxsize=None) +def _bessel_lowpass_coeffs(order, filt_coeff): + """Design (and cache) a low-pass Bessel filter. + + Coefficients depend only on the order and the normalized cutoff, which are + fixed across sweeps sharing a sampling rate, so caching avoids redesigning + the filter on every call. + """ + return signal.bessel(order, filt_coeff, "low") + + def calculate_dvdt(v, t, filter=None): """Low-pass filters (if requested) and differentiates voltage by time. @@ -43,7 +60,7 @@ def calculate_dvdt(v, t, filter=None): 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") + b, a = _bessel_lowpass_coeffs(4, filt_coeff) v_filt = signal.filtfilt(b, a, v, axis=0) dv = np.diff(v_filt) else: diff --git a/tests/data/feature_vector/features_T301.csv b/tests/data/feature_vector/features_T301.csv index 3bf53163..a38bb9c7 100644 --- a/tests/data/feature_vector/features_T301.csv +++ b/tests/data/feature_vector/features_T301.csv @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:d8e63f82a504a72780b35785947772552c7a39705f01375e7e82d62454a81df2 -size 7738 +oid sha256:1f7f174bee6eb7139d7eb329688ef600b05ac13ba37370b8f7a8f2d60720057c +size 7749 diff --git a/tests/data/feature_vector/fv_inst_freq_TEMP.npy b/tests/data/feature_vector/fv_inst_freq_TEMP.npy index 37939215..0c0d3b7f 100644 --- a/tests/data/feature_vector/fv_inst_freq_TEMP.npy +++ b/tests/data/feature_vector/fv_inst_freq_TEMP.npy @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:070c61ec83138f6369bade59da292c0ec2da2e313a0f118290f4308cb5d48c13 +oid sha256:87d63d408a5ed18e719941839992bc81ed03e6f415f3cac70c67bcf6c654c363 size 4928 diff --git a/tests/data/feature_vector/fv_isi_shape_TEMP.npy b/tests/data/feature_vector/fv_isi_shape_TEMP.npy index d03e6241..c98a06c8 100644 --- a/tests/data/feature_vector/fv_isi_shape_TEMP.npy +++ b/tests/data/feature_vector/fv_isi_shape_TEMP.npy @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:3ec22cc42846ecf06a2b761b311fe72abbffed07d137f5e7d452789b961903a6 -size 928 +oid sha256:e3e92a0e6a8b6c39b29fbec5a36dbe9a85c207b4c86ceef6853fb8d7bf9353e9 +size 1728 diff --git a/tests/data/feature_vector/fv_spiking_upstroke_downstroke_ratio_TEMP.npy b/tests/data/feature_vector/fv_spiking_upstroke_downstroke_ratio_TEMP.npy index f88b012f..77e58b13 100644 --- a/tests/data/feature_vector/fv_spiking_upstroke_downstroke_ratio_TEMP.npy +++ b/tests/data/feature_vector/fv_spiking_upstroke_downstroke_ratio_TEMP.npy @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:7716d3f24c2ebcf5c2027f353853abc4d1ef3c6cd8cf62806bce4a5461b21549 +oid sha256:196111c5ac4fd45d0da243b5b595f961e3560fec965dab6e68d56a0a5b414493 size 4928 diff --git a/tests/data/feature_vector/fv_subthresh_norm_TEMP.npy b/tests/data/feature_vector/fv_subthresh_norm_TEMP.npy index d3c20cff..746b5271 100644 --- a/tests/data/feature_vector/fv_subthresh_norm_TEMP.npy +++ b/tests/data/feature_vector/fv_subthresh_norm_TEMP.npy @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:84f129e00341e7f4d1353b2ba87b6047245599c2688000d81e9c039d208a6d48 +oid sha256:149b23e025be656b4abab2a5dc16b3fbf24026e2b42848d084daa0e8ba87a46a size 1248 diff --git a/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/pipeline_output.json b/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/pipeline_output.json index 0c2df9cc..22459fbd 100644 --- a/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/pipeline_output.json +++ b/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/pipeline_output.json @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:368bcc976f11437ccea47c7533b3b0669ce7471d3c7a0f51f1b051d74f9c6524 -size 380671 +oid sha256:5c420091696210ce2a19161a36d51af7c7ac7b148c11d984e5d07a19c92a72ac +size 382792 diff --git a/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/stimulus_ontology.json b/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/stimulus_ontology.json index 650c318a..b158829b 100644 --- a/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/stimulus_ontology.json +++ b/tests/data/specimens/Ctgf-T2A-dgCre;Ai14-495723.05.02.01/stimulus_ontology.json @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:62ce5ed148c99b343e672d06364c433ec29046c58ea5602cbac5679650bb3af9 -size 24526 +oid sha256:32f47b15b6f3c9a829473d57631dd28e0316ad746f183b83b8592a7488e62c13 +size 27660 diff --git a/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/pipeline_output.json b/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/pipeline_output.json index 3b5a965e..c7bd65a4 100644 --- a/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/pipeline_output.json +++ b/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/pipeline_output.json @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:99982b633670a238baa4758c9026fb2091321800c473e1f3590ce3d108ff15b6 -size 1652494 +oid sha256:c28cf799bad67ed76960d58da670b0f6126cd71e6a01f2111f70d0f1da0a553d +size 1676647 diff --git a/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/stimulus_ontology.json b/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/stimulus_ontology.json index 650c318a..b158829b 100644 --- a/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/stimulus_ontology.json +++ b/tests/data/specimens/Pvalb-IRES-Cre;Ai14(IVSCC)-165172.05.02/stimulus_ontology.json @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:62ce5ed148c99b343e672d06364c433ec29046c58ea5602cbac5679650bb3af9 -size 24526 +oid sha256:32f47b15b6f3c9a829473d57631dd28e0316ad746f183b83b8592a7488e62c13 +size 27660 diff --git a/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/pipeline_output.json b/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/pipeline_output.json index 2fceb01d..d20cce3a 100644 --- a/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/pipeline_output.json +++ b/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/pipeline_output.json @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:022ab74919eff869aa9893c241d736909ef4cdda762fd9e7feb12c729fc630fa -size 441882 +oid sha256:7414055b5bb3a99ccdbc194c817ae5b5dae0709d2470c071ea2ec6db159194ea +size 444272 diff --git a/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/stimulus_ontology.json b/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/stimulus_ontology.json index 650c318a..b158829b 100644 --- a/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/stimulus_ontology.json +++ b/tests/data/specimens/Vip-IRES-Cre;Ai14-331294.04.01.01/stimulus_ontology.json @@ -1,3 +1,3 @@ version https://git-lfs.github.com/spec/v1 -oid sha256:62ce5ed148c99b343e672d06364c433ec29046c58ea5602cbac5679650bb3af9 -size 24526 +oid sha256:32f47b15b6f3c9a829473d57631dd28e0316ad746f183b83b8592a7488e62c13 +size 27660 diff --git a/tests/test_feature_vector.py b/tests/test_feature_vector.py index b922b00e..9e889c9d 100755 --- a/tests/test_feature_vector.py +++ b/tests/test_feature_vector.py @@ -5,6 +5,7 @@ import ipfx.feature_vectors as fv from ipfx.stimulus import StimulusOntology from ipfx.sweep import Sweep, SweepSet +from ipfx.script_utils import StimulusTiming import ipfx.json_utilities as ju import pytest @@ -239,21 +240,32 @@ def test_identify_sweep_for_isi_shape(): end = 1.5 start = 0.5 + ndata = len(sweeps["data"][0]) + nsweeps = len(sweeps["index"]) + + class FakeSweep: + def __init__(self, sweep_number, v): + self.sweep_number = sweep_number + self.v = v + class SomeSweeps: @property def sweeps(self): - ndata = len(sweeps["data"][0]) - nsweeps = len(sweeps["index"]) - return np.arange(nsweeps * ndata).reshape(nsweeps, ndata) + arr = np.arange(nsweeps * ndata).reshape(nsweeps, ndata) + return [FakeSweep(i, arr[i]) for i in range(nsweeps)] + + stim_timing_dict = { + i: StimulusTiming(start, end, end - start) for i in range(nsweeps) + } isi_sweep, isi_sweep_spike_info = fv.identify_sweep_for_isi_shape( SomeSweeps(), {"sweeps": pd.DataFrame(**sweeps), "spikes_set": {5: "foo"}}, - end - start, + stim_timing_dict, ) assert isi_sweep_spike_info == "foo" - assert np.allclose(isi_sweep, np.arange(55, 66)) + assert np.allclose(isi_sweep.v, np.arange(55, 66)) def test_isi_shape(): @@ -262,9 +274,12 @@ def test_isi_shape(): "fast_trough_index": [0, 10, -10000], "threshold_index": [-10000, 10, 20], "threshold_v": [1, 2, -10000], + "clipped": [False, False, False], } class Sweep: + sweep_number = 0 + @property def v(self): return np.arange(20) @@ -276,7 +291,7 @@ def t(self): obtained = fv.isi_shape( Sweep(), pd.DataFrame(sweep_spike_info), - 50, + {0: StimulusTiming(0, 50, 50)}, n_points=10 ) assert np.allclose(np.arange(3.5, 13.5, 1.0), obtained) @@ -309,6 +324,8 @@ def sweeps(self): def test_step_subthreshold(): class Sweep: + sweep_number = 0 + @property def v(self): return np.arange(10) @@ -322,7 +339,9 @@ def t(self): } obtained = fv.step_subthreshold( - subthresh_hyperpol_dict, [-30], 4, 6, subsample_interval=1) + subthresh_hyperpol_dict, [-30], {0: StimulusTiming(4, 6, 2)}, + extend_duration_before=0, extend_duration_after=0, + subsample_interval=1, remove_transients=False) assert np.allclose(obtained, [4, 5]) @@ -340,18 +359,23 @@ def test_step_subthreshold_interpolation(): } sampling_rate = 1 clamp_mode = "CurrentClamp" - for a in test_amps: + for idx, a in enumerate(test_amps): v = np.hstack([np.zeros(2), np.ones(2) * a, np.zeros(2)]) - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, + sweep_number=idx, epochs=epochs) test_sweep_list.append(test_sweep) amp_sweep_dict = dict(zip(test_amps, test_sweep_list)) + stim_timing_dict = { + swp.sweep_number: StimulusTiming(2, 4, 2) for swp in test_sweep_list + } output = fv.step_subthreshold( amp_sweep_dict, target_amps, - start=2, - end=4, - extend_duration=1, + stim_timing_dict, + extend_duration_before=1, + extend_duration_after=1, subsample_interval=1, + remove_transients=False, ) assert np.all(output[1:3] == -90) assert np.array_equal(output[4:8], test_sweep_list[0].v[1:-1]) @@ -373,7 +397,8 @@ def test_subthresh_norm_normalization(): } sampling_rate = 1 clamp_mode = "CurrentClamp" - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, sweep_number=0, + epochs=epochs) base = v[0] deflect_v = np.min(v) @@ -383,25 +408,41 @@ def test_subthresh_norm_normalization(): output = fv.subthresh_norm( amp_sweep_dict, deflect_dict, - start=t[0], - end=t[-1], + {0: StimulusTiming(t[0], t[-1], t[-1] - t[0])}, target_amp=-10, - extend_duration=0, + extend_duration_before=0, + extend_duration_after=0, subsample_interval=1, + remove_transients=False, ) assert np.isclose(output[0], 0) assert np.isclose(np.min(output), -1) def test_subthresh_depol_norm_bad_steady_state_interval(): + class Sweep: + sweep_number = 0 + + @property + def v(self): + return np.arange(10) + + @property + def t(self): + return np.arange(10) + + amp_sweep_dict = {50: Sweep()} + deflect_dict = {50: (1, 2)} with pytest.raises(ValueError): fv.subthresh_depol_norm( - {}, {}, start=0, end=1, steady_state_interval=2 + amp_sweep_dict, deflect_dict, {0: StimulusTiming(0, 1, 1)}, + steady_state_interval=2 ) def test_subthresh_depol_norm_empty_result(): - output = fv.subthresh_depol_norm({}, {}, start=1.02, end=2.02) + output = fv.subthresh_depol_norm( + {}, {}, {0: StimulusTiming(1.02, 2.02, 1.0)}) assert len(output) == 140 assert np.all(np.isnan(output)) @@ -419,7 +460,8 @@ def test_subthresh_depol_norm_normalization(): } sampling_rate = 1 clamp_mode = "CurrentClamp" - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, sweep_number=0, + epochs=epochs) base = v[0] deflect_v = np.max(v) @@ -430,8 +472,7 @@ def test_subthresh_depol_norm_normalization(): output = fv.subthresh_depol_norm( amp_sweep_dict, deflect_dict, - start=t[0], - end=t[-1], + {0: StimulusTiming(t[0], t[-1], t[-1] - t[0])}, steady_state_interval=1, subsample_interval=1, extend_duration=0, @@ -473,6 +514,8 @@ def sweeps(self): def test_subthresh_depol_norm(): class Sweep: + sweep_number = 0 + @property def v(self): return np.arange(10) @@ -485,7 +528,8 @@ def t(self): deflect_dict = {50: (1, 2)} obtained = fv.subthresh_depol_norm( - amp_sweep_dict, deflect_dict, 4, 7, subsample_interval=1, + amp_sweep_dict, deflect_dict, {0: StimulusTiming(4, 7, 3)}, + extend_duration=0, subsample_interval=1, steady_state_interval=2 ) assert np.allclose([2/3, 8/9, 10/9], obtained) @@ -564,7 +608,9 @@ def test_psth_sparse_firing(): width = 50 # All spikes are in own bins - output = fv.psth_vector([spike_info], start=start, end=end, width=width) + output = fv.psth_vector( + [spike_info], [StimulusTiming(start, end, end - start)], + width=width, duration=end - start) assert np.sum(output > 0) == len(test_spike_times) assert np.max(output) == 1 / (width * 0.001) @@ -577,7 +623,9 @@ def test_psth_compressed_firing(): width = 50 # All spikes are within one bin - output = fv.psth_vector([spike_info], start=start, end=end, width=width) + output = fv.psth_vector( + [spike_info], [StimulusTiming(start, end, end - start)], + width=width, duration=end - start) assert np.sum(output > 0) == 1 assert np.max(output) == len(test_spike_times) / (width * 0.001) @@ -591,7 +639,9 @@ def test_psth_number_of_spikes(): test_spike_times = np.random.random(n_spikes) * (end - start) + start spike_info = pd.DataFrame({"threshold_t": test_spike_times}) - output = fv.psth_vector([spike_info], start=start, end=end, width=width) + output = fv.psth_vector( + [spike_info], [StimulusTiming(start, end, end - start)], + width=width, duration=end - start) assert np.isclose(output.mean(), n_spikes) @@ -608,7 +658,9 @@ def test_psth_between_sweep_interpolation(): n_bins = int((end - start) / (width * 0.001)) si_list = [None, available_list[0], None, available_list[1], None] - output = fv.psth_vector(si_list, start=start, end=end, width=width) + stim_timing_list = [StimulusTiming(start, end, end - start)] * len(si_list) + output = fv.psth_vector( + si_list, stim_timing_list, width=width, duration=end - start) assert np.array_equal(output[:n_bins], output[n_bins : 2 * n_bins]) assert np.array_equal( @@ -633,9 +685,14 @@ def test_psth_duration_rounding(): spike_info = pd.DataFrame({"threshold_t": test_spike_times}) width = 50 + duration = 1.0 - output_a = fv.psth_vector([spike_info], start=start_a, end=end_a, width=width) - output_b = fv.psth_vector([spike_info], start=start_b, end=end_b, width=width) + output_a = fv.psth_vector( + [spike_info], [StimulusTiming(start_a, end_a, end_a - start_a)], + width=width, duration=duration) + output_b = fv.psth_vector( + [spike_info], [StimulusTiming(start_b, end_b, end_b - start_b)], + width=width, duration=duration) assert output_a.shape == output_b.shape @@ -647,7 +704,8 @@ def test_inst_freq_one_spike(): end = 1 width = 20 output = fv.inst_freq_vector( - [spike_info], start=start, end=end, width=width + [spike_info], [StimulusTiming(start, end, end - start)], + width=width, duration=end - start ) assert np.all(output >= 1 / (end - start)) @@ -658,8 +716,17 @@ def test_inst_freq_initial_freq(): start = 0 end = 1 width = 20 - output = fv.inst_freq_vector([spike_info], start=start, end=end, width=width) - assert output[0] == 1.0 / (test_spike_times[0] - start) + output = fv.inst_freq_vector( + [spike_info], [StimulusTiming(start, end, end - start)], + width=width, duration=end - start) + + # The instantaneous rate is now only estimated at spike times (the rate at + # each spike being the inverse of the average of its adjacent ISIs), and + # bins before the first spike extrapolate to that first estimate. + first_isi = test_spike_times[0] - start + second_isi = test_spike_times[1] - test_spike_times[0] + first_spike_rate = 1.0 / ((first_isi + second_isi) / 2) + assert output[0] == first_spike_rate def test_inst_freq_between_sweep_interpolation(): @@ -675,7 +742,9 @@ def test_inst_freq_between_sweep_interpolation(): n_bins = int((end - start) / (width * 0.001)) si_list = [None, available_list[0], None, available_list[1], None] - output = fv.inst_freq_vector(si_list, start=start, end=end, width=width) + stim_timing_list = [StimulusTiming(start, end, end - start)] * len(si_list) + output = fv.inst_freq_vector( + si_list, stim_timing_list, width=width, duration=end - start) assert np.array_equal(output[:n_bins], output[n_bins : 2 * n_bins]) assert np.array_equal( @@ -697,9 +766,14 @@ def test_inst_freq_duration_rounding(): spike_info = pd.DataFrame({"threshold_t": test_spike_times}) width = 20 + duration = 1.0 - output_a = fv.inst_freq_vector([spike_info], start=start_a, end=end_a, width=width) - output_b = fv.inst_freq_vector([spike_info], start=start_b, end=end_b, width=width) + output_a = fv.inst_freq_vector( + [spike_info], [StimulusTiming(start_a, end_a, end_a - start_a)], + width=width, duration=duration) + output_b = fv.inst_freq_vector( + [spike_info], [StimulusTiming(start_b, end_b, end_b - start_b)], + width=width, duration=duration) assert output_a.shape == output_b.shape @@ -719,7 +793,8 @@ def test_spike_feature_within_sweep_interpolation(): end = 1 width = 20 output = fv.spike_feature_vector( - feature, [spike_info], start=start, end=end, width=width + feature, [spike_info], [StimulusTiming(start, end, end - start)], + width=width, duration=end - start ) assert output[0] == test_feature_values[0] assert output[len(output) // 2] > test_feature_values[0] @@ -747,8 +822,9 @@ def test_spike_feature_between_sweep_interpolation(): width = 20 n_bins = int((end - start) / (width * 0.001)) si_list = [None, available_list[0], None, available_list[1], None] + stim_timing_list = [StimulusTiming(start, end, end - start)] * len(si_list) output = fv.spike_feature_vector( - feature, si_list, start=start, end=end, width=width + feature, si_list, stim_timing_list, width=width, duration=end - start ) assert np.all(output[:n_bins] == test_feature_values[0]) assert np.all(output[n_bins : 2 * n_bins] == test_feature_values[0]) @@ -773,11 +849,14 @@ def test_spike_feature_duration_rounding(): }) width = 20 + duration = 1.0 output_a = fv.spike_feature_vector( - feature, [spike_info], start=start_a, end=end_a, width=width) + feature, [spike_info], [StimulusTiming(start_a, end_a, end_a - start_a)], + width=width, duration=duration) output_b = fv.spike_feature_vector( - feature, [spike_info], start=start_b, end=end_b, width=width) + feature, [spike_info], [StimulusTiming(start_b, end_b, end_b - start_b)], + width=width, duration=duration) assert output_a.shape == output_b.shape @@ -981,14 +1060,18 @@ def test_identify_isi_shape_min_spike(): sampling_rate = 1 clamp_mode = "CurrentClamp" sweep_list = [] - for a in test_input_amplitudes: + for idx, a in enumerate(test_input_amplitudes): v = np.random.randn(n_points) - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, + sweep_number=idx, epochs=epochs) sweep_list.append(test_sweep) test_sweep_set = SweepSet(sweep_list) + stim_timing_dict = { + swp.sweep_number: StimulusTiming(0, 1, 1) for swp in sweep_list + } selected_sweep, _ = fv.identify_sweep_for_isi_shape( - test_sweep_set, test_features, duration=1, min_spike=min_spike + test_sweep_set, test_features, stim_timing_dict, min_spike=min_spike ) assert np.array_equal(selected_sweep.v, sweep_list[2].v) @@ -1020,14 +1103,18 @@ def test_identify_isi_shape_largest_below_min_spike(): sampling_rate = 1 clamp_mode = "CurrentClamp" sweep_list = [] - for a in test_input_amplitudes: + for idx, a in enumerate(test_input_amplitudes): v = np.random.randn(n_points) - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, + sweep_number=idx, epochs=epochs) sweep_list.append(test_sweep) test_sweep_set = SweepSet(sweep_list) + stim_timing_dict = { + swp.sweep_number: StimulusTiming(0, 1, 1) for swp in sweep_list + } selected_sweep, _ = fv.identify_sweep_for_isi_shape( - test_sweep_set, test_features, duration=1, min_spike=min_spike + test_sweep_set, test_features, stim_timing_dict, min_spike=min_spike ) assert np.array_equal(selected_sweep.v, sweep_list[-1].v) @@ -1059,14 +1146,18 @@ def test_identify_isi_shape_one_spike(): sampling_rate = 1 clamp_mode = "CurrentClamp" sweep_list = [] - for a in test_input_amplitudes: + for idx, a in enumerate(test_input_amplitudes): v = np.random.randn(n_points) - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, + sweep_number=idx, epochs=epochs) sweep_list.append(test_sweep) test_sweep_set = SweepSet(sweep_list) + stim_timing_dict = { + swp.sweep_number: StimulusTiming(0, 1, 1) for swp in sweep_list + } selected_sweep, _ = fv.identify_sweep_for_isi_shape( - test_sweep_set, test_features, duration=1, min_spike=min_spike + test_sweep_set, test_features, stim_timing_dict, min_spike=min_spike ) assert np.array_equal(selected_sweep.v, sweep_list[1].v) @@ -1140,7 +1231,8 @@ def test_isi_shape_aligned(): } sampling_rate = 1 clamp_mode = "CurrentClamp" - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, sweep_number=0, + epochs=epochs) end = t[-100] test_threshold_index = np.array([100, 220, 340]) @@ -1155,12 +1247,14 @@ def test_isi_shape_aligned(): "fast_trough_index": test_fast_trough_index, "threshold_v": test_threshold_v, "fast_trough_t": test_fast_trough_index, + "clipped": np.zeros(len(test_threshold_index), dtype=bool), } ) n_points = 100 isi_norm = fv.isi_shape( - test_sweep, test_spike_info, end, n_points=n_points + test_sweep, test_spike_info, {0: StimulusTiming(0, end, end)}, + n_points=n_points ) assert len(isi_norm) == n_points assert isi_norm[0] == np.mean( @@ -1182,7 +1276,8 @@ def test_isi_shape_skip_short(): } sampling_rate = 1 clamp_mode = "CurrentClamp" - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, sweep_number=0, + epochs=epochs) end = t[-100] test_subsample = 3 @@ -1198,12 +1293,14 @@ def test_isi_shape_skip_short(): "fast_trough_index": test_fast_trough_index, "threshold_v": test_threshold_v, "fast_trough_t": test_fast_trough_index, + "clipped": np.zeros(len(test_threshold_index), dtype=bool), } ) n_points = 100 isi_norm = fv.isi_shape( - test_sweep, test_spike_info, end, n_points=n_points + test_sweep, test_spike_info, {0: StimulusTiming(0, end, end)}, + n_points=n_points ) assert len(isi_norm) == n_points @@ -1232,7 +1329,8 @@ def test_isi_shape_one_spike(): } sampling_rate = 1 clamp_mode = "CurrentClamp" - test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, epochs=epochs) + test_sweep = Sweep(t, v, i, clamp_mode, sampling_rate, sweep_number=0, + epochs=epochs) end = t[-100] test_threshold_index = [80] @@ -1244,11 +1342,13 @@ def test_isi_shape_one_spike(): "fast_trough_index": test_fast_trough_index, "threshold_v": test_threshold_v, "fast_trough_t": test_fast_trough_index, + "clipped": [False], }) n_points = 100 isi_norm = fv.isi_shape( - test_sweep, test_spike_info, end, n_points=n_points, + test_sweep, test_spike_info, {0: StimulusTiming(0, end, end)}, + n_points=n_points, steady_state_interval=10, single_max_duration=500 ) assert len(isi_norm) == n_points diff --git a/tests/test_lims_queries.py b/tests/test_lims_queries.py index a15fe053..e3d95718 100644 --- a/tests/test_lims_queries.py +++ b/tests/test_lims_queries.py @@ -15,7 +15,7 @@ def test_get_nwb_path_from_lims(): ephys_roi_result = 500844779 result = lq.get_nwb_path_from_lims(ephys_roi_result) - assert result == "/allen/programs/celltypes/production/mousecelltypes/prod589/Ephys_Roi_Result_500844779/500844779.nwb" + assert result == "/allen/programs/celltypes/production/mousecelltypes/prod589/Ephys_Roi_Result_500844779/nwb2_Vip-IRES-Cre;Ai14(IVSCC)-226110.03.01.nwb" @pytest.mark.requires_lims def test_get_igorh5_path_from_lims(): diff --git a/tests/test_qc_features.py b/tests/test_qc_features.py index 302c7598..e4f0337c 100644 --- a/tests/test_qc_features.py +++ b/tests/test_qc_features.py @@ -22,33 +22,66 @@ def test_measure_electrode_0(): def test_measure_seal(): - i = np.array([0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0]) - v = np.array([0, 1, 1, 0, 0, 0, 0, 1, 1, 1, 0]) - t = np.arange(len(v)) * 1E-3 - b = qcf.measure_seal(v, i, t) - assert np.allclose([b], [1.0]) + # measure_seal now averages a series of square voltage test pulses (skipping + # the first) and fits the capacitive transient to recover the steady-state + # resistance. Build a cell-attached-style recording: v is a clean square + # command (mV), curr is a step with a decaying capacitive transient (pA). + dt = 1e-5 + pulse_dur_pts = 500 # 5 ms + gap_pts = 400 # 4 ms between pulses + lead_pts = 400 + n_pulses = 4 # first is treated as the test pulse and skipped + + delta_v_mV = 5.0 + r_seal = 1e9 # 1 GOhm + i_ss_pA = (delta_v_mV * 1e-3 / r_seal) * 1e12 # steady current, pA + peak_extra_pA = 200.0 + tau = 5e-5 + + total = lead_pts + n_pulses * (pulse_dur_pts + gap_pts) + v = np.zeros(total) + curr = np.zeros(total) + t_rel = np.arange(pulse_dur_pts) * dt + idx = lead_pts + for _ in range(n_pulses): + up = idx + down = idx + pulse_dur_pts + v[up:down] = delta_v_mV + curr[up:down] = i_ss_pA + peak_extra_pA * np.exp(-t_rel / tau) + idx = down + gap_pts + t = dt * np.arange(total) + + b = qcf.measure_seal(v, curr, t) + assert np.allclose([b], [1.0], rtol=1e-3) def test_measure_input_resistance(): - - ir = 50.0 - dt = 1E-4 - time_range = [0, 0.5] - time = np.arange(time_range[0], time_range[1], dt) - current = np.zeros(time.shape) - - pulse_intervals = [(0.1, 0.2), (0.24, 0.3), (0.31, 0.32)] - current_magnitudes = [1, 2, 3] - - for pulse_interval, current_magnitude in zip(pulse_intervals, current_magnitudes): - - ix = np.where((time > pulse_interval[0]) & (time < pulse_interval[1])) - current[ix] = current_magnitude - - voltage = ir*current - - ir_tested = qcf.get_r_from_stable_pulse_response(voltage, current, time) - assert np.isclose(ir_tested, ir) + # get_r_from_stable_pulse_response_fit now operates on a single averaged + # pulse (avg_v in V, avg_i in A) plus the relative up/down indices, and fits + # the capacitive transient to estimate the steady-state resistance. + dt = 1e-5 + n = 1000 + up_ind = 200 + down_ind = 800 + + delta_v = 5e-3 # V + r_expected = 50e6 # Ohm (50 MOhm) + i_ss = delta_v / r_expected # steady-state current, A + peak_extra = 200e-12 + tau = 5e-5 + + avg_v = np.zeros(n) + avg_v[up_ind:down_ind] = delta_v + + avg_i = np.zeros(n) + t_rel = np.arange(down_ind - up_ind) * dt + avg_i[up_ind:down_ind] = i_ss + peak_extra * np.exp(-t_rel / tau) + + t = dt * np.arange(n) + + r = qcf.get_r_from_stable_pulse_response_fit( + avg_v, avg_i, t, up_ind, down_ind, post_transient_shift_ms=1.0) + assert np.isclose(r, r_expected, rtol=1e-3) def test_get_square_pulse_idx(): diff --git a/tests/test_subthresh_features.py b/tests/test_subthresh_features.py index 791d2248..2ded28d1 100644 --- a/tests/test_subthresh_features.py +++ b/tests/test_subthresh_features.py @@ -3,13 +3,26 @@ def test_input_resistance(): - t = np.arange(0, 1.0, 5e-6) - v1 = np.ones_like(t) * -5. - v2 = np.ones_like(t) * -10. - i1 = np.ones_like(t) * -50. - i2 = np.ones_like(t) * -100. - - ri = subf.input_resistance([t, t], [i1, i2], [v1, v2], 0, t[-1]) + # input_resistance now measures the deflection relative to a pre-stimulus + # baseline, so the sweeps need a real baseline period before `start` and a + # (gently-settling, to avoid transient rejection) step during the stimulus. + dt = 5e-6 + t = np.arange(0, 1.0, dt) + start = 0.3 + end = 0.8 + ramp_dur = 0.1 + + def make_sweep(deflection, current): + ramp = np.clip((t - start) / ramp_dur, 0.0, 1.0) + v = np.where(t < start, 0.0, deflection * ramp) + i = np.ones_like(t) * current + return v, i + + # 100 MOhm: 0.1 mV/pA -> -50 pA gives -5 mV, -100 pA gives -10 mV + v1, i1 = make_sweep(-5., -50.) + v2, i2 = make_sweep(-10., -100.) + + ri = subf.input_resistance([t, t], [i1, i2], [v1, v2], start, end) assert np.allclose(ri, 100.)