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: JUND
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/JUND_multitask_profile_fold2/JUND_multitask_profile_fold2_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/JUND_multitask_profile_fold2/JUND_multitask_profile_fold2_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/JUND_multitask_profile_fold2/JUND_multitask_profile_fold2_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%|██████████| 350/350 [04:12<00:00,  1.38it/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/8

6653 seqlets

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

Pattern 2/8

1596 seqlets

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

Pattern 3/8

972 seqlets

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

Pattern 4/8

139 seqlets

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

Pattern 5/8

89 seqlets

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

Pattern 6/8

73 seqlets

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

Pattern 7/8

68 seqlets

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

Pattern 8/8

49 seqlets

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

Metacluster 2/2

Pattern 1/8

252 seqlets

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

Pattern 2/8

163 seqlets

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

Pattern 3/8

135 seqlets

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

Pattern 4/8

65 seqlets

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

Pattern 5/8

59 seqlets

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

Pattern 6/8

57 seqlets

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

Pattern 7/8

52 seqlets

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

Pattern 8/8

46 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

#	Seqlets	Forward	Reverse
1	6653
2	1596
3	972
4	139
5	89
6	73
7	68
8	49

Metacluster 2/2

#	Seqlets	Forward	Reverse
1	252
2	163
3	135
4	65
5	59
6	57
7	52
8	46

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

Motif ID	q-val	PWM
JUN_HUMAN.H11MO.0.A	3.7376400000000005e-08
FOSL2_HUMAN.H11MO.0.A	3.7376400000000005e-08
FOSB_HUMAN.H11MO.0.A	5.00775e-07
FOSL1_HUMAN.H11MO.0.A	1.3773599999999999e-06
MA0099.3_FOS::JUN	1.3773599999999999e-06
JUND_HUMAN.H11MO.0.A	1.7217e-06	Not shown
MA1130.1_FOSL2::JUN	1.9676499999999997e-06	Not shown
MA1622.1_Smad2::Smad3	2.58254e-06	Not shown
MA1141.1_FOS::JUND	3.9001999999999995e-06	Not shown
MA1128.1_FOSL1::JUN	5.660740000000001e-06	Not shown

Motif 2/8

Motif ID	q-val	PWM
ATF2_HUMAN.H11MO.0.B	7.2941499999999994e-06
MA1133.1_JUN::JUNB(var.2)	7.2941499999999994e-06
MA0656.1_JDP2(var.2)	2.0112e-05
MA1136.1_FOSB::JUNB(var.2)	2.0112e-05
MA0840.1_Creb5	2.21767e-05
JDP2_HUMAN.H11MO.0.D	2.21767e-05	Not shown
MA0605.2_ATF3	2.21767e-05	Not shown
ATF7_HUMAN.H11MO.0.D	3.12823e-05	Not shown
MA1129.1_FOSL1::JUN(var.2)	3.53748e-05	Not shown
MA1145.1_FOSL2::JUND(var.2)	3.84061e-05	Not shown

Motif 3/8

Motif ID	q-val	PWM
MA1636.1_CEBPG(var.2)	0.000539998
ATF4_HUMAN.H11MO.0.A	0.000539998
MA0833.2_ATF4	0.000539998
MA0466.2_CEBPB	0.000539998
MA0837.1_CEBPE	0.000539998
BATF_HUMAN.H11MO.1.A	0.000539998	Not shown
CEBPG_HUMAN.H11MO.0.B	0.000539998	Not shown
CEBPD_HUMAN.H11MO.0.C	0.000668805	Not shown
DDIT3_HUMAN.H11MO.0.D	0.0009080119999999999	Not shown
CEBPB_HUMAN.H11MO.0.A	0.000952852	Not shown

Motif 4/8

Motif ID	q-val	PWM
CPEB1_HUMAN.H11MO.0.D	9.97135e-06
MA1125.1_ZNF384	0.017242700000000003
PRDM6_HUMAN.H11MO.0.C	0.0198267
FOXL1_HUMAN.H11MO.0.D	0.0198267
FOXG1_HUMAN.H11MO.0.D	0.0362718
MA0679.2_ONECUT1	0.0622139	Not shown
FOXJ3_HUMAN.H11MO.0.A	0.0760227	Not shown
ANDR_HUMAN.H11MO.0.A	0.077457	Not shown
FUBP1_HUMAN.H11MO.0.D	0.118801	Not shown
FOXJ3_HUMAN.H11MO.1.B	0.13378099999999998	Not shown

Motif 5/8

No TOMTOM matches passing threshold

Motif 6/8

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	0.000261975
SP1_HUMAN.H11MO.0.A	0.00187169
USF2_HUMAN.H11MO.0.A	0.00189819
SP3_HUMAN.H11MO.0.B	0.00247564
AP2D_HUMAN.H11MO.0.D	0.00594311
CTCFL_HUMAN.H11MO.0.A	0.00675133	Not shown
ZFX_HUMAN.H11MO.1.A	0.016271	Not shown
KLF3_HUMAN.H11MO.0.B	0.016271	Not shown
MA0146.2_Zfx	0.0165169	Not shown
SP1_HUMAN.H11MO.1.A	0.0181601	Not shown

Motif 7/8

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	8.5138e-05
SP1_HUMAN.H11MO.0.A	0.000547364
SP3_HUMAN.H11MO.0.B	0.00068174
ZFX_HUMAN.H11MO.1.A	0.000715395
MA0146.2_Zfx	0.00544758
SP1_HUMAN.H11MO.1.A	0.00677153	Not shown
KLF3_HUMAN.H11MO.0.B	0.009348100000000002	Not shown
MA1513.1_KLF15	0.0116137	Not shown
USF2_HUMAN.H11MO.0.A	0.0116137	Not shown
WT1_HUMAN.H11MO.0.C	0.0204811	Not shown

Motif 8/8

Motif ID	q-val	PWM
MAFK_HUMAN.H11MO.1.A	0.0385864
BACH1_HUMAN.H11MO.0.A	0.0385864
MAFF_HUMAN.H11MO.0.B	0.0385864
BACH2_HUMAN.H11MO.0.A	0.0385864
MAFG_HUMAN.H11MO.0.A	0.0385864
MA1633.1_BACH1	0.0385864	Not shown
MA0591.1_Bach1::Mafk	0.0385864	Not shown
JUND_HUMAN.H11MO.0.A	0.045124300000000006	Not shown
MAFB_HUMAN.H11MO.0.B	0.045124300000000006	Not shown
MAFF_HUMAN.H11MO.1.B	0.0465835	Not shown

Metacluster 2/2

Motif 1/8

No TOMTOM matches passing threshold

Motif 2/8

Motif ID	q-val	PWM
MA1132.1_JUN::JUNB	0.06871239999999999
MA0835.2_BATF3	0.06871239999999999
MA1634.1_BATF	0.06871239999999999
JUNB_HUMAN.H11MO.0.A	0.06871239999999999
MA1135.1_FOSB::JUNB	0.06871239999999999
MA1138.1_FOSL2::JUNB	0.06871239999999999	Not shown
MA1622.1_Smad2::Smad3	0.06871239999999999	Not shown
MA0462.2_BATF::JUN	0.06871239999999999	Not shown
FOSL1_HUMAN.H11MO.0.A	0.06871239999999999	Not shown
MA1134.1_FOS::JUNB	0.06871239999999999	Not shown

Motif 3/8

No TOMTOM matches passing threshold

Motif 4/8

Motif ID	q-val	PWM
COT1_HUMAN.H11MO.1.C	0.266405
MA1550.1_PPARD	0.266405
RXRA_HUMAN.H11MO.0.A	0.266405
RXRA_HUMAN.H11MO.1.A	0.266405
COT1_HUMAN.H11MO.0.C	0.266405
PPARG_HUMAN.H11MO.1.A	0.266405	Not shown
MA1149.1_RARA::RXRG	0.266405	Not shown
NR1H3_HUMAN.H11MO.0.B	0.266405	Not shown
MA1578.1_VEZF1	0.266405	Not shown
MA0504.1_NR2C2	0.266405	Not shown

Motif 5/8

Motif ID	q-val	PWM
MA0476.1_FOS	0.050478300000000004
MA1521.1_MAFA	0.050478300000000004
BATF_HUMAN.H11MO.0.A	0.050478300000000004
MA1134.1_FOS::JUNB	0.050478300000000004
MA0781.1_PAX9	0.050478300000000004
MA1144.1_FOSL2::JUND	0.050478300000000004	Not shown
MA1101.2_BACH2	0.050478300000000004	Not shown
MAFG_HUMAN.H11MO.0.A	0.050478300000000004	Not shown
MA0655.1_JDP2	0.050478300000000004	Not shown
MAFK_HUMAN.H11MO.0.A	0.050478300000000004	Not shown

Motif 6/8

Motif ID	q-val	PWM
MA0489.1_JUN(var.2)	0.038253800000000004
MAZ_HUMAN.H11MO.0.A	0.038253800000000004
MA0099.3_FOS::JUN	0.038253800000000004
JUND_HUMAN.H11MO.0.A	0.038253800000000004
FOSB_HUMAN.H11MO.0.A	0.038253800000000004
MA1144.1_FOSL2::JUND	0.038253800000000004	Not shown
MA1142.1_FOSL1::JUND	0.038253800000000004	Not shown
MA0478.1_FOSL2	0.038253800000000004	Not shown
PATZ1_HUMAN.H11MO.0.C	0.038253800000000004	Not shown
FOS_HUMAN.H11MO.0.A	0.038253800000000004	Not shown

Motif 7/8

Motif ID	q-val	PWM
PAX5_HUMAN.H11MO.0.A	0.345837
TFCP2_HUMAN.H11MO.0.D	0.345837

Motif 8/8

Motif ID	q-val	PWM
SP2_HUMAN.H11MO.0.A	0.012809000000000001
KLF3_HUMAN.H11MO.0.B	0.012809000000000001
SP1_HUMAN.H11MO.0.A	0.0147393
KLF6_HUMAN.H11MO.0.A	0.0147393
KLF1_HUMAN.H11MO.0.A	0.0229318
THAP1_HUMAN.H11MO.0.C	0.0229318	Not shown
SP3_HUMAN.H11MO.0.B	0.0446143	Not shown
ZF64A_HUMAN.H11MO.0.D	0.0452856	Not shown
P73_HUMAN.H11MO.1.A	0.0623568	Not shown
KLF5_HUMAN.H11MO.0.A	0.0623568	Not shown