def plot_example_hit(chrom, start, end, prof_center_size=1000, score_center_size=50, save_path=None):
"""
For a given putative motif hit, plots the ChIP-nexus and ChIP-seq
profiles, as well as the importance scores for the task around the
hit.
"""
mid = (start + end) // 2
prof_start = mid - (prof_center_size // 2)
prof_end = prof_start + prof_center_size
score_start = mid - (score_center_size // 2)
score_end = score_start + score_center_size
with h5py.File(true_profiles_hdf5_path, "r") as f:
true_profs = f[chrom][prof_start:prof_end][:, 0]
with h5py.File(pred_profiles_hdf5_path, "r") as f:
# Need to use the coordinates of the profiles themselves
prof_len = f["predictions"]["log_pred_profs"].shape[2]
prof_coords_chrom = f["coords"]["coords_chrom"][:].astype(str)
prof_coords_start = f["coords"]["coords_start"][:]
prof_coords_end = f["coords"]["coords_end"][:]
mid = (prof_coords_start + prof_coords_end) // 2
prof_coords_start = mid - (prof_len // 2)
prof_coords_end = prof_coords_start + prof_len
match_inds = np.where(
(prof_coords_chrom == chrom) &
(prof_coords_start <= prof_start) &
(prof_coords_end >= prof_end)
)[0]
if not match_inds.size:
print("Warning: did not find sufficiently large prediction track for %s:%d-%s" % (chrom, prof_start, prof_end))
return
match_ind = match_inds[0]
coord_start, coord_end = prof_coords_start[match_ind], prof_coords_end[match_ind]
cut_start = prof_start - coord_start
cut_end = cut_start + prof_center_size
pred_profs = np.exp(f["predictions"]["log_pred_profs"][match_ind][0][cut_start:cut_end])
with h5py.File(imp_scores_hdf5_path, "r") as f:
match_inds = np.where(
(f["coords_chrom"][:].astype(str) == chrom) &
(f["coords_start"][:] <= score_start) &
(f["coords_end"][:] >= score_end)
)[0]
if not match_inds.size:
print("Warning: did not find sufficiently large importance score track for %s:%d-%s" % (chrom, score_start, score_end))
return
match_ind = match_inds[0]
coord_start, coord_end = f["coords_start"][match_ind], f["coords_end"][match_ind]
hyp_scores = f[hyp_score_key][match_ind]
one_hot_seq = f["input_seqs"][match_ind]
cut_start = score_start - coord_start
cut_end = cut_start + score_center_size
hyp_scores = hyp_scores[cut_start:cut_end]
one_hot_seq = one_hot_seq[cut_start:cut_end]
repeats_matches = repeats_table[
(repeats_table["genoName"] == chrom) &
(repeats_table["genoStart"] < prof_end) &
(repeats_table["genoEnd"] > prof_start)
]
prof_fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(20, 10))
# Draw profiles
ax[0].plot(true_profs[:, 0], color="darkslateblue")
ax[0].plot(-true_profs[:, 1], color="darkorange")
ax[0].set_title("True ChIP-exo profiles")
ax[1].plot(pred_profs[:, 0], color="darkslateblue")
ax[1].plot(-pred_profs[:, 1], color="darkorange")
ax[1].set_title("Predicted ChIP-exo profiles")
# Draw repeat annotation(s)
for _, row in repeats_matches.iterrows():
start_pos = max(row["genoStart"] - prof_start, 0)
end_pos = min(row["genoEnd"] - prof_start, prof_center_size)
ax[2].add_patch(patches.Rectangle(
xy=(start_pos, 0), width=(end_pos - start_pos), height=1, color="red", fill=False
))
ax[2].text(
x=((end_pos + start_pos) // 2), y=0.5, s=row["repName"]
)
# Draw vertical lines that denote the portion with importance scores
for i in range(3):
ax[i].axvline(score_start - prof_start, color="gray")
ax[i].axvline(score_end - prof_start, color="gray")
if save_path:
plt.savefig(
os.path.join(save_path + "_profiles.svg"), format="svg"
)
plt.show()
score_fig = viz_sequence.plot_weights(
hyp_scores * one_hot_seq, figsize=(20, 4), subticks_frequency=score_center_size, return_fig=True
)
score_fig.tight_layout()
if save_path:
plt.savefig(
os.path.join(save_path + "_impscores.svg"), format="svg"
)
plt.show()
print(one_hot_to_seq(one_hot_seq))
