<jupyter_start><jupyter_text>Classification Attention Analysis<jupyter_code>%load_ext autoreload
%autoreload 2

from __future__ import print_function
import keras
from keras import backend as K
import tensorflow as tf

import sys
import os
import numpy as np
import pandas as pd
import scipy

import Bio
from Bio import SeqIO
import Bio.motifs

%matplotlib inline
from matplotlib import pyplot as plt
import seaborn as sns
from sklearn import model_selection
import sklearn

from IPython.display import SVG
from keras.utils.vis_utils import model_to_dot

from signal_models import *
from utility_functions import *

sys.path.insert(0, '/home/jtao/code/tba/motif_tools')
from motif_utilities import *

sns.set_context('notebook')
working_directory = '/home/jtao/analysis/genomic_grammar_analysis/'<jupyter_output><empty_output><jupyter_text>Read Data<jupyter_code>treatment_positive_fasta_dict = {}
treatment_negative_fasta_dict = {}
treatment_positive_seqArray_dict = {}
treatment_positive_rc_seqArray_dict = {}
treatment_negative_seqArray_dict = {}
treatment_negative_rc_seqArray_dict = {}

for ps in ['c57bl6_kla-1h_peaks.fasta', 'c57bl6_veh_peaks.fasta', 'c57bl6_il4-24h_peaks.fasta']:
    treatment = ps.split('_')[1]
    print('***',treatment, '***')
    positive_seqRecords = list(SeqIO.parse(working_directory + '/peak_sequences/' + ps, 'fasta'))
    negative_seqRecords = list(SeqIO.parse(working_directory + '/background_files/' + ps.replace('_peaks', '_background'), 'fasta'))[:len(positive_seqRecords)]
    
    positive_fasta_seq = [str(x.seq[:200]) for x in positive_seqRecords]
    negative_fasta_seq = [str(x.seq[:200]) for x in negative_seqRecords]
    
    positive_rc_fasta_seq = [str(x[:200].reverse_complement().seq) for x in positive_seqRecords]
    negative_rc_fasta_seq = [str(x[:200].reverse_complement().seq) for x in negative_seqRecords]
    
    positive_sequence_arrays = convert_sequences_to_array(positive_fasta_seq)
    positive_rc_sequence_arrays = convert_sequences_to_array(positive_rc_fasta_seq)
    
    negative_sequence_arrays = convert_sequences_to_array(negative_fasta_seq)
    negative_rc_sequence_arrays = convert_sequences_to_array(negative_rc_fasta_seq)
    
    treatment_positive_fasta_dict[treatment] = positive_fasta_seq
    treatment_negative_fasta_dict[treatment] = negative_fasta_seq
    treatment_positive_seqArray_dict[treatment] = positive_sequence_arrays
    treatment_positive_rc_seqArray_dict[treatment] = positive_rc_sequence_arrays
    treatment_negative_seqArray_dict[treatment] = negative_sequence_arrays
    treatment_negative_rc_seqArray_dict[treatment] = negative_rc_sequence_arrays<jupyter_output>*** kla-1h ***
*** veh ***
*** il4-24h ***<jupyter_text>Set Model Settings<jupyter_code>### settings for attention model ###
attention_num_motifs=150
attention_motif_size=16
adjacent_bp_pool_size=1
attention_num_dense_neurons = 1
attention_dropout_rate=0.1<jupyter_output><empty_output><jupyter_text>Read in Models<jupyter_code>condition_model_dict = {}
for f in os.listdir(working_directory + '/models_classification/'):
    treatment_modelType = f.replace('_weights.h5','')
    modelType = treatment_modelType.split('_')[0]
    if modelType == 'DotAttention':
        current_model = get_dotProductAttention_model(200,
            mode='classification',
            num_motifs=attention_num_motifs, 
            motif_size=attention_motif_size,
            adjacent_bp_pool_size=adjacent_bp_pool_size,
            num_dense_neurons = attention_num_dense_neurons,
            dropout_rate=attention_dropout_rate,
            num_classes=2)
        current_model.load_weights(working_directory + '/models_classification/' + treatment_modelType + '_weights.h5')
        condition_model_dict[treatment_modelType] = current_model
        print(f, current_model.count_params())<jupyter_output>DotAttention_kla-1h_weights.h5 10753
DotAttention_veh_weights.h5 10753
DotAttention_il4-24h_weights.h5 10753<jupyter_text>Get correctly classified sequences<jupyter_code>treatment_truePositive_seqArray_dict = {}
for treatment in sorted(treatment_positive_seqArray_dict.keys()):
    positive_sequence_arrays = treatment_positive_seqArray_dict[treatment] 
    
    current_model = condition_model_dict['DotAttention_' + treatment]
    
    probs = current_model.predict([positive_sequence_arrays])
    predictions = np.array([0 if x[0] > x[1] else 1 for x in probs])
    truePositve_sequence_arrays = positive_sequence_arrays[predictions == 1]
    treatment_truePositive_seqArray_dict[treatment] = truePositve_sequence_arrays

    print(treatment, len(truePositve_sequence_arrays)/len(positive_sequence_arrays))<jupyter_output>il4-24h 0.9252640372337252
kla-1h 0.9226212343864806
veh 0.8527099841521395<jupyter_text>Weight Extraction<jupyter_code>def get_sequence_attention(model, sequence_arrays, sequence_size):
    get_attention = K.function([model.get_layer('input_fwd').input, 
                                K.learning_phase()
                               ], 
                               [model.get_layer('attention_softmax_layer').output])
    trimmed_sequences = [x[:sequence_size] for x in sequence_arrays]
    
    attention = get_attention((trimmed_sequences, 0))[0]

    return attention
treatment_attention_dict = {}
for treatment in sorted(treatment_truePositive_seqArray_dict.keys()):
    print('***',treatment, '***')

    sequence_arrays = treatment_truePositive_seqArray_dict[treatment]
    if 'DotAttention_'+treatment in condition_model_dict:
        current_model = condition_model_dict['DotAttention_'+treatment]

        attention = get_sequence_attention(current_model,
            sequence_arrays,
            200)

        treatment_attention_dict[treatment] = attention<jupyter_output>*** il4-24h ***
*** kla-1h ***
*** veh ***<jupyter_text>Visualize Attention<jupyter_code>attention = treatment_attention_dict['kla-1h']<jupyter_output><empty_output><jupyter_text>Overall Matrix<jupyter_code>index = 1811
plt.figure(figsize=(15,15))
hm = sns.heatmap(attention[index], 
            yticklabels = False, 
            xticklabels = False, 
            square = True,
           vmin=0,
#            vmax=0.4
           )
cb = hm.collections[0].colorbar
cb.ax.set_visible(False)
index = 19999
# index = 19996
plt.figure(figsize=(15,15))
hm = sns.heatmap(attention[index], 
            yticklabels = False, 
            xticklabels = False, 
            square = True,
#            vmin=0,
#            vmax=0.1
           )
cb = hm.collections[0].colorbar
cb.ax.set_visible(False)
for treatment in treatment_attention_dict:
    treatment_attention = treatment_attention_dict[treatment]
    sns.distplot(treatment_attention[:100].flatten(), label = treatment)
plt.xlabel('Attention')
plt.ylabel('Frequency')
plt.legend()<jupyter_output><empty_output><jupyter_text>Per sequence frequencies<jupyter_code>num_positions = 200
matrix_size = num_positions**2
count_list = []
treatment_list = []
threshold_list = []
for treatment in treatment_attention_dict:
    treatment_attention = treatment_attention_dict[treatment]
    print(treatment)
    
    for threshold in [np.round(x,2) for x in np.arange(0.1,1,0.1)]:
        counter = 0
        print(threshold)
        for am in treatment_attention:
            count = np.sum(am.flatten() > threshold)
            count_list.append(count)
#             normalized_count_list.append(count/matrix_size)
            treatment_list.append(treatment)
            threshold_list.append(threshold)
#             counter+=1
#             if counter % 10000 == 0:
#                 print(counter)
frame = pd.DataFrame({'Treatment': treatment_list, 
                      'Threshold': threshold_list,
                      'Count':count_list})
sns.factorplot(data = frame, 
    x = 'Treatment',
    y = 'Count',
    hue='Threshold',
    kind = 'box',
    palette = 'Blues',
    showfliers=False,
    size = 10)<jupyter_output><empty_output><jupyter_text>Max per Position<jupyter_code>for treatment in treatment_attention_dict:
    treatment_attention = treatment_attention_dict[treatment]
    max_attention = np.array([x.max(axis=0) for x in treatment_attention[:100]])
    sns.distplot(max_attention.flatten(), label = treatment)
plt.xlabel('Max Attention Per Position')
plt.ylabel('Frequency')
plt.legend()<jupyter_output><empty_output><jupyter_text>Per sequence frequencies<jupyter_code>count_list = []
treatment_list = []
threshold_list = []
normalized_count_list = []
for treatment in treatment_attention_dict:
    treatment_attention = treatment_attention_dict[treatment]
    max_attention = np.array([x.max(axis=0) for x in treatment_attention])
    print(treatment)
    for threshold in [np.round(x,2) for x in np.arange(0.1,1,0.1)]:
        for fam in max_attention:
            count = np.sum(fam.flatten() >= threshold)
            count_list.append(count)
            treatment_list.append(treatment)
            threshold_list.append(threshold)
        

frame = pd.DataFrame({'Treatment': treatment_list, 
                      'Threshold': threshold_list,
                      'Count':count_list})
sns.factorplot(data = frame, 
    x = 'Treatment',
    y = 'Count',
    hue='Threshold',
    kind = 'box',
    palette = 'Greens',
    showfliers=False,
    size = 10)<jupyter_output><empty_output><jupyter_text>Motif Extraction<jupyter_code>def get_sequence_motif_scores(model, sequence_arrays, sequence_size):
    get_convolution_kernels = K.function([model.get_layer('input_fwd').input, 
                                K.learning_phase()
                               ], 
                               [model.get_layer('convolution_layer').output])
    trimmed_sequences = [x[:sequence_size] for x in sequence_arrays]
    
    motif_scores = np.array(get_convolution_kernels((trimmed_sequences, 0)))

    return np.squeeze(motif_scores)
treatment_motifScores_dict = {}
for treatment in sorted(treatment_truePositive_seqArray_dict.keys()):
    print('***',treatment, '***')

    sequence_arrays = treatment_truePositive_seqArray_dict[treatment]
    if 'DotAttention_'+treatment in condition_model_dict:
        current_model = condition_model_dict['DotAttention_'+treatment]

        scores = get_sequence_motif_scores(current_model,
            sequence_arrays,
            200)

        treatment_motifScores_dict[treatment] = scores<jupyter_output>*** il4-24h ***
*** kla-1h ***
*** veh ***<jupyter_text>Compare Convolution kernels to known motifs<jupyter_code>def normalize_motif(motif):
    normed_motif = np.zeros(motif.shape)
    mins = motif.min(axis=1)
    for i in range(4):
        normed_motif[:,i] = motif[:,i] - mins
    sums = normed_motif.sum(axis = 1)
    for i in range(4):
        normed_motif[:,i] = normed_motif[:,i]/sums
    return normed_motif
motifs = current_model.get_layer('convolution_layer').get_weights()[0]
motifs.shape
index_seq_dict = dict(zip(range(4), ['A','C', 'G', 'T']))
seq_list = []
for i in range(motifs.shape[-1]):
    motif = motifs[:,:,i]
    normed_motif = normalize_motif(motif)
    
    sequence = ''.join([index_seq_dict[np.argmax(x)] for x in motif])
    seq_list.append(sequence)
    motif_to_plot = pd.DataFrame(motif, index = list(sequence), columns=['A','C','G','T'])
    normed_motif_to_plot = pd.DataFrame(normed_motif, index = list(sequence), columns=['A','C','G','T'])
    print(sequence)
    cg = sns.clustermap(motif_to_plot,
              row_cluster = False,
              col_cluster=False,
              cmap='bwr',
              square = True
              )
    plt.setp(cg.ax_heatmap.yaxis.get_majorticklabels(), rotation=0)
    plt.show()
    
    cg = sns.clustermap(normed_motif_to_plot,
                  row_cluster = False,
                  col_cluster=False,
                  cmap='bwr',
                  square = True
                  )
    plt.setp(cg.ax_heatmap.yaxis.get_majorticklabels(), rotation=0)
    plt.show()
    if i == 1:
        break<jupyter_output>CGAAACTGCAACTCGT<jupyter_text>Score motifs against library<jupyter_code>def readMotifFile(motifPath):
    with open(motifPath) as f:
        m = Bio.motifs.read(f, 'jaspar')
    name = m.name
    matrix = [m.pwm['A'], m.pwm['C'], m.pwm['G'], m.pwm['T']]
    matrix = np.array(matrix).T
    
    return (name,np.array(matrix), m.matrix_id)
def find_best_matching_motif(motif, motif_library):
    cleaned_motif_1 = (motif[0], cleanMatrix(motif[1]))
    similarity_scores = []
    counter = 0
    for other_motif in motif_library:
        cleaned_motif_2 = (other_motif[0], cleanMatrix(other_motif[1]))
        alignment_fwd, alignScore_fwd = local_align_motifs(cleaned_motif_1, cleaned_motif_2)
        r_fwd = calcCorrelation(alignment_fwd[0], alignment_fwd[1])
        alignment_rev, alignScore_rev = local_align_motifs(cleaned_motif_1, revCompMotif(cleaned_motif_2))
        r_rev = calcCorrelation(alignment_rev[0], alignment_rev[1])
        r = np.max([r_fwd, r_rev])
        similarity_scores.append(r)
    best_index = np.argmax(similarity_scores)
    best_score = similarity_scores[best_index]
    best_motif = motif_library[best_index]
    
    return best_score, best_motif, similarity_scores
def minmax_normalize_motif(motif):
    normed_motif = np.zeros(motif.shape)
    mins = motif.min(axis=1)
    for i in range(4):
        normed_motif[:,i] = motif[:,i] - mins
    distances = motif.max(axis = 1) - motif.min(axis=1)
    for i in range(4):
        normed_motif[:,i] = normed_motif[:,i]/distances
    return normed_motif
# read in motif library
motif_dir = '/home/jtao/code/tba/default_motifs/'
motif_library = []
for f in os.listdir(motif_dir):
    motif = readMotifFile(motif_dir + f)
    motif_library.append(motif)
treatment_conv_motif_dict = {}
for condition in sorted(condition_model_dict.keys()):
    treatment = condition.split('_')[1]
    print('***', treatment, '***')
    conv_motif_dict = {}
    current_model = condition_model_dict[condition]
    motifs = current_model.get_layer('convolution_layer').get_weights()[0]
    all_scores = []
    for i in range(motifs.shape[-1]):
        motif = motifs[:,:,i]

        kernel_name = treatment +  '_' + str(i)
        normed_motif = (kernel_name, normalize_motif(motif))
#         best_score, best_motif, similarity_scores = find_best_matching_motif(normed_motif, motif_library)
        best_score, best_motif, similarity_scores = find_best_matching_motif((kernel_name, 
                                                                              minmax_normalize_motif(motif)), motif_library)
        all_scores.append(best_score)
        print(i, best_motif[0], best_score)
        conv_motif_dict[kernel_name] = best_motif[0]
    sns.distplot(all_scores)
    plt.title(treatment)
    plt.xlabel('Pearson Corrlation of Best Matching Motif')
    treatment_conv_motif_dict[treatment] = conv_motif_dict
sns.distplot(current_model.get_layer('dense_layer').get_weights()[0])
import pickle
if len(treatment_conv_motif_dict) > 0:
    with open(working_directory + '/treatment_conv_motif_dict.pickle', 'wb') as handle:
        pickle.dump(treatment_conv_motif_dict, handle)
with open(working_directory + '/treatment_conv_motif_dict.pickle', 'rb') as handle:
    treatment_conv_motif_dict = pickle.load(handle)
motifs = current_model.get_layer('convolution_layer').get_weights()[0]

best_score, best_motif, similarity_scores = find_best_matching_motif(('motif_0', normalized_motifs[0]), motif_library)
print(best_motif[0], best_score)<jupyter_output>/home/jtao/software/anaconda3/lib/python3.5/site-packages/scipy/stats/stats.py:3021: RuntimeWarning: invalid value encountered in double_scalars
  r = r_num / r_den<jupyter_text>Visualize Motif Scores  Summarize mean and variance<jupyter_code>for treatment in sorted(treatment_motifScores_dict.keys()):
    variances = []
    means = []
    for i in range(scores.shape[-1]):
        scores = treatment_motifScores_dict[treatment]
        flattened_scores = scores[:,:,i].flatten()
        score_mean = np.mean(flattened_scores)
        score_var = np.var(flattened_scores)
        means.append(score_mean)
        variances.append(score_var)
    with sns.axes_style('ticks'):
        plt.scatter(means, 
                    variances,
                    edgecolor=None,
                    s=10)
        plt.xlim(0,0.04)
        plt.ylim(0,0.02)
        ax = plt.gca()
        ax.set_aspect(ax.get_xlim()[1]/ax.get_ylim()[1])

        plt.title(treatment + ' motif score summary')
        plt.xlabel('Mean Score')
        plt.ylabel('Score Variance')
        sns.despine()
        plt.show()<jupyter_output><empty_output><jupyter_text>Test for normality<jupyter_code>threshold = 0.01
pval_lists = []
for treatment in sorted(treatment_motifScores_dict.keys()):
    print(treatment,'***')
    log_pvals = []
    for i in range(scores.shape[-1]):
        scores = treatment_motifScores_dict[treatment]
        flattened_scores = scores[:,:,i].flatten()
        
        stat, pval = scipy.stats.normaltest(flattened_scores)
        log_pvals.append(pval)
        if pval > threshold:
            print(i, 'mean:', np.mean(flattened_scores), 'var:', np.var(flattened_scores))
    
    print(treatment, np.sum(np.array(log_pvals) > threshold))<jupyter_output>il4-24h ***
62 mean: 0.0 var: 0.0
il4-24h 1
kla-1h ***
2 mean: 0.0 var: 0.0
39 mean: 0.0 var: 0.0
42 mean: 0.0 var: 0.0
45 mean: 0.0 var: 0.0
50 mean: 0.0 var: 0.0
70 mean: 0.0 var: 0.0
114 mean: 0.0 var: 0.0
116 mean: 0.0 var: 0.0
kla-1h 8
veh ***
138 mean: 0.0 var: 0.0
139 mean: 0.0 var: 0.0
veh 2<jupyter_text>Fraction of positions with significant score<jupyter_code>threshold = 2
pval_lists = []
for treatment in sorted(treatment_motifScores_dict.keys()):
    print(treatment,'***')
    fractions = []
    for i in range(scores.shape[-1]):
        scores = treatment_motifScores_dict[treatment]
        flattened_scores = scores[:,:,i].flatten()
        if np.var(flattened_scores) > 0:
            zscores = scipy.stats.zscore(flattened_scores)
            sig_fraction = np.sum(zscores > threshold)

            fractions.append(sig_fraction/len(flattened_scores))
    sns.distplot(fractions)
    plt.title(treatment)
    plt.xlabel('Fraction of Loci with Significant Match')
    plt.ylabel('Frequency')
    plt.show()<jupyter_output>il4-24h ***<jupyter_text>Fraction of sequences with significant score<jupyter_code>threshold = 2
pval_lists = []
for treatment in sorted(treatment_motifScores_dict.keys()):
    print(treatment,'***')
    fractions = []
    for i in range(scores.shape[-1]):
        scores = treatment_motifScores_dict[treatment]
        current_scores = scores[:,:,i]
        if np.var(current_scores) > 0:
            zscores = scipy.stats.zscore(current_scores)
            sig_fraction = np.sum(zscores.max(axis=1) > threshold)

            fractions.append(sig_fraction/len(scores))
    sns.distplot(fractions)
    plt.title(treatment)
    plt.xlabel('Fraction of Sequences with Significant Match')
    plt.ylabel('Frequency')
    plt.show()<jupyter_output>il4-24h ***<jupyter_text>Convert Motif Scores to p-values<jupyter_code>len(flattened_scores)
np.sum( > 2)/len(flattened_scores)
sns.distplot(scipy.stats.zscore(flattened_scores))
scipy.stats.ttest_1samp([1,2,3],2f)<jupyter_output><empty_output><jupyter_text>Scratch<jupyter_code>plt.figure(figsize=(15,15))
sns.heatmap(treatment_attention_dict['kla-1h'][0], xticklabels=False, 
            yticklabels=False,
            square = True)
treatment_motifScores_dict['kla-1h'][0].shape
plt.figure(figsize=(15,15))
sns.clustermap(treatment_motifScores_dict['kla-1h'][0],
           xticklabels=False, yticklabels=False,
           square=True,
           z_score=1,
#            cmap='Greens',
           row_cluster = False,
           col_cluster = False)
attention = treatment_attention_dict['kla-1h'][0]
zscores = scipy.stats.zscore(treatment_motifScores_dict['kla-1h'][0], axis=0)
np.where(attention.max(axis=0) > 0.3)
treatment_conv_motif_dict['kla-1h']['kla-1h_'+str(np.argmax(zscores[94]))]
treatment_conv_motif_dict['kla-1h']['kla-1h_'+str(np.argmax(zscores[182]))]
zscores[attention.max(axis=0) > 0.3].shape
zscores = 
len(set(treatment_conv_motif_dict['kla-1h'].values()))
k<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11587685921826357, 'fraction_numerical': 0.020457577704205167, 'mean_word_length': 3.894316807738815, 'pattern_counts': {'":': 0, '<': 71, '<?xml version=': 0, '>': 84, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '343878', 'license': 'Apache-2.0', 'n_tokens_mistral': 7402, 'n_tokens_neox': 6791, 'n_words': 1073, 'path': 'Classification Attention Analysis.ipynb', 'repo_name': 'jenhantao/genomic_grammar', 'size': '20237'}