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

from matlas.genome_data import *
from matlas.pwms import shorten_tfname, get_motifDB_id_maps
from matlas.matches import fig2vdom, vdom_pssm
from basepair.modisco.results import ModiscoResult
from vdom.helpers import (h1, p, li, img, div, b, br, ul, a,
                          details, summary,
                          table, thead, th, tr, tbody, td, ol)
from IPython.display import display


project_root = "/mnt/lab_data/kundaje/msharmin/mouse_hem/with_tfd"
oak_root = "/oak/stanford/groups/akundaje/msharmin/mouse_hem/from_labcluster/with_tfd"
srv_root = "/srv/scratch/msharmin/mouse_hem/with_tfd"
data_name = "full_mouse50"
modisco_root = "{0}/{1}/Naive_modisco2019_msharmin".format(oak_root, data_name)
modisco_root = "{0}/{1}/Naive_modisco2019".format(srv_root, data_name)

homer_root = "{0}/{1}/Naive_scans".format(srv_root, data_name)
model_root = "{0}/{1}/fineFactorized".format(srv_root, data_name)
basset_root = "{0}/{1}/Naive_basset_motifs".format(srv_root, data_name)
reportfile= "/mnt/lab_data/kundaje/msharmin/mouse_hem/with_tfd/full_mouse50/filtering samples_MS2.xlsx"
sheetname = "filter23"

load data from labcluster

Using TensorFlow backend.
2019-09-04 16:52:56,824 [WARNING] git-lfs not installed

from matlas.matches import DenovoHomer, DenovoModisco
from matlas.utils import rolling_window
from matlas.dlutils import load_deeplift_data

task_idx = 273
modisco_dir = "{0}/task_{1}-naivegw".format(modisco_root, task_idx)
ob = DenovoModisco(modisco_dir)
ob.load_matched_motifs()
ob.get_motif_per_celltype()

TF-MoDISco is using the TensorFlow backend.

from matlas.utils import get_calibrated_peaks

summit_file = "{0}/{1}/fineFactorized/task_{2}-naivegw/predictions/task_{2}-{3}Summit.calibrated.{1}.tab.gz".format(project_root, data_name, task_idx, "Naive")
perf_file = "{0}/{1}/fineFactorized/task_{2}-naivegw/{3}auPRC.txt".format(project_root, data_name, task_idx, "Naive")
calibrated_data = get_calibrated_peaks(summit_file, perf_file)
selected_summits = calibrated_data['selected_summits']
indices = calibrated_data['indices']

deeplift_hdf = "{0}/{1}/Naive_deeplift2019/task_{2}-naivegw/summit.h5".format(srv_root, data_name, task_idx)
deeplift_data = load_deeplift_data(deeplift_hdf, keys=['one_hot', 'hyp_scores', 'scores', 
                                                       'shuffled_onehot', 'shuffled_scores', 'peaks'])
deeplift_data.keys()

odict_keys(['one_hot', 'shuffled_onehot', 'peaks', 'hyp_scores', 'scores', 'shuffled_scores'])

def get_probs(arr_weights):
    """
    convert the array weights to probability
    """
    arr_weights = arr_weights / arr_weights.sum(1, keepdims=True)
    arr_weights = arr_weights + 0.01  # add pseudo-counts
    arr_weights = arr_weights / arr_weights.sum(1, keepdims=True)
    return arr_weights

def get_per_base_ic(arr_weights, arr_type='weights'):
    """Get information content per each position per base pair
    """
    if(arr_type=='weights'):
        probs = get_probs(arr_weights)
    else:
        probs = arr_weights
    bg_ic = sum(-np.array(BP_BACKGROUND)*np.log2(np.array(BP_BACKGROUND)))
    ic = bg_ic + np.sum(probs * np.log2(probs), axis=1)

    return ic

def idx_to_chomp(arr_weights, arr_type='probabilities'):
    """
    leading and trailing bases with low IC is droppped
    """
    ic_values = get_per_base_ic(arr_weights, arr_type)
    
    true_indices = np.argwhere(ic_values > 0.3)
    start_idx = true_indices[0][0]
    end_idx = true_indices[-1][0]
    min_pwm_len = sum(ic_values>1)
    
    return start_idx, end_idx

def chomp_pwm(arr_weights, arr_type='probabilities'):
    """
    leading and trailing bases with low IC is droppped
    """
    start_idx, end_idx = idx_to_chomp(arr_weights, arr_type)
    if(arr_type=='weights'):
        return arr_weights[start_idx:end_idx+1, :]
    else:
        return probs[start_idx:end_idx+1, :]
    
    
def pwm2pssm(arr_weights, arr_type='probabilities', background_probs=BP_BACKGROUND, ):
    """Convert the probanility array to pssm array
    In probanility array: rows sum to one
    """
    if(arr_type=='weights'):
        probs = get_probs(arr_weights)
    else:
        probs = arr_weights
    b = np.array(background_probs)[np.newaxis]
    return np.log(arr_weights / b).astype(arr_weights.dtype)

def load_modisco_pattern(tf_idx):
    weights = ob.denovo_pwms[pattern_names[tf_idx]]['PPM']
    #weights = weights[::-1, ::-1]
    start_idx, end_idx = idx_to_chomp(weights, 'weights')
    probs = get_probs(weights)
    modisco_weights_to_scan = weights[start_idx:end_idx+1, :]
    modisco_pfm_to_scan = probs[start_idx:end_idx+1, :]
    modisco_pwm_to_scan = pwm2pssm(modisco_pfm_to_scan)

    display(summary('modisco_pattern_weights', vdom_pssm(weights), 
            'modisco_pattern_probability', vdom_pssm(probs), 
            'chomped modisco_pattern_weights', vdom_pssm(modisco_weights_to_scan),
            'chomped modisco_pattern_probability', vdom_pssm(modisco_pfm_to_scan),
            'chomped modisco_pattern_pwm', vdom_pssm(modisco_pwm_to_scan),
           ))
    return modisco_weights_to_scan, modisco_pfm_to_scan, modisco_pwm_to_scan

pattern_names = ['metacluster_1/pattern_7', 'metacluster_1/pattern_0', 'metacluster_1/pattern_1']
tfnames = ['GATA', 'CTCF', 'SPI1']
weights_query, pfm_query, pwm_query = load_modisco_pattern(tf_idx = 0)

from basepair.plot.profiles import extract_signal
def get_seqlet_scores(pattern_name):
    
    seqlets_per_pattern = ob.mr.seqlets()
    seqlets_ = seqlets_per_pattern[pattern_name]
    contrib_per_pattern = extract_signal(deeplift_data['scores'][indices], seqlets_)
    seqs_per_pattern = extract_signal(deeplift_data['one_hot'][indices], seqlets_)
    return contrib_per_pattern, seqs_per_pattern

#get the scores for seqlets
seqlet_scores, seqlet_onehot = get_seqlet_scores(pattern_names[0])

from modisco.visualization import viz_sequence

#check the data
for seqno in range(5):
    viz_sequence.plot_weights(seqlet_scores[seqno],
                              #highlight={"black": [(start_motif, end_motif)]},
                              subticks_frequency=10, figsize=(20,1))

from matlas.sliding_similarities import sliding_similarity
def pwm_scan(pwm, array):
    strided_arr = rolling_window(array.transpose((0,2,1)), window=len(pwm)).transpose((0,2,3,1))
    pwm_scores = np.sum(strided_arr*pwm[None,None,:,:],axis=(2,3))
    return pwm_scores

method_to_scores = OrderedDict()
# method_to_scores['pwm_multisum'] = pwm_scan(pwm=pwm_query, array=deeplift_data['one_hot']) 
# method_to_scores['pwm_convolution'] = sliding_similarity(pwm_query, deeplift_data['one_hot'], 
#                                                         'dotproduct', None, n_jobs=10, verbose=False)

# In sliding_similarity, convolution between pwm and one_hot encoding was used
# testing convolution of pwm*one_hot product with av's pwm_scan. Both produces same scores
# np.max(np.abs(method_to_scores['pwm_multisum']-method_to_scores['pwm_convolution']))

def sliding_sum(query, target_scores):
    query = query.swapaxes(-2, -1)  #
    target_scores = target_scores.swapaxes(-2, -1)

    assert target_scores.shape[-1] >= query.shape[-1]  # target needs to be longer than the query
    window_len = query.shape[-1]

    ta_strided = rolling_window(target_scores, window_len).swapaxes(-2, -1)
    ta_L1_norm = np.sum((ta_strided), axis=(-3, -2))
    
    return ta_L1_norm

# method_to_scores['sum_score'] = sliding_sum(pwm_query, deeplift_data['scores'])

# method_to_scores['pfm_convolution'] = sliding_similarity(pfm_query, deeplift_data['one_hot'], 
#                                                         'dotproduct', None, n_jobs=10, verbose=False)
# method_to_scores['weight_convolution'] = sliding_similarity(weights_query, deeplift_data['one_hot'], 
#                                                         'dotproduct', None, n_jobs=10, verbose=False)

def sliding_masked_convolution(query, target_scores):#imp_scores, pwm):
    strided_arr = rolling_window(target_scores.transpose((0,2,1)),
                                 window=len(query)).transpose((0,2,3,1))
    # this finds the cosine similarity with a masked version of the pwm
    # where only the positions that are nonzero in the deeplift scores are
    # considered
    dot_product_pwm = np.sum(strided_arr*query[None,None,:,:], axis=(2,3)) #doproduct(pwm, score)
    
    norm_deeplift_scores = np.sqrt(np.sum(np.square(strided_arr), axis=(2,3))) #norm(sum_score)
    norm_masked_pwm = np.sqrt(np.sum(np.square((np.abs(strided_arr)>0)*query[None,None,:,:]), axis=(2,3))) #masked_dot
    cosine_sim = dot_product_pwm/(norm_deeplift_scores*norm_masked_pwm)
    return cosine_sim

method_to_seqlet_scores = OrderedDict()
#get sum of deeplift scores
method_to_seqlet_scores['sum_score'] = sliding_sum(pwm_query, seqlet_scores)

#get convolution of pwm with deeplift
method_to_seqlet_scores['pwm_convolution'] = sliding_similarity(pwm_query, seqlet_scores, 
                                                        'dotproduct', None, n_jobs=10, verbose=False)
method_to_seqlet_scores['pfm_convolution'] = sliding_similarity(pfm_query, seqlet_scores, 
                                                        'dotproduct', None, n_jobs=10, verbose=False)
method_to_seqlet_scores['weight_convolution'] = sliding_similarity(weights_query, seqlet_scores, 
                                                        'dotproduct', None, n_jobs=10, verbose=False)

#get convolution of pwm with sequence
method_to_seqlet_scores['pwm_scan'] = sliding_similarity(pwm_query, seqlet_onehot, 
                                                        'dotproduct', None, n_jobs=10, verbose=False)
method_to_seqlet_scores['pfm_scan'] = sliding_similarity(pfm_query, seqlet_onehot, 
                                                        'dotproduct', None, n_jobs=10, verbose=False)
method_to_seqlet_scores['weight_scan'] = sliding_similarity(weights_query, seqlet_onehot, 
                                                        'dotproduct', None, n_jobs=10, verbose=False)
#get jaccard simillarity of pwm with deeplift
method_to_seqlet_scores['pwm_jaccard'] = sliding_similarity(pwm_query, seqlet_scores, 
                                                        'continousjaccard', None, n_jobs=10, verbose=False)
method_to_seqlet_scores['pfm_jaccard'] = sliding_similarity(pfm_query, seqlet_scores, 
                                                        'continousjaccard', None, n_jobs=10, verbose=False)
method_to_seqlet_scores['weight_jaccard'] = sliding_similarity(weights_query, seqlet_scores, 
                                                        'continousjaccard', None, n_jobs=10, verbose=False)


#get masked convolution of pwm with deeplift
# method_to_seqlet_scores['pwm_masked_convolution'] = sliding_similarity(pwm_query, seqlet_scores, 
#                                                         'masked_dotproduct', None, n_jobs=10, verbose=False)
# method_to_seqlet_scores['pfm_masked_convolution'] = sliding_similarity(pfm_query, seqlet_scores, 
#                                                         'masked_dotproduct', None, n_jobs=10, verbose=False)
# method_to_seqlet_scores['weight_masked_convolution'] = sliding_similarity(weights_query, seqlet_scores, 
#                                                         'masked_dotproduct', None, n_jobs=10, verbose=False)

method_to_seqlet_scores['pwm_masked_convolution'] = sliding_masked_convolution(pwm_query, seqlet_scores)
method_to_seqlet_scores['pfm_masked_convolution'] = sliding_masked_convolution(pfm_query, seqlet_scores)
method_to_seqlet_scores['weight_masked_convolution'] = sliding_masked_convolution(weights_query, seqlet_scores)

def get_hits(method_to_scores, min_tot_imp, keys=None):
    if keys is None:
        keys = list(method_to_scores.keys())
    sum_scores = method_to_scores['sum_score']
    sumscores_mask_idx1, sumscores_mask_idx2 = np.nonzero(sum_scores > min_tot_imp)
    
    matches = [(sum_scores[idx1, idx2], idx1, idx2) for (idx1, idx2) in zip(sumscores_mask_idx1, sumscores_mask_idx2)]
    match_df = pd.DataFrame(np.array(matches))
    match_df.columns = ['sum_score', 'seqno', 'start_idx']
    match_df.seqno = match_df.seqno.values.astype(int)
    match_df.start_idx = match_df.start_idx.values.astype(int)
    for key in keys:
        if(key=='sum_score'):
            continue
        match_df[key] = [method_to_scores[key][idx1, idx2] for 
                         (idx1, idx2) in zip(sumscores_mask_idx1, sumscores_mask_idx2)]
    
    return match_df

# hits = get_hits(method_to_scores, min_tot_imp=0.125*len(pwm_query),
#                 keys=set(method_scores.keys()).difference(set(['seqno', 'start_idx'])))
# hits.shape

seqlet_hits = get_hits(method_to_seqlet_scores, min_tot_imp=0.125*len(pwm_query),
                       keys=set(method_to_seqlet_scores.keys()).difference(set(['seqno', 'start_idx'])))
seqlet_hits.shape, seqlet_hits.keys()

((6810, 15), Index(['sum_score', 'seqno', 'start_idx', 'weight_convolution',
        'pfm_convolution', 'pwm_convolution', 'pwm_scan',
        'pfm_masked_convolution', 'pwm_jaccard', 'pfm_scan', 'pfm_jaccard',
        'weight_jaccard', 'pwm_masked_convolution', 'weight_masked_convolution',
        'weight_scan'],
       dtype='object'))

import seaborn as sns
def plots(df):
#     dfm = df.melt(var_name='columns')
#     g = sns.FacetGrid(dfm, col='columns')
#     g = (g.map(sns.distplot, 'value'))
    g = sns.pairplot(df, diag_kind="kde", kind='scatter')
    return None

plots(seqlet_hits[['sum_score', 
                   'weight_convolution', 'pfm_convolution', 'pwm_convolution' 
                ]])

/users/msharmin/anaconda2/envs/basepair/lib/python3.6/site-packages/scipy/stats/stats.py:1713: 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.

plots(seqlet_hits[[
                   'weight_scan', 'pfm_scan', 'pwm_scan',
                ]])

/users/msharmin/anaconda2/envs/basepair/lib/python3.6/site-packages/scipy/stats/stats.py:1713: 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.

plots(seqlet_hits[[
                   'weight_jaccard', 'pfm_jaccard', 'pwm_jaccard', 
                ]])
plots(seqlet_hits[[
    'weight_masked_convolution', 'pfm_masked_convolution', 'pwm_masked_convolution'
                ]])

/users/msharmin/anaconda2/envs/basepair/lib/python3.6/site-packages/scipy/stats/stats.py:1713: 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.

import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
import seaborn as sns; sns.set()
from matlas.seaborn_utils import SeabornFig2Grid
gs = gridspec.GridSpec(2, 2)
g0 = sns.jointplot("pwm_convolution", "pfm_convolution", data=seqlet_hits, kind="kde", color="skyblue")
g1 = sns.jointplot("pwm_convolution", "weight_convolution", data=seqlet_hits, kind="kde", color="skyblue")
g2 = sns.jointplot("pwm_convolution", "pwm_scan", data=seqlet_hits, kind="kde", color="skyblue")
g3 = sns.jointplot("pwm_convolution", "pwm_jaccard", data=seqlet_hits, kind="kde", color="skyblue")
# g0 = sns.jointplot("pwm_convolution", "pfm_jaccard", data=seqlet_hits, kind="kde", color="skyblue")
# g1 = sns.jointplot("pwm_convolution", "weight_jaccard", data=seqlet_hits, kind="kde", color="skyblue")
# g2 = sns.jointplot("pwm_convolution", "pwm_masked_convolution", data=seqlet_hits, kind="kde", color="skyblue")
# g3 = sns.jointplot("pwm_convolution", "pfm_masked_convolution", data=seqlet_hits, kind="kde", color="skyblue")

fig = plt.figure(figsize=(10,8))
mg0 = SeabornFig2Grid(g0, fig, gs[0])
mg1 = SeabornFig2Grid(g1, fig, gs[1])
mg2 = SeabornFig2Grid(g2, fig, gs[2])
mg3 = SeabornFig2Grid(g3, fig, gs[3])

gs.tight_layout(fig)
plt.show()

/users/msharmin/anaconda2/envs/basepair/lib/python3.6/site-packages/scipy/stats/stats.py:1713: 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.

gs = gridspec.GridSpec(2, 2)
g0 = sns.jointplot("pwm_convolution", "pfm_jaccard", data=seqlet_hits, kind="kde", color="skyblue")
g1 = sns.jointplot("pwm_convolution", "weight_jaccard", data=seqlet_hits, kind="kde", color="skyblue")
g2 = sns.jointplot("pwm_convolution", "pwm_masked_convolution", data=seqlet_hits, kind="kde", color="skyblue")
g3 = sns.jointplot("pwm_convolution", "pfm_masked_convolution", data=seqlet_hits, kind="kde", color="skyblue")

fig = plt.figure(figsize=(10,8))
mg0 = SeabornFig2Grid(g0, fig, gs[0])
mg1 = SeabornFig2Grid(g1, fig, gs[1])
mg2 = SeabornFig2Grid(g2, fig, gs[2])
mg3 = SeabornFig2Grid(g3, fig, gs[3])

gs.tight_layout(fig)
plt.show()

/users/msharmin/anaconda2/envs/basepair/lib/python3.6/site-packages/scipy/stats/stats.py:1713: 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.

# include really low scoring examples
seqlet_hits = get_hits(method_to_seqlet_scores, min_tot_imp=0.0625*len(pwm_query),
                       keys=['sum_score', 
                             'pwm_scan', 'pfm_scan', 'weight_scan', 
                             'pwm_jaccard', 'pfm_jaccard', 'weight_jaccard', 
                             'pwm_convolution', 'pfm_convolution', 'weight_convolution',
                             'pwm_masked_convolution', 'pfm_masked_convolution', 'weight_masked_convolution'
                            ])
seqlet_hits.shape

(13587, 15)

from sklearn.cluster import SpectralClustering

def select_hits(df):
    x = df[['sum_score', 
            'pwm_convolution', 
            'pfm_convolution', 
            'weight_convolution',
            'pwm_masked_convolution', 
            'pfm_masked_convolution', 
            'weight_masked_convolution',
            'pwm_jaccard', 
            'pfm_jaccard', 
            'weight_jaccard', 
            'pwm_scan', 
            'pfm_scan', 
            'weight_scan'
           ]]
    #TODO: try umap clustering
    clustering = SpectralClustering(n_clusters=2, assign_labels="discretize", random_state=0).fit(x.values)

    mean0 = np.mean(df['pwm_convolution'][clustering.labels_==0])
    mean1 = np.mean(df['pwm_convolution'][clustering.labels_==1])
    cno = 1
    if(mean0 > mean1):
        cno = 0
    df = df[clustering.labels_==cno]
    df = df[df['pwm_jaccard']>0]
    return df

filtered_seqlet_hits = select_hits(seqlet_hits)
filtered_seqlet_hits.shape

(953, 15)

#exclude the overlapping ones
def exclude_overlapping_hits(df, pwm):
    from scipy.cluster.hierarchy import fcluster, linkage
    keep = []
    for seqno, seq_df in df.groupby('seqno'):
        if(seq_df.shape[0]==1):
            keep.append(seq_df.index[0])
            continue
        starts_to_merge = seq_df.start_idx.values
        starts_before_cluster = starts_to_merge.reshape((-1,1)).tolist()
        z = linkage(starts_before_cluster, 'complete')
        max_len_per_motif = len(pwm)
        clusters = fcluster(z, max_len_per_motif, criterion='distance')
        seq_df = seq_df.join(pd.DataFrame({'cluster': clusters}, index=seq_df.index).astype(int))
        #print(seq_df)
        for seqno, seq_df in seq_df.groupby('cluster'):
            if(seq_df.shape[0]==1):
                keep.append(seq_df.index[0])
                continue
            candidate_indx = seq_df.index[np.argmax(seq_df['pwm_scan'].values)]
            keep.append(candidate_indx)
    df = df.loc[keep, :]
    df = df.sort_values(['seqno', 'start_idx'])
    return df

#exclude_overlapping_hits(filtered_seqlet_hits[filtered_seqlet_hits.seqno==0], pwm_query)
filtered_seqlet_hits2 = exclude_overlapping_hits(filtered_seqlet_hits, pwm_query)
filtered_seqlet_hits2.shape

(917, 15)

filtered_seqlet_hits2[filtered_seqlet_hits2.seqno==714]

# check any seqlet is missing
def check_summary(df):
    print(df.shape)
    print(df.columns)
    for key in df.columns:
        print(key, np.min(df[key]))
    
    missing_seqs = set(range(seqlet_onehot.shape[0])).difference(set(df.seqno))
    print(missing_seqs)
    return None
check_summary(filtered_seqlet_hits)
check_summary(filtered_seqlet_hits2)
#0.791953	578	20	1.521076	0.898678	1.814853

(953, 15)
Index(['sum_score', 'seqno', 'start_idx', 'pwm_scan', 'pfm_scan',
       'weight_scan', 'pwm_jaccard', 'pfm_jaccard', 'weight_jaccard',
       'pwm_convolution', 'pfm_convolution', 'weight_convolution',
       'pwm_masked_convolution', 'pfm_masked_convolution',
       'weight_masked_convolution'],
      dtype='object')
sum_score 0.5294555804808624
seqno 0
start_idx 2
pwm_scan 1.2664481252093567
pfm_scan 4.629802247582542
weight_scan 8.156709088644678
pwm_jaccard 0.0005829842385300647
pfm_jaccard 0.06520915295226179
weight_jaccard 0.004442358236267636
pwm_convolution -0.9119619644006709
pfm_convolution 0.2714236631550992
weight_convolution 0.041829735601686874
pwm_masked_convolution -0.18962783250233317
pfm_masked_convolution 0.26093742514812696
weight_masked_convolution 0.021984330015196046
set()
(917, 15)
Index(['sum_score', 'seqno', 'start_idx', 'pwm_scan', 'pfm_scan',
       'weight_scan', 'pwm_jaccard', 'pfm_jaccard', 'weight_jaccard',
       'pwm_convolution', 'pfm_convolution', 'weight_convolution',
       'pwm_masked_convolution', 'pfm_masked_convolution',
       'weight_masked_convolution'],
      dtype='object')
sum_score 0.5294555804808624
seqno 0
start_idx 2
pwm_scan 1.335839857548609
pfm_scan 4.629802247582542
weight_scan 8.267396133835705
pwm_jaccard 0.009196646206622246
pfm_jaccard 0.2071279446085702
weight_jaccard 0.21956038132732028
pwm_convolution -0.10803838253418951
pfm_convolution 0.49349497595970093
weight_convolution 0.6961941237791255
pwm_masked_convolution -0.0410643408158796
pfm_masked_convolution 0.43906283124540557
weight_masked_convolution 0.36458707939615875
set()

def view_test_case(scores, one_hot, df, pwm, seqno, start=400, end=600):
    df = df[df.seqno==seqno]
    if(df.shape[0]==0):
        print('Not part of the hits')
        viz_sequence.plot_weights(scores[seqno, start:end,:],
                                  highlight={"black": [(21, 29)]},
                                  subticks_frequency=10, figsize=(20,1))

        viz_sequence.plot_weights(one_hot[seqno, start:end, :], 
                                  highlight={"black": [(21, 29)]},
                                  subticks_frequency=10, figsize=(20,1))
        return None
    keys = set(seqlet_hits.columns).difference(set(['seqno', 'start_idx']))
    for i, row in df.iterrows():
        printables = ["{}: {}".format(key, round(row[key], 3)) for key in keys]
        seqno = int(row['seqno'])
        start_idx = int(row['start_idx'])
        print(i,
              printables,
              "coordinates:",(seqno, start_idx))
        start_motif = start_idx - start
        end_motif = start_motif+len(pwm)
        print(start_idx, start_motif, end_motif)
        if(end_motif > (end-start) or end_motif<0):
            continue
        viz_sequence.plot_weights(scores[seqno, start:end,:],
                                  highlight={"black": [(start_motif, end_motif)]},
                                  subticks_frequency=10, figsize=(20,1))

        viz_sequence.plot_weights(one_hot[seqno, start:end, :], 
                                  highlight={"black": [(start_motif, end_motif)]},
                                  subticks_frequency=10, figsize=(20,1))
    return None
view_test_case(seqlet_scores, seqlet_onehot, filtered_seqlet_hits2, pwm_query, seqno=411, start=0, end=70)

8146 ['weight_convolution: 2.24', 'pfm_convolution: 1.164', 'sum_score: 1.297', 'pwm_convolution: 1.591', 'pwm_scan: 8.58', 'pfm_masked_convolution: 0.917', 'pwm_jaccard: 0.068', 'pfm_scan: 6.335', 'pfm_jaccard: 0.506', 'weight_jaccard: 0.652', 'pwm_masked_convolution: 0.909', 'weight_masked_convolution: 0.946', 'weight_scan: 10.689'] coordinates: (411, 20)
20 20 28

def view_test_cases(scores, one_hot, df, pwm, start=400, end=600):
    max_examples_to_view=5
    keys = set(seqlet_hits.columns).difference(set(['seqno', 'start_idx']))
    for i, row in df.iterrows():
        printables = ["{}: {}".format(key, round(row[key], 3)) for key in keys]
        seqno = int(row['seqno'])
        start_idx = int(row['start_idx'])
        print(i,
              printables,
              "coordinates:",(seqno, start_idx))
        start_motif = start_idx - start
        end_motif = start_motif+len(pwm)
        print(start_idx, start_motif, end_motif)
        if(end_motif > (end-start) or end_motif<0):
            continue
        viz_sequence.plot_weights(scores[seqno, start:end,:],
                                  highlight={"black": [(start_motif, end_motif)]},
                                  subticks_frequency=10, figsize=(20,1))

        viz_sequence.plot_weights(one_hot[seqno, start:end, :], 
                                  highlight={"black": [(start_motif, end_motif)]},
                                  subticks_frequency=10, figsize=(20,1))
        max_examples_to_view -= 1
#         if(max_examples_to_view==0):
#             break
            
    return None

view_test_cases(seqlet_scores, seqlet_onehot, filtered_seqlet_hits2, pwm_query, start=0, end=70)

4 ['weight_convolution: 4.215', 'pfm_convolution: 2.205', 'sum_score: 2.382', 'pwm_convolution: 3.044', 'pwm_scan: 9.189', 'pfm_masked_convolution: 0.933', 'pwm_jaccard: 0.073', 'pfm_scan: 6.621', 'pfm_jaccard: 0.547', 'weight_jaccard: 0.748', 'pwm_masked_convolution: 0.927', 'weight_masked_convolution: 0.974', 'weight_scan: 10.975'] coordinates: (0, 20)
20 20 28

24 ['weight_convolution: 1.331', 'pfm_convolution: 0.758', 'sum_score: 1.325', 'pwm_convolution: -0.108', 'pwm_scan: 1.773', 'pfm_masked_convolution: 0.641', 'pwm_jaccard: 0.009', 'pfm_scan: 4.862', 'pfm_jaccard: 0.232', 'weight_jaccard: 0.26', 'pwm_masked_convolution: -0.041', 'weight_masked_convolution: 0.595', 'weight_scan: 8.39'] coordinates: (0, 50)
50 50 58

33 ['weight_convolution: 5.405', 'pfm_convolution: 2.879', 'sum_score: 3.144', 'pwm_convolution: 4.006', 'pwm_scan: 9.189', 'pfm_masked_convolution: 0.89', 'pwm_jaccard: 0.073', 'pfm_scan: 6.621', 'pfm_jaccard: 0.543', 'weight_jaccard: 0.635', 'pwm_masked_convolution: 0.891', 'weight_masked_convolution: 0.912', 'weight_scan: 10.975'] coordinates: (1, 20)
20 20 28

62 ['weight_convolution: 5.591', 'pfm_convolution: 2.901', 'sum_score: 3.177', 'pwm_convolution: 4.004', 'pwm_scan: 8.58', 'pfm_masked_convolution: 0.926', 'pwm_jaccard: 0.068', 'pfm_scan: 6.335', 'pfm_jaccard: 0.552', 'weight_jaccard: 0.708', 'pwm_masked_convolution: 0.927', 'weight_masked_convolution: 0.958', 'weight_scan: 10.689'] coordinates: (2, 20)
20 20 28

100 ['weight_convolution: 7.166', 'pfm_convolution: 3.73', 'sum_score: 4.044', 'pwm_convolution: 5.181', 'pwm_scan: 8.871', 'pfm_masked_convolution: 0.899', 'pwm_jaccard: 0.071', 'pfm_scan: 6.47', 'pfm_jaccard: 0.532', 'weight_jaccard: 0.642', 'pwm_masked_convolution: 0.9', 'weight_masked_convolution: 0.934', 'weight_scan: 10.923'] coordinates: (3, 20)
20 20 28

131 ['weight_convolution: 3.916', 'pfm_convolution: 2.051', 'sum_score: 2.248', 'pwm_convolution: 2.824', 'pwm_scan: 8.871', 'pfm_masked_convolution: 0.939', 'pwm_jaccard: 0.071', 'pfm_scan: 6.47', 'pfm_jaccard: 0.564', 'weight_jaccard: 0.745', 'pwm_masked_convolution: 0.933', 'weight_masked_convolution: 0.97', 'weight_scan: 10.923'] coordinates: (4, 20)
20 20 28

174 ['weight_convolution: 6.193', 'pfm_convolution: 3.244', 'sum_score: 3.775', 'pwm_convolution: 4.362', 'pwm_scan: 8.263', 'pfm_masked_convolution: 0.901', 'pwm_jaccard: 0.065', 'pfm_scan: 6.184', 'pfm_jaccard: 0.522', 'weight_jaccard: 0.61', 'pwm_masked_convolution: 0.881', 'weight_masked_convolution: 0.913', 'weight_scan: 10.636'] coordinates: (5, 20)
20 20 28

196 ['weight_convolution: 3.782', 'pfm_convolution: 1.964', 'sum_score: 2.184', 'pwm_convolution: 2.69', 'pwm_scan: 7.616', 'pfm_masked_convolution: 0.853', 'pwm_jaccard: 0.06', 'pfm_scan: 6.033', 'pfm_jaccard: 0.467', 'weight_jaccard: 0.514', 'pwm_masked_convolution: 0.841', 'weight_masked_convolution: 0.866', 'weight_scan: 10.584'] coordinates: (6, 20)
20 20 28

235 ['weight_convolution: 6.409', 'pfm_convolution: 3.276', 'sum_score: 3.503', 'pwm_convolution: 4.523', 'pwm_scan: 8.58', 'pfm_masked_convolution: 0.902', 'pwm_jaccard: 0.068', 'pfm_scan: 6.335', 'pfm_jaccard: 0.473', 'weight_jaccard: 0.642', 'pwm_masked_convolution: 0.904', 'weight_masked_convolution: 0.947', 'weight_scan: 10.689'] coordinates: (7, 20)
20 20 28

257 ['weight_convolution: 3.529', 'pfm_convolution: 1.822', 'sum_score: 1.957', 'pwm_convolution: 2.51', 'pwm_scan: 8.225', 'pfm_masked_convolution: 0.932', 'pwm_jaccard: 0.065', 'pfm_scan: 6.32', 'pfm_jaccard: 0.535', 'weight_jaccard: 0.711', 'pwm_masked_convolution: 0.918', 'weight_masked_convolution: 0.971', 'weight_scan: 10.87'] coordinates: (8, 20)
20 20 28

# tmp_hits = seqlet_hits[(seqlet_hits['pfm_convolution']>=1) & (seqlet_hits['pwm_convolution']>=0)]
# tmp_hits[tmp_hits.seqno==589]
# seqlet_hits[seqlet_hits.seqno==589]
def print_hits(seqno, df):
    df = df[df.seqno==seqno]
    return df

seqno = 411
print_hits(seqno, seqlet_hits)
#print_hits(seqno, filtered_seqlet_hits)
#print_hits(seqno, filtered_seqlet_hits2)

def print_positional_metadata(seqno, position):
    print("seqno: ", seqno, ", position: ", position)
    for key in method_to_seqlet_scores.keys():
        print(key, ": ", method_to_seqlet_scores[key][seqno, position])
    return None

print_positional_metadata(411, 50)

seqno:  411 , position:  50
sum_score :  0.49528664289973673
pwm_convolution :  -0.9707458526143764
pfm_convolution :  0.11525304126095902
weight_convolution :  0.19439835836030966
pwm_scan :  -15.951820489446279
pfm_scan :  1.2235821935708682
weight_scan :  1.7877666263372463
pwm_jaccard :  -0.13061477707862001
pfm_jaccard :  0.08280641557377288
weight_jaccard :  0.08207471727288053
pwm_masked_convolution :  -0.6653442622974229
pfm_masked_convolution :  0.5990337238706755
weight_masked_convolution :  0.5733796309913389

#query_weights on deeplift
#query pfm on deeplift
#query pwm on deeplift
#query_hyp on deeplift

#query_weights on deeplift
#query pfm on deeplift
#query pwm on deeplift
#query_hyp on deeplift

import modisco.visualization
from modisco.visualization import viz_sequence
for i in range(10):
    print(i)
    print('regression model trained with naive peaks')
    viz_sequence.plot_weights(naive_scores[i, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
    print('regression model trained with idr peaks')
    viz_sequence.plot_weights(idr_scores[i, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
    print('binary model trained with idr peaks as 1 and naive peaks as -1')
    viz_sequence.plot_weights(binary_scores[i, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
    

#print(hyp_scores[0, :4, :])
#print(scores[0, :4, :])
#viz_sequence.plot_weights(naive_scores[0, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
#print(0)
#viz_sequence.plot_weights(idr_scores[89214, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
#viz_sequence.plot_weights(binary_scores[0, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
#print(3)
#viz_sequence.plot_weights(idr_scores[90511, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})
#viz_sequence.plot_weights(binary_scores[3, :, :], colors={0:'red', 1:'blue', 2:'orange', 3:'green'})

0
regression model trained with naive peaks

regression model trained with idr peaks

binary model trained with idr peaks as 1 and naive peaks as -1

1
regression model trained with naive peaks

regression model trained with idr peaks

binary model trained with idr peaks as 1 and naive peaks as -1

2
regression model trained with naive peaks

regression model trained with idr peaks

binary model trained with idr peaks as 1 and naive peaks as -1

3
regression model trained with naive peaks

	sum_score	seqno	start_idx	pwm_scan	pfm_scan	weight_scan	pwm_jaccard	pfm_jaccard	weight_jaccard	pwm_convolution	pfm_convolution	weight_convolution	pwm_masked_convolution	pfm_masked_convolution	weight_masked_convolution
8143	0.731279	411	17	-4.610036	3.220707	5.001780	-0.037773	0.152309	0.161939	-0.301083	0.292727	0.454521	-0.151113	0.497255	0.425820
8144	1.008944	411	18	-4.971968	3.587682	4.977184	-0.041432	0.162063	0.133626	-1.182558	0.361936	0.486159	-0.406227	0.453763	0.373196
8145	1.233264	411	19	-14.423033	1.433749	2.192847	-0.123225	0.084110	0.069878	-3.097933	0.158248	0.282611	-0.831991	0.288853	0.266391
8146	1.296638	411	20	8.579902	6.335090	10.688832	0.068057	0.505689	0.652191	1.591254	1.164406	2.240044	0.909064	0.916588	0.946240
8147	1.364044	411	21	-15.305287	1.373857	1.840224	-0.131670	0.075323	0.052173	-3.355454	0.147573	0.208874	-0.845576	0.271204	0.226026
8148	1.370298	411	22	-15.136421	1.555710	2.192695	-0.130826	0.107704	0.102571	-2.227585	0.372027	0.622283	-0.554834	0.625662	0.569018
8149	1.293506	411	23	-6.011615	3.429785	5.845349	-0.050320	0.182729	0.164286	-1.122411	0.508080	0.846865	-0.351520	0.548245	0.481960
8150	1.230326	411	24	-19.920621	0.220664	0.136918	-0.173633	0.013984	0.005843	-2.488539	0.050517	0.034333	-0.641842	0.981765	0.829766
8151	0.938041	411	25	-17.122085	1.242094	0.756873	-0.134464	0.084165	0.033177	-0.785398	0.319872	0.193643	-0.237843	0.884913	0.749474
8152	0.649387	411	26	-12.027125	1.998911	2.733568	-0.102334	0.051859	0.031008	-0.946346	0.111807	0.076402	-0.457633	0.309741	0.122291