TF-MoDISco results

Summary of motifs

TOMTOM matches to motifs

In [1]:
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

In [2]:
# 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: E2F6
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/E2F6_multitask_profile_fold2/E2F6_multitask_profile_fold2_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/E2F6_multitask_profile_fold2/E2F6_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/E2F6_multitask_profile_fold2/E2F6_multitask_profile_fold2_profile
In [3]:
# Define paths and constants
input_length = 2114
shap_score_center_size = 400
In [4]:
if tfm_motifs_cache_dir:
    os.makedirs(tfm_motifs_cache_dir, exist_ok=True)

Import SHAP scores and TF-MoDISco results

In [5]:
# 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%|██████████| 52/52 [00:36<00:00,  1.43it/s]
In [6]:
# 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

In [7]:
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

In [8]:
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")

Metacluster 1/2

Pattern 1/10

5766 seqlets

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

Pattern 2/10

3464 seqlets

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

Pattern 3/10

1880 seqlets

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

Pattern 4/10

913 seqlets

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

Pattern 5/10

325 seqlets

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

Pattern 6/10

291 seqlets

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

Pattern 7/10

160 seqlets

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

Pattern 8/10

88 seqlets

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

Pattern 9/10

64 seqlets

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

Pattern 10/10

42 seqlets

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

Metacluster 2/2

Pattern 1/1

34 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 [9]:
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

#SeqletsForwardReverse
15766
23464
31880
4913
5325
6291
7160
888
964
1042

Metacluster 2/2

#SeqletsForwardReverse
134

Top TOMTOM matches for each motif

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

In [10]:
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/10

Motif IDq-valPWM
MA0059.1_MAX::MYC0.0009135219999999999
MA0147.3_MYC0.00121075
MXI1_HUMAN.H11MO.1.A0.00121075
MAX_HUMAN.H11MO.0.A0.00121075
MXI1_HUMAN.H11MO.0.A0.00121075
MYCN_HUMAN.H11MO.0.A0.0012835999999999998Not shown
MA0058.3_MAX0.0012835999999999998Not shown
MYC_HUMAN.H11MO.0.A0.0025098000000000004Not shown
MA0825.1_MNT0.00286036Not shown
MA0104.4_MYCN0.00543002Not shown

Motif 2/10

Motif IDq-valPWM
E2F1_HUMAN.H11MO.0.A3.1025900000000004e-05
MA0471.2_E2F60.00031386900000000004
E2F3_HUMAN.H11MO.0.A0.000327077
E2F6_HUMAN.H11MO.0.A0.000367961
TFDP1_HUMAN.H11MO.0.C0.000707585
E2F4_HUMAN.H11MO.1.A0.000760058Not shown
MA0758.1_E2F70.000760058Not shown
MA0865.1_E2F80.000760058Not shown
E2F4_HUMAN.H11MO.0.A0.000760058Not shown
E2F7_HUMAN.H11MO.0.B0.00147234Not shown

Motif 3/10

No TOMTOM matches passing threshold

Motif 4/10

No TOMTOM matches passing threshold

Motif 5/10

Motif IDq-valPWM
SP1_HUMAN.H11MO.0.A0.00419218
SP2_HUMAN.H11MO.0.A0.00419218
SP3_HUMAN.H11MO.0.B0.0122296
USF2_HUMAN.H11MO.0.A0.015671
MA1650.1_ZBTB140.015671
MXI1_HUMAN.H11MO.0.A0.037386300000000004Not shown
AP2D_HUMAN.H11MO.0.D0.059696400000000004Not shown
MA1513.1_KLF150.0677144Not shown
SP1_HUMAN.H11MO.1.A0.0677144Not shown
KLF16_HUMAN.H11MO.0.D0.06842949999999999Not shown

Motif 6/10

Motif IDq-valPWM
MXI1_HUMAN.H11MO.0.A0.0111205
SP1_HUMAN.H11MO.0.A0.037448300000000004
SP2_HUMAN.H11MO.0.A0.0724522
CTCFL_HUMAN.H11MO.0.A0.0774709
MA1513.1_KLF150.0774709
USF2_HUMAN.H11MO.0.A0.08229Not shown
AP2D_HUMAN.H11MO.0.D0.08229Not shown
SP3_HUMAN.H11MO.0.B0.08229Not shown
EGR4_HUMAN.H11MO.0.D0.08229Not shown
AP2B_HUMAN.H11MO.0.B0.08229Not shown

Motif 7/10

Motif IDq-valPWM
MA1565.1_TBX183.83004e-05
MA0801.1_MGA3.83004e-05
MA0803.1_TBX153.83004e-05
MA0805.1_TBX13.83004e-05
MA1566.1_TBX30.000425627
MA1567.1_TBX60.000425627Not shown
MA0806.1_TBX40.000425627Not shown
MA0807.1_TBX50.000542603Not shown
MA0689.1_TBX200.000691998Not shown
TBX21_HUMAN.H11MO.0.A0.000863987Not shown

Motif 8/10

Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A0.00332001
THAP1_HUMAN.H11MO.0.C0.00723798
SP3_HUMAN.H11MO.0.B0.00723798
MXI1_HUMAN.H11MO.0.A0.00723798
MA1099.2_HES10.0106569
AP2D_HUMAN.H11MO.0.D0.014267700000000003Not shown
SP1_HUMAN.H11MO.0.A0.0174703Not shown
USF2_HUMAN.H11MO.0.A0.0267389Not shown
MA1650.1_ZBTB140.041357099999999994Not shown
MA0616.2_HES20.0539332Not shown

Motif 9/10

Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A0.00128491
SP1_HUMAN.H11MO.0.A0.00676468
SP3_HUMAN.H11MO.0.B0.011275
KLF3_HUMAN.H11MO.0.B0.011792
USF2_HUMAN.H11MO.0.A0.0192105
MXI1_HUMAN.H11MO.0.A0.020061199999999998Not shown
SP1_HUMAN.H11MO.1.A0.04221369999999999Not shown
MA1650.1_ZBTB140.0431712Not shown
AP2D_HUMAN.H11MO.0.D0.056854199999999994Not shown
ZFX_HUMAN.H11MO.1.A0.057782799999999995Not shown

Motif 10/10

Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A2.39978e-07
TBX15_HUMAN.H11MO.0.D2.39978e-07
SP1_HUMAN.H11MO.0.A2.39978e-07
KLF16_HUMAN.H11MO.0.D2.39978e-07
MAZ_HUMAN.H11MO.0.A2.39978e-07
SP3_HUMAN.H11MO.0.B2.39978e-07Not shown
ZN467_HUMAN.H11MO.0.C3.00163e-07Not shown
PATZ1_HUMAN.H11MO.0.C1.5985e-06Not shown
VEZF1_HUMAN.H11MO.0.C1.7018400000000002e-06Not shown
WT1_HUMAN.H11MO.0.C2.34902e-06Not shown

Metacluster 2/2

Motif 1/1

Motif IDq-valPWM
HES7_HUMAN.H11MO.0.D0.00853714
MAX_HUMAN.H11MO.0.A0.00863933
MA0058.3_MAX0.00863933
MA1648.1_TCF12(var.2)0.00863933
MA1099.2_HES10.00863933
MXI1_HUMAN.H11MO.1.A0.00899483Not shown
MA0823.1_HEY10.012939700000000002Not shown
MA0616.2_HES20.012939700000000002Not shown
MYC_HUMAN.H11MO.0.A0.012939700000000002Not shown
MA0821.1_HES50.012939700000000002Not shown