import basepair
import modisco

Using TensorFlow backend.
Can not use cuDNN on context None: cannot compile with cuDNN. We got this error:
b'/usr/bin/ld: cannot find -lcudnn\ncollect2: error: ld returned 1 exit status\n'
Mapped name None to device cuda: GeForce GTX TITAN X (0000:09:00.0)

from keras.models import Model, load_model
from basepair.losses import twochannel_multinomial_nll

# Use gpus 3, 5
import os
os.environ["CUDA_VISIBLE_DEVICES"] = "3, 5"

model = load_model("model.h5", custom_objects={"twochannel_multinomial_nll": twochannel_multinomial_nll})

WARNING:tensorflow:From /users/amr1/miniconda3/envs/basepair/lib/python3.6/site-packages/tensorflow/python/util/deprecation.py:497: calling conv1d (from tensorflow.python.ops.nn_ops) with data_format=NHWC is deprecated and will be removed in a future version.
Instructions for updating:
`NHWC` for data_format is deprecated, use `NWC` instead

2018-08-19 16:37:10,889 [WARNING] From /users/amr1/miniconda3/envs/basepair/lib/python3.6/site-packages/tensorflow/python/util/deprecation.py:497: calling conv1d (from tensorflow.python.ops.nn_ops) with data_format=NHWC is deprecated and will be removed in a future version.
Instructions for updating:
`NHWC` for data_format is deprecated, use `NWC` instead

WARNING:tensorflow:From /users/amr1/miniconda3/envs/basepair/lib/python3.6/site-packages/tensorflow/contrib/learn/python/learn/datasets/base.py:198: retry (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.
Instructions for updating:
Use the retry module or similar alternatives.

2018-08-19 16:37:19,198 [WARNING] From /users/amr1/miniconda3/envs/basepair/lib/python3.6/site-packages/tensorflow/contrib/learn/python/learn/datasets/base.py:198: retry (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.
Instructions for updating:
Use the retry module or similar alternatives.

from basepair.utils import read_pkl
train,valid,test = read_pkl("/users/avsec/workspace/basepair-workflow/models/0/data.pkl")

from basepair.eval import evaluate
evaluate(model, valid[0], valid[1])

{'loss': 953.9045610385498,
 'profile/Oct4_loss': 271.63486264911495,
 'profile/Sox2_loss': 162.5192502521077,
 'profile/Klf4_loss': 228.8282409395277,
 'profile/Nanog_loss': 265.30063137148727,
 'counts/Oct4_loss': 0.612707477322428,
 'counts/Sox2_loss': 0.5658578240037332,
 'counts/Klf4_loss': 0.8421652589296793,
 'counts/Nanog_loss': 0.5414273575545541}

task_names = ["profile/Oct4", "profile/Sox2", "profile/Klf4", "profile/Nanog",
              "counts/Oct4", "counts/Sox2", "counts/Klf4", "counts/Nanog"]

import keras.backend as K
inp = model.inputs[0]
fn_pos = {}
fn_neg = {}
for task_id, task_name in enumerate(task_names):
    if "counts" in task_name:
        fn_pos[task_name] = K.function([inp], K.gradients(model.outputs[task_id][:, 0], inp))
        fn_neg[task_name] = K.function([inp], K.gradients(model.outputs[task_id][:, 1], inp))
    else:
        fn_pos[task_name] = K.function([inp], K.gradients(K.sum(K.stop_gradient(K.softmax(
            model.outputs[task_id][:, :, 0])) * model.outputs[task_id][:, :, 0], axis=-1), inp))
        fn_neg[task_name] = K.function([inp], K.gradients(K.sum(K.stop_gradient(K.softmax(
            model.outputs[task_id][:, :, 1])) * model.outputs[task_id][:, :, 1], axis=-1), inp))

import numpy as np
from basepair.data import numpy_minibatch

grads_pos = {}
grads_neg = {}

for task_name in task_names:
    grads_pos[task_name] = np.concatenate([np.array(fn_pos[task_name]([batch])).squeeze()
                                           for batch in numpy_minibatch(valid[0], 512)])
    grads_neg[task_name] = np.concatenate([np.array(fn_neg[task_name]([batch])).squeeze()
                                           for batch in numpy_minibatch(valid[0], 512)])

2018-08-19 16:37:42,441 [WARNING] git-lfs not installed
2018-08-19 16:37:42,616 [WARNING] git-lfs not installed

# Setup different scores
hyp_scores = {}
scores = {}
for task_name in task_names:
    hyp_scores[task_name] = grads_pos[task_name] + grads_neg[task_name]
    hyp_scores[task_name] = hyp_scores[task_name] - hyp_scores[task_name].mean(-1, keepdims=True)
    scores[task_name] = hyp_scores[task_name] * valid[0]

onehot_data = valid[0]
task_to_scores = scores
task_to_hyp_scores = hyp_scores

import h5py
import modisco
import modisco.util
import modisco.core
import modisco.tfmodisco_workflow.seqlets_to_patterns
from modisco.tfmodisco_workflow import workflow

tfmodisco_results = modisco.tfmodisco_workflow.workflow.TfModiscoWorkflow(
                    #Slight modifications from the default settings
                    sliding_window_size=15,
                    flank_size=5,
                    target_seqlet_fdr=0.01,
                    min_seqlets_per_task=1000,
                    seqlets_to_patterns_factory=
                     modisco.tfmodisco_workflow.seqlets_to_patterns.TfModiscoSeqletsToPatternsFactory(
                        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=task_names,
                contrib_scores=task_to_scores,
                hypothetical_contribs=task_to_hyp_scores,
                one_hot=onehot_data)

On task profile/Oct4
Computing windowed sums
Computing threshold
peak(mu)= 0.06276555971801281
Mu: 6.276556e-02 +/- 5.909070e-06
Lablace_b: 0.139662450995146 and 0.26768297362318977
Thresholds: -2.282341041997075 and 2.3123230548352005
#fdrs pass: 0 and 52224
CDFs: 0.9999999489901056 and 0.9997759883296964
Est. FDRs: 0.0626664590409689 and 0.009998536975671766

Got 5362 coords
On task profile/Sox2
Computing windowed sums
Computing threshold
peak(mu)= 0.04190081917122006
Mu: 4.190082e-02 +/- 8.923907e-06
Lablace_b: 0.18997234100246468 and 0.38309747398215255
Thresholds: -2.999906593929976 and 3.62982863746956
#fdrs pass: 0 and 20393
CDFs: 0.9999998887897856 and 0.999914378505272
Est. FDRs: 0.13096704643032808 and 0.00999979911681125

Got 2483 coords
On task profile/Klf4
Computing windowed sums
Computing threshold
peak(mu)= 0.044450335949659345
Mu: 4.445034e-02 +/- 5.139410e-06
Lablace_b: 0.129423190632036 and 0.24263794268149524
Thresholds: -1.2058626884354116 and 3.8841040663421156
#fdrs pass: 7907 and 33
CDFs: 0.9999362573568685 and 0.9999998658922182
Est. FDRs: 0.00999604251532824 and 0.009426219409136976

Got 804 coords
On task profile/Nanog
Computing windowed sums
Computing threshold
peak(mu)= 0.05894492487795651
Manually adjusting thresholds to get desired num seqlets
Mu: 5.894492e-02 +/- 9.355104e-06
Lablace_b: 0.2250517944452598 and 0.6103430685661245
Thresholds: -6.540546671343967 and 6.658436521099881
#fdrs pass: 3 and 87
CDFs: 0.999999999999816 and 0.9999798590903745
Est. FDRs: 0.0015856588303443016 and 0.009948188832590832

Got 240 coords
On task counts/Oct4
Computing windowed sums
Computing threshold
peak(mu)= 0.012566710790619252
Mu: 1.256671e-02 +/- 1.162811e-05
Lablace_b: 0.3693071719980157 and 0.6606773467999582
Thresholds: -4.142623280206696 and 5.241926337560638
#fdrs pass: 0 and 83286
CDFs: 0.9999870097488164 and 0.9996348305182882
Est. FDRs: 1.0 and 0.009999506220365736

Got 7480 coords
On task counts/Sox2
Computing windowed sums
Computing threshold
peak(mu)= 0.012341707596555352
Mu: 1.234171e-02 +/- 1.008116e-05
Lablace_b: 0.3565778686238186 and 0.7367579119032016
Thresholds: -4.614141255965456 and 6.935839861283079
#fdrs pass: 0 and 19954
CDFs: 0.999997681703296 and 0.9999170481150947
Est. FDRs: 1.0 and 0.009999303991058353

Got 2908 coords
On task counts/Klf4
Computing windowed sums
Computing threshold
peak(mu)= 0.10090510392561555
Manually adjusting thresholds to get desired num seqlets
Mu: 1.009051e-01 +/- 7.118348e-06
Lablace_b: 0.43305931247688123 and 0.6413692939038103
Thresholds: -4.394236402269453 and 4.596046610120684
#fdrs pass: 0 and 0
CDFs: 0.9999689517250595 and 0.9990959839251266
Est. FDRs: 0.280650626313502 and 1.0

Got 180 coords
On task counts/Nanog
Computing windowed sums
Computing threshold
peak(mu)= 0.3062464974075556
Manually adjusting thresholds to get desired num seqlets
Mu: 3.062465e-01 +/- 8.409873e-06
Lablace_b: 0.39656970316185786 and 0.9808883116430085
Thresholds: -7.42988374017179 and 8.0423767349869
#fdrs pass: 0 and 0
CDFs: 0.9999999966275215 and 0.9996243532970951
Est. FDRs: 1.0 and 1.0

Got 252 coords
After resolving overlaps, got 9885 seqlets
Across all tasks, the weakest laplace threshold used was: 0.9990959839251266
9885 identified in total
13 activity patterns with support >= 100 out of 6561 possible patterns
Metacluster sizes:  [4500, 1107, 791, 501, 454, 400, 386, 224, 218, 206, 196, 154, 114]
Idx to activities:  {0: '1,1,0,0,1,1,0,0', 1: '1,0,0,0,1,1,0,0', 2: '1,0,0,0,1,0,0,0', 3: '1,1,0,0,1,1,0,1', 4: '-1,-1,-1,0,0,0,0,0', 5: '1,1,1,1,0,1,0,1', 6: '1,0,0,0,1,0,-1,0', 7: '1,1,0,1,1,1,0,0', 8: '0,1,0,0,1,1,0,1', 9: '0,1,0,0,0,1,0,1', 10: '0,0,-1,0,0,0,0,0', 11: '0,1,0,0,0,1,0,0', 12: '-1,0,-1,0,0,0,0,0'}
On metacluster 12
Metacluster size 114
Relevant tasks:  ('profile/Oct4', 'profile/Klf4')
Relevant signs:  (-1, -1)
(Round 1) num seqlets: 114
(Round 1) Computing coarse affmat
Beginning embedding computation
Computing embeddings
Finished embedding computation in 0.03 s
Starting affinity matrix computations
Normalization computed in 0.0 s
Cosine similarity mat computed in 0.0 s
Normalization computed in 0.0 s
Cosine similarity mat computed in 0.0 s
Finished affinity matrix computations in 0.01 s
(Round 1) Compute nearest neighbors from coarse affmat
Computed nearest neighbors in 0.0 s
(Round 1) Computing affinity matrix on nearest neighbors
Launching nearest neighbors affmat calculation job
Job completed in: 9.36 s
Launching nearest neighbors affmat calculation job
Job completed in: 0.44 s
(Round 1) Computed affinity matrix on nearest neighbors in 9.87 s
Filtered down to 102 of 114
(Round 1) Retained 102 rows out of 114 after filtering
(Round 1) Computing density adapted affmat
[t-SNE] Computing 31 nearest neighbors...
[t-SNE] Indexed 102 samples in 0.000s...
[t-SNE] Computed neighbors for 102 samples in 0.001s...
[t-SNE] Computed conditional probabilities for sample 102 / 102
[t-SNE] Mean sigma: 0.234632
(Round 1) Computing clustering
Beginning preprocessing + Louvain
Wrote graph to binary file in 0.02069568634033203 seconds
Running Louvain modularity optimization

[Parallel(n_jobs=4)]: Using backend LokyBackend with 4 concurrent workers.

from importlib import reload
import h5py
import modisco.util
reload(modisco.util)
!rm modisco_results_on_valid.hdf5
grp = h5py.File("modisco_results_on_valid.hdf5")
tfmodisco_results.save_hdf5(grp)
grp.close()

rm: cannot remove 'modisco_results_on_valid.hdf5': No such file or directory