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: NR3C1-reddytime
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/NR3C1-reddytime_multitask_profile_fold3/NR3C1-reddytime_multitask_profile_fold3_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/NR3C1-reddytime_multitask_profile_fold3/NR3C1-reddytime_multitask_profile_fold3_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/NR3C1-reddytime_multitask_profile_fold3/NR3C1-reddytime_multitask_profile_fold3_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%|██████████| 186/186 [01:31<00:00,  2.04it/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")
if tfm_motifs_cache_dir:
    motif_hdf5.close()

Metacluster 1/2

Pattern 1/11

5110 seqlets

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

Pattern 2/11

2688 seqlets

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

Pattern 3/11

1771 seqlets

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

Pattern 4/11

1383 seqlets

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

Pattern 5/11

803 seqlets

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

Pattern 6/11

717 seqlets

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

Pattern 7/11

239 seqlets

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

Pattern 8/11

56 seqlets

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

Pattern 9/11

35 seqlets

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

Pattern 10/11

34 seqlets

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

Pattern 11/11

30 seqlets

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

Metacluster 2/2

Pattern 1/2

47 seqlets

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

Pattern 2/2

43 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

/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)
#SeqletsForwardReverse
15110
22688
31771
41383
5803
6717
7239
856
935
1034
1130

Metacluster 2/2

#SeqletsForwardReverse
147
243

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

Motif IDq-valPWM
GCR_HUMAN.H11MO.0.A2.30565e-08
PRGR_HUMAN.H11MO.0.A2.39988e-07
ANDR_HUMAN.H11MO.1.A3.22972e-06
MA0727.1_NR3C28.65976e-05
MA0113.3_NR3C10.000177774
MA0007.3_Ar0.00810513Not shown
PRGR_HUMAN.H11MO.1.A0.00812799Not shown
GCR_HUMAN.H11MO.1.A0.0234287Not shown
ESR2_HUMAN.H11MO.0.A0.30753800000000003Not shown
MA1508.1_IKZF10.30753800000000003Not shown

Motif 2/11

Motif IDq-valPWM
MA0489.1_JUN(var.2)1.21254e-07
MA1135.1_FOSB::JUNB4.9435400000000004e-06
MA1138.1_FOSL2::JUNB6.62401e-06
MA1144.1_FOSL2::JUND9.01907e-06
MA0478.1_FOSL21.08229e-05
MA1134.1_FOS::JUNB1.20254e-05Not shown
MA0099.3_FOS::JUN1.57758e-05Not shown
MA0476.1_FOS1.57758e-05Not shown
FOSB_HUMAN.H11MO.0.A1.71575e-05Not shown
FOSL1_HUMAN.H11MO.0.A1.81668e-05Not shown

Motif 3/11

Motif IDq-valPWM
CEBPB_HUMAN.H11MO.0.A3.5945100000000002e-09
CEBPA_HUMAN.H11MO.0.A1.03536e-07
CEBPD_HUMAN.H11MO.0.C2.30146e-06
MA0836.2_CEBPD3.45219e-06
MA0102.4_CEBPA4.20228e-05
MA0837.1_CEBPE0.00117285Not shown
MA0466.2_CEBPB0.0013141Not shown
MA0025.2_NFIL30.0013141Not shown
MA0838.1_CEBPG0.0019908Not shown
DBP_HUMAN.H11MO.0.B0.0023640999999999996Not shown

Motif 4/11

Motif IDq-valPWM
FOXA1_HUMAN.H11MO.0.A4.59154e-07
FOXA2_HUMAN.H11MO.0.A1.01512e-06
FOXM1_HUMAN.H11MO.0.A2.11617e-06
FOXF2_HUMAN.H11MO.0.D8.632539999999998e-06
FOXA3_HUMAN.H11MO.0.B1.72651e-05
FOXD3_HUMAN.H11MO.0.D1.72651e-05Not shown
MA0846.1_FOXC22.1410900000000002e-05Not shown
MA0847.2_FOXD23.4694200000000004e-05Not shown
FOXC1_HUMAN.H11MO.0.C4.6259e-05Not shown
FOXJ2_HUMAN.H11MO.0.C0.00108032Not shown

Motif 5/11

Motif IDq-valPWM
MA1127.1_FOSB::JUN1.2052e-05
MA1136.1_FOSB::JUNB(var.2)3.61849e-05
ATF2_HUMAN.H11MO.2.C3.61849e-05
ATF7_HUMAN.H11MO.0.D3.61849e-05
JDP2_HUMAN.H11MO.0.D3.61849e-05
MA1139.1_FOSL2::JUNB(var.2)3.61849e-05Not shown
MA1133.1_JUN::JUNB(var.2)3.61849e-05Not shown
MA0834.1_ATF74.50301e-05Not shown
MA1145.1_FOSL2::JUND(var.2)4.63167e-05Not shown
MA1140.2_JUNB(var.2)5.16975e-05Not shown

Motif 6/11

Motif IDq-valPWM
HSF2_HUMAN.H11MO.0.A0.1381
HSF1_HUMAN.H11MO.0.A0.1381
PRGR_HUMAN.H11MO.0.A0.251298
MA0808.1_TEAD30.251298
MA0771.1_HSF40.251298
HXB2_HUMAN.H11MO.0.D0.251298Not shown
HSF4_HUMAN.H11MO.0.D0.251298Not shown
MA0486.2_HSF10.251298Not shown
MA0770.1_HSF20.251298Not shown
P73_HUMAN.H11MO.1.A0.251298Not shown

Motif 7/11

Motif IDq-valPWM
FOXD2_HUMAN.H11MO.0.D0.00348996
FOXB1_HUMAN.H11MO.0.D0.110642
MA0032.2_FOXC10.237871
MA0845.1_FOXB10.237871
PO4F3_HUMAN.H11MO.0.D0.384282

Motif 8/11

Motif IDq-valPWM
CPEB1_HUMAN.H11MO.0.D8.60681e-05
PRDM6_HUMAN.H11MO.0.C0.013070400000000001
FOXL1_HUMAN.H11MO.0.D0.013070400000000001
MA1125.1_ZNF3840.013070400000000001
FOXG1_HUMAN.H11MO.0.D0.032309
FOXJ3_HUMAN.H11MO.0.A0.0567692Not shown
ANDR_HUMAN.H11MO.0.A0.06869789999999999Not shown
FUBP1_HUMAN.H11MO.0.D0.0818431Not shown
FOXJ3_HUMAN.H11MO.1.B0.0818431Not shown
MA0679.2_ONECUT10.0946926Not shown

Motif 9/11

Motif IDq-valPWM
MA0809.2_TEAD40.00639766
TEAD2_HUMAN.H11MO.0.D0.00639766
MA0090.3_TEAD10.00639766
MA1121.1_TEAD20.011840799999999999
NFIC_HUMAN.H11MO.0.A0.0241861
MA0119.1_NFIC::TLX10.027212599999999997Not shown
TEAD1_HUMAN.H11MO.0.A0.0281535Not shown
MA0808.1_TEAD30.0281535Not shown
TEAD4_HUMAN.H11MO.0.A0.0497142Not shown
MA1643.1_NFIB0.051105300000000006Not shown

Motif 10/11

Motif IDq-valPWM
MA0808.1_TEAD30.10241099999999999
MA1121.1_TEAD20.1776
MA0090.3_TEAD10.192465
TEAD1_HUMAN.H11MO.0.A0.192465
TEAD4_HUMAN.H11MO.0.A0.23768499999999998
TEAD2_HUMAN.H11MO.0.D0.389563Not shown
ZF64A_HUMAN.H11MO.0.D0.409403Not shown
NRF1_HUMAN.H11MO.0.A0.409403Not shown
MA0162.4_EGR10.409403Not shown
MA0506.1_NRF10.409403Not shown

Motif 11/11

Motif IDq-valPWM
ZN467_HUMAN.H11MO.0.C0.000322048
MAZ_HUMAN.H11MO.0.A0.000322048
ZN341_HUMAN.H11MO.0.C0.000322048
WT1_HUMAN.H11MO.0.C0.000322048
PATZ1_HUMAN.H11MO.0.C0.000322048
VEZF1_HUMAN.H11MO.0.C0.000322048Not shown
KLF15_HUMAN.H11MO.0.A0.000916191Not shown
SP3_HUMAN.H11MO.0.B0.00135164Not shown
SP1_HUMAN.H11MO.0.A0.00197712Not shown
ZN263_HUMAN.H11MO.0.A0.00354936Not shown

Metacluster 2/2

Motif 1/2

Motif IDq-valPWM
MA0478.1_FOSL20.022644099999999997
MA1144.1_FOSL2::JUND0.022644099999999997
MA1135.1_FOSB::JUNB0.022644099999999997
MA0099.3_FOS::JUN0.022644099999999997
MA1141.1_FOS::JUND0.022644099999999997
MA1138.1_FOSL2::JUNB0.022644099999999997Not shown
MA1130.1_FOSL2::JUN0.022644099999999997Not shown
MA1128.1_FOSL1::JUN0.022644099999999997Not shown
MAFB_HUMAN.H11MO.0.B0.022644099999999997Not shown
MA1633.1_BACH10.022644099999999997Not shown

Motif 2/2

Motif IDq-valPWM
PRGR_HUMAN.H11MO.0.A8.41604e-06
GCR_HUMAN.H11MO.0.A3.6594099999999996e-05
ANDR_HUMAN.H11MO.1.A5.70945e-05
MA0727.1_NR3C20.000804155
MA0113.3_NR3C10.000804155
MA0007.3_Ar0.00763375Not shown
GCR_HUMAN.H11MO.1.A0.0248251Not shown
MA1623.1_Stat20.0688146Not shown
PRGR_HUMAN.H11MO.1.A0.107423Not shown
SOX18_HUMAN.H11MO.0.D0.174873Not shown