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_fold7/REST_multitask_profile_fold7_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/REST_multitask_profile_fold7/REST_multitask_profile_fold7_profile_tfm.h5
Importance score key: profile_hyp_scores
Saved TF-MoDISco-derived motifs cache: /users/amtseng/tfmodisco/results/reports/tfmodisco_results//cache/multitask_profile/REST_multitask_profile_fold7/REST_multitask_profile_fold7_profile

# 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 [03:01<00:00,  1.58it/s]

# 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/21

3742 seqlets

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

Pattern 2/21

2573 seqlets

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

Pattern 3/21

472 seqlets

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

Pattern 4/21

387 seqlets

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

Pattern 5/21

383 seqlets

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

Pattern 6/21

346 seqlets

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

Pattern 7/21

246 seqlets

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

Pattern 8/21

239 seqlets

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

Pattern 9/21

192 seqlets

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

Pattern 10/21

161 seqlets

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

Pattern 11/21

116 seqlets

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

Pattern 12/21

79 seqlets

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

Pattern 13/21

67 seqlets

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

Pattern 14/21

51 seqlets

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

Pattern 15/21

50 seqlets

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

Pattern 16/21

44 seqlets

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

Pattern 17/21

38 seqlets

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

Pattern 18/21

34 seqlets

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

Pattern 19/21

32 seqlets

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

Pattern 20/21

31 seqlets

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

Pattern 21/21

30 seqlets

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

Metacluster 2/2

Pattern 1/4

51 seqlets

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

Pattern 2/4

42 seqlets

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

Pattern 3/4

36 seqlets

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

Pattern 4/4

33 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	3742
2	2573
3	472
4	387
5	383
6	346
7	246
8	239
9	192
10	161
11	116
12	79
13	67
14	51
15	50
16	44
17	38
18	34
19	32
20	31
21	30

Metacluster 2/2

#	Seqlets	Forward	Reverse
1	51
2	42
3	36
4	33

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

Motif ID	q-val	PWM
MA0138.2_REST	6.44492e-07
REST_HUMAN.H11MO.0.A	2.6035200000000003e-06

Motif 2/21

Motif ID	q-val	PWM
MA0138.2_REST	1.4204400000000001e-19
REST_HUMAN.H11MO.0.A	1.2278599999999999e-17

Motif 3/21

No TOMTOM matches passing threshold

Motif 4/21

No TOMTOM matches passing threshold

Motif 5/21

Motif ID	q-val	PWM
MA0139.1_CTCF	1.2968100000000002e-10
CTCF_HUMAN.H11MO.0.A	1.2968100000000002e-10
CTCFL_HUMAN.H11MO.0.A	1.7104499999999998e-10
MA1102.2_CTCFL	0.000152604
MA1568.1_TCF21(var.2)	0.11778599999999999
MA1638.1_HAND2	0.141547	Not shown
SNAI1_HUMAN.H11MO.0.C	0.250092	Not shown

Motif 6/21

Motif ID	q-val	PWM
MA0138.2_REST	0.16053399999999998
REST_HUMAN.H11MO.0.A	0.184363

Motif 7/21

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

Motif 8/21

Motif ID	q-val	PWM
REST_HUMAN.H11MO.0.A	0.27535

Motif 9/21

Motif ID	q-val	PWM
TBX20_HUMAN.H11MO.0.D	0.0610922
REST_HUMAN.H11MO.0.A	0.0610922
MA0138.2_REST	0.0610922

Motif 10/21

No TOMTOM matches passing threshold

Motif 11/21

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	0.000127348
CTCFL_HUMAN.H11MO.0.A	0.000312947
SP3_HUMAN.H11MO.0.B	0.00119845
SP1_HUMAN.H11MO.0.A	0.00130775
CTCF_HUMAN.H11MO.0.A	0.00214741
SP1_HUMAN.H11MO.1.A	0.00365823	Not shown
MA1102.2_CTCFL	0.00451772	Not shown
THAP1_HUMAN.H11MO.0.C	0.00475894	Not shown
MA0830.2_TCF4	0.00769479	Not shown
RFX1_HUMAN.H11MO.0.B	0.00781938	Not shown

Motif 12/21

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

Motif 13/21

No TOMTOM matches passing threshold

Motif 14/21

Motif ID	q-val	PWM
MA1631.1_ASCL1(var.2)	0.0807602
MA0830.2_TCF4	0.228584
MA0146.2_Zfx	0.43635799999999997
MA0139.1_CTCF	0.43635799999999997
MYOD1_HUMAN.H11MO.0.A	0.43635799999999997
SP1_HUMAN.H11MO.1.A	0.43635799999999997	Not shown

Motif 15/21

Motif ID	q-val	PWM
CPEB1_HUMAN.H11MO.0.D	9.16806e-07
MA1125.1_ZNF384	0.00194064
MA0679.2_ONECUT1	0.00385698
FOXG1_HUMAN.H11MO.0.D	0.00839897
HXC10_HUMAN.H11MO.0.D	0.00873881
FOXL1_HUMAN.H11MO.0.D	0.017896799999999997	Not shown
LMX1A_HUMAN.H11MO.0.D	0.017896799999999997	Not shown
PRDM6_HUMAN.H11MO.0.C	0.0186954	Not shown
MA0757.1_ONECUT3	0.0374379	Not shown
ONEC2_HUMAN.H11MO.0.D	0.0374379	Not shown

Motif 16/21

Motif ID	q-val	PWM
CPEB1_HUMAN.H11MO.0.D	2.8087199999999998e-05
PRDM6_HUMAN.H11MO.0.C	0.0234668
MA1125.1_ZNF384	0.0234668
FOXL1_HUMAN.H11MO.0.D	0.0234668
FOXG1_HUMAN.H11MO.0.D	0.057834300000000005
ANDR_HUMAN.H11MO.0.A	0.107599	Not shown
FOXJ3_HUMAN.H11MO.0.A	0.109215	Not shown
MA0679.2_ONECUT1	0.109215	Not shown
ONEC2_HUMAN.H11MO.0.D	0.133878	Not shown
MA0080.5_SPI1	0.15697	Not shown

Motif 17/21

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	6.193180000000001e-09
SP3_HUMAN.H11MO.0.B	1.529e-07
SP1_HUMAN.H11MO.0.A	1.59719e-07
KLF16_HUMAN.H11MO.0.D	3.35111e-05
PATZ1_HUMAN.H11MO.0.C	8.87259e-05
KLF3_HUMAN.H11MO.0.B	9.940530000000001e-05	Not shown
WT1_HUMAN.H11MO.0.C	0.00017087599999999998	Not shown
TBX15_HUMAN.H11MO.0.D	0.00021587400000000001	Not shown
SP1_HUMAN.H11MO.1.A	0.000344245	Not shown
MA1650.1_ZBTB14	0.000344245	Not shown

Motif 18/21

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

Motif 19/21

Motif ID	q-val	PWM
ZN770_HUMAN.H11MO.0.C	1.03053e-06
MA1596.1_ZNF460	0.00631478
ZN770_HUMAN.H11MO.1.C	0.0159578
MA1587.1_ZNF135	0.021223500000000003
ZN219_HUMAN.H11MO.0.D	0.141774
PITX2_HUMAN.H11MO.0.D	0.187867	Not shown
KLF6_HUMAN.H11MO.0.A	0.44492200000000004	Not shown
ZFX_HUMAN.H11MO.1.A	0.44492200000000004	Not shown
IKZF1_HUMAN.H11MO.0.C	0.44492200000000004	Not shown
ZN341_HUMAN.H11MO.0.C	0.44492200000000004	Not shown

Motif 20/21

No TOMTOM matches passing threshold

Motif 21/21

Motif ID	q-val	PWM
ZN549_HUMAN.H11MO.0.C	0.044506800000000006
REST_HUMAN.H11MO.0.A	0.288672
MA0138.2_REST	0.365429

Metacluster 2/2

Motif 1/4

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	0.00682785
SP1_HUMAN.H11MO.1.A	0.00925742
THAP1_HUMAN.H11MO.0.C	0.015552700000000001
SP3_HUMAN.H11MO.0.B	0.015552700000000001
MA1513.1_KLF15	0.015552700000000001
KLF3_HUMAN.H11MO.0.B	0.015552700000000001	Not shown
ZN335_HUMAN.H11MO.0.A	0.0287013	Not shown
KLF16_HUMAN.H11MO.0.D	0.0329694	Not shown
SP1_HUMAN.H11MO.0.A	0.0349652	Not shown
MA1650.1_ZBTB14	0.0373496	Not shown

Motif 2/4

Motif ID	q-val	PWM
MA0138.2_REST	3.53006e-05
REST_HUMAN.H11MO.0.A	3.6754699999999996e-05

Motif 3/4

Motif ID	q-val	PWM
MA0753.2_ZNF740	0.0227575
MA1650.1_ZBTB14	0.06712689999999999
MA1513.1_KLF15	0.0678762
MA0146.2_Zfx	0.0731752
KLF16_HUMAN.H11MO.0.D	0.0743157
MBD2_HUMAN.H11MO.0.B	0.12258	Not shown
SP1_HUMAN.H11MO.0.A	0.128128	Not shown
AP2B_HUMAN.H11MO.0.B	0.128128	Not shown
MA1653.1_ZNF148	0.173869	Not shown
MA1522.1_MAZ	0.173869	Not shown

Motif 4/4

Motif ID	q-val	PWM
REST_HUMAN.H11MO.0.A	6.38215e-06
MA0138.2_REST	6.38215e-06
ZN335_HUMAN.H11MO.0.A	0.0511709