<jupyter_start><jupyter_text>Application of deep learning approach to classify SARS-CoV-2 variants Coding Notebook Install and load the following libraries and packages<jupyter_code>import pandas as pd
import numpy as np
import os
import joblib
from keras import utils as np_utils
from sklearn.model_selection import train_test_split, GridSearchCV
from keras.models import Sequential
from keras.layers import Dense, Dropout, Flatten, Conv1D, MaxPooling1D, MaxPooling2D, Conv2D, LSTM, GRU, Bidirectional
from keras import regularizers
from keras.callbacks import EarlyStopping, ModelCheckpoint
from keras.wrappers.scikit_learn import KerasClassifier
import keras
from keras import utils as np_utils
from numpy import array
import seaborn as sns
from sklearn.metrics import confusion_matrix, classification_report, f1_score,ConfusionMatrixDisplay
print(keras.__version__)<jupyter_output>2.8.0<jupyter_text>Loading datasets and preprocessingdataset_2 - sequence dataset with over 26000 used for training and validating the model dataset_3 - sequences to validate the model<jupyter_code>dataset_2 = '../data/sequences.aln.fas' # huge dataset with 26k sequences
dataset_3 = '../data/validation.sequences.fas'
#Converting sequence data into data frame
#Ensure the fasta sequences are in one-liner and not blocks
#The sequence file had been pre-processed to include varinat names in the sequnce headers. Check file in fasta format.

#inputfile = [dataset_2,dataset_3]

def preprocess(inputfile):
    lines = []
    header = []
    with open (inputfile, 'r') as fasta_file:
    
        for line in fasta_file:
            if line.startswith(">"):
            #extract fasta headers
                lines.append(line.replace('\w', '').replace(">", "").replace("\n", ''))
            
            next
            #extract sequence lines
            if not line.startswith(">"):
                header.append(line.replace("\n",''))

            
    #Create dictionary of header and sequences and create a dataframe
    data = {lines[i]:header[i] for i in range(len(header))}
    data = pd.DataFrame(list(data.items()),columns = ['Seq_ID','sequence']) 

    data["sequence"] = data["sequence"].str.upper() #ensure the sequence is un uppercase. Important for encoding
    data["Variant_Id"] = data["Seq_ID"].apply(lambda i: i.split('/')[0]) #extracting variant names from sequnce file
    data = data.drop("Seq_ID", axis=1)
    
    return data



dataset_2 = preprocess(dataset_2)
dataset_3 = preprocess(dataset_3)<jupyter_output><empty_output><jupyter_text>Plot number of sequences for each variant<jupyter_code>dataset_2["Variant_Id"].value_counts().plot(kind="bar", 
                                            title = "Variant data", 
                                            fontsize=18,
                                           color = "blue")
dataset_3["Variant_Id"].value_counts().plot(kind="bar", 
                                            title = "Variant data", 
                                            fontsize=18,
                                           color = "blue")<jupyter_output><empty_output><jupyter_text>Assign variant name to numerical formgamma - 0 \delta - 1 \beta - 2 \alpha - 3 \omicron - 4<jupyter_code>def variant_to_numerical(sequence_data):
    """
    Function to map convetional variant names to numeric
    pass a sequence data with column named "Variant_ID".
    
    """
    
    variant = {"gamma" : 0, "delta" : 1, "beta" : 2, "alpha" : 3, "omicron":4}
    
    
    sequence_data["Variant_Id"] = sequence_data["Variant_Id"].map(variant)
    
    return sequence_data
    
variant_to_numerical(dataset_2)

dataset_2.head()
def encode_seq(seq_data):
    """
    Function to convert nucleotides to binary form and hot-encoded

    Arguments :
    seq_data = sequence dataframe with "sequence" and "Variant_ID" columns

    char : total number of unique characters in the whole dataset. This is used to 
    bin : array list of zeros corresponding to the number of characters

    """

    ntds = []
    for i in seq_data['sequence']:
        for seq in i:
            ntd = seq.strip()
            ntds.append(ntd.upper())
        
    char = list(set(ntds))
    
    bin = [0] * len(char)
    bin0  = [0] * len(char)
    bin1  = [0] * len(char)
    bin2  = [0] * len(char)
    bin3  = [0] * len(char)
    
    #Create a dictionary of arrays to the seqeunce characters
    Encode = {}

    for i in char:
        if i == 'A':
            bin0[0] = 1
            Encode[i] = bin0
        elif i == 'C':
            bin1[1] = 1
            Encode[i] = bin1
        elif i == 'G':
            bin2[2] = 1
            Encode[i] = bin2
        elif i == 'T':
            bin3[3] = 1
            Encode[i] = bin3
        else:
            Encode[i] = bin
    
      #convert the sequence characters to encoded form and append to the list

    encoded_list = []

    for seq in seq_data["sequence"]:
        encoded_list.append([Encode[x] for x in seq])

    return encoded_list<jupyter_output><empty_output><jupyter_text>Encoding sequences<jupyter_code>#call function to encode the sequences and convert to array form
X_mine_data_kmer_2 = np.array(encode_seq(dataset_2))

#np.save("../../my_array", X_mine_data_kmer) #just setting DIR where to store this. Next time no need to run the process again

#X_mine_data_kmer_2 = np.load("../../my_array.npy")

#set variable for 
y_mine_kmer_2 = dataset_2['Variant_Id']

#y_mine_kmer_2.shape
dataset_2.head()<jupyter_output><empty_output><jupyter_text>**Model Fitting**Split the dataset into training and validation (80:20) ration.<jupyter_code>#Split dataset into training and validation
X_train, X_test, y_train, y_test = train_test_split(X_mine_data_kmer_2, y_mine_kmer_2, random_state = 42, stratify = y_mine_kmer_2)

#have the data in tensor-readable form. I also don't know what that means.
y_train = keras.utils.np_utils.to_categorical(y_train)
y_test = keras.utils.np_utils.to_categorical(y_test)
X_train.shape
X_test.shape<jupyter_output><empty_output><jupyter_text>Set the model<jupyter_code>#Model parameter setting

cnn_model_c = Sequential()

#Layer 1 - feature extraction
#Performs dimension reduction

cnn_model_c.add(Conv1D(filters = 54, kernel_size = (6), activation = 'relu', 
                     input_shape = X_mine_data_kmer_2.shape[1:])) #ensure you change the shape

cnn_model_c.add(MaxPooling1D(pool_size= (5)))
cnn_model_c.add(Dropout(0.5))

#layer2
cnn_model_c.add(Conv1D(filters = 27, kernel_size = (3), activation = 'relu', padding = 'same'))
cnn_model_c.add(MaxPooling1D(pool_size= (3)))
cnn_model_c.add(Dropout(0.2))

#layer3
cnn_model_c.add(Conv1D(filters = 14, kernel_size = (2), activation = 'relu', padding = 'same'))
cnn_model_c.add(MaxPooling1D(pool_size= (3)))
cnn_model_c.add(Dropout(0.1))

#layer4
cnn_model_c.add(Conv1D(filters = 7, kernel_size = (2), activation = 'relu', padding = 'same'))#
cnn_model_c.add(MaxPooling1D(pool_size= (3)))
cnn_model_c.add(Bidirectional(GRU(512, activation = 'relu')))
cnn_model_c.add(Dropout(0.01))

#Flatten
cnn_model_c.add(Flatten()) #right place for flattening

#Hidden dense layer
cnn_model_c.add(Dense(256, activation = 'relu'))
cnn_model_c.add(Dense(128))#, activation = 'relu'))
cnn_model_c.add(Dense(64))#, activation = 'relu'))
cnn_model_c.add(Dense(32))#, activation = 'relu'))
cnn_model_c.add(Dense(16))#, activation = 'relu', kernel_regularizer = regularizers.l2(0.01)))
cnn_model_c.add(Dense(5, activation = 'softmax'))

#Compile model
cnn_model_c.compile(loss='categorical_crossentropy', 
              optimizer='adam',
              metrics=['accuracy'])

cnn_model_c.summary()
#Fit the model on the dataset

early_stop = keras.callbacks.EarlyStopping(monitor = 'val_accuracy', min_delta = 0.0005, patience=8, 
                                           restore_best_weights=True )

history = cnn_model_c.fit(X_train, y_train, batch_size = 1000, validation_data=(X_test, y_test), 
                        epochs=30)
# save the model to disk
#filename = 'cnn_model.sav'
#joblib.dump(cnn_model, 'cnn_model.sav')
 
# load the model from disk
#cnn_model_c = joblib.load(filename)<jupyter_output><empty_output><jupyter_text>Explore Results  a. Plot accuracy and loss values of the training dataset<jupyter_code># Plot the loss and accuracy curves for training

import matplotlib.pyplot as plt
plt.plot(history.history['loss'], color='b', label="train_loss")
plt.title("Training Loss")
plt.xlabel("Number of Epochs")
plt.ylabel("Percent(%)")
plt.xlim(0,30)
plt.legend()
plt.savefig("training_loss.png")
# Plot the test accuracy

plt.plot(history.history['accuracy'], color='g', label="train_accuracy")
plt.title("Train Accuracy")
plt.xlabel("Number of Epochs")
plt.ylabel("Percent(%)")
plt.legend()
plt.savefig("test_accuracy.png")<jupyter_output><empty_output><jupyter_text>b. Plot accuracy and loss values of the validation dataset<jupyter_code># Plot the accuracy for validation

plt.plot(history.history['val_accuracy'], color='g', label="validation_accuracy")
plt.title("Validation Accuracy")
plt.xlabel("Number of Epochs")
plt.ylabel("Percent(%)")
plt.legend()
plt.savefig("validation_accuracy.png")
# Plot the loss and accuracy curves for validation 

import matplotlib.pyplot as plt
plt.plot(history.history['val_loss'], color='b', label="validation loss")
plt.title("Validation Loss")
plt.xlabel("Number of Epochs")
plt.ylabel("Loss")
plt.legend()
plt.savefig("validation_loss.png")<jupyter_output><empty_output><jupyter_text>c. Confusion matrix of the validation dataset<jupyter_code># Predict the values from the validation dataset
Y_pred = cnn_model_c.predict(X_test)

# Convert predictions classes to one hot vectors 
Y_pred_classes = np.argmax(y_test,axis = 1) 

# Convert validation observations to one hot vectors
Y_true = np.argmax(y_test,axis = 1) 
# compute the confusion matrix
confusion_mtx = confusion_matrix(Y_true, Y_pred_classes) 
# plot the confusion matrix
f,ax = plt.subplots(figsize=(8, 8))


sns.heatmap(confusion_mtx, annot=True, 
            linewidths=0.01,cmap="Greens",linecolor="gray", 
            fmt= '.1f',ax=ax,
            xticklabels=["gamma","delta","beta","alpha","omicron"], 
            yticklabels=["gamma","delta","beta","alpha","omicron"])

plt.xlabel("Predicted Label")
plt.ylabel("True Label")
plt.title("Confusion Matrix")
plt.savefig("confusion.png")<jupyter_output><empty_output><jupyter_text>c. Other metrices<jupyter_code>print(confusion_mtx)<jupyter_output>[[1399    0    0    0    0]
 [   0 1374    0    0    0]
 [   0    0 1068    0    0]
 [   0    0    0 1378    0]
 [   0    0    0    0 1290]]<jupyter_text>Section 2: Evaluating Model on External dataset A summary of the dataset<jupyter_code>dataset_3["Variant_Id"].value_counts().plot(kind="bar", 
                                            title = "Variant data", 
                                            fontsize=18,
                                           color = "blue")
variant_to_numerical(dataset_3)

dataset_3.head()<jupyter_output><empty_output><jupyter_text>Perform Hot-encoding on the dataset<jupyter_code>#call function to encode the sequences and convert to array form
X_mine_data_kmer_3 = np.array(encode_seq(dataset_3))

#np.save("../../my_array_3", X_mine_data_kmer) #just setting DIR where to store this. Next time no need to run the process again

#X_mine_data_kmer_3 = np.load("../../my_array_3.npy")<jupyter_output><empty_output><jupyter_text>Load saved model and test on the new dataset<jupyter_code># load the model from disk
#cnn_model_loaded = joblib.load('cnn_model.sav')

#X_mine_data_kmer_3 = np.load("../../my_array_3.npy")

#X_mine_data_kmer_3.shape


# Predict the values from the validation dataset
Y_pred = cnn_model_c.predict(X_mine_data_kmer_3)


# Convert predictions classes to one hot vectors 
#set variable for

y_mine_kmer_3 = dataset_3['Variant_Id'] 
y_mine_kmer_3 = np.array(y_mine_kmer_3)
y_mine_kmer_3

Y_pred_values = (Y_pred > 0.5)
Y_pred_values = np.argmax(Y_pred, axis=1)


val = f1_score(Y_pred_values, y_mine_kmer_3,average='weighted' )
print("F1 score for the model is", val)<jupyter_output>F1 score for the model is 0.9840722740824901<jupyter_text>Compute and plot Confusion Matrix<jupyter_code># compute the confusion matrix
confusion_mtx_testing = confusion_matrix(y_mine_kmer_3, Y_pred_values)#_classes)
#confusion_mtx = ConfusionMatrixDisplay(confusion_mtx, display_labels = ["gamma", "delta", "beta", "alpha", "omicron"])

confusion_mtx_testing
# plot the confusion matrix
f,ax = plt.subplots(figsize=(8, 8))
sns.heatmap(confusion_mtx_testing, annot=True, 
            xticklabels=["gamma","delta","beta","alpha","omicron"], 
            yticklabels=["gamma","delta","beta","alpha","omicron"],
            linewidths=0.01,
            cmap="Greens",
            linecolor="gray", fmt= '.1f',ax=ax)
            
plt.xlabel("Predicted Label")
plt.ylabel("True Label")
plt.title("Confusion Matrix")
plt.savefig("confusion.png")

#gamma" : 0, "delta" : 1, "beta" : 2, "alpha" : 3, "omicron":4
print(confusion_mtx_testing)
cmx = confusion_matrix(y_mine_kmer_3, Y_pred_values, normalize='all')

print(classification_report(y_mine_kmer_3, Y_pred_values))<jupyter_output>precision    recall  f1-score   support

           0       1.00      0.99      0.99       437
           1       1.00      0.94      0.97       494
           2       0.97      1.00      0.98      1403
           3       0.99      0.99      0.99       625
           4       1.00      1.00      1.00       103

    accuracy                           0.98      3062
   macro avg       0.99      0.98      0.99      3062
weighted avg       0.98      0.98      0.98      3062
{'fraction_non_alphanumeric': 0.11605118829981718, 'fraction_numerical': 0.02676416819012797, 'mean_word_length': 3.4547231270358307, 'pattern_counts': {'":': 2, '<': 57, '<?xml version=': 0, '>': 61, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 27, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '15', 'dir': 'jupyter-structured-clean-dedup', 'id': '102868', 'license': 'MIT', 'n_tokens_mistral': 4788, 'n_tokens_neox': 4519, 'n_words': 1080, 'path': 'notebook/.ipynb_checkpoints/notebook-checkpoint.ipynb', 'repo_name': 'omicscodeathon/ml_sarscov2', 'size': '13675'}