Skip to content

Commit

Permalink
Merge pull request #24 from CardiacModelling/js_dev_plots
Browse files Browse the repository at this point in the history 
Apply QC to leak-corrected traces. Fix subtraction plots.
  • Loading branch information
@joeyshuttleworth
joeyshuttleworth authored Aug 8, 2024
2 parents 40ada3d + 5e6ba3e commit e520c88
Show file tree
Hide file tree
Showing 5 changed files with 66 additions and 48 deletions.
2 changes: 1 addition & 1 deletion .github/workflows/pytest.yml
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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/
Expand Down
2 changes: 1 addition & 1 deletion pcpostprocess/hergQC.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
4 changes: 2 additions & 2 deletions pcpostprocess/leak_correct.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
53 changes: 28 additions & 25 deletions pcpostprocess/scripts/run_herg_qc.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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]
Expand All @@ -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]
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down
53 changes: 34 additions & 19 deletions pcpostprocess/subtraction_plots.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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


Expand All @@ -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 = []
Expand All @@ -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():
Expand All @@ -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()
Expand All @@ -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'))

Expand All @@ -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)

0 comments on commit e520c88

Please sign in to comment.