{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\n\n# Load just 20000 rows - the entire input data is massive!\ndf = pd.read_csv('/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv', nrows = 20000)\ndf.head()","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-10-06T21:17:44.554828Z","iopub.execute_input":"2023-10-06T21:17:44.555170Z","iopub.status.idle":"2023-10-06T21:17:46.131772Z","shell.execute_reply.started":"2023-10-06T21:17:44.555142Z","shell.execute_reply":"2023-10-06T21:17:46.130677Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's take a look at the data, shall we?\n# Don't worry about understanding this code - it's just generally helpful to visualize your data when\n# starting a Kaggle competition.\nimport collections\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nchar_frequencies = [sum(df[\"sequence\"].apply(lambda x: x.count(char))) for char in ['G','A','U','C']] \n\nsns.set_theme()\nplt.pie(char_frequencies, labels = ['G','A','U','C'], autopct='%1.1f%%')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:17:46.133612Z","iopub.execute_input":"2023-10-06T21:17:46.134657Z","iopub.status.idle":"2023-10-06T21:17:47.023677Z","shell.execute_reply.started":"2023-10-06T21:17:46.134618Z","shell.execute_reply":"2023-10-06T21:17:47.022377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To understand what we have here, we actually have **timeseries data.** \n### For each sequence of mRNA, we want to predict its reactivity at each position:\n\n```\n    C        A        U         G        A    ...\n  0.61     0.42     0.01      0.19     0.99   ...\n```\n\nOur data is formatted with many columns, but only a few of them are informative. \n### In fact, we only have one piece of data to condition our prediction on:\n1. The mRNA sequence itself.\n\nAside from that, we also have a few pieces of metadata to help us decide how to train our model.\n\n2. Reads: Number of reads in the sequencing experiment that were assigned to the RNA sequence, and whose mutations were tabulated to compile the reactivity profile. (higher = more likely that this mRNA sequence data is correct).\n3. signal_to_noise: mean(measurement value)/mean(statistical error in measurement value). Higher = smaller spread for the reactivity target.\n4. SN_filter: a boolean of whether the sequence has > 100 reads and > 1 signal to noise. Basically if the data is good quality or not.\n5. reactivity error: a measure of (I'm assuming expected) error for the reactivity at each position. High error at a position = large spread = the reactivity at that position is less likely to be good training data.\n\nOur target is a sequence of the reactivities at each position in the mRNA sequence. The columns [reactivity_0001... reactivity_n], where n is the length of your sequence, contains the training target. \n\nUnfortunately the reactivities of the positions near the starts and ends of all sequences are unable to be scanned for technical reasons, so we'll have to figure out some way to impute data for those targets.","metadata":{}},{"cell_type":"markdown","source":"## Sample solution: Naive RNN\nWe'll build a recurrent neural network that will take in mRNA sequences and spit out a sequence of the same length that denotes the reactivity at each position in the sequence.\n\nThis sucks for multiple reasons and is strictly worse than attention (since the network can't find position-relative patterns in the sequence).\n\nNote that this model is for example purposes only - please don't use this model as a blueprint when submitting competitive code 👀.","metadata":{}},{"cell_type":"code","source":"import torch as t\ndf = df.fillna(0)\nmapping = {'G': [1, 0, 0, 0],\n           'A': [0, 1, 0, 0],\n           'U': [0, 0, 1, 0],\n           'C': [0, 0, 0, 1]}\ninput_sequences = np.stack(df[\"sequence\"].apply(lambda x: np.array([mapping[i] for i in x])))\noutput_sequences = df[[\"reactivity_\"+str(i).zfill(4) for i in range(1, 171)]].values[:, :, np.newaxis]\n\n# Convert numpy arrays to PyTorch tensors\ninput_sequences = t.from_numpy(input_sequences).float()\noutput_sequences = t.from_numpy(output_sequences).float()","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:17:52.931141Z","iopub.execute_input":"2023-10-06T21:17:52.931470Z","iopub.status.idle":"2023-10-06T21:18:01.166213Z","shell.execute_reply.started":"2023-10-06T21:17:52.931442Z","shell.execute_reply":"2023-10-06T21:18:01.165139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Let's replace all the NaNs with 0s. This is bad practice - it's much better to impute missing values\n# with the sample mean or a regression-imputed value instead. Try improving this on your own time!\n\nimport torch as t\nimport torch.nn as nn\nimport torch.nn.functional as F\n\nclass NaiveRNN(nn.Module):\n    def __init__(self, input_size, hidden_size, output_size, num_layers):\n        super(NaiveRNN, self).__init__()\n        self.hidden_size = hidden_size\n        self.num_layers = num_layers\n\n        self.rnn = nn.RNN(input_size, hidden_size, num_layers, batch_first=True)\n        self.fc = nn.Linear(hidden_size, output_size)\n\n    def forward(self, x):\n        h0 = t.zeros(self.num_layers, x.size(0), self.hidden_size).to(x.device) \n\n        out, _ = self.rnn(x, h0)  \n        out = self.fc(out)  \n        return out\n\n# Initialize the network\nmodel = NaiveRNN(4, 512, 1, 3)","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:18:03.632615Z","iopub.execute_input":"2023-10-06T21:18:03.633101Z","iopub.status.idle":"2023-10-06T21:18:03.667154Z","shell.execute_reply.started":"2023-10-06T21:18:03.633072Z","shell.execute_reply":"2023-10-06T21:18:03.666291Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Using GPUs with Pytorch\nWe heavily recommend using GPU acceleration, which should be on by default in this notebook.","metadata":{}},{"cell_type":"code","source":"if t.cuda.is_available():    \n    device = t.device(\"cuda:0\")\n    print('There are %d GPU(s) available.' % t.cuda.device_count())\nelse:\n    print('No GPU found.')\n    device = t.device(\"cpu\")\n    \nmodel = nn.DataParallel(model, device_ids=[0, 1]) # Parallelize across 2 gpus\nmodel = model.to(device)","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:18:08.760916Z","iopub.execute_input":"2023-10-06T21:18:08.761251Z","iopub.status.idle":"2023-10-06T21:18:15.936949Z","shell.execute_reply.started":"2023-10-06T21:18:08.761226Z","shell.execute_reply":"2023-10-06T21:18:15.935957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Loss + Optimizer\ncriterion = nn.MSELoss()\noptimizer = t.optim.Adam(model.parameters(), lr=5e-4)\n\nfrom sklearn.model_selection import train_test_split\nX_train, X_test, y_train, y_test = train_test_split(input_sequences, output_sequences, test_size = 0.2)\n\nclass sequence_dataset(t.utils.data.Dataset):\n    def __init__(self, input_sequences, output_sequences):\n        self.input_sequences = input_sequences; self.output_sequences = output_sequences\n\n    def __len__(self):\n        assert len(self.input_sequences) == len(self.output_sequences)\n        return len(self.input_sequences)\n\n    def __getitem__(self, idx):\n        return self.input_sequences[idx], self.output_sequences[idx]\n    \nBATCH_SIZE = 32 * 2\ntrain_dl=t.utils.data.DataLoader(sequence_dataset(X_train, y_train), batch_size=BATCH_SIZE, shuffle=True)\ntest_dl=t.utils.data.DataLoader(sequence_dataset(X_test, y_test), batch_size=BATCH_SIZE, shuffle=True)\n\n# Training loop\nfrom tqdm import tqdm, trange\nepochs = 20; ema_loss = None\ntrain_history, test_history = [], []\nfor epoch in range(epochs):\n    print(\"Epoch %d\"%(epoch+1))\n    with tqdm(train_dl) as pbar:\n        for input_sequence, output_sequence in pbar:\n            input_sequence = input_sequence.to(device)\n            output_sequence = output_sequence.to(device)\n            # Zero the parameter gradients\n            optimizer.zero_grad()\n\n            # Forward pass\n            outputs = model(input_sequence, )\n            loss = criterion(outputs, output_sequence)\n\n            # Backward pass and optimization\n            loss.backward()\n            optimizer.step()\n            if ema_loss is None: ema_loss = loss.item()\n            else: ema_loss = ema_loss * 0.98 + loss.item() * 0.02\n            pbar.set_description(\"Training Loss: %f\"%ema_loss)\n            train_history.append(loss.item())\n    with tqdm(test_dl) as pbar:\n        for input_sequence, output_sequence in pbar:\n            input_sequence = input_sequence.to(device)\n            output_sequence = output_sequence.to(device)\n            # Forward pass\n            with t.no_grad():\n                outputs = model(input_sequence)\n                loss = criterion(outputs, output_sequence)\n\n                pbar.set_description(\"Validation Loss: %f\"%loss)\n            test_history.append(loss.item())","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:18:15.938850Z","iopub.execute_input":"2023-10-06T21:18:15.939222Z","iopub.status.idle":"2023-10-06T21:21:32.462096Z","shell.execute_reply.started":"2023-10-06T21:18:15.939187Z","shell.execute_reply":"2023-10-06T21:21:32.461098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import seaborn as sns\ndef ema(data, smoothing = 0.98):\n    smoothed_data = []\n    val = None\n    for item in data:\n        if val is None: val = item\n        else: val = val * smoothing + item * (1 - smoothing)\n        smoothed_data.append(val)\n    return smoothed_data\n\nfig, axes = plt.subplots(2, 1, figsize = (8, 5))\nsns.lineplot(ema(train_history), ax = axes[0], label = \"Train loss\")\nsns.lineplot(ema(test_history), ax = axes[1], label = \"Val loss\")","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:21:39.620188Z","iopub.execute_input":"2023-10-06T21:21:39.620537Z","iopub.status.idle":"2023-10-06T21:21:40.307660Z","shell.execute_reply.started":"2023-10-06T21:21:39.620509Z","shell.execute_reply":"2023-10-06T21:21:40.306710Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now that we've got a model, let's make predictions!","metadata":{}},{"cell_type":"code","source":"test_input_df = pd.read_csv(\"/kaggle/input/stanford-ribonanza-rna-folding/test_sequences.csv\")\n\ntest_input_df = test_input_df.fillna(0)\nmapping = {'G': [1, 0, 0, 0],\n           'A': [0, 1, 0, 0],\n           'U': [0, 0, 1, 0],\n           'C': [0, 0, 0, 1]}\ninput_sequences = test_input_df[\"sequence\"].apply(lambda x: t.Tensor([mapping[i] for i in x]).to(device))","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:21:58.522548Z","iopub.execute_input":"2023-10-06T21:21:58.523267Z","iopub.status.idle":"2023-10-06T21:25:28.957366Z","shell.execute_reply.started":"2023-10-06T21:21:58.523233Z","shell.execute_reply":"2023-10-06T21:25:28.956265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import csv\nwith open(\"output.csv\", \"a\") as csv_file:\n    writer = csv.writer(csv_file, delimiter=',')\n    writer.writerow(['id','reactivity_DMS_MaP','reactivity_2A3_MaP'])\n    start_idx = 0\n    for input_sequence in tqdm(input_sequences):\n        # Forward pass\n        with t.no_grad():\n            output = model(input_sequence.unsqueeze(dim = 0))\n            output_list = output.squeeze().cpu()\n            writer.writerows([[start_idx + index, i.item(), i.item()] for index, i in enumerate(output_list)])\n            start_idx += len(output_list)","metadata":{"execution":{"iopub.status.busy":"2023-10-06T21:37:23.009229Z","iopub.execute_input":"2023-10-06T21:37:23.009632Z","iopub.status.idle":"2023-10-06T21:37:26.521329Z","shell.execute_reply.started":"2023-10-06T21:37:23.009599Z","shell.execute_reply":"2023-10-06T21:37:26.519858Z"},"trusted":true},"execution_count":null,"outputs":[]}]}