{"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":"markdown","source":"# Ribonanza RNA - Simple EDA 🧬🔬\n\nA very simple notebook I did for myself to get started in the competion and decide to share.\n\n## Initial Hypotesis\n\nRNA lengths differ across train and test samples (🚨 spoiler alert: they do).\n\n\n## Exploratory Data Analysis\n\nLoading required Python libraries.","metadata":{}},{"cell_type":"code","source":"import matplotlib.patches as mpatches\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nimport numpy as np \nimport pandas as pd \nimport os, gc, re\n\nCMAP = plt.cm.Pastel2","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-09-19T13:22:17.393467Z","iopub.execute_input":"2023-09-19T13:22:17.393844Z","iopub.status.idle":"2023-09-19T13:22:19.823380Z","shell.execute_reply.started":"2023-09-19T13:22:17.393815Z","shell.execute_reply":"2023-09-19T13:22:19.822116Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Peeking training data - first 10 rows:","metadata":{}},{"cell_type":"code","source":"train_file = '/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv'\ntrain = pd.read_csv(train_file,\n                    nrows = 10,\n                   )\ntrain.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:19.825438Z","iopub.execute_input":"2023-09-19T13:22:19.825941Z","iopub.status.idle":"2023-09-19T13:22:19.926672Z","shell.execute_reply.started":"2023-09-19T13:22:19.825908Z","shell.execute_reply":"2023-09-19T13:22:19.925343Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As expected there are many columns, `reactivity_*` are the ones the target come from, also `signal_to_noise` but here the focous will be upon `sequence`. For test set only the later column was loaded.","metadata":{}},{"cell_type":"code","source":"test_file = '/kaggle/input/stanford-ribonanza-rna-folding/test_sequences.csv'\ntest = pd.read_csv(test_file,\n                    usecols = ['sequence'],\n                   )\ntest.head()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:19.928286Z","iopub.execute_input":"2023-09-19T13:22:19.928990Z","iopub.status.idle":"2023-09-19T13:22:27.238307Z","shell.execute_reply.started":"2023-09-19T13:22:19.928943Z","shell.execute_reply":"2023-09-19T13:22:27.236913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Organizing columns by category:","metadata":{}},{"cell_type":"code","source":"# regex was very useful to separate reactivity_* and reactivity_error_* cols\nreactivity_cols = [\n    col for col in train.columns \n    if bool(re.search('^reactivity_\\d{4}$', col))\n]\n\nerror_cols = [\n    col for col in train.columns \n    if bool(re.search('^reactivity_error_\\d{4}$', col))\n]\n\nqualify_cols = [\n    'sequence_id',\n    'sequence',\n    'experiment_type'\n]\n\nquality_cols = [\n    'reads',\n    'signal_to_noise',\n    'SN_filter',\n]","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:27.241578Z","iopub.execute_input":"2023-09-19T13:22:27.242085Z","iopub.status.idle":"2023-09-19T13:22:27.252036Z","shell.execute_reply.started":"2023-09-19T13:22:27.242024Z","shell.execute_reply":"2023-09-19T13:22:27.250590Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Describing the columns later used as target:","metadata":{}},{"cell_type":"code","source":"train[reactivity_cols].describe()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:27.253964Z","iopub.execute_input":"2023-09-19T13:22:27.254439Z","iopub.status.idle":"2023-09-19T13:22:27.687716Z","shell.execute_reply.started":"2023-09-19T13:22:27.254397Z","shell.execute_reply":"2023-09-19T13:22:27.686574Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"At frist glance it may seen as they are all empty, but as the sponsor stated in the data tab (and as I show later) [\"several positions early and late in the sequence also cannot be probed due to technical reasons, and their reactivity values are `null.`\"](https://www.kaggle.com/competitions/stanford-ribonanza-rna-folding/data)","metadata":{}},{"cell_type":"code","source":"train[error_cols].describe()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:27.689189Z","iopub.execute_input":"2023-09-19T13:22:27.690232Z","iopub.status.idle":"2023-09-19T13:22:28.113636Z","shell.execute_reply.started":"2023-09-19T13:22:27.690192Z","shell.execute_reply":"2023-09-19T13:22:28.112345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The same thing can be seen for the error related columns.\n\nDescribing quality related columns:","metadata":{}},{"cell_type":"code","source":"train[quality_cols].describe()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.115207Z","iopub.execute_input":"2023-09-19T13:22:28.115749Z","iopub.status.idle":"2023-09-19T13:22:28.138305Z","shell.execute_reply.started":"2023-09-19T13:22:28.115718Z","shell.execute_reply":"2023-09-19T13:22:28.137094Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looking at how many `reactivity_*` columns there are:","metadata":{}},{"cell_type":"code","source":"train[reactivity_cols].shape","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.139741Z","iopub.execute_input":"2023-09-19T13:22:28.140105Z","iopub.status.idle":"2023-09-19T13:22:28.149249Z","shell.execute_reply.started":"2023-09-19T13:22:28.140076Z","shell.execute_reply":"2023-09-19T13:22:28.148090Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The reason we saw only nulls while using [`pandas.DataFrame.describe`](https://pandas.pydata.org/docs/reference/api/pandas.DataFrame.describe.html) is due to Pandas don't showing the in-between columns as default (for good, there are 206 columns).\n\nBut following we can see a list with booleans wheter any of the first 10 rows has at list a non null value:","metadata":{}},{"cell_type":"code","source":"(~train[reactivity_cols].isnull()).any().values","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.151060Z","iopub.execute_input":"2023-09-19T13:22:28.151508Z","iopub.status.idle":"2023-09-19T13:22:28.168331Z","shell.execute_reply.started":"2023-09-19T13:22:28.151476Z","shell.execute_reply":"2023-09-19T13:22:28.166996Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The non-empty columns are in between. These are the indexes for which we may have at least a value:","metadata":{}},{"cell_type":"code","source":"(~train[reactivity_cols].isnull()).any().values.nonzero()[0]","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.172675Z","iopub.execute_input":"2023-09-19T13:22:28.173095Z","iopub.status.idle":"2023-09-19T13:22:28.185242Z","shell.execute_reply.started":"2023-09-19T13:22:28.173049Z","shell.execute_reply":"2023-09-19T13:22:28.184075Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"If the corresponding column is all null `train[_].isnull().mean()` should be 1. We could do the following or just attaching `reactivity_` at the begining of each index above plus one.","metadata":{}},{"cell_type":"code","source":"non_null_targets = [\n    _ for _ in reactivity_cols\n    if train[_].isnull().mean() != 1\n]","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.187232Z","iopub.execute_input":"2023-09-19T13:22:28.187579Z","iopub.status.idle":"2023-09-19T13:22:28.239650Z","shell.execute_reply.started":"2023-09-19T13:22:28.187545Z","shell.execute_reply":"2023-09-19T13:22:28.238341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now we can describe these values only to see the obvius - their values differ accross diferent RNA sequences (as expected).","metadata":{}},{"cell_type":"code","source":"train[non_null_targets].describe()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.241215Z","iopub.execute_input":"2023-09-19T13:22:28.241572Z","iopub.status.idle":"2023-09-19T13:22:28.474476Z","shell.execute_reply.started":"2023-09-19T13:22:28.241542Z","shell.execute_reply":"2023-09-19T13:22:28.473264Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Descriptive are very useful while comparing two distributions:","metadata":{}},{"cell_type":"code","source":"train_rnas = pd.read_csv(train_file, usecols = ['sequence'])['sequence']\ntest_rnas = test['sequence']\n\ntrain_rna_lens = train_rnas.apply(len)\ntest_rna_lens = test_rnas.apply(len)\n\nprint(f'''== RNA lenghts ==\n\nTrain\n___\nMax: {train_rna_lens.max()}\nMin: {train_rna_lens.min()}\nMean: {train_rna_lens.mean():.2f}\nMedian: {train_rna_lens.median()}\nStd: {train_rna_lens.std():.2f}\n\nTest\n___\nMax: {test_rna_lens.max()}\nMin: {test_rna_lens.min()}\nMean: {test_rna_lens.mean():.2f}\nMedian: {test_rna_lens.median()}\nStd: {test_rna_lens.std():.2f}\n''')","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:22:28.476286Z","iopub.execute_input":"2023-09-19T13:22:28.476611Z","iopub.status.idle":"2023-09-19T13:23:20.374179Z","shell.execute_reply.started":"2023-09-19T13:22:28.476583Z","shell.execute_reply":"2023-09-19T13:23:20.373066Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"They clearly differ - here is what called my attetion: the smaller length for the test sample is the median value for the training sample, and the bigger is way bigger than any length from trainig. \n\nPlotting both distributions.","metadata":{}},{"cell_type":"code","source":"sns.kdeplot(train_rna_lens, label='train', alpha=.3)\nsns.kdeplot(test_rna_lens, label='test', alpha=.3)\nsns.despine()\nplt.xlabel('RNA lenghts')\nplt.legend()","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:20.375387Z","iopub.execute_input":"2023-09-19T13:23:20.375684Z","iopub.status.idle":"2023-09-19T13:23:35.513607Z","shell.execute_reply.started":"2023-09-19T13:23:20.375658Z","shell.execute_reply":"2023-09-19T13:23:35.512447Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Again, which are the possible nucleobases?","metadata":{}},{"cell_type":"code","source":"avaliable_bases = np.unique(list(train_rnas[0]))\nprint(avaliable_bases)","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:35.514919Z","iopub.execute_input":"2023-09-19T13:23:35.515237Z","iopub.status.idle":"2023-09-19T13:23:35.521149Z","shell.execute_reply.started":"2023-09-19T13:23:35.515208Z","shell.execute_reply":"2023-09-19T13:23:35.519936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looking at a complete RNA sequence:","metadata":{}},{"cell_type":"code","source":"train_rnas[0]","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:35.522692Z","iopub.execute_input":"2023-09-19T13:23:35.523060Z","iopub.status.idle":"2023-09-19T13:23:35.539580Z","shell.execute_reply.started":"2023-09-19T13:23:35.522999Z","shell.execute_reply":"2023-09-19T13:23:35.538472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The following function receives an RNA sequence plus the targeted nucleobase and returns a tuple with the actual positions for it in the sequence and their relative positions - that's the position controlled by its size:","metadata":{}},{"cell_type":"code","source":"def get_base_distribution(rna_sequence, nitrogen_base='A'):\n    positions = [\n        current_base_index for current_base_index,current_base \n        in enumerate(rna_sequence) if current_base == nitrogen_base\n    ]\n    positions = np.array(positions)\n    return (positions, positions / len(rna_sequence) )\n\nget_base_distribution(rna_sequence = train_rnas[0],\n                      nitrogen_base = 'A')","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:35.540974Z","iopub.execute_input":"2023-09-19T13:23:35.541380Z","iopub.status.idle":"2023-09-19T13:23:35.555142Z","shell.execute_reply.started":"2023-09-19T13:23:35.541350Z","shell.execute_reply":"2023-09-19T13:23:35.553873Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Visualizing a sequence:","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(16,2))\nplt.axis('off')\nplt.title(train_rnas[0], size=7)\n\nax = plt.subplot(111)\n\nfor base in avaliable_bases:\n    x = get_base_distribution(rna_sequence = train_rnas[0],\n                              nitrogen_base = base)[0]\n    y = np.repeat(1, len(x))\n    plt.bar(x=x, height=y, label=base, width=1)\n\nsns.despine()\n\nlegend = plt.legend()\nlegend.get_frame().set_facecolor('none')\n\nbox = ax.get_position()\nax.set_position([box.x0, box.y0 + box.height * 0.1,\n                 box.width, box.height * 0.9])\nax.legend(loc='upper center',\n          bbox_to_anchor=(0.5, -0.05),\n          ncol=5,\n          frameon=False\n         )","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:35.556531Z","iopub.execute_input":"2023-09-19T13:23:35.557158Z","iopub.status.idle":"2023-09-19T13:23:36.129060Z","shell.execute_reply.started":"2023-09-19T13:23:35.557118Z","shell.execute_reply":"2023-09-19T13:23:36.127884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Creating a matrix where each row is a sequence and the columns represents the nucleobases:","metadata":{}},{"cell_type":"code","source":"paading = train_rna_lens.max()\nn_rnas = 500\nrna_map = np.zeros((len(train_rnas[:n_rnas]), paading))\nrna_color = {\n    'A' : 1,\n    'C' : 2,\n    'G' : 3,\n    'U' : 4\n}\n\nfor r, rna_seq in enumerate(train_rnas[:n_rnas]):\n    rna_map[r] = np.array([rna_color[_] for _ in list(rna_seq)] + [0 for __ in range(paading - len(rna_seq))])\n\nrna_map","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:36.130630Z","iopub.execute_input":"2023-09-19T13:23:36.130993Z","iopub.status.idle":"2023-09-19T13:23:36.174176Z","shell.execute_reply.started":"2023-09-19T13:23:36.130961Z","shell.execute_reply":"2023-09-19T13:23:36.173083Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Checking if it's correct:","metadata":{}},{"cell_type":"code","source":"print(rna_map[9])\nprint(train_rnas[9])","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:36.175646Z","iopub.execute_input":"2023-09-19T13:23:36.176317Z","iopub.status.idle":"2023-09-19T13:23:36.185271Z","shell.execute_reply.started":"2023-09-19T13:23:36.176275Z","shell.execute_reply":"2023-09-19T13:23:36.183926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Plotting 500 RNAs (TRAINING Sample)\n\n### first 500","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize=(16,4))\n\nim = plt.imshow(rna_map, aspect='auto', cmap=CMAP, vmin = 0, vmax = 4, interpolation='nearest')\n\nrna_values = {\n    1 : 'A',\n    2 : 'C',\n    3 : 'G',\n    4 : 'U',\n    0 : 'Pad'\n}\n\ncolor_values = [1, 2, 3, 4, 0]\ncolors = [ im.cmap(im.norm(color_value)) for color_value in color_values]\npatches = [ mpatches.Patch(color=colors[i], label=rna_values[color_value]) for i,color_value in enumerate(color_values) ]\nplt.legend(handles=patches, bbox_to_anchor=(0.65, -0.1), ncol=5, borderaxespad=0. )","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:36.187112Z","iopub.execute_input":"2023-09-19T13:23:36.187510Z","iopub.status.idle":"2023-09-19T13:23:36.731383Z","shell.execute_reply.started":"2023-09-19T13:23:36.187469Z","shell.execute_reply":"2023-09-19T13:23:36.730073Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For these first 125 or so _nucleobases_ from these 500 sequences we can see that it differs very little from each other. The ones in between \\#125 - \\#150 seens to differ more. We're looking into the first 500, which is not that much if you think about the remaining 1 Million. I suspect that there's nothing to it, might be the way rows were organized. \n\n### last 500","metadata":{}},{"cell_type":"code","source":"rna_map = np.zeros((len(train_rnas[-n_rnas:]), paading))\n\nfor r, rna_seq in enumerate(train_rnas[-n_rnas:]):\n    rna_map[r] = np.array([rna_color[_] for _ in list(rna_seq)] + [0 for __ in range(paading - len(rna_seq))])\n\nplt.figure(figsize=(16,4))\nim = plt.imshow(rna_map, aspect='auto', cmap=CMAP, vmin = 0, vmax = 4, interpolation='nearest')\nplt.legend(handles=patches, bbox_to_anchor=(0.65, -0.1), ncol=5, borderaxespad=0. )","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:36.733408Z","iopub.execute_input":"2023-09-19T13:23:36.734100Z","iopub.status.idle":"2023-09-19T13:23:37.301519Z","shell.execute_reply.started":"2023-09-19T13:23:36.734065Z","shell.execute_reply":"2023-09-19T13:23:37.300267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For the last 500 training sequences they very much differ from the first 500 but are similar accross one another.\n\n## Plotting 500 RNAs (TEST Sample)\n\n### first 500","metadata":{}},{"cell_type":"code","source":"paading = test_rna_lens.max()\nrna_map = np.zeros((len(test_rnas[:n_rnas]), paading))\n\nfor r, rna_seq in enumerate(test_rnas[:n_rnas]):\n    rna_map[r] = np.array([rna_color[_] for _ in list(rna_seq)] + [0 for __ in range(paading - len(rna_seq))])\n\n\nplt.figure(figsize=(16,4))\nim = plt.imshow(rna_map, aspect='auto', cmap=CMAP, vmin = 0, vmax = 4, interpolation='nearest')\nplt.legend(handles=patches, bbox_to_anchor=(0.65, -0.1), ncol=5, borderaxespad=0. )","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:37.303437Z","iopub.execute_input":"2023-09-19T13:23:37.304239Z","iopub.status.idle":"2023-09-19T13:23:37.997775Z","shell.execute_reply.started":"2023-09-19T13:23:37.304201Z","shell.execute_reply":"2023-09-19T13:23:37.996742Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"For the first 500 test sequences it's hard to see it because the padding is too big.\n\n### last 500","metadata":{}},{"cell_type":"code","source":"rna_map = np.zeros((len(test_rnas[-n_rnas:]), paading))\n\nfor r, rna_seq in enumerate(test_rnas[-n_rnas:]):\n    rna_map[r] = np.array([rna_color[_] for _ in list(rna_seq)] + [0 for __ in range(paading - len(rna_seq))])\n\n\nplt.figure(figsize=(16,4))\nim = plt.imshow(rna_map, aspect='auto', cmap=CMAP, vmin = 0, vmax = 4, interpolation='nearest')\nplt.legend(handles=patches, bbox_to_anchor=(0.65, -0.1), ncol=5, borderaxespad=0. )","metadata":{"execution":{"iopub.status.busy":"2023-09-19T13:23:37.999603Z","iopub.execute_input":"2023-09-19T13:23:37.999933Z","iopub.status.idle":"2023-09-19T13:23:38.575765Z","shell.execute_reply.started":"2023-09-19T13:23:37.999904Z","shell.execute_reply":"2023-09-19T13:23:38.574588Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"I've liked to see the patterns that showed up, hope you've too. Let me know how to improve my analysis and what you think the implications should be for this competition.","metadata":{}}]}