Merge branch 'main' into hergQC-docstring

Author: MichaelClerx

README.md

@@ -5,7 +5,7 @@ This repository contains a python package and scripts for handling time-series d
 The package has been tested with data from a SyncroPatch 384, but may be adapted to work with data in other formats.
 It can also be used to perform quality control (QC) as described in [Lei et al. (2019)](https://doi.org/10.1016%2Fj.bpj.2019.07.029).
 
-This package is tested on Ubuntu with Python 3.8, 3.9, 3.10, 3.11, 3.12 and 3.13.
+This package is tested on Ubuntu with Python 3.10, 3.11, 3.12, 3.13 and 3.14.
 
 ## Getting Started
 

pcpostprocess/leak_correct.py

@@ -149,7 +149,7 @@ def fit_linear_leak(current, voltage, times, ramp_start_index, ramp_end_index,
         ax4.plot(times, I_obs, label=r'$I_\mathrm{obs}$')
         ax4.plot(times, I_leak, linestyle='--', label=r'$I_\mathrm{L}$')
         ax4.plot(times, I_obs - I_leak,
-                 linestyle='--', alpha=0.5, label=r'$I_\mathrm{obs} - I_\mathrm{L}$')
+                 alpha=0.5, label=r'$I_\mathrm{obs} - I_\mathrm{L}$')
         ax4.legend(frameon=False)
 
         if not os.path.exists(output_dir):

pcpostprocess/scripts/run_herg_qc.py

@@ -51,7 +51,7 @@ def main():
     parser.add_argument('--reversal_spread_threshold', type=float, default=10)
     parser.add_argument('--export_failed', action='store_true')
     parser.add_argument('--selection_file')
-    parser.add_argument('--figsize', nargs=2, type=int, default=[5, 8])
+    parser.add_argument('--figsize', nargs=2, type=int, default=[16, 18])
     parser.add_argument('--debug', action='store_true')
     parser.add_argument('--log_level', default='INFO')
     parser.add_argument('--Erev', default=-90.71, type=float)
@@ -517,7 +517,6 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args):
     #  Find start of leak section
     desc = voltage_protocol.get_all_sections()
     ramp_bounds = detect_ramp_bounds(times, desc)
-    tstart, tend = ramp_bounds
 
     nsweeps_before = before_trace.NofSweeps = 2
     nsweeps_after = after_trace.NofSweeps = 2

pcpostprocess/subtraction_plots.py

@@ -5,9 +5,9 @@
 
 
 def setup_subtraction_grid(fig, nsweeps):
-    # Use 5 x 2 grid when there are 2 sweeps
-    gs = GridSpec(6, nsweeps, figure=fig)
-
+    # Use 6 x 2 grid when there are 2 sweeps
+    gs = GridSpec(6, nsweeps, figure=fig, height_ratios=[0.7, 2, 2, 2, 0.7, 3])
+    # gs.subplots(sharey='row')
     # plot protocol at the top
     protocol_axs = [fig.add_subplot(gs[0, i]) for i in range(nsweeps)]
 
@@ -20,18 +20,20 @@ def setup_subtraction_grid(fig, nsweeps):
     # Leak corrected traces
     corrected_axs = [fig.add_subplot(gs[3, i]) for i in range(nsweeps)]
 
-    # Subtracted traces on one axis
-    subtracted_ax = fig.add_subplot(gs[4, :])
-
     # Long axis for protocol on the bottom (full width)
-    long_protocol_ax = fig.add_subplot(gs[5, :])
+    long_protocol_ax = fig.add_subplot(gs[4, :])
+
+    # Subtracted traces on one axis
+    subtracted_ax = fig.add_subplot(gs[5, :])
 
     for ax, cap in zip(list(protocol_axs) + list(before_axs)
                        + list(after_axs) + list(corrected_axs)
-                       + [subtracted_ax] + [long_protocol_ax],
+                       + [long_protocol_ax] + [subtracted_ax],
                        'abcdefghijklm'):
         ax.spines[['top', 'right']].set_visible(False)
         ax.set_title(cap, loc='left', fontweight='bold')
+        if cap != 'a':
+            ax.sharex(protocol_axs[0])
 
     return protocol_axs, before_axs, after_axs, corrected_axs, subtracted_ax, long_protocol_ax
 
@@ -46,100 +48,91 @@ def do_subtraction_plot(fig, times, sweeps, before_currents, after_currents,
     protocol_axs, before_axs, after_axs, corrected_axs, \
         subtracted_ax, long_protocol_ax = axs
 
-    for ax in protocol_axs:
+    for i, ax in enumerate(protocol_axs):
         ax.plot(times*1e-3, voltages, color='black')
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$V_\mathrm{cmd}$ (mV)')
+        ax.set_title(f'Well {well}, sweep {sweeps[i]}', fontweight='bold')
+        ax.tick_params(axis='x', labelbottom=False)
+
+    protocol_axs[0].set_ylabel(r'$V_\mathrm{cmd}$ (mV)', fontsize=16)
 
     all_leak_params_before = []
     all_leak_params_after = []
-    for i in range(len(sweeps)):
-        before_params, _ = fit_linear_leak(before_currents, voltages, times,
-                                           *ramp_bounds)
-        all_leak_params_before.append(before_params)
 
-        after_params, _ = fit_linear_leak(after_currents, voltages, times,
-                                          *ramp_bounds)
-        all_leak_params_after.append(after_params)
+    alpha_of_zero = 0.2
+    style_of_zero = '-'
+    range_of_zero = [times[0]*1e-3, times[-1]*1e-3]
 
     # Compute and store leak currents
     before_leak_currents = np.full((nsweeps, voltages.shape[0]),
                                    np.nan)
     after_leak_currents = np.full((nsweeps, voltages.shape[0]),
                                   np.nan)
-    for i, sweep in enumerate(sweeps):
-
-        b0, b1 = all_leak_params_before[i]
-        gleak = b1
-        Eleak = -b1/b0
-        before_leak_currents[i, :] = gleak * (voltages - Eleak)
 
-        b0, b1 = all_leak_params_after[i]
-        gleak = b1
-        Eleak = -b1/b0
+    for i in range(len(sweeps)):
+        before_params, before_leak_current = fit_linear_leak(before_currents[i], voltages, times,
+                                                             *ramp_bounds)
+        before_leak_currents[i, :] = before_leak_current
+        all_leak_params_before.append(before_params)
 
-        after_leak_currents[i, :] = gleak * (voltages - Eleak)
+        after_params, after_leak_current = fit_linear_leak(after_currents[i], voltages, times,
+                                                           *ramp_bounds)
+        all_leak_params_after.append(after_params)
+        after_leak_currents[i, :] = after_leak_current
 
     for i, (sweep, ax) in enumerate(zip(sweeps, before_axs)):
-        gleak, Eleak = all_leak_params_before[i]
-        ax.plot(times*1e-3, before_currents[i, :], label=f"pre-drug raw, sweep {sweep}")
+        b0, b1 = all_leak_params_before[i]
+        gleak = b1
+        Eleak = -b0/b1
+
+        ax.plot(times*1e-3, before_currents[i, :], label="Pre-drug raw")
         ax.plot(times*1e-3, before_leak_currents[i, :],
-                label=r'$I_\mathrm{L}$.' f"g={gleak:1E}, E={Eleak:.1e}")
-        # ax.legend()
+                label=f"Fitted leak g={gleak:7.5g}, E={Eleak:7.4g} mV")
+        ax.plot(range_of_zero, [0, 0], color='black', linestyle=style_of_zero, alpha=alpha_of_zero)
+        ax.legend()
+        ax.tick_params(axis='x', labelbottom=False)
 
-        if ax.get_legend():
-            ax.get_legend().remove()
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'pre-drug trace')
-        # ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
-        # ax.tick_params(axis='y', rotation=90)
+    before_axs[0].set_ylabel(r'Pre-drug trace', fontsize=16)
 
     for i, (sweep, ax) in enumerate(zip(sweeps, after_axs)):
-        gleak, Eleak = all_leak_params_before[i]
-        ax.plot(times*1e-3, after_currents[i, :], label=f"post-drug raw, sweep {sweep}")
+        b0, b1 = all_leak_params_after[i]
+        gleak = b1
+        Eleak = -b0/b1
+
+        ax.plot(times*1e-3, after_currents[i, :], label="Post-drug raw")
         ax.plot(times*1e-3, after_leak_currents[i, :],
-                label=r"$I_\mathrm{L}$." f"g={gleak:1E}, E={Eleak:.1e}")
-        # ax.legend()
-        if ax.get_legend():
-            ax.get_legend().remove()
-        ax.set_xlabel('$t$ (s)')
-        ax.set_ylabel(r'post-drug trace')
-        # ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
-        # ax.tick_params(axis='y', rotation=90)
+                label=f"Fitted leak g={gleak:7.5g}, E={Eleak:7.4g} mV")
+        ax.plot(range_of_zero, [0, 0], color='black', linestyle=style_of_zero, alpha=alpha_of_zero)
+        ax.legend()
+        ax.tick_params(axis='x', labelbottom=False)
+    after_axs[0].set_ylabel(r'Post-drug trace', fontsize=16)
 
     for i, (sweep, ax) in enumerate(zip(sweeps, corrected_axs)):
         corrected_before_currents = before_currents[i, :] - before_leak_currents[i, :]
         corrected_after_currents = after_currents[i, :] - after_leak_currents[i, :]
         ax.plot(times*1e-3, corrected_before_currents,
-                label=f"leak-corrected pre-drug trace, sweep {sweep}")
+                label="Leak-corrected pre-drug trace")
         ax.plot(times*1e-3, corrected_after_currents,
-                label=f"leak-corrected post-drug trace, sweep {sweep}")
-        ax.set_xlabel(r'$t$ (s)')
-        ax.set_ylabel(r'leak-corrected traces')
-        # ax.tick_params(axis='y', rotation=90)
-        # ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
+                label="Leak-corrected post-drug trace")
+        ax.plot(range_of_zero, [0, 0], color='black', linestyle=style_of_zero, alpha=alpha_of_zero)
+        ax.set_xlabel(r'Time (s)')
+        ax.legend()
+    corrected_axs[0].set_ylabel(r'Leak-corrected traces', fontsize=16)
 
-    ax = subtracted_ax
     for i, sweep in enumerate(sweeps):
-        before_trace = before_currents[i, :].flatten()
-        after_trace = after_currents[i, :].flatten()
-        before_params, before_leak = fit_linear_leak(before_trace, voltages, times,
-                                                     *ramp_bounds)
-        after_params, after_leak = fit_linear_leak(after_trace, voltages, times,
-                                                   *ramp_bounds)
-
         subtracted_currents = before_currents[i, :] - before_leak_currents[i, :] - \
             (after_currents[i, :] - after_leak_currents[i, :])
-        ax.plot(times*1e-3, subtracted_currents, label=f"sweep {sweep}", alpha=.5)
 
-        #  Cycle to next colour
-        ax.plot([np.nan], [np.nan], label=f"sweep {sweep}", alpha=.5)
+        subtracted_ax.plot(times*1e-3, subtracted_currents, label=f"sweep {sweep}", alpha=.5)
 
-    ax.set_ylabel(r'$I_\mathrm{obs} - I_\mathrm{L}$ (mV)')
-    ax.set_xlabel('$t$ (s)')
+    subtracted_ax.legend()
+    subtracted_ax.plot(range_of_zero, [0, 0], color='black',
+                       linestyle=style_of_zero, alpha=alpha_of_zero)
+    subtracted_ax.set_ylabel('Final subtracted traces', fontsize=16)
+    subtracted_ax.set_xlabel('Time (s)', fontsize=16)
 
     long_protocol_ax.plot(times*1e-3, voltages, color='black')
-    long_protocol_ax.set_xlabel('time (s)')
-    long_protocol_ax.set_ylabel(r'$V_\mathrm{cmd}$ (mV)')
-    long_protocol_ax.tick_params(axis='y', rotation=90)
+    long_protocol_ax.set_ylabel(r'$V_\mathrm{cmd}$ (mV)', fontsize=16)
+    long_protocol_ax.tick_params(axis='x', labelbottom=False)
+
+    fig.align_ylabels()
 

setup.py

@@ -13,18 +13,18 @@
     # Module name (lowercase)
     name='pcpostprocess',
     version='0.0.1',
-    description='Export high-throughput patch-clamp data from the Nanion SyncroPatch',
+    description='Post-process patch clamp recordings with the staircase protocol',
     long_description=readme,
     long_description_content_type="text/markdown",
-    author='Frankie Patten-Elliot, Joseph Shuttleworth, Chon Lok Lei',
+    author='Frankie Patten-Elliot, Joseph Shuttleworth, Chon Lok Lei, Michael Clerx',
     author_email='joseph.shuttleworth@nottingham.ac.uk',
     maintainer='Joseph Shuttleworth',
     maintainer_email='joseph.shuttleworth@nottingham.ac.uk',
-    # url='https://github.com/CardiacModelling/markov-builder',
+    url='https://github.com/CardiacModelling/pcpostprocess',
     classifiers=[
-        "Programming Language :: Python :: 3",
-        "License :: OSI Approved :: BSD License",
-        "Operating System :: OS Independent",
+        'Programming Language :: Python :: 3',
+        'License :: OSI Approved :: BSD License',
+        'Operating System :: OS Independent',
     ],
 
     # Packages to include
@@ -35,7 +35,7 @@
     include_package_data=True,
 
     # Required Python version
-    python_requires='>=3.7',
+    python_requires='>=3.10',
 
     # List of dependencies
     install_requires=[
@@ -44,25 +44,23 @@
         'matplotlib>=3.4',
         'pandas>=1.3',
         'regex>=2023.12.25',
-        'openpyxl>=3.1.2',
-        'jinja2>=3.1.0',
-        'seaborn>=0.12.2'
+        'seaborn>=0.12.2',
+        'openpyxl>=3.1.2',      # Used via pandas (to create excel doc)
+        'jinja2>=3.1.0',        # Used via pandas (to create latex doc)
     ],
     extras_require={
         'test': [
             'pytest-cov>=2.10',     # For coverage checking
             'pytest>=4.6',          # For unit tests
             'flake8>=3',            # For code style checking
             'isort',
-            'mock>=3.0.5',         # For mocking command line args etc.
             'codecov>=2.1.3',
             'syncropatch_export @ git+https://github.com/CardiacModelling/syncropatch_export.git'
         ],
     },
     entry_points={
         'console_scripts': [
-            'pcpostprocess='
-            'pcpostprocess.scripts.__main__:main',
+            'pcpostprocess=pcpostprocess.scripts.__main__:main',
         ],
     },
 )

tests/test_leak_correct.py

@@ -5,6 +5,7 @@
 from syncropatch_export.trace import Trace
 
 from pcpostprocess import leak_correct
+from pcpostprocess.detect_ramp_bounds import detect_ramp_bounds
 
 
 class TestLeakCorrect(unittest.TestCase):
@@ -18,15 +19,20 @@ def setUp(self):
         if not os.path.exists(self.output_dir):
             os.makedirs(self.output_dir)
 
-        self.ramp_bounds = [1700, 2500]
         self.test_trace = Trace(test_data_dir, json_file)
 
         # get currents and QC from trace object
         self.currents = self.test_trace.get_all_traces(leakcorrect=False)
         self.currents['times'] = self.test_trace.get_times()
         self.currents['voltages'] = self.test_trace.get_voltage()
+
         self.QC = self.test_trace.get_onboard_QC_values()
 
+        # Find first times ahead of these times
+        voltage_protocol = self.test_trace.get_voltage_protocol().get_all_sections()
+        times = self.currents['times'].flatten()
+        self.ramp_bound_indices = detect_ramp_bounds(times, voltage_protocol, ramp_no=0)
+
     def test_plot_leak_fit(self):
         well = 'A01'
         sweep = 0
@@ -37,7 +43,7 @@ def test_plot_leak_fit(self):
         current = self.test_trace.get_trace_sweeps(sweeps=[sweep])[well][0, :]
 
         leak_correct.fit_linear_leak(current, voltage, times,
-                                     *self.ramp_bounds,
+                                     *self.ramp_bound_indices,
                                      output_dir=self.output_dir,
                                      save_fname=f"{well}_sweep{sweep}_leak_correction")
 
@@ -50,7 +56,8 @@ def test_get_leak_correct(self):
         times = trace.get_times()
 
         current = currents[well][sweep, :]
-        x = leak_correct.get_leak_corrected(current, voltage, times, *self.ramp_bounds)
+        x = leak_correct.get_leak_corrected(current, voltage, times,
+                                            *self.ramp_bound_indices)
         self.assertEqual(x.shape, (30784,))
 
 

tests/test_subtraction_plots.py

@@ -11,23 +11,26 @@
 
 class TestSubtractionPlots(unittest.TestCase):
     def setUp(self):
-        test_data_dir = os.path.join('tests', 'test_data', '13112023_MW2_FF',
-                                     "staircaseramp (2)_2kHz_15.01.07")
-        json_file = "staircaseramp (2)_2kHz_15.01.07.json"
+        test_data_dir_before = os.path.join('tests', 'test_data', '13112023_MW2_FF',
+                                            "staircaseramp (2)_2kHz_15.01.07")
+
+        test_data_dir_after = os.path.join('tests', 'test_data', '13112023_MW2_FF',
+                                           "staircaseramp (2)_2kHz_15.11.33")
+
+        json_file_before = "staircaseramp (2)_2kHz_15.01.07.json"
+        json_file_after = "staircaseramp (2)_2kHz_15.11.33.json"
 
         self.output_dir = os.path.join('test_output', 'test_trace_class')
 
         if not os.path.exists(self.output_dir):
             os.makedirs(self.output_dir)
 
-        self.ramp_bounds = [1700, 2500]
-
         # Use identical traces for purpose of the test
-        self.before_trace = Trace(test_data_dir, json_file)
-        self.after_trace = Trace(test_data_dir, json_file)
+        self.before_trace = Trace(test_data_dir_before, json_file_before)
+        self.after_trace = Trace(test_data_dir_after, json_file_after)
 
     def test_do_subtraction_plot(self):
-        fig = plt.figure(layout='constrained')
+        fig = plt.figure(figsize=(16, 18), layout='constrained')
         times = self.before_trace.get_times()
 
         well = 'A01'
@@ -44,6 +47,8 @@ def test_do_subtraction_plot(self):
         do_subtraction_plot(fig, times, sweeps, before_current, after_current,
                             voltages, ramp_bounds, well=well)
 
+        fig.savefig(os.path.join(self.output_dir, f"subtraction_plot_{well}"))
+
 
 if __name__ == "__main__":
     pass