<jupyter_start><jupyter_text>Imports<jupyter_code>from __future__ import division
import pandas as pd
import numpy as np
from numpy import argmax
from scipy import constants
import random

import os
import sys
import re
import pdb
import glob

import suftware

from sklearn.preprocessing import LabelEncoder
from sklearn.preprocessing import OneHotEncoder
from sklearn.model_selection import train_test_split

import keras
from keras.models import Sequential, Model
from keras.layers import Dense, Dropout, Activation, Input, Lambda, Concatenate
from keras.optimizers import Adam
from keras import metrics
from keras import regularizers
from keras import callbacks
import keras.backend as K

import tensorflow as tf

from tensorflow.python.framework import ops
from tensorflow.python.ops import math_ops

import seaborn as sns
import matplotlib.pyplot as plt
import matplotlib.patches as patches
from matplotlib import cm
from matplotlib.colors import ListedColormap, LinearSegmentedColormap
%matplotlib inline
plt.ion()<jupyter_output>Using TensorFlow backend.<jupyter_text>Load the energy matrices from Kinney et al. 2010. These will be used later for comparison<jupyter_code>crp_df = pd.read_csv('data/energy_matrices/crp_tau_final_all.26.txt', sep = "\t", usecols = range(1,5))
rnap_df = pd.read_csv('data/energy_matrices/rnap_tau_final_all.41.txt', sep = "\t", usecols = range(1,5))

# Convert the dataframe to numpy array for downstream processing
crp_array = pd.DataFrame.to_numpy(crp_df)
rnap_array = pd.DataFrame.to_numpy(rnap_df)<jupyter_output><empty_output><jupyter_text>Helper function for one-hot encoding<jupyter_code># Fit a label encoder and a onehot encoder
bases = ["A","C","G","T"]
label_encoder = LabelEncoder()
label_encoder.fit(bases)
tmp = label_encoder.transform(bases)
tmp = tmp.reshape(len(tmp), 1)
onehot_encoder = OneHotEncoder(sparse = False)
onehot_encoder.fit(tmp)

# Encode sequence into onehot
def onehot_sequence(sequence, lab_encoder = label_encoder, one_encoder = onehot_encoder):
    """Sequence as a string"""
    tmp = lab_encoder.transform(list(sequence))
    tmp = tmp.reshape(len(tmp),1)
    tmp = one_encoder.transform(tmp)
    return tmp<jupyter_output><empty_output><jupyter_text>Load Input Data<jupyter_code>input_data_full_wt = pd.read_csv('data/full-wt/data.txt',delim_whitespace=True)<jupyter_output><empty_output><jupyter_text>Put data into lists and feed them to sklearn's train_test_split method<jupyter_code>sequences = []
bin_number = []

for _ in range(len(input_data_full_wt)):
    sequences.append(input_data_full_wt['seq'][_])
    bin_number.append(list(input_data_full_wt.iloc[_,1:len(input_data_full_wt.columns)].values))
X_train, X_test, y_train, y_test = train_test_split(sequences, bin_number, test_size=0.2)
######
## Define custom function for keras
######

def log_poisson_loss(y_true, y_pred):
    return tf.nn.log_poisson_loss(y_true,y_pred)

def poisson_loss(y_true, y_pred):
    y_pred = ops.convert_to_tensor(y_pred)
    y_true = math_ops.cast(y_true, y_pred.dtype)
    return K.mean(y_pred - y_true * math_ops.log(y_pred + K.epsilon()), axis=-1)<jupyter_output><empty_output><jupyter_text>One-hot encode the sequences at the location of CRP and RNAP binding sites<jupyter_code>input_seqs_ohe = []
for _ in range(len(X_train)):
    input_seqs_ohe.append(onehot_sequence(X_train[_][1:27]+X_train[_][34:75]))
    
test_seqs_ohe = []    
for _ in range(len(X_test)):
    test_seqs_ohe.append(onehot_sequence(X_test[_][1:27]+X_test[_][34:75]))<jupyter_output><empty_output><jupyter_text>Shapes of the training and test sets:<jupyter_code>print(np.shape(input_seqs_ohe))
print(np.shape(test_seqs_ohe))
from keras.layers import multiply

class CustomConnected(Dense):

    def __init__(self,units,connections, output_dim,**kwargs):
        
        self.output_dim = output_dim
        #this is matrix of connections. I.e., this should be [0,0,trainable,trainable]
        self.connections = connections                        
        #initalize the original Dense with all the usual arguments   
        super(CustomConnected,self).__init__(units,**kwargs)  

    def call(self, inputs):
        #print(self.connections)
        
        kernel_times_connections = multiply([self.kernel, self.connections])
        #output = K.dot(inputs, self.kernel * self.connections)
        
        # for implementing custom constraint, print out kernel. Set last two elements equal to each other.
        output = K.dot(inputs,kernel_times_connections)
        
        if self.use_bias:
            output = K.bias_add(output, self.bias)
        if self.activation is not None:
            output = self.activation(output)
        return output
    
    def compute_output_shape(self, input_shape):
        return (input_shape[0], self.output_dim)<jupyter_output><empty_output><jupyter_text>Model definitions<jupyter_code>number_input_layer_nodes = len(input_seqs_ohe[0])*4 +1   
inputTensor = Input((number_input_layer_nodes,),name='Sequence')

# define two groups which represent binding sites for the two proteins crp and rnap
group1 = Lambda(lambda x: x[:,0:104], output_shape=((104,)))(inputTensor)
group2 = Lambda(lambda x: x[:,104:268], output_shape=((164,)))(inputTensor)
group3 = Lambda(lambda x: x[:,268:269], output_shape=((1,)))(inputTensor)


# From input to first hidden layer. The groups are connected to two nodes
# node group1 represents energy of crp and group2 represents rnap energy.

# rename group1 to crp pwm

group1 = Dense(1,use_bias=True,name='CRP_pwm')(group1)
group2 = Dense(1,use_bias=True,name='RNAP_pwm')(group2)
group3 = Dense(1,use_bias=False,name='one-feature')(group3)

# the groups from the previous layer have to be concatenated:
concatenateLayer = Concatenate()([group1,group2,group3])

# gammaInput = Input(shape=(1,),name='gamma_Input')   # i.e. weights going into softmax: [0,0,0,1]

softmaxLayer = Dense(4,activation='softmax', trainable=False, use_bias=False,name='softmax_layer')(concatenateLayer)

tauTensor = CustomConnected(1,tf.constant([0.0,0.0,1.0,1.0]),1,use_bias=False,name='custom_tau_layer')(softmaxLayer)
intermediateTensor = Dense(10,activation='softmax', trainable=True, use_bias=False,name='softmax_layer_intermediate')(tauTensor)
outputTensor = Dense(np.shape(y_train[1])[0],activation='softmax', trainable=True, use_bias=False,name='softmax_layer_bins')(intermediateTensor)


#create the model:
model = Model(inputTensor,outputTensor)
model.summary()
model.layers<jupyter_output><empty_output><jupyter_text>Model layer 6 represents the value of $\Delta G_I$, check its value before training<jupyter_code>print(model.layers[6].get_weights())<jupyter_output>[array([[-0.6022438]], dtype=float32)]<jupyter_text>Set weights for $\tau$ layer...... so that only the two microstates (with bound RNAP) contribute to transcription<jupyter_code>tau_layer_weights = model.layers[9].get_weights()[0].copy()
tau_layer_weights.shape
tau_layer_initial_weights = np.array([0,0,1,1])
tau_layer_initial_weights = tau_layer_initial_weights.reshape(4,1)
model.layers[9].set_weights([tau_layer_initial_weights])
model.layers[9].get_weights()[0]<jupyter_output><empty_output><jupyter_text>The following are the non-trainable and trainable weights in the softmax layerand their values before training.<jupyter_code>print(model.layers[8].get_weights())<jupyter_output>[array([[-0.6904056 ,  0.5450921 ,  0.64685476, -0.87900007],
       [ 0.78645945,  0.36611283,  0.79684687,  0.9018415 ],
       [ 0.8125454 ,  0.10814333, -0.5667167 ,  0.81788135]],
      dtype=float32)]<jupyter_text>Set softmax layer weights manually, will be kept frozen during trainingNote that the bias weights here may be interpreted as gamma<jupyter_code># the following are the input weights to the softmax array.
# the first two arrays indicate 4 weights each coming out of the e_c and 
# e_r nodes. The last node indicates the gamma input, connecting only to
# microstate 4.

softmax_weights = np.array([[0,1,0,1],[0,0,1,1],[0,0,0,1]]) # this works

# syntax for set_weights expects weights in a list
softmax_weights = [softmax_weights]

# set weights
model.layers[8].set_weights(softmax_weights)
# see set weights
model.layers[8].get_weights()<jupyter_output><empty_output><jupyter_text>Appenindg a constant input of 1 to the input tensor. The weight of this tensor will correspond to $\Delta G_I$<jupyter_code>input_sequences_tensor2 = []
test_input_sequences_tensor2 = []

for _ in range(len(input_seqs_ohe)):
    temp = input_seqs_ohe[_].flatten()
    temp = np.append(temp,1)
    input_sequences_tensor2.append(temp)
    
for _ in range(len(test_seqs_ohe)):
    
    temp = test_seqs_ohe[_].flatten()
    temp = np.append(temp,1)
    test_input_sequences_tensor2.append(temp)<jupyter_output><empty_output><jupyter_text>Change inputs to TF/Keras to numpy arrays because that's the fit function expects.<jupyter_code>input_sequences_tensor2 = np.array(input_sequences_tensor2)
test_input_sequences_tensor2 = np.array(test_input_sequences_tensor2)
input_sequences_tensor2[0].shape
labels = np.array(y_train)
test_labels = np.array(y_test)
# Define custom function as follows

def poisson_loss(y_true, y_pred):
    y_pred = ops.convert_to_tensor(y_pred)
    y_true = math_ops.cast(y_true, y_pred.dtype)
    return K.mean(y_pred - y_true * math_ops.log(y_pred + K.epsilon()), axis=-1)

def log_poisson_loss(y_true, y_pred):
    return tf.nn.log_poisson_loss(y_true,y_pred)
model.compile(loss=log_poisson_loss,optimizer=Adam(lr=0.0005))
# fit model with early stopping call back
# esCallBack = keras.callbacks.EarlyStopping(monitor='val_loss',mode='auto')
# history = model.fit(input_sequences_tensor2, labels, 
#                     validation_split=0.2, epochs=300,verbose=1,callbacks=[esCallBack])  

history = model.fit(input_sequences_tensor2, labels, validation_split=0.2, epochs=300,verbose=1)<jupyter_output>Train on 32331 samples, validate on 8083 samples
Epoch 1/300
32331/32331 [==============================] - 1s 39us/step - loss: 1.0952 - categorical_accuracy: 0.1114 - val_loss: 1.0952 - val_categorical_accuracy: 0.1181
Epoch 2/300
32331/32331 [==============================] - 1s 31us/step - loss: 1.0952 - categorical_accuracy: 0.1115 - val_loss: 1.0951 - val_categorical_accuracy: 0.1181
Epoch 3/300
32331/32331 [==============================] - 1s 31us/step - loss: 1.0951 - categorical_accuracy: 0.1115 - val_loss: 1.0951 - val_categorical_accuracy: 0.1181
Epoch 4/300
32331/32331 [==============================] - 1s 31us/step - loss: 1.0950 - categorical_accuracy: 0.1165 - val_loss: 1.0948 - val_categorical_accuracy: 0.1636
Epoch 5/300
32331/32331 [==============================] - 1s 31us/step - loss: 1.0944 - categorical_accuracy: 0.1801 - val_loss: 1.0939 - val_categorical_accuracy: 0.1732
Epoch 6/300
32331/32331 [==============================] - 1s 31us/step - loss: 1.0934 - ca[...]<jupyter_text>plot losses<jupyter_code>plt.figure()
plt.plot(history.history['loss'],color='blue')
plt.plot(history.history['val_loss'],color='orange')
plt.title('Model loss',fontsize=12)
plt.ylabel('loss',fontsize=12)
plt.xlabel('epoch',fontsize=12)
plt.legend(['train', 'validation'])
#plt.savefig('Write_up/model_loss.png')
plt.show()
predictions_on_test = model.predict(test_input_sequences_tensor2)
model.evaluate(test_input_sequences_tensor2,test_labels)<jupyter_output>10104/10104 [==============================] - 0s 12us/step<jupyter_text>Visualize inferred parameters (PWMs)<jupyter_code>crp_weights = model.layers[4].get_weights()[0].copy()
rnap_weights = model.layers[5].get_weights()[0].copy()
crp_weights = crp_weights.reshape(26,4)
rnap_weights = rnap_weights.reshape(41,4)


plt.figure(figsize=(12,3))
sns.heatmap(crp_weights.T,cmap='RdBu_r')
plt.show()
plt.figure(figsize=(16,3))
sns.heatmap(rnap_weights.T,cmap='RdBu_r')
#sns.heatmap(rnap_weights.T,cmap='RdBu_r')
plt.show()<jupyter_output><empty_output><jupyter_text>visualize parameters as sequences logos using logomaker<jupyter_code>import logomaker
crp_logo_df = pd.DataFrame(crp_weights,columns=['A','C','G','T'])

crp_df.columns = ['A','C','G','T']
logomaker.Logo(-crp_df,center_values=True, font_name='Arial Rounded MT Bold')
plt.tick_params(labelsize=14)
plt.title('CRP PSAM Reported in Kinney et al. 2010',fontsize=14)
plt.ylabel("$-\Delta \Delta G$ ($k_{B}T$)", labelpad=-1)
logomaker.Logo(-crp_logo_df,center_values=True, font_name='Arial Rounded MT Bold')
plt.title('CRP PSAM learned by neural network (full-wt data).',fontsize=14)
plt.tick_params(labelsize=14)
plt.ylabel("$-\Delta \Delta G$ ($k_{B}T$)", labelpad=-1)

plt.show()
#logomaker.Logo(rnap_logo_df,center_values=True, font_name='Arial Rounded MT Bold')
rnap_logo_df = pd.DataFrame(rnap_weights,columns=['A','C','G','T'])
rnap_df.columns = ['A','C','G','T']
logomaker.Logo(-rnap_df,center_values=True, font_name='Arial Rounded MT Bold')
plt.tick_params(labelsize=14)
plt.title('RNAP PSAM Reported in Kinney et al. 2010',fontsize=14)
plt.ylabel("$-\Delta \Delta G$ ($k_{B}T$)", labelpad=-1)
logomaker.Logo(rnap_logo_df,center_values=True, font_name='Arial Rounded MT Bold')
plt.tick_params(labelsize=14)
plt.title('RNAP PSAM learned by neural network (full-wt data).',fontsize=14)
plt.ylabel("$-\Delta \Delta G$ ($k_{B}T$)", labelpad=-1)
plt.tick_params(labelsize=14)<jupyter_output><empty_output><jupyter_text>Inferred $\Delta G_I$<jupyter_code>model.layers[6].get_weights()<jupyter_output><empty_output><jupyter_text>value in kcal/mol<jupyter_code>-4.6800294*0.62<jupyter_output><empty_output><jupyter_text>Ensure that weights connecting $\Delta G$ nodes to microstate nodes stay fixed after training<jupyter_code>model.layers[8].get_weights()<jupyter_output><empty_output><jupyter_text>Weights of the $\tau$ layer, with the last two elements corresponding to $\tau_{max}$<jupyter_code>model.layers[9].get_weights()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.138504753673293, 'fraction_numerical': 0.041342552578507635, 'mean_word_length': 6.311743022643497, 'pattern_counts': {'":': 0, '<': 80, '<?xml version=': 0, '>': 80, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '478984', 'license': 'MIT', 'n_tokens_mistral': 5297, 'n_tokens_neox': 4873, 'n_words': 960, 'path': 'tareen_kinney_mlcb_2019.ipynb', 'repo_name': 'jbkinney/19_mlcb', 'size': '13884'}