Direct links to results¶
In [1]:
import os
import util
from tomtom import match_motifs_to_database
import viz_sequence
import numpy as np
import pandas as pd
import sklearn.cluster
import scipy.cluster.hierarchy
import matplotlib.pyplot as plt
import matplotlib.font_manager as font_manager
import vdom.helpers as vdomh
from IPython.display import display
import tqdm
tqdm.tqdm_notebook()
Out[1]:
In [2]:
# Plotting defaults
font_manager.fontManager.ttflist.extend(
font_manager.createFontList(
font_manager.findSystemFonts(fontpaths="/users/amtseng/modules/fonts")
)
)
plot_params = {
"figure.titlesize": 22,
"axes.titlesize": 22,
"axes.labelsize": 20,
"legend.fontsize": 18,
"xtick.labelsize": 16,
"ytick.labelsize": 16,
"font.family": "Roboto",
"font.weight": "bold"
}
plt.rcParams.update(plot_params)
Define constants and paths¶
In [3]:
# Define parameters/fetch arguments
preds_path = os.environ["TFM_PRED_PATH"]
shap_scores_path = os.environ["TFM_SHAP_PATH"]
tfm_results_path = os.environ["TFM_TFM_PATH"]
peak_bed_paths = [os.environ["TFM_PEAKS_PATH"]]
print("Predictions path: %s" % preds_path)
print("DeepSHAP scores path: %s" % shap_scores_path)
print("TF-MoDISco results path: %s" % tfm_results_path)
print("Peaks path: %s" % peak_bed_paths[0])
In [4]:
# Define constants
input_length, profile_length = 2114, 1000
shap_score_center_size = 400
profile_display_center_size = 400
hyp_score_key = "hyp_scores"
task_index = None
Helper functions¶
In [5]:
def extract_profiles_and_coords(
seqlets_arr, one_hot_seqs, hyp_scores, true_profs, pred_profs, pred_coords,
input_length, profile_length, input_center_cut_size, profile_center_cut_size,
task_index=None
):
"""
From the seqlets object of a TF-MoDISco pattern's seqlets and alignments,
extracts the predicted and observed profiles of the model, as well as the
set of coordinates for the seqlets.
Arguments:
`seqlets_arr`: a TF-MoDISco pattern's seqlets object array (N-array)
`one_hot_seqs`: an N x R x 4 array of input sequences, where R is
the cut centered size
`hyp_scores`: an N x R x 4 array of hypothetical importance scores
`true_profs`: an N x T x O x 2 array of true profile counts
`pred_profs`: an N x T x O x 2 array of predicted profile probabilities
`pred_coords`: an N x 3 object array of coordinates for the input sequence
underlying the predictions
`input_length`: length of original input sequences, I
`profile_length`: length of profile predictions, O
`input_center_cut_size`: centered cut size of SHAP scores used
`profile_center_cut_size`: size to cut profiles to when returning them, P
`task_index`: index of task to focus on for profiles; if None, returns
profiles for all tasks
Returns an N x (T or 1) x P x 2 array of true profile counts, an
N x (T or 1) x P x 2 array of predicted profile probabilities, an N x Q x 4
array of one-hot seqlet sequences, an N x Q x 4 array of hypothetical seqlet
importance scores, and an N x 3 object array of seqlet coordinates, where P
is the profile cut size and Q is the seqlet length. Returned profiles are
centered at the same center as the seqlets.
Note that it is important that the seqlet indices match exactly with the indices
out of the N. This should be the exact sequences in the original SHAP scores.
"""
true_seqlet_profs, pred_seqlet_profs, seqlet_seqs, seqlet_hyps, seqlet_coords = [], [], [], [], []
def seqlet_coord_to_profile_coord(seqlet_coord):
return seqlet_coord + ((input_length - input_center_cut_size) // 2) - ((input_length - profile_length) // 2)
def seqlet_coord_to_input_coord(seqlet_coord):
return seqlet_coord + ((input_length - input_center_cut_size) // 2)
# For each seqlet, fetch the true/predicted profiles
for seqlet in seqlets_arr:
coord_index = seqlet.coor.example_idx
seqlet_start = seqlet.coor.start
seqlet_end = seqlet.coor.end
seqlet_rc = seqlet.coor.is_revcomp
# Get indices of profile to cut out
seqlet_center = (seqlet_start + seqlet_end) // 2
prof_center = seqlet_coord_to_profile_coord(seqlet_center)
prof_start = prof_center - (profile_center_cut_size // 2)
prof_end = prof_start + profile_center_cut_size
if task_index is None or true_profs.shape[1] == 1:
# Use all tasks if the predictions only have 1 task to begin with
task_start, task_end = None, None
else:
task_start, task_end = task_index, task_index + 1
true_prof = true_profs[coord_index, task_start:task_end, prof_start:prof_end] # (T or 1) x P x 2
pred_prof = pred_profs[coord_index, task_start:task_end, prof_start:prof_end] # (T or 1) x P x 2
true_seqlet_profs.append(true_prof)
pred_seqlet_profs.append(pred_prof)
# The one-hot-sequences and hypothetical scores are assumed to already by cut/centered,
# so the indices match the seqlet indices
if seqlet_rc:
seqlet_seqs.append(np.flip(one_hot_seqs[coord_index, seqlet_start:seqlet_end], axis=(0, 1)))
seqlet_hyps.append(np.flip(hyp_scores[coord_index, seqlet_start:seqlet_end], axis=(0, 1)))
else:
seqlet_seqs.append(one_hot_seqs[coord_index, seqlet_start:seqlet_end])
seqlet_hyps.append(hyp_scores[coord_index, seqlet_start:seqlet_end])
# Get the coordinates of the seqlet based on the input coordinates
inp_start = seqlet_coord_to_input_coord(seqlet_start)
inp_end = seqlet_coord_to_input_coord(seqlet_end)
chrom, start, _ = pred_coords[coord_index]
seqlet_coords.append([chrom, start + inp_start, start + inp_end])
return np.stack(true_seqlet_profs), np.stack(pred_seqlet_profs), np.stack(seqlet_seqs), np.stack(seqlet_hyps), np.array(seqlet_coords, dtype=object)
In [6]:
def plot_profiles(seqlet_true_profs, seqlet_pred_profs, kmeans_clusters=5):
"""
Plots the given profiles with a heatmap.
Arguments:
`seqlet_true_profs`: an N x O x 2 NumPy array of true profiles, either as raw
counts or probabilities (they will be normalized)
`seqlet_pred_profs`: an N x O x 2 NumPy array of predicted profiles, either as
raw counts or probabilities (they will be normalized)
`kmeans_cluster`: when displaying profile heatmaps, there will be this
many clusters
"""
assert len(seqlet_true_profs.shape) == 3
assert seqlet_true_profs.shape == seqlet_pred_profs.shape
num_profs, width, _ = seqlet_true_profs.shape
# First, normalize the profiles along the output profile dimension
def normalize(arr, axis=0):
arr_sum = np.sum(arr, axis=axis, keepdims=True)
arr_sum[arr_sum == 0] = 1 # If 0, keep 0 as the quotient instead of dividing by 0
return arr / arr_sum
true_profs_norm = normalize(seqlet_true_profs, axis=1)
pred_profs_norm = normalize(seqlet_pred_profs, axis=1)
# Compute the mean profiles across all examples
true_profs_mean = np.mean(true_profs_norm, axis=0)
pred_profs_mean = np.mean(pred_profs_norm, axis=0)
# Perform k-means clustering on the predicted profiles, with the strands pooled
kmeans_clusters = max(5, num_profs // 50) # Set number of clusters based on number of profiles, with minimum
kmeans = sklearn.cluster.KMeans(n_clusters=kmeans_clusters)
cluster_assignments = kmeans.fit_predict(
np.reshape(pred_profs_norm, (pred_profs_norm.shape[0], -1))
)
# Perform hierarchical clustering on the cluster centers to determine optimal ordering
kmeans_centers = kmeans.cluster_centers_
cluster_order = scipy.cluster.hierarchy.leaves_list(
scipy.cluster.hierarchy.optimal_leaf_ordering(
scipy.cluster.hierarchy.linkage(kmeans_centers, method="centroid"), kmeans_centers
)
)
# Order the profiles so that the cluster assignments follow the optimal ordering
cluster_inds = []
for cluster_id in cluster_order:
cluster_inds.append(np.where(cluster_assignments == cluster_id)[0])
cluster_inds = np.concatenate(cluster_inds)
# Compute a matrix of profiles, normalized to the maximum height, ordered by clusters
def make_profile_matrix(flat_profs, order_inds):
matrix = flat_profs[order_inds]
maxes = np.max(matrix, axis=1, keepdims=True)
maxes[maxes == 0] = 1 # If 0, keep 0 as the quotient instead of dividing by 0
return matrix / maxes
true_matrix = make_profile_matrix(true_profs_norm, cluster_inds)
pred_matrix = make_profile_matrix(pred_profs_norm, cluster_inds)
# Create a figure with the right dimensions
mean_height = 4
heatmap_height = min(num_profs * 0.004, 8)
fig_height = mean_height + (2 * heatmap_height)
fig, ax = plt.subplots(
3, 2, figsize=(16, fig_height), sharex=True,
gridspec_kw={
"width_ratios": [1, 1],
"height_ratios": [mean_height / fig_height, heatmap_height / fig_height, heatmap_height / fig_height]
}
)
# Plot the average predictions
ax[0, 0].plot(true_profs_mean[:, 0], color="darkslateblue")
ax[0, 0].plot(-true_profs_mean[:, 1], color="darkorange")
ax[0, 1].plot(pred_profs_mean[:, 0], color="darkslateblue")
ax[0, 1].plot(-pred_profs_mean[:, 1], color="darkorange")
# Set axes on average predictions
max_mean_val = max(np.max(true_profs_mean), np.max(pred_profs_mean))
mean_ylim = max_mean_val * 1.05 # Make 5% higher
ax[0, 0].set_title("True profiles")
ax[0, 0].set_ylabel("Average probability")
ax[0, 1].set_title("Predicted profiles")
for j in (0, 1):
ax[0, j].set_ylim(-mean_ylim, mean_ylim)
ax[0, j].label_outer()
# Plot the heatmaps
ax[1, 0].imshow(true_matrix[:, :, 0], interpolation="nearest", aspect="auto", cmap="Blues")
ax[1, 1].imshow(pred_matrix[:, :, 0], interpolation="nearest", aspect="auto", cmap="Blues")
ax[2, 0].imshow(true_matrix[:, :, 1], interpolation="nearest", aspect="auto", cmap="Oranges")
ax[2, 1].imshow(pred_matrix[:, :, 1], interpolation="nearest", aspect="auto", cmap="Oranges")
# Set axes on heatmaps
for i in (1, 2):
for j in (0, 1):
ax[i, j].set_yticks([])
ax[i, j].set_yticklabels([])
ax[i, j].label_outer()
width = true_matrix.shape[1]
delta = 100
num_deltas = (width // 2) // delta
labels = list(range(max(-width // 2, -num_deltas * delta), min(width // 2, num_deltas * delta) + 1, delta))
tick_locs = [label + max(width // 2, num_deltas * delta) for label in labels]
for j in (0, 1):
ax[2, j].set_xticks(tick_locs)
ax[2, j].set_xticklabels(labels)
ax[2, j].set_xlabel("Distance from peak summit (bp)")
fig.tight_layout()
plt.show()
In [7]:
def get_summit_distances(coords, peak_table):
"""
Given a set of coordinates, computes the distance of the center of each
coordinate to the nearest summit.
Arguments:
`coords`: an N x 3 object array of coordinates
`peak_table`: a 6-column table of peak data, as imported by
`import_peak_table`
Returns and N-array of integers, which is the distance of each coordinate
midpoint to the nearest coordinate.
"""
chroms = coords[:, 0]
midpoints = (coords[:, 1] + coords[:, 2]) // 2
dists = []
for i in range(len(coords)):
chrom = chroms[i]
midpoint = midpoints[i]
rows = peak_table[peak_table["chrom"] == chrom]
dist_arr = (midpoint - rows["summit"]).values
min_dist = dist_arr[np.argmin(np.abs(dist_arr))]
dists.append(min_dist)
return np.array(dists)
In [8]:
def plot_summit_dists(summit_dists):
"""
Plots the distribution of seqlet distances to summits.
Arguments:
`summit_dists`: the array of distances as returned by
`get_summit_distances`
"""
plt.figure(figsize=(8, 6))
num_bins = max(len(summit_dists) // 30, 20)
plt.hist(summit_dists, bins=num_bins, color="purple")
plt.title("Histogram of distance of seqlets to peak summits")
plt.xlabel("Signed distance from seqlet center to nearest peak summit (bp)")
plt.show()
Import SHAP scores, profile predictions, and TF-MoDISco results¶
In [9]:
# Import SHAP coordinates and one-hot sequences
hyp_scores, _, one_hot_seqs, shap_coords = util.import_shap_scores(shap_scores_path, hyp_score_key, center_cut_size=shap_score_center_size, remove_non_acgt=False)
# This cuts the sequences/scores off just as how TF-MoDISco saw them, but the coordinates are uncut
In [10]:
# Import the set of all profiles and their coordinates
true_profs, pred_profs, all_pred_coords = util.import_profiles(preds_path)
In [11]:
# Import the set of peaks
peak_table = util.import_peak_table(peak_bed_paths)
In [12]:
# Subset the predicted profiles/coordinates to the task-specific SHAP coordinates/scores
shap_coords_table = pd.DataFrame(shap_coords, columns=["chrom", "start", "end"])
pred_coords_table = pd.DataFrame(all_pred_coords, columns=["chrom", "start", "end"])
subset_inds = pred_coords_table.reset_index().drop_duplicates(["chrom", "start", "end"]).merge(
shap_coords_table.reset_index(), on=["chrom", "start", "end"]
).sort_values("index_y")["index_x"].values
true_profs = true_profs[subset_inds]
pred_profs = pred_profs[subset_inds]
pred_coords = all_pred_coords[subset_inds]
# Make sure the coordinates all match
assert np.all(pred_coords == shap_coords)
In [13]:
# Import the TF-MoDISco results object
tfm_obj = util.import_tfmodisco_results(tfm_results_path, hyp_scores, one_hot_seqs, shap_score_center_size)
Plot some SHAP score tracks¶
Plot the central region of some randomly selected actual importance scores
In [14]:
for index in np.random.choice(hyp_scores.shape[0], size=5, replace=False):
viz_sequence.plot_weights((hyp_scores[index] * one_hot_seqs[index])[100:300], subticks_frequency=100)
In [15]:
motif_pfms, motif_hcwms, motif_cwms = [], [], [] # Save the trimmed PFMs, hCWMs, and CWMs
motif_pfms_short = [] # PFMs that are even more trimmed (for TOMTOM)
num_seqlets = [] # Number of seqlets for each motif
motif_seqlets = [] # Save seqlets of each motif
metaclusters = tfm_obj.metacluster_idx_to_submetacluster_results
num_metaclusters = len(metaclusters.keys())
for metacluster_i, metacluster_key in enumerate(metaclusters.keys()):
metacluster = metaclusters[metacluster_key]
display(vdomh.h3("Metacluster %d/%d" % (metacluster_i + 1, num_metaclusters)))
patterns = metacluster.seqlets_to_patterns_result.patterns
if not patterns:
break
motif_pfms.append([])
motif_hcwms.append([])
motif_cwms.append([])
motif_pfms_short.append([])
num_seqlets.append([])
motif_seqlets.append([])
num_patterns = len(patterns)
for pattern_i, pattern in enumerate(patterns):
seqlets = pattern.seqlets
display(vdomh.h4("Pattern %d/%d" % (pattern_i + 1, num_patterns)))
display(vdomh.p("%d seqlets" % len(seqlets)))
pfm = pattern["sequence"].fwd
hcwm = pattern["task0_hypothetical_contribs"].fwd
cwm = pattern["task0_contrib_scores"].fwd
pfm_fig = viz_sequence.plot_weights(pfm, subticks_frequency=10, return_fig=True)
hcwm_fig = viz_sequence.plot_weights(hcwm, subticks_frequency=10, return_fig=True)
cwm_fig = viz_sequence.plot_weights(cwm, subticks_frequency=10, return_fig=True)
pfm_fig.tight_layout()
hcwm_fig.tight_layout()
cwm_fig.tight_layout()
motif_table = vdomh.table(
vdomh.tr(
vdomh.td("Sequence (PFM)"),
vdomh.td(util.figure_to_vdom_image(pfm_fig))
),
vdomh.tr(
vdomh.td("Hypothetical contributions (hCWM)"),
vdomh.td(util.figure_to_vdom_image(hcwm_fig))
),
vdomh.tr(
vdomh.td("Actual contributions (CWM)"),
vdomh.td(util.figure_to_vdom_image(cwm_fig))
)
)
display(motif_table)
plt.close("all") # Remove all standing figures
# Trim motif based on total contribution
score = np.sum(np.abs(cwm), axis=1)
trim_thresh = np.max(score) * 0.2 # Cut off anything less than 20% of max score
pass_inds = np.where(score >= trim_thresh)[0]
short_trimmed_pfm = pfm[np.min(pass_inds): np.max(pass_inds) + 1]
motif_pfms_short[-1].append(short_trimmed_pfm)
# Expand trimming to +/- 4bp on either side
start, end = max(0, np.min(pass_inds) - 4), min(len(cwm), np.max(pass_inds) + 4 + 1)
trimmed_pfm = pfm[start:end]
trimmed_hcwm = hcwm[start:end]
trimmed_cwm = cwm[start:end]
motif_pfms[-1].append(trimmed_pfm)
motif_hcwms[-1].append(trimmed_hcwm)
motif_cwms[-1].append(trimmed_cwm)
num_seqlets[-1].append(len(seqlets))
seqlet_true_profs, seqlet_pred_profs, seqlet_seqs, seqlet_hyps, seqlet_coords = extract_profiles_and_coords(
seqlets, one_hot_seqs, hyp_scores, true_profs, pred_profs, pred_coords,
input_length, profile_length, shap_score_center_size,
profile_display_center_size, task_index=task_index
)
motif_seqlets[-1].append((seqlet_seqs, seqlet_hyps))
assert np.allclose(np.sum(seqlet_seqs, axis=0) / len(seqlet_seqs), pattern["sequence"].fwd)
# ^Sanity check: PFM derived from seqlets match the PFM stored in the pattern
plot_profiles(
# Flatten to NT x O x 2
np.reshape(seqlet_true_profs, (-1, seqlet_true_profs.shape[2], seqlet_true_profs.shape[3])),
np.reshape(seqlet_pred_profs, (-1, seqlet_pred_profs.shape[2], seqlet_pred_profs.shape[3]))
)
summit_dists = get_summit_distances(seqlet_coords, peak_table)
plot_summit_dists(summit_dists)
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
| Sequence (PFM) |  |
| Hypothetical contributions (hCWM) |  |
| Actual contributions (CWM) |  |
Summary of motifs¶
Motifs are trimmed based on information content, and presented in descending order by number of supporting seqlets. The motifs are separated by metacluster. The motifs are presented as hCWMs. The forward orientation is defined as the orientation that is richer in purines.
In [16]:
colgroup = vdomh.colgroup(
vdomh.col(style={"width": "5%"}),
vdomh.col(style={"width": "5%"}),
vdomh.col(style={"width": "45%"}),
vdomh.col(style={"width": "45%"})
)
header = vdomh.thead(
vdomh.tr(
vdomh.th("#", style={"text-align": "center"}),
vdomh.th("Seqlets", style={"text-align": "center"}),
vdomh.th("Forward", style={"text-align": "center"}),
vdomh.th("Reverse", style={"text-align": "center"})
)
)
for i in range(len(motif_hcwms)):
display(vdomh.h3("Metacluster %d/%d" % (i + 1, num_metaclusters)))
body = []
for j in range(len(motif_hcwms[i])):
motif = motif_hcwms[i][j]
if np.sum(motif[:, [0, 2]]) > 0.5 * np.sum(motif):
# Forward is purine-rich, reverse-complement is pyrimidine-rich
f, rc = motif, np.flip(motif, axis=(0, 1))
else:
f, rc = np.flip(motif, axis=(0, 1)), motif
f_fig = viz_sequence.plot_weights(f, figsize=(20, 4), return_fig=True)
f_fig.tight_layout()
rc_fig = viz_sequence.plot_weights(rc, figsize=(20, 4), return_fig=True)
rc_fig.tight_layout()
body.append(
vdomh.tr(
vdomh.td(str(j + 1)),
vdomh.td(str(num_seqlets[i][j])),
vdomh.td(util.figure_to_vdom_image(f_fig)),
vdomh.td(util.figure_to_vdom_image(rc_fig))
)
)
display(vdomh.table(colgroup, header, vdomh.tbody(*body)))
plt.close("all")
| # | Seqlets | Forward | Reverse |
|---|---|---|---|
| 1 | 9059 |  |  |
| 2 | 695 |  |  |
| 3 | 71 |  |  |
| 4 | 67 |  |  |
| 5 | 65 |  |  |
| 6 | 81 |  |  |
| 7 | 55 |  |  |
| 8 | 57 |  |  |
| 9 | 38 |  |  |
Top TOMTOM matches for each motif¶
Here, the TF-MoDISco motifs are plotted as hCWMs, but the TOMTOM matches are shown as PWMs.
In [17]:
num_matches_to_keep = 10
num_matches_to_show = 5
header = vdomh.thead(
vdomh.tr(
vdomh.th("Motif ID", style={"text-align": "center"}),
vdomh.th("q-val", style={"text-align": "center"}),
vdomh.th("PWM", style={"text-align": "center"})
)
)
for i in range(len(motif_pfms)):
display(vdomh.h3("Metacluster %d/%d" % (i + 1, num_metaclusters)))
# Compute TOMTOM matches for all motifs in the metacluster at once
tomtom_matches = match_motifs_to_database(motif_pfms_short[i], top_k=num_matches_to_keep)
for j in range(len(motif_pfms[i])):
display(vdomh.h4("Motif %d/%d" % (j + 1, len(motif_pfms[i]))))
viz_sequence.plot_weights(motif_hcwms[i][j])
body = []
for k, (match_name, match_pfm, match_qval) in enumerate(tomtom_matches[j]):
fig = viz_sequence.plot_weights(util.pfm_to_pwm(match_pfm), return_fig=True)
fig.tight_layout()
if k < num_matches_to_show:
body.append(
vdomh.tr(
vdomh.td(match_name),
vdomh.td(str(match_qval)),
vdomh.td(util.figure_to_vdom_image(fig))
)
)
else:
body.append(
vdomh.tr(
vdomh.td(match_name),
vdomh.td(str(match_qval)),
vdomh.td("Not shown")
)
)
if not body:
display(vdomh.p("No TOMTOM matches passing threshold"))
else:
display(vdomh.table(header, vdomh.tbody(*body)))
plt.close("all")
| Motif ID | q-val | PWM |
|---|---|---|
| GATA2_HUMAN.H11MO.0.A | 1.10049e-09 |  |
| GATA1_HUMAN.H11MO.0.A | 0.0009263489999999999 |  |
| TAL1_HUMAN.H11MO.0.A | 0.0009263489999999999 |  |
| GATA4_HUMAN.H11MO.0.A | 0.012155899999999999 |  |
| GATA1_HUMAN.H11MO.1.A | 0.0201947 |  |
| GATA2_HUMAN.H11MO.1.A | 0.0201947 | Not shown |
| GATA3_HUMAN.H11MO.0.A | 0.0475029 | Not shown |
| MA0036.3_GATA2 | 0.0505214 | Not shown |
| GATA6_HUMAN.H11MO.0.A | 0.0658582 | Not shown |
| MA0482.2_GATA4 | 0.0732721 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| FOSL1_HUMAN.H11MO.0.A | 6.451390000000001e-07 |  |
| FOSL2_HUMAN.H11MO.0.A | 6.680100000000001e-07 |  |
| JUN_HUMAN.H11MO.0.A | 7.34134e-07 |  |
| MA0489.1_JUN(var.2) | 9.2411e-07 |  |
| JUND_HUMAN.H11MO.0.A | 1.08289e-06 |  |
| MA1144.1_FOSL2::JUND | 6.40837e-06 | Not shown |
| MA1135.1_FOSB::JUNB | 7.40885e-06 | Not shown |
| MA1138.1_FOSL2::JUNB | 7.40885e-06 | Not shown |
| MA0478.1_FOSL2 | 8.70471e-06 | Not shown |
| MA0477.2_FOSL1 | 1.3088599999999999e-05 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| MA0037.3_GATA3 | 0.000177706 |  |
| GATA2_HUMAN.H11MO.1.A | 0.000890688 |  |
| MA0766.2_GATA5 | 0.000890688 |  |
| MA0036.3_GATA2 | 0.00555745 |  |
| GATA1_HUMAN.H11MO.1.A | 0.00578997 | |
| GATA4_HUMAN.H11MO.0.A | 0.00830561 | Not shown |
| GATA3_HUMAN.H11MO.0.A | 0.00871587 | Not shown |
| TAL1_HUMAN.H11MO.0.A | 0.00958375 | Not shown |
| GATA2_HUMAN.H11MO.0.A | 0.0114812 | Not shown |
| GATA1_HUMAN.H11MO.0.A | 0.014897299999999999 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| OLIG2_HUMAN.H11MO.1.B | 0.0109375 |  |
| MA1472.1_BHLHA15(var.2) | 0.09438730000000001 |  |
| ATOH1_HUMAN.H11MO.0.B | 0.09438730000000001 |  |
| PTF1A_HUMAN.H11MO.1.B | 0.09438730000000001 |  |
| LYL1_HUMAN.H11MO.0.A | 0.09438730000000001 |  |
| NDF2_HUMAN.H11MO.0.B | 0.0945005 | Not shown |
| MA1467.1_ATOH1(var.2) | 0.0945005 | Not shown |
| MA0665.1_MSC | 0.125325 | Not shown |
| MA0816.1_Ascl2 | 0.125325 | Not shown |
| FOXC2_HUMAN.H11MO.0.D | 0.131924 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| OLIG2_HUMAN.H11MO.1.B | 0.0035861999999999995 | |
| LYL1_HUMAN.H11MO.0.A | 0.019397400000000002 | |
| NDF2_HUMAN.H11MO.0.B | 0.019397400000000002 |  |
| NDF1_HUMAN.H11MO.0.A | 0.019397400000000002 |  |
| ATOH1_HUMAN.H11MO.0.B | 0.019397400000000002 | |
| MA0665.1_MSC | 0.0330346 | Not shown |
| MA1467.1_ATOH1(var.2) | 0.0330346 | Not shown |
| MA1472.1_BHLHA15(var.2) | 0.0330834 | Not shown |
| MA0091.1_TAL1::TCF3 | 0.0347904 | Not shown |
| MA0667.1_MYF6 | 0.043989 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| MA1472.1_BHLHA15(var.2) | 0.00601984 | |
| MA0816.1_Ascl2 | 0.00601984 |  |
| LYL1_HUMAN.H11MO.0.A | 0.00601984 | |
| OLIG1_HUMAN.H11MO.0.D | 0.00601984 |  |
| OLIG2_HUMAN.H11MO.1.B | 0.00836597 | |
| MA1467.1_ATOH1(var.2) | 0.0114053 | Not shown |
| BHE22_HUMAN.H11MO.0.D | 0.0114053 | Not shown |
| MA0665.1_MSC | 0.0114053 | Not shown |
| ZN317_HUMAN.H11MO.0.C | 0.0119941 | Not shown |
| ASCL2_HUMAN.H11MO.0.D | 0.0121946 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| GATA6_HUMAN.H11MO.0.A | 0.0372414 |  |
| TAL1_HUMAN.H11MO.0.A | 0.054278999999999994 | |
| MA1104.2_GATA6 | 0.0547539 |  |
| GATA2_HUMAN.H11MO.0.A | 0.0547539 | |
| GATA4_HUMAN.H11MO.0.A | 0.130679 | |
| MA0037.3_GATA3 | 0.16400399999999998 | Not shown |
| GATA1_HUMAN.H11MO.0.A | 0.16400399999999998 | Not shown |
| MA0029.1_Mecom | 0.16400399999999998 | Not shown |
| GATA5_HUMAN.H11MO.0.D | 0.200975 | Not shown |
| GATA3_HUMAN.H11MO.0.A | 0.367261 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| GATA4_HUMAN.H11MO.0.A | 0.20747800000000002 | |
| GATA2_HUMAN.H11MO.1.A | 0.346831 | |
| GATA1_HUMAN.H11MO.1.A | 0.47356800000000004 | |
| MA0766.2_GATA5 | 0.47356800000000004 | |
| MA0036.3_GATA2 | 0.47356800000000004 | |
| GATA3_HUMAN.H11MO.0.A | 0.47356800000000004 | Not shown |
| MA0037.3_GATA3 | 0.47356800000000004 | Not shown |
| Motif ID | q-val | PWM |
|---|---|---|
| MA0036.3_GATA2 | 0.027977600000000002 | |
| MA0037.3_GATA3 | 0.027977600000000002 | |
| GATA6_HUMAN.H11MO.0.A | 0.027977600000000002 | |
| GATA2_HUMAN.H11MO.1.A | 0.0361672 | |
| Z324A_HUMAN.H11MO.0.C | 0.037268199999999994 |  |
| MA0482.2_GATA4 | 0.037268199999999994 | Not shown |
| GATA3_HUMAN.H11MO.0.A | 0.037268199999999994 | Not shown |
| MA0766.2_GATA5 | 0.037268199999999994 | Not shown |
| GATA4_HUMAN.H11MO.0.A | 0.06606210000000001 | Not shown |
| MA1104.2_GATA6 | 0.06702430000000001 | Not shown |
Sample of seqlets supporting each motif¶
Here, the motifs are presented as hCWMs, along with the actual importance scores of a random sample of seqlets that support the motif.
In [18]:
num_seqlets_to_show = 10
colgroup = vdomh.colgroup(
vdomh.col(style={"width": "50%"}),
vdomh.col(style={"width": "50%"})
)
header = vdomh.thead(
vdomh.tr(
vdomh.th("Motif hCWM", style={"text-align": "center"}),
vdomh.th("Seqlets", style={"text-align": "center"})
)
)
for i in range(len(motif_hcwms)):
display(vdomh.h3("Metacluster %d/%d" % (i + 1, num_metaclusters)))
for j in range(len(motif_hcwms[i])):
display(vdomh.h4("Motif %d/%d" % (j + 1, len(motif_hcwms[i]))))
motif_fig = viz_sequence.plot_weights(motif_hcwms[i][j], figsize=(20, 4), return_fig=True)
motif_fig.tight_layout()
seqlet_seqs, seqlet_hyps = motif_seqlets[i][j]
sample_size = min(num_seqlets_to_show, len(seqlet_seqs))
sample_inds = np.random.choice(len(seqlet_seqs), size=sample_size, replace=False)
sample = []
for k in sample_inds:
fig = viz_sequence.plot_weights(seqlet_hyps[k] * seqlet_seqs[k], subticks_frequency=10, return_fig=True)
fig.tight_layout()
sample.append(util.figure_to_vdom_image(fig))
body = vdomh.tbody(vdomh.tr(vdomh.td(util.figure_to_vdom_image(motif_fig)), vdomh.td(*sample)))
display(vdomh.table(colgroup, header, body))
plt.close("all")
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |
| Motif hCWM | Seqlets |
|---|---|
|           |