Link to results

DiChIPMunk: on peaks and on multi-task seqlets and on single-task seqlets

HOMER: on peaks and on multi-task seqlets and on single-task seqlets

MEME: on peaks and on multitask seqlets and on single-task seqlets

import sys
import os
sys.path.append(os.path.abspath("/users/amtseng/tfmodisco/src/"))
from util import figure_to_vdom_image
import motif.read_motifs as read_motifs
from motif.read_motifs import pfm_to_pwm
import plot.viz_sequence as viz_sequence
import numpy as np
import matplotlib.pyplot as plt
import vdom.helpers as vdomh
from IPython.display import display

# Define parameters/fetch arguments
tf_name = os.environ["TFM_TF_NAME"]
multitask_fold = int(os.environ["TFM_MULTITASK_FOLD"])

if "TFM_TASK_INDEX" in os.environ:
    task_index = int(os.environ["TFM_TASK_INDEX"])
    singletask_fold = int(os.environ["TFM_SINGLETASK_FOLD"])
else:
    task_index = None
    singletask_fold = None
    
print("TF name: %s" % tf_name)
print("Multi-task fold: %s" % multitask_fold)
print("Task index: %s" % task_index)
print("Single-task fold: %s" % singletask_fold)

TF name: CEBPB
Multi-task fold: 7
Task index: 0
Single-task fold: 7

# Define paths and constants
base_path = "/users/amtseng/tfmodisco/results/classic_motifs/"

multitask_seqlets_dir = os.path.join(
    base_path, "seqlets", "multitask_profile_finetune",
    "%s_multitask_profile_finetune_fold%s" % (tf_name, multitask_fold)
)

if task_index is None:
    peaks_path = os.path.join(base_path, "peaks", tf_name, "%s_peaks_taskall" % tf_name)
    multitask_profile_seqlets_path = os.path.join(
        multitask_seqlets_dir,
        "%s_seqlets_profile_taskall" % tf_name
    )
    multitask_count_seqlets_path = os.path.join(
        multitask_seqlets_dir,
        "%s_seqlets_count_taskall" % tf_name
    )
else:
    peaks_path = os.path.join(base_path, "peaks", tf_name, "%s_peaks_task%d" % (tf_name, task_index))
    multitask_profile_seqlets_path = os.path.join(
        multitask_seqlets_dir,
        "%s_seqlets_profile_task%d" % (tf_name, task_index)
    )
    multitask_count_seqlets_path = os.path.join(
        multitask_seqlets_dir,
        "%s_seqlets_count_task%d" % (tf_name, task_index)
    )
    
    singletask_seqlets_dir = os.path.join(
        base_path, "seqlets", "singletask_profile_finetune",
        "%s_singletask_profile_finetune_fold%s" % (tf_name, singletask_fold),
        "task_%d" % task_index
    )
    singletask_profile_seqlets_path = os.path.join(
        singletask_seqlets_dir,
        "%s_seqlets_profile_task%d" % (tf_name, task_index)
    )
    singletask_count_seqlets_path = os.path.join(
        singletask_seqlets_dir,
        "%s_seqlets_count_task%d" % (tf_name, task_index)
    )

Helper functions

def show_peaks_motif_table(results_path, mode):
    """
    Shows a table of motifs from the given results path.
    `mode` is either `dichipmunk`, `homer`, `meme`, or `memechip`.
    """
    assert mode in ("dichipmunk", "homer", "meme", "memechip")
    if mode == "dichipmunk":
        score_name = "Supporting sequences"
        pfms, score_vals = read_motifs.import_dichipmunk_pfms(results_path)
    elif mode == "homer":
        score_name = "Log enrichment"
        pfms, score_vals = read_motifs.import_homer_pfms(results_path)
    elif mode == "meme":
        score_name = "E-value"
        pfms, score_vals = read_motifs.import_meme_pfms(results_path)
    else:
        score_name = "E-value"
        pfms, score_vals = read_motifs.import_meme_pfms(
            os.path.join(results_path, "meme_out")
        )
        
    colgroup = vdomh.colgroup(
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "40%"})
    )
    header = vdomh.thead(
        vdomh.tr(
            vdomh.th("Motif", style={"text-align": "center"}),
            vdomh.th(score_name, style={"text-align": "center"}),
            vdomh.th("PWM", style={"text-align": "center"})
        )
    )

    body = []
    for i, pfm in enumerate(pfms):
        pwm = pfm_to_pwm(pfm)
        if np.sum(pwm[:, [0, 2]]) < 0.5 * np.sum(pwm):
            # Flip to purine-rich version
            pwm = np.flip(pwm, axis=(0, 1))
        fig = viz_sequence.plot_weights(pwm, figsize=(20, 4), return_fig=True)
        fig.tight_layout()

        body.append(
            vdomh.tr(
                vdomh.td(str(i + 1)),
                vdomh.td(str(score_vals[i])),
                vdomh.td(figure_to_vdom_image(fig))
            )
        )

    display(vdomh.table(colgroup, header, vdomh.tbody(*body)))
    plt.close("all")

def show_seqlets_motif_table(profile_results_path, count_results_path, mode):
    """
    Shows a table of motifs from the given results path.
    `mode` is either `dichipmunk`, `homer`, `meme`, or `memechip`
    """
    assert mode in ("dichipmunk", "homer", "meme", "memechip")
    if mode == "dichipmunk":
        score_name = "Supporting sequences"
        p_pfms, p_score_vals = read_motifs.import_dichipmunk_pfms(profile_results_path)
        c_pfms, c_score_vals = read_motifs.import_dichipmunk_pfms(count_results_path)
    elif mode == "homer":
        score_name = "Log enrichment"
        p_pfms, p_score_vals = read_motifs.import_homer_pfms(profile_results_path)
        c_pfms, c_score_vals = read_motifs.import_homer_pfms(count_results_path)
    elif mode == "meme":
        score_name = "E-value"
        p_pfms, p_score_vals = read_motifs.import_meme_pfms(profile_results_path)
        c_pfms, c_score_vals = read_motifs.import_meme_pfms(count_results_path)
    else:
        score_name = "E-value"
        p_pfms, p_score_vals = read_motifs.import_meme_pfms(
            os.path.join(profile_results_path, "meme_out")
        )
        c_pfms, c_score_vals = read_motifs.import_meme_pfms(
            os.path.join(count_results_path, "meme_out")
        )
        
    colgroup = vdomh.colgroup(
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "40%"}),
        vdomh.col(style={"width": "5%"}),
        vdomh.col(style={"width": "40%"})
    )
    header = vdomh.thead(
        vdomh.tr(
            vdomh.th("Motif", style={"text-align": "center"}),
            vdomh.th(score_name + " (profile)", style={"text-align": "center"}),
            vdomh.th("PWM (profile)", style={"text-align": "center"}),
            vdomh.th(score_name + " (count)", style={"text-align": "center"}),
            vdomh.th("PWM (count)", style={"text-align": "center"})
        )
    )

    body = []
    for i in range(max(len(p_pfms), len(c_pfms))):
        rows = [vdomh.td(str(i + 1))]
        if i < len(p_pfms):
            pwm = pfm_to_pwm(p_pfms[i])
            if np.sum(pwm[:, [0, 2]]) < 0.5 * np.sum(pwm):
                # Flip to purine-rich version
                pwm = np.flip(pwm, axis=(0, 1))
            fig = viz_sequence.plot_weights(pwm, figsize=(20, 4), return_fig=True)
            fig.tight_layout()
            rows.extend([
                vdomh.td(str(p_score_vals[i])),
                vdomh.td(figure_to_vdom_image(fig))
            ])
        else:
            rows.extend([vdomh.td(), vdomh.td()])
            
        if i < len(c_pfms):
            pwm = pfm_to_pwm(c_pfms[i])
            if np.sum(pwm[:, [0, 2]]) < 0.5 * np.sum(pwm):
                # Flip to purine-rich version
                pwm = np.flip(pwm, axis=(0, 1))
            fig = viz_sequence.plot_weights(pwm, figsize=(20, 4), return_fig=True)
            fig.tight_layout()
            rows.extend([
                vdomh.td(str(c_score_vals[i])),
                vdomh.td(figure_to_vdom_image(fig))
            ])
        else:
            rows.extend([vdomh.td(), vdomh.td()])
            

        body.append(vdomh.tr(*rows))

    display(vdomh.table(colgroup, header, vdomh.tbody(*body)))
    plt.close("all")

Show benchmark motifs

DiChIPMunk on peaks

show_peaks_motif_table(os.path.join(peaks_path, "dichipmunk"), "dichipmunk")

Motif	Supporting sequences	PWM
1	1800
2	2000
3	2000
4	2000
5	1628
6	913
7	335
8	86
9	25
10	18

DiChIPMunk on multi-task seqlets

show_seqlets_motif_table(
    os.path.join(multitask_profile_seqlets_path, "dichipmunk"),
    os.path.join(multitask_count_seqlets_path, "dichipmunk"),
    "dichipmunk"
)

Motif	Supporting sequences (profile)	PWM (profile)	Supporting sequences (count)	PWM (count)
1	13810		15317
2	2079		1322
3	659		175
4	413		120
5			78
6			44

DiChIPMunk on single-task seqlets

if task_index is not None:
    show_seqlets_motif_table(
        os.path.join(singletask_profile_seqlets_path, "dichipmunk"),
        os.path.join(singletask_count_seqlets_path, "dichipmunk"),
        "dichipmunk"
    )

Motif	Supporting sequences (profile)	PWM (profile)	Supporting sequences (count)	PWM (count)
1	14846		16036
2	2017		1462
3	194		116
4	60		7
5	1

HOMER on peaks

show_peaks_motif_table(os.path.join(peaks_path, "homer"), "homer")

Motif	Log enrichment	PWM
1	-23381.782555
2	-6440.771162
3	-3708.969408
4	-3588.026123
5	-3510.919406
6	-3429.225319
7	-3416.102803
8	-2948.841
9	-2829.056889
10	-2386.236895
11	-2265.972397
12	-1428.377979
13	-1363.737305
14	-1306.367475
15	-1077.45807
16	-1065.271271
17	-1032.047176
18	-263.213401
19	-250.175692

HOMER on multi-task seqlets

show_seqlets_motif_table(
    os.path.join(multitask_profile_seqlets_path, "homer"),
    os.path.join(multitask_count_seqlets_path, "homer"),
    "homer"
)

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

Motif	Log enrichment (profile)	PWM (profile)	Log enrichment (count)	PWM (count)
1	-19855.090039		-30308.553771
2	-3844.158413		-3770.46649
3	-3248.946282		-3463.110153
4	-2513.719948		-3146.785216
5	-2286.655204		-2057.0444
6	-2078.766947		-1943.708211
7	-1430.88208		-1937.171771
8	-1195.335244		-1117.527814
9	-923.838882		-954.364358
10	-744.774094		-688.308241
11	-705.82604		-630.544162
12	-694.621665		-274.660338
13	-663.488113		-253.668958
14	-601.862331		-218.411604
15	-499.378094
16	-440.609689
17	-422.461764
18	-347.915285
19	-190.115833
20	-103.957809
21	-87.414862

HOMER on single-task seqlets

if task_index is not None:
    show_seqlets_motif_table(
        os.path.join(singletask_profile_seqlets_path, "homer"),
        os.path.join(singletask_count_seqlets_path, "homer"),
        "homer"
    )

Motif	Log enrichment (profile)	PWM (profile)	Log enrichment (count)	PWM (count)
1	-29538.28702		-34742.471769
2	-3550.748299		-3575.800943
3	-2307.374426		-3525.369837
4	-2126.663084		-2724.118078
5	-1971.73812		-2575.222865
6	-1634.300594		-2315.901311
7	-1417.379139		-382.733068
8	-618.322023		-232.623177
9	-606.766173		-159.816938
10	-445.740357		-150.988834
11	-215.41867		-108.951224
12	-189.390414		-90.560808
13	-179.130107		-67.541734
14	-169.360276		-52.08057
15	-73.623228		-38.054462
16	-19.95657

MEMEChIP on peaks

show_peaks_motif_table(os.path.join(peaks_path, "memechip"), "memechip")

Motif	E-value	PWM
1	0.0
2	5.5e-294
3	2.7e-222
4	5.7e-75
5	5e-53
6	2.9e-39
7	4.2e-24
8	1.8e-21
9	2.3e-12
10	1.2e-11

MEME on multi-task seqlets

show_seqlets_motif_table(
    os.path.join(multitask_profile_seqlets_path, "meme"),
    os.path.join(multitask_count_seqlets_path, "meme"),
    "meme"
)

Motif	E-value (profile)	PWM (profile)	E-value (count)	PWM (count)
1	0.0		0.0
2	2.3e-93		8.2e-107
3	4.2e-42		7.6e-34
4	2.2e-27		1.9e-13
5	5.9e-26		4.1e-11
6	5.7e-16		1.5
7	5.8e-08		0.028
8	4.7e-05		250.0
9	0.0013		730.0
10	0.011		78.0

MEME on single-task seqlets

if task_index is not None:
    show_seqlets_motif_table(
        os.path.join(singletask_profile_seqlets_path, "meme"),
        os.path.join(singletask_count_seqlets_path, "meme"),
        "meme"
    )

Motif	E-value (profile)	PWM (profile)	E-value (count)	PWM (count)
1	0.0		0.0
2	1.2e-83		4.6e-68
3	1.6e-20		1.2e-34
4	2.7e-19		7.6e-27
5	5.8e-08		1.3e-20
6	0.28		1.9e-06
7	110.0		13.0
8	2200.0		43.0
9	2900.0		210.0
10	32000.0		270.0