{"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":"Coming at this project from the chemical side, it seems obvious that the chemical nature of the RNA molecule should have a big effect on which sites are reactive. I see previous notebooks that computed baselines using medians for each experiment across all data, as well as medians by position number in the RNA sequence. Unsurprisingly, position number makes basically no difference, since the reactivity is only available for RNA nucleotides reasonably far from the ends. (https://www.kaggle.com/code/alexandervc/ribonanza-baselines-medians-onehot-ridge-lb0-229). I computed means for each experiment type and got very similar error values (see table below)\n\nLooking at the actual experiments performed, the 2A3 experiment reacts with the RNA backbone, which is the same for all RNA nucleotides (A, C, G, and U). The DMS experiment reacts with the nucleotide itself. Traditional DMS experiments more or less only reacted with A and C, but the modified current procedure at least gives some reactivity for all 4 bases. Based on this, I expect the nucleotide type to have a big effect on the reactivity.\n\nThis notebook computes the mean reactivity for each nucleotide for each of the two experiment type (8 total values), and assigns each value in the submitted table to one of those 8 values. The results summarized in the table below show that just knowing which nucleotide you're looking at is enough to have a fairly big effect on the model accuracy. It's definitely not enough to be competitive with the top models, but it's a useful baseline to have.\n\n| Submission Type | Score |\n| --- | --- |\nSubmit all zeros | 0.33975\nMedians for each experiment type (DMS=0.123, 2A3=0.214) | 0.28255\nMedians for each nucleotide position | 0.28208\nMeans for each experiment type (DMS=0.334, 2A3=0.348) | 0.29931\nMeans for each nucleotide+experiment type | 0.23305\n","metadata":{}},{"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-10-09T16:48:53.138236Z","iopub.execute_input":"2023-10-09T16:48:53.138591Z","iopub.status.idle":"2023-10-09T16:48:53.503215Z","shell.execute_reply.started":"2023-10-09T16:48:53.138561Z","shell.execute_reply":"2023-10-09T16:48:53.502242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the training data\nfn = '/kaggle/input/stanford-ribonanza-rna-folding-converted/train_data.parquet'\ndf = pd.read_parquet(fn)\nprint(df.shape)","metadata":{"execution":{"iopub.status.busy":"2023-10-09T16:48:53.505150Z","iopub.execute_input":"2023-10-09T16:48:53.505674Z","iopub.status.idle":"2023-10-09T16:49:04.415019Z","shell.execute_reply.started":"2023-10-09T16:48:53.505638Z","shell.execute_reply":"2023-10-09T16:49:04.413981Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check for and drop duplicate rows\n# To count as a duplicate, both the sequence and the experiment type must be the same\ndf = df.drop_duplicates(subset=['sequence', 'experiment_type'])\nprint(df.shape)","metadata":{"execution":{"iopub.status.busy":"2023-10-09T16:49:04.416184Z","iopub.execute_input":"2023-10-09T16:49:04.416787Z","iopub.status.idle":"2023-10-09T16:49:07.376631Z","shell.execute_reply.started":"2023-10-09T16:49:04.416761Z","shell.execute_reply":"2023-10-09T16:49:07.375585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Add length data to the dataframe\ndf['len'] = [len(t) for t in df['sequence']]\ndf['len'].value_counts()","metadata":{"execution":{"iopub.status.busy":"2023-10-09T16:49:07.378907Z","iopub.execute_input":"2023-10-09T16:49:07.379925Z","iopub.status.idle":"2023-10-09T16:49:08.005410Z","shell.execute_reply.started":"2023-10-09T16:49:07.379886Z","shell.execute_reply":"2023-10-09T16:49:08.004724Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Compute the average values and the counts for each nucleotide, as a function of position\nmeans = pd.DataFrame({'A_DMS': [], 'A_2A3': [], 'A_count': [],\n                          'C_DMS': [], 'C_2A3': [], 'C_count': [],\n                          'G_DMS': [], 'G_2A3': [], 'G_count': [],\n                          'U_DMS': [], 'U_2A3': [], 'U_count': []})\n\nfor i in range(206):\n    if i%10 == 9: print(i+1)\n    row = {}\n    counts = df[(df['experiment_type'] == 'DMS_MaP') \n                & (~df['reactivity_'+str(i+1).zfill(4)].isna())]['sequence'].str[i].value_counts()\n    for b, val in counts.items():\n        row[b+'_DMS'] = df.loc[(df['experiment_type'] == 'DMS_MaP') & (df['sequence'].str[i] == b), \n                            'reactivity_'+str(i+1).zfill(4)].mean()\n        row[b+'_2A3'] = df.loc[(df['experiment_type'] == '2A3_MaP') & (df['sequence'].str[i] == b), \n                            'reactivity_'+str(i+1).zfill(4)].mean()\n        row[b+'_count'] = val\n    means.loc[len(means.index)] = row","metadata":{"execution":{"iopub.status.busy":"2023-10-09T16:49:09.183770Z","iopub.execute_input":"2023-10-09T16:49:09.184536Z","iopub.status.idle":"2023-10-09T17:08:09.606962Z","shell.execute_reply.started":"2023-10-09T16:49:09.184499Z","shell.execute_reply":"2023-10-09T17:08:09.605720Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# From the counts and averages for each position, compute the overall average\nbases = ['A', 'C', 'G', 'U']\nfor b in bases:\n    means[b+'_DMS_prod'] = means[b+'_DMS'] * means[b+'_count']\n    means[b+'_2A3_prod'] = means[b+'_2A3'] * means[b+'_count']\navg_DMS = {}\navg_2A3 = {}\nprod_DMS, prod_2A3, count = 0, 0, 0\nfor b in bases:\n    avg_DMS[b] = means[b+'_DMS_prod'].sum() / means[b+'_count'].sum()\n    avg_2A3[b] = means[b+'_2A3_prod'].sum() / means[b+'_count'].sum()\n    prod_DMS += means[b+'_DMS_prod'].sum()\n    prod_2A3 += means[b+'_2A3_prod'].sum()\n    count += means[b+'_count'].sum()\navg_DMS['All'] = prod_DMS / count\navg_2A3['All'] = prod_2A3 / count\nprint('DMS: ', avg_DMS)\nprint('2A3: ', avg_2A3)","metadata":{"execution":{"iopub.status.busy":"2023-10-09T17:09:06.165213Z","iopub.execute_input":"2023-10-09T17:09:06.165548Z","iopub.status.idle":"2023-10-09T17:09:06.181398Z","shell.execute_reply.started":"2023-10-09T17:09:06.165522Z","shell.execute_reply":"2023-10-09T17:09:06.180002Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The important thing to note here is that the different nucleotides have very different reactivity. DMS has nearly 20x the reactivity with A as it does with U. The two highly reactive bases are A and C as expected.\n\nThe 2A3 experiment has a smaller range, but still larger than I expected: a 2.6x difference between U (most reactive) and C (least reactive). The two more reactive bases are A and U. That makes sense because A and U are the bases that tend to pair more weakly, and so it's easier to get more flexibility in the RNA backbone at those sites, increasing their reactivity.","metadata":{}},{"cell_type":"code","source":"# Outputs from earlier stored so I don't have to re-run the time-consuming step unless necessary\n'''\navg_DMS = {'A': 0.7289046806161731, \n           'C': 0.4686539396568741, \n           'G': 0.07791616790852551, \n           'U': 0.04017014799018061, \n           'All': 0.3346237541194974}\navg_2A3 = {'A': 0.4098115008601963, \n           'C': 0.19420878305815112, \n           'G': 0.24506126619913013, \n           'U': 0.504049637303388, \n           'All': 0.3478078818875119}\n'''","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read in the test sequences\nfn_test = '/kaggle/input/stanford-ribonanza-rna-folding-converted/test_sequences.parquet'\ndf_test = pd.read_parquet(fn_test)\nprint(df_test.shape)\nfn_sub = '/kaggle/input/stanford-ribonanza-rna-folding/sample_submission.csv'\ndf_sub = pd.read_csv(fn_sub)\ndf_test.head()","metadata":{"execution":{"iopub.status.busy":"2023-10-09T17:09:52.690619Z","iopub.execute_input":"2023-10-09T17:09:52.690950Z","iopub.status.idle":"2023-10-09T17:11:38.650965Z","shell.execute_reply.started":"2023-10-09T17:09:52.690924Z","shell.execute_reply":"2023-10-09T17:11:38.649899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in range(df_test.shape[0]):\n    if (i+1) % 10000 == 0: print(i+1)\n    if df_test['future'].iloc[i] == 0:\n        seq = df_test['sequence'].iloc[i]\n        pred_DMS = list(map(avg_DMS.get, list(seq)))\n        pred_2A3 = list(map(avg_2A3.get, list(seq)))\n        #print(pred_DMS)\n        df_sub['reactivity_DMS_MaP'].iloc[df_test['id_min'].iloc[i]:df_test['id_max'].iloc[i]+1] = pred_DMS\n        df_sub['reactivity_2A3_MaP'].iloc[df_test['id_min'].iloc[i]:df_test['id_max'].iloc[i]+1] = pred_2A3\ndf_sub.to_parquet('predict_means.parquet', index=False)","metadata":{"execution":{"iopub.status.busy":"2023-10-09T17:43:52.857671Z","iopub.execute_input":"2023-10-09T17:43:52.858600Z","iopub.status.idle":"2023-10-09T17:49:06.386134Z","shell.execute_reply.started":"2023-10-09T17:43:52.858557Z","shell.execute_reply":"2023-10-09T17:49:06.384876Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}