Direct links to results

TF-MoDISco results

Summary of motifs

TOMTOM matches to motifs

import os
import sys
sys.path.append(os.path.abspath("/users/amtseng/tfmodisco/src/"))
from tfmodisco.run_tfmodisco import import_shap_scores, import_tfmodisco_results
from motif.read_motifs import pfm_info_content, pfm_to_pwm, trim_motif_by_ic
from motif.match_motifs import match_motifs_to_database
from util import figure_to_vdom_image
import plot.viz_sequence as viz_sequence
import numpy as np
import h5py
import matplotlib.pyplot as plt
import vdom.helpers as vdomh
from IPython.display import display

Define constants and paths

# Define parameters/fetch arguments
tf_name = os.environ["TFM_TF_NAME"]
shap_scores_path = os.environ["TFM_SHAP_PATH"]
tfm_results_path = os.environ["TFM_TFM_PATH"]
hyp_score_key = os.environ["TFM_HYP_SCORE_KEY"]
if "TFM_MOTIF_CACHE" in os.environ:
    tfm_motifs_cache_dir = os.environ["TFM_MOTIF_CACHE"]
else:
    tfm_motifs_cache_dir = None

print("TF name: %s" % tf_name)
print("DeepSHAP scores path: %s" % shap_scores_path)
print("TF-MoDISco results path: %s" % tfm_results_path)
print("Importance score key: %s" % hyp_score_key)
print("Saved TF-MoDISco-derived motifs cache: %s" % tfm_motifs_cache_dir)

TF name: REST
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/REST_multitask_profile_fold9/REST_multitask_profile_fold9_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/REST_multitask_profile_fold9/REST_multitask_profile_fold9_count_tfm.h5
Importance score key: count_hyp_scores
Saved TF-MoDISco-derived motifs cache: /users/amtseng/tfmodisco/results/reports/tfmodisco_results//cache/multitask_profile/REST_multitask_profile_fold9/REST_multitask_profile_fold9_count

# Define paths and constants
input_length = 2114
shap_score_center_size = 400

if tfm_motifs_cache_dir:
    os.makedirs(tfm_motifs_cache_dir, exist_ok=True)

Import SHAP scores and TF-MoDISco results

# Import SHAP coordinates and one-hot sequences
hyp_scores, _, one_hot_seqs, shap_coords = import_shap_scores(shap_scores_path, hyp_score_key, center_cut_size=shap_score_center_size)
# This cuts the sequences/scores off just as how TF-MoDISco saw them, but the coordinates are uncut

Importing SHAP scores: 100%|██████████| 287/287 [05:14<00:00,  1.09s/it]

# Import the TF-MoDISco results object
tfm_obj = 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

plot_slice = slice(int(shap_score_center_size / 4), int(3 * shap_score_center_size / 4))
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])[plot_slice], subticks_frequency=100)

Plot TF-MoDISco results

Plot all motifs by metacluster

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())
if tfm_motifs_cache_dir:
    motif_hdf5 = h5py.File(os.path.join(tfm_motifs_cache_dir, "all_motifs.h5"), "w")
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(figure_to_vdom_image(pfm_fig))
            ),
            vdomh.tr(
                vdomh.td("Hypothetical contributions (hCWM)"),
                vdomh.td(figure_to_vdom_image(hcwm_fig))
            ),
            vdomh.tr(
                vdomh.td("Actual contributions (CWM)"),
                vdomh.td(figure_to_vdom_image(cwm_fig))
            )
        )
        display(motif_table)
        plt.close("all")  # Remove all standing figures
        
        # Trim motif based on information content
        short_trimmed_pfm = trim_motif_by_ic(pfm, pfm)
        motif_pfms_short[-1].append(short_trimmed_pfm)
        
        # Expand trimming to +/- 4bp on either side
        trimmed_pfm = trim_motif_by_ic(pfm, pfm, pad=4)
        trimmed_hcwm = trim_motif_by_ic(pfm, hcwm, pad=4)
        trimmed_cwm = trim_motif_by_ic(pfm, cwm, pad=4)
        
        motif_pfms[-1].append(trimmed_pfm)
        motif_hcwms[-1].append(trimmed_hcwm)
        motif_cwms[-1].append(trimmed_cwm)
        
        num_seqlets[-1].append(len(seqlets))
        
        if tfm_motifs_cache_dir:
            # Save results and figures
            motif_id = "%d_%d" % (metacluster_i, pattern_i)
            pfm_fig.savefig(os.path.join(tfm_motifs_cache_dir, motif_id + "_pfm_full.png"))
            hcwm_fig.savefig(os.path.join(tfm_motifs_cache_dir, motif_id + "_hcwm_full.png"))
            cwm_fig.savefig(os.path.join(tfm_motifs_cache_dir, motif_id + "_cwm_full.png"))
            motif_dset = motif_hdf5.create_group(motif_id)
            motif_dset.create_dataset("pfm_full", data=pfm, compression="gzip")
            motif_dset.create_dataset("hcwm_full", data=hcwm, compression="gzip")
            motif_dset.create_dataset("cwm_full", data=cwm, compression="gzip")
            motif_dset.create_dataset("pfm_trimmed", data=trimmed_pfm, compression="gzip")
            motif_dset.create_dataset("hcwm_trimmed", data=trimmed_hcwm, compression="gzip")
            motif_dset.create_dataset("cwm_trimmed", data=trimmed_cwm, compression="gzip")
            motif_dset.create_dataset("pfm_short_trimmed", data=short_trimmed_pfm, compression="gzip")
if tfm_motifs_cache_dir:
    motif_hdf5.close()

Metacluster 1/2

Pattern 1/15

6445 seqlets

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

Pattern 2/15

3144 seqlets

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

Pattern 3/15

693 seqlets

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

Pattern 4/15

657 seqlets

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

Pattern 5/15

499 seqlets

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

Pattern 6/15

343 seqlets

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

Pattern 7/15

312 seqlets

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

Pattern 8/15

169 seqlets

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

Pattern 9/15

155 seqlets

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

Pattern 10/15

149 seqlets

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

Pattern 11/15

139 seqlets

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

Pattern 12/15

63 seqlets

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

Pattern 13/15

52 seqlets

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

Pattern 14/15

47 seqlets

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

Pattern 15/15

47 seqlets

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

Metacluster 2/2

Pattern 1/7

52 seqlets

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

Pattern 2/7

44 seqlets

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

Pattern 3/7

44 seqlets

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

Pattern 4/7

40 seqlets

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

Pattern 5/7

36 seqlets

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

Pattern 6/7

36 seqlets

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

Pattern 7/7

32 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.

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()
        
        if tfm_motifs_cache_dir:
            # Save results and figures
            motif_id = "%d_%d" % (i, j)
            f_fig.savefig(os.path.join(tfm_motifs_cache_dir, motif_id + "_hcwm_trimmed_fwd.png"))
            rc_fig.savefig(os.path.join(tfm_motifs_cache_dir, motif_id + "_hcwm_trimmed_rev.png"))

        body.append(
            vdomh.tr(
                vdomh.td(str(j + 1)),
                vdomh.td(str(num_seqlets[i][j])),
                vdomh.td(figure_to_vdom_image(f_fig)),
                vdomh.td(figure_to_vdom_image(rc_fig))
            )
        )
    display(vdomh.table(colgroup, header, vdomh.tbody(*body)))
    plt.close("all")

Metacluster 1/2

/users/amtseng/tfmodisco/src/plot/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)

#	Seqlets	Forward	Reverse
1	6445
2	3144
3	693
4	657
5	499
6	343
7	312
8	169
9	155
10	149
11	139
12	63
13	52
14	47
15	47

Metacluster 2/2

#	Seqlets	Forward	Reverse
1	52
2	44
3	44
4	40
5	36
6	36
7	32

Top TOMTOM matches for each motif

Here, the TF-MoDISco motifs are plotted as hCWMs, but the TOMTOM matches are shown as PWMs.

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
    out_dir = os.path.join(tfm_motifs_cache_dir, "tomtom", "metacluster_%d" % i) if tfm_motifs_cache_dir else None
    tomtom_matches = match_motifs_to_database(motif_pfms_short[i], top_k=num_matches_to_keep, temp_dir=out_dir)
    
    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(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(figure_to_vdom_image(fig))
                    )
                )
                if tfm_motifs_cache_dir:
                    # Save results and figures
                    motif_id = "%d_%d" % (i, j)
                    fig.savefig(os.path.join(out_dir, motif_id + ("_hit-%d.png" % (k + 1))))
            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/15

Motif ID	q-val	PWM
MA0138.2_REST	1.0142700000000001e-19
REST_HUMAN.H11MO.0.A	8.50033e-16

Motif 2/15

Motif ID	q-val	PWM
MA0138.2_REST	0.061813400000000004
REST_HUMAN.H11MO.0.A	0.0710657

Motif 3/15

Motif ID	q-val	PWM
REST_HUMAN.H11MO.0.A	4.20642e-10
MA0138.2_REST	7.758580000000001e-10

Motif 4/15

Motif ID	q-val	PWM
CTCF_HUMAN.H11MO.0.A	1.2001700000000002e-22
MA0139.1_CTCF	1.2001700000000002e-22
CTCFL_HUMAN.H11MO.0.A	9.186190000000001e-09
MA1102.2_CTCFL	0.000103559
MA1568.1_TCF21(var.2)	0.248086
SNAI1_HUMAN.H11MO.0.C	0.27015100000000003	Not shown
MA1638.1_HAND2	0.27015100000000003	Not shown
ZIC3_HUMAN.H11MO.0.B	0.415437	Not shown
MYCN_HUMAN.H11MO.0.A	0.415437	Not shown
MA1629.1_Zic2	0.415437	Not shown

Motif 5/15

No TOMTOM matches passing threshold

Motif 6/15

Motif ID	q-val	PWM
REST_HUMAN.H11MO.0.A	0.21811
MA0138.2_REST	0.492205

Motif 7/15

No TOMTOM matches passing threshold

Motif 8/15

Motif ID	q-val	PWM
MA0138.2_REST	0.0330365
REST_HUMAN.H11MO.0.A	0.0494635
MA1631.1_ASCL1(var.2)	0.217671

Motif 9/15

Motif ID	q-val	PWM
CTCF_HUMAN.H11MO.0.A	6.070949999999999e-09
MA0139.1_CTCF	6.070949999999999e-09
CTCFL_HUMAN.H11MO.0.A	1.0120200000000002e-08
MA1102.2_CTCFL	1.84858e-06
AP2B_HUMAN.H11MO.0.B	0.030818
MA1648.1_TCF12(var.2)	0.07910299999999999	Not shown
ZIC3_HUMAN.H11MO.0.B	0.101258	Not shown
MA0155.1_INSM1	0.11938900000000001	Not shown
MA1629.1_Zic2	0.11938900000000001	Not shown
BHA15_HUMAN.H11MO.0.B	0.11938900000000001	Not shown

Motif 10/15

Motif ID	q-val	PWM
MA0138.2_REST	0.24928899999999998
REST_HUMAN.H11MO.0.A	0.333723

Motif 11/15

Motif ID	q-val	PWM
CTCFL_HUMAN.H11MO.0.A	7.057e-05
CTCF_HUMAN.H11MO.0.A	0.0010845999999999998
MA0139.1_CTCF	0.0010845999999999998
MA1102.2_CTCFL	0.0010845999999999998
KLF8_HUMAN.H11MO.0.C	0.115729
AP2B_HUMAN.H11MO.0.B	0.136819	Not shown
MA1628.1_Zic1::Zic2	0.144025	Not shown
MXI1_HUMAN.H11MO.0.A	0.144025	Not shown
ZIC3_HUMAN.H11MO.0.B	0.144025	Not shown
MA1629.1_Zic2	0.144025	Not shown

Motif 12/15

Motif ID	q-val	PWM
ZN549_HUMAN.H11MO.0.C	0.07886810000000001

Motif 13/15

Motif ID	q-val	PWM
TBX20_HUMAN.H11MO.0.D	0.178121

Motif 14/15

Motif ID	q-val	PWM
TBX20_HUMAN.H11MO.0.D	0.454092
MA0499.2_MYOD1	0.454092
MYOD1_HUMAN.H11MO.0.A	0.455266
REST_HUMAN.H11MO.0.A	0.455266
MYOG_HUMAN.H11MO.0.B	0.455266
MAFA_HUMAN.H11MO.0.D	0.455266	Not shown
MA0138.2_REST	0.455266	Not shown
ZN708_HUMAN.H11MO.1.D	0.455266	Not shown

Motif 15/15

Motif ID	q-val	PWM
MA0527.1_ZBTB33	0.00768516
KAISO_HUMAN.H11MO.1.A	0.00768516
KAISO_HUMAN.H11MO.0.A	0.012855000000000002
MA0749.1_ZBED1	0.27570900000000004

Metacluster 2/2

Motif 1/7

Motif ID	q-val	PWM
TBX15_HUMAN.H11MO.0.D	4.4775200000000005e-05
KLF16_HUMAN.H11MO.0.D	0.00043737400000000003
SP2_HUMAN.H11MO.0.A	0.0006769710000000001
SP1_HUMAN.H11MO.0.A	0.00100319
SP3_HUMAN.H11MO.0.B	0.00155064
ZN148_HUMAN.H11MO.0.D	0.00228384	Not shown
MA1650.1_ZBTB14	0.00272305	Not shown
MAZ_HUMAN.H11MO.0.A	0.00302853	Not shown
PATZ1_HUMAN.H11MO.0.C	0.00302853	Not shown
KLF3_HUMAN.H11MO.0.B	0.00506225	Not shown

Motif 2/7

No TOMTOM matches passing threshold

Motif 3/7

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	0.00299075
SP1_HUMAN.H11MO.0.A	0.00299075
SP3_HUMAN.H11MO.0.B	0.00308246
ZFX_HUMAN.H11MO.1.A	0.0107729
KLF3_HUMAN.H11MO.0.B	0.09597860000000001
MXI1_HUMAN.H11MO.0.A	0.09597860000000001	Not shown
USF2_HUMAN.H11MO.0.A	0.09597860000000001	Not shown
MA1513.1_KLF15	0.09597860000000001	Not shown
SP1_HUMAN.H11MO.1.A	0.10232100000000001	Not shown
KLF16_HUMAN.H11MO.0.D	0.10232100000000001	Not shown

Motif 4/7

Motif ID	q-val	PWM
SP3_HUMAN.H11MO.0.B	0.00081129
SP1_HUMAN.H11MO.0.A	0.00308451
SP2_HUMAN.H11MO.0.A	0.00346994
TBX15_HUMAN.H11MO.0.D	0.00346994
ZN148_HUMAN.H11MO.0.D	0.00385196
WT1_HUMAN.H11MO.0.C	0.00521073	Not shown
MA1627.1_Wt1	0.00619415	Not shown
KLF16_HUMAN.H11MO.0.D	0.00619415	Not shown
WT1_HUMAN.H11MO.1.B	0.00678258	Not shown
KLF3_HUMAN.H11MO.0.B	0.00678258	Not shown

Motif 5/7

Motif ID	q-val	PWM
MAFA_HUMAN.H11MO.0.D	0.336894

Motif 6/7

Motif ID	q-val	PWM
KLF16_HUMAN.H11MO.0.D	0.114025
RARA_HUMAN.H11MO.0.A	0.114025
SP1_HUMAN.H11MO.0.A	0.114025
SRBP2_HUMAN.H11MO.0.B	0.115645
EGR2_HUMAN.H11MO.1.A	0.115645
MA0728.1_Nr2f6(var.2)	0.132548	Not shown
VEZF1_HUMAN.H11MO.0.C	0.132548	Not shown
MA0163.1_PLAG1	0.212207	Not shown
EGR1_HUMAN.H11MO.0.A	0.212207	Not shown
MA1531.1_NR1D1	0.212207	Not shown

Motif 7/7

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	0.0165939
ZN263_HUMAN.H11MO.0.A	0.0165939
MAZ_HUMAN.H11MO.0.A	0.0165939
KLF15_HUMAN.H11MO.0.A	0.027971299999999998
SP3_HUMAN.H11MO.0.B	0.027971299999999998
ZN467_HUMAN.H11MO.0.C	0.027971299999999998	Not shown
MA1652.1_ZKSCAN5	0.029420299999999996	Not shown
SP1_HUMAN.H11MO.0.A	0.030938900000000002	Not shown
WT1_HUMAN.H11MO.0.C	0.0339524	Not shown
ELF5_HUMAN.H11MO.0.A	0.037288699999999994	Not shown