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_fold10/NR3C1-reddytime_multitask_profile_fold10_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/NR3C1-reddytime_multitask_profile_fold10/NR3C1-reddytime_multitask_profile_fold10_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_fold10/NR3C1-reddytime_multitask_profile_fold10_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.03it/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

4926 seqlets

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

Pattern 2/10

2742 seqlets

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

Pattern 3/10

1959 seqlets

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

Pattern 4/10

1123 seqlets

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

Pattern 5/10

661 seqlets

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

Pattern 6/10

656 seqlets

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

Pattern 7/10

269 seqlets

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

Pattern 8/10

89 seqlets

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

Pattern 9/10

74 seqlets

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

Pattern 10/10

34 seqlets

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

Metacluster 2/2

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
14926
22742
31959
41123
5661
6656
7269
889
974
1034

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.A2.4626900000000004e-08
PRGR_HUMAN.H11MO.0.A4.1351100000000003e-07
ANDR_HUMAN.H11MO.1.A2.72731e-06
MA0727.1_NR3C20.000124712
MA0113.3_NR3C10.00026489999999999993
PRGR_HUMAN.H11MO.1.A0.00668969Not shown
MA0007.3_Ar0.012939500000000001Not shown
GCR_HUMAN.H11MO.1.A0.0342133Not shown
MA1508.1_IKZF10.338292Not shown
RARG_HUMAN.H11MO.0.B0.416871Not shown

Motif 2/10

Motif IDq-valPWM
JUN_HUMAN.H11MO.0.A1.45819e-05
FOSL2_HUMAN.H11MO.0.A1.45819e-05
MA0099.3_FOS::JUN3.5890300000000004e-05
MA1130.1_FOSL2::JUN3.5890300000000004e-05
MA1141.1_FOS::JUND3.5890300000000004e-05
NFE2_HUMAN.H11MO.0.A3.5890300000000004e-05Not shown
MA1622.1_Smad2::Smad33.5890300000000004e-05Not shown
MA0591.1_Bach1::Mafk3.5890300000000004e-05Not shown
MA1144.1_FOSL2::JUND5.27034e-05Not shown
MA1138.1_FOSL2::JUNB5.27034e-05Not shown

Motif 3/10

Motif IDq-valPWM
CEBPB_HUMAN.H11MO.0.A2.50106e-09
CEBPA_HUMAN.H11MO.0.A8.942639999999999e-08
CEBPD_HUMAN.H11MO.0.C1.5148500000000002e-06
MA0836.2_CEBPD7.821860000000001e-06
MA0102.4_CEBPA3.41003e-05
MA0025.2_NFIL30.00138365Not shown
MA0466.2_CEBPB0.00206258Not shown
MA0837.1_CEBPE0.00206258Not shown
DBP_HUMAN.H11MO.0.B0.00228343Not shown
NFIL3_HUMAN.H11MO.0.D0.00228343Not shown

Motif 4/10

Motif IDq-valPWM
FOXA1_HUMAN.H11MO.0.A7.72962e-07
FOXM1_HUMAN.H11MO.0.A8.18394e-07
FOXA2_HUMAN.H11MO.0.A1.05487e-06
FOXF2_HUMAN.H11MO.0.D6.28971e-06
FOXD3_HUMAN.H11MO.0.D1.9128900000000002e-05
MA0846.1_FOXC21.9128900000000002e-05Not shown
FOXA3_HUMAN.H11MO.0.B2.1861599999999998e-05Not shown
MA0847.2_FOXD22.6705e-05Not shown
FOXC1_HUMAN.H11MO.0.C4.90653e-05Not shown
FOXJ2_HUMAN.H11MO.0.C0.0004961480000000001Not shown

Motif 5/10

Motif IDq-valPWM
ATF2_HUMAN.H11MO.2.C2.01082e-05
ATF7_HUMAN.H11MO.0.D2.01082e-05
MA1127.1_FOSB::JUN2.01082e-05
JDP2_HUMAN.H11MO.0.D2.01082e-05
MA1133.1_JUN::JUNB(var.2)2.01082e-05
MA1145.1_FOSL2::JUND(var.2)2.68109e-05Not shown
MA1140.2_JUNB(var.2)2.9433e-05Not shown
MA1139.1_FOSL2::JUNB(var.2)2.9433e-05Not shown
MA0840.1_Creb52.9433e-05Not shown
MA0656.1_JDP2(var.2)2.9433e-05Not shown

Motif 6/10

Motif IDq-valPWM
MA0808.1_TEAD30.13655799999999998
HSF1_HUMAN.H11MO.0.A0.183821
HSF2_HUMAN.H11MO.0.A0.183821
MA0770.1_HSF20.187131
MA0486.2_HSF10.187131
TEAD4_HUMAN.H11MO.0.A0.187131Not shown
P63_HUMAN.H11MO.1.A0.187131Not shown
MA0771.1_HSF40.187131Not shown
HXB2_HUMAN.H11MO.0.D0.187131Not shown
P73_HUMAN.H11MO.1.A0.187131Not shown

Motif 7/10

Motif IDq-valPWM
FOXD2_HUMAN.H11MO.0.D0.00555494
FOXB1_HUMAN.H11MO.0.D0.0312853
MA0032.2_FOXC10.352607
MA0845.1_FOXB10.352607
MA0791.1_POU4F30.352607
MA0683.1_POU4F20.352607Not shown
PO4F3_HUMAN.H11MO.0.D0.389219Not shown
MA0847.2_FOXD20.389219Not shown
FOXJ2_HUMAN.H11MO.0.C0.389219Not shown
FOXA1_HUMAN.H11MO.0.A0.389219Not shown

Motif 8/10

Motif IDq-valPWM
CPEB1_HUMAN.H11MO.0.D9.43867e-06
FOXL1_HUMAN.H11MO.0.D0.00287734
MA1125.1_ZNF3840.00301386
FOXG1_HUMAN.H11MO.0.D0.00301386
PRDM6_HUMAN.H11MO.0.C0.010825399999999999
FOXJ3_HUMAN.H11MO.0.A0.0151388Not shown
MA0679.2_ONECUT10.0169666Not shown
HXC10_HUMAN.H11MO.0.D0.0366091Not shown
FOXJ3_HUMAN.H11MO.1.B0.038109699999999996Not shown
ANDR_HUMAN.H11MO.0.A0.038447300000000004Not shown

Motif 9/10

Motif IDq-valPWM
TBX15_HUMAN.H11MO.0.D8.1923e-05
SP3_HUMAN.H11MO.0.B8.1923e-05
ZN467_HUMAN.H11MO.0.C8.1923e-05
SP2_HUMAN.H11MO.0.A8.1923e-05
PATZ1_HUMAN.H11MO.0.C0.000130484
WT1_HUMAN.H11MO.0.C0.000135508Not shown
MAZ_HUMAN.H11MO.0.A0.000135508Not shown
KLF16_HUMAN.H11MO.0.D0.000346004Not shown
ZN263_HUMAN.H11MO.0.A0.00041926699999999997Not shown
ZN281_HUMAN.H11MO.0.A0.00041926699999999997Not shown

Motif 10/10

Motif IDq-valPWM
MA0841.1_NFE24.08909e-06
JUND_HUMAN.H11MO.0.A8.93534e-05
FOSB_HUMAN.H11MO.0.A0.00017068400000000002
MA0655.1_JDP20.000221996
NFE2_HUMAN.H11MO.0.A0.000450224
MA0478.1_FOSL20.000450224Not shown
MA1130.1_FOSL2::JUN0.000489261Not shown
JUN_HUMAN.H11MO.0.A0.000489261Not shown
MA0591.1_Bach1::Mafk0.000489261Not shown
MA1141.1_FOS::JUND0.000489261Not shown