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: FOXA2
DeepSHAP scores path: /users/amtseng/tfmodisco/results/importance_scores/multitask_profile/FOXA2_multitask_profile_fold1/FOXA2_multitask_profile_fold1_imp_scores.h5
TF-MoDISco results path: /users/amtseng/tfmodisco/results/tfmodisco/multitask_profile/FOXA2_multitask_profile_fold1/FOXA2_multitask_profile_fold1_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/FOXA2_multitask_profile_fold1/FOXA2_multitask_profile_fold1_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%|██████████| 174/174 [02:25<00:00,  1.19it/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")

Metacluster 1/2

Pattern 1/10

9470 seqlets

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

Pattern 2/10

1284 seqlets

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

Pattern 3/10

909 seqlets

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

Pattern 4/10

536 seqlets

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

Pattern 5/10

302 seqlets

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

Pattern 6/10

197 seqlets

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

Pattern 7/10

153 seqlets

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

Pattern 8/10

143 seqlets

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

Pattern 9/10

131 seqlets

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

Pattern 10/10

54 seqlets

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

Metacluster 2/2

Pattern 1/6

64 seqlets

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

Pattern 2/6

51 seqlets

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

Pattern 3/6

47 seqlets

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

Pattern 4/6

46 seqlets

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

Pattern 5/6

42 seqlets

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

Pattern 6/6

35 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
19470
21284
3909
4536
5302
6197
7153
8143
9131
1054

Metacluster 2/2

#SeqletsForwardReverse
164
251
347
446
542
635

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
FOXA1_HUMAN.H11MO.0.A4.9867e-10
FOXA2_HUMAN.H11MO.0.A9.59807e-09
FOXA3_HUMAN.H11MO.0.B3.1536299999999995e-08
MA0846.1_FOXC25.3656499999999995e-08
FOXD3_HUMAN.H11MO.0.D1.00896e-06
FOXF2_HUMAN.H11MO.0.D1.00896e-06Not shown
MA0847.2_FOXD21.00896e-06Not shown
MA0032.2_FOXC13.5341900000000005e-06Not shown
FOXC1_HUMAN.H11MO.0.C3.6594499999999998e-06Not shown
FOXM1_HUMAN.H11MO.0.A4.08081e-05Not shown

Motif 2/10

Motif IDq-valPWM
HNF4G_HUMAN.H11MO.0.B1.7960000000000002e-07
HNF4A_HUMAN.H11MO.0.A1.9045400000000002e-07
MA0677.1_Nr2f69.401180000000001e-05
MA0856.1_RXRG9.401180000000001e-05
MA0512.2_Rxra9.401180000000001e-05
MA1550.1_PPARD9.401180000000001e-05Not shown
MA1537.1_NR2F1(var.2)9.401180000000001e-05Not shown
MA1574.1_THRB9.401180000000001e-05Not shown
MA0855.1_RXRB0.00012073Not shown
MA0115.1_NR1H2::RXRA0.00012073Not shown

Motif 3/10

Motif IDq-valPWM
CEBPB_HUMAN.H11MO.0.A1.4822799999999999e-06
MA0837.1_CEBPE3.43904e-06
CEBPD_HUMAN.H11MO.0.C3.62526e-06
MA0466.2_CEBPB3.62526e-06
MA0838.1_CEBPG1.98859e-05
CEBPA_HUMAN.H11MO.0.A4.97149e-05Not shown
MA0836.2_CEBPD5.6816999999999995e-05Not shown
MA0043.3_HLF0.00034921800000000004Not shown
MA0102.4_CEBPA0.000494575Not shown
DBP_HUMAN.H11MO.0.B0.000581814Not shown

Motif 4/10

Motif IDq-valPWM
RXRG_HUMAN.H11MO.0.B6.68087e-05
MA0856.1_RXRG6.68087e-05
MA0512.2_Rxra6.68087e-05
MA0855.1_RXRB6.68087e-05
MA1550.1_PPARD6.68087e-05
MA0677.1_Nr2f66.68087e-05Not shown
MA1537.1_NR2F1(var.2)9.07814e-05Not shown
MA1574.1_THRB9.38121e-05Not shown
HNF4G_HUMAN.H11MO.0.B0.00010234Not shown
MA1148.1_PPARA::RXRA0.00010234Not shown

Motif 5/10

Motif IDq-valPWM
MA0489.1_JUN(var.2)1.10589e-06
MA1135.1_FOSB::JUNB5.27217e-06
MA1138.1_FOSL2::JUNB5.44951e-06
MA1144.1_FOSL2::JUND5.44951e-06
MA0099.3_FOS::JUN1.32948e-05
MA1134.1_FOS::JUNB1.32948e-05Not shown
FOSB_HUMAN.H11MO.0.A2.09275e-05Not shown
MA1622.1_Smad2::Smad32.09275e-05Not shown
MA0476.1_FOS2.09275e-05Not shown
MA0478.1_FOSL22.09275e-05Not shown

Motif 6/10

Motif IDq-valPWM
FOXA2_HUMAN.H11MO.0.A0.20843000000000003
FOXA3_HUMAN.H11MO.0.B0.20843000000000003
MA1683.1_FOXA30.20843000000000003
MA0848.1_FOXO40.20843000000000003
MA0042.2_FOXI10.20843000000000003
MA0849.1_FOXO60.20843000000000003Not shown
MA0795.1_SMAD30.20843000000000003Not shown
MA0846.1_FOXC20.20843000000000003Not shown
MA1602.1_ZSCAN290.20843000000000003Not shown
FOXA1_HUMAN.H11MO.0.A0.20843000000000003Not shown

Motif 7/10

Motif IDq-valPWM
HNF4G_HUMAN.H11MO.0.B4.5780599999999996e-08
HNF4A_HUMAN.H11MO.0.A2.0248899999999998e-07
MA1550.1_PPARD4.79134e-07
MA0856.1_RXRG4.79134e-07
MA0677.1_Nr2f64.84373e-07
MA0512.2_Rxra6.317159999999999e-07Not shown
MA0855.1_RXRB1.22984e-06Not shown
MA0504.1_NR2C23.60895e-06Not shown
MA1574.1_THRB4.7264e-06Not shown
MA1537.1_NR2F1(var.2)7.0257e-06Not shown

Motif 8/10

Motif IDq-valPWM
HNF4A_HUMAN.H11MO.0.A5.43947e-07
HNF4G_HUMAN.H11MO.0.B5.43947e-07
MA0512.2_Rxra4.37911e-06
MA0677.1_Nr2f64.37911e-06
MA0856.1_RXRG4.37911e-06
MA1550.1_PPARD5.42814e-06Not shown
MA0855.1_RXRB5.42814e-06Not shown
MA1537.1_NR2F1(var.2)3.9874e-05Not shown
MA1574.1_THRB3.9874e-05Not shown
MA0115.1_NR1H2::RXRA0.000141917Not shown

Motif 9/10

Motif IDq-valPWM
HNF4G_HUMAN.H11MO.0.B1.72302e-08
HNF4A_HUMAN.H11MO.0.A3.1357600000000003e-07
MA0677.1_Nr2f60.000178403
MA1537.1_NR2F1(var.2)0.000281475
MA0856.1_RXRG0.000281475
MA1550.1_PPARD0.000281475Not shown
MA0115.1_NR1H2::RXRA0.000281475Not shown
MA0512.2_Rxra0.000281475Not shown
MA0855.1_RXRB0.000422984Not shown
MA1574.1_THRB0.000422984Not shown

Motif 10/10

Motif IDq-valPWM
HNF1B_HUMAN.H11MO.0.A4.64844e-09
HNF1A_HUMAN.H11MO.0.C2.20822e-07
MA0153.2_HNF1B5.33434e-07
MA0046.2_HNF1A1.73137e-06
HNF1B_HUMAN.H11MO.1.A1.33257e-05
ZFHX3_HUMAN.H11MO.0.D0.0315857Not shown
MEOX2_HUMAN.H11MO.0.D0.13455899999999998Not shown
PAX4_HUMAN.H11MO.0.D0.155851Not shown
MA0853.1_Alx40.230223Not shown
MA0601.1_Arid3b0.241585Not shown

Metacluster 2/2

Motif 1/6

No TOMTOM matches passing threshold

Motif 2/6

No TOMTOM matches passing threshold

Motif 3/6

Motif IDq-valPWM
PLAG1_HUMAN.H11MO.0.D0.299646
MA1644.1_NFYC0.299646
FOXI1_HUMAN.H11MO.0.B0.299646
MA0060.3_NFYA0.299646
MA0503.1_Nkx2-5(var.2)0.463887

Motif 4/6

Motif IDq-valPWM
EOMES_HUMAN.H11MO.0.D0.100097
MA0850.1_FOXP30.100097
MA0042.2_FOXI10.100097
MA1567.1_TBX60.100097
MA0848.1_FOXO40.100097
MA0849.1_FOXO60.100097Not shown
MA1565.1_TBX180.100097Not shown
MA1566.1_TBX30.114508Not shown
TBX21_HUMAN.H11MO.0.A0.114508Not shown
MA0033.2_FOXL10.114508Not shown

Motif 5/6

No TOMTOM matches passing threshold

Motif 6/6

No TOMTOM matches passing threshold