Apply QC to leak-corrected traces. Fix subtraction plots.

.github/workflows/pytest.yml

@@ -27,7 +27,7 @@ jobs:
         eval `ssh-agent -s`
         ssh-add - <<< '${{ secrets.syncropatch_export_key }}'
         python -m pip install --upgrade pip
-        python -m pip install .[test]
+        python -m pip install -e .[test]
     - name: Extract test data
       run: |
         wget https://cardiac.nottingham.ac.uk/syncropatch_export/test_data.tar.xz -P tests/

pcpostprocess/hergQC.py

@@ -18,7 +18,7 @@ class hERGQC(object):
     no_QC = len(QCnames)
 
     def __init__(self, sampling_rate=5, plot_dir=None, voltage=np.array([]),
-                 n_sweeps=None, removal_time=2):
+                 n_sweeps=None, removal_time=5):
         # TODO docstring
 
         if plot_dir is not None:

pcpostprocess/leak_correct.py

@@ -168,9 +168,9 @@ def fit_linear_leak(current, voltage, times, ramp_start_index, ramp_end_index,
                  I_leak[ramp_start_index:ramp_end_index+1], '--')
 
         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_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}$')
+                 linestyle='--', 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

@@ -229,10 +229,6 @@ def main():
             fout.write(well)
             fout.write('\n')
 
-    logfile = os.path.join(savedir, 'table-%s.txt' % saveID)
-    with open(logfile, 'a') as f:
-        f.write('\\end{table}\n')
-
     # Export all protocols
     savenames, readnames, times_list = [], [], []
     for protocol in res_dict:
@@ -670,7 +666,8 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args):
             row_dict['R_leftover'] =\
                 np.sqrt(np.sum((after_corrected)**2)/(np.sum(before_corrected**2)))
 
-            row_dict['QC.R_leftover'] = row_dict['R_leftover'] < 0.5
+            QC_R_leftover = np.all(row_dict['R_leftover'] < 0.5)
+            row_dict['QC.R_leftover'] = QC_R_leftover
 
             row_dict['E_rev'] = E_rev
             row_dict['E_rev_before'] = E_rev_before
@@ -735,7 +732,7 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args):
                                                                           'debug',
                                                                           '-120mV time constant',
                                                                           f"{savename}-{well}-sweep"
-                                                                          "{sweep}-time-constant-fit.png"))
+                                                                          "{sweep}-time-constant-fit.pdf"))
 
             row_dict['-120mV decay time constant 1'] = res[0][0]
             row_dict['-120mV decay time constant 2'] = res[0][1]
@@ -754,13 +751,6 @@ def extract_protocol(readname, savename, time_strs, selected_wells, args):
     before_current_all = before_trace.get_trace_sweeps()
     after_current_all = after_trace.get_trace_sweeps()
 
-    # Convert everything to nA...
-    before_current_all = {key: value * 1e-3 for key, value in before_current_all.items()}
-    after_current_all = {key: value * 1e-3 for key, value in after_current_all.items()}
-
-    before_leak_current_dict = {key: value * 1e-3 for key, value in before_leak_current_dict.items()}
-    after_leak_current_dict = {key: value * 1e-3 for key, value in after_leak_current_dict.items()}
-
     for well in selected_wells:
         before_current = before_current_all[well]
         after_current = after_current_all[well]
@@ -865,12 +855,11 @@ def run_qc_for_protocol(readname, savename, time_strs, args):
 
         # Check if any cell first!
         if (None in qc_before[well][0]) or (None in qc_after[well][0]):
-            # no_cell = True
+            no_cell = True
             continue
 
         else:
-            # no_cell = False
-            pass
+            no_cell = False
 
         nsweeps = before_trace.NofSweeps
         assert after_trace.NofSweeps == nsweeps
@@ -924,19 +913,21 @@ def run_qc_for_protocol(readname, savename, time_strs, args):
         logging.info(f"{well} {savename}\n----------")
         logging.info(f"sampling_rate is {sampling_rate}")
 
-        voltage_steps = [tstart
+        voltage_steps = [tend
                          for tstart, tend, vstart, vend in
                          voltage_protocol.get_all_sections() if vend == vstart]
 
         # Run QC with raw currents
-        selected, QC = hergqc.run_qc(voltage_steps, times,
-                                     before_currents,
-                                     after_currents,
-                                     np.array(qc_before[well])[0, :],
-                                     np.array(qc_after[well])[0, :], nsweeps)
+        _, QC = hergqc.run_qc(voltage_steps, times,
+                              before_currents_corrected,
+                              after_currents_corrected,
+                              np.array(qc_before[well])[0, :],
+                              np.array(qc_after[well])[0, :], nsweeps)
 
+        QC = list(QC)
         df_rows.append([well] + list(QC))
 
+        selected = np.all(QC) and not no_cell
         if selected:
             selected_wells.append(well)
 
@@ -1044,9 +1035,20 @@ def qc3_bookend(readname, savename, time_strs, args):
         last_before_current = last_before_current_dict[well][-1, :]
         last_after_current = last_after_current_dict[well][-1, :]
 
+        output_directory = os.path.join(args.output_dir, "leak_correction")
+        save_fname = f"{well}_{savename}_before0.pdf"
+
+        #  Plot subtraction
+        get_leak_corrected(first_before_current,
+                           voltage, times,
+                           *ramp_bounds,
+                           save_fname=save_fname,
+                           output_dir=output_directory)
+
         before_traces_first[well] = get_leak_corrected(first_before_current,
                                                        voltage, times,
                                                        *ramp_bounds)
+
         before_traces_last[well] = get_leak_corrected(last_before_current,
                                                       voltage, times,
                                                       *ramp_bounds)
@@ -1197,11 +1199,12 @@ def fit_func(x, args=None):
         for ax in axs:
             ax.spines[['top', 'right']].set_visible(False)
             ax.set_ylabel(r'$I_\mathrm{obs}$ (pA)')
-            ax.set_xlabel(r'$t$ (ms)')
+
+        axs[-1].set_xlabel(r'$t$ (ms)')
 
         protocol_ax, fit_ax = axs
-        protocol_ax.set_title('a', fontweight='bold')
-        fit_ax.set_title('b', fontweight='bold')
+        protocol_ax.set_title('a', fontweight='bold', loc='left')
+        fit_ax.set_title('b', fontweight='bold', loc='left')
         fit_ax.plot(peak_time, peak_current, marker='x', color='red')
 
         a, b, c, d = res1.x

pcpostprocess/subtraction_plots.py

@@ -26,6 +26,13 @@ def setup_subtraction_grid(fig, nsweeps):
     # Long axis for protocol on the bottom (full width)
     long_protocol_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],
+                       'abcdefghijklm'):
+        ax.spines[['top', 'right']].set_visible(False)
+        ax.set_title(cap, loc='left', fontweight='bold')
+
     return protocol_axs, before_axs, after_axs, corrected_axs, subtracted_ax, long_protocol_ax
 
 
@@ -40,9 +47,9 @@ def do_subtraction_plot(fig, times, sweeps, before_currents, after_currents,
         subtracted_ax, long_protocol_ax = axs
 
     for ax in protocol_axs:
-        ax.plot(times, voltages, color='black')
+        ax.plot(times*1e-3, voltages, color='black')
         ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$V_\mathrm{command}$ (mV)')
+        ax.set_ylabel(r'$V_\mathrm{cmd}$ (mV)')
 
     all_leak_params_before = []
     all_leak_params_after = []
@@ -62,17 +69,22 @@ def do_subtraction_plot(fig, times, sweeps, before_currents, after_currents,
                                   np.nan)
     for i, sweep in enumerate(sweeps):
 
-        gleak, Eleak = all_leak_params_before[i]
+        b0, b1 = all_leak_params_before[i]
+        gleak = b1
+        Eleak = -b1/b0
         before_leak_currents[i, :] = gleak * (voltages - Eleak)
 
-        gleak, Eleak = all_leak_params_after[i]
+        b0, b1 = all_leak_params_after[i]
+        gleak = b1
+        Eleak = -b1/b0
+
         after_leak_currents[i, :] = gleak * (voltages - Eleak)
 
     for i, (sweep, ax) in enumerate(zip(sweeps, before_axs)):
         gleak, Eleak = all_leak_params_before[i]
-        ax.plot(times, before_currents[i, :], label=f"pre-drug raw, sweep {sweep}")
-        ax.plot(times, before_leak_currents[i, :],
-                label=r'$I_\mathrm{leak}$.' f"g={gleak:1E}, E={Eleak:.1e}")
+        ax.plot(times*1e-3, before_currents[i, :], label=f"pre-drug raw, sweep {sweep}")
+        ax.plot(times*1e-3, before_leak_currents[i, :],
+                label=r'$I_\mathrm{L}$.' f"g={gleak:1E}, E={Eleak:.1e}")
         # ax.legend()
 
         if ax.get_legend():
@@ -84,9 +96,9 @@ def do_subtraction_plot(fig, times, sweeps, before_currents, after_currents,
 
     for i, (sweep, ax) in enumerate(zip(sweeps, after_axs)):
         gleak, Eleak = all_leak_params_before[i]
-        ax.plot(times, after_currents[i, :], label=f"post-drug raw, sweep {sweep}")
-        ax.plot(times, after_leak_currents[i, :],
-                label=r"$I_\mathrm{leak}$." f"g={gleak:1E}, E={Eleak:.1e}")
+        ax.plot(times*1e-3, after_currents[i, :], label=f"post-drug raw, sweep {sweep}")
+        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()
@@ -98,12 +110,12 @@ def do_subtraction_plot(fig, times, sweeps, before_currents, after_currents,
     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, corrected_before_currents,
-                label=f"leak corrected before drug trace, sweep {sweep}")
-        ax.plot(times, corrected_after_currents,
-                label=f"leak corrected after drug trace, sweep {sweep}")
+        ax.plot(times*1e-3, corrected_before_currents,
+                label=f"leak-corrected pre-drug trace, sweep {sweep}")
+        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.set_ylabel(r'leak-corrected traces')
         # ax.tick_params(axis='y', rotation=90)
         # ax.yaxis.set_major_formatter(mtick.FormatStrFormatter('%.1e'))
 
@@ -118,13 +130,16 @@ def do_subtraction_plot(fig, times, sweeps, before_currents, after_currents,
 
         subtracted_currents = before_currents[i, :] - before_leak_currents[i, :] - \
             (after_currents[i, :] - after_leak_currents[i, :])
-        ax.plot(times, subtracted_currents, label=f"sweep {sweep}")
+        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)
 
-    ax.set_ylabel(r'$I_\mathrm{obs} - I_\mathrm{l}$ (mV)')
+    ax.set_ylabel(r'$I_\mathrm{obs} - I_\mathrm{L}$ (mV)')
     ax.set_xlabel('$t$ (s)')
 
-    long_protocol_ax.plot(times, voltages, color='black')
+    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{command}$ (mV)')
+    long_protocol_ax.set_ylabel(r'$V_\mathrm{cmd}$ (mV)')
     long_protocol_ax.tick_params(axis='y', rotation=90)