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/ENCSR000DTO/profile_predictions.h5
DeepSHAP scores path: /users/zahoor/TF-Atlas/02-16-2021/shap/ENCSR000DTO/counts_scores_alex_format.h5
TF-MoDISco results path: /users/zahoor/TF-Atlas/02-16-2021/modisco/ENCSR000DTO/counts/modisco_results.hd5
Peaks path: /users/zahoor/TF-Atlas/data/idr_peaks//ENCFF832INR.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%|██████████| 53/53 [00:10<00:00,  5.07it/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/2










    







Pattern 1/2










    







7740 seqlets










    







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










    

















    

















    







Pattern 2/2










    







314 seqlets










    







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










    

















    

















    







Metacluster 2/2










    







Pattern 1/16










    







915 seqlets










    







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










    

















    

















    







Pattern 2/16










    







228 seqlets










    







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










    

















    

















    







Pattern 3/16










    







215 seqlets










    







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










    

















    

















    







Pattern 4/16










    







143 seqlets










    







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










    

















    

















    







Pattern 5/16










    







150 seqlets










    







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










    

















    

















    







Pattern 6/16










    







142 seqlets










    







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










    

















    

















    







Pattern 7/16










    







138 seqlets










    







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










    

















    

















    







Pattern 8/16










    







126 seqlets










    







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










    

















    

















    







Pattern 9/16










    







132 seqlets










    







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










    

















    

















    







Pattern 10/16










    







140 seqlets










    







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










    

















    

















    







Pattern 11/16










    







116 seqlets










    







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










    

















    

















    







Pattern 12/16










    







75 seqlets










    







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










    

















    

















    







Pattern 13/16










    







71 seqlets










    







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










    

















    

















    







Pattern 14/16










    







63 seqlets










    







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










    

















    

















    







Pattern 15/16










    







42 seqlets










    







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










    

















    

















    







Pattern 16/16










    







41 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/2










    







#SeqletsForwardReverse
17740
2314










    







Metacluster 2/2










    






/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
1915
2228
3215
4143
5150
6142
7138
8126
9132
10140
11116
1275
1371
1463
1542
1641























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/2










    







Motif 1/2










    

















    







Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A8.514799999999999e-17
MA0139.1_CTCF1.12174e-15
CTCFL_HUMAN.H11MO.0.A1.81263e-08
MA1102.2_CTCFL5.55329e-05
MA1568.1_TCF21(var.2)0.21925999999999998
MA1638.1_HAND20.23626Not shown
SNAI1_HUMAN.H11MO.0.C0.23626Not shown
KLF8_HUMAN.H11MO.0.C0.354638Not shown
ZIC3_HUMAN.H11MO.0.B0.354638Not shown
MA0155.1_INSM10.354638Not shown










    







Motif 2/2










    

















    







Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A0.4772
MA0139.1_CTCF0.4772
SUH_HUMAN.H11MO.0.A0.4772
CTCFL_HUMAN.H11MO.0.A0.4772
MYBA_HUMAN.H11MO.0.D0.4772










    







Metacluster 2/2










    







Motif 1/16










    

















    







Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A1.05261e-11
MA0139.1_CTCF2.85483e-10
CTCFL_HUMAN.H11MO.0.A4.49902e-08
MA1102.2_CTCFL1.10767e-05
SNAI1_HUMAN.H11MO.0.C0.21650500000000003
PLAL1_HUMAN.H11MO.0.D0.21650500000000003Not shown
KLF8_HUMAN.H11MO.0.C0.21650500000000003Not shown
MA1628.1_Zic1::Zic20.220675Not shown
MA1629.1_Zic20.22628600000000001Not shown
ZIC3_HUMAN.H11MO.0.B0.386974Not shown










    







Motif 2/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A3.24968e-05
SP3_HUMAN.H11MO.0.B0.000129701
SP1_HUMAN.H11MO.0.A0.000896146
WT1_HUMAN.H11MO.0.C0.0014397
MA0146.2_Zfx0.00194424
PATZ1_HUMAN.H11MO.0.C0.00194424Not shown
ZFX_HUMAN.H11MO.1.A0.00194424Not shown
AP2B_HUMAN.H11MO.0.B0.00254714Not shown
MA0528.2_ZNF2630.00254714Not shown
ZN263_HUMAN.H11MO.0.A0.00648162Not shown










    







Motif 3/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A1.58084e-05
SP1_HUMAN.H11MO.0.A0.00019015599999999997
SP3_HUMAN.H11MO.0.B0.00021013099999999998
THAP1_HUMAN.H11MO.0.C0.0005578790000000001
KLF3_HUMAN.H11MO.0.B0.00177451
MA1513.1_KLF150.016814Not shown
RFX1_HUMAN.H11MO.0.B0.016814Not shown
MXI1_HUMAN.H11MO.0.A0.016814Not shown
USF2_HUMAN.H11MO.0.A0.0172948Not shown
SP1_HUMAN.H11MO.1.A0.0172948Not shown










    







Motif 4/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A0.00286641
SP3_HUMAN.H11MO.0.B0.00376617
AP2D_HUMAN.H11MO.0.D0.00723837
SP1_HUMAN.H11MO.0.A0.0143343
MA1650.1_ZBTB140.0158501
USF2_HUMAN.H11MO.0.A0.0250834Not shown
ZFX_HUMAN.H11MO.1.A0.027973400000000002Not shown
THAP1_HUMAN.H11MO.0.C0.027973400000000002Not shown
HAND1_HUMAN.H11MO.1.D0.0318278Not shown
WT1_HUMAN.H11MO.0.C0.034431300000000005Not shown










    







Motif 5/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A3.88584e-08
ZN467_HUMAN.H11MO.0.C6.08063e-07
SP3_HUMAN.H11MO.0.B1.02446e-06
WT1_HUMAN.H11MO.0.C3.05397e-06
PATZ1_HUMAN.H11MO.0.C3.05397e-06
SP1_HUMAN.H11MO.0.A3.8369e-06Not shown
VEZF1_HUMAN.H11MO.0.C8.05555e-06Not shown
MAZ_HUMAN.H11MO.0.A1.37993e-05Not shown
ZN263_HUMAN.H11MO.0.A3.48542e-05Not shown
ZN341_HUMAN.H11MO.0.C4.72402e-05Not shown










    







Motif 6/16










    

















    







Motif IDq-valPWM
TBX15_HUMAN.H11MO.0.D1.4802500000000001e-06
SP2_HUMAN.H11MO.0.A1.4802500000000001e-06
SP3_HUMAN.H11MO.0.B1.4802500000000001e-06
SP1_HUMAN.H11MO.0.A1.98179e-06
MA1513.1_KLF150.000155186
ZN341_HUMAN.H11MO.0.C0.000211397Not shown
KLF16_HUMAN.H11MO.0.D0.000317736Not shown
SP4_HUMAN.H11MO.1.A0.000317736Not shown
SP1_HUMAN.H11MO.1.A0.000317736Not shown
PATZ1_HUMAN.H11MO.0.C0.0005518809999999999Not shown










    







Motif 7/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A5.37967e-06
SP1_HUMAN.H11MO.0.A4.98697e-05
SP3_HUMAN.H11MO.0.B7.44685e-05
KLF3_HUMAN.H11MO.0.B0.000851355
MA0146.2_Zfx0.00270894
MA1513.1_KLF150.00300371Not shown
SP1_HUMAN.H11MO.1.A0.00431349Not shown
NR1H4_HUMAN.H11MO.0.B0.0130644Not shown
MA1650.1_ZBTB140.0160258Not shown
PATZ1_HUMAN.H11MO.0.C0.0170095Not shown










    







Motif 8/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A0.00100312
SP1_HUMAN.H11MO.0.A0.00103263
MA0814.2_TFAP2C(var.2)0.00588679
MA1615.1_Plagl10.00588679
SP3_HUMAN.H11MO.0.B0.0103651
AP2D_HUMAN.H11MO.0.D0.0103651Not shown
MA1513.1_KLF150.0103651Not shown
RFX1_HUMAN.H11MO.0.B0.0146257Not shown
USF2_HUMAN.H11MO.0.A0.0146257Not shown
THAP1_HUMAN.H11MO.0.C0.020112Not shown










    







Motif 9/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A0.00184202
SP1_HUMAN.H11MO.0.A0.00291615
SP3_HUMAN.H11MO.0.B0.00557688
MBD2_HUMAN.H11MO.0.B0.0305682
SMAD3_HUMAN.H11MO.0.B0.0305682
RFX1_HUMAN.H11MO.0.B0.0307995Not shown
AP2D_HUMAN.H11MO.0.D0.0358041Not shown
USF2_HUMAN.H11MO.0.A0.0358041Not shown
THAP1_HUMAN.H11MO.0.C0.0358041Not shown
ZN770_HUMAN.H11MO.0.C0.0358041Not shown










    







Motif 10/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A5.2164199999999994e-06
SP3_HUMAN.H11MO.0.B0.00084602
BHA15_HUMAN.H11MO.0.B0.0015079000000000002
SP1_HUMAN.H11MO.0.A0.0015079000000000002
KLF3_HUMAN.H11MO.0.B0.00196009
MYOD1_HUMAN.H11MO.0.A0.00196009Not shown
MA1631.1_ASCL1(var.2)0.00376796Not shown
THAP1_HUMAN.H11MO.0.C0.00384747Not shown
MA1513.1_KLF150.00384747Not shown
MA0500.2_MYOG0.00384747Not shown










    







Motif 11/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A9.72001e-09
SP1_HUMAN.H11MO.0.A4.30397e-08
SP3_HUMAN.H11MO.0.B1.7779800000000001e-07
PATZ1_HUMAN.H11MO.0.C1.67235e-05
KLF16_HUMAN.H11MO.0.D1.82735e-05
TBX15_HUMAN.H11MO.0.D6.17493e-05Not shown
KLF3_HUMAN.H11MO.0.B0.000136308Not shown
SP1_HUMAN.H11MO.1.A0.00026488200000000003Not shown
WT1_HUMAN.H11MO.0.C0.000284972Not shown
ZN467_HUMAN.H11MO.0.C0.000594247Not shown










    







Motif 12/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A2.3184000000000002e-07
SP1_HUMAN.H11MO.0.A3.25995e-06
MA1513.1_KLF154.77752e-05
MA1522.1_MAZ5.5257700000000005e-05
SP3_HUMAN.H11MO.0.B5.5257700000000005e-05
SP1_HUMAN.H11MO.1.A0.000156336Not shown
MA1653.1_ZNF1480.000156336Not shown
KLF3_HUMAN.H11MO.0.B0.00050872Not shown
KLF6_HUMAN.H11MO.0.A0.000774481Not shown
ZFX_HUMAN.H11MO.1.A0.0008970139999999999Not shown










    







Motif 13/16










    

















    







Motif IDq-valPWM
RFX1_HUMAN.H11MO.0.B0.034597199999999995
SP2_HUMAN.H11MO.0.A0.034597199999999995
MBD2_HUMAN.H11MO.0.B0.08154589999999999
SP1_HUMAN.H11MO.0.A0.0842097
MA1631.1_ASCL1(var.2)0.0842097
SP3_HUMAN.H11MO.0.B0.0842097Not shown
MA0830.2_TCF40.0842097Not shown
ZBT7B_HUMAN.H11MO.0.D0.0842097Not shown
ZN331_HUMAN.H11MO.0.C0.14146199999999998Not shown
THAP1_HUMAN.H11MO.0.C0.143057Not shown










    







Motif 14/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A8.431529999999999e-07
SP3_HUMAN.H11MO.0.B9.116030000000001e-06
SP1_HUMAN.H11MO.0.A5.0752200000000005e-05
KLF16_HUMAN.H11MO.0.D0.00297567
MA1513.1_KLF150.00297567
PATZ1_HUMAN.H11MO.0.C0.00309836Not shown
MA1650.1_ZBTB140.00622493Not shown
KLF3_HUMAN.H11MO.0.B0.00645018Not shown
TBX15_HUMAN.H11MO.0.D0.00817637Not shown
WT1_HUMAN.H11MO.0.C0.00817637Not shown










    







Motif 15/16










    

















    







Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A2.2015400000000003e-06
SP3_HUMAN.H11MO.0.B3.67171e-05
SP1_HUMAN.H11MO.0.A8.37111e-05
SP1_HUMAN.H11MO.1.A0.000512659
MA1650.1_ZBTB140.000550873
MA1513.1_KLF150.0008462389999999999Not shown
THAP1_HUMAN.H11MO.0.C0.00179472Not shown
KLF16_HUMAN.H11MO.0.D0.00179472Not shown
KLF3_HUMAN.H11MO.0.B0.00983265Not shown
AP2D_HUMAN.H11MO.0.D0.010665700000000002Not shown










    







Motif 16/16










    

















    







Motif IDq-valPWM
SP1_HUMAN.H11MO.0.A0.0403729
SP2_HUMAN.H11MO.0.A0.0403729
THAP1_HUMAN.H11MO.0.C0.0652132
ZN770_HUMAN.H11MO.0.C0.0878221
MA0163.1_PLAG10.0878221
SP3_HUMAN.H11MO.0.B0.151504Not shown
GLIS3_HUMAN.H11MO.0.D0.205719Not shown
MA0528.2_ZNF2630.21796100000000002Not shown
MA1513.1_KLF150.21796100000000002Not shown
NR0B1_HUMAN.H11MO.0.D0.222372Not 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/2










    







Motif 1/2










    







Motif hCWMSeqlets










    







Motif 2/2










    







Motif hCWMSeqlets










    







Metacluster 2/2










    







Motif 1/16










    







Motif hCWMSeqlets










    







Motif 2/16










    







Motif hCWMSeqlets










    







Motif 3/16










    







Motif hCWMSeqlets










    







Motif 4/16










    







Motif hCWMSeqlets










    







Motif 5/16










    







Motif hCWMSeqlets










    







Motif 6/16










    







Motif hCWMSeqlets










    







Motif 7/16










    







Motif hCWMSeqlets










    







Motif 8/16










    







Motif hCWMSeqlets










    







Motif 9/16










    







Motif hCWMSeqlets










    







Motif 10/16










    







Motif hCWMSeqlets










    







Motif 11/16










    







Motif hCWMSeqlets










    







Motif 12/16










    







Motif hCWMSeqlets










    







Motif 13/16










    







Motif hCWMSeqlets










    







Motif 14/16










    







Motif hCWMSeqlets










    







Motif 15/16










    







Motif hCWMSeqlets










    







Motif 16/16










    







Motif hCWMSeqlets