Skip to content

Commit

Permalink
add peptide plotting example
Browse files Browse the repository at this point in the history
  • Loading branch information
Jhsmit committed Oct 14, 2024
1 parent 40ff584 commit a46afaf
Showing 1 changed file with 40 additions and 4 deletions.
44 changes: 40 additions & 4 deletions templates/05_SecB_fit_dG.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,17 +3,20 @@
# %%
from pathlib import Path

import numpy as np
import proplot as pplt
import yaml

from pyhdx.datasets import HDXDataSet
from pyhdx.models import HDXMeasurement
from pyhdx.fitting import fit_gibbs_global, fit_rates_weighted_average
from pyhdx.fileIO import csv_to_dataframe, save_fitresult
from pyhdx.fitting import fit_gibbs_global, fit_rates_weighted_average
from pyhdx.local_cluster import default_client
import yaml
from pyhdx.models import HDXMeasurement

# %%

guess = False # Set to True to redo initial guesses
fit_kwargs = {"epochs": 10000, "lr": 1e4, "stop_loss": 1e-6}
fit_kwargs = {"epochs": 200000, "lr": 1e4, "stop_loss": 1e-6}

# %%
current_dir = Path(__file__).parent
Expand Down Expand Up @@ -55,3 +58,36 @@
fr_torch.to_file(output_dir / "SecB_fit_result.csv", fmt="csv")

save_fitresult(output_dir / "SecB_fit", fr_torch)

# %%
# evaluate the fit, plot measured and fitted D-uptake for a few peptides
NUM_EVAL_POINTS = 100
all_timepoints = np.concatenate([hdxm.timepoints for hdxm in fr_torch.hdxm_set])

elem = all_timepoints[np.nonzero(all_timepoints)]
start = np.log10(elem.min())
end = np.log10(elem.max())
pad = (end - start) * 0.1
time = np.logspace(start - pad, end + pad, num=NUM_EVAL_POINTS, endpoint=True)


# %%
# evaluate the fitted model at timepoints
d_calc = fr_torch(time) # Ns x Np x Nt
d_calc
# %%
# choose which sample to plot
# here we plot the first one (= SecB tetramer)
sample_idx = 0
hdxm = fr_torch.hdxm_set.hdxm_list[sample_idx]
d_calc_s = d_calc[sample_idx]

# %%
# make a subplot grid to plot the first 24 peptides
fig, axes = pplt.subplots(ncols=4, nrows=6)
for i, ax in enumerate(axes):
ax.plot(time, d_calc_s[i], color="r")
ax.scatter(hdxm.timepoints, hdxm.d_exp.iloc[i], color="k")

axes.format(xscale="log", xlabel="Time (s)", ylabel="Corrected D-uptake", ylim=(0, None))
# %%

0 comments on commit a46afaf

Please sign in to comment.