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: CEBPB
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/CEBPB_multitask_profile_fold4/CEBPB_multitask_profile_fold4_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/CEBPB_multitask_profile_fold4/CEBPB_multitask_profile_fold4_count_tfm.h5
Importance score key: count_hyp_scores
Saved TF-MoDISco-derived motifs cache: /users/amtseng/tfmodisco/results/reports/tfmodisco_results//cache/multitask_profile/CEBPB_multitask_profile_fold4/CEBPB_multitask_profile_fold4_count
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%|██████████| 273/273 [06:23<00:00,  1.41s/it]
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/9

12796 seqlets

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

Pattern 2/9

937 seqlets

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

Pattern 3/9

626 seqlets

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

Pattern 4/9

275 seqlets

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

Pattern 5/9

179 seqlets

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

Pattern 6/9

173 seqlets

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

Pattern 7/9

67 seqlets

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

Pattern 8/9

67 seqlets

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

Pattern 9/9

55 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
112796
2937
3626
4275
5179
6173
767
867
955

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

Motif IDq-valPWM
CEBPB_HUMAN.H11MO.0.A2.0903799999999997e-10
CEBPD_HUMAN.H11MO.0.C9.243450000000002e-09
CEBPA_HUMAN.H11MO.0.A1.10237e-06
MA0836.2_CEBPD1.68858e-05
MA0102.4_CEBPA0.000167612
MA0837.1_CEBPE0.000167612Not shown
MA0466.2_CEBPB0.000217813Not shown
MA0838.1_CEBPG0.000659141Not shown
MA0025.2_NFIL30.0009923430000000001Not shown
MA1636.1_CEBPG(var.2)0.00199618Not shown

Motif 2/9

Motif IDq-valPWM
FOSL1_HUMAN.H11MO.0.A1.5796999999999999e-06
JUND_HUMAN.H11MO.0.A1.8688900000000002e-06
FOSB_HUMAN.H11MO.0.A2.0275e-06
JUN_HUMAN.H11MO.0.A3.76477e-06
FOS_HUMAN.H11MO.0.A3.76477e-06
FOSL2_HUMAN.H11MO.0.A4.05499e-06Not shown
MA1130.1_FOSL2::JUN2.4601399999999998e-05Not shown
MA0478.1_FOSL22.4601399999999998e-05Not shown
MA0099.3_FOS::JUN2.4601399999999998e-05Not shown
MA1622.1_Smad2::Smad32.86675e-05Not shown

Motif 3/9

Motif IDq-valPWM
MA0139.1_CTCF5.86627e-18
CTCF_HUMAN.H11MO.0.A1.58918e-13
CTCFL_HUMAN.H11MO.0.A2.0604299999999997e-08
MA1102.2_CTCFL8.49391e-05
MA1568.1_TCF21(var.2)0.132437
MA1638.1_HAND20.171031Not shown
SNAI1_HUMAN.H11MO.0.C0.293695Not shown
ZIC3_HUMAN.H11MO.0.B0.293695Not shown
ZIC2_HUMAN.H11MO.0.D0.402451Not shown
PLAG1_HUMAN.H11MO.0.D0.402451Not shown

Motif 4/9

Motif IDq-valPWM
CTCFL_HUMAN.H11MO.0.A6.19502e-06
MA1102.2_CTCFL5.6159700000000004e-05
CTCF_HUMAN.H11MO.0.A5.63254e-05
MA0139.1_CTCF0.000373756
SP2_HUMAN.H11MO.0.A0.000926763
SP4_HUMAN.H11MO.0.A0.00477052Not shown
SP3_HUMAN.H11MO.0.B0.0052216Not shown
PLAL1_HUMAN.H11MO.0.D0.0278475Not shown
KLF15_HUMAN.H11MO.0.A0.038519599999999994Not shown
MA0830.2_TCF40.038519599999999994Not shown

Motif 5/9

Motif IDq-valPWM
GATA2_HUMAN.H11MO.0.A0.000923032
MA0036.3_GATA20.00262682
GATA1_HUMAN.H11MO.1.A0.00262682
MA0482.2_GATA40.00262682
MA1104.2_GATA60.00262682
GATA6_HUMAN.H11MO.0.A0.00262682Not shown
GATA1_HUMAN.H11MO.0.A0.00365877Not shown
TAL1_HUMAN.H11MO.0.A0.00365877Not shown
GATA4_HUMAN.H11MO.0.A0.00459859Not shown
GATA2_HUMAN.H11MO.1.A0.00929108Not shown

Motif 6/9

Motif IDq-valPWM
FOXA1_HUMAN.H11MO.0.A1.62995e-06
FOXF2_HUMAN.H11MO.0.D3.3636699999999997e-06
FOXA2_HUMAN.H11MO.0.A3.3636699999999997e-06
FOXM1_HUMAN.H11MO.0.A6.88988e-06
FOXA3_HUMAN.H11MO.0.B1.15978e-05
MA0846.1_FOXC21.65905e-05Not shown
MA0847.2_FOXD21.66746e-05Not shown
FOXD3_HUMAN.H11MO.0.D1.66746e-05Not shown
FOXC1_HUMAN.H11MO.0.C3.21743e-05Not shown
FOXD1_HUMAN.H11MO.0.D0.000251013Not shown

Motif 7/9

Motif IDq-valPWM
HNF4G_HUMAN.H11MO.0.B7.54379e-12
HNF4A_HUMAN.H11MO.0.A2.21062e-09
MA1494.1_HNF4A(var.2)0.000699863
MA0856.1_RXRG0.00103119
PPARA_HUMAN.H11MO.0.B0.00103119
MA0065.2_Pparg::Rxra0.00103119Not shown
MA0512.2_Rxra0.00103119Not shown
MA0855.1_RXRB0.00103119Not shown
MA1574.1_THRB0.00117698Not shown
MA1550.1_PPARD0.00124092Not shown

Motif 8/9

Motif IDq-valPWM
CEBPA_HUMAN.H11MO.0.A0.0007793569999999999
CEBPB_HUMAN.H11MO.0.A0.00238286
MA0102.4_CEBPA0.0045485
CEBPD_HUMAN.H11MO.0.C0.0045485
MA0836.2_CEBPD0.0045485
DBP_HUMAN.H11MO.0.B0.0539579Not shown
NFIL3_HUMAN.H11MO.0.D0.0539579Not shown
MA0025.2_NFIL30.056839Not shown
MA0043.3_HLF0.06578289999999999Not shown
CEBPG_HUMAN.H11MO.0.B0.154569Not shown

Motif 9/9

Motif IDq-valPWM
GATA6_HUMAN.H11MO.0.A0.00557649
MA1104.2_GATA60.00557649
MA0029.1_Mecom0.034264800000000005
EVI1_HUMAN.H11MO.0.B0.0425733
MA0035.4_GATA10.0960545
MA0025.2_NFIL30.10557799999999999Not shown
MA0037.3_GATA30.129746Not shown
MA0468.1_DUX40.28730900000000004Not shown
GATA3_HUMAN.H11MO.0.A0.28730900000000004Not shown
LMX1A_HUMAN.H11MO.0.D0.28730900000000004Not shown