%load_ext autoreload
%autoreload 2

import pandas as pd
import numpy as np 
import glob
import os
from collections import OrderedDict
import pickle
import h5py

Generating modisco motifs¶

def load_deeplift_data(deeplift_hdf, keys=['scores', 'one_hot']):
    fp = h5py.File(deeplift_hdf, "r")
    hyp_scores = fp['deeplift_scores'][:]
    one_hot = fp['inputs'][:]
    shuffled_onehot = fp['shuffled_inputs'][:]

    deeplift_data = OrderedDict()
    if 'one_hot' in keys:
        deeplift_data['one_hot'] = one_hot
    if 'peaks' in keys :
        df = OrderedDict()
        for key in list(fp['metadata/range']):
            df[key] = fp['metadata/range/{}'.format(key)][:]
        df = pd.DataFrame(df)
        df.chr = np.array([df.chr.values[i].decode('utf-8') for i in range(df.shape[0])])
        deeplift_data['peaks'] = df
        
    if 'hyp_scores' in keys :
        deeplift_data['hyp_scores'] = hyp_scores
    if 'scores' in keys :
        deeplift_data['scores'] = np.multiply(hyp_scores, one_hot)#no need to sum before mult, taken care of in deeplift, not in deepshap though
    
    
    return deeplift_data

deeplift_hdf = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/deeplift_out/summit.h5"
deeplift_data = load_deeplift_data(deeplift_hdf, keys=['hyp_scores', 'scores', 'one_hot'])
deeplift_data['scores'].shape, deeplift_data['hyp_scores'].shape, deeplift_data['one_hot'].shape

((35024, 1000, 4), (35024, 1000, 4), (35024, 1000, 4))

import modisco
null_per_pos_scores = modisco.coordproducers.LaplaceNullDist(num_to_samp=5000)

ggrfile = "/mnt/lab_data3/dskim89/ggr/nn/2019-03-12.freeze/motifs.input_x_grad.early/ggr.scanmotifs.h5"
with h5py.File(ggrfile, "r") as fp:
    traj_lab = fp['TRAJ_LABELS'][:]
traj_lab.shape
# #0, 7, 8-10-11, and 9
# indices = np.logical_or(np.logical_or(traj_lab[:, 8]==1,  traj_lab[:, 10]==1), traj_lab[:, 11]==1)
# indices.shape, deeplift_data['scores'][indices,:,:].shape, np.sum(indices)

(35024, 15)

traj_no = -1
if traj_no==8:
    indices = np.logical_or(np.logical_or(traj_lab[:, 8]==1,  traj_lab[:, 10]==1), traj_lab[:, 11]==1)
else:
    indices = traj_lab[:, traj_no]==1

if traj_no<0:
    _score = deeplift_data['scores']
    _hyp_score = deeplift_data['hyp_scores']
    _one_hot = deeplift_data['one_hot']
else:
    _score = deeplift_data['scores'][indices,:,:]
    _hyp_score = deeplift_data['hyp_scores'][indices,:,:]
    _one_hot = deeplift_data['one_hot'][indices,:,:]

tfmodisco_results = modisco.tfmodisco_workflow.workflow.TfModiscoWorkflow(
                    #Slight modifications from the default settings
                    sliding_window_size=[5,9,13,17],
                    flank_size=5,
                    target_seqlet_fdr=0.15,
                    
                    min_passing_windows_frac=0.03,
                    max_seqlets_per_metacluster=20000,
                    
                    seqlets_to_patterns_factory=
                     modisco.tfmodisco_workflow.seqlets_to_patterns.TfModiscoSeqletsToPatternsFactory(
                        #Note: as of version 0.5.6.0, it's possible to use the results of a motif discovery
                        # software like MEME to improve the TF-MoDISco clustering. To use the meme-based
                        # initialization, you would specify the initclusterer_factory as shown in the
                        # commented-out code below:
                        #initclusterer_factory=modisco.clusterinit.memeinit.MemeInitClustererFactory(    
                        #    meme_command="meme", base_outdir="meme_out",            
                        #    max_num_seqlets_to_use=10000, nmotifs=10, n_jobs=1),
                        
                        embedder_factory=(modisco.seqlet_embedding
                          .advanced_gapped_kmer.AdvancedGappedKmerEmbedderFactory(max_entries=500)),
                        
                        trim_to_window_size=15,
                        initial_flank_to_add=5,
                        #kmer_len=5, num_gaps=1,
                        #num_mismatches=0,
                        final_min_cluster_size=60)
                )(
                 task_names=["early0"],# "task1", "task2"],
                 contrib_scores={'early0': _score}, #[:,420:580,:]
                 hypothetical_contribs={'early0': _hyp_score }, #[:,420:580,:]
                 one_hot=_one_hot, #[:,420:580,:]
                 null_per_pos_scores = null_per_pos_scores)

MEMORY 7.065616384
On task early0
Fitting - on window size 5
peak(mu)= -1.6098946790181717e-07
Computing window sums
Done computing window sums
Subsampling!
For increasing = True , the minimum IR precision was 0.4475815087661602 occurring at 0.0 implying a frac_neg of 0.8102218080471498
To be conservative, adjusted frac neg is 0.95
For increasing = False , the minimum IR precision was 0.37161335387877126 occurring at -1.149310263626191e-08 implying a frac_neg of 0.5913769112895493
To be conservative, adjusted frac neg is 0.95

Fitting - on window size 9
peak(mu)= -9.305308456869672e-07
Computing window sums
Done computing window sums
Subsampling!
For increasing = True , the minimum IR precision was 0.44210022527471415 occurring at 0.0 implying a frac_neg of 0.7924366441847153
To be conservative, adjusted frac neg is 0.95
For increasing = False , the minimum IR precision was 0.4003950993610411 occurring at -4.536996130177773e-08 implying a frac_neg of 0.667764888069405
To be conservative, adjusted frac neg is 0.95

Fitting - on window size 13
peak(mu)= 6.035615486338266e-07
Computing window sums
Done computing window sums
Subsampling!
For increasing = True , the minimum IR precision was 0.4397550756148556 occurring at 0.0 implying a frac_neg of 0.7849336182696817
To be conservative, adjusted frac neg is 0.95
For increasing = False , the minimum IR precision was 0.3506084713838218 occurring at -1.0496441937846157e-07 implying a frac_neg of 0.5399030568368383
To be conservative, adjusted frac neg is 0.95

Fitting - on window size 17
peak(mu)= 2.842042537608714e-06
Computing window sums
Done computing window sums
Subsampling!
For increasing = True , the minimum IR precision was 0.44183136128423395 occurring at 0.0 implying a frac_neg of 0.7915732462160525
To be conservative, adjusted frac neg is 0.95
For increasing = False , the minimum IR precision was 0.2544225791333979 occurring at -2.1372761692395326e-08 implying a frac_neg of 0.34124233381112395
To be conservative, adjusted frac neg is 0.95

Subsampling!

Thresholds from null dist were -0.85  and  0.85 with frac passing 0.100903
Got 332134 coords
After resolving overlaps, got 332134 seqlets
Across all tasks, the weakest transformed threshold used was: 0.8499
MEMORY 9.551650816
332134 identified in total
min_metacluster_size_frac * len(seqlets) = 3321 is more than min_metacluster_size=100.
Using it as a new min_metacluster_size
2 activity patterns with support >= 3321 out of 2 possible patterns
Metacluster sizes:  [290277, 41857]
Idx to activities:  {0: '1', 1: '-1'}
MEMORY 9.551650816
On metacluster 1
Metacluster size 41857 limited to 20000
Relevant tasks:  ('early0',)
Relevant signs:  (-1,)
TfModiscoSeqletsToPatternsFactory: seed=1234
(Round 1) num seqlets: 20000
(Round 1) Computing coarse affmat
MEMORY 9.551650816
Beginning embedding computation

[Parallel(n_jobs=10)]: Using backend LokyBackend with 10 concurrent workers.
[Parallel(n_jobs=10)]: Done  30 tasks      | elapsed:    4.7s
[Parallel(n_jobs=10)]: Done 180 tasks      | elapsed:    9.5s
[Parallel(n_jobs=10)]: Done 430 tasks      | elapsed:   18.2s
[Parallel(n_jobs=10)]: Done 780 tasks      | elapsed:   29.6s
[Parallel(n_jobs=10)]: Done 1230 tasks      | elapsed:   42.2s
[Parallel(n_jobs=10)]: Done 1780 tasks      | elapsed:   57.3s
[Parallel(n_jobs=10)]: Done 2430 tasks      | elapsed:  1.2min
[Parallel(n_jobs=10)]: Done 3180 tasks      | elapsed:  1.6min
[Parallel(n_jobs=10)]: Done 4030 tasks      | elapsed:  1.9min
[Parallel(n_jobs=10)]: Done 4980 tasks      | elapsed:  2.2min
[Parallel(n_jobs=10)]: Done 6030 tasks      | elapsed:  2.6min
[Parallel(n_jobs=10)]: Done 7180 tasks      | elapsed:  3.1min
[Parallel(n_jobs=10)]: Done 8430 tasks      | elapsed:  3.5min
[Parallel(n_jobs=10)]: Done 9780 tasks      | elapsed:  4.0min
[Parallel(n_jobs=10)]: Done 11230 tasks      | elapsed:  4.6min
[Parallel(n_jobs=10)]: Done 12780 tasks      | elapsed:  5.2min
[Parallel(n_jobs=10)]: Done 14430 tasks      | elapsed:  5.8min
[Parallel(n_jobs=10)]: Done 16180 tasks      | elapsed:  6.5min
[Parallel(n_jobs=10)]: Done 18030 tasks      | elapsed:  7.1min
[Parallel(n_jobs=10)]: Done 19980 tasks      | elapsed:  7.9min
[Parallel(n_jobs=10)]: Done 20000 out of 20000 | elapsed:  7.9min finished
[Parallel(n_jobs=10)]: Using backend LokyBackend with 10 concurrent workers.
[Parallel(n_jobs=10)]: Done  40 tasks      | elapsed:    1.6s
[Parallel(n_jobs=10)]: Done 340 tasks      | elapsed:   11.7s
[Parallel(n_jobs=10)]: Done 840 tasks      | elapsed:   27.1s
[Parallel(n_jobs=10)]: Done 1540 tasks      | elapsed:   47.1s
[Parallel(n_jobs=10)]: Done 2440 tasks      | elapsed:  1.2min
[Parallel(n_jobs=10)]: Done 3540 tasks      | elapsed:  1.6min
[Parallel(n_jobs=10)]: Done 4840 tasks      | elapsed:  2.2min
[Parallel(n_jobs=10)]: Done 6340 tasks      | elapsed:  2.7min
[Parallel(n_jobs=10)]: Done 8040 tasks      | elapsed:  3.4min

modisco_hdf = '/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_out_1k/results.hdf5'
grp = h5py.File(modisco_hdf)
tfmodisco_results.save_hdf5(grp)

print('..')

..

from matlas.matches import DenovoModisco, DenovoHomer
from vdom.helpers import (b, summary, details)
from IPython.display import display
import numpy as np


def show_patterns_using_hoccomocco_db(sample_name, modiscodir):
    ob = DenovoModisco(modiscodir)
#     ob.fetch_tomtom_matches(
#                 meme_db="/mnt/lab_data/kundaje/users/msharmin/annotations/HOCOMOCOv11_core_pwms_HUMAN_mono.renamed.nonredundant.annotated.meme",
#                 database_name="HOCOMOCO.nonredundant.annotated",
#                 save_report=True, tomtom_dir= "{0}/{1}_tomtomout".format(modiscodir, "HOCOMOCO.nonredundant.annotated"))
    ob.load_matched_motifs(database_name="HOCOMOCO.nonredundant.annotated")
    ob.get_motif_per_celltype(match_threshold=0.03, database_name="HOCOMOCO.nonredundant.annotated")
    #ob.display_individual_table()
    
    pattern_tab, pattern_dict = ob.visualize_pattern_table()
    display(details(summary('Click here for ', b('Denovo Patterns'), ' by ', b('{}'.format('MoDISco')),
                        ' in ', b(sample_name),
                        ": #{}".format(len(pattern_dict)),
                       ), pattern_tab))
    return None



def display_denovo_patterns(sample_name, modiscodir, match_threshold=0.05):
    display(summary(b(sample_name)))
    
    ob = DenovoModisco(modiscodir)
#     ob.fetch_tomtom_matches(save_report=True, 
#                                   tomtom_dir= "{0}/{1}_tomtomout".format(modiscodir, "CISBP_2.00"))
    
    ob.load_matched_motifs()
    ob.get_motif_per_celltype(match_threshold=match_threshold)
    pattern_tab, pattern_dict = ob.visualize_pattern_table()
    display(details(summary('Click here for ', b('Denovo Patterns'), ' by ', b('{}'.format('MoDISco')),
                        ' in ', b(sample_name),
                        ": #{}".format(len(pattern_dict)),
                       ), pattern_tab))
    #ob.display_individual_table()
    
    return None

Early timepoint¶

With all peaks¶

using central 160bp¶

# sample_name = 'early_fold0'

# display_denovo_patterns(
#     sample_name,
#     modiscodir="/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_out_1k"
# )

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_out"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

using central 1kb¶

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_out_1k"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

With time-point based peaks¶

using central 160bp¶

# sample_name = 'early_fold0'
# modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_out"

# show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_0"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_7"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_8"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_9"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

using central 1kbp¶

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_0_1k"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_7_1k"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_8_1k"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

sample_name = 'early_fold0'
modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/modisco_9_1k"

show_patterns_using_hoccomocco_db(sample_name, modiscodir)

Late timepoint¶

# sample_name = 'late_fold0'

# display_denovo_patterns(
#     sample_name,
#     modiscodir="/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_late/modisco_out"
# )

# sample_name = 'late_fold0'
# modiscodir = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_late/modisco_out"

# show_patterns_using_hoccomocco_db(sample_name, modiscodir)

Loading keras model¶

from matlas.model_test import getSkinModel
from matlas.model_test import setup_keras_session
setup_keras_session('4')
init_weights = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/weights_from_raw_tf.p"
model = getSkinModel(init_weights, 19, classification=False)
model_h5 = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/model.h5"
model.save(model_h5)
model.summary()

channels_last
compiling!
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input (InputLayer)           (None, 1000, 4)           0         
_________________________________________________________________
conv1d_7 (Conv1D)            (None, 1000, 300)         23100     
_________________________________________________________________
batch_normalization_11 (Batc (None, 1000, 300)         1200      
_________________________________________________________________
activation_14 (Activation)   (None, 1000, 300)         0         
_________________________________________________________________
max_pooling1d_7 (MaxPooling1 (None, 333, 300)          0         
_________________________________________________________________
conv1d_8 (Conv1D)            (None, 333, 200)          660200    
_________________________________________________________________
batch_normalization_12 (Batc (None, 333, 200)          800       
_________________________________________________________________
activation_15 (Activation)   (None, 333, 200)          0         
_________________________________________________________________
max_pooling1d_8 (MaxPooling1 (None, 83, 200)           0         
_________________________________________________________________
conv1d_9 (Conv1D)            (None, 83, 200)           280200    
_________________________________________________________________
batch_normalization_13 (Batc (None, 83, 200)           800       
_________________________________________________________________
activation_16 (Activation)   (None, 83, 200)           0         
_________________________________________________________________
max_pooling1d_9 (MaxPooling1 (None, 20, 200)           0         
_________________________________________________________________
flatten_3 (Flatten)          (None, 4000)              0         
_________________________________________________________________
dense_5 (Dense)              (None, 1000)              4001000   
_________________________________________________________________
batch_normalization_14 (Batc (None, 1000)              4000      
_________________________________________________________________
activation_17 (Activation)   (None, 1000)              0         
_________________________________________________________________
dropout_5 (Dropout)          (None, 1000)              0         
_________________________________________________________________
dense_6 (Dense)              (None, 1000)              1001000   
_________________________________________________________________
batch_normalization_15 (Batc (None, 1000)              4000      
_________________________________________________________________
activation_18 (Activation)   (None, 1000)              0         
_________________________________________________________________
dropout_6 (Dropout)          (None, 1000)              0         
_________________________________________________________________
final_dense19 (Dense)        (None, 19)                19019     
=================================================================
Total params: 5,995,319
Trainable params: 5,989,919
Non-trainable params: 5,400
_________________________________________________________________

ggrfile = "/mnt/lab_data3/dskim89/ggr/nn/2019-03-12.freeze/motifs.input_x_grad.early/ggr.scanmotifs.h5"
with h5py.File(ggrfile, "r") as fp:
    labels = fp['labels'][:]
    logits = fp['logits'][:]
    seqs = fp['sequence'][:]
labels.shape, logits.shape, seqs.shape

((35024, 19), (35024, 19), (35024, 1, 1000, 4))

from matlas.generators import EmbeddingsGenerator
def get_predictions(cur_seqs, model):
    e_generator = EmbeddingsGenerator(cur_seqs, batch_size=1000, num_rows=cur_seqs.shape[0])
    #batch = e_generator.get_batch(i)
    #e = model.predict_on_batch(batch[0])
    e = model.predict_generator(
                e_generator,
                max_queue_size=100,
                workers=1,
                use_multiprocessing=False,
                verbose=1
            )
    return e
keras_op = get_predictions(np.squeeze(seqs[:1000]), model)

1/1 [==============================] - 1s 664ms/step

from matplotlib import pylab as plt
# plt.scatter(activations_all['activation_2/Relu:0'][:1000,0,0,0], cnv1[:1000,0,0,0])
plt.scatter(logits[:1000,0], keras_op[:1000,0])
plt.xlabel('raw tensorflow prediction')
plt.ylabel('keras predictions')

Text(0, 0.5, 'keras predictions')

import scipy.stats
print(scipy.stats.pearsonr(logits[:1000,0], keras_op[:1000,0]))
print(scipy.stats.spearmanr(logits[:1000,0], keras_op[:1000,0]))

(0.7495659591048396, 4.981212730424334e-181)
SpearmanrResult(correlation=0.7251773091773093, pvalue=6.437351695707496e-164)

Deeplifting the keras model¶

from matlas.deeplift_run import *
contrib_funcs, input_layer_shape = retrieve_func_from_model(
    model_h5, 
    algorithm="rescale_conv_revealcancel_fc", 
    regression=True,
    sequential=False, 
    w0=None, w1=None, logger=None)
input_layer_shape

load data from labcluster

TF-MoDISco is using the TensorFlow backend.

nonlinear_mxts_mode is set to: DeepLIFT_GenomicsDefault
For layer activation_4_0 the preceding linear layer is conv1d_8_0 of type Conv1D;
In accordance with nonlinear_mxts_mode=DeepLIFT_GenomicsDefault we are setting the NonlinearMxtsMode to Rescale
Heads-up: current implementation assumes maxpool layer is followed by a linear transformation (conv/dense layer)
For layer activation_5_0 the preceding linear layer is conv1d_9_0 of type Conv1D;
In accordance with nonlinear_mxts_mode=DeepLIFT_GenomicsDefault we are setting the NonlinearMxtsMode to Rescale
Heads-up: current implementation assumes maxpool layer is followed by a linear transformation (conv/dense layer)
For layer activation_6_0 the preceding linear layer is conv1d_10_0 of type Conv1D;
In accordance with nonlinear_mxts_mode=DeepLIFT_GenomicsDefault we are setting the NonlinearMxtsMode to Rescale
Heads-up: current implementation assumes maxpool layer is followed by a linear transformation (conv/dense layer)
For layer activation_7_0 the preceding linear layer is dense_1_0 of type Dense;
In accordance with nonlinear_mxts_modeDeepLIFT_GenomicsDefault we are setting the NonlinearMxtsMode to RevealCancel
For layer activation_8_0 the preceding linear layer is dense_2_0 of type Dense;
In accordance with nonlinear_mxts_modeDeepLIFT_GenomicsDefault we are setting the NonlinearMxtsMode to RevealCancel

[None, 1000, 4]

#provide list of strings to run deeplift
# def read_ggr_active_sequences(ggr_h5):
#     with h5py.File(ggr_h5, "r") as fp:
#         seqs = fp['sequence.active.string'][:]
#     sequences = []
#     for seq in seqs:
#         sequences.append(seq[0].decode('utf-8'))
    
#     return sequences

# ggrfile = "/mnt/lab_data3/dskim89/ggr/nn/2019-03-12.freeze/motifs.input_x_grad.early/ggr.scanmotifs.h5"
# sequences = read_ggr_sequences(ggrfile)
# type(sequences), len(sequences) #sequences[0], sequences[1], seqs[0,0], seqs[1,0]

(list, 35024)

def get_genome_coordinates(ggr_h5, bed_file):
    with h5py.File(ggr_h5, "r") as fp:
        regions = fp['example_metadata'][:]
    
    chroms = []
    starts = []
    ends = []
    for region in regions[:,0]:
        region = region.decode("utf-8")
        if region!='':
            region = region.split("features=")[1]
        else:
            continue
        chroms.append(region.split(":")[0])
        starts.append(region.split(":")[1].split("-")[0])
        ends.append(region.split(":")[1].split("-")[1])
    df = pd.DataFrame({'chrom': chroms, 'start':starts, 'end': ends})
    df.to_csv(bed_file, header=False, index=False, sep="\t", compression="gzip")
    return None
ggrfile = "/mnt/lab_data3/dskim89/ggr/nn/2019-03-12.freeze/motifs.input_x_grad.early/ggr.scanmotifs.h5"
bed_file = "/mnt/lab_data2/msharmin/oc-atlas/DanSkinData/fold_0_ggr/result_early/regions.bed.gz"
get_genome_coordinates(ggrfile, bed_file)

from matlas.model_layer import retrieve_sequences
sequences, intervals_wo_flanks = retrieve_sequences(
    bed_file, 
    fasta_file="/mnt/lab_data3/dskim89/ggr/annotations/hg19.genome.fa", flank_size=0)

num_refs_per_seq = 10
from deeplift.dinuc_shuffle import dinuc_shuffle
from matlas.model_layer import one_hot_encode_along_col_axis
from matlas.dlutils import get_shuffled_seqs
input_data_list, input_references_list = get_shuffled_seqs(sequences[:45], num_refs_per_seq, shuffle_func=dinuc_shuffle,
                                                                one_hot_func=lambda x: np.array([one_hot_encode_along_col_axis(seq) for seq in x]),
                                                                progress_update=10000)
input_data_list[0].shape, len(sequences[0])
# input_data_list = [np.expand_dims(input_data_list[0], axis=1)]
# input_references_list = [np.expand_dims(input_references_list[0], axis=1)]

One hot encoding sequences...
One hot encoding done...

((450, 1000, 4), 1000)

from matlas.dlutils import get_given_seq_ref_function
shuffled_score_funcs = {input_name: get_given_seq_ref_function(score_computation_function=score_func)
                        for input_name, score_func in contrib_funcs.items()}

task_idx = 0
batch_size = 256
num_refs_per_seq = 10
for input_name, score_func in shuffled_score_funcs.items():
    hyp_scores = None
    b = 10000
    c = int(np.ceil(1.0*len(input_data_list[0])/b))
    for si in range(c):
        if(si==c-1):
            tmp = score_func(task_idx=int(task_idx), input_data_list=[input_data_list[0][si*b:len(input_data_list[0])]],
                               input_references_list=[input_references_list[0][si*b:len(input_data_list[0])]],
                               num_refs_per_seq=num_refs_per_seq, batch_size=batch_size,
                               progress_update=10000)
        else:
            #print('batch: ', si, si*b, (si+1)*b) 
            tmp = score_func(task_idx=int(task_idx), input_data_list=[input_data_list[0][si*b:(si+1)*b]],
                               input_references_list=[input_references_list[0][si*b:(si+1)*b]],
                               num_refs_per_seq=num_refs_per_seq, 
                               batch_size=batch_size,
                               progress_update=10000)
        if(hyp_scores is None):
            hyp_scores = tmp
        else:
            hyp_scores = np.vstack((hyp_scores, tmp))
    input_data_list[0] = np.squeeze(input_data_list[0])
    input_references_list[0] = np.squeeze(input_references_list[0])
    one_hot = input_data_list[0][[range(0, len(input_data_list[0]), num_refs_per_seq)]]
    shuffled_onehot = input_references_list[0].reshape((one_hot.shape[0], num_refs_per_seq, 
                                                       input_references_list[0].shape[-2], #seq_len
                                                        input_references_list[0].shape[-1]))#alphabet 
    scores = np.multiply(hyp_scores, one_hot)
       
hyp_scores.shape, one_hot.shape, scores.shape

Done 0

/users/msharmin/anaconda2/envs/basepair13/lib/python3.6/site-packages/ipykernel_launcher.py:27: FutureWarning: Using a non-tuple sequence for multidimensional indexing is deprecated; use `arr[tuple(seq)]` instead of `arr[seq]`. In the future this will be interpreted as an array index, `arr[np.array(seq)]`, which will result either in an error or a different result.

# create_deeplift_h5(bed_file, score_hdf, hyp_scores, one_hot, shuffled_onehot)

ggrfile = "/mnt/lab_data3/dskim89/ggr/nn/2019-03-12.freeze/motifs.input_x_grad.early/ggr.scanmotifs.h5"
with h5py.File(ggrfile, "r") as fp:
    print(list(fp))
    scores_ggr = fp['sequence-weighted'][:]
    scores_ggr_active = fp['sequence-weighted.active'][:]
scores_ggr.shape, scores_ggr_active.shape

['ATAC_LABELS', 'ATAC_SIGNALS', 'ATAC_SIGNALS.NORM', 'CTCF_LABELS', 'CTCF_SIGNALS', 'CTCF_SIGNALS.NORM', 'DYNAMIC_MARK_LABELS', 'DYNAMIC_STATE_LABELS', 'H3K27ac_LABELS', 'H3K27ac_SIGNALS', 'H3K27ac_SIGNALS.NORM', 'H3K27me3_LABELS', 'H3K27me3_SIGNALS', 'H3K27me3_SIGNALS.NORM', 'H3K4me1_LABELS', 'H3K4me1_SIGNALS', 'H3K4me1_SIGNALS.NORM', 'KLF4_LABELS', 'POL2_LABELS', 'STABLE_MARK_LABELS', 'STABLE_STATE_LABELS', 'TP63_LABELS', 'TRAJ_LABELS', 'ZNF750_LABELS', 'example_metadata', 'gradients', 'labels', 'logits', 'logits.ci', 'logits.ci.thresh', 'logits.multimodel', 'logits.multimodel.norm', 'logits.norm', 'positive_importance_bp_sum', 'probs', 'pwm-scores.null.idx', 'sequence', 'sequence-weighted', 'sequence-weighted.active', 'sequence-weighted.active.ci', 'sequence-weighted.active.ci.thresh', 'sequence-weighted.active.pwm-scores.thresh', 'sequence-weighted.active.pwm-scores.thresh.max.idx', 'sequence-weighted.active.pwm-scores.thresh.max.val', 'sequence-weighted.active.pwm-scores.thresh.sum', 'sequence-weighted.thresholds', 'sequence.active', 'sequence.active.gc_fract', 'sequence.active.pwm-hits', 'sequence.active.pwm-hits.densities', 'sequence.active.pwm-hits.densities.max', 'sequence.active.pwm-scores.thresh', 'sequence.active.pwm-scores.thresh.sum', 'sequence.active.string']

((35024, 10, 1000, 4), (35024, 10, 160, 4))

import modisco.visualization
from modisco.visualization import viz_sequence
viz_sequence.plot_weights(scores_ggr[0,0][500:600], subticks_frequency=20)
viz_sequence.plot_weights(scores_ggr_active[0,0], subticks_frequency=20)

-0.4860307276248932 1.3380632400512695

-0.0440836176276207 0.1500825583934784

scores_ggr.shape

(35024, 10, 1000, 4)

import modisco.visualization
from modisco.visualization import viz_sequence

viz_sequence.plot_weights(scores[0][500:600], subticks_frequency=20)
viz_sequence.plot_weights(hyp_scores[0][500:600], subticks_frequency=20)
viz_sequence.plot_weights(one_hot[0][500:600], subticks_frequency=20)
viz_sequence.plot_weights(scores_ggr[0][500:600], subticks_frequency=20)

-0.0688343504909426 0.35437235310673715

-0.976342553505674 0.35437235310673715

0.0 1.0

Pattern Name	TF Name(s)	Modisco
metacluster_0/pattern_0 # seqlets: 10140 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-101_NFE2.UNK.0.A, HCLUST-179_BACH1.UNK.0.A
metacluster_0/pattern_1 # seqlets: 643 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-170_TP53.UNK.0.A
metacluster_0/pattern_2 # seqlets: 598 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_3 # seqlets: 536 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A
metacluster_0/pattern_4 # seqlets: 304 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-176_CBFB.UNK.0.A
metacluster_0/pattern_5 # seqlets: 200 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-184_KLF12.UNK.0.A
metacluster_0/pattern_6 # seqlets: 176 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_7 # seqlets: 169 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-149_CTCFL.UNK.0.A, HCLUST-134_INSM1.UNK.0.A
metacluster_0/pattern_8 # seqlets: 165 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_9 # seqlets: 154 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-156_TEAD1.UNK.0.A, HCLUST-101_NFE2.UNK.0.A
metacluster_0/pattern_10 # seqlets: 153 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-159_EHF.UNK.0.A, HCLUST-130_ERG.UNK.0.A, HCLUST-129_ELF1.UNK.0.A
metacluster_0/pattern_11 # seqlets: 91 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-169_TFAP2A.UNK.0.A
metacluster_0/pattern_12 # seqlets: 76 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-167_NFKB1.UNK.0.A, HCLUST-145_RELA.UNK.0.A
metacluster_0/pattern_13 # seqlets: 77 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-137_NFIA.UNK.0.A

Pattern Name	TF Name(s)	Modisco
metacluster_0/pattern_0 # seqlets: 4530 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-101_NFE2.UNK.0.A, HCLUST-179_BACH1.UNK.0.A
metacluster_0/pattern_1 # seqlets: 885 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-170_TP53.UNK.0.A
metacluster_0/pattern_2 # seqlets: 693 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_3 # seqlets: 565 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-176_CBFB.UNK.0.A
metacluster_0/pattern_4 # seqlets: 555 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_5 # seqlets: 504 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_6 # seqlets: 452 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-137_NFIA.UNK.0.A
metacluster_0/pattern_7 # seqlets: 428 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-169_TFAP2A.UNK.0.A
metacluster_0/pattern_8 # seqlets: 366 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-175_ATF1.UNK.0.A
metacluster_0/pattern_9 # seqlets: 320 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-3_ZBTB18.UNK.0.A, HCLUST-162_ASCL1.UNK.0.A, HCLUST-138_MEIS2.UNK.0.A, HCLUST-4_ATOH1.UNK.0.A
metacluster_0/pattern_10 # seqlets: 283 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_11 # seqlets: 264 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-106_ATF6.UNK.0.A
metacluster_0/pattern_12 # seqlets: 254 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_13 # seqlets: 225 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-10_IRF1.UNK.0.A, HCLUST-9_IRF7.UNK.0.A
metacluster_0/pattern_14 # seqlets: 191 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-101_NFE2.UNK.0.A
metacluster_0/pattern_15 # seqlets: 158 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-167_NFKB1.UNK.0.A, HCLUST-145_RELA.UNK.0.A
metacluster_0/pattern_16 # seqlets: 138 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-110_GRHL2.UNK.0.A
metacluster_0/pattern_17 # seqlets: 130 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_18 # seqlets: 170 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-149_CTCFL.UNK.0.A
metacluster_0/pattern_19 # seqlets: 112 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_20 # seqlets: 103 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-16_RFX1.UNK.0.A, HCLUST-15_RFX5.UNK.0.A

Pattern Name	TF Name(s)	Modisco
metacluster_0/pattern_0 # seqlets: 6819 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-101_NFE2.UNK.0.A, HCLUST-179_BACH1.UNK.0.A
metacluster_0/pattern_1 # seqlets: 353 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_2 # seqlets: 311 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-170_TP53.UNK.0.A
metacluster_0/pattern_3 # seqlets: 292 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A
metacluster_0/pattern_4 # seqlets: 250 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-176_CBFB.UNK.0.A
metacluster_0/pattern_5 # seqlets: 125 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-175_ATF1.UNK.0.A
metacluster_0/pattern_6 # seqlets: 103 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-184_KLF12.UNK.0.A
metacluster_0/pattern_7 # seqlets: 103 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_8 # seqlets: 111 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-159_EHF.UNK.0.A, HCLUST-130_ERG.UNK.0.A, HCLUST-129_ELF1.UNK.0.A
metacluster_0/pattern_9 # seqlets: 83 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-149_CTCFL.UNK.0.A, HCLUST-134_INSM1.UNK.0.A
metacluster_0/pattern_10 # seqlets: 75 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-169_TFAP2A.UNK.0.A

Pattern Name	TF Name(s)	Modisco
metacluster_0/pattern_0 # seqlets: 2962 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-101_NFE2.UNK.0.A, HCLUST-124_FOSB.UNK.0.A, HCLUST-179_BACH1.UNK.0.A
metacluster_0/pattern_1 # seqlets: 181 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-170_TP53.UNK.0.A
metacluster_0/pattern_2 # seqlets: 113 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_3 # seqlets: 73 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-179_BACH1.UNK.0.A, HCLUST-124_FOSB.UNK.0.A, HCLUST-101_NFE2.UNK.0.A, HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_4 # seqlets: 111 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-184_KLF12.UNK.0.A
metacluster_0/pattern_5 # seqlets: 80 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-159_EHF.UNK.0.A, HCLUST-130_ERG.UNK.0.A, HCLUST-129_ELF1.UNK.0.A
metacluster_0/pattern_6 # seqlets: 105 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-175_ATF1.UNK.0.A

Pattern Name	TF Name(s)	Modisco
metacluster_0/pattern_0 # seqlets: 5497 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-101_NFE2.UNK.0.A, HCLUST-179_BACH1.UNK.0.A
metacluster_0/pattern_1 # seqlets: 500 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-170_TP53.UNK.0.A
metacluster_0/pattern_2 # seqlets: 304 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_3 # seqlets: 290 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A
metacluster_0/pattern_4 # seqlets: 261 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-176_CBFB.UNK.0.A
metacluster_0/pattern_5 # seqlets: 212 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A, HCLUST-96_NKX6-1.UNK.0.A, HCLUST-29_ZNF547.UNK.0.A
metacluster_0/pattern_6 # seqlets: 160 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-149_CTCFL.UNK.0.A, HCLUST-134_INSM1.UNK.0.A
metacluster_0/pattern_7 # seqlets: 143 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_8 # seqlets: 138 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-184_KLF12.UNK.0.A
metacluster_0/pattern_9 # seqlets: 139 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A
metacluster_0/pattern_10 # seqlets: 92 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_11 # seqlets: 90 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-101_NFE2.UNK.0.A
metacluster_0/pattern_12 # seqlets: 106 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-159_EHF.UNK.0.A, HCLUST-130_ERG.UNK.0.A, HCLUST-129_ELF1.UNK.0.A
metacluster_0/pattern_13 # seqlets: 67 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-179_BACH1.UNK.0.A, HCLUST-101_NFE2.UNK.0.A
metacluster_0/pattern_14 # seqlets: 68 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-169_TFAP2A.UNK.0.A

Pattern Name	TF Name(s)	Modisco
metacluster_0/pattern_0 # seqlets: 4723 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-101_NFE2.UNK.0.A, HCLUST-124_FOSB.UNK.0.A, HCLUST-179_BACH1.UNK.0.A
metacluster_0/pattern_1 # seqlets: 769 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-170_TP53.UNK.0.A
metacluster_0/pattern_2 # seqlets: 651 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-176_CBFB.UNK.0.A
metacluster_0/pattern_3 # seqlets: 573 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-159_EHF.UNK.0.A, HCLUST-130_ERG.UNK.0.A
metacluster_0/pattern_4 # seqlets: 525 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-156_TEAD1.UNK.0.A
metacluster_0/pattern_5 # seqlets: 512 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_6 # seqlets: 508 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-175_ATF1.UNK.0.A
metacluster_0/pattern_7 # seqlets: 489 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-137_NFIA.UNK.0.A
metacluster_0/pattern_8 # seqlets: 441 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-124_FOSB.UNK.0.A, HCLUST-156_TEAD1.UNK.0.A, HCLUST-101_NFE2.UNK.0.A, HCLUST-141_BATF3.UNK.0.A
metacluster_0/pattern_9 # seqlets: 388 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-169_TFAP2A.UNK.0.A
metacluster_0/pattern_10 # seqlets: 263 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-162_ASCL1.UNK.0.A
metacluster_0/pattern_11 # seqlets: 235 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-167_NFKB1.UNK.0.A, HCLUST-145_RELA.UNK.0.A
metacluster_0/pattern_12 # seqlets: 232 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_13 # seqlets: 210 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-10_IRF1.UNK.0.A, HCLUST-9_IRF7.UNK.0.A, HCLUST-158_BCL11A.UNK.0.A
metacluster_0/pattern_14 # seqlets: 210 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_15 # seqlets: 194 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_16 # seqlets: 90 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-149_CTCFL.UNK.0.A
metacluster_0/pattern_17 # seqlets: 90 SequenceContrib ScoresHyp_Contrib Scores
metacluster_0/pattern_18 # seqlets: 81 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-149_CTCFL.UNK.0.A
metacluster_0/pattern_19 # seqlets: 116 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-16_RFX1.UNK.0.A
metacluster_0/pattern_20 # seqlets: 67 SequenceContrib ScoresHyp_Contrib Scores	HCLUST-110_GRHL2.UNK.0.A