diff --git a/.github/workflows/pytest.yml b/.github/workflows/pytest.yml index 76f5fdd..5878ab6 100644 --- a/.github/workflows/pytest.yml +++ b/.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/ diff --git a/pcpostprocess/hergQC.py b/pcpostprocess/hergQC.py index 049e010..06d8c17 100644 --- a/pcpostprocess/hergQC.py +++ b/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: diff --git a/pcpostprocess/leak_correct.py b/pcpostprocess/leak_correct.py index 39dad18..ac542c3 100644 --- a/pcpostprocess/leak_correct.py +++ b/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): diff --git a/pcpostprocess/scripts/run_herg_qc.py b/pcpostprocess/scripts/run_herg_qc.py index 237c08e..496c9f2 100644 --- a/pcpostprocess/scripts/run_herg_qc.py +++ b/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 diff --git a/pcpostprocess/subtraction_plots.py b/pcpostprocess/subtraction_plots.py index 941c658..2f90141 100644 --- a/pcpostprocess/subtraction_plots.py +++ b/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)