{"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"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":51294,"databundleVersionId":6923401,"sourceType":"competition"}],"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Testing different sequence encodings","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd\nimport matplotlib.pyplot as plt\n\nfrom sklearn.model_selection import KFold\nfrom sklearn.ensemble import RandomForestRegressor\nfrom sklearn.model_selection import train_test_split\n\ndef MakeSequenceEncoding(Sequence,tokenDict=None):\n    \n    stringFrags = [val for val in Sequence]\n    nToAdd = (206) - len(stringFrags)\n    \n    encoded = [tokenDict[val] for val in stringFrags] \n\n    toaddvec = [0 for _ in range(len(encoded[0]))]\n    toAdd = [toaddvec for k in range(nToAdd)]\n    \n    encoded = np.array(encoded + toAdd).ravel()\n    \n    return encoded","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Loading the data ","metadata":{}},{"cell_type":"code","source":"\ndata = pd.read_csv('/kaggle/input/stanford-ribonanza-rna-folding/train_data_QUICK_START.csv')\ndata = data.fillna(0)\n\ndataseqs = data.groupby('sequence_id')['sequence'].unique().to_frame()\ndataseqs['sequence'] = [val[0] for val in dataseqs['sequence']]\n\ncats = ['2A3_MaP', 'DMS_MaP']\n\nexp0 = data[data['experiment_type']==cats[0]]\nexp0 = exp0.set_index('sequence_id')\n\nexp1 = data[data['experiment_type']==cats[1]]\nexp1 = exp1.set_index('sequence_id')\n\ndatacols = exp0.columns[3:209]\n","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data splitting ","metadata":{}},{"cell_type":"code","source":"X_train, X_test, _, _ = train_test_split(dataseqs.index,dataseqs.index,test_size=0.1, random_state=42)\n\nX_train = X_train[0:50000]\nX_test = X_test[0:25000]\n\naydata0 = exp0[datacols].loc[X_train].values\naydata1= exp1[datacols].loc[X_train].values\n\nayval0 = exp0[datacols].loc[X_test].values\nayval1 = exp1[datacols].loc[X_test].values\n\ngroup_y = [aydata0,ayval0,aydata1,ayval1]\n\ngroup_x = [dataseqs['sequence'].loc[X_train].values, dataseqs['sequence'].loc[X_test].values]","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Testing different sequence encodings","metadata":{}},{"cell_type":"code","source":"def TestEncoding(data_x,data_y,encoderdict):\n    \n    dst_xtr,dst_xval = data_x\n    ydata0,yval0,ydata1,yval1 = data_y\n    \n    xdata = [MakeSequenceEncoding(val,tokenDict=encoderdict) for val in dst_xtr]\n    xdata = np.stack(xdata)\n    \n    xval = [MakeSequenceEncoding(val,tokenDict=encoderdict) for val in dst_xval]\n    xval = np.stack(xval)\n    \n    kf = KFold(n_splits=6,random_state=1,shuffle=True)\n    \n    fig,axs = plt.subplots(6,2,figsize=(15,15))\n    \n    for i, (train_index, test_index) in enumerate(kf.split(xdata)):\n        \n        xtrain,xtest = xdata[train_index],xdata[test_index]\n        ytrain0,ytest0 = ydata0[train_index],ydata0[test_index]\n        ytrain1,ytest1 = ydata1[train_index],ydata1[test_index]\n        \n        model0 = RandomForestRegressor(50,random_state=0,n_jobs=-1)\n        model0.fit(xtrain, ytrain0)\n        preds0 = model0.predict(xval)\n        preds0f = model0.predict(xtest)\n        \n        model1 = RandomForestRegressor(50,random_state=0,n_jobs=-1)\n        model1.fit(xtrain, ytrain1)\n        preds1 = model1.predict(xval)\n        preds1f = model1.predict(xtest)\n        \n        serror0 = np.mean((yval0.ravel()-preds0.ravel())**2)\n        serror0f = np.mean((ytest0.ravel()-preds0f.ravel())**2)\n        serror1 = np.mean((yval1.ravel()-preds1.ravel())**2)\n        serror1f = np.mean((ytest1.ravel()-preds1f.ravel())**2)\n        \n        error0 = np.mean(np.abs(yval0.ravel()-preds0.ravel()))\n        error0f = np.mean(np.abs(ytest0.ravel()-preds0f.ravel()))\n        error1 = np.mean(np.abs(yval1.ravel()-preds1.ravel()))\n        error1f = np.mean(np.abs(ytest1.ravel()-preds1f.ravel()))\n        \n        axs[i,0].scatter(yval0.ravel(),preds0.ravel(),alpha=0.1)\n        axs[i,0].text(0.05,0.85,'val MSE = '+str(round(serror0,3)),transform = axs[i,0].transAxes)\n        axs[i,0].text(0.05,0.65,'val MAE = '+str(round(error0,3)),transform = axs[i,0].transAxes)\n        axs[i,0].text(0.05,0.45,'fold MSE = '+str(round(serror0f,3)),transform = axs[i,0].transAxes)\n        axs[i,0].text(0.05,0.25,'fold MAE = '+str(round(error0f,3)),transform = axs[i,0].transAxes)\n        axs[i,0].set_title('Fold = '+str(i))\n        \n        axs[i,1].scatter(yval1.ravel(),preds1.ravel(),alpha=0.1)\n        axs[i,1].text(0.05,0.85,'val MSE = '+str(round(serror1,3)),transform = axs[i,1].transAxes)\n        axs[i,1].text(0.05,0.65,'val MAE = '+str(round(error1,3)),transform = axs[i,1].transAxes)\n        axs[i,1].text(0.05,0.45,'fold MSE = '+str(round(serror1f,3)),transform = axs[i,1].transAxes)\n        axs[i,1].text(0.05,0.25,'fold MAE = '+str(round(error1f,3)),transform = axs[i,1].transAxes)\n        axs[i,1].set_title('Fold = '+str(i))\n    plt.tight_layout()","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true}},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## One hot encoded sequences","metadata":{}},{"cell_type":"code","source":"Alphabet = ['A','C','U','G']\n\nTokenDictionary = {}\n\nfor k,val in enumerate(Alphabet):\n    currentVec = [0 for j in range(len(Alphabet))]\n    currentVec[k] = 1\n    TokenDictionary[val]=currentVec\n\nTestEncoding(group_x,group_y,TokenDictionary)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Basic chemical properties\n\nFirst two element lists encodes for chemical structure(purines or pyrimidines), second one number bonds","metadata":{}},{"cell_type":"code","source":"TokenDictionary2 = {}\n\nTokenDictionary2['A'] = [1,0] + [0,1]\nTokenDictionary2['C'] = [0,1] + [1,0]\nTokenDictionary2['U'] = [0,1] + [0,1]\nTokenDictionary2['G'] = [1,0] + [1,0]\n\nTestEncoding(group_x,group_y,TokenDictionary2)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Nucleotide pka\npKa constants taken from https://pubs.acs.org/doi/10.1021/acs.jpca.1c10905","metadata":{}},{"cell_type":"code","source":"TokenDictionary3 = {}\n\nTokenDictionary3['A'] = [3.5]\nTokenDictionary3['C'] = [4.2]\nTokenDictionary3['U'] = [9.25]\nTokenDictionary3['G'] = [9.2]\n\nTestEncoding(group_x,group_y,TokenDictionary3)","metadata":{},"execution_count":null,"outputs":[]}]}