{"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":"Uses code from: https://kaggle.com/code/gusthema/cafa-5-protein-function-with-tensorflow/notebook\n\nCross validation code from: https://www.section.io/engineering-education/how-to-implement-k-fold-cross-validation/","metadata":{}},{"cell_type":"code","source":"# UTILITARIES\nimport time\nimport tensorflow as tf\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\n# for progressbar widget\nimport progressbar","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-08-09T07:15:12.989998Z","iopub.execute_input":"2023-08-09T07:15:12.990792Z","iopub.status.idle":"2023-08-09T07:15:12.997579Z","shell.execute_reply.started":"2023-08-09T07:15:12.990752Z","shell.execute_reply":"2023-08-09T07:15:12.996327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# for improved reproducibilty\nimport random as rn\n\nnp.random.seed(1)\nrn.seed(2)\ntf.random.set_seed(3)","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:13.000292Z","iopub.execute_input":"2023-08-09T07:15:13.000681Z","iopub.status.idle":"2023-08-09T07:15:13.011163Z","shell.execute_reply.started":"2023-08-09T07:15:13.000647Z","shell.execute_reply":"2023-08-09T07:15:13.010074Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# This cell outputs the intended submission format\n\nsample_submission = pd.read_csv(\"/kaggle/input/cafa-5-protein-function-prediction/sample_submission.tsv\", sep = \"\\t\", header = None)\nsample_submission.columns = [\"The Protein ID\", \"The Gene Ontology term (GO) ID\", \"Predicted link probability that GO appear in Protein\"]\n\nsample_submission.head(5)","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:13.012349Z","iopub.execute_input":"2023-08-09T07:15:13.012785Z","iopub.status.idle":"2023-08-09T07:15:13.33404Z","shell.execute_reply.started":"2023-08-09T07:15:13.012752Z","shell.execute_reply":"2023-08-09T07:15:13.33284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"MAIN_DIR = \"/kaggle/input/cafa-5-protein-function-prediction\"","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:13.335817Z","iopub.execute_input":"2023-08-09T07:15:13.336297Z","iopub.status.idle":"2023-08-09T07:15:13.342185Z","shell.execute_reply.started":"2023-08-09T07:15:13.336255Z","shell.execute_reply":"2023-08-09T07:15:13.341014Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_sequences_path = MAIN_DIR + \"/Train/train_sequences.fasta\"\ntrain_labels_path = MAIN_DIR + \"/Train/train_terms.tsv\"\ntest_sequences_path = MAIN_DIR + \"/Test (Targets)/testsuperset.fasta\"","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:13.345049Z","iopub.execute_input":"2023-08-09T07:15:13.345622Z","iopub.status.idle":"2023-08-09T07:15:13.354702Z","shell.execute_reply.started":"2023-08-09T07:15:13.345587Z","shell.execute_reply":"2023-08-09T07:15:13.353768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# load the data set\ntrain_terms = pd.read_csv(train_labels_path, sep = \"\\t\")\ntrain_terms.shape","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:13.356245Z","iopub.execute_input":"2023-08-09T07:15:13.356669Z","iopub.status.idle":"2023-08-09T07:15:17.384157Z","shell.execute_reply.started":"2023-08-09T07:15:13.356609Z","shell.execute_reply":"2023-08-09T07:15:17.382841Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms.head()","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:17.385789Z","iopub.execute_input":"2023-08-09T07:15:17.386233Z","iopub.status.idle":"2023-08-09T07:15:17.399996Z","shell.execute_reply.started":"2023-08-09T07:15:17.386189Z","shell.execute_reply":"2023-08-09T07:15:17.398538Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# load precalculated protein embeddings\ntrain_protein_id = np.load('/kaggle/input/t5embeds/train_ids.npy')\ntrain_protein_id[:5]","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:17.401956Z","iopub.execute_input":"2023-08-09T07:15:17.402407Z","iopub.status.idle":"2023-08-09T07:15:17.461662Z","shell.execute_reply.started":"2023-08-09T07:15:17.402372Z","shell.execute_reply":"2023-08-09T07:15:17.460354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Each protein embedding is a vector of length 1024.\n# Protein embeddings are a way to encode the functional and structural properties of a protein\n# into a machine-freindly format.\n\ntrain_embeddings = np.load('/kaggle/input/t5embeds/train_embeds.npy')\n#train_embeddings[:5]","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:17.462956Z","iopub.execute_input":"2023-08-09T07:15:17.463292Z","iopub.status.idle":"2023-08-09T07:15:28.297699Z","shell.execute_reply.started":"2023-08-09T07:15:17.463262Z","shell.execute_reply":"2023-08-09T07:15:28.296748Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create a pandas dataframe with 1024 columns to represent each value in the vector's 1024 \n# places.\ncol_num = train_embeddings.shape[1]\ntrain_df = pd.DataFrame(train_embeddings, columns = [\"col\"+str(i) for i in range (1, col_num + 1)])\n\n# add indexing\ntrain_df.index = train_protein_id\ntrain_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:28.299125Z","iopub.execute_input":"2023-08-09T07:15:28.299437Z","iopub.status.idle":"2023-08-09T07:15:28.349919Z","shell.execute_reply.started":"2023-08-09T07:15:28.299409Z","shell.execute_reply":"2023-08-09T07:15:28.348601Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('train_protein_id shape: {0}\\ntrain_df shape:{1}'.format(train_protein_id.shape, train_df.shape))","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:28.351753Z","iopub.execute_input":"2023-08-09T07:15:28.352093Z","iopub.status.idle":"2023-08-09T07:15:28.358042Z","shell.execute_reply.started":"2023-08-09T07:15:28.352063Z","shell.execute_reply":"2023-08-09T07:15:28.356913Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparing the dataset","metadata":{}},{"cell_type":"code","source":"# First, extract all the needed labels (GO Term ID's) from the train_terms.tsv file. To \n# simplify the model, we choose the most frequent 1500 GO Term Ids as labels.\n\nnum_labels = 1000\n\nlabels = train_terms['term'].value_counts().index[:num_labels].tolist()\nlabels[:5]","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:28.359652Z","iopub.execute_input":"2023-08-09T07:15:28.360434Z","iopub.status.idle":"2023-08-09T07:15:29.408545Z","shell.execute_reply.started":"2023-08-09T07:15:28.360388Z","shell.execute_reply":"2023-08-09T07:15:29.40768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create a new datafram by filtering train_terms with the selected labels\n\n# fetch the train_terms data for relevant labels only!\ntrain_terms_updated = train_terms.loc[train_terms['term'].isin(labels)]\ntrain_terms_updated.head()","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:29.413207Z","iopub.execute_input":"2023-08-09T07:15:29.414088Z","iopub.status.idle":"2023-08-09T07:15:30.198832Z","shell.execute_reply.started":"2023-08-09T07:15:29.414049Z","shell.execute_reply":"2023-08-09T07:15:30.198007Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Since this is a mutli label classification problem, in the array we denote the presence\n# or absence of each Go Term Id using a 1 or 0.\n\n# Setting up progressbar (aesthetic)\nbar = progressbar.ProgressBar(maxval = num_labels, \\\n                             widgets = [progressbar.Bar('=', '[', ']'), ' ', progressbar.Percentage()])\n\n# Create empty dataframe of size (train_size x num_labels) for storing labels\ntrain_size = train_protein_id.shape[0]\ntrain_labels = np.zeros((train_size, num_labels))\n\n#train_labels.shape\n\n# Convert from numpy to pandas series for better handling\ntrain_protein_id_series = pd.Series(train_protein_id)\n\ntrain_protein_id_series","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:30.200116Z","iopub.execute_input":"2023-08-09T07:15:30.200795Z","iopub.status.idle":"2023-08-09T07:15:30.230123Z","shell.execute_reply.started":"2023-08-09T07:15:30.200761Z","shell.execute_reply":"2023-08-09T07:15:30.228808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Loop through each label to update the train_labels array with the appropriate values.\n\nfor i in range(num_labels):\n    \n    # for each label, fetch train_term data that has the term label[i]\n    n_train_terms = train_terms_updated[train_terms_updated['term'] == labels[i]]\n    \n    # fetch all unique related proteins (under EntryID) related to current label (GO 'term' ID)\n    label_related_proteins = n_train_terms['EntryID'].unique()\n    \n    # If a protein in train_protein_id_series is related to the current label, mark it as 1 \n    # (else 0). Replace the i'th column with the updated version of train_protein_id_series.\n    train_labels[:,i] = train_protein_id_series.isin(label_related_proteins).astype(float)\n    \n    # Increase progress bar %\n    bar.update(i+1)\n    \n\n# Notify the end of the progress bar\nbar.finish()\n","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:15:30.231297Z","iopub.execute_input":"2023-08-09T07:15:30.231647Z","iopub.status.idle":"2023-08-09T07:28:23.089886Z","shell.execute_reply.started":"2023-08-09T07:15:30.231597Z","shell.execute_reply":"2023-08-09T07:28:23.088647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert train_Y numby into pandas dataframe\nlabels_df = pd.DataFrame(data = train_labels, columns = labels)\nlabels_df.index = train_protein_id_series\nlabels_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:23.091486Z","iopub.execute_input":"2023-08-09T07:28:23.092385Z","iopub.status.idle":"2023-08-09T07:28:23.130889Z","shell.execute_reply.started":"2023-08-09T07:28:23.09235Z","shell.execute_reply":"2023-08-09T07:28:23.129712Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.index = train_protein_id_series","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:23.132413Z","iopub.execute_input":"2023-08-09T07:28:23.13289Z","iopub.status.idle":"2023-08-09T07:28:23.13999Z","shell.execute_reply.started":"2023-08-09T07:28:23.132845Z","shell.execute_reply":"2023-08-09T07:28:23.138777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:23.141319Z","iopub.execute_input":"2023-08-09T07:28:23.141679Z","iopub.status.idle":"2023-08-09T07:28:23.323569Z","shell.execute_reply.started":"2023-08-09T07:28:23.141648Z","shell.execute_reply":"2023-08-09T07:28:23.322237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Training\nK-fold cross validation from scratch.","metadata":{}},{"cell_type":"code","source":"# Merge the labels and training data\nmerged_df = pd.merge(train_df, labels_df, left_index=True, right_index=True, how='left')\n\n# randomize the data\n#merged_random_df = merged_random_df.reindex(np.random.permutation(merged_random_df.index))\n\n# get all test data\ntest_df = pd.DataFrame(np.load('/kaggle/input/t5embeds/test_embeds.npy'))\n\nmerged_df.head(4)\n\n# split `merged_df` into data and labels\ndata_df = merged_df.iloc[:, : train_df.shape[1]]\nlabels_df = merged_df.iloc[:, -labels_df.shape[1]:]\n\ndata_df","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:23.325371Z","iopub.execute_input":"2023-08-09T07:28:23.325735Z","iopub.status.idle":"2023-08-09T07:28:39.630585Z","shell.execute_reply.started":"2023-08-09T07:28:23.325704Z","shell.execute_reply":"2023-08-09T07:28:39.629354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Early stopping\n\nfrom tensorflow import keras\nfrom tensorflow.keras import layers, callbacks\nearly_stopping = callbacks.EarlyStopping(\n    min_delta = 0.0005,\n    patience = 5,\n    restore_best_weights = True,\n)\n","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:39.632538Z","iopub.execute_input":"2023-08-09T07:28:39.633016Z","iopub.status.idle":"2023-08-09T07:28:39.640618Z","shell.execute_reply.started":"2023-08-09T07:28:39.63297Z","shell.execute_reply":"2023-08-09T07:28:39.639303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def training(train_data, train_labels, val_data, val_labels, test_data):\n    \n    \"\"\"\n    Helper function that does the training \n    \"\"\"\n    \n    INPUT_SHAPE = [train_df.shape[1]]\n    BATCH_SIZE = 500\n\n    model = tf.keras.Sequential([\n        tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),\n        layers.Dropout(rate = 0.3),\n        tf.keras.layers.Dense(units = 512, activation = 'relu'),\n        tf.keras.layers.Dense(units = 512, activation = 'relu'),\n        tf.keras.layers.Dense(units = 512, activation = 'relu'),\n        tf.keras.layers.Dense(units = num_labels, activation = 'sigmoid')\n    ])\n\n    # Compile model\n    model.compile(\n        optimizer=tf.keras.optimizers.Adam(learning_rate=0.001),\n        loss='binary_crossentropy',\n        metrics=['binary_accuracy', tf.keras.metrics.AUC()],\n    )\n\n    # train\n    history = model.fit(\n        train_data, train_labels,\n        batch_size=BATCH_SIZE,\n        validation_data = (val_data, val_labels),\n        shuffle = True,\n        epochs= 20,\n        callbacks = [early_stopping],\n        verbose = 0,\n    )\n    \n    \n    history_df = pd.DataFrame(history.history)\n    history_df.loc[:, ['loss', 'val_loss']].plot(title = \"Loss/Validation Loss\")\n    history_df.loc[:, ['binary_accuracy']].plot(title=\"Accuracy\")\n    history_df.loc[:, ['loss']].plot(title=\"Accuracy\")\n\n    plt.show()\n\n    return model.predict(test_df)","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:39.641945Z","iopub.execute_input":"2023-08-09T07:28:39.642873Z","iopub.status.idle":"2023-08-09T07:28:39.658212Z","shell.execute_reply.started":"2023-08-09T07:28:39.642828Z","shell.execute_reply":"2023-08-09T07:28:39.657161Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def k_fold_cross_validation_prediction(data, k, labels, test_data):\n    \n    \"\"\"\n    Preform k-fold cross validation:\n    \n    Parameters:\n        k - number of rows\n        \n    \n    \"\"\"\n    num_rows = data.shape[0]\n    fold_size = num_rows // k\n    prediction = np.empty((test_data.shape[0], labels.shape[1]))\n\n    for fold in range(k):\n        \n        start = data.index[fold * fold_size]\n        end = data.index[(fold + 1) * fold_size]\n    \n        # Split data into train and test sets\n        val_data = data.loc[start:end]\n        val_labels = labels.loc[start:end]\n        train_data = np.concatenate([data.loc[:start], data.loc[end:]])\n        train_labels = np.concatenate([labels.loc[:start], labels.loc[end:]])\n    \n        \n        # Train the Model\n        prediction += training(train_data, train_labels, val_data, val_labels, test_data)\n        \n    # return the average predictions\n    return (prediction / k)","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:39.659707Z","iopub.execute_input":"2023-08-09T07:28:39.660509Z","iopub.status.idle":"2023-08-09T07:28:39.679188Z","shell.execute_reply.started":"2023-08-09T07:28:39.660466Z","shell.execute_reply":"2023-08-09T07:28:39.677767Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prediction = k_fold_cross_validation_prediction(data_df, 5, labels_df, test_df)\n\nprediction","metadata":{"execution":{"iopub.status.busy":"2023-08-09T07:28:39.681537Z","iopub.execute_input":"2023-08-09T07:28:39.681928Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"type(prediction)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prediction.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission\nfrom https://www.kaggle.com/code/gusthema/cafa-5-protein-function-with-tensorflow/notebook !","metadata":{}},{"cell_type":"code","source":"df_submission = pd.DataFrame(columns = ['Protein Id', 'GO Term Id','Prediction'])\ntest_protein_ids = np.load('/kaggle/input/t5embeds/test_ids.npy')\nl = []\nfor k in list(test_protein_ids):\n    l += [ k] * prediction.shape[1]   \n\ndf_submission['Protein Id'] = l\ndf_submission['GO Term Id'] = labels * prediction.shape[0]\ndf_submission['Prediction'] = prediction.ravel()\ndf_submission.to_csv(\"submission.tsv\",header=False, index=False, sep=\"\\t\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_submission","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Combining Scores\n#this submission was obtained by training models on BlstP and Foldseek offline\nsubmission2 = pd.read_csv('/kaggle/input/cafa5-using-protbert-embeds/submission.tsv',\n    sep='\\t', header=None, names=['Protein Id', 'GO Term Id', 'Prediction']) ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submissions = submission2.merge(df_submission, left_on=['Protein Id', 'GO Term Id'], \n                                                  right_on=['Protein Id', 'GO Term Id'], how='outer')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submissions.drop(['Id', 'GO term'], axis=1, inplace=True)\nsubmissions['confidence_combined'] = submissions.apply(lambda row: row['Confidence2'] if not np.isnan(row['Confidence2']) else row['Confidence'], axis=1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submissions[['Id2', 'GO term2', 'confidence_combined']].to_csv('submission.tsv', sep='\\t', header=False, index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submissions","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}