def plot_example_hits(
chrom, start, end, tfm_hits_table, moods_hits_table, profiles_hdf5_path,
imp_scores_hdf5_path, task_index, pfm_dict, prof_center_size=700, score_center_size=100,
hyp_score_key="profile_hyp_scores", save_path=None
):
"""
For a given region, plots the true/predicted profiles, importance scores, and the
TFM and MOODS hits at that region on each.
"""
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(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
if f["predictions"]["log_pred_profs"].shape[1] == 1:
task_index = 0
pred_profs = np.exp(f["predictions"]["log_pred_profs"][match_ind][task_index][cut_start:cut_end])
true_profs = f["predictions"]["true_profs"][match_ind][task_index][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]
tfm_hits = tfm_hits_table[
(tfm_hits_table["chrom"] == chrom) &
(tfm_hits_table["start"] < end) &
(tfm_hits_table["end"] > start)
]
moods_hits = moods_hits_table[
(moods_hits_table["chrom"] == chrom) &
(moods_hits_table["start"] < end) &
(moods_hits_table["end"] > start)
]
prof_fig, ax = plt.subplots(nrows=3, sharex=True, figsize=(20, 8))
# Draw profiles
ax[0].plot(true_profs[:, 0], color="darkslateblue")
ax[0].plot(-true_profs[:, 1], color="darkorange")
ax[0].set_title("True ChIP-seq profiles")
ax[1].plot(pred_profs[:, 0], color="darkslateblue")
ax[1].plot(-pred_profs[:, 1], color="darkorange")
ax[1].set_title("Predicted ChIP-seq profiles")
# Draw motif hits
for i, (hit_table, color) in enumerate([(tfm_hits, "blue"), (moods_hits, "red")]):
for _, row in hit_table.iterrows():
start_pos = max(row["start"] - prof_start, 0)
end_pos = min(row["end"] - prof_start, prof_center_size)
ax[2].add_patch(patches.Rectangle(
xy=(start_pos, i * 0.5), width=(end_pos - start_pos), height=0.5, color=color, fill=False
))
ax[2].text(
x=((end_pos + start_pos) // 2), y=((i * 0.5) + 0.25), s=row["key"]
)
# Draw vertical lines that denote the portion with importance scores
for i in range(2):
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()
motif_keys = np.unique(np.concatenate([tfm_hits["key"], moods_hits["key"]]))
fig, ax = plt.subplots(nrows=len(motif_keys), figsize=(20, 2 * len(motif_keys)))
if len(motif_keys) == 1:
ax = [ax]
for i, motif_key in enumerate(motif_keys):
pfm = pfm_dict[motif_key]
pwm = np.log2((pfm + 0.0001) / read_motifs.BACKGROUND_FREQS[None])
match_track = pwm_score_track(pwm, one_hot_seq)
match_track[match_track < 0] = 0
ax[i].plot(match_track)
ax[i].set_title(motif_key)
fig.suptitle("PWM log-odds scores across the sequence")
fig.tight_layout()
if save_path:
plt.savefig(
os.path.join(save_path + "_pwm_logodds.svg"), format="svg"
)
plt.show()
highlights = {}
for i, (hit_table, color) in enumerate([(tfm_hits, "blue"), (moods_hits, "red")]):
highlights[color] = []
for _, row in hit_table.iterrows():
if row["start"] >= score_end or row["end"] <= score_start:
continue
start_pos = max(row["start"] - score_start, 0)
end_pos = min(row["end"] - score_start, score_center_size)
highlights[color].append((start_pos, end_pos))
print(row["key"], start_pos, end_pos)
score_fig = viz_sequence.plot_weights(
hyp_scores * one_hot_seq, figsize=(20, 4), subticks_frequency=score_center_size, highlight=highlights, 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))
