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_fold8/MAFK_multitask_profile_fold8_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/MAFK_multitask_profile_fold8/MAFK_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/MAFK_multitask_profile_fold8/MAFK_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%|██████████| 311/311 [10:16<00:00,  1.98s/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/12

10884 seqlets

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

Pattern 2/12

747 seqlets

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

Pattern 3/12

280 seqlets

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

Pattern 4/12

175 seqlets

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

Pattern 5/12

170 seqlets

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

Pattern 6/12

117 seqlets

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

Pattern 7/12

89 seqlets

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

Pattern 8/12

41 seqlets

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

Pattern 9/12

36 seqlets

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

Pattern 10/12

35 seqlets

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

Pattern 11/12

35 seqlets

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

Pattern 12/12

33 seqlets

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

Metacluster 2/2

Pattern 1/1

37 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
110884
2747
3280
4175
5170
6117
789
841
936
1035
1135
1233

Metacluster 2/2

#SeqletsForwardReverse
137

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

Motif IDq-valPWM
MAFK_HUMAN.H11MO.0.A8.27531e-11
MA1520.1_MAF8.27531e-11
MAFG_HUMAN.H11MO.0.A3.41836e-10
MA0496.3_MAFK3.2461100000000005e-09
MA1521.1_MAFA7.729319999999999e-09
MAFB_HUMAN.H11MO.0.B9.995579999999999e-09Not shown
MAF_HUMAN.H11MO.0.A2.39141e-08Not shown
MAFK_HUMAN.H11MO.1.A7.0301e-08Not shown
MAFF_HUMAN.H11MO.0.B7.329160000000001e-08Not shown
MAFF_HUMAN.H11MO.1.B7.85165e-06Not shown

Motif 2/12

Motif IDq-valPWM
MA0139.1_CTCF1.6981200000000001e-18
CTCF_HUMAN.H11MO.0.A3.37681e-14
CTCFL_HUMAN.H11MO.0.A8.59287e-08
MA1102.2_CTCFL0.000165425
MA1568.1_TCF21(var.2)0.11198599999999999
MA1638.1_HAND20.11198599999999999Not shown
SNAI1_HUMAN.H11MO.0.C0.24634499999999998Not shown

Motif 3/12

Motif IDq-valPWM
MAFB_HUMAN.H11MO.0.B0.00136275
MAF_HUMAN.H11MO.1.B0.00201262
MA0659.2_MAFG0.00201262
MA0495.3_MAFF0.00201262
MAF_HUMAN.H11MO.0.A0.00201262
NRL_HUMAN.H11MO.0.D0.00201262Not shown
MA0117.2_Mafb0.00207012Not shown
MA0842.2_NRL0.00211325Not shown
MAFG_HUMAN.H11MO.0.A0.00797719Not shown
MAFK_HUMAN.H11MO.1.A0.00797719Not shown

Motif 4/12

Motif IDq-valPWM
MA0139.1_CTCF0.00251898
CTCF_HUMAN.H11MO.0.A0.00690663
CTCFL_HUMAN.H11MO.0.A0.197765
MA1568.1_TCF21(var.2)0.218982

Motif 5/12

Motif IDq-valPWM
TBX15_HUMAN.H11MO.0.D0.00575022
ZN467_HUMAN.H11MO.0.C0.00575022
PATZ1_HUMAN.H11MO.0.C0.007326999999999999
SP2_HUMAN.H11MO.0.A0.007326999999999999
ZN770_HUMAN.H11MO.0.C0.012332899999999999
WT1_HUMAN.H11MO.0.C0.012332899999999999Not shown
SP3_HUMAN.H11MO.0.B0.0156558Not shown
MA1587.1_ZNF1350.0156558Not shown
ZN341_HUMAN.H11MO.0.C0.016859299999999997Not shown
VEZF1_HUMAN.H11MO.0.C0.016859299999999997Not shown

Motif 6/12

Motif IDq-valPWM
MA0501.1_MAF::NFE21.06604e-05
MA1633.1_BACH11.06604e-05
NFE2_HUMAN.H11MO.0.A1.06604e-05
MAFK_HUMAN.H11MO.1.A1.06604e-05
BACH2_HUMAN.H11MO.0.A1.06604e-05
NF2L2_HUMAN.H11MO.0.A1.06604e-05Not shown
MA0150.2_Nfe2l21.0965e-05Not shown
MA0089.2_NFE2L11.1193399999999999e-05Not shown
MA0591.1_Bach1::Mafk1.97606e-05Not shown
BACH1_HUMAN.H11MO.0.A2.1297199999999997e-05Not shown

Motif 7/12

Motif IDq-valPWM
PAX5_HUMAN.H11MO.0.A5.07186e-05
ZN121_HUMAN.H11MO.0.C0.000266215

Motif 8/12

Motif IDq-valPWM
CPEB1_HUMAN.H11MO.0.D4.1245900000000004e-05
MA1125.1_ZNF3840.0177785
PRDM6_HUMAN.H11MO.0.C0.0201721
FOXL1_HUMAN.H11MO.0.D0.0201721
FOXG1_HUMAN.H11MO.0.D0.0201721
FOXJ3_HUMAN.H11MO.0.A0.0789709Not shown
ZFP28_HUMAN.H11MO.0.C0.0789709Not shown
MA0679.2_ONECUT10.0789709Not shown
ANDR_HUMAN.H11MO.0.A0.090932Not shown
HXC10_HUMAN.H11MO.0.D0.10733399999999998Not shown

Motif 9/12

No TOMTOM matches passing threshold

Motif 10/12

Motif IDq-valPWM
TBX15_HUMAN.H11MO.0.D1.3778800000000002e-07
MAZ_HUMAN.H11MO.0.A1.3778800000000002e-07
ZN467_HUMAN.H11MO.0.C1.3778800000000002e-07
KLF16_HUMAN.H11MO.0.D2.36253e-07
SP1_HUMAN.H11MO.0.A4.94243e-07
PATZ1_HUMAN.H11MO.0.C4.94243e-07Not shown
VEZF1_HUMAN.H11MO.0.C4.94243e-07Not shown
SP2_HUMAN.H11MO.0.A8.7154e-07Not shown
WT1_HUMAN.H11MO.0.C1.20987e-06Not shown
SP3_HUMAN.H11MO.0.B1.4359799999999998e-06Not shown

Motif 11/12

Motif IDq-valPWM
CPEB1_HUMAN.H11MO.0.D4.4106e-05
MA0679.2_ONECUT10.000742085
LMX1A_HUMAN.H11MO.0.D0.0012485999999999999
ARI3A_HUMAN.H11MO.0.D0.00209387
FOXG1_HUMAN.H11MO.0.D0.00209387
HXC10_HUMAN.H11MO.0.D0.00222469Not shown
MA1125.1_ZNF3840.00276967Not shown
FOXL1_HUMAN.H11MO.0.D0.00293222Not shown
MA0757.1_ONECUT30.019372999999999998Not shown
PO3F3_HUMAN.H11MO.0.D0.019372999999999998Not shown

Motif 12/12

Motif IDq-valPWM
CPEB1_HUMAN.H11MO.0.D5.5934300000000004e-05
MA1125.1_ZNF3840.018474900000000002
PRDM6_HUMAN.H11MO.0.C0.018474900000000002
FOXL1_HUMAN.H11MO.0.D0.0226054
FOXG1_HUMAN.H11MO.0.D0.08654939999999998
ANDR_HUMAN.H11MO.0.A0.100776Not shown
MA0679.2_ONECUT10.100776Not shown
FOXJ3_HUMAN.H11MO.0.A0.100776Not shown
ONEC2_HUMAN.H11MO.0.D0.149592Not shown
FUBP1_HUMAN.H11MO.0.D0.176137Not shown

Metacluster 2/2

Motif 1/1

Motif IDq-valPWM
MA0117.2_Mafb0.00650182
MAF_HUMAN.H11MO.1.B0.00650182
MAF_HUMAN.H11MO.0.A0.0088238
MA0659.2_MAFG0.0088238
MAFG_HUMAN.H11MO.0.A0.0088238
MAFK_HUMAN.H11MO.0.A0.0180264Not shown
MA0501.1_MAF::NFE20.03439280000000001Not shown
MA0842.2_NRL0.0558302Not shown
MAFG_HUMAN.H11MO.1.A0.0558302Not shown
MA0495.3_MAFF0.0716683Not shown