# ---
# jupyter:
#   jupytext:
#     text_representation:
#       extension: .py
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.14.4
#   kernelspec:
#     display_name: Python 3
#     language: python
#     name: python3
# ---

# # MLP cross validation
# Classify pc vs nc RNA.
# Use K=2,3,4.
#
# Set aside the 20% test set, stratified by length.
# On the remaining 80%,
# perform 5-fold cross validation.

# +
import numpy as np
import pandas as pd
import pandas as pd
import matplotlib.pyplot as plt
import tensorflow as tf
from tensorflow import keras
from sklearn.model_selection import cross_val_score
from sklearn.model_selection import RepeatedKFold
from tensorflow.keras.wrappers.scikit_learn import KerasRegressor

tf.keras.backend.set_floatx('float32')


# -

# ## K-mer frequency, K=2

def read_features(nc_file,pc_file):
    nc = pd.read_csv (nc_file)
    pc = pd.read_csv (pc_file)
    nc['class']=0
    pc['class']=1
    rna_mer=pd.concat((nc,pc),axis=0)
    return rna_mer
rna_mer = read_features('ncRNA.2mer.features.csv','pcRNA.2mer.features.csv')
rna_mer


# Split into train/test stratified by sequence length.
def sizebin(df):
    return pd.cut(df["seqlen"],
                              bins=[0,1000,2000,4000,8000,16000,np.inf],
                              labels=[0,1,2,3,4,5])
def make_train_test(data):
    bin_labels= sizebin(data)
    from sklearn.model_selection import StratifiedShuffleSplit
    splitter = StratifiedShuffleSplit(n_splits=1, test_size=0.2, random_state=37863)
    # split(x,y) expects that y is the labels. 
    # Trick: Instead of y, give it it the bin labels that we generated.
    for train_index,test_index in splitter.split(data,bin_labels):
        train_set = rna_mer.iloc[train_index]
        test_set = rna_mer.iloc[test_index]
    return (train_set,test_set)
(train_set,test_set)=make_train_test(rna_mer)
print("train,test")
train_set.shape,test_set.shape


# +
def prepare_test_set(test_set):
    y_test=   test_set[['class']].copy()
    X_test=   test_set.div(test_set['seqlen'],axis=0)
    X_test=   X_test.drop(columns=['class','seqnum','seqlen'])
    return (X_test,y_test)
(X_test,y_test)=prepare_test_set(test_set)

def prepare_train_set(train_set):
    y_train_all=  train_set[['class']].copy()
    X_train_all=  train_set.div(train_set['seqlen'],axis=0)
    X_train_all=  X_train_all.drop(columns=['class','seqnum','seqlen'])
    return (X_train_all,y_train_all)

(X_train_all,y_train_all)=prepare_train_set(train_set)
X_train_all.shape,y_train_all.shape

# +
act="sigmoid"
def builder():
    mlp = keras.models.Sequential([
        keras.layers.LayerNormalization(trainable=False),
        keras.layers.Dense(32, activation=act,dtype='float32'),
        keras.layers.Dense(32, activation=act,dtype='float32'),
        keras.layers.Dense(1,  activation=act,dtype='float32')
    ])
    bc=tf.keras.losses.BinaryCrossentropy(from_logits=False)
    mlp.compile(loss=bc, optimizer="Adam",metrics=["accuracy"])
    return mlp

estimator= KerasRegressor( #batch_size=10
    build_fn=builder, epochs=50, verbose=1
)
kfold= RepeatedKFold(n_splits=5, n_repeats=1)
results= cross_val_score(
    estimator, X_train_all, y_train_all, cv=kfold, n_jobs=2)
results

# +
# See page 302 for explanation of these parameters.
# See also the keras docs e.g. 
# https://www.tensorflow.org/api_docs/python/tf/keras/losses/sparse_categorical_crossentropy
# Note keras can take parameters for % train vs % validation.

# It seems the BinaryCrossentropy assumes labels are probabilities.
# Instead of loss="binary_crossentropy",
bc=tf.keras.losses.BinaryCrossentropy(from_logits=False)
# Tried optimizers SGD, Adam
mlp2mer.compile(loss=bc, optimizer="Adam",metrics=["accuracy"])

# With one dense layer and Adam optimizer, accuracy increases slowly.
history2mer = mlp2mer.fit(X_train,y_train,epochs=100,validation_data=(X_valid,y_valid))
# -

pd.DataFrame(history2mer.history).plot(figsize=(8,5))
plt.grid(True)
plt.gca().set_ylim(0,1)
plt.show()

# ## K-mer frequency, K=3

rna_mer = read_features('ncRNA.3mer.features.csv','pcRNA.3mer.features.csv')
(train_set,test_set)=make_train_test(rna_mer)
(X_train,y_train,X_valid,y_valid)=prepare_train_set(train_set)

act="sigmoid"
mlp3mer = keras.models.Sequential([
    keras.layers.LayerNormalization(trainable=False),
    keras.layers.Dense(32, activation=act,dtype='float32'),
    keras.layers.Dense(32, activation=act,dtype='float32'),
    keras.layers.Dense(1,  activation=act,dtype='float32')
])
bc=tf.keras.losses.BinaryCrossentropy(from_logits=False)
mlp3mer.compile(loss=bc, optimizer="Adam",metrics=["accuracy"])
history3mer = mlp3mer.fit(X_train,y_train,epochs=100,validation_data=(X_valid,y_valid))

pd.DataFrame(history3mer.history).plot(figsize=(8,5))
plt.grid(True)
plt.gca().set_ylim(0,1)
plt.show()

# ## K-mer frequency, K=4

rna_mer = read_features('ncRNA.4mer.features.csv','pcRNA.4mer.features.csv')
(train_set,test_set)=make_train_test(rna_mer)
(X_train,y_train,X_valid,y_valid)=prepare_train_set(train_set)
act="sigmoid"
mlp4mer = keras.models.Sequential([
    keras.layers.LayerNormalization(trainable=False),
    keras.layers.Dense(32, activation=act,dtype='float32'),
    keras.layers.Dense(32, activation=act,dtype='float32'),
    keras.layers.Dense(1,  activation=act,dtype='float32')
])
bc=tf.keras.losses.BinaryCrossentropy(from_logits=False)
mlp4mer.compile(loss=bc, optimizer="Adam",metrics=["accuracy"])
history4mer = mlp4mer.fit(X_train,y_train,epochs=100,validation_data=(X_valid,y_valid))

pd.DataFrame(history4mer.history).plot(figsize=(8,5))
plt.grid(True)
plt.gca().set_ylim(0,1)
plt.show()



StarcoderdataJupyterScriptsDedupFiltered