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-05 11:06:02,273 [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 = 2)

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[2])

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_sum, sliding_similarity
def collect_various_scoring(pwm_query, scores, one_hot):
    method_to_scores = OrderedDict()
    #get sum of deeplift scores
    method_to_scores['sum_score'] = sliding_sum(pwm_query, scores)

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

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


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

#     method_to_scores['pwm_masked_convolution'] = sliding_masked_convolution(pwm_query, scores)
#     method_to_scores['pfm_masked_convolution'] = sliding_masked_convolution(pfm_query, scores)
#     method_to_scores['weight_masked_convolution'] = sliding_masked_convolution(weights_query, scores)


    return method_to_scores
method_to_seqlet_scores = collect_various_scoring(pwm_query, seqlet_scores, seqlet_onehot)
#method_to_scores = collect_various_scoring(pwm_query, deeplift_data['scores'], deeplift_data['one_hot'])

#collet all the instances with greater than min_tot_imp deeplift_some_score
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

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

((68563, 15),
 Index(['sum_score', 'seqno', 'start_idx', 'pwm_jaccard', 'pfm_convolution',
        'pwm_scan', 'weight_jaccard', 'pfm_jaccard',
        'weight_masked_convolution', 'weight_scan', 'pfm_masked_convolution',
        'pwm_masked_convolution', 'pwm_convolution', 'weight_convolution',
        'pfm_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.

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

(114461, 15)

from sklearn.cluster import SpectralClustering
# cluster the insances into two groups and select the higher score group
def select_hits(df, clustering=None, cno=None):
    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
    if clustering is None:
        clustering = SpectralClustering(n_clusters=2, assign_labels="discretize", random_state=0).fit(x.values)
    else:
        clustering = clustering.fit_predict(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, (clustering, cno)

filtered_seqlet_hits, (clustering, cno) = select_hits(seqlet_hits)
filtered_seqlet_hits.shape

(7803, 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

(6974, 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)

(7803, 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.6875157222239069
seqno 0
start_idx 0
pwm_scan 1.534107442506396
pfm_scan 5.071605340673137
weight_scan 7.753818386101301
pwm_jaccard 9.963956813823062e-05
pfm_jaccard 0.008519029397890848
weight_jaccard -0.010971767952219543
pwm_convolution -2.7326738123957672
pfm_convolution 0.14043241161255088
weight_convolution 0.1180979073190949
pwm_masked_convolution -0.407800768709883
pfm_masked_convolution 0.1605796616485947
weight_masked_convolution 0.04252452831016403
{5055, 1635, 2437, 4616, 1353, 2185, 6253, 1550, 2446, 2257, 2353, 2325, 2038, 1945, 5019, 2303}
(6974, 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.6875157222239069
seqno 0
start_idx 0
pwm_scan 1.534107442506396
pfm_scan 5.071605340673137
weight_scan 7.814415591311394
pwm_jaccard 0.0006804774557692602
pfm_jaccard 0.03658666368648715
weight_jaccard 0.009186136872677716
pwm_convolution -0.44046356489329536
pfm_convolution 0.30598300167084314
weight_convolution 0.1180979073190949
pwm_masked_convolution -0.2788103773203126
pfm_masked_convolution 0.24785115301053345
weight_masked_convolution 0.04252452831016403
{5055, 1635, 2437, 4616, 1353, 2185, 6253, 1550, 2446, 2257, 2353, 2325, 2038, 1945, 5019, 2303}

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('Seqno', seqno, 'is NOT part of any hits')
        viz_sequence.plot_weights(scores[seqno, start:end,:],
                                  highlight={"black": [(37, 37+len(pwm))]},
                                  subticks_frequency=10, figsize=(20,1))

        viz_sequence.plot_weights(one_hot[seqno, start:end, :], 
                                  highlight={"black": [(37, 37+len(pwm))]},
                                  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=5055, start=0, end=70)
for seqno in np.sort(np.array([5055, 1635, 2437, 4616, 1353, 2185, 6253, 1550, 2446, 2257, 2353, 2325, 2038, 1945, 5019, 2303])):
    view_test_case(seqlet_scores, seqlet_onehot, filtered_seqlet_hits2, pwm_query, seqno=seqno, start=0, end=70)

Seqno 1353 is NOT part of any hits

Seqno 1550 is NOT part of any hits

Seqno 1635 is NOT part of any hits

Seqno 1945 is NOT part of any hits

Seqno 2038 is NOT part of any hits

Seqno 2185 is NOT part of any hits

Seqno 2257 is NOT part of any hits

Seqno 2303 is NOT part of any hits

Seqno 2325 is NOT part of any hits

Seqno 2353 is NOT part of any hits

def view_test_cases(scores, one_hot, df, pwm, start=400, end=600, max_examples_to_view=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))
        if max_examples_to_view is None:
            continue
        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, 
               max_examples_to_view=1000
               )

8 ['pwm_jaccard: 0.063', 'pfm_convolution: 3.353', 'pwm_scan: 9.753', 'sum_score: 4.121', 'weight_jaccard: 0.591', 'pfm_jaccard: 0.423', 'weight_masked_convolution: 0.964', 'weight_scan: 10.593', 'pfm_masked_convolution: 0.915', 'pwm_masked_convolution: 0.853', 'pwm_convolution: 4.512', 'weight_convolution: 5.716', 'pfm_scan: 7.628'] coordinates: (0, 23)
23 23 34

22 ['pwm_jaccard: 0.087', 'pfm_convolution: 6.069', 'pwm_scan: 12.42', 'sum_score: 7.098', 'weight_jaccard: 0.671', 'pfm_jaccard: 0.542', 'weight_masked_convolution: 0.956', 'weight_scan: 11.712', 'pfm_masked_convolution: 0.931', 'pwm_masked_convolution: 0.927', 'pwm_convolution: 8.616', 'weight_convolution: 9.656', 'pfm_scan: 8.703'] coordinates: (0, 37)
37 37 48

55 ['pwm_jaccard: 0.058', 'pfm_convolution: 5.343', 'pwm_scan: 9.621', 'sum_score: 6.51', 'weight_jaccard: 0.481', 'pfm_jaccard: 0.401', 'weight_masked_convolution: 0.889', 'weight_scan: 10.747', 'pfm_masked_convolution: 0.864', 'pwm_masked_convolution: 0.823', 'pwm_convolution: 7.351', 'weight_convolution: 8.925', 'pfm_scan: 7.647'] coordinates: (1, 37)
37 37 48

88 ['pwm_jaccard: 0.085', 'pfm_convolution: 6.392', 'pwm_scan: 11.034', 'sum_score: 7.291', 'weight_jaccard: 0.574', 'pfm_jaccard: 0.464', 'weight_masked_convolution: 0.914', 'weight_scan: 11.421', 'pfm_masked_convolution: 0.908', 'pwm_masked_convolution: 0.918', 'pwm_convolution: 9.264', 'weight_convolution: 10.218', 'pfm_scan: 8.208'] coordinates: (2, 37)
37 37 48

101 ['pwm_jaccard: 0.095', 'pfm_convolution: 1.805', 'pwm_scan: 9.39', 'sum_score: 1.974', 'weight_jaccard: 0.613', 'pfm_jaccard: 0.51', 'weight_masked_convolution: 0.955', 'weight_scan: 10.961', 'pfm_masked_convolution: 0.921', 'pwm_masked_convolution: 0.89', 'pwm_convolution: 2.797', 'weight_convolution: 2.96', 'pfm_scan: 8.011'] coordinates: (2, 50)
50 50 61

117 ['pwm_jaccard: 0.089', 'pfm_convolution: 7.955', 'pwm_scan: 12.42', 'sum_score: 9.376', 'weight_jaccard: 0.61', 'pfm_jaccard: 0.49', 'weight_masked_convolution: 0.912', 'weight_scan: 11.712', 'pfm_masked_convolution: 0.9', 'pwm_masked_convolution: 0.897', 'pwm_convolution: 11.292', 'weight_convolution: 12.482', 'pfm_scan: 8.703'] coordinates: (3, 37)
37 37 48

138 ['pwm_jaccard: 0.089', 'pfm_convolution: 6.631', 'pwm_scan: 12.42', 'sum_score: 7.83', 'weight_jaccard: 0.552', 'pfm_jaccard: 0.532', 'weight_masked_convolution: 0.887', 'weight_scan: 11.712', 'pfm_masked_convolution: 0.879', 'pwm_masked_convolution: 0.883', 'pwm_convolution: 9.479', 'weight_convolution: 10.356', 'pfm_scan: 8.703'] coordinates: (4, 37)
37 37 48

152 ['pwm_jaccard: 0.04', 'pfm_convolution: 1.09', 'pwm_scan: 5.649', 'sum_score: 1.662', 'weight_jaccard: 0.451', 'pfm_jaccard: 0.352', 'weight_masked_convolution: 0.842', 'weight_scan: 9.499', 'pfm_masked_convolution: 0.871', 'pwm_masked_convolution: 0.663', 'pwm_convolution: 1.211', 'weight_convolution: 1.864', 'pfm_scan: 6.158'] coordinates: (4, 51)
51 51 62

169 ['pwm_jaccard: 0.089', 'pfm_convolution: 6.118', 'pwm_scan: 12.42', 'sum_score: 7.175', 'weight_jaccard: 0.674', 'pfm_jaccard: 0.514', 'weight_masked_convolution: 0.947', 'weight_scan: 11.712', 'pfm_masked_convolution: 0.928', 'pwm_masked_convolution: 0.922', 'pwm_convolution: 8.662', 'weight_convolution: 9.678', 'pfm_scan: 8.703'] coordinates: (5, 37)
37 37 48

187 ['pwm_jaccard: 0.084', 'pfm_convolution: 2.257', 'pwm_scan: 9.179', 'sum_score: 2.493', 'weight_jaccard: 0.537', 'pfm_jaccard: 0.415', 'weight_masked_convolution: 0.861', 'weight_scan: 10.903', 'pfm_masked_convolution: 0.879', 'pwm_masked_convolution: 0.887', 'pwm_convolution: 3.335', 'weight_convolution: 3.602', 'pfm_scan: 7.642'] coordinates: (5, 55)
55 55 66