<jupyter_start><jupyter_text>![mrna](Images/mrna.jpg) [![stanford](https://img.shields.io/badge/OpenVaccine-COVID_19_mRNA_Vaccine_Degradation_Prediction-00c0ff?logo=Kaggle&logoColor=white&style=for-the-badge)](https://www.kaggle.com/c/stanford-covid-vaccine)  Introduction Winning the fight against the COVID-19 pandemic will require an effective vaccine that can be equitably and widely distributed. Building upon decades of research has allowed scientists to accelerate the search for a vaccine against COVID-19, but every day that goes by without a vaccine has enormous costs for the world nonetheless. We need new, fresh ideas from all corners of the world. Could online gaming and crowdsourcing help solve a worldwide pandemic? Pairing scientific and crowdsourced intelligence could help computational biochemists make measurable progress. mRNA VaccinemRNA vaccines have taken the lead as the fastest vaccine candidates for COVID-19, but currently, they face key potential limitations. One of the biggest challenges right now is how to design super stable messenger RNA molecules (mRNA). Conventional vaccines (like your seasonal flu shots) are packaged in disposable syringes and shipped under refrigeration around the world, but that is not currently possible for mRNA vaccines. mRNA Stability IssueResearchers have observed that RNA molecules have the tendency to spontaneously degrade. This is a serious limitation--a single cut can render the mRNA vaccine useless. Currently, little is known on the details of where in the backbone of a given RNA is most prone to being affected. Without this knowledge, current mRNA vaccines against COVID-19 must be prepared and shipped under intense refrigeration, and are unlikely to reach more than a tiny fraction of human beings on the planet unless they can be stabilized.  Building Predictive ModelWe need to design the model that will predict likely degradation rates at each base of an RNA molecule, trained on a subset of an Eterna dataset comprising over 3000 RNA molecules (which span a panoply of sequences and structures) and their degradation rates at each position. Improving Stability of mRNAImproving the stability of mRNA vaccines was a problem that was being explored before the pandemic but was expected to take many years to solve.Now, we must solve this deep scientific challenge in months, if not weeks, to accelerate mRNA vaccine research and deliver a refrigerator-stable vaccine against SARS-CoV-2, the virus behind COVID-19. ![Picture title](Images/image-20210621-125546.png)  Getting Ready<jupyter_code># Dataframe
import json
import pandas as pd
import numpy as np
# Visualization
import plotly.express as px
# Deeplearning
import tensorflow.keras.layers as L
import tensorflow as tf
# Sklearn
from sklearn.model_selection import train_test_split
#Setting seeds
tf.random.set_seed(2021)
np.random.seed(2021)
# This will tell us the columns we are predicting
target_cols = ['reactivity', 'deg_Mg_pH10', 'deg_Mg_50C', 'deg_pH10', 'deg_50C']
Model_Train = True # True if you want to Train model which take 1 hour to train.<jupyter_output><empty_output><jupyter_text>Metric ![Picture title](Images/image-20210621-122600.png)![Picture title](Images/image-20210621-122723.png)<jupyter_code>def MCRMSE(y_true, y_pred):## Monte Carlo root mean squared errors 
    colwise_mse = tf.reduce_mean(tf.square(y_true - y_pred), axis=1)
    return tf.reduce_mean(tf.sqrt(colwise_mse), axis=1)<jupyter_output><empty_output><jupyter_text>Loading Data [![stanford](https://img.shields.io/badge/Stanford_University-OpenVaccine-800080?logo=Kaggle&logoColor=white&style=for-the-badge)](https://www.kaggle.com/c/stanford-covid-vaccine/data)  Columns detail explaination - **id** - An arbitrary identifier for each sample.- **seq_scored** - (68 in Train and Public Test, 91 in Private Test) Integer value denoting the number of positions used in scoring with predicted values. This should match the length of reactivity, deg_* and *_error_* columns. Note that molecules used for the Private Test will be longer than those in the Train and Public Test data, so the size of this vector will be different.- **seq_length** - (107 in Train and Public Test, 130 in Private Test) Integer values, denotes the length of sequence. Note that molecules used for the Private Test will be longer than those in the Train and Public Test data, so the size of this vector will be different.- **sequence** - (1x107 string in Train and Public Test, 130 in Private Test) Describes the RNA sequence, a combination of A, G, U, and C for each sample. Should be 107 characters long, and the first 68 bases should correspond to the 68 positions specified in seq_scored (note: indexed starting at 0).- **structure** - (1x107 string in Train and Public Test, 130 in Private Test) An array of (, ), and . characters that describe whether a base is estimated to be paired or unpaired. Paired bases are denoted by opening and closing parentheses e.g. (....) means that base 0 is paired to base 5, and bases 1-4 are unpaired.- **reactivity** - (1x68 vector in Train and Public Test, 1x91 in Private Test) An array of floating point numbers, should have the same length as seq_scored. These numbers are reactivity values for the first 68 bases as denoted in sequence, and used to determine the likely secondary structure of the RNA sample.- **deg_pH10** - (1x68 vector in Train and Public Test, 1x91 in Private Test) An array of floating point numbers, should have the same length as seq_scored. These numbers are reactivity values for the first 68 bases as denoted in sequence, and used to determine the likelihood of degradation at the base/linkage after incubating without magnesium at high pH (pH 10).- **deg_Mg_pH10** - (1x68 vector in Train and Public Test, 1x91 in Private Test) An array of floating point numbers, should have the same length as seq_scored. These numbers are reactivity values for the first 68 bases as denoted in sequence, and used to determine the likelihood of degradation at the base/linkage after incubating with magnesium in high pH (pH 10).- **deg_50C** - (1x68 vector in Train and Public Test, 1x91 in Private Test) An array of floating point numbers, should have the same length as seq_scored. These numbers are reactivity values for the first 68 bases as denoted in sequence, and used to determine the likelihood of degradation at the base/linkage after incubating without magnesium at high temperature (50 degrees Celsius).- **deg_Mg_50C** - (1x68 vector in Train and Public Test, 1x91 in Private Test) An array of floating point numbers, should have the same length as seq_scored. These numbers are reactivity values for the first 68 bases as denoted in sequence, and used to determine the likelihood of degradation at the base/linkage after incubating with magnesium at high temperature (50 degrees Celsius).*_error_* - An array of floating point numbers, should have the same length as the corresponding reactivity or deg_* columns, calculated errors in experimental values obtained in reactivity and deg_* columns.- **predicted_loop_type** - (1x107 string) Describes the structural context (also referred to as 'loop type')of each character in sequence. Loop types assigned by bpRNA from Vienna RNAfold 2 structure. From the bpRNA_documentation: S: paired "Stem" M: Multiloop I: Internal loop B: Bulge H: Hairpin loop E: dangling End X: eXternal loopS/N filter Indicates if the sample passed filters described below in Additional Notes.  Observing Data<jupyter_code>data_dir = "Data/"
train = pd.read_json(data_dir + "train.json", lines=True)
test = pd.read_json(data_dir + "test.json", lines=True)
sample_df = pd.read_csv(data_dir + "sample_submission.csv")<jupyter_output><empty_output><jupyter_text>we have sequence, structure and predicted loop type which are text formate and we have to convert them into token to train the model. Then we have arrays with in columns from reactivity_error to deg_50C.<jupyter_code>train.head(2)
print('Train shapes: ', train.shape)
print('Test shapes: ', test.shape)<jupyter_output>Train shapes:  (2400, 19)
Test shapes:  (3634, 7)<jupyter_text>Test have sequence, structure, predicted_loop_type, seq_length, seq_scored and id    Signal to noise distribution<jupyter_code>fig = px.histogram(
    train, 
    "signal_to_noise", 
    nbins=25, 
    title='signal_to_noise distribution', 
    width=800,
    height=400
)
fig.show()<jupyter_output><empty_output><jupyter_text>Removing negetive values<jupyter_code>train = train.query("signal_to_noise >= 1")<jupyter_output><empty_output><jupyter_text>Sequence Test lengthWe need to divide our test dataset into two different data frames based on different sequence lengths. Splitting them into Public and private data frames.<jupyter_code>fig = px.histogram(
    test, 
    "seq_length", 
    nbins=25, 
    title='sequence_length distribution', 
    width=800,
    height=400
)
fig.show()<jupyter_output><empty_output><jupyter_text>Spliting Test into Public and Private DataframeSplitting test datasets into public and private data frames due to variation in sequence length. This has improved the overall result.<jupyter_code>public_df = test.query("seq_length == 107")
private_df = test.query("seq_length == 130")<jupyter_output><empty_output><jupyter_text>Pre-processing Data We will use this dictionary to map each character to an integer<jupyter_code>token2int = {x: i for i, x in enumerate("().<KEY>")}

token2int<jupyter_output><empty_output><jupyter_text>Converting DataFrame to 3D ArrayIt takes dataframe of shape (x, y), containing list of length `l` and returning numpy 3D array `(x,l,y)`.<jupyter_code>def dataframe_to_array(df):
   return np.transpose(np.array(df.values.tolist()), (0, 2, 1))<jupyter_output><empty_output><jupyter_text>Tokeniziation of SequenceIt takes dataframe then convert it into 3D numpy array and then tokenize it using Label Encoding.<jupyter_code>def dataframe_label_encoding(
    df, token2int, cols=["sequence", "structure", "predicted_loop_type"]
):
    return dataframe_to_array(
        df[cols].applymap(lambda seq: [token2int[x] for x in seq])
    )  ## tokenization of Sequence, Structure, Predicted loop<jupyter_output><empty_output><jupyter_text>Preprocessing Features and Labels Train input contains `["sequence", "structure", "predicted_loop_type"]` which need to tokenize and before that converting dataframe into 3D numpy array.<jupyter_code>train_inputs = dataframe_label_encoding(train, token2int)  ## Label encoding
train_labels = dataframe_to_array(train[target_cols])  ## dataframe to 3D array to<jupyter_output><empty_output><jupyter_text>Train & Validation split<jupyter_code>x_train, x_val, y_train, y_val = train_test_split(
    train_inputs, train_labels, test_size=0.1, random_state=34, stratify=train.SN_filter
)<jupyter_output><empty_output><jupyter_text>Preprocessing Public and Private Dataframe<jupyter_code>public_inputs = dataframe_label_encoding(public_df, token2int)
private_inputs = dataframe_label_encoding(private_df, token2int)<jupyter_output><empty_output><jupyter_text>Training / Evaluating Model  Build ModelI have experimented with multiple traditional models including Light GBM, Catboost, and BiLSTM, but the result was quite bad as compare to triple GRU layers. Using simple 3 Bidirectional GRU layer with linear activation. This model is quite simple and drived from [xhlulu](https://www.kaggle.com/xhlulu/openvaccine-simple-gru-model) inital model. I am big fan of his work and its just blew my mind on how he produced lower MCRMSE score with no feature engineering or data augumentation.  To learn more about RNNs, LSTM and GRU, please see [this blog post](https://colah.github.io/posts/2015-08-Understanding-LSTMs/).<jupyter_code>def build_model(
    embed_size,  # Length of uniqe tokens
    seq_len=107,  # public dataset seq_len
    pred_len=68,  # pred_len for public data
    dropout=0.5,  # trying best dropout (general)
    sp_dropout=0.2,  # Spatial Droput
    embed_dim=200,  # embeding dimention
    hidden_dim=256,  # hidden layer units
):
    inputs = L.Input(shape=(seq_len, 3))
    embed = L.Embedding(input_dim=embed_size, output_dim=embed_dim)(inputs)

    reshaped = tf.reshape(
        embed, shape=(-1, embed.shape[1], embed.shape[2] * embed.shape[3])
    )
    hidden = L.SpatialDropout1D(sp_dropout)(reshaped)

    # 3X BiGRU layers
    hidden = L.Bidirectional(
        L.GRU(
            hidden_dim,
            dropout=dropout,
            return_sequences=True,
            kernel_initializer="orthogonal",
        )
    )(hidden)

    hidden = L.Bidirectional(
        L.GRU(
            hidden_dim,
            dropout=dropout,
            return_sequences=True,
            kernel_initializer="orthogonal",
        )
    )(hidden)

    hidden = L.Bidirectional(
        L.GRU(
            hidden_dim,
            dropout=dropout,
            return_sequences=True,
            kernel_initializer="orthogonal",
        )
    )(hidden)

    # Since we are only making predictions on the first part of each sequence,
    # we have to truncate it
    truncated = hidden[:, :pred_len]
    out = L.Dense(5, activation="linear")(truncated)

    model = tf.keras.Model(inputs=inputs, outputs=out)
    model.compile(optimizer="Adam", loss=MCRMSE)  # loss function as of Eval Metric

    return model<jupyter_output><empty_output><jupyter_text>Building model<jupyter_code>model = build_model(
    embed_size=len(token2int) ## embed_size = 14
)  ## uniquie token in sequence, structure, predicted_loop_type
model.summary()<jupyter_output>Model: "model"
_________________________________________________________________
Layer (type)                 Output Shape              Param #   
=================================================================
input_1 (InputLayer)         [(None, 107, 3)]          0         
_________________________________________________________________
embedding (Embedding)        (None, 107, 3, 200)       2800      
_________________________________________________________________
tf.reshape (TFOpLambda)      (None, 107, 600)          0         
_________________________________________________________________
spatial_dropout1d (SpatialDr (None, 107, 600)          0         
_________________________________________________________________
bidirectional (Bidirectional (None, 107, 512)          1317888   
_________________________________________________________________
bidirectional_1 (Bidirection (None, 107, 512)          1182720   
_____________________________________________________________[...]<jupyter_text>Training ModelFitting our model on train dataset and saving best model checkpoint in Model folder.- batch size = 64- epochs = 40<jupyter_code>if Model_Train:
    history = model.fit(
        x_train,
        y_train,
        validation_data=(x_val, y_val),
        batch_size=64,
        epochs=40,
        verbose=2,
        callbacks=[
            tf.keras.callbacks.ReduceLROnPlateau(patience=5),
            tf.keras.callbacks.ModelCheckpoint("Model/model.h5"),
        ],
    )<jupyter_output>Epoch 1/40
30/30 - 69s - loss: 0.4536 - val_loss: 0.3796
Epoch 2/40
30/30 - 57s - loss: 0.3856 - val_loss: 0.3601
Epoch 3/40
30/30 - 57s - loss: 0.3637 - val_loss: 0.3410
Epoch 4/40
30/30 - 57s - loss: 0.3488 - val_loss: 0.3255
Epoch 5/40
30/30 - 57s - loss: 0.3357 - val_loss: 0.3188
Epoch 6/40
30/30 - 57s - loss: 0.3295 - val_loss: 0.3163
Epoch 7/40
30/30 - 57s - loss: 0.3200 - val_loss: 0.3098
Epoch 8/40
30/30 - 57s - loss: 0.3117 - val_loss: 0.2997
Epoch 9/40
30/30 - 57s - loss: 0.3046 - val_loss: 0.2899
Epoch 10/40
30/30 - 57s - loss: 0.2993 - val_loss: 0.2875
Epoch 11/40
30/30 - 57s - loss: 0.2919 - val_loss: 0.2786
Epoch 12/40
30/30 - 57s - loss: 0.2830 - val_loss: 0.2711
Epoch 13/40
30/30 - 57s - loss: 0.2777 - val_loss: 0.2710
Epoch 14/40
30/30 - 57s - loss: 0.2712 - val_loss: 0.2584
Epoch 15/40
30/30 - 57s - loss: 0.2640 - val_loss: 0.2580
Epoch 16/40
30/30 - 57s - loss: 0.2592 - val_loss: 0.2518
Epoch 17/40
30/30 - 57s - loss: 0.2540 - val_loss: 0.2512
Epoch 18/40
30/30 - 57s[...]<jupyter_text>Evaluate training historyThe validation MCRMSE score became flat on 35 but training MCRMSE kept reducing so I have to limit my training to 40 ephchs at most<jupyter_code>if Model_Train:
    fig = px.line(
        history.history,
        y=["loss", "val_loss"],
        labels={"index": "epoch", "value": "MCRMSE"},
        title="History",
    )
    fig.show()<jupyter_output><empty_output><jupyter_text>Loading models and making predictions Public and private sets have different sequence lengths, so we will preprocess them separately and load models of different tensor shapes. This is possible because RNN models can accept sequences of varying lengths as inputs. If the length of the sequence of an id is, e.g., 107, then you should make 107 predictions. Positions greater than the seq_scored value of a sample are not scored, but still need a value in the solution file.| id_seqpos      | reactivity | deg_Mg_pH10 | deg_pH10 | deg_Mg_50C | deg_50C ||----------------|------------|-------------|----------|------------|---------|| id_00073f8be_0 | 0.1        | 0.3         | 0.2      | 0.5        | 0.4     || id_00073f8be_1 | 0.3        | 0.2         | 0.5      | 0.4        | 0.2     || id_00073f8be_2 | 0.5        | 0.4         | 0.2      | 0.1        | 0.2     |  Loading weights Prediction formate required by competition organizers, `seq_length` is similar to `pred_len` so we are also creating two models with different `seq_len` and `pred_len`. After that loading our saved weights into both models.<jupyter_code>model_public = build_model(seq_len=107, pred_len=107, embed_size=len(token2int))
model_private = build_model(seq_len=130, pred_len=130, embed_size=len(token2int))

model_public.load_weights("Model/model.h5")
model_private.load_weights("Model/model.h5")<jupyter_output><empty_output><jupyter_text>Prediction<jupyter_code>public_preds = model_public.predict(public_inputs)
private_preds = model_private.predict(private_inputs)<jupyter_output><empty_output><jupyter_text>Checking the output shape<jupyter_code>private_preds.shape<jupyter_output><empty_output><jupyter_text>Post-processing and submit **Converting 3D numpy array into Dataframe:**- combining both **private(df,prediction )** and **public(df,prediction)** dataframe- Dividing unique id into 130 and 107 diferent samples `[id_00073f8be_0,id_00073f8be_1,id_00073f8be_2 ..]`- concating all of the data into 2D Dataframe and preparing for submission.<jupyter_code>preds_ls = []

for df, preds in [(public_df, public_preds), (private_df, private_preds)]:
    for i, uid in enumerate(df.id):

        single_pred = preds[i]

        single_df = pd.DataFrame(single_pred, columns=target_cols)
        single_df["id_seqpos"] = [f"{uid}_{x}" for x in range(single_df.shape[0])]

        preds_ls.append(single_df)

preds_df = pd.concat(preds_ls)
preds_df.head()<jupyter_output><empty_output><jupyter_text>Submissionmergining our predicted dataframe on `id_seqpos` so that all the predicted values follows the same id. The Creating submission file<jupyter_code>submission = sample_df[["id_seqpos"]].merge(preds_df, on=["id_seqpos"])
submission.to_csv("Submission/submission.csv", index=False)<jupyter_output><empty_output><jupyter_text>Observing submission file<jupyter_code>submission.head()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08923330415529582, 'fraction_numerical': 0.04114103290252819, 'mean_word_length': 4.428172163219676, 'pattern_counts': {'":': 2, '<': 95, '<?xml version=': 0, '>': 96, 'https://': 6, 'lorem ipsum': 0, 'www.': 3, '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': '24', 'dir': 'jupyter-structured-clean-dedup', 'id': '239362', 'license': 'Apache-2.0', 'n_tokens_mistral': 6409, 'n_tokens_neox': 5662, 'n_words': 2033, 'path': 'openvaccine-simple-gru-model.ipynb', 'repo_name': 'kingabzpro/mRNA-Vaccine-Degradation-Prediction', 'size': '19427'}