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: MAFK
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/MAFK_multitask_profile_fold4/MAFK_multitask_profile_fold4_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/MAFK_multitask_profile_fold4/MAFK_multitask_profile_fold4_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/MAFK_multitask_profile_fold4/MAFK_multitask_profile_fold4_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%|██████████| 311/311 [06:32<00:00,  1.26s/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

10569 seqlets

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

Pattern 2/9

856 seqlets

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

Pattern 3/9

337 seqlets

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

Pattern 4/9

236 seqlets

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

Pattern 5/9

77 seqlets

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

Pattern 6/9

61 seqlets

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

Pattern 7/9

49 seqlets

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

Pattern 8/9

35 seqlets

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

Pattern 9/9

33 seqlets

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

Metacluster 2/2

Pattern 1/1

64 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
110569
2856
3337
4236
577
661
749
835
933

Metacluster 2/2

#SeqletsForwardReverse
164

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
MA0496.3_MAFK3.3744600000000004e-09
MAFK_HUMAN.H11MO.0.A3.49464e-09
MAFG_HUMAN.H11MO.0.A5.26031e-09
MA1520.1_MAF6.024390000000001e-09
MAFB_HUMAN.H11MO.0.B1.12617e-08
MAFK_HUMAN.H11MO.1.A1.3956999999999999e-08Not shown
MAF_HUMAN.H11MO.0.A4.57935e-08Not shown
MA1521.1_MAFA6.7919e-08Not shown
MAFF_HUMAN.H11MO.0.B1.59208e-07Not shown
MAF_HUMAN.H11MO.1.B1.17725e-05Not shown

Motif 2/9

Motif IDq-valPWM
MA0139.1_CTCF1.69182e-16
CTCF_HUMAN.H11MO.0.A5.07102e-12
CTCFL_HUMAN.H11MO.0.A5.06601e-07
MA1102.2_CTCFL0.000273948
MA1568.1_TCF21(var.2)0.0804931
MA1638.1_HAND20.09571110000000001Not shown
SNAI1_HUMAN.H11MO.0.C0.0958145Not shown
ZIC3_HUMAN.H11MO.0.B0.288598Not shown
MA0155.1_INSM10.379029Not shown
BHA15_HUMAN.H11MO.0.B0.379029Not shown

Motif 3/9

Motif IDq-valPWM
MA0117.2_Mafb1.34017e-05
MAFK_HUMAN.H11MO.0.A1.34017e-05
MAFG_HUMAN.H11MO.0.A0.00158915
MA0659.2_MAFG0.00158915
MA0495.3_MAFF0.00158915
MA0501.1_MAF::NFE20.00740178Not shown
MAFF_HUMAN.H11MO.0.B0.012844Not shown
NF2L2_HUMAN.H11MO.0.A0.017558Not shown
MA0842.2_NRL0.017558Not shown
MA0496.3_MAFK0.038749599999999995Not shown

Motif 4/9

Motif IDq-valPWM
ZN770_HUMAN.H11MO.0.C3.30125e-05
ZFX_HUMAN.H11MO.1.A0.00629324
EGR1_HUMAN.H11MO.0.A0.029957099999999997
SP1_HUMAN.H11MO.0.A0.029957099999999997
TBX15_HUMAN.H11MO.0.D0.0498532
ZN341_HUMAN.H11MO.0.C0.0557474Not shown
SP2_HUMAN.H11MO.0.A0.0651404Not shown
KLF15_HUMAN.H11MO.0.A0.06630010000000001Not shown
EGR2_HUMAN.H11MO.1.A0.0874312Not shown
WT1_HUMAN.H11MO.0.C0.0874312Not shown

Motif 5/9

Motif IDq-valPWM
MAFK_HUMAN.H11MO.1.A0.0113589
MAF_HUMAN.H11MO.0.A0.011506
MAF_HUMAN.H11MO.1.B0.0121694
MA1520.1_MAF0.0161153
MA0496.3_MAFK0.0161153
MA1521.1_MAFA0.0161153Not shown
MAFG_HUMAN.H11MO.0.A0.0161153Not shown
MAFF_HUMAN.H11MO.0.B0.0161153Not shown
MAFK_HUMAN.H11MO.0.A0.018639700000000002Not shown
MAFB_HUMAN.H11MO.0.B0.018639700000000002Not shown

Motif 6/9

Motif IDq-valPWM
MA0139.1_CTCF0.0201361
CTCF_HUMAN.H11MO.0.A0.0201361
CTCFL_HUMAN.H11MO.0.A0.103852
ZN436_HUMAN.H11MO.0.C0.18434
SNAI1_HUMAN.H11MO.0.C0.195724
MA1638.1_HAND20.24803899999999998Not shown
PLAG1_HUMAN.H11MO.0.D0.48714700000000005Not shown
ZIC3_HUMAN.H11MO.0.B0.48714700000000005Not shown
SCRT1_HUMAN.H11MO.0.D0.48714700000000005Not shown
ZIC2_HUMAN.H11MO.0.D0.48714700000000005Not shown

Motif 7/9

Motif IDq-valPWM
MA0139.1_CTCF0.17189200000000002
SCRT1_HUMAN.H11MO.0.D0.22181599999999999
MA0743.2_SCRT10.22181599999999999
MA0820.1_FIGLA0.22181599999999999
MA0744.2_SCRT20.22181599999999999
CTCF_HUMAN.H11MO.0.A0.23686999999999997Not shown

Motif 8/9

Motif IDq-valPWM
BACH2_HUMAN.H11MO.0.A0.00019177900000000002
MAFK_HUMAN.H11MO.0.A0.00026073
MAFG_HUMAN.H11MO.0.A0.00026073
MA0501.1_MAF::NFE20.00026073
NF2L2_HUMAN.H11MO.0.A0.00026073
MAFF_HUMAN.H11MO.0.B0.00029179Not shown
BACH1_HUMAN.H11MO.0.A0.0008593030000000001Not shown
JUNB_HUMAN.H11MO.0.A0.0010939Not shown
MAFK_HUMAN.H11MO.1.A0.0013861Not shown
MA0659.2_MAFG0.0013861Not shown

Motif 9/9

Motif IDq-valPWM
LMX1A_HUMAN.H11MO.0.D0.00304686
FOXD2_HUMAN.H11MO.0.D0.00304686
PO3F3_HUMAN.H11MO.0.D0.0186408
HXC10_HUMAN.H11MO.0.D0.038788199999999995
FOXG1_HUMAN.H11MO.0.D0.0403047
MA0845.1_FOXB10.0445086Not shown
LMX1B_HUMAN.H11MO.0.D0.0534067Not shown
MA0679.2_ONECUT10.0534067Not shown
MA0681.2_PHOX2B0.0534067Not shown
CPEB1_HUMAN.H11MO.0.D0.0551131Not shown

Metacluster 2/2

Motif 1/1

Motif IDq-valPWM
MAFG_HUMAN.H11MO.0.A2.5690300000000003e-06
MAFK_HUMAN.H11MO.0.A2.5690300000000003e-06
MA0659.2_MAFG7.096360000000001e-06
MA0495.3_MAFF5.12943e-05
MAFF_HUMAN.H11MO.0.B0.000168145
MAF_HUMAN.H11MO.0.A0.000168145Not shown
MA0117.2_Mafb0.000210015Not shown
MA0496.3_MAFK0.000766066Not shown
MAF_HUMAN.H11MO.1.B0.0008327780000000001Not shown
MAFK_HUMAN.H11MO.1.A0.00107331Not shown