Link to results

Within-motif heterogeneity

Clustering of peaks

import os
import h5py
import util
import moods
import viz_sequence
import numpy as np
import pandas as pd
import modisco
import sklearn.decomposition
import umap
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()

/opt/conda/lib/python3.7/site-packages/ipykernel_launcher.py:16: TqdmDeprecationWarning: This function will be removed in tqdm==5.0.0
Please use `tqdm.notebook.tqdm` instead of `tqdm.tqdm_notebook`
  app.launch_new_instance()

|<bar/>| 0/? [00:00<?, ?it/s]

# Plotting defaults
plot_params = {
    "figure.titlesize": 22,
    "axes.titlesize": 22,
    "axes.labelsize": 20,
    "legend.fontsize": 18,
    "xtick.labelsize": 16,
    "ytick.labelsize": 16,
    "font.weight": "bold"
}
plt.rcParams.update(plot_params)

Define constants and paths

# Define parameters/fetch arguments
tfm_results_path = os.environ["TFM_TFM_PATH"]
shap_scores_path = os.environ["TFM_SHAP_PATH"]
embeddings_path = os.environ["TFM_EMB_PATH"]
moods_dir = os.environ["TFM_MOODS_DIR"]
if "TFM_CLUSTER_CACHE" in os.environ:
    cluster_cache_dir = os.environ["TFM_CLUSTER_CACHE"]
else:
    cluster_cache_dir = None

print("TF-MoDISco results path: %s" % tfm_results_path)
print("DeepSHAP scores path: %s" % shap_scores_path)
print("Embeddings path: %s" % embeddings_path)
print("MOODS directory: %s" % moods_dir)
print("Saved clusters cache: %s" % cluster_cache_dir)

TF-MoDISco results path: /mydata/modisco/profile/modisco_results.h5
DeepSHAP scores path: /mydata/shap/profile_scores.h5
Embeddings path: /mydata/embeddings/embeddings.h5
MOODS directory: /mydata/reports_output/moods/profile
Saved clusters cache: /mydata/reports_output/cluster_cache_profile

# Define constants
shap_score_center_size = 400
hyp_score_key = "hyp_scores"
task_index = None

Helper functions

For plotting and organizing things

def check_tfmodisco_motif_subcluster(tfm_results):
    """
    From an imported TF-MoDISco results object, returns whether or not
    the results contain the subclustering of each motif/pattern.
    """
    metaclusters = tfm_results.metacluster_idx_to_submetacluster_results
    # Take an arbitrary metacluster
    metacluster = next(iter(metaclusters.values()))
    
    patterns = metacluster.seqlets_to_patterns_result.patterns
    # Take an arbitrary pattern
    pattern = patterns[0]
    return pattern.subclusters is not None

def compute_tfmodisco_motif_subclusters(tfm_results):
    """
    From an imported TF-MoDISco results object, computes the subclustering
    of heterogeneity within each motif/pattern.
    """
    metaclusters = tfm_results.metacluster_idx_to_submetacluster_results
    num_metaclusters = len(metaclusters.keys())
    for metacluster_i, metacluster_key in enumerate(metaclusters.keys()):
        metacluster = metaclusters[metacluster_key]
        patterns = metacluster.seqlets_to_patterns_result.patterns
        if not patterns:
            break
        num_patterns = len(patterns)
        for pattern_i, pattern in enumerate(patterns):
            # Compute subclustering for each pattern (motif)
            pattern.compute_subclusters_and_embedding(
                pattern_comparison_settings=modisco.affinitymat.core.PatternComparisonSettings(
                    track_names=["task0_hypothetical_contribs", "task0_contrib_scores"],
                    track_transformer=modisco.affinitymat.L1Normalizer(),
                    min_overlap=None  # This argument is irrelevant here
                ),
                perplexity=30, n_jobs=4, verbose=True
            )

def plot_motif_heterogeneity(tfm_results, save_dir=None):
    """
    Plots subclusters of motifs. If `save_dir` is provided, saves the
    results and figures there.
    """
    colgroup = vdomh.colgroup(
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "50%"}),
        vdomh.col(style={"width": "40%"})
    )
    header = vdomh.thead(
        vdomh.tr(
            vdomh.th("Subpattern", style={"text-align": "center"}),
            vdomh.th("Seqlets", style={"text-align": "center"}),
            vdomh.th("Embeddings", style={"text-align": "center"}),
            vdomh.th("hCWM", style={"text-align": "center"})
        )
    )
    
    if save_dir:
        motif_hdf5 = h5py.File(os.path.join(save_dir, "all_motif_subclusters.h5"), "w")
        
    metaclusters = tfm_results.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
        num_patterns = len(patterns)
        for pattern_i, pattern in enumerate(patterns):
            display(vdomh.h4("Pattern %d/%d" % (pattern_i + 1, num_patterns)))

            embedding = pattern.twod_embedding
            subpattern_clusters = pattern.subclusters

            # Aggregate motif
            pfm = pattern["sequence"].fwd
            hcwm = pattern["task0_hypothetical_contribs"].fwd
            trimmed_hcwm = util.trim_motif(pfm, hcwm, pad=4)
            hcwm_fig = viz_sequence.plot_weights(
                trimmed_hcwm, subticks_frequency=(len(trimmed_hcwm) + 1), return_fig=True
            )
            emb_fig, ax = plt.subplots()
            ax.scatter(
                embedding[:,0], embedding[:,1], c=subpattern_clusters, cmap="tab20", alpha=0.3
            )

            table_rows = [vdomh.tr(
                vdomh.td("Agg."),
                vdomh.td(str(len(pattern.seqlets))),
                vdomh.td(util.figure_to_vdom_image(emb_fig)),
                vdomh.td(util.figure_to_vdom_image(hcwm_fig))
            )]
            
            if save_dir:
                # Save aggregate embedding plot
                motif_id = "%d_%d" % (metacluster_i, pattern_i)
                emb_fig.savefig(os.path.join(save_dir, motif_id + "_subcluster_agg.png"))
                
                # Create dictionaries for subclusters
                sc_emb_figs, sc_hcwm_figs, sc_motifs = {}, {}, {}

            for subpattern_key, subpattern in pattern.subcluster_to_subpattern.items():
                pfm = subpattern["sequence"].fwd
                cwm = subpattern["task0_contrib_scores"].fwd
                hcwm = subpattern["task0_hypothetical_contribs"].fwd
                trimmed_hcwm = util.trim_motif(pfm, hcwm, pad=4)
                hcwm_fig = viz_sequence.plot_weights(
                    trimmed_hcwm, subticks_frequency=(len(trimmed_hcwm) + 1), return_fig=True
                )
                emb_fig, ax = plt.subplots()
                ax.scatter(
                    embedding[:,0], embedding[:,1], c=(subpattern_clusters == subpattern_key), alpha=0.3
                )

                table_rows.append(vdomh.tr(
                    vdomh.td(str(subpattern_key)),
                    vdomh.td(str(len(subpattern.seqlets))),
                    vdomh.td(util.figure_to_vdom_image(emb_fig)),
                    vdomh.td(util.figure_to_vdom_image(hcwm_fig))
                ))
                
                if save_dir:
                    sc_emb_figs[subpattern_key] = emb_fig
                    sc_hcwm_figs[subpattern_key] = hcwm_fig
                    sc_motifs[subpattern_key] = (pfm, cwm, hcwm, trimmed_hcwm)
                
            if save_dir:
                # Save embedding plots, hCWM figure, and motifs
                for sc_key, emb_fig in sc_emb_figs.items():
                    emb_fig.savefig(os.path.join(save_dir, motif_id + ("_subcluster_%s.png" % sc_key)))
                for sc_key, hcwm_fig in sc_hcwm_figs.items():
                    hcwm_fig.savefig(os.path.join(save_dir, motif_id + ("_subcluster_%s_hcwm_trimmed.png" % sc_key)))

                motif_dset = motif_hdf5.create_group(motif_id)
                for sc_key, (pfm, cwm, hcwm, trimmed_hcwm) in sc_motifs.items():
                    sc_dset = motif_dset.create_group("subcluster_%s" % sc_key)
                    sc_dset.create_dataset("pfm_full", data=pfm, compression="gzip")
                    sc_dset.create_dataset("hcwm_full", data=hcwm, compression="gzip")
                    sc_dset.create_dataset("cwm_full", data=cwm, compression="gzip")
                    sc_dset.create_dataset("hcwm_trimmed", data=trimmed_hcwm, compression="gzip")

            table = vdomh.table(header, vdomh.tbody(*table_rows))
            display(table)
            plt.close("all")  # Remove all standing figures

def import_tfmodisco_motifs(tfm_results_path, trim=True, only_pos=True):
    """
    Imports hCWMs to into a dictionary, mapping `(x, y)` to the hCWM,
    where `x` is the metacluster index and `y` is the pattern index.
    Arguments:
        `tfm_results_path`: path to HDF5 containing TF-MoDISco results
        `out_dir`: where to save motifs
        `trim`: if True, trim the motif flanks based on total importance
        `only_pos`: if True, only return motifs with positive contributions
    Returns the dictionary of hCWM.
    """ 
    hcwms = {}
    with h5py.File(tfm_results_path, "r") as f:
        metaclusters = f["metacluster_idx_to_submetacluster_results"]
        num_metaclusters = len(metaclusters.keys())
        for metacluster_i, metacluster_key in enumerate(metaclusters.keys()):
            metacluster = metaclusters[metacluster_key]
            if "patterns" not in metacluster["seqlets_to_patterns_result"]:
                continue
            patterns = metacluster["seqlets_to_patterns_result"]["patterns"]
            num_patterns = len(patterns["all_pattern_names"][:])
            for pattern_i, pattern_name in enumerate(patterns["all_pattern_names"][:]):
                pattern_name = pattern_name.decode()
                pattern = patterns[pattern_name]
                pfm = pattern["sequence"]["fwd"][:]
                hcwm = pattern["task0_hypothetical_contribs"]["fwd"][:]
                cwm = pattern["task0_contrib_scores"]["fwd"][:]
                
                # Check that the contribution scores are overall positive
                if only_pos and np.sum(cwm) < 0:
                    continue
                    
                if trim:
                    hcwm = util.trim_motif(pfm, hcwm, pad=4)
                    
                hcwms["%d_%d" % (metacluster_i,pattern_i)] = hcwm
    return hcwms

def get_hit_peak_indices(hit_table, motif_keys):
    """
    Returns a dictionary of NumPy arrays, mapping each motif key to
    the set of peak indices that contain that motif.
    """
    hit_peak_indices = {}
    for motif_key in motif_keys:
        hit_peak_indices[motif_key] = hit_table[hit_table["key"] == motif_key]["peak_index"].values
    return hit_peak_indices

def umap_transform(matrix):
    """
    Converts N x D matrix into transformed N x 2 matrix using
    UMAP. First projects down to 50 components using PCA.
    """
    # First reduce using PCA
    centered = matrix - np.mean(matrix, axis=0, keepdims=True)
    pca = sklearn.decomposition.PCA(n_components=50)
    reduced = pca.fit_transform(centered)

    # Run UMAP
    um = umap.UMAP(verbose=False)
    return um.fit_transform(centered)

def plot_peak_clustering(embeddings_path, motif_keys, hcwms, hit_emb_indices, save_dir=None):
    """
    Plots clustering of peaks based on motif hits. If `save_dir` is provided,
    saves the results and figures there.
    """
    colgroup = vdomh.colgroup(
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "55"}),
        vdomh.col(style={"width": "40%"})
    )
    header = vdomh.thead(
        vdomh.tr(
            vdomh.th("Motif key", style={"text-align": "center"}),
            vdomh.th("Embeddings", style={"text-align": "center"}),
            vdomh.th("hCWM", style={"text-align": "center"})
        )
    )
    
    if save_dir:
        emb_hdf5 = h5py.File(os.path.join(save_dir, "transformed_embeddings.h5"), "w")
    
    embeddings_reader = h5py.File(embeddings_path, "r")
    num_layers = embeddings_reader["embeddings"]["mean"].shape[1]
    
    def plot_embeddings(embeddings, layer_id):
        umap_trans = umap_transform(embeddings)
        
        if save_dir:
            emb_hdf5.create_dataset(
                "layer_%s" % layer_id, data=umap_trans, compression="gzip"
            )
        
        table_rows = []
        for motif_key in motif_keys:
            hcwm = hcwms[motif_key]
            hcwm_fig = viz_sequence.plot_weights(
                hcwm, subticks_frequency=(len(hcwm) + 1), return_fig=True
            )
            emb_fig, ax = plt.subplots()
            subset = np.zeros(len(embeddings), dtype=int)
            subset[hit_emb_indices[motif_key]] = 1
            ax.scatter(
                umap_trans[:,0], umap_trans[:,1], c=subset, alpha=0.1
            )

            table_rows.append(vdomh.tr(
                vdomh.td(motif_key),
                vdomh.td(util.figure_to_vdom_image(emb_fig)),
                vdomh.td(util.figure_to_vdom_image(hcwm_fig))
            ))
            
            if save_dir:
                emb_fig.savefig(
                    os.path.join(save_dir, "layer_%s_motif_%s_peak_embeddings.png" % (layer_id, motif_key))
                )

        table = vdomh.table(header, vdomh.tbody(*table_rows))
        display(table)
        plt.close("all")  # Remove all standing figures
    
    for i in range(num_layers):
        display(vdomh.h3("Layer %d/%d" % (i + 1, num_layers)))
        
        embeddings = np.concatenate([
            embeddings_reader["embeddings"]["mean"][:, i],
            embeddings_reader["embeddings"]["std"][:, i],
            embeddings_reader["embeddings"]["max"][:, i],
            embeddings_reader["embeddings"]["min"][:, i]
        ], axis=1)  # Shape: N x (F * 4)
        
        plot_embeddings(embeddings, i)
        
    # Concatenate of all layers
    embeddings = np.concatenate([
        embeddings_reader["embeddings"]["mean"][:],
        embeddings_reader["embeddings"]["std"][:],
        embeddings_reader["embeddings"]["max"][:],
        embeddings_reader["embeddings"]["min"][:]
    ], axis=2)  # Shape: N x C x (F * 4)
    embeddings = np.reshape(embeddings, (embeddings.shape[0], -1))
    display(vdomh.h3("All layers"))
    plot_embeddings(embeddings, "all")
    
    embeddings_reader.close()

Import TF-MoDISco results

Run motif subclustering

# 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%|██████████| 11/11 [00:01<00:00,  7.08it/s]

# Import the TF-MoDISco results object
tfm_obj = util.import_tfmodisco_results(tfm_results_path, hyp_scores, one_hot_seqs, shap_score_center_size)

# Compute subclusters (needed for older versions of TF-MoDISco); this takes awhile!
if not check_tfmodisco_motif_subcluster(tfm_obj):
    compute_tfmodisco_motif_subclusters(tfm_obj)

Import motif hits

For each motif, determine the peaks that contain it

# Import the hCWMs
hcwms = import_tfmodisco_motifs(tfm_results_path)
motif_keys = list(hcwms.keys())

# Import the motif hits
hit_table = moods.import_moods_hits(os.path.join(moods_dir, "moods_filtered_collapsed_scored_thresholded.bed"))

hit_peak_indices = get_hit_peak_indices(hit_table, motif_keys)

if cluster_cache_dir:
    # Save hit embedding indices
    with h5py.File(os.path.join(cluster_cache_dir, "motif_to_embedding_inds.h5"), "w") as f:
        for motif_key, emb_inds in hit_peak_indices.items():
            f.create_dataset(motif_key, data=emb_inds, compression="gzip")

Within-motif heterogeneity

For each motif, show the subclusters that exist within the TF-MoDISco-identified subpatterns

plot_motif_heterogeneity(tfm_obj, cluster_cache_dir)

Metacluster 1/2

Pattern 1/5

/scratch/TF-Atlas/pipeline/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)

Subpattern	Seqlets	Embeddings	hCWM
Agg.	9250
0	1825
1	1824
2	1314
3	1192
4	1037
5	767
6	735
7	221
8	202
9	122
10	8
11	3

Pattern 2/5

Subpattern	Seqlets	Embeddings	hCWM
Agg.	1364
0	546
1	508
2	110
3	105
4	95

Pattern 3/5

Subpattern	Seqlets	Embeddings	hCWM
Agg.	208
0	69
1	43
2	31
3	20
4	17
5	17
6	11

Pattern 4/5

Subpattern	Seqlets	Embeddings	hCWM
Agg.	166
0	44
1	42
2	29
3	25
5	13
4	13

Pattern 5/5

Subpattern	Seqlets	Embeddings	hCWM
Agg.	133
0	33
1	29
2	22
3	20
4	18
5	7
6	4

Metacluster 2/2

Pattern 1/1

Subpattern	Seqlets	Embeddings	hCWM
Agg.	234
0	55
1	48
2	27
3	26
4	23
5	19
7	18
6	18

Peak clustering

For each peak, cluster the peaks by embeddings to highlight the structure of different peaks and different motifs

plot_peak_clustering(embeddings_path, motif_keys, hcwms, hit_peak_indices, cluster_cache_dir)

Layer 1/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 2/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 3/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 4/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 5/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 6/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 7/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

Layer 8/8

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0

All layers

Motif key	Embeddings	hCWM
0_0
0_1
0_2
0_3
0_4
1_0