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: MAX
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/MAX_multitask_profile_fold2/MAX_multitask_profile_fold2_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/MAX_multitask_profile_fold2/MAX_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/MAX_multitask_profile_fold2/MAX_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%|██████████| 204/204 [01:46<00:00,  1.92it/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

9810 seqlets

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

Pattern 2/10

827 seqlets

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

Pattern 3/10

669 seqlets

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

Pattern 4/10

443 seqlets

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

Pattern 5/10

373 seqlets

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

Pattern 6/10

314 seqlets

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

Pattern 7/10

42 seqlets

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

Pattern 8/10

39 seqlets

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

Pattern 9/10

35 seqlets

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

Pattern 10/10

31 seqlets

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

Metacluster 2/2

Pattern 1/1

201 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
19810
2827
3669
4443
5373
6314
742
839
935
1031

Metacluster 2/2

#SeqletsForwardReverse
1201

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
BMAL1_HUMAN.H11MO.0.A0.00198904
MYCN_HUMAN.H11MO.0.A0.00362214
MA0059.1_MAX::MYC0.00385158
MA0871.2_TFEC0.00409706
MYC_HUMAN.H11MO.0.A0.00556513
MA0626.1_Npas20.00556513Not shown
MXI1_HUMAN.H11MO.0.A0.00591773Not shown
MA0825.1_MNT0.00790926Not shown
MA0058.3_MAX0.00790926Not shown
MA0104.4_MYCN0.00790926Not shown

Motif 2/10

Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A3.4335500000000004e-09
MA0139.1_CTCF1.52092e-08
CTCFL_HUMAN.H11MO.0.A1.52092e-08
MA1102.2_CTCFL1.1483800000000002e-05
PLAL1_HUMAN.H11MO.0.D0.18528
KLF8_HUMAN.H11MO.0.C0.18528Not shown
MA1628.1_Zic1::Zic20.20540799999999998Not shown
ZIC3_HUMAN.H11MO.0.B0.20540799999999998Not shown
MA1548.1_PLAGL20.20540799999999998Not shown
SNAI1_HUMAN.H11MO.0.C0.20540799999999998Not shown

Motif 3/10

Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A5.4130600000000007e-05
MA0139.1_CTCF5.50594e-05
CTCFL_HUMAN.H11MO.0.A0.000145633
SNAI1_HUMAN.H11MO.0.C0.0123101
MA1102.2_CTCFL0.023805700000000003
MA1648.1_TCF12(var.2)0.0376963Not shown
MYC_HUMAN.H11MO.0.A0.0454098Not shown
BHA15_HUMAN.H11MO.0.B0.0454098Not shown
MA0830.2_TCF40.0454098Not shown
MA1568.1_TCF21(var.2)0.051375699999999996Not shown

Motif 4/10

Motif IDq-valPWM
CTCFL_HUMAN.H11MO.0.A0.00048493099999999997
CTCF_HUMAN.H11MO.0.A0.00191876
AP2B_HUMAN.H11MO.0.B0.0106125
MA0139.1_CTCF0.0106125
MA1102.2_CTCFL0.0106125
MA0155.1_INSM10.022870500000000002Not shown
MA0830.2_TCF40.0243819Not shown
INSM1_HUMAN.H11MO.0.C0.07244830000000001Not shown
MA0104.4_MYCN0.130907Not shown
MXI1_HUMAN.H11MO.0.A0.155642Not shown

Motif 5/10

Motif IDq-valPWM
CTCF_HUMAN.H11MO.0.A0.0409416
MA0139.1_CTCF0.0409416
MA1568.1_TCF21(var.2)0.0809269
CTCFL_HUMAN.H11MO.0.A0.0809269
MA0668.1_NEUROD20.258245
MA0095.2_YY10.258245Not shown
TYY1_HUMAN.H11MO.0.A0.258245Not shown
MYC_HUMAN.H11MO.0.A0.263884Not shown
MA1524.1_MSGN10.263884Not shown
SNAI1_HUMAN.H11MO.0.C0.34778400000000004Not shown

Motif 6/10

Motif IDq-valPWM
MA0139.1_CTCF6.00073e-08
CTCF_HUMAN.H11MO.0.A1.28069e-07
CTCFL_HUMAN.H11MO.0.A0.000923255
MA1102.2_CTCFL0.015489700000000002
ZIC2_HUMAN.H11MO.0.D0.41261499999999995
ZIC3_HUMAN.H11MO.0.B0.41261499999999995Not shown
MA0154.4_EBF10.45536400000000005Not shown
MA1568.1_TCF21(var.2)0.45536400000000005Not shown
SNAI1_HUMAN.H11MO.0.C0.45536400000000005Not shown
MA1638.1_HAND20.45536400000000005Not shown

Motif 7/10

Motif IDq-valPWM
CTCFL_HUMAN.H11MO.0.A0.00727025
MA0748.2_YY20.04165
CTCF_HUMAN.H11MO.0.A0.11068099999999999
MA0139.1_CTCF0.18776500000000002
USF2_HUMAN.H11MO.0.A0.229758
MA0095.2_YY10.315637Not shown
TYY1_HUMAN.H11MO.0.A0.315637Not shown
MXI1_HUMAN.H11MO.0.A0.327392Not shown
TAF1_HUMAN.H11MO.0.A0.327392Not shown
MA1108.2_MXI10.385797Not shown

Motif 8/10

Motif IDq-valPWM
AP2D_HUMAN.H11MO.0.D0.028177999999999998
MA1650.1_ZBTB140.028177999999999998
KLF1_HUMAN.H11MO.0.A0.028177999999999998
SP3_HUMAN.H11MO.0.B0.028177999999999998
SP2_HUMAN.H11MO.0.A0.028177999999999998
SP1_HUMAN.H11MO.0.A0.028177999999999998Not shown
MXI1_HUMAN.H11MO.0.A0.0492114Not shown
KLF3_HUMAN.H11MO.0.B0.051349900000000004Not shown
MA0753.2_ZNF7400.06871089999999999Not shown
EGR1_HUMAN.H11MO.0.A0.0721838Not shown

Motif 9/10

Motif IDq-valPWM
MA1108.2_MXI10.0589112
MA0104.4_MYCN0.101699
MA0668.1_NEUROD20.101699
MA1642.1_NEUROG2(var.2)0.101699
USF2_HUMAN.H11MO.0.A0.101699
CR3L1_HUMAN.H11MO.0.D0.101699Not shown
MA0093.3_USF10.126303Not shown
HEY2_HUMAN.H11MO.0.D0.13350499999999998Not shown
MA0059.1_MAX::MYC0.13350499999999998Not shown
MA0147.3_MYC0.13350499999999998Not shown

Motif 10/10

Motif IDq-valPWM
SP2_HUMAN.H11MO.0.A1.4772299999999998e-08
SP1_HUMAN.H11MO.0.A1.92912e-08
KLF16_HUMAN.H11MO.0.D1.92912e-08
SP3_HUMAN.H11MO.0.B1.9545200000000002e-07
PATZ1_HUMAN.H11MO.0.C5.135620000000001e-07
TBX15_HUMAN.H11MO.0.D8.56708e-07Not shown
MAZ_HUMAN.H11MO.0.A8.56708e-07Not shown
ZN467_HUMAN.H11MO.0.C8.56708e-07Not shown
WT1_HUMAN.H11MO.0.C1.6919999999999999e-06Not shown
KLF15_HUMAN.H11MO.0.A1.96131e-05Not shown

Metacluster 2/2

Motif 1/1

Motif IDq-valPWM
MITF_HUMAN.H11MO.0.A0.00353487
MA0825.1_MNT0.00353487
MA1464.1_ARNT20.00353487
MA1560.1_SOHLH20.00353487
BHE40_HUMAN.H11MO.0.A0.00353487
MA0626.1_Npas20.00353487Not shown
MA0692.1_TFEB0.00353487Not shown
MA1493.1_HES60.00353487Not shown
USF2_HUMAN.H11MO.0.A0.00353487Not shown
MA0664.1_MLXIPL0.00353487Not shown