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_fold8/NR3C1-reddytime_multitask_profile_fold8_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/NR3C1-reddytime_multitask_profile_fold8/NR3C1-reddytime_multitask_profile_fold8_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_fold8/NR3C1-reddytime_multitask_profile_fold8_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:30<00:00,  2.05it/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/10

5461 seqlets

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

Pattern 2/10

2631 seqlets

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

Pattern 3/10

1780 seqlets

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

Pattern 4/10

1187 seqlets

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

Pattern 5/10

691 seqlets

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

Pattern 6/10

512 seqlets

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

Pattern 7/10

246 seqlets

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

Pattern 8/10

198 seqlets

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

Pattern 9/10

46 seqlets

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

Pattern 10/10

39 seqlets

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

Metacluster 2/2

Pattern 1/2

46 seqlets

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

Pattern 2/2

45 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
15461
22631
31780
41187
5691
6512
7246
8198
946
1039

Metacluster 2/2

#SeqletsForwardReverse
146
245

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
GCR_HUMAN.H11MO.0.A4.242220000000001e-09
PRGR_HUMAN.H11MO.0.A4.8267e-07
ANDR_HUMAN.H11MO.1.A1.4495600000000002e-06
MA0727.1_NR3C20.000127595
MA0113.3_NR3C10.000250938
PRGR_HUMAN.H11MO.1.A0.00773315Not shown
MA0007.3_Ar0.0149531Not shown
GCR_HUMAN.H11MO.1.A0.032844Not shown
MA1508.1_IKZF10.356104Not shown
RARG_HUMAN.H11MO.0.B0.404974Not shown

Motif 2/10

Motif IDq-valPWM
JUN_HUMAN.H11MO.0.A2.13451e-06
FOSL1_HUMAN.H11MO.0.A2.13451e-06
MA1622.1_Smad2::Smad32.13451e-06
FOSB_HUMAN.H11MO.0.A2.13451e-06
FOSL2_HUMAN.H11MO.0.A2.33883e-06
MA1130.1_FOSL2::JUN2.33883e-06Not shown
MA0099.3_FOS::JUN2.33883e-06Not shown
MA1137.1_FOSL1::JUNB2.33883e-06Not shown
MA1141.1_FOS::JUND3.88093e-06Not shown
MA1144.1_FOSL2::JUND3.88093e-06Not shown

Motif 3/10

Motif IDq-valPWM
CEBPB_HUMAN.H11MO.0.A3.22725e-09
CEBPA_HUMAN.H11MO.0.A1.91587e-06
CEBPD_HUMAN.H11MO.0.C1.91587e-06
MA0836.2_CEBPD1.37901e-05
MA0102.4_CEBPA4.8054e-05
MA0025.2_NFIL30.0008963660000000001Not shown
MA0837.1_CEBPE0.00123493Not shown
MA0466.2_CEBPB0.00158858Not shown
MA0838.1_CEBPG0.00198173Not shown
DBP_HUMAN.H11MO.0.B0.00198173Not shown

Motif 4/10

Motif IDq-valPWM
FOXA1_HUMAN.H11MO.0.A9.36936e-07
FOXM1_HUMAN.H11MO.0.A9.36936e-07
FOXA2_HUMAN.H11MO.0.A4.66609e-06
FOXF2_HUMAN.H11MO.0.D4.66609e-06
FOXA3_HUMAN.H11MO.0.B1.67636e-05
FOXD3_HUMAN.H11MO.0.D2.0954499999999997e-05Not shown
MA0846.1_FOXC22.54851e-05Not shown
MA0847.2_FOXD23.96435e-05Not shown
FOXC1_HUMAN.H11MO.0.C3.96435e-05Not shown
FOXJ2_HUMAN.H11MO.0.C0.0008636669999999999Not shown

Motif 5/10

Motif IDq-valPWM
MA1133.1_JUN::JUNB(var.2)1.41849e-05
ATF7_HUMAN.H11MO.0.D4.5117299999999996e-05
MA1127.1_FOSB::JUN4.5117299999999996e-05
MA0656.1_JDP2(var.2)4.5117299999999996e-05
JDP2_HUMAN.H11MO.0.D4.5117299999999996e-05
MA0834.1_ATF74.73731e-05Not shown
MA1145.1_FOSL2::JUND(var.2)4.8125100000000004e-05Not shown
MA1140.2_JUNB(var.2)9.17789e-05Not shown
MA1139.1_FOSL2::JUNB(var.2)9.17789e-05Not shown
MA0840.1_Creb59.17789e-05Not shown

Motif 6/10

Motif IDq-valPWM
TEAD4_HUMAN.H11MO.0.A9.72439e-05
TEAD2_HUMAN.H11MO.0.D9.72439e-05
TEAD1_HUMAN.H11MO.0.A0.000324358
MA1121.1_TEAD20.0006829910000000001
MA0809.2_TEAD40.00373963
MA0090.3_TEAD10.00456789Not shown
MA0808.1_TEAD30.00477391Not shown
P63_HUMAN.H11MO.1.A0.10257899999999999Not shown
P73_HUMAN.H11MO.1.A0.189448Not shown
P53_HUMAN.H11MO.1.A0.217227Not shown

Motif 7/10

Motif IDq-valPWM
HSF1_HUMAN.H11MO.0.A0.165242
HSF2_HUMAN.H11MO.0.A0.165242
MA0808.1_TEAD30.165242
MA0770.1_HSF20.334488
P63_HUMAN.H11MO.1.A0.334488
MA0486.2_HSF10.334488Not shown
GCR_HUMAN.H11MO.0.A0.334488Not shown
HSF4_HUMAN.H11MO.0.D0.334488Not shown
MA0771.1_HSF40.334488Not shown
P73_HUMAN.H11MO.1.A0.334488Not shown

Motif 8/10

Motif IDq-valPWM
FOXD2_HUMAN.H11MO.0.D0.00655635
FOXB1_HUMAN.H11MO.0.D0.0536967
MA0791.1_POU4F30.255284
MA0683.1_POU4F20.255284
MA0032.2_FOXC10.255284
MA0845.1_FOXB10.259339Not shown
PO4F3_HUMAN.H11MO.0.D0.268669Not shown
FOXJ2_HUMAN.H11MO.0.C0.3754Not shown
MA0847.2_FOXD20.450721Not shown

Motif 9/10

Motif IDq-valPWM
FOS_HUMAN.H11MO.0.A0.00260271
FOSB_HUMAN.H11MO.0.A0.00334244
JUNB_HUMAN.H11MO.0.A0.00334244
FOSL1_HUMAN.H11MO.0.A0.01126
MA1137.1_FOSL1::JUNB0.01126
FOSL2_HUMAN.H11MO.0.A0.011679799999999999Not shown
MA0476.1_FOS0.011679799999999999Not shown
MA0478.1_FOSL20.0148514Not shown
MA1128.1_FOSL1::JUN0.0148514Not shown
MA1144.1_FOSL2::JUND0.0148514Not shown

Motif 10/10

Motif IDq-valPWM
MA0466.2_CEBPB0.00048093800000000004
MA0837.1_CEBPE0.00048093800000000004
MA0838.1_CEBPG0.000647584
CEBPE_HUMAN.H11MO.0.A0.00086351
MA0836.2_CEBPD0.00120099
CEBPB_HUMAN.H11MO.0.A0.00239367Not shown
MA0639.1_DBP0.00769095Not shown
MA0102.4_CEBPA0.00769095Not shown
MA0843.1_TEF0.0147944Not shown
MA0025.2_NFIL30.0160975Not shown

Metacluster 2/2

Motif 1/2

Motif IDq-valPWM
GCR_HUMAN.H11MO.0.A5.632130000000001e-06
PRGR_HUMAN.H11MO.0.A5.632130000000001e-06
ANDR_HUMAN.H11MO.1.A8.2643e-06
MA0113.3_NR3C11.27603e-05
MA0727.1_NR3C22.05008e-05
MA0007.3_Ar0.000114249Not shown
GCR_HUMAN.H11MO.1.A0.00412836Not shown
ANDR_HUMAN.H11MO.2.A0.0641511Not shown
MA1623.1_Stat20.086546Not shown
SOX18_HUMAN.H11MO.0.D0.26982399999999995Not shown

Motif 2/2

Motif IDq-valPWM
BATF_HUMAN.H11MO.1.A0.0105934
MA1143.1_FOSL1::JUND(var.2)0.034674
MA0489.1_JUN(var.2)0.034674
MA0466.2_CEBPB0.034674
FOSL2_HUMAN.H11MO.0.A0.034674
MA0837.1_CEBPE0.034674Not shown
MA0655.1_JDP20.034674Not shown
BACH2_HUMAN.H11MO.0.A0.034674Not shown
MA1132.1_JUN::JUNB0.034674Not shown
FOSL1_HUMAN.H11MO.0.A0.034674Not shown