{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Ribonanza with Kmers","metadata":{}},{"cell_type":"code","source":"import numpy as np \nimport pandas as pd \n\nfrom itertools import product\n\nAlphabet = ['A','C','G','U']\n\nBlocks = []\n\nmaxSize = 5\nfor k in range(1,maxSize):    \n    Blocks.append([''.join(i) for i in product(Alphabet, repeat = k)])\n\ndef SplitString(String,ChunkSize):\n    return [String[k:k+ChunkSize] for k in range(len(String)-ChunkSize+1)]\n\ndef DictionaryLibrary(ElementsBlocks):\n    \n    maindict = {}\n    \n    for k,val in enumerate(ElementsBlocks):\n        innerDict = {}\n        for kk,sal in enumerate(val):\n            innerDict[sal] = kk\n        \n        maindict[k+1] = innerDict\n    \n    return maindict\n\nDicts = DictionaryLibrary(Blocks)\n\ndef CountUniqueElements(String,dictionary):\n\n    localCounter = [0 for k in range(dictionary.__len__())]\n    ProcessedString = SplitString(String,len(list(dictionary.keys())[0]))    \n    \n    for val in ProcessedString:\n        \n        if val in dictionary.keys():\n            \n            localPosition=dictionary[val]\n            localCounter[localPosition]=localCounter[localPosition]+1\n            \n    localCounter=[val/len(ProcessedString) for val in localCounter]\n        \n    return localCounter","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T00:20:20.617581Z","iopub.execute_input":"2023-09-14T00:20:20.617968Z","iopub.status.idle":"2023-09-14T00:20:21.023592Z","shell.execute_reply.started":"2023-09-14T00:20:20.617936Z","shell.execute_reply":"2023-09-14T00:20:21.022607Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Making the dataset","metadata":{}},{"cell_type":"code","source":"data = pd.read_csv('/kaggle/input/stanford-ribonanza-rna-folding/train_data.csv')\ndata = data.query('signal_to_noise < 0.7 and signal_to_noise > 0')\ndata = data.fillna(0)\n\nuniquedata = data.groupby('sequence_id')['sequence'].value_counts().index\n\ncontainer = []\n\nfor idx in uniquedata:\n    seq = idx[1]\n    vec = []\n    for kk in range(len(Dicts)):\n        minivec = CountUniqueElements(seq, Dicts[kk+1])\n        minv,maxv = np.min(minivec),np.max(minivec)\n        normvec = [(val-minv)/(maxv-minv) for val in minivec]\n        vec = vec + normvec\n    container.append(vec)\n\ncontainer = np.array(container).astype(np.float16)\n\nKmerDF = pd.DataFrame()\nheaders = [val for li in Blocks for val in li]\nKmerDF = pd.DataFrame(container,columns=headers)\nKmerDF['id'] = [val[0] for val in uniquedata]\nKmerDF = KmerDF.set_index('id')","metadata":{"execution":{"iopub.status.busy":"2023-09-14T00:20:21.025332Z","iopub.execute_input":"2023-09-14T00:20:21.026345Z","iopub.status.idle":"2023-09-14T00:28:06.602646Z","shell.execute_reply.started":"2023-09-14T00:20:21.026306Z","shell.execute_reply":"2023-09-14T00:28:06.601345Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# MLP","metadata":{}},{"cell_type":"code","source":"import time\nfrom sklearn import preprocessing as pr\n\nfrom typing import Sequence\n\nimport jax\nimport optax\nimport jax.numpy as jnp\nfrom jax import random\nfrom flax import linen as nn\n\nclass Coder(nn.Module):\n    \n    Units: Sequence[int]\n    Name: str \n    train: bool = True \n    \n    @nn.compact\n    def __call__(self,inputs):\n        x = inputs\n        for k,feat in enumerate(self.Units):\n            x = nn.Dense(feat,use_bias=False,name = self.Name+' layer_'+str(k))(x)\n            x = nn.BatchNorm(use_running_average=not self.train,name = self.Name+' norm_'+str(k))(x)\n            x = nn.leaky_relu(x)\n        return x\n\nclass Regressor(nn.Module):\n    \n    Units: Sequence[int]\n    train: bool = True \n    \n    def setup(self):\n        self.Body = Coder(self.Units,'body',train=self.train)\n        self.out = nn.Dense(self.Units[-1],name='regressorOut')\n        self.norm = nn.BatchNorm(use_running_average=not self.train,name='outnorm')\n    \n    @nn.compact\n    def __call__(self,inputs):\n        \n        x = self.Body(inputs)\n        x = self.out(x)\n        x = self.norm(x)\n        \n        return x\n\ndef MainLoss(Model,params,batchStats,batch,target):\n  \n    block, newbatchst = Model().apply({'params': params, 'batch_stats': batchStats}, batch,mutable=['batch_stats'])\n    y_target = block\n    \n    loss_value = optax.l2_loss(y_target, target).mean()\n    total_loss = loss_value\n    \n    return total_loss,newbatchst['batch_stats']\n\ndef TrainModel(TrainDataSet,TestDataSet,Loss,params,batchStats,rng,epochs=10,batch_size=64,lr=0.005):\n    \n    TrainData,TrainTargets = TrainDataSet\n    TestData, TestTargets = TestDataSet\n    \n    totalSteps = epochs*(TrainData.shape[0]//batch_size) + epochs\n    stepsPerCycle = totalSteps//4\n    \n    esp = [{\"init_value\":lr/10, \n            \"peak_value\":(lr)/((k+1)), \n            \"decay_steps\":int(stepsPerCycle*0.75), \n            \"warmup_steps\":int(stepsPerCycle*0.25), \n            \"end_value\":lr/10} for k in range(8)]\n    \n    Scheduler = optax.sgdr_schedule(esp)\n    localOptimizer = optax.adam(learning_rate=Scheduler)\n    optState = localOptimizer.init(params)\n    \n    @jax.jit\n    def step(params,batchStats ,optState, batch,target):\n        \n        (loss_value,batchStats), grads = jax.value_and_grad(Loss,has_aux=True)(params,batchStats, batch,target)\n        updates, optState = localOptimizer.update(grads, optState, params)\n        params = optax.apply_updates(params, updates)\n        \n        return params,batchStats, optState, loss_value\n    \n    @jax.jit\n    def getloss(params,batchStats, batch,target):\n        (loss_value,_), _ = jax.value_and_grad(Loss,has_aux=True)(params,batchStats,batch,target)\n        return loss_value\n    \n    trainloss = []\n    testloss = []\n    \n    for epoch in range(epochs):\n        \n        st = time.time()\n        batchtime = []\n        losses = []\n        \n        for k in range(0,TrainData.shape[0],batch_size):\n    \n            stb = time.time()\n            batch = TrainData[k:k+batch_size]\n            b_targets = TrainTargets[k:k+batch_size]\n            \n            params,batchStats ,optState, lossval = step(params,batchStats,optState,batch,b_targets)\n            losses.append(lossval)\n            batchtime.append(time.time()-stb)\n        \n        valloss = []\n        for i in range(0,TestData.shape[0],batch_size):\n            val_batch = TestData[i:i+batch_size]\n            val_targets = TestTargets[i:i+batch_size]\n            valloss.append(getloss(params,batchStats,val_batch,val_targets))\n        \n        mbatch = 1000*np.mean(batchtime)\n        meanloss = np.mean(losses)\n        meanvalloss = np.mean(valloss)\n        \n        trainloss.append(meanloss)\n        testloss.append(meanvalloss)\n        \n        localIndex = np.arange(len(TrainData))\n        np.random.shuffle(localIndex)\n        \n        TrainData = TrainData[localIndex]\n        TrainTargets = TrainTargets[localIndex]\n    \n        end = time.time()\n        output = 'Epoch = '+str(epoch) + ' Time per epoch = ' + str(round(end-st,3)) + 's  Time per batch = ' + str(round(mbatch,3)) + 'ms' + ' Train Loss = ' + str(meanloss) +' Test Loss = ' + str(meanvalloss)\n        print(output)\n        \n    return trainloss,testloss,params,batchStats\n","metadata":{"_kg_hide-input":true,"jupyter":{"source_hidden":true},"execution":{"iopub.status.busy":"2023-09-14T00:28:06.604417Z","iopub.execute_input":"2023-09-14T00:28:06.604800Z","iopub.status.idle":"2023-09-14T00:28:09.099490Z","shell.execute_reply.started":"2023-09-14T00:28:06.604769Z","shell.execute_reply":"2023-09-14T00:28:09.098177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\nfrom sklearn.metrics import mean_absolute_error\nfrom sklearn.model_selection import train_test_split\n\nmainUnits  = [340,300,250,170]\n\nbatchSize = 1024\nInputShape = 340\n\nepochs = 15\ncolumns = data.columns[7:177]","metadata":{"execution":{"iopub.status.busy":"2023-09-14T00:28:09.102172Z","iopub.execute_input":"2023-09-14T00:28:09.102874Z","iopub.status.idle":"2023-09-14T00:28:09.230850Z","shell.execute_reply.started":"2023-09-14T00:28:09.102830Z","shell.execute_reply":"2023-09-14T00:28:09.229753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for exps in data['experiment_type'].unique():\n    \n    exp0 = data[data['experiment_type']==exps]\n    \n    Inxtrain,Inxtest, _, _ = train_test_split(exp0['sequence_id'], exp0['sequence_id'], test_size=0.10, random_state=42)\n    \n    Inxtrain = Inxtrain[0:batchSize*(len(Inxtrain)//batchSize)]\n    Inxtest = Inxtest[0:batchSize*(len(Inxtest)//batchSize)]\n    \n    trainIn = exp0[exp0['sequence_id'].isin(Inxtrain)]\n    testIn = exp0[exp0['sequence_id'].isin(Inxtest)]\n    \n    trainTarget = trainIn[columns].values\n    testTarget = testIn[columns].values\n    \n    trainData = KmerDF.loc[trainIn['sequence_id']].values\n    testData = KmerDF.loc[testIn['sequence_id']].values\n\n    scaler = pr.MinMaxScaler()\n    scaler.fit(trainData)\n    \n    scalerTarget = pr.MinMaxScaler()\n    scalerTarget.fit(trainTarget)\n    \n    trainDta = scaler.transform(trainData)\n    testDta = scaler.transform(testData)\n    \n    trainTarg = scalerTarget.transform(trainTarget)\n    testTarg = scalerTarget.transform(testTarget)\n    \n    np.random.seed(128)\n    \n    def RegressorModel():\n        return Regressor(mainUnits)\n    \n    def loss(params,batchStats,batch,target):\n        return MainLoss(RegressorModel,params,batchStats,batch,target)\n    \n    rng = random.PRNGKey(0)\n    rng, key = random.split(rng)\n    \n    init_data = jnp.ones((batchSize, InputShape), jnp.float32)\n    initModel = RegressorModel().init(key, init_data)\n    \n    params0 = initModel['params']\n    batchStats = initModel['batch_stats']\n    \n    trloss,tstloss,params0,batchStats = TrainModel([trainDta,trainTarg],[testDta,testTarg],loss,params0,\n                                        batchStats,rng,lr=0.001,epochs=epochs,\n                                        batch_size=batchSize)\n    \n    state = {'params':params0,'batch_stats':batchStats}\n    def RegressorModelPreds(batch):\n        return Regressor(mainUnits,train=False).apply(state,batch)\n    \n    preds = RegressorModelPreds(testData)\n    preds = scalerTarget.inverse_transform(preds)\n    \n    plt.figure()\n    plt.scatter(testTarget.ravel(),preds.ravel(),alpha=0.25)\n    ax = plt.gca()\n    ax.set_xlabel('test')\n    ax.set_ylabel('predictions')\n    print(mean_absolute_error(testTarget.ravel(),preds.ravel()))\n","metadata":{"execution":{"iopub.status.busy":"2023-09-14T00:28:09.232489Z","iopub.execute_input":"2023-09-14T00:28:09.232921Z","iopub.status.idle":"2023-09-14T00:38:03.052404Z","shell.execute_reply.started":"2023-09-14T00:28:09.232880Z","shell.execute_reply":"2023-09-14T00:38:03.051070Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Random Forest regressor ","metadata":{}},{"cell_type":"code","source":"from sklearn.ensemble import RandomForestRegressor\n\nfor exps in data['experiment_type'].unique():\n    \n    exp0 = data[data['experiment_type']==exps]\n    \n    Inxtrain,Inxtest, _, _ = train_test_split(exp0['sequence_id'], exp0['sequence_id'], test_size=0.10, random_state=42)\n    \n    Inxtrain = Inxtrain[0:batchSize*(len(Inxtrain)//batchSize)]\n    Inxtest = Inxtest[0:batchSize*(len(Inxtest)//batchSize)]\n    \n    trainIn = exp0[exp0['sequence_id'].isin(Inxtrain)]\n    testIn = exp0[exp0['sequence_id'].isin(Inxtest)]\n    \n    trainTarget = trainIn[columns].values\n    testTarget = testIn[columns].values\n    \n    trainData = KmerDF.loc[trainIn['sequence_id']].values\n    testData = KmerDF.loc[testIn['sequence_id']].values\n\n    scaler = pr.MinMaxScaler()\n    scaler.fit(trainData)\n    \n    scalerTarget = pr.MinMaxScaler()\n    scalerTarget.fit(trainTarget)\n    \n    trainDta = scaler.transform(trainData)\n    testDta = scaler.transform(testData)\n    \n    trainTarg = scalerTarget.transform(trainTarget)\n    testTarg = scalerTarget.transform(testTarget)\n    \n    model = RandomForestRegressor(n_estimators=5,n_jobs=-1)\n    model.fit(trainDta,trainTarg)\n    \n    preds = model.predict(testDta)\n    \n    preds = scalerTarget.inverse_transform(preds)\n    \n    plt.figure()\n    plt.scatter(testTarget.ravel(),preds.ravel(),alpha=0.25)\n    ax = plt.gca()\n    ax.set_xlabel('test')\n    ax.set_ylabel('predictions')\n    print(mean_absolute_error(testTarget.ravel(),preds.ravel()))\n","metadata":{"execution":{"iopub.status.busy":"2023-09-14T00:38:03.054075Z","iopub.execute_input":"2023-09-14T00:38:03.055329Z","iopub.status.idle":"2023-09-14T00:39:36.667552Z","shell.execute_reply.started":"2023-09-14T00:38:03.055276Z","shell.execute_reply":"2023-09-14T00:39:36.665571Z"},"trusted":true},"execution_count":null,"outputs":[]}],"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"}}