From 527f6b68fafda469300f4565e83c598a4e03bf6d Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Thu, 16 Oct 2025 15:08:25 +0100 Subject: [PATCH 01/13] Slightly tidying hergQC.py --- pcpostprocess/hergQC.py | 29 +++++++++++++---------------- 1 file changed, 13 insertions(+), 16 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index e248a3ca..28a5b370 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -13,14 +13,14 @@ class QCDict(UserDict): """ Stores the results from QC checks. - Each entry is a label -> [(bool, value),...] mapping. + Each entry is a ``label -> [(bool, value),...]`` mapping. The bool in each tuple indicates whether the QC passed, and the value is the result that was checked (e.g. the SNR value). The list can contain multiple tuples if the QC checks multiple values, as QC1 does, or if the checks are run multiple times e.g. once per sweep. """ - labels = [ + labels = ( "qc1.rseal", "qc1.cm", "qc1.rseries", @@ -37,7 +37,7 @@ class QCDict(UserDict): "qc6.subtracted", "qc6.1.subtracted", "qc6.2.subtracted", - ] + ) def __init__(self): super().__init__([(label, [(False, None)]) for label in QCDict.labels]) @@ -49,15 +49,15 @@ def __setitem__(self, key, value): def qc_passed(self, label): """Return whether a single QC passed.""" - return all([x for x, _ in self[label]]) + return all(x for x, _ in self[label]) def passed_list(self): """Return a list of booleans indicating whether each QC passed.""" - return [self.qc_passed(label) for label in QCDict.labels] + return [all(x[0] for x in tuples) for tuples in self.values()] def all_passed(self): """Return whether all QC passed.""" - return all(self.passed_list()) + return all(all(x[0] for x in tuples) for tuples in self.values()) class hERGQC: @@ -68,19 +68,24 @@ def __init__(self, sampling_rate=5, plot_dir=None, voltage=np.array([]), self._plot_dir = plot_dir - self._n_qc = 16 + self._n_qc = len(QCDict.labels) self.removal_time = removal_time self.voltage = np.array(voltage) self.logger = logging.getLogger(__name__) - self.logger.setLevel(logging.DEBUG) + self._debug = True + if self._debug: + self.logger.setLevel(logging.DEBUG) + # https://github.com/CardiacModelling/pcpostprocess/issues/42 self.sampling_rate = sampling_rate # Define all thresholds + # TODO: These should be args? Or perhaps this is good so that this QC + # class can be extended # qc1 self.rsealc = [1e8, 1e12] # in Ohm, converted from [0.1, 1000] GOhm self.cmc = [1e-12, 1e-10] # in F, converted from [1, 100] pF @@ -122,8 +127,6 @@ def __init__(self, sampling_rate=5, plot_dir=None, voltage=np.array([]), self.qc6_1_win = self.qc6_1_win.astype(int) self.qc6_2_win = self.qc6_2_win.astype(int) - self._debug = True - @property def plot_dir(self): return self._plot_dir @@ -132,12 +135,6 @@ def plot_dir(self): def plot_dir(self, path): self._plot_dir = path - def set_trace(self, before, after, qc_vals_before, qc_vals_after, n_sweeps): - self._before = before - self._qc_vals_before = qc_vals_before - self._after = after - self._qc_vals_after = qc_vals_after - def set_debug(self, debug): self._debug = debug From e7de3a2a168b0b54afa53346f3f7e55d2aa0559f Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Fri, 17 Oct 2025 16:51:00 +0100 Subject: [PATCH 02/13] Removed unused qc2 and qc3 options --- pcpostprocess/hergQC.py | 49 ++++++++++++++--------------------------- 1 file changed, 16 insertions(+), 33 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index 28a5b370..0aec19d1 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -21,11 +21,11 @@ class QCDict(UserDict): """ labels = ( - "qc1.rseal", - "qc1.cm", - "qc1.rseries", - "qc2.raw", - "qc2.subtracted", + "qc1.rseal", # [before, after] + "qc1.cm", # [before, after] + "qc1.rseries", # [before, after] + "qc2.raw", # [sweep[0], sweep[1], ...] + "qc2.subtracted", # [sweep[0], sweep[1], ...] "qc3.raw", "qc3.E4031", "qc3.subtracted", @@ -301,18 +301,10 @@ def qc1(self, rseal, cm, rseries): return [(qc11, rseal), (qc12, cm), (qc13, rseries)] - def qc2(self, recording, method=3): + def qc2(self, recording): # Check SNR is good - if method == 1: - # Not sure if this is good... - snr = scipy.stats.signaltonoise(recording) - elif method == 2: - noise = np.std(recording[:NOISE_LEN]) - snr = (np.max(recording) - np.min(recording) - 2 * noise) / noise - elif method == 3: - noise = np.std(recording[:NOISE_LEN]) - snr = (np.std(recording) / noise) ** 2 - + noise = np.std(recording[:NOISE_LEN]) + snr = (np.std(recording) / noise) ** 2 if snr < self.snrc or not np.isfinite(snr): self.logger.debug(f"snr: {snr}") result = False @@ -321,24 +313,14 @@ def qc2(self, recording, method=3): return (result, snr) - def qc3(self, recording1, recording2, method=3): + def qc3(self, recording1, recording2): # Check 2 sweeps similar - if method == 1: - rmsdc = 2 # A/F * F - elif method == 2: - noise_1 = np.std(recording1[:NOISE_LEN]) - peak_1 = (np.max(recording1) - noise_1) - noise_2 = np.std(recording2[:NOISE_LEN]) - peak_2 = (np.max(recording2) - noise_2) - rmsdc = max(np.mean([peak_1, peak_2]) * 0.1, - np.mean([noise_1, noise_2]) * 5) - elif method == 3: - noise_1 = np.std(recording1[:NOISE_LEN]) - noise_2 = np.std(recording2[:NOISE_LEN]) - rmsd0_1 = np.sqrt(np.mean((recording1) ** 2)) - rmsd0_2 = np.sqrt(np.mean((recording2) ** 2)) - rmsdc = max(np.mean([rmsd0_1, rmsd0_2]) * self.rmsd0c, - np.mean([noise_1, noise_2]) * 6) + noise_1 = np.std(recording1[:NOISE_LEN]) + noise_2 = np.std(recording2[:NOISE_LEN]) + rmsd0_1 = np.sqrt(np.mean((recording1) ** 2)) + rmsd0_2 = np.sqrt(np.mean((recording2) ** 2)) + rmsdc = max(np.mean([rmsd0_1, rmsd0_2]) * self.rmsd0c, + np.mean([noise_1, noise_2]) * 6) rmsd = np.sqrt(np.mean((recording1 - recording2) ** 2)) if rmsd > rmsdc or not (np.isfinite(rmsd) and np.isfinite(rmsdc)): self.logger.debug(f"rmsd: {rmsd}, rmsdc: {rmsdc}") @@ -497,6 +479,7 @@ def filter_capacitive_spikes(self, current, times, voltage_step_times): if i_end == 0: break if len(current.shape) == 2: + # When would this happen? Is one dimension always 1 ? current[:, i_start:i_end] = 0 elif len(current.shape) == 1: current[i_start:i_end] = 0 From 1acc6a905f239e4d5cf3062086bb34780625b269 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Mon, 3 Nov 2025 14:13:54 +0000 Subject: [PATCH 03/13] Working on docstrings for hergQC module --- pcpostprocess/hergQC.py | 235 ++++++++++++++++++++-------------------- tests/test_herg_qc.py | 84 +++++++------- 2 files changed, 161 insertions(+), 158 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index 0aec19d1..4e213e47 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -4,9 +4,6 @@ import matplotlib.pyplot as plt import numpy as np -import scipy.stats - -NOISE_LEN = 200 class QCDict(UserDict): @@ -21,22 +18,22 @@ class QCDict(UserDict): """ labels = ( - "qc1.rseal", # [before, after] - "qc1.cm", # [before, after] - "qc1.rseries", # [before, after] - "qc2.raw", # [sweep[0], sweep[1], ...] - "qc2.subtracted", # [sweep[0], sweep[1], ...] - "qc3.raw", - "qc3.E4031", - "qc3.subtracted", - "qc4.rseal", - "qc4.cm", - "qc4.rseries", - "qc5.staircase", - "qc5.1.staircase", - "qc6.subtracted", - "qc6.1.subtracted", - "qc6.2.subtracted", + 'qc1.rseal', # [before, after] + 'qc1.cm', # [before, after] + 'qc1.rseries', # [before, after] + 'qc2.raw', # [sweep[0], sweep[1], ...] + 'qc2.subtracted', # [sweep[0], sweep[1], ...] + 'qc3.raw', + 'qc3.E4031', + 'qc3.subtracted', + 'qc4.rseal', + 'qc4.cm', + 'qc4.rseries', + 'qc5.staircase', + 'qc5.1.staircase', + 'qc6.subtracted', + 'qc6.1.subtracted', + 'qc6.2.subtracted', ) def __init__(self): @@ -44,7 +41,7 @@ def __init__(self): def __setitem__(self, key, value): if key not in QCDict.labels: - raise KeyError(f"Invalid QC key: {key}") + raise KeyError(f'Invalid QC key: {key}') super().__setitem__(key, value) def qc_passed(self, label): @@ -61,18 +58,26 @@ def all_passed(self): class hERGQC: + """ + Performs hERG staircase quality control similar to Lei 2019. - def __init__(self, sampling_rate=5, plot_dir=None, voltage=np.array([]), - n_sweeps=None, removal_time=5): - # TODO docstring - - self._plot_dir = plot_dir + @param voltage: A voltage trace, sampled at ``sampling_rate``. + @param sampling_rate: The number of samples per time unit. + @param removal_time: Number of time units to remove after each step in the protocol + @param noise_len: Number of initial samples during which the signal is flat, to use for noise estimates. + @param plot_dir: An optional directory to store plots in - self._n_qc = len(QCDict.labels) + """ + def __init__(self, voltage, sampling_rate=5, removal_time=5, noise_len=200, + plot_dir=None): + self.voltage = np.array(voltage) + self.sampling_rate = sampling_rate self.removal_time = removal_time + self.noise_len = int(noise_len) + self._plot_dir = plot_dir - self.voltage = np.array(voltage) + self._n_qc = len(QCDict.labels) self.logger = logging.getLogger(__name__) self._debug = True @@ -80,8 +85,6 @@ def __init__(self, sampling_rate=5, plot_dir=None, voltage=np.array([]), self.logger.setLevel(logging.DEBUG) # https://github.com/CardiacModelling/pcpostprocess/issues/42 - self.sampling_rate = sampling_rate - # Define all thresholds # TODO: These should be args? Or perhaps this is good so that this QC @@ -100,32 +103,32 @@ def __init__(self, sampling_rate=5, plot_dir=None, voltage=np.array([]), self.rseriessc = 0.5 # qc5 self.max_diffc = 0.75 - # self.qc5_win = [3275, 5750] # indices with hERG screening peak! + # self._qc5_win = [3275, 5750] # indices with hERG screening peak! # indices where hERG could peak (for different temperatures) - self.qc5_win = np.array([8550 + 400, 10950 + 400]) * sampling_rate + self._qc5_win = np.array([8550 + 400, 10950 + 400]) * sampling_rate # qc5_1 self.rmsd0_diffc = 0.5 # qc6 self.negative_tolc = -2 ''' # These are for `whole` (just staircase) protocol at 5kHz - self.qc6_win = [3000, 7000] # indices for first +40 mV - self.qc6_1_win = [35250, 37250] # indices for second +40 mV - self.qc6_2_win = [40250, 42250] # indices for third +40 mV + self._qc6_win = [3000, 7000] # indices for first +40 mV + self._qc6_1_win = [35250, 37250] # indices for second +40 mV + self._qc6_2_win = [40250, 42250] # indices for third +40 mV ''' # These are for `staircaseramp` protocol # Firstly, indices for 1st +40 mV - self.qc6_win = np.array([1000, 1800]) * sampling_rate + self._qc6_win = np.array([1000, 1800]) * sampling_rate # indices for 2nd +40 mV - self.qc6_1_win = np.array([7450, 7850]) * sampling_rate + self._qc6_1_win = np.array([7450, 7850]) * sampling_rate # indices for 3rd +40 mV - self.qc6_2_win = np.array([8450, 8850]) * sampling_rate + self._qc6_2_win = np.array([8450, 8850]) * sampling_rate # Ensure these indices are integers - self.qc5_win = self.qc5_win.astype(int) - self.qc6_win = self.qc6_win.astype(int) - self.qc6_1_win = self.qc6_1_win.astype(int) - self.qc6_2_win = self.qc6_2_win.astype(int) + self._qc5_win = self._qc5_win.astype(int) + self._qc6_win = self._qc6_win.astype(int) + self._qc6_1_win = self._qc6_1_win.astype(int) + self._qc6_2_win = self._qc6_2_win.astype(int) @property def plot_dir(self): @@ -138,28 +141,26 @@ def plot_dir(self, path): def set_debug(self, debug): self._debug = debug - def run_qc(self, voltage_steps, times, - before, after, qc_vals_before, + def run_qc(self, voltage_steps, times, before, after, qc_vals_before, qc_vals_after, n_sweeps=None): """Run each QC criteria on a single (before-trace, after-trace) pair. @param voltage_steps is a list of times at which there are discontinuities in Vcmd @param times is the array of observation times - @param before is the before-drug current trace - @param after is the post-drug current trace - @param qc_vals_before is an array of values for each pre-drug sweep where each row is (Rseal, Cm, Rseries) - @param qc_vals_after is an array of values for each post-drug sweep where each row is (Rseal, Cm, Rseries) - @n_sweeps is the number of sweeps to be considered + @param before are the before-drug current traces, in an array with sweeps on the first index and values on the second. + @param after is the post-drug current traces, in an array with sweeps on the first index and values on the second. + @param qc_vals_before is an array of values for each before-drug sweep where each row is (Rseal, Cm, Rseries) + @param qc_vals_after is an array of values for each after-drug sweep where each row is (Rseal, Cm, Rseries) + @n_sweeps is the number of sweeps to be run QC on. + + @returns A :class:`QCDict` with the test results. """ if not n_sweeps: n_sweeps = len(before) - before = np.array(before) - after = np.array(after) - - before = self.filter_capacitive_spikes(before, times, voltage_steps) - after = self.filter_capacitive_spikes(after, times, voltage_steps) + before = self._filter_capacitive_spikes(np.array(before), times, voltage_steps) + after = self._filter_capacitive_spikes(np.array(after), times, voltage_steps) QC = QCDict() @@ -169,8 +170,8 @@ def run_qc(self, voltage_steps, times, if (None in qc_vals_before) or (None in qc_vals_after): return QC - qc1_1 = self.qc1(*qc_vals_before) - qc1_2 = self.qc1(*qc_vals_after) + qc1_1 = self._qc1(*qc_vals_before) + qc1_2 = self._qc1(*qc_vals_after) QC['qc1.rseal'] = [qc1_1[0], qc1_2[0]] QC['qc1.cm'] = [qc1_1[1], qc1_2[1]] @@ -179,73 +180,73 @@ def run_qc(self, voltage_steps, times, QC['qc2.raw'] = [] QC['qc2.subtracted'] = [] for i in range(n_sweeps): - qc2_1 = self.qc2(before[i]) - qc2_2 = self.qc2(before[i] - after[i]) + qc2_1 = self._qc2(before[i]) + qc2_2 = self._qc2(before[i] - after[i]) QC['qc2.raw'].append(qc2_1) QC['qc2.subtracted'].append(qc2_2) - qc3_1 = self.qc3(before[0, :], before[1, :]) - qc3_2 = self.qc3(after[0, :], after[1, :]) - qc3_3 = self.qc3(before[0, :] - after[0, :], + qc3_1 = self._qc3(before[0, :], before[1, :]) + qc3_2 = self._qc3(after[0, :], after[1, :]) + qc3_3 = self._qc3(before[0, :] - after[0, :], before[1, :] - after[1, :]) QC['qc3.raw'] = [qc3_1] - QC['qc3.E4031'] = [qc3_2] + QC['qc3.E4031'] = [qc3_2] # Change to 'control' and 'blocker' ? QC['qc3.subtracted'] = [qc3_3] rseals = [qc_vals_before[0], qc_vals_after[0]] cms = [qc_vals_before[1], qc_vals_after[1]] rseriess = [qc_vals_before[2], qc_vals_after[2]] - qc4 = self.qc4(rseals, cms, rseriess) + qc4 = self._qc4(rseals, cms, rseriess) QC['qc4.rseal'] = [qc4[0]] QC['qc4.cm'] = [qc4[1]] QC['qc4.rseries'] = [qc4[2]] # indices where hERG peaks - qc5 = self.qc5(before[0, :], after[0, :], - self.qc5_win) + qc5 = self._qc5(before[0, :], after[0, :], + self._qc5_win) - qc5_1 = self.qc5_1(before[0, :], after[0, :], label='1') + qc5_1 = self._qc5_1(before[0, :], after[0, :], label='1') QC['qc5.staircase'] = [qc5] QC['qc5.1.staircase'] = [qc5_1] # Ensure thats the windows are correct by checking the voltage trace assert np.all( - np.abs(self.voltage[self.qc6_win[0]: self.qc6_win[1]] - 40.0))\ + np.abs(self.voltage[self._qc6_win[0]: self._qc6_win[1]] - 40.0))\ < 1e-8 assert np.all( - np.abs(self.voltage[self.qc6_1_win[0]: self.qc6_1_win[1]] - 40.0))\ + np.abs(self.voltage[self._qc6_1_win[0]: self._qc6_1_win[1]] - 40.0))\ < 1e-8 assert np.all( - np.abs(self.voltage[self.qc6_2_win[0]: self.qc6_2_win[1]] - 40.0))\ + np.abs(self.voltage[self._qc6_2_win[0]: self._qc6_2_win[1]] - 40.0))\ < 1e-8 QC['qc6.subtracted'] = [] QC['qc6.1.subtracted'] = [] QC['qc6.2.subtracted'] = [] for i in range(before.shape[0]): - qc6 = self.qc6((before[i, :] - after[i, :]), - self.qc6_win, label="0") - qc6_1 = self.qc6((before[i, :] - after[i, :]), - self.qc6_1_win, label="1") - qc6_2 = self.qc6((before[i, :] - after[i, :]), - self.qc6_2_win, label="2") + qc6 = self._qc6((before[i, :] - after[i, :]), + self._qc6_win, label='0') + qc6_1 = self._qc6((before[i, :] - after[i, :]), + self._qc6_1_win, label='1') + qc6_2 = self._qc6((before[i, :] - after[i, :]), + self._qc6_2_win, label='2') QC['qc6.subtracted'].append(qc6) QC['qc6.1.subtracted'].append(qc6_1) QC['qc6.2.subtracted'].append(qc6_2) - if self.plot_dir and self._debug: + if self._plot_dir and self._debug: fig = plt.figure(figsize=(8, 5)) ax = fig.subplots() ax.plot(times, (before - after).T, label='subtracted') ax.plot(times, (before).T, label='before') ax.plot(times, (after).T, label='after') - for l1, l2, l3, l4 in zip(self.qc5_win, self.qc6_win, - self.qc6_1_win, self.qc6_2_win): + for l1, l2, l3, l4 in zip(self._qc5_win, self._qc6_win, + self._qc6_1_win, self._qc6_2_win): plt.axvline(times[l1], c='#7f7f7f', label='qc5') plt.axvline(times[l2], c='#ff7f0e', ls='--', label='qc6') plt.axvline(times[l3], c='#2ca02c', ls='-.', label='qc6_1') @@ -259,12 +260,12 @@ def run_qc(self, voltage_steps, times, by_label = OrderedDict(zip(labels, handles)) fig.legend(by_label.values(), by_label.keys()) - fig.savefig(os.path.join(self.plot_dir, 'qc_debug.png')) + fig.savefig(os.path.join(self._plot_dir, 'qc_debug.png')) plt.close(fig) return QC - def qc1(self, rseal, cm, rseries): + def _qc1(self, rseal, cm, rseries): # Check R_seal, C_m, R_series within desired range if ( rseal is None @@ -272,7 +273,7 @@ def qc1(self, rseal, cm, rseries): or rseal > self.rsealc[1] or not np.isfinite(rseal) ): - self.logger.debug(f"rseal: {rseal}") + self.logger.debug(f'rseal: {rseal}') qc11 = False else: qc11 = True @@ -283,7 +284,7 @@ def qc1(self, rseal, cm, rseries): or cm > self.cmc[1] or not np.isfinite(cm) ): - self.logger.debug(f"cm: {cm}") + self.logger.debug(f'cm: {cm}') qc12 = False else: qc12 = True @@ -294,43 +295,43 @@ def qc1(self, rseal, cm, rseries): or rseries > self.rseriesc[1] or not np.isfinite(rseries) ): - self.logger.debug(f"rseries: {rseries}") + self.logger.debug(f'rseries: {rseries}') qc13 = False else: qc13 = True return [(qc11, rseal), (qc12, cm), (qc13, rseries)] - def qc2(self, recording): + def _qc2(self, recording): # Check SNR is good - noise = np.std(recording[:NOISE_LEN]) + noise = np.std(recording[:self.noise_len]) snr = (np.std(recording) / noise) ** 2 if snr < self.snrc or not np.isfinite(snr): - self.logger.debug(f"snr: {snr}") + self.logger.debug(f'snr: {snr}') result = False else: result = True return (result, snr) - def qc3(self, recording1, recording2): + def _qc3(self, recording1, recording2): # Check 2 sweeps similar - noise_1 = np.std(recording1[:NOISE_LEN]) - noise_2 = np.std(recording2[:NOISE_LEN]) + noise_1 = np.std(recording1[:self.noise_len]) + noise_2 = np.std(recording2[:self.noise_len]) rmsd0_1 = np.sqrt(np.mean((recording1) ** 2)) rmsd0_2 = np.sqrt(np.mean((recording2) ** 2)) rmsdc = max(np.mean([rmsd0_1, rmsd0_2]) * self.rmsd0c, np.mean([noise_1, noise_2]) * 6) rmsd = np.sqrt(np.mean((recording1 - recording2) ** 2)) if rmsd > rmsdc or not (np.isfinite(rmsd) and np.isfinite(rmsdc)): - self.logger.debug(f"rmsd: {rmsd}, rmsdc: {rmsdc}") + self.logger.debug(f'rmsd: {rmsd}, rmsdc: {rmsdc}') result = False else: result = True return (result, rmsdc - rmsd) - def qc4(self, rseals, cms, rseriess): + def _qc4(self, rseals, cms, rseriess): # Check R_seal, C_m, R_series stability # Require std/mean < x% qc41 = True @@ -344,7 +345,7 @@ def qc4(self, rseals, cms, rseriess): d_rseal = np.std(rseals) / np.mean(rseals) if d_rseal > self.rsealsc or not ( np.isfinite(np.mean(rseals)) and np.isfinite(np.std(rseals))): - self.logger.debug(f"d_rseal: {d_rseal}") + self.logger.debug(f'd_rseal: {d_rseal}') qc41 = False if None in list(cms): @@ -354,7 +355,7 @@ def qc4(self, rseals, cms, rseriess): d_cm = np.std(cms) / np.mean(cms) if d_cm > self.cmsc or not ( np.isfinite(np.mean(cms)) and np.isfinite(np.std(cms))): - self.logger.debug(f"d_cm: {d_cm}") + self.logger.debug(f'd_cm: {d_cm}') qc42 = False if None in list(rseriess): @@ -365,12 +366,12 @@ def qc4(self, rseals, cms, rseriess): if d_rseries > self.rseriessc or not ( np.isfinite(np.mean(rseriess)) and np.isfinite(np.std(rseriess))): - self.logger.debug(f"d_rseries: {d_rseries}") + self.logger.debug(f'd_rseries: {d_rseries}') qc43 = False return [(qc41, d_rseal), (qc42, d_cm), (qc43, d_rseries)] - def qc5(self, recording1, recording2, win=None, label=''): + def _qc5(self, recording1, recording2, win=None, label=''): # Check pharma peak value drops after E-4031 application # Require subtracted peak > 70% of the original peak if win is not None: @@ -378,12 +379,12 @@ def qc5(self, recording1, recording2, win=None, label=''): else: i, f = 0, None - if self.plot_dir and self._debug: + if self._plot_dir and self._debug: if win is not None: plt.axvspan(win[0], win[1], color='grey', alpha=.1) plt.plot(recording1, label='recording1') plt.plot(recording2, label='recording2') - plt.savefig(os.path.join(self.plot_dir, f"qc5_{label}")) + plt.savefig(os.path.join(self._plot_dir, f'qc5_{label}')) plt.clf() wherepeak = np.argmax(recording1[i:f]) @@ -392,14 +393,14 @@ def qc5(self, recording1, recording2, win=None, label=''): if (max_diff < max_diffc) or not (np.isfinite(max_diff) and np.isfinite(max_diffc)): - self.logger.debug(f"max_diff: {max_diff}, max_diffc: {max_diffc}") + self.logger.debug(f'max_diff: {max_diff}, max_diffc: {max_diffc}') result = False else: result = True return (result, max_diffc - max_diff) - def qc5_1(self, recording1, recording2, win=None, label=''): + def _qc5_1(self, recording1, recording2, win=None, label=''): # Check RMSD_0 drops after E-4031 application # Require RMSD_0 (after E-4031 / before) diff > 50% of RMSD_0 before if win is not None: @@ -407,14 +408,14 @@ def qc5_1(self, recording1, recording2, win=None, label=''): else: i, f = 0, -1 - if self.plot_dir and self._debug: + if self._plot_dir and self._debug: if win is not None: plt.axvspan(win[0], win[1], color='grey', alpha=.1) fig = plt.figure() ax = fig.subplots() ax.plot(recording1, label='recording1') ax.plot(recording2, label='recording2') - fig.savefig(os.path.join(self.plot_dir, f"qc5_{label}")) + fig.savefig(os.path.join(self._plot_dir, f'qc5_{label}')) plt.close(fig) rmsd0_diff = np.sqrt(np.mean(recording1[i:f] ** 2)) \ @@ -426,64 +427,58 @@ def qc5_1(self, recording1, recording2, win=None, label=''): if (rmsd0_diff < rmsd0_diffc) or not (np.isfinite(rmsd0_diff) and np.isfinite(rmsd0_diffc)): self.logger.debug( - f"rmsd0_diff: {rmsd0_diff}, rmsd0c: {rmsd0_diffc}") + f'rmsd0_diff: {rmsd0_diff}, rmsd0c: {rmsd0_diffc}') result = False else: result = True return (result, rmsd0_diffc - rmsd0_diff) - def qc6(self, recording1, win=None, label=''): + def _qc6(self, recording1, win=None, label=''): # Check subtracted staircase +40mV step up is non negative if win is not None: i, f = win else: i, f = 0, -1 - if self.plot_dir and self._debug: + if self._plot_dir and self._debug: if win is not None: plt.axvspan(win[0], win[1], color='grey', alpha=.1) plt.plot(recording1, label='recording1') - plt.savefig(os.path.join(self.plot_dir, f"qc6_{label}")) + plt.savefig(os.path.join(self._plot_dir, f'qc6_{label}')) plt.clf() val = np.mean(recording1[i:f]) - - # valc = -0.005 * np.abs(np.sqrt(np.mean((recording1) ** 2))) # or just 0 - valc = self.negative_tolc * np.std(recording1[:NOISE_LEN]) + valc = self.negative_tolc * np.std(recording1[:self.noise_len]) if (val < valc) or not (np.isfinite(val) and np.isfinite(valc)): - self.logger.debug(f"qc6_{label} val: {val}, valc: {valc}") + self.logger.debug(f'qc6_{label} val: {val}, valc: {valc}') result = False else: result = True return (result, valc - val) - def filter_capacitive_spikes(self, current, times, voltage_step_times): + def _filter_capacitive_spikes(self, current, times, voltage_step_times): """ - Set values to 0 where they lie less that self.removal time after a change in voltage + Set values to 0 where they lie less than ``removal_time`` after a change in voltage. @param current: The observed current @param times: the times at which the current was observed @param voltage_step_times: the times at which there are discontinuities in Vcmd """ + if len(current.shape) != 2: + raise ValueError('Current must have 2 dimensions (sweep, values)') - for tstart, tend in zip(voltage_step_times, - np.append(voltage_step_times[1:], np.inf)): + voltage_step_ends = np.append(voltage_step_times[1:], np.inf) + for tstart, tend in zip(voltage_step_times, voltage_step_ends): win_end = tstart + self.removal_time win_end = min(tend, win_end) i_start = np.argmax(times >= tstart) i_end = np.argmax(times > win_end) - if i_end == 0: break - if len(current.shape) == 2: - # When would this happen? Is one dimension always 1 ? - current[:, i_start:i_end] = 0 - elif len(current.shape) == 1: - current[i_start:i_end] = 0 - else: - raise ValueError("Current must have 1 or 2 dimensions") + + current[:, i_start:i_end] = 0 return current diff --git a/tests/test_herg_qc.py b/tests/test_herg_qc.py index ebf64ef0..10460f5e 100755 --- a/tests/test_herg_qc.py +++ b/tests/test_herg_qc.py @@ -1,5 +1,4 @@ #!/usr/bin/env python3 -import copy import os import string import unittest @@ -7,7 +6,7 @@ import numpy as np from syncropatch_export.trace import Trace -from pcpostprocess.hergQC import NOISE_LEN, hERGQC +from pcpostprocess.hergQC import hERGQC def all_passed(result): @@ -45,7 +44,8 @@ def setUp(self): for i in range(1, 25) ] - sampling_rate = int(1.0 / (self.times[1] - self.times[0])) # in kHz + # in kHz + self.sampling_rate = int(1 / (self.times[1] - self.times[0])) #  Assume that there are no discontinuities at the start or end of ramps voltage_protocol = trace_before.get_voltage_protocol() @@ -55,29 +55,37 @@ def setUp(self): if vend == vstart ] - plot_dir = "test_output" - os.makedirs(plot_dir, exist_ok=True) + def create_hergqc(self, plot_dir=None): + """ + Creates and returns a hERGQC object. - self.hergqc = hERGQC( - sampling_rate=sampling_rate, + If a ``plot_dir`` is set this will be used for debug output + """ + # Set this to True to generate plot directories + debug = False + + if debug and plot_dir is not None: + plot_dir = os.path.join('test_output', plot_dir) + os.makedirs(plot_dir, exist_ok=True) + else: + plot_dir = None + + return hERGQC( + sampling_rate=self.sampling_rate, plot_dir=plot_dir, voltage=self.voltage, ) - self.hergqc._debug = False # Turn this on to output plots - - def clone_hergqc(self, plot_dir): - hergqc = copy.deepcopy(self.hergqc) - plot_dir = os.path.join(hergqc.plot_dir, plot_dir) - os.makedirs(plot_dir, exist_ok=True) - hergqc.plot_dir = plot_dir - return hergqc def test_qc_inputs(self): + hergqc = self.create_hergqc() times = self.times - voltage = self.hergqc.voltage - qc6_win = self.hergqc.qc6_win - qc6_1_win = self.hergqc.qc6_1_win - qc6_2_win = self.hergqc.qc6_2_win + + # TODO: This should work some nicer way, without accessing private properties + # But first we probably need to stop hardcoding these windows + voltage = hergqc._voltage + qc6_win = hergqc._qc6_win + qc6_1_win = hergqc._qc6_1_win + qc6_2_win = hergqc._qc6_2_win self.assertTrue(np.all(np.isfinite(times))) self.assertTrue(np.all(np.isfinite(voltage))) @@ -89,7 +97,7 @@ def test_qc_inputs(self): def test_qc1(self): # qc1 checks that rseal, cm, rseries are within range - hergqc = self.clone_hergqc("test_qc1") + hergqc = self.create_hergqc('qc1') rseal_lo, rseal_hi = 1e8, 1e12 rseal_mid = (rseal_lo + rseal_hi) / 2 @@ -297,7 +305,7 @@ def test_qc1(self): def test_qc2(self): # qc2 checks that raw and subtracted SNR are above a minimum threshold - hergqc = self.clone_hergqc("test_qc2") + hergqc = self.create_hergqc('qc2') test_matrix = [ (10, True, 8082.1), @@ -310,7 +318,7 @@ def test_qc2(self): ] for i, ex_pass, ex_snr in test_matrix: - recording = np.asarray([0, 0.1] * (NOISE_LEN // 2) + [i] * 500) + recording = np.asarray([0, 0.1] * 100 + [i] * 500) pass_, snr = hergqc.qc2(recording) self.assertAlmostEqual( snr, ex_snr, 1, f"QC2: ({i}) {snr} != {ex_snr}" @@ -349,7 +357,7 @@ def test_qc2(self): def test_qc3(self): # qc3 checks that rmsd of two sweeps are similar - hergqc = self.clone_hergqc("test_qc3") + hergqc = self.create_hergqc('qc3') # Test with same noise, different signal test_matrix = [ @@ -363,10 +371,10 @@ def test_qc3(self): (10, False, -0.8), ] - recording1 = np.asarray([0, 0.1] * (NOISE_LEN // 2) + [40] * 500) + recording1 = np.asarray([0, 0.1] * 100 + [40] * 500) for i, ex_pass, ex_d_rmsd in test_matrix: recording2 = np.asarray( - [0, 0.1] * (NOISE_LEN // 2) + [40 + i] * 500 + [0, 0.1] * 100 + [40 + i] * 500 ) pass_, d_rmsd = hergqc.qc3(recording1, recording2) self.assertAlmostEqual( @@ -388,10 +396,10 @@ def test_qc3(self): (100, True, 11.4), ] - recording1 = np.asarray([0, 0.1] * (NOISE_LEN // 2) + [40] * 500) + recording1 = np.asarray([0, 0.1] * 100 + [40] * 500) for i, ex_pass, ex_d_rmsd in test_matrix: recording2 = np.asarray( - [0, 0.1 * i] * (NOISE_LEN // 2) + [40] * 500 + [0, 0.1 * i] * 100 + [40] * 500 ) pass_, d_rmsd = hergqc.qc3(recording1, recording2) self.assertAlmostEqual( @@ -456,7 +464,7 @@ def test_qc3(self): def test_qc4(self): # qc4 checks that rseal, cm, rseries are similar before/after E-4031 change - hergqc = self.clone_hergqc("test_qc4") + hergqc = self.create_hergqc('qc4') r_lo, r_hi = 1e6, 3e7 c_lo, c_hi = 1e-12, 1e-10 @@ -617,7 +625,7 @@ def test_qc4(self): def test_qc5(self): # qc5 checks that the maximum current during the second half of the # staircase changes by at least 75% of the raw trace after E-4031 addition - hergqc = self.clone_hergqc("test_qc5") + hergqc = self.create_hergqc('qc5') test_matrix = [ (-1.0, True, -12.5), @@ -631,10 +639,10 @@ def test_qc5(self): (1.0, False, 7.5), ] - recording1 = np.asarray([0, 0.1] * (NOISE_LEN // 2) + [10] * 500) + recording1 = np.asarray([0, 0.1] * 100 + [10] * 500) for i, ex_pass, ex_d_max_diff in test_matrix: recording2 = np.asarray( - [0, 0.1] * (NOISE_LEN // 2) + [10 * i] * 500 + [0, 0.1] * 100 + [10 * i] * 500 ) pass_, d_max_diff = hergqc.qc5(recording1, recording2) self.assertAlmostEqual( @@ -689,7 +697,7 @@ def test_qc5(self): def test_qc5_1(self): # qc5_1 checks that the RMSD to zero of staircase protocol changes # by at least 50% of the raw trace after E-4031 addition. - hergqc = self.clone_hergqc("test_qc5_1") + hergqc = self.create_hergqc('qc5_1') test_matrix = [ (-1.0, False, 4.22), @@ -704,10 +712,10 @@ def test_qc5_1(self): (1.0, False, 4.22), ] - recording1 = np.asarray([0, 0.1] * (NOISE_LEN // 2) + [10] * 500) + recording1 = np.asarray([0, 0.1] * 100 + [10] * 500) for i, ex_pass, ex_d_max_diff in test_matrix: recording2 = np.asarray( - [0, 0.1] * (NOISE_LEN // 2) + [10 * i] * 500 + [0, 0.1] * 100 + [10 * i] * 500 ) pass_, d_max_diff = hergqc.qc5_1(recording1, recording2) self.assertAlmostEqual( @@ -765,7 +773,7 @@ def test_qc5_1(self): def test_qc6(self): # qc6 checks that the first step up to +40 mV, before the staircase, in # the subtracted trace is bigger than -2 x estimated noise level. - hergqc = self.clone_hergqc("test_qc6") + hergqc = self.create_hergqc('qc6') test_matrix = [ (-100, False, 9.9), @@ -779,9 +787,9 @@ def test_qc6(self): for i, ex_pass, ex_d_val in test_matrix: recording = np.asarray( - [0, 0.1] * (NOISE_LEN // 2) + [0.1 * i] * 500 + [0, 0.1] * 100 + [0.1 * i] * 500 ) - pass_, d_val = hergqc.qc6(recording, win=[NOISE_LEN, -1]) + pass_, d_val = hergqc.qc6(recording, win=[200, -1]) self.assertAlmostEqual( d_val, ex_d_val, @@ -888,7 +896,7 @@ def test_qc6(self): def test_run_qc(self): # Test all wells - hergqc = self.clone_hergqc("test_run_qc") + hergqc = self.create_hergqc('run_qc') failed_wells = [ 'A04', 'A05', 'A06', 'A07', 'A08', 'A10', 'A11', 'A12', 'A13', 'A15', From f1b1ba24dd1a37faed64c7b2cccdadf6bfb72c64 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Tue, 4 Nov 2025 17:01:16 +0000 Subject: [PATCH 04/13] Added docstrings to hergQC --- pcpostprocess/hergQC.py | 342 +++++++++++++++++++++++----------------- tests/test_herg_qc.py | 17 +- 2 files changed, 208 insertions(+), 151 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index 4e213e47..bdf9a2d1 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -66,7 +66,6 @@ class hERGQC: @param removal_time: Number of time units to remove after each step in the protocol @param noise_len: Number of initial samples during which the signal is flat, to use for noise estimates. @param plot_dir: An optional directory to store plots in - """ def __init__(self, voltage, sampling_rate=5, removal_time=5, noise_len=200, plot_dir=None): @@ -75,20 +74,19 @@ def __init__(self, voltage, sampling_rate=5, removal_time=5, noise_len=200, self.sampling_rate = sampling_rate self.removal_time = removal_time self.noise_len = int(noise_len) - self._plot_dir = plot_dir - - self._n_qc = len(QCDict.labels) + # Passing in a plot dir enables debug mode + self._plot_dir = plot_dir self.logger = logging.getLogger(__name__) - self._debug = True - if self._debug: + if self._plot_dir is not None: self.logger.setLevel(logging.DEBUG) - # https://github.com/CardiacModelling/pcpostprocess/issues/42 + # https://github.com/CardiacModelling/pcpostprocess/issues/42 + self._plot_dir = plot_dir # Define all thresholds # TODO: These should be args? Or perhaps this is good so that this QC - # class can be extended + # class can be extended? # qc1 self.rsealc = [1e8, 1e12] # in Ohm, converted from [0.1, 1000] GOhm self.cmc = [1e-12, 1e-10] # in F, converted from [1, 100] pF @@ -103,43 +101,32 @@ def __init__(self, voltage, sampling_rate=5, removal_time=5, noise_len=200, self.rseriessc = 0.5 # qc5 self.max_diffc = 0.75 - # self._qc5_win = [3275, 5750] # indices with hERG screening peak! + # self.qc5_win = [3275, 5750] # indices with hERG screening peak! # indices where hERG could peak (for different temperatures) - self._qc5_win = np.array([8550 + 400, 10950 + 400]) * sampling_rate + self.qc5_win = np.array([8550 + 400, 10950 + 400]) * sampling_rate # qc5_1 self.rmsd0_diffc = 0.5 # qc6 self.negative_tolc = -2 ''' # These are for `whole` (just staircase) protocol at 5kHz - self._qc6_win = [3000, 7000] # indices for first +40 mV - self._qc6_1_win = [35250, 37250] # indices for second +40 mV - self._qc6_2_win = [40250, 42250] # indices for third +40 mV + self.qc6_win = [3000, 7000] # indices for first +40 mV + self.qc6_1_win = [35250, 37250] # indices for second +40 mV + self.qc6_2_win = [40250, 42250] # indices for third +40 mV ''' # These are for `staircaseramp` protocol # Firstly, indices for 1st +40 mV - self._qc6_win = np.array([1000, 1800]) * sampling_rate + self.qc6_win = np.array([1000, 1800]) * sampling_rate # indices for 2nd +40 mV - self._qc6_1_win = np.array([7450, 7850]) * sampling_rate + self.qc6_1_win = np.array([7450, 7850]) * sampling_rate # indices for 3rd +40 mV - self._qc6_2_win = np.array([8450, 8850]) * sampling_rate + self.qc6_2_win = np.array([8450, 8850]) * sampling_rate # Ensure these indices are integers - self._qc5_win = self._qc5_win.astype(int) - self._qc6_win = self._qc6_win.astype(int) - self._qc6_1_win = self._qc6_1_win.astype(int) - self._qc6_2_win = self._qc6_2_win.astype(int) - - @property - def plot_dir(self): - return self._plot_dir - - @plot_dir.setter - def plot_dir(self, path): - self._plot_dir = path - - def set_debug(self, debug): - self._debug = debug + self.qc5_win = self.qc5_win.astype(int) + self.qc6_win = self.qc6_win.astype(int) + self.qc6_1_win = self.qc6_1_win.astype(int) + self.qc6_2_win = self.qc6_2_win.astype(int) def run_qc(self, voltage_steps, times, before, after, qc_vals_before, qc_vals_after, n_sweeps=None): @@ -147,14 +134,16 @@ def run_qc(self, voltage_steps, times, before, after, qc_vals_before, @param voltage_steps is a list of times at which there are discontinuities in Vcmd @param times is the array of observation times - @param before are the before-drug current traces, in an array with sweeps on the first index and values on the second. - @param after is the post-drug current traces, in an array with sweeps on the first index and values on the second. - @param qc_vals_before is an array of values for each before-drug sweep where each row is (Rseal, Cm, Rseries) - @param qc_vals_after is an array of values for each after-drug sweep where each row is (Rseal, Cm, Rseries) + @param before are the before-drug current traces, in an array with + sweeps on the first index and values on the second. + @param after is the post-drug current traces, in an array with sweeps + on the first index and values on the second. + @param qc_vals_before is a sequence (Rseal, Cm, Rseries) + @param qc_vals_after is a sequence (Rseal, Cm, Rseries) @n_sweeps is the number of sweeps to be run QC on. - @returns A :class:`QCDict` with the test results. """ + # TODO: Why doesn't each sweep have its own "qc_vals" ? if not n_sweeps: n_sweeps = len(before) @@ -170,8 +159,8 @@ def run_qc(self, voltage_steps, times, before, after, qc_vals_before, if (None in qc_vals_before) or (None in qc_vals_after): return QC - qc1_1 = self._qc1(*qc_vals_before) - qc1_2 = self._qc1(*qc_vals_after) + qc1_1 = self.qc1(*qc_vals_before) + qc1_2 = self.qc1(*qc_vals_after) QC['qc1.rseal'] = [qc1_1[0], qc1_2[0]] QC['qc1.cm'] = [qc1_1[1], qc1_2[1]] @@ -180,16 +169,15 @@ def run_qc(self, voltage_steps, times, before, after, qc_vals_before, QC['qc2.raw'] = [] QC['qc2.subtracted'] = [] for i in range(n_sweeps): - qc2_1 = self._qc2(before[i]) - qc2_2 = self._qc2(before[i] - after[i]) + qc2_1 = self.qc2(before[i]) + qc2_2 = self.qc2(before[i] - after[i]) QC['qc2.raw'].append(qc2_1) QC['qc2.subtracted'].append(qc2_2) - qc3_1 = self._qc3(before[0, :], before[1, :]) - qc3_2 = self._qc3(after[0, :], after[1, :]) - qc3_3 = self._qc3(before[0, :] - after[0, :], - before[1, :] - after[1, :]) + qc3_1 = self.qc3(before[0, :], before[1, :]) + qc3_2 = self.qc3(after[0, :], after[1, :]) + qc3_3 = self.qc3(before[0, :] - after[0, :], before[1, :] - after[1, :]) QC['qc3.raw'] = [qc3_1] QC['qc3.E4031'] = [qc3_2] # Change to 'control' and 'blocker' ? @@ -198,55 +186,53 @@ def run_qc(self, voltage_steps, times, before, after, qc_vals_before, rseals = [qc_vals_before[0], qc_vals_after[0]] cms = [qc_vals_before[1], qc_vals_after[1]] rseriess = [qc_vals_before[2], qc_vals_after[2]] - qc4 = self._qc4(rseals, cms, rseriess) + qc4 = self.qc4(rseals, cms, rseriess) QC['qc4.rseal'] = [qc4[0]] QC['qc4.cm'] = [qc4[1]] QC['qc4.rseries'] = [qc4[2]] # indices where hERG peaks - qc5 = self._qc5(before[0, :], after[0, :], - self._qc5_win) - - qc5_1 = self._qc5_1(before[0, :], after[0, :], label='1') + qc5 = self.qc5(before[0, :], after[0, :], self.qc5_win) + qc5_1 = self.qc5_1(before[0, :], after[0, :], label='1') QC['qc5.staircase'] = [qc5] QC['qc5.1.staircase'] = [qc5_1] # Ensure thats the windows are correct by checking the voltage trace assert np.all( - np.abs(self.voltage[self._qc6_win[0]: self._qc6_win[1]] - 40.0))\ + np.abs(self.voltage[self.qc6_win[0]: self.qc6_win[1]] - 40.0))\ < 1e-8 assert np.all( - np.abs(self.voltage[self._qc6_1_win[0]: self._qc6_1_win[1]] - 40.0))\ + np.abs(self.voltage[self.qc6_1_win[0]: self.qc6_1_win[1]] - 40.0))\ < 1e-8 assert np.all( - np.abs(self.voltage[self._qc6_2_win[0]: self._qc6_2_win[1]] - 40.0))\ + np.abs(self.voltage[self.qc6_2_win[0]: self.qc6_2_win[1]] - 40.0))\ < 1e-8 QC['qc6.subtracted'] = [] QC['qc6.1.subtracted'] = [] QC['qc6.2.subtracted'] = [] for i in range(before.shape[0]): - qc6 = self._qc6((before[i, :] - after[i, :]), - self._qc6_win, label='0') - qc6_1 = self._qc6((before[i, :] - after[i, :]), - self._qc6_1_win, label='1') - qc6_2 = self._qc6((before[i, :] - after[i, :]), - self._qc6_2_win, label='2') + qc6 = self.qc6( + before[i, :] - after[i, :], self.qc6_win, label='0') + qc6_1 = self.qc6( + before[i, :] - after[i, :], self.qc6_1_win, label='1') + qc6_2 = self.qc6( + before[i, :] - after[i, :], self.qc6_2_win, label='2') QC['qc6.subtracted'].append(qc6) QC['qc6.1.subtracted'].append(qc6_1) QC['qc6.2.subtracted'].append(qc6_2) - if self._plot_dir and self._debug: + if self._plot_dir is not None: fig = plt.figure(figsize=(8, 5)) ax = fig.subplots() ax.plot(times, (before - after).T, label='subtracted') ax.plot(times, (before).T, label='before') ax.plot(times, (after).T, label='after') - for l1, l2, l3, l4 in zip(self._qc5_win, self._qc6_win, - self._qc6_1_win, self._qc6_2_win): + for l1, l2, l3, l4 in zip(self.qc5_win, self.qc6_win, + self.qc6_1_win, self.qc6_2_win): plt.axvline(times[l1], c='#7f7f7f', label='qc5') plt.axvline(times[l2], c='#ff7f0e', ls='--', label='qc6') plt.axvline(times[l3], c='#2ca02c', ls='-.', label='qc6_1') @@ -265,75 +251,113 @@ def run_qc(self, voltage_steps, times, before, after, qc_vals_before, return QC - def _qc1(self, rseal, cm, rseries): + def qc1(self, rseal, cm, rseries): + """ + Checks that the given ``rseal``, ``cm`` and ``rseries`` are within the + desired range. + + @param rseal A scalar indicating a seal resistance in Ohm + @param cm A scalar cell capacitance in Farad + @param rseries A scalar series resistance in Ohm + @returns A tuple ``((qc11_passed, rseal), (qc12_passed, cm), (qc13_passed, rseries))``. + """ + # TODO Is there any reason these are public, other than that it was + # convenient for testing? + # TODO Could just be 1 method called 3 times + # TODO Should boundaries be arguments? + # Check R_seal, C_m, R_series within desired range - if ( - rseal is None - or rseal < self.rsealc[0] - or rseal > self.rsealc[1] - or not np.isfinite(rseal) - ): + qc11 = not (rseal is None + or rseal < self.rsealc[0] + or rseal > self.rsealc[1] + or not np.isfinite(rseal)) + qc12 = not (cm is None + or cm < self.cmc[0] + or cm > self.cmc[1] + or not np.isfinite(cm)) + qc13 = not (rseries is None + or rseries < self.rseriesc[0] + or rseries > self.rseriesc[1] + or not np.isfinite(rseries)) + + if not qc11: self.logger.debug(f'rseal: {rseal}') - qc11 = False - else: - qc11 = True - - if ( - cm is None - or cm < self.cmc[0] - or cm > self.cmc[1] - or not np.isfinite(cm) - ): + if not qc12: self.logger.debug(f'cm: {cm}') - qc12 = False - else: - qc12 = True - - if ( - rseries is None - or rseries < self.rseriesc[0] - or rseries > self.rseriesc[1] - or not np.isfinite(rseries) - ): + if not qc13: self.logger.debug(f'rseries: {rseries}') - qc13 = False - else: - qc13 = True return [(qc11, rseal), (qc12, cm), (qc13, rseries)] - def _qc2(self, recording): - # Check SNR is good - noise = np.std(recording[:self.noise_len]) - snr = (np.std(recording) / noise) ** 2 - if snr < self.snrc or not np.isfinite(snr): + def qc2(self, recording): + """ + Checks that the signal-to-noise ratio is below a certain threshold. + + Here the signal-to-noise ratio is defined as + ``(std(recording) / std(noise))``, where ``noise`` is the initial part + of ``recording``. + + @param recording A 1-d numpy array containing recorded currents. + @returns a tuple ``(passed, signal-to-noise-ratio)``. + """ + snr = (np.std(recording) / np.std(recording[:self.noise_len])) ** 2 + passed = snr >= self.snrc and np.isfinite(snr) + if not passed: self.logger.debug(f'snr: {snr}') - result = False - else: - result = True + return (passed, snr) - return (result, snr) - - def _qc3(self, recording1, recording2): - # Check 2 sweeps similar - noise_1 = np.std(recording1[:self.noise_len]) - noise_2 = np.std(recording2[:self.noise_len]) - rmsd0_1 = np.sqrt(np.mean((recording1) ** 2)) - rmsd0_2 = np.sqrt(np.mean((recording2) ** 2)) - rmsdc = max(np.mean([rmsd0_1, rmsd0_2]) * self.rmsd0c, - np.mean([noise_1, noise_2]) * 6) - rmsd = np.sqrt(np.mean((recording1 - recording2) ** 2)) - if rmsd > rmsdc or not (np.isfinite(rmsd) and np.isfinite(rmsdc)): - self.logger.debug(f'rmsd: {rmsd}, rmsdc: {rmsdc}') - result = False - else: - result = True + def qc3(self, recording1, recording2): + """ + Checks that ``recording1`` and ``recording2`` are similar, using a + measure based on root-mean-square-distance. + + Particularly, it checks whether + ``RMSD(recording1, recording2) < threshold``, + where ``threshold`` is determined as the maximum of + ``w * (RMSD(recording1, 0) + RMSD(recording2, 0))`` and + ``6 * (std(noise1) + std(noise2))``. + + @param recording1 A 1-d numpy array containing recorded currents + @param recording2 The same, but for a repeat measurement + @returns a tuple ``(passed, RMSD - threshold)`` + """ + def rmsd(x, y=0): + return np.sqrt(np.mean((x - y)**2)) + + t1 = self.rmsd0c * np.mean((rmsd(recording1), rmsd(recording2))) + t2 = 6 * np.mean(( + np.std(recording1[:self.noise_len]), + np.std(recording2[:self.noise_len]))) + # TODO: This second cut-off is not documented in the paper + t = max(t1, t2) + r = rmsd(recording1, recording2) + passed = r <= t and np.isfinite(t) and np.isfinite(r) + if not passed: + self.logger.debug(f'rmsd: {r}, rmsdc: {t}') + + # TODO: Is this the best thing to return? Why not t _and_ r? + return (passed, t - r) + + def qc4(self, rseals, cms, rseriess): + """ + Checks that Rseal, Cm, and Rseries values remain steady during an + experiment. + + For each quantity, it checks if ``std(x) / mean(x) <= cutoff``. + + @param rseals A list of Rseal values, for repeated experiments, in Ohm. + @param cms A list of Cm values, in Farad. + @param rseriess A list of Rseries values, in Ohm + @returns a tuple ``((qc41_passed, ex41), (qc42_passed, ex42), (qc43_passed, ex43))`` + where ``exY`` is the amount by which the criterion exceeded the threshold. + """ + # TODO: Run this on values for all sweeps, not 1 before and after ? + # Or, if using only 2, write as just |a-b|/(a+b) - return (result, rmsdc - rmsd) + # TODO: This uses numpy's standard `std` call without the bias + # correction. Do we want to use the bias corrected one instead? + # matters for above, as becomes 2|a-b|/(a+b) ! - def _qc4(self, rseals, cms, rseriess): - # Check R_seal, C_m, R_series stability - # Require std/mean < x% qc41 = True qc42 = True qc43 = True @@ -371,15 +395,29 @@ def _qc4(self, rseals, cms, rseriess): return [(qc41, d_rseal), (qc42, d_cm), (qc43, d_rseries)] - def _qc5(self, recording1, recording2, win=None, label=''): - # Check pharma peak value drops after E-4031 application - # Require subtracted peak > 70% of the original peak + def qc5(self, recording1, recording2, win=None, label=None): + """ + Checks whether the peak current in ``recording1`` has been blocked in + ``recording2``. + + First, find an ``i`` such that ``recording1[i]`` is the largest + (positive) current (if a window is given, only this window is + searched). Then, checks whether + ``(recording1[i] - recording2[i]) / recording1[i] >= 0.75`` (in other + words whether at least 75% of the peak current has been blocked). + + @param recording1 A staircase before blocker application + @param recording2 A staircase with a strong blocker, e.g. E-4031 + @param win An optional window (lower and upper indices) within which to search for peak current + @param label An optional label for a plot + @returns a tuple ``(passed, 0.75 * recording1[i] - (recording1[i] - recording2[i]))``. + """ if win is not None: i, f = win else: i, f = 0, None - if self._plot_dir and self._debug: + if self._plot_dir is not None and label is not None: if win is not None: plt.axvspan(win[0], win[1], color='grey', alpha=.1) plt.plot(recording1, label='recording1') @@ -398,19 +436,25 @@ def _qc5(self, recording1, recording2, win=None, label=''): else: result = True + # TODO: More sensible to return max_diff / recording1[i:f][wherepeak] here? return (result, max_diffc - max_diff) - def _qc5_1(self, recording1, recording2, win=None, label=''): - # Check RMSD_0 drops after E-4031 application - # Require RMSD_0 (after E-4031 / before) diff > 50% of RMSD_0 before - if win is not None: - i, f = win - else: - i, f = 0, -1 + def qc5_1(self, recording1, recording2, label=None): + """ + Checks whether root-mean-squared current in ``recording1`` has been + blocked in ``recording2``, - if self._plot_dir and self._debug: - if win is not None: - plt.axvspan(win[0], win[1], color='grey', alpha=.1) + The test passes if + ``(RMSD(recording2, 0) - RMSD(recording1, 0)) / RMSD(recording1, 0) <= 0.5`` + + @param recording1 A staircase before blocker application + @param recording2 A staircase with a strong blocker, e.g. E-4031 + @param label An optional label for a plot + @returns a tuple ``(passed, 0.5 * recording1[i] - (recording1[i] - recording2[i]))``. + """ + # TODO: Just rephrase this as RMSD(I2) / RMSD(I1) <= 0.5 + + if self._plot_dir is not None and label is not None: fig = plt.figure() ax = fig.subplots() ax.plot(recording1, label='recording1') @@ -418,11 +462,11 @@ def _qc5_1(self, recording1, recording2, win=None, label=''): fig.savefig(os.path.join(self._plot_dir, f'qc5_{label}')) plt.close(fig) - rmsd0_diff = np.sqrt(np.mean(recording1[i:f] ** 2)) \ - - np.sqrt(np.mean(recording2[i:f] ** 2)) + def rmsd(x, y=0): + return np.sqrt(np.mean((x - y)**2)) - rmsd0_diffc = self.rmsd0_diffc *\ - np.sqrt(np.mean(recording1[i:f] ** 2)) + rmsd0_diff = rmsd(recording1) - rmsd(recording2) + rmsd0_diffc = self.rmsd0_diffc * rmsd(recording1) if (rmsd0_diff < rmsd0_diffc) or not (np.isfinite(rmsd0_diff) and np.isfinite(rmsd0_diffc)): @@ -432,22 +476,32 @@ def _qc5_1(self, recording1, recording2, win=None, label=''): else: result = True + # TODO Just return rmsd(I2) / RMSD(I1) return (result, rmsd0_diffc - rmsd0_diff) - def _qc6(self, recording1, win=None, label=''): + def qc6(self, recording1, win, label=None): + """ + Checks that the current in a particular window is non-negative. + + Instead of comparing the whole current to 0, we compare the mean + current to -2 times the noise level (estimed from the start of the + trace). + + @param recording1 The full current trace + @param win A tuple ``(i, j)`` such that ``recording[i:j]`` is the window to check + @param label A label for an optional plot + @returns a tuple ``(passed, -2 * noise - + """ # Check subtracted staircase +40mV step up is non negative - if win is not None: - i, f = win - else: - i, f = 0, -1 - if self._plot_dir and self._debug: + if self._plot_dir is not None and label is not None: if win is not None: plt.axvspan(win[0], win[1], color='grey', alpha=.1) plt.plot(recording1, label='recording1') plt.savefig(os.path.join(self._plot_dir, f'qc6_{label}')) plt.clf() + i, f = win val = np.mean(recording1[i:f]) valc = self.negative_tolc * np.std(recording1[:self.noise_len]) if (val < valc) or not (np.isfinite(val) @@ -457,6 +511,7 @@ def _qc6(self, recording1, win=None, label=''): else: result = True + # TODO Return val / valc here? return (result, valc - val) def _filter_capacitive_spikes(self, current, times, voltage_step_times): @@ -466,6 +521,7 @@ def _filter_capacitive_spikes(self, current, times, voltage_step_times): @param current: The observed current @param times: the times at which the current was observed @param voltage_step_times: the times at which there are discontinuities in Vcmd + @returns the ``current`` with some samples set to 0 """ if len(current.shape) != 2: raise ValueError('Current must have 2 dimensions (sweep, values)') diff --git a/tests/test_herg_qc.py b/tests/test_herg_qc.py index 10460f5e..f6217472 100755 --- a/tests/test_herg_qc.py +++ b/tests/test_herg_qc.py @@ -80,12 +80,13 @@ def test_qc_inputs(self): hergqc = self.create_hergqc() times = self.times - # TODO: This should work some nicer way, without accessing private properties - # But first we probably need to stop hardcoding these windows - voltage = hergqc._voltage - qc6_win = hergqc._qc6_win - qc6_1_win = hergqc._qc6_1_win - qc6_2_win = hergqc._qc6_2_win + # TODO: This should work some nicer way, without accessing what should + # really be private properties. But first we probably need to + # stop hardcoding these windows + voltage = hergqc.voltage + qc6_win = hergqc.qc6_win + qc6_1_win = hergqc.qc6_1_win + qc6_2_win = hergqc.qc6_2_win self.assertTrue(np.all(np.isfinite(times))) self.assertTrue(np.all(np.isfinite(voltage))) @@ -700,7 +701,7 @@ def test_qc5_1(self): hergqc = self.create_hergqc('qc5_1') test_matrix = [ - (-1.0, False, 4.22), + (-1.0, False, 4.23), (-0.5, False, 0), (-0.412, True, -0.74), (-0.4, True, -0.84), @@ -709,7 +710,7 @@ def test_qc5_1(self): (0.4, True, -0.84), (0.412, True, -0.74), (0.5, False, 0), - (1.0, False, 4.22), + (1.0, False, 4.23), ] recording1 = np.asarray([0, 0.1] * 100 + [10] * 500) From c403a3a24212a259ab911cfa51811e0aa24c705f Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Tue, 4 Nov 2025 17:11:18 +0000 Subject: [PATCH 05/13] Made filter_capacitive_spikes public again --- pcpostprocess/hergQC.py | 23 +++++++++-------------- tests/test_herg_qc.py | 2 +- 2 files changed, 10 insertions(+), 15 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index bdf9a2d1..3fbd3f48 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -148,8 +148,8 @@ def run_qc(self, voltage_steps, times, before, after, qc_vals_before, if not n_sweeps: n_sweeps = len(before) - before = self._filter_capacitive_spikes(np.array(before), times, voltage_steps) - after = self._filter_capacitive_spikes(np.array(after), times, voltage_steps) + before = self.filter_capacitive_spikes(np.array(before), times, voltage_steps) + after = self.filter_capacitive_spikes(np.array(after), times, voltage_steps) QC = QCDict() @@ -395,10 +395,10 @@ def qc4(self, rseals, cms, rseriess): return [(qc41, d_rseal), (qc42, d_cm), (qc43, d_rseries)] - def qc5(self, recording1, recording2, win=None, label=None): + def qc5(self, recording1, recording2, win, label=None): """ - Checks whether the peak current in ``recording1`` has been blocked in - ``recording2``. + Checks whether the peak current in ``recording1`` during a given window + has been blocked in ``recording2``. First, find an ``i`` such that ``recording1[i]`` is the largest (positive) current (if a window is given, only this window is @@ -408,23 +408,18 @@ def qc5(self, recording1, recording2, win=None, label=None): @param recording1 A staircase before blocker application @param recording2 A staircase with a strong blocker, e.g. E-4031 - @param win An optional window (lower and upper indices) within which to search for peak current + @param win A tuple ``(i, j)`` such that ``recording[i:j]`` is the window to search in @param label An optional label for a plot @returns a tuple ``(passed, 0.75 * recording1[i] - (recording1[i] - recording2[i]))``. """ - if win is not None: - i, f = win - else: - i, f = 0, None - if self._plot_dir is not None and label is not None: - if win is not None: - plt.axvspan(win[0], win[1], color='grey', alpha=.1) + plt.axvspan(win[0], win[1], color='grey', alpha=.1) plt.plot(recording1, label='recording1') plt.plot(recording2, label='recording2') plt.savefig(os.path.join(self._plot_dir, f'qc5_{label}')) plt.clf() + i, f = win wherepeak = np.argmax(recording1[i:f]) max_diff = recording1[i:f][wherepeak] - recording2[i:f][wherepeak] max_diffc = self.max_diffc * recording1[i:f][wherepeak] @@ -514,7 +509,7 @@ def qc6(self, recording1, win, label=None): # TODO Return val / valc here? return (result, valc - val) - def _filter_capacitive_spikes(self, current, times, voltage_step_times): + def filter_capacitive_spikes(self, current, times, voltage_step_times): """ Set values to 0 where they lie less than ``removal_time`` after a change in voltage. diff --git a/tests/test_herg_qc.py b/tests/test_herg_qc.py index f6217472..e2a4da54 100755 --- a/tests/test_herg_qc.py +++ b/tests/test_herg_qc.py @@ -645,7 +645,7 @@ def test_qc5(self): recording2 = np.asarray( [0, 0.1] * 100 + [10 * i] * 500 ) - pass_, d_max_diff = hergqc.qc5(recording1, recording2) + pass_, d_max_diff = hergqc.qc5(recording1, recording2, (0, -1)) self.assertAlmostEqual( d_max_diff, ex_d_max_diff, From b3a93cbe82c998793a03d0d8540ce027ee2aab4b Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Tue, 4 Nov 2025 17:25:06 +0000 Subject: [PATCH 06/13] Made filter_capacitive_spikes accept 1-d arrays again. --- pcpostprocess/hergQC.py | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index 3fbd3f48..a5e7b694 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -513,14 +513,13 @@ def filter_capacitive_spikes(self, current, times, voltage_step_times): """ Set values to 0 where they lie less than ``removal_time`` after a change in voltage. - @param current: The observed current + @param current: The observed current, as a 1 or multi-dimensional array. + If a multi-dimensional array is used, repeats must be + on the first axis, and time series values on the 2nd. @param times: the times at which the current was observed @param voltage_step_times: the times at which there are discontinuities in Vcmd @returns the ``current`` with some samples set to 0 """ - if len(current.shape) != 2: - raise ValueError('Current must have 2 dimensions (sweep, values)') - voltage_step_ends = np.append(voltage_step_times[1:], np.inf) for tstart, tend in zip(voltage_step_times, voltage_step_ends): win_end = tstart + self.removal_time From 50360f6d22f3befcb4b81ae958eb40dc68bb6d03 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Tue, 4 Nov 2025 18:34:22 +0000 Subject: [PATCH 07/13] Made filter_capacitive_spikes accept 1-d arrays again. --- pcpostprocess/hergQC.py | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index a5e7b694..9f4f9311 100644 --- a/pcpostprocess/hergQC.py +++ b/pcpostprocess/hergQC.py @@ -511,10 +511,10 @@ def qc6(self, recording1, win, label=None): def filter_capacitive_spikes(self, current, times, voltage_step_times): """ - Set values to 0 where they lie less than ``removal_time`` after a change in voltage. + Set currents to 0 where they lie less than ``removal_time`` after a change in voltage. - @param current: The observed current, as a 1 or multi-dimensional array. - If a multi-dimensional array is used, repeats must be + @param current: The observed current, as a 1 or 2-dimensional array. + If a 2-dimensional array is used, repeats must be on the first axis, and time series values on the 2nd. @param times: the times at which the current was observed @param voltage_step_times: the times at which there are discontinuities in Vcmd @@ -529,6 +529,11 @@ def filter_capacitive_spikes(self, current, times, voltage_step_times): if i_end == 0: break - current[:, i_start:i_end] = 0 + if len(current.shape) == 2: + current[:, i_start:i_end] = 0 + elif len(current.shape) == 1: + current[i_start:i_end] = 0 + else: + raise ValueError('Current array must be 1 or 2-dimensional') return current From 9d22edbe83e0e9609545915b86586d9652ed62b9 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 5 Nov 2025 10:46:26 +0000 Subject: [PATCH 08/13] Updated run_herg_qc to not use unnecessary n_sweeps constructor arg --- pcpostprocess/scripts/run_herg_qc.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/pcpostprocess/scripts/run_herg_qc.py b/pcpostprocess/scripts/run_herg_qc.py index 530a55a5..ae2bc9d9 100644 --- a/pcpostprocess/scripts/run_herg_qc.py +++ b/pcpostprocess/scripts/run_herg_qc.py @@ -693,8 +693,7 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args): os.makedirs(plot_dir) hergqc = hERGQC(sampling_rate=before_trace.sampling_rate, - plot_dir=plot_dir, - n_sweeps=before_trace.NofSweeps) + plot_dir=plot_dir) times = before_trace.get_times() voltage = before_trace.get_voltage() From f305a42577821d12abe62b39f47c292e70acdf57 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 5 Nov 2025 11:26:03 +0000 Subject: [PATCH 09/13] Added voltage to hergqc construction in run_herg_qc call --- .gitignore | 1 + pcpostprocess/scripts/run_herg_qc.py | 1 + 2 files changed, 2 insertions(+) diff --git a/.gitignore b/.gitignore index f1ff992d..c02b1046 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ /tests/test_data /test_output +/output *__pycache__* *.egg-info *.DS_Store diff --git a/pcpostprocess/scripts/run_herg_qc.py b/pcpostprocess/scripts/run_herg_qc.py index ae2bc9d9..1d0f1197 100644 --- a/pcpostprocess/scripts/run_herg_qc.py +++ b/pcpostprocess/scripts/run_herg_qc.py @@ -693,6 +693,7 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args): os.makedirs(plot_dir) hergqc = hERGQC(sampling_rate=before_trace.sampling_rate, + voltage=voltage, plot_dir=plot_dir) times = before_trace.get_times() From e6b480d98d4e3641f0fb008dfb6de793736bfa72 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 5 Nov 2025 13:29:46 +0000 Subject: [PATCH 10/13] Added new generated version file to gitignore --- .gitignore | 1 + 1 file changed, 1 insertion(+) diff --git a/.gitignore b/.gitignore index c02b1046..0a39bdfa 100644 --- a/.gitignore +++ b/.gitignore @@ -1,3 +1,4 @@ +pcpostprocess/_version.py /tests/test_data /test_output /output From 4ce60be4eecf4683f0f1a67e88eabfdca713a903 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Wed, 5 Nov 2025 13:29:55 +0000 Subject: [PATCH 11/13] Style fix --- tests/test_herg_qc.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_herg_qc.py b/tests/test_herg_qc.py index 95a2333a..e2a4da54 100755 --- a/tests/test_herg_qc.py +++ b/tests/test_herg_qc.py @@ -69,7 +69,7 @@ def create_hergqc(self, plot_dir=None): os.makedirs(plot_dir, exist_ok=True) else: plot_dir = None - + return hERGQC( sampling_rate=self.sampling_rate, plot_dir=plot_dir, From 430fd64f49329bbeee08cc27713657e5826e87c7 Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Tue, 11 Nov 2025 18:12:50 +0000 Subject: [PATCH 12/13] Clarified changelog use in contributing.md --- CONTRIBUTING.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index a8770f83..15e1b12c 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -149,6 +149,9 @@ For example, the first entry in our Changelog was: - [#104](https://github.com/CardiacModelling/pcpostprocess/pull/104) Added a CHANGELOG.md and CONTRIBUTING.md ``` +Changelog entries are intended for _users_, and so should focus on changes to the public API or command line interface. +Changes to internals are less likely to need reporting. + ## Packaging This project uses a minimal [`setup.py`](./setup.py), and instead uses [pyproject.toml](./pyproject.toml). From 7300038a25c11639866534915cf7ee0b1c788b5a Mon Sep 17 00:00:00 2001 From: Michael Clerx Date: Tue, 11 Nov 2025 18:12:56 +0000 Subject: [PATCH 13/13] Updated changelog --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 7d3fe562..685a6de4 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -4,10 +4,13 @@ This page lists the main changes made to pcpostprocess in each release. ## Unreleased - Added + - [#81](https://github.com/CardiacModelling/pcpostprocess/pull/81) Added docstrings to the `hERGQC` class. - [#104](https://github.com/CardiacModelling/pcpostprocess/pull/104) Added a CHANGELOG.md and CONTRIBUTING.md - Changed + - [#81](https://github.com/CardiacModelling/pcpostprocess/pull/81) Changed the constructor arguments for `hERGQC`. - Deprecated - Removed + - [#81](https://github.com/CardiacModelling/pcpostprocess/pull/81) Removed `hERGQC.plot_dir`, `hERGQC.set_trace` and `hERGQC.set_debug`. - Fixed ## [0.2.0] - 2025-11-11