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-10-07 09:58:06,657 [WARNING] git-lfs not installed
TF-MoDISco is using the TensorFlow backend.

from matlas.matches import DenovoModisco
task_idx = 273
modisco_dir = "{0}/task_{1}-naivegw".format(modisco_root, task_idx)
ob = DenovoModisco(modisco_dir)

task_idx = 273
modisco_dir = "{0}/task_{1}-naivegw".format(modisco_root, task_idx)
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}/task_{1}-naivegw/NaiveauPRC.txt".format(model_root, task_idx)
deeplift_hdf = "{0}/{1}/Naive_deeplift2019/task_{2}-naivegw/summit.h5".format(srv_root, data_name, task_idx)
ob = DenovoModisco(modisco_dir)
pattern_name = "metacluster_1/pattern_1"
query = ob.load_chomped_pattern(pattern_name, plot=True)

ob.initialize_with_importance_scores(summit_file, perf_file, deeplift_hdf, query.keys()) #loads all the important scores
deeplift_data = ob.deeplift_data

from matlas.utils import get_calibrated_peaks

calibrated_data = get_calibrated_peaks(summit_file, perf_file)
selected_summits = calibrated_data['selected_summits']
indices = calibrated_data['indices']

from basepair.plot.profiles import extract_signal
pattern_name = 'metacluster_1/pattern_1'
one_hot = ob.seqs_per_pattern[pattern_name]
scores = ob.contrib_per_pattern[pattern_name]
seqlets_ = ob.seqlets_per_pattern[pattern_name]
shuffled_onehot = extract_signal(ob.deeplift_data['shuffled_onehot'][indices], seqlets_)
shuffled_scores = extract_signal(ob.deeplift_data['shuffled_scores'][indices], seqlets_)

one_hot.shape, scores.shape, shuffled_onehot.shape, shuffled_scores.shape

((6372, 70, 4), (6372, 70, 4), (6372, 70, 4), (6372, 70, 4))

from matlas.matches import get_null_sequences
seq_null_onehot, seq_null_scores = get_null_sequences(one_hot, scores,
                                                shuffle_seq=True, shuffle_weight=False
                                                     )

weight_null_onehot, weight_null_scores = get_null_sequences(one_hot, scores,
                                                shuffle_seq=False, shuffle_weight=True
                                                     )

weightedseq_null_onehot, weightedseq_null_scores = get_null_sequences(one_hot, scores,
                                                shuffle_seq=True, shuffle_weight=True
                                                     )

done multiprocessing
(6372, 70, 4)
(6372, 70, 4)
done multiprocessing
(6372, 70, 4)
(6372, 70, 4)
done multiprocessing
(6372, 70, 4)
(6372, 70, 4)

from matlas.matches import collect_scoring_measures
seqlet_scores = ob.contrib_per_pattern[pattern_name]
seqlet_onehot = ob.seqs_per_pattern[pattern_name]
method_to_seqlet_scores = collect_scoring_measures(query, seqlet_scores, seqlet_onehot)
query.keys(), np.sort(np.array(list(method_to_seqlet_scores.keys())))

(dict_keys(['ppm', 'pwm', 'log10pwm', 'pssm', 'cwm', 'hyp_cwm', 'consensus']),
 array(['cwm_convolution', 'cwm_jaccard', 'cwm_masked_convolution',
        'cwm_scan', 'hyp_cwm_convolution', 'hyp_cwm_jaccard',
        'hyp_cwm_masked_convolution', 'hyp_cwm_scan',
        'log10pwm_convolution', 'log10pwm_jaccard',
        'log10pwm_masked_convolution', 'log10pwm_scan', 'ppm_convolution',
        'ppm_jaccard', 'ppm_masked_convolution', 'ppm_scan',
        'pssm_convolution', 'pssm_jaccard', 'pssm_masked_convolution',
        'pssm_scan', 'pwm_convolution', 'pwm_jaccard',
        'pwm_masked_convolution', 'pwm_scan', 'sum_score'], dtype='<U27'))

method_to_shuffled_scores = collect_scoring_measures(query, shuffled_scores, shuffled_onehot)
method_to_seq_null_scores = collect_scoring_measures(query, seq_null_scores, seq_null_onehot)
method_to_weight_null_scores = collect_scoring_measures(query, weight_null_scores, weight_null_onehot)
method_to_weightedseq_null_scores = collect_scoring_measures(query, weightedseq_null_scores, weightedseq_null_onehot)

from matlas.matches import get_any_high_scoring_instances
per_base_imp=0.0625

keys_for_scoring = ['sum_score'] 
keys_for_scoring += ["{}_scan".format(key) for key in query.keys() if key!= 'consensus']
keys_for_scoring += ["{}_jaccard".format(key) for key in query.keys() if key!= 'consensus']
keys_for_scoring += ["{}_convolution".format(key) for key in query.keys() if key!= 'consensus']
keys_for_scoring += ["{}_masked_convolution".format(key) for key in query.keys() if key!= 'consensus']
        
seqlet_hits = get_any_high_scoring_instances(method_to_seqlet_scores,
                                                min_tot_imp=per_base_imp*len(query['ppm']),
                                                keys=keys_for_scoring)

shuffled_seqlet_hits = get_any_high_scoring_instances(method_to_shuffled_scores,
                                                min_tot_imp=per_base_imp*len(query['ppm']),
                                                keys=keys_for_scoring)
seq_null_seqlet_hits = get_any_high_scoring_instances(method_to_seq_null_scores,
                                                min_tot_imp=per_base_imp*len(query['ppm']),
                                                keys=keys_for_scoring)
weight_null_seqlet_hits = get_any_high_scoring_instances(method_to_weight_null_scores,
                                                min_tot_imp=per_base_imp*len(query['ppm']),
                                                keys=keys_for_scoring)
weightedseq_null_seqlet_hits = get_any_high_scoring_instances(method_to_weightedseq_null_scores,
                                                min_tot_imp=per_base_imp*len(query['ppm']),
                                                keys=keys_for_scoring)

seqlet_hits['marker'] = 'foreground'
shuffled_seqlet_hits['marker'] = 'suffled'
seq_null_seqlet_hits['marker'] = 'seq_null'
weight_null_seqlet_hits['marker'] = 'weight_null'
weightedseq_null_seqlet_hits['marker'] = 'weightedseq'
test_hits = pd.concat([seqlet_hits,
                       shuffled_seqlet_hits,
                       seq_null_seqlet_hits,
                       weight_null_seqlet_hits,
                       weightedseq_null_seqlet_hits
                      ])

from sklearn.preprocessing import StandardScaler
from sklearn.decomposition import PCA

def transform_data(_hits, scaler=None, pca=None):
    if scaler is None:
        # should be used only during training
        scaler = StandardScaler()
        scaler.fit(_hits)
        
    # Apply transform to both the training set and testing set
    _hits = scaler.transform(_hits)

    if pca is None:
        # Make an instance of the Model
        pca = PCA(.98)
        #fitting PCA on the training set only.
        pca.fit(_hits)
        
    #Apply the mapping (transform) to both the training set and the test set.
    _hits = pca.transform(_hits)
    
    return scaler, pca, _hits

scaler, pca, train_hits = transform_data(seqlet_hits[keys_for_scoring])
pca.explained_variance_ratio_

array([0.69987347, 0.09950445, 0.08237547, 0.04125603, 0.02750585,
       0.01404074, 0.01112795, 0.00870734])

#read the real truth
def load_true_positives(tp_file):
    coordinates = []
    with open(tp_file, "r") as fp:
        for line in fp:
            line = line.strip().split(",")
            seqno = line[0]
            pos = int(line[1])
            if "-" in seqno:
                seqno = seqno.split("-")
                for i in range(int(seqno[0]), int(seqno[1])+1):
                    coordinates.append((i, pos))
            else:
                coordinates.append((int(seqno), pos))
    return coordinates

def get_tp_indices(seqlet_hits, verbose=True):
    keep = []
    for cord in coordinates:
        seqno, pos = cord
        df = seqlet_hits[(seqlet_hits['seqno']==seqno) & (seqlet_hits['start_idx']==pos)]
        if df.shape[0]>0:
            keep.append(df.index[0])
        else:
            if verbose: print(cord)

    return keep

tp_file = "{}/web_mouse_heme_data/modisco_tracks/binary_model/task_{}-naivegw/mm11_tp_spectral.txt".format(project_root, task_idx)
coordinates = load_true_positives(tp_file)
true_indices = get_tp_indices(seqlet_hits, False)

import matplotlib.pyplot as plt
def plot_pca(df, true_ids=None):
    print(df.shape)
    fig = plt.figure(figsize = (8,8))
    ax = fig.add_subplot(1,1,1) 
    ax.set_xlabel('PC1', fontsize = 15)
    ax.set_ylabel('PC2', fontsize = 15)
    ax.set_title('2 component PCA', fontsize = 20)
    if true_ids is not None:
        true_df = df.loc[np.array(true_ids),:]
        df = df.loc[np.array(list(set(list(df.index)).difference(set(true_ids)))),:]
        ax.scatter(true_df['PC1'], true_df['PC2'], s = 1, c='r', marker='+')
    ax.scatter(df['PC1'], df['PC2'], s = 1, c='b')
    #ax.axvline(9)
    #ax.axhline(-1)
    ax.grid()
    return None
pc_df = pd.DataFrame(data = train_hits, columns=["PC{}".format(i+1) for i in range(train_hits.shape[1])])
plot_pca(pc_df, true_indices)

(114461, 8)

len(true_indices)

7312

from matlas.matches import view_test_case
from modisco.visualization import viz_sequence
#from matlas.matches import view_test_cases

def view_test_cases(scores, one_hot, df, pwm, start=400, end=600, examples_to_view=None):
    
    seqnos = df.seqno.values
    positions = df.start_idx.values
    if examples_to_view is not None:
        seqnos = seqnos[examples_to_view]
        positions = positions[examples_to_view]
    
    for seqno, pos in zip(seqnos, positions):
        view_test_case(scores, one_hot, df, pwm, seqno, start=start, end=end, itemno=pos)
        
    return None
tmp_df = seqlet_hits.loc[np.sort((np.array(true_indices))),:]
# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=tmp_df.start_idx.values!=37)

# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=tmp_df.seqno.values<=1000 
#                )

# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=np.logical_and(tmp_df.seqno.values>1000, tmp_df.seqno.values<=2000) 
#                )

# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=np.logical_and(tmp_df.seqno.values>2000, tmp_df.seqno.values<=3000) 
#                )

view_test_cases(seqlet_scores, 
               seqlet_onehot, 
               tmp_df, 
               query['pwm'], 
               start=0, end=70, 
               examples_to_view=np.logical_and(tmp_df.seqno.values>3000, tmp_df.seqno.values<=4000) 
               )

56051 ['sum_score: 2.305', 'hyp_cwm_convolution: 0.596', 'pwm_masked_convolution: 0.839', 'pwm_jaccard: 0.076', 'pwm_scan: 11.691', 'log10pwm_convolution: 3.044', 'cwm_convolution: 0.85', 'ppm_jaccard: 0.379', 'ppm_convolution: 2.102', 'cwm_jaccard: 0.516', 'hyp_cwm_scan: 1.769', 'cwm_scan: 2.482', 'pssm_convolution: 3.517', 'ppm_masked_convolution: 0.853', 'log10pwm_jaccard: 0.076', 'log10pwm_scan: 8.104', 'pssm_jaccard: 0.478', 'cwm_masked_convolution: 0.888', 'hyp_cwm_jaccard: 0.249', 'hyp_cwm_masked_convolution: 0.887', 'pssm_scan: 10.25', 'ppm_scan: 7.323', 'pssm_masked_convolution: 0.881', 'pwm_convolution: 4.391', 'log10pwm_masked_convolution: 0.839'] coordinates: (3001, 37)
37 37 48

56074 ['sum_score: 2.731', 'hyp_cwm_convolution: 0.622', 'pwm_masked_convolution: 0.931', 'pwm_jaccard: 0.061', 'pwm_scan: 12.492', 'log10pwm_convolution: 3.112', 'cwm_convolution: 0.912', 'ppm_jaccard: 0.487', 'ppm_convolution: 2.325', 'cwm_jaccard: 0.688', 'hyp_cwm_scan: 1.718', 'cwm_scan: 2.477', 'pssm_convolution: 3.801', 'ppm_masked_convolution: 0.967', 'log10pwm_jaccard: 0.061', 'log10pwm_scan: 8.659', 'pssm_jaccard: 0.616', 'cwm_masked_convolution: 0.98', 'hyp_cwm_jaccard: 0.218', 'hyp_cwm_masked_convolution: 0.967', 'pssm_scan: 10.576', 'ppm_scan: 7.434', 'pssm_masked_convolution: 0.974', 'pwm_convolution: 4.489', 'log10pwm_masked_convolution: 0.931'] coordinates: (3002, 37)
37 37 48

56104 ['sum_score: 3.358', 'hyp_cwm_convolution: 0.725', 'pwm_masked_convolution: 0.885', 'pwm_jaccard: 0.071', 'pwm_scan: 14.492', 'log10pwm_convolution: 3.607', 'cwm_convolution: 1.048', 'ppm_jaccard: 0.495', 'ppm_convolution: 2.744', 'cwm_jaccard: 0.684', 'hyp_cwm_scan: 1.891', 'cwm_scan: 2.615', 'pssm_convolution: 4.454', 'ppm_masked_convolution: 0.928', 'log10pwm_jaccard: 0.071', 'log10pwm_scan: 10.045', 'pssm_jaccard: 0.659', 'cwm_masked_convolution: 0.938', 'hyp_cwm_jaccard: 0.242', 'hyp_cwm_masked_convolution: 0.933', 'pssm_scan: 10.869', 'ppm_scan: 7.949', 'pssm_masked_convolution: 0.955', 'pwm_convolution: 5.203', 'log10pwm_masked_convolution: 0.885'] coordinates: (3003, 37)
37 37 48

56119 ['sum_score: 3.196', 'hyp_cwm_convolution: 0.619', 'pwm_masked_convolution: 0.825', 'pwm_jaccard: 0.062', 'pwm_scan: 12.599', 'log10pwm_convolution: 3.261', 'cwm_convolution: 0.888', 'ppm_jaccard: 0.504', 'ppm_convolution: 2.507', 'cwm_jaccard: 0.628', 'hyp_cwm_scan: 1.628', 'cwm_scan: 2.335', 'pssm_convolution: 3.736', 'ppm_masked_convolution: 0.947', 'log10pwm_jaccard: 0.062', 'log10pwm_scan: 8.733', 'pssm_jaccard: 0.54', 'cwm_masked_convolution: 0.944', 'hyp_cwm_jaccard: 0.208', 'hyp_cwm_masked_convolution: 0.926', 'pssm_scan: 9.774', 'ppm_scan: 7.597', 'pssm_masked_convolution: 0.933', 'pwm_convolution: 4.704', 'log10pwm_masked_convolution: 0.825'] coordinates: (3004, 37)
37 37 48

56143 ['sum_score: 0.702', 'hyp_cwm_convolution: 0.095', 'pwm_masked_convolution: 0.183', 'pwm_jaccard: 0.036', 'pwm_scan: 0.663', 'log10pwm_convolution: 0.225', 'cwm_convolution: 0.141', 'ppm_jaccard: 0.296', 'ppm_convolution: 0.405', 'cwm_jaccard: 0.335', 'hyp_cwm_scan: 1.248', 'cwm_scan: 1.884', 'pssm_convolution: 0.59', 'ppm_masked_convolution: 0.776', 'log10pwm_jaccard: 0.036', 'log10pwm_scan: 0.459', 'pssm_jaccard: 0.308', 'cwm_masked_convolution: 0.716', 'hyp_cwm_jaccard: 0.179', 'hyp_cwm_masked_convolution: 0.641', 'pssm_scan: 7.816', 'ppm_scan: 5.705', 'pssm_masked_convolution: 0.711', 'pwm_convolution: 0.324', 'log10pwm_masked_convolution: 0.183'] coordinates: (3005, 18)
18 18 29

56144 ['sum_score: 0.7', 'hyp_cwm_convolution: 0.122', 'pwm_masked_convolution: 0.529', 'pwm_jaccard: 0.056', 'pwm_scan: 3.837', 'log10pwm_convolution: 0.617', 'cwm_convolution: 0.193', 'ppm_jaccard: 0.392', 'ppm_convolution: 0.524', 'cwm_jaccard: 0.454', 'hyp_cwm_scan: 1.236', 'cwm_scan: 2.229', 'pssm_convolution: 0.827', 'ppm_masked_convolution: 0.917', 'log10pwm_jaccard: 0.056', 'log10pwm_scan: 2.66', 'pssm_jaccard: 0.457', 'cwm_masked_convolution: 0.853', 'hyp_cwm_jaccard: 0.18', 'hyp_cwm_masked_convolution: 0.804', 'pssm_scan: 9.581', 'ppm_scan: 6.434', 'pssm_masked_convolution: 0.861', 'pwm_convolution: 0.891', 'log10pwm_masked_convolution: 0.529'] coordinates: (3005, 19)
19 19 30

56153 ['sum_score: 1.577', 'hyp_cwm_convolution: 0.378', 'pwm_masked_convolution: 0.785', 'pwm_jaccard: 0.064', 'pwm_scan: 11.34', 'log10pwm_convolution: 1.872', 'cwm_convolution: 0.531', 'ppm_jaccard: 0.35', 'ppm_convolution: 1.377', 'cwm_jaccard: 0.522', 'hyp_cwm_scan: 1.451', 'cwm_scan: 2.218', 'pssm_convolution: 2.177', 'ppm_masked_convolution: 0.877', 'log10pwm_jaccard: 0.064', 'log10pwm_scan: 7.86', 'pssm_jaccard: 0.441', 'cwm_masked_convolution: 0.961', 'hyp_cwm_jaccard: 0.215', 'hyp_cwm_masked_convolution: 0.97', 'pssm_scan: 9.77', 'ppm_scan: 7.453', 'pssm_masked_convolution: 0.909', 'pwm_convolution: 2.701', 'log10pwm_masked_convolution: 0.785'] coordinates: (3005, 37)
37 37 48

56163 ['sum_score: 1.126', 'hyp_cwm_convolution: 0.231', 'pwm_masked_convolution: 0.795', 'pwm_jaccard: 0.061', 'pwm_scan: 7.967', 'log10pwm_convolution: 1.238', 'cwm_convolution: 0.341', 'ppm_jaccard: 0.349', 'ppm_convolution: 0.94', 'cwm_jaccard: 0.448', 'hyp_cwm_scan: 1.313', 'cwm_scan: 1.995', 'pssm_convolution: 1.475', 'ppm_masked_convolution: 0.933', 'log10pwm_jaccard: 0.061', 'log10pwm_scan: 5.523', 'pssm_jaccard: 0.44', 'cwm_masked_convolution: 0.913', 'hyp_cwm_jaccard: 0.173', 'hyp_cwm_masked_convolution: 0.918', 'pssm_scan: 8.915', 'ppm_scan: 6.439', 'pssm_masked_convolution: 0.906', 'pwm_convolution: 1.786', 'log10pwm_masked_convolution: 0.795'] coordinates: (3006, 23)
23 23 34

56177 ['sum_score: 2.293', 'hyp_cwm_convolution: 0.552', 'pwm_masked_convolution: 0.883', 'pwm_jaccard: 0.066', 'pwm_scan: 12.061', 'log10pwm_convolution: 2.645', 'cwm_convolution: 0.773', 'ppm_jaccard: 0.433', 'ppm_convolution: 1.964', 'cwm_jaccard: 0.63', 'hyp_cwm_scan: 1.775', 'cwm_scan: 2.483', 'pssm_convolution: 3.173', 'ppm_masked_convolution: 0.925', 'log10pwm_jaccard: 0.066', 'log10pwm_scan: 8.36', 'pssm_jaccard: 0.574', 'cwm_masked_convolution: 0.941', 'hyp_cwm_jaccard: 0.237', 'hyp_cwm_masked_convolution: 0.972', 'pssm_scan: 10.534', 'ppm_scan: 7.379', 'pssm_masked_convolution: 0.92', 'pwm_convolution: 3.815', 'log10pwm_masked_convolution: 0.883'] coordinates: (3006, 37)
37 37 48

56205 ['sum_score: 2.019', 'hyp_cwm_convolution: 0.406', 'pwm_masked_convolution: 0.785', 'pwm_jaccard: 0.047', 'pwm_scan: 8.981', 'log10pwm_convolution: 1.922', 'cwm_convolution: 0.576', 'ppm_jaccard: 0.41', 'ppm_convolution: 1.548', 'cwm_jaccard: 0.543', 'hyp_cwm_scan: 1.478', 'cwm_scan: 2.117', 'pssm_convolution: 2.434', 'ppm_masked_convolution: 0.951', 'log10pwm_jaccard: 0.047', 'log10pwm_scan: 6.225', 'pssm_jaccard: 0.507', 'cwm_masked_convolution: 0.942', 'hyp_cwm_jaccard: 0.191', 'hyp_cwm_masked_convolution: 0.93', 'pssm_scan: 9.053', 'ppm_scan: 6.624', 'pssm_masked_convolution: 0.933', 'pwm_convolution: 2.773', 'log10pwm_masked_convolution: 0.785'] coordinates: (3007, 37)
37 37 48

# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=np.logical_and(tmp_df.seqno.values>4000, tmp_df.seqno.values<=5000) 
#                )

# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=np.logical_and(tmp_df.seqno.values>5000, tmp_df.seqno.values<=6000) 
#                )

# view_test_cases(seqlet_scores, 
#                seqlet_onehot, 
#                tmp_df, 
#                query['pwm'], 
#                start=0, end=70, 
#                examples_to_view=tmp_df.seq.values>6000 
#                )

from sklearn.mixture import GaussianMixture
#import sklearn.mixture.GaussianMixture
gmm = GaussianMixture(n_components=2).fit(pc_df[['PC1', 'PC2']])
labels = gmm.predict(pc_df[['PC1', 'PC2']])
plt.scatter(pc_df['PC1'], pc_df['PC2'], c=labels, s=1, cmap='viridis');

probs = gmm.predict_proba(pc_df[['PC1', 'PC2']])
plt.scatter(pc_df['PC1'], pc_df['PC2'], c=probs[:,1], s=1, cmap='seismic');

_, _, test_hits = transform_data(shuffled_seqlet_hits[keys_for_scoring], scaler, pca)
shuffled_pc_df = pd.DataFrame(data = test_hits, columns=["PC{}".format(i+1) for i in range(test_hits.shape[1])])
#plot_pca(shuffled_pc_df, true_ids=None)

probs = gmm.predict_proba(shuffled_pc_df[['PC1', 'PC2']])
plt.scatter(shuffled_pc_df['PC1'], shuffled_pc_df['PC2'], c=probs[:,1], s=1, cmap='seismic');

_, _, test_hits = transform_data(seq_null_seqlet_hits[keys_for_scoring], scaler, pca)
seq_null_pc_df = pd.DataFrame(data = test_hits, columns=["PC{}".format(i+1) for i in range(test_hits.shape[1])])
#plot_pca(shuffled_pc_df, true_ids=None)

probs = gmm.predict_proba(seq_null_pc_df[['PC1', 'PC2']])
plt.scatter(seq_null_pc_df['PC1'], seq_null_pc_df['PC2'], c=probs[:,1], s=1, cmap='seismic');

_, _, test_hits = transform_data(weight_null_seqlet_hits[keys_for_scoring], scaler, pca)
weight_null_pc_df = pd.DataFrame(data = test_hits, columns=["PC{}".format(i+1) for i in range(test_hits.shape[1])])
#plot_pca(weight_null_pc_df, true_ids=None)

probs = gmm.predict_proba(weight_null_pc_df[['PC1', 'PC2']])
plt.scatter(weight_null_pc_df['PC1'], weight_null_pc_df['PC2'], c=probs[:,1], s=1, cmap='seismic');

_, _, test_hits = transform_data(weightedseq_null_seqlet_hits[keys_for_scoring], scaler, pca)
weightedseq_null_pc_df = pd.DataFrame(data = test_hits, columns=["PC{}".format(i+1) for i in range(test_hits.shape[1])])
#plot_pca(weightedseq_null_pc_df, true_ids=None)

probs = gmm.predict_proba(weightedseq_null_pc_df[['PC1', 'PC2']])
plt.scatter(weightedseq_null_pc_df['PC1'], weightedseq_null_pc_df['PC2'], c=probs[:,1], s=1, cmap='seismic');

pca_seqlet_hits = pc_df[['PC1', 'PC2', 'PC3']]
pca_shuffled_seqlet_hits = shuffled_pc_df[['PC1', 'PC2', 'PC3']]
pca_seq_null_seqlet_hits = seq_null_pc_df[['PC1', 'PC2', 'PC3']]
pca_weight_null_seqlet_hits = weight_null_pc_df[['PC1', 'PC2', 'PC3']]
pca_weightedseq_null_seqlet_hits = weightedseq_null_pc_df[['PC1', 'PC2', 'PC3']]
pca_seqlet_hits['marker'] = 'foreground'
pca_shuffled_seqlet_hits['marker'] = 'suffled'
pca_seq_null_seqlet_hits['marker'] = 'seq_null'
pca_weight_null_seqlet_hits['marker'] = 'weight_null'
pca_weightedseq_null_seqlet_hits['marker'] = 'weightedseq'
pca_test_hits = pd.concat([pca_seqlet_hits,
                       pca_shuffled_seqlet_hits,
                       pca_seq_null_seqlet_hits,
                       pca_weight_null_seqlet_hits,
                       pca_weightedseq_null_seqlet_hits
                      ])

# import seaborn as sns
# sns.pairplot(pca_test_hits, 
#              diag_kind="kde", kind='scatter', hue="marker")

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

<seaborn.axisgrid.PairGrid at 0x7f0669a8beb8>

# sns.pairplot(pd.concat([pca_seqlet_hits, pca_shuffled_seqlet_hits]), 
#              diag_kind="kde", kind='scatter', hue="marker")

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

<seaborn.axisgrid.PairGrid at 0x7f0683430780>

# sns.pairplot(pd.concat([pca_seqlet_hits, pca_seq_null_seqlet_hits]), 
#              diag_kind="kde", kind='scatter', hue="marker")

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

<seaborn.axisgrid.PairGrid at 0x7f013adc8278>

# sns.pairplot(pd.concat([pca_seqlet_hits, pca_weight_null_seqlet_hits]), 
#              diag_kind="kde", kind='scatter', hue="marker")

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

<seaborn.axisgrid.PairGrid at 0x7f06833e0978>

# sns.pairplot(pd.concat([pca_seqlet_hits, pca_weightedseq_null_seqlet_hits]), 
#              diag_kind="kde", kind='scatter', hue="marker")

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

<seaborn.axisgrid.PairGrid at 0x7f0680a80cf8>

from sklearn.cluster import SpectralClustering
from sklearn.svm import LinearSVC, SVC
from sklearn.cluster import AffinityPropagation
def get_spectral_high_scoring_instances(train_df, keys_to_cluster, test_df=None, clustering=None, cno=None):
    train = train_df[keys_to_cluster]
    
    if clustering is None:
        clustering = SpectralClustering(n_clusters=2, assign_labels="discretize", random_state=0).fit(train.values)
        mean0 = np.mean(train_df[keys_to_cluster[0]][clustering.labels_==0])
        mean1 = np.mean(train_df[keys_to_cluster[0]][clustering.labels_==1])
        cno = 1
        if(mean0 > mean1):
            cno = 0
        df = train_df[clustering.labels_==cno]
    else:
        test = test_df[keys_to_cluster]
        y_tr = clustering.labels_
        #clf = LinearSVC(dual=False).fit(train.values, y_tr)
        clf = SVC().fit(train.values, y_tr)
        batch_size = 1000
        batch_count = int(np.ceil(1.0*test_df.shape[0]/batch_size))
        predicted_labels = clf.predict(test.values)
        predicted_labels = np.array(predicted_labels)
        assert len(predicted_labels)==test_df.shape[0], "not all the values have been predicted"
        df = test_df.loc[predicted_labels==cno, :]
        
    #df = df[df['pwm_jaccard']>0]
    
    return df, (clustering, cno)

# import sys
# from matlas.matches import clustering_method_to_functions#, get_spectral_high_scoring_instances
# from random import sample

# max_instances_to_cluster = min(seqlet_hits.shape[0], 20000) #similar size as the gata motif
# pca_seqlet_hits = pc_df[['PC1', 'PC2']]
# train_hits_forcluster = pca_seqlet_hits.loc[sample(list(seqlet_hits.index), max_instances_to_cluster),:]
# #train_hits_forcluster = pc_df[['PC1', 'PC2']]

# clustering_method = 'spectral'
# thismodule = sys.modules[__name__]
# method_to_call = getattr(thismodule, clustering_method_to_functions[clustering_method])

# #use a subset of the hits to build the cluster
# clustered_seqlet_hits, (clustering, cno) = method_to_call(train_hits_forcluster,
#                                                           keys_to_cluster=['PC1', 'PC2']
#                                                          )

# #use the cluster to predict on the full set of seqlets
# clustered_seqlet_hits2, (_, _) = method_to_call(train_hits_forcluster, #train_df
#                                         ['PC1', 'PC2'], #keys_to_cluster
#                                         pca_seqlet_hits, #test_df
#                                         clustering, #clustering
#                                         cno)

# clustered_seqlet_hits.shape, clustered_seqlet_hits2.shape, clustered_seqlet_hits3.shape

((1275, 2), (7307, 4), (605, 4))

# clustered_seqlet_hits3, (_, _) = method_to_call(train_hits_forcluster, #train_df
#                                         ['PC1', 'PC2'], #keys_to_cluster
#                                         pca_shuffled_seqlet_hits, #test_df
#                                         clustering, #clustering
#                                         cno)

# tp_file = "{}/web_mouse_heme_data/modisco_tracks/binary_model/task_{}-naivegw/mm11_tp_random.txt".format(project_root, task_idx)
# write_tp(shuffled_seqlet_hits.loc[clustered_seqlet_hits3.index,:], tp_file)
# coordinates = load_true_positives(tp_file)
# true_indices = get_tp_indices(shuffled_seqlet_hits, verbose=False)
# plot_pca(shuffled_pc_df, true_indices)

(48842, 8)

# clustered_seqlet_hits3, (_, _) = method_to_call(train_hits_forcluster, #train_df
#                                         ['PC1', 'PC2'], #keys_to_cluster
#                                         pca_seq_null_seqlet_hits, #test_df
#                                         clustering, #clustering
#                                         cno)

# tp_file = "{}/web_mouse_heme_data/modisco_tracks/binary_model/task_{}-naivegw/mm11_tp_random.txt".format(project_root, task_idx)
# write_tp(seqlet_hits.loc[clustered_seqlet_hits3.index,:], tp_file)
# coordinates = load_true_positives(tp_file)
# true_indices = get_tp_indices(seq_null_seqlet_hits, verbose=False)
# plot_pca(seq_null_pc_df, true_indices)

(114461, 8)

# clustered_seqlet_hits3, (_, _) = method_to_call(train_hits_forcluster, #train_df
#                                         ['PC1', 'PC2'], #keys_to_cluster
#                                         pca_weight_null_seqlet_hits, #test_df
#                                         clustering, #clustering
#                                         cno)

# tp_file = "{}/web_mouse_heme_data/modisco_tracks/binary_model/task_{}-naivegw/mm11_tp_random.txt".format(project_root, task_idx)
# write_tp(weight_null_seqlet_hits.loc[clustered_seqlet_hits3.index,:], tp_file)
# coordinates = load_true_positives(tp_file)
# true_indices = get_tp_indices(weight_null_seqlet_hits, verbose=False)
# plot_pca(weight_null_pc_df, true_indices)

(124266, 8)

# clustered_seqlet_hits3, (_, _) = method_to_call(train_hits_forcluster, #train_df
#                                         ['PC1', 'PC2'], #keys_to_cluster
#                                         pca_weightedseq_null_seqlet_hits, #test_df
#                                         clustering, #clustering
#                                         cno)

# tp_file = "{}/web_mouse_heme_data/modisco_tracks/binary_model/task_{}-naivegw/mm11_tp_random.txt".format(project_root, task_idx)
# write_tp(weightedseq_null_seqlet_hits.loc[clustered_seqlet_hits3.index,:], tp_file)
# coordinates = load_true_positives(tp_file)
# true_indices = get_tp_indices(weightedseq_null_seqlet_hits, verbose=False)
# plot_pca(weightedseq_null_pc_df, true_indices)

(123596, 8)