TF-MoDISco results

Summary of motifs

TOMTOM matches to motifs

Sample of seqlets for each motif

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()

/users/amtseng/miniconda3/envs/tfmodisco-mini/lib/python3.7/site-packages/ipykernel_launcher.py:14: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0
Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`
Out[1]:
0it [00:00, ?it/s]
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)

/users/amtseng/miniconda3/envs/tfmodisco-mini/lib/python3.7/site-packages/ipykernel_launcher.py:4: MatplotlibDeprecationWarning: 
The createFontList function was deprecated in Matplotlib 3.2 and will be removed two minor releases later. Use FontManager.addfont instead.
  after removing the cwd from sys.path.

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])

Predictions path: /users/zahoor/TF-Atlas/02-16-2021/predictions/ENCSR000BGZ/profile_predictions.h5
DeepSHAP scores path: /users/zahoor/TF-Atlas/02-16-2021/shap/ENCSR000BGZ/counts_scores_alex_format.h5
TF-MoDISco results path: /users/zahoor/TF-Atlas/02-16-2021/modisco/ENCSR000BGZ/counts/modisco_results.hd5
Peaks path: /users/zahoor/TF-Atlas/data/idr_peaks//ENCFF068YYR.bed.gz
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

Importing SHAP scores: 100%|██████████| 26/26 [00:05<00:00,  4.73it/s]
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)



    















    

















    

















    

















    

















    






























Plot TF-MoDISco results
Plot all motifs by metacluster













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)



    















    







Metacluster 1/1










    







Pattern 1/19










    







6514 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 2/19










    







1517 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 3/19










    







1427 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 4/19










    







1328 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 5/19










    







927 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 6/19










    







779 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 7/19










    







445 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 8/19










    







236 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 9/19










    







217 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 10/19










    







182 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 11/19










    







165 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 12/19










    







171 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 13/19










    







68 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 14/19










    







72 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 15/19










    







63 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 16/19










    







44 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 17/19










    







42 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 18/19










    







40 seqlets










    







Sequence (PFM)
Hypothetical contributions (hCWM)
Actual contributions (CWM)










    

















    

















    







Pattern 19/19










    







35 seqlets










    







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")



    















    







Metacluster 1/1










    






/users/amtseng/TF-Atlas/3M/reports/viz_sequence.py:152: RuntimeWarning: More than 20 figures have been opened. Figures created through the pyplot interface (`matplotlib.pyplot.figure`) are retained until explicitly closed and may consume too much memory. (To control this warning, see the rcParam `figure.max_open_warning`).
  fig = plt.figure(figsize=figsize)










    







#SeqletsForwardReverse
16514
21517
31427
41328
5927
6779
7445
8236
9217
10182
11165
12171
1368
1472
1563
1644
1742
1840
1935























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")



    















    







Metacluster 1/1










    







Motif 1/19










    

















    







Motif IDq-valPWM
ASCL2_HUMAN.H11MO.0.D1.28145e-05
MA1635.1_BHLHE22(var.2)0.000106217
MYF6_HUMAN.H11MO.0.C0.000106217
TFE2_HUMAN.H11MO.0.A0.000106217
ITF2_HUMAN.H11MO.0.C0.000176032
ASCL1_HUMAN.H11MO.0.A0.000185829Not shown
BHA15_HUMAN.H11MO.0.B0.00042569699999999995Not shown
MA1648.1_TCF12(var.2)0.00042569699999999995Not shown
HTF4_HUMAN.H11MO.0.A0.00042569699999999995Not shown
MA0522.3_TCF30.00042569699999999995Not shown










    







Motif 2/19










    

















    







Motif IDq-valPWM
IRF4_HUMAN.H11MO.0.A7.59907e-11
IRF8_HUMAN.H11MO.0.B1.3847000000000002e-07
BC11A_HUMAN.H11MO.0.A0.000682792
SPIB_HUMAN.H11MO.0.A0.000973653
MA0081.2_SPIB0.00147881
SPI1_HUMAN.H11MO.0.A0.00227877Not shown
MA0080.5_SPI10.00741069Not shown
MA0517.1_STAT1::STAT20.011507700000000001Not shown
MA1623.1_Stat20.0117134Not shown
IRF7_HUMAN.H11MO.0.C0.012849600000000001Not shown










    







Motif 3/19










    

















    







Motif IDq-valPWM
BATF_HUMAN.H11MO.0.A5.78604e-05
MA0489.1_JUN(var.2)0.00278567
MA0655.1_JDP20.00594256
MA0490.2_JUNB0.00839013
MA0462.2_BATF::JUN0.010504600000000001
MA0835.2_BATF30.010504600000000001Not shown
MA1634.1_BATF0.010504600000000001Not shown
MA0491.2_JUND0.0114833Not shown
MA0477.2_FOSL10.0139609Not shown
FOSB_HUMAN.H11MO.0.A0.017259200000000002Not shown










    







Motif 4/19










    

















    







Motif IDq-valPWM
RUNX3_HUMAN.H11MO.0.A6.43059e-05
RUNX1_HUMAN.H11MO.0.A6.43059e-05
MA0002.2_RUNX10.000466383
RUNX2_HUMAN.H11MO.0.A0.00192346
PEBB_HUMAN.H11MO.0.C0.00713931
MA0511.2_RUNX20.027234500000000002Not shown
MA0684.2_RUNX30.12984600000000002Not shown
FOXH1_HUMAN.H11MO.0.A0.21936599999999998Not shown
TCF7_HUMAN.H11MO.0.A0.337663Not shown
MA0122.3_Nkx3-20.39828Not shown










    







Motif 5/19










    

















    







Motif IDq-valPWM
NFKB1_HUMAN.H11MO.1.B1.6530899999999998e-05
NFKB2_HUMAN.H11MO.0.B0.000327737
MA0107.1_RELA0.000327737
TF65_HUMAN.H11MO.0.A0.000360759
RELB_HUMAN.H11MO.0.C0.00143517
ZEP1_HUMAN.H11MO.0.D0.00143517Not shown
MA0101.1_REL0.00191432Not shown
MA0778.1_NFKB20.00224795Not shown
MA1117.1_RELB0.00474916Not shown
MA0105.4_NFKB10.00475783Not shown










    







Motif 6/19










    

















    







Motif IDq-valPWM
MA0598.3_EHF0.00025341299999999997
MA0473.3_ELF10.00031239
MA0062.3_GABPA0.00031239
ERG_HUMAN.H11MO.0.A0.00033702300000000004
ETS1_HUMAN.H11MO.0.A0.0004113
ETV5_HUMAN.H11MO.0.C0.00041852199999999995Not shown
ETV2_HUMAN.H11MO.0.B0.00046749599999999997Not shown
SPI1_HUMAN.H11MO.0.A0.00046749599999999997Not shown
FLI1_HUMAN.H11MO.1.A0.00046749599999999997Not shown
ELF2_HUMAN.H11MO.0.C0.0005229880000000001Not shown










    







Motif 7/19










    

















    







Motif IDq-valPWM
BATF_HUMAN.H11MO.0.A2.8555400000000003e-06
BATF3_HUMAN.H11MO.0.B2.8555400000000003e-06
IRF2_HUMAN.H11MO.0.A0.0520754
MA0050.2_IRF10.0520754
IRF1_HUMAN.H11MO.0.A0.059466099999999994
MA1420.1_IRF50.11798299999999999Not shown
MA1509.1_IRF60.120402Not shown
MA0653.1_IRF90.22984200000000002Not shown
IRF7_HUMAN.H11MO.0.C0.25727300000000003Not shown
STAT2_HUMAN.H11MO.0.A0.25727300000000003Not shown










    







Motif 8/19










    

















    







Motif IDq-valPWM
ITF2_HUMAN.H11MO.0.C0.00442095
MA1635.1_BHLHE22(var.2)0.00442095
HTF4_HUMAN.H11MO.0.A0.00849501
MYF6_HUMAN.H11MO.0.C0.010810799999999999
TFE2_HUMAN.H11MO.0.A0.021396099999999998
PTF1A_HUMAN.H11MO.1.B0.021396099999999998Not shown
ASCL2_HUMAN.H11MO.0.D0.021396099999999998Not shown
MA0048.2_NHLH10.0456454Not shown
ZIC3_HUMAN.H11MO.0.B0.10461500000000001Not shown
MA1648.1_TCF12(var.2)0.10461500000000001Not shown










    







Motif 9/19










    

















    







Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A1.20437e-06
MA0139.1_CTCF4.38988e-06
CTCFL_HUMAN.H11MO.0.A7.84471e-06
MA1102.2_CTCFL3.5781399999999996e-05
SP1_HUMAN.H11MO.0.A0.06030219999999999
RFX1_HUMAN.H11MO.0.B0.0701706Not shown
SP2_HUMAN.H11MO.0.A0.0982163Not shown
SP3_HUMAN.H11MO.0.B0.11988299999999999Not shown
THAP1_HUMAN.H11MO.0.C0.119956Not shown
MA0155.1_INSM10.119956Not shown










    







Motif 10/19










    

















    







Motif IDq-valPWM
OLIG2_HUMAN.H11MO.0.B0.00221995
BHA15_HUMAN.H11MO.0.B0.014627600000000001
MYF6_HUMAN.H11MO.0.C0.0208393
ASCL2_HUMAN.H11MO.0.D0.0247384
MYOG_HUMAN.H11MO.0.B0.0247384
MA1472.1_BHLHA15(var.2)0.0247384Not shown
MA1635.1_BHLHE22(var.2)0.0264278Not shown
MA0816.1_Ascl20.0264278Not shown
PTF1A_HUMAN.H11MO.1.B0.0264278Not shown
MA0521.1_Tcf120.0264278Not shown










    







Motif 11/19










    

















    







Motif IDq-valPWM
COE1_HUMAN.H11MO.0.A1.89971e-07
MA1604.1_Ebf20.00028571900000000003
MA1637.1_EBF30.000895921
MA0154.4_EBF10.00350973
ZBT48_HUMAN.H11MO.0.C0.138535
SP2_HUMAN.H11MO.0.A0.379284Not shown










    







Motif 12/19










    

















    







Motif IDq-valPWM
MYOD1_HUMAN.H11MO.0.A0.00174359
MA1472.1_BHLHA15(var.2)0.0217317
MA1100.2_ASCL10.0217317
BHA15_HUMAN.H11MO.0.B0.0217317
ASCL2_HUMAN.H11MO.0.D0.0217317
MA1635.1_BHLHE22(var.2)0.022677799999999998Not shown
MA0816.1_Ascl20.022677799999999998Not shown
MA1467.1_ATOH1(var.2)0.022677799999999998Not shown
PTF1A_HUMAN.H11MO.1.B0.022677799999999998Not shown
MYF6_HUMAN.H11MO.0.C0.022677799999999998Not shown










    







Motif 13/19










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A9.845780000000001e-05
ATF3_HUMAN.H11MO.0.A9.845780000000001e-05
USF2_HUMAN.H11MO.0.A0.000134626
SP3_HUMAN.H11MO.0.B0.000134626
SP1_HUMAN.H11MO.0.A0.00027638799999999997
KLF16_HUMAN.H11MO.0.D0.00118939Not shown
MITF_HUMAN.H11MO.0.A0.0021250999999999996Not shown
MA0609.2_CREM0.00278122Not shown
TFE3_HUMAN.H11MO.0.B0.00320723Not shown
CREM_HUMAN.H11MO.0.C0.005821100000000001Not shown










    







Motif 14/19










    

















    







Motif IDq-valPWM
LYL1_HUMAN.H11MO.0.A0.100458
MYF6_HUMAN.H11MO.0.C0.100458
TAL1_HUMAN.H11MO.1.A0.100458
MEIS1_HUMAN.H11MO.1.B0.100458
MA1635.1_BHLHE22(var.2)0.13356300000000002
MA1109.1_NEUROD10.13356300000000002Not shown
MA0797.1_TGIF20.13356300000000002Not shown
ITF2_HUMAN.H11MO.0.C0.13356300000000002Not shown
MA0796.1_TGIF10.13356300000000002Not shown
ASCL2_HUMAN.H11MO.0.D0.13356300000000002Not shown










    







Motif 15/19










    

















    







Motif IDq-valPWM
BHA15_HUMAN.H11MO.0.B0.0114651
MYF6_HUMAN.H11MO.0.C0.0190421
MYOD1_HUMAN.H11MO.0.A0.0511216
ASCL2_HUMAN.H11MO.0.D0.0511216
MA0521.1_Tcf120.0511216
TFE2_HUMAN.H11MO.0.A0.0511216Not shown
MYOD1_HUMAN.H11MO.1.A0.0511216Not shown
PTF1A_HUMAN.H11MO.1.B0.0511216Not shown
MA0500.2_MYOG0.0511216Not shown
MA0816.1_Ascl20.09327669999999999Not shown










    







Motif 16/19










    

















    







Motif IDq-valPWM
MA0006.1_Ahr::Arnt0.0855202
NR0B1_HUMAN.H11MO.0.D0.0855202
MA0733.1_EGR40.109795
MA0472.2_EGR20.109795
MA0732.1_EGR30.109795
SUH_HUMAN.H11MO.0.A0.139311Not shown
EGR2_HUMAN.H11MO.1.A0.161148Not shown
ARNT2_HUMAN.H11MO.0.D0.194101Not shown
MA1627.1_Wt10.25656799999999996Not shown
EGR1_HUMAN.H11MO.0.A0.274515Not shown










    







Motif 17/19










    

















    







Motif IDq-valPWM
MA1135.1_FOSB::JUNB6.06775e-07
MA1138.1_FOSL2::JUNB6.06775e-07
MA1144.1_FOSL2::JUND1.15759e-06
MA0099.3_FOS::JUN1.4283600000000001e-06
MA0655.1_JDP20.000180116
FOSB_HUMAN.H11MO.0.A0.000286202Not shown
MA1132.1_JUN::JUNB0.000333507Not shown
MA1142.1_FOSL1::JUND0.000574694Not shown
MA1134.1_FOS::JUNB0.000789906Not shown
FOS_HUMAN.H11MO.0.A0.00150024Not shown










    







Motif 18/19










    

















    







Motif IDq-valPWM
MYF6_HUMAN.H11MO.0.C0.0158303
MA1635.1_BHLHE22(var.2)0.108474
MESP1_HUMAN.H11MO.0.D0.108474
ITF2_HUMAN.H11MO.0.C0.150598
HTF4_HUMAN.H11MO.0.A0.178479
ASCL2_HUMAN.H11MO.0.D0.178479Not shown
MA1619.1_Ptf1a(var.2)0.209894Not shown
ZSC31_HUMAN.H11MO.0.C0.23211700000000002Not shown
PTF1A_HUMAN.H11MO.1.B0.24179499999999998Not shown
TFE2_HUMAN.H11MO.0.A0.24179499999999998Not shown










    







Motif 19/19










    

















    







Motif IDq-valPWM
BACH2_HUMAN.H11MO.0.A0.00119261
MAFB_HUMAN.H11MO.0.B0.00208129
MA1138.1_FOSL2::JUNB0.00208129
MA1144.1_FOSL2::JUND0.00208129
JUNB_HUMAN.H11MO.0.A0.00208129
MA0099.3_FOS::JUN0.00208129Not shown
MA1135.1_FOSB::JUNB0.00208129Not shown
MA1142.1_FOSL1::JUND0.00286984Not shown
MA0462.2_BATF::JUN0.00286984Not shown
MA1634.1_BATF0.00286984Not 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")



    















    







Metacluster 1/1










    







Motif 1/19










    







Motif hCWMSeqlets










    







Motif 2/19










    







Motif hCWMSeqlets










    







Motif 3/19










    







Motif hCWMSeqlets










    







Motif 4/19










    







Motif hCWMSeqlets










    







Motif 5/19










    







Motif hCWMSeqlets










    







Motif 6/19










    







Motif hCWMSeqlets










    







Motif 7/19










    







Motif hCWMSeqlets










    







Motif 8/19










    







Motif hCWMSeqlets










    







Motif 9/19










    







Motif hCWMSeqlets










    







Motif 10/19










    







Motif hCWMSeqlets










    







Motif 11/19










    







Motif hCWMSeqlets










    







Motif 12/19










    







Motif hCWMSeqlets










    







Motif 13/19










    







Motif hCWMSeqlets










    







Motif 14/19










    







Motif hCWMSeqlets










    







Motif 15/19










    







Motif hCWMSeqlets










    







Motif 16/19










    







Motif hCWMSeqlets










    







Motif 17/19










    







Motif hCWMSeqlets










    







Motif 18/19










    







Motif hCWMSeqlets










    







Motif 19/19










    







Motif hCWMSeqlets