{"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":"*This notebook is based on [Gusthema's starter notebook](https://www.kaggle.com/code/gusthema/cafa-5-protein-function-with-tensorflow).*\n\n# Competition Objective\nThe objective of the competition is to predict [Gene Ontology (GO)](http://geneontology.org/) terms associated with a set of test protein IDs to provide indicators for protein function. The GO terms need to be predicted at the level of Biological Process (BP), Molecular Function (MF) and Cellular Component (CC) in the GO tree.","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-08-21T08:42:39.177005Z","iopub.execute_input":"2023-08-21T08:42:39.177389Z","iopub.status.idle":"2023-08-21T08:42:39.184407Z","shell.execute_reply.started":"2023-08-21T08:42:39.177359Z","shell.execute_reply":"2023-08-21T08:42:39.183377Z"}}},{"cell_type":"markdown","source":"# Available input data used for analysis\n\nThe train and test datasets provided include protein IDs (Uniprot accession ID) and their respective protein sequences in fasta format. Corresponding taxon IDs for each protein ID have been provided in separate files. This gives an idea of the different species associated with the protein IDs.  \nThe trainset also includes the ground truth, experimentally verified or inferred GO term IDs that each protein has been annotated with. This is the prediction target.\n\nSince the input feature includes protein sequences, Natural Language Processing (NLP) methods can potentially provide insights into structural and functional information contained in the sequences. This is because each protein sequence is represented as a chain of different combinations of 22 known amino acids. Each amino acid is an organic compound represented by one letter of the alphabet, e.g. 'R' represents the amino acid arginine and 'S' represents serine.","metadata":{}},{"cell_type":"markdown","source":"## Protein embeddings\n\nDifferent embeddings have been generated for the train and test sets in the competition using models pretrained on large protein sequence corpuses. These include [EMS2](https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/406168), [ProtBert](https://www.kaggle.com/datasets/henriupton/protbert-embeddings-for-cafa5) and [ProtT5](https://www.kaggle.com/code/siddhvr/cafa-5-t5-embeds). \n\nThis paper cites that ProtT5 was more informative than ProtBert for downstream tasks: *Elnaggar A, Heinzinger M, Dallago C, Rehawi G, Wang Y, Jones L, Gibbs T, Feher T, Angerer C, Steinegger M, Bhowmik D, Rost B. ProtTrans: Toward Understanding the Language of Life Through Self-Supervised Learning. IEEE Trans Pattern Anal Mach Intell. 2022 Oct;44(10)*","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport os","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load the T5 embeds file.\ntrain_embeds_file = '/kaggle/input/t5embeds/train_embeds.npy'\ntrain_embeds_np = np.load(train_embeds_file)\ntrain_embeds_df = pd.DataFrame(train_embeds_np, columns=['Dim_'+ str(i) for i in np.arange(1, train_embeds_np.shape[1]+1)])\nprint(f\"The shape of the trainset embeddings dataframe: {train_embeds_df.shape}\")\n\n#load corresponding train IDs file\ntrain_ids_file = \"/kaggle/input/t5embeds/train_ids.npy\"\ntrain_ids = np.load(train_ids_file)\ntrain_ids_df = pd.DataFrame(train_ids, columns=[\"EntryID\"])\nprint(f\"The shape of the trainset ids dataframe: {train_ids_df.shape}\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_ids_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_embeds_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_embeds_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Encode the target variable \nThis is a multilabel classification problem since one protein can be associated with multiple GO term IDs. So, the GO term ID classes need to be encoded appropriately.\n\nSome GO term IDs appear very few times in the trainset. Training on these would not provide sufficient information to the model. So, only GO term IDs that appear 300 or more times are considered here.","metadata":{}},{"cell_type":"code","source":"#load the train terms file with ground truth target GO term IDs for supervised learning.\ntrain_terms_file = \"/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv\"\ntrain_terms = pd.read_csv(train_terms_file, sep=\"\\t\")\nprint(f\"The shape of the trainset GO term ID mappings dataframe: {train_terms.shape}\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.preprocessing import MultiLabelBinarizer","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"The number of unique GO term IDs associated with trainset proteins: {train_terms['term'].nunique()}\")\ntrain_term_cnts = train_terms['term'].value_counts()\n\ntrain_term_grt_300 = train_term_cnts[train_term_cnts >= 300]\nprint(f\"The number of GO term IDs that occur 300 or more times in the target: {len(train_term_grt_300)}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"target_all = train_terms.groupby('EntryID')['term'].apply(tuple)\nprint(f\"The shape of the target: {target_all.shape}\")\n\ntarget_list = list()\n\nfor val in target_all.values:\n    tmp = tuple([go_val for go_val in val if go_val in train_term_grt_300])\n    target_list.append(tmp)\n\ntarget = pd.Series(target_list, index=target_all.index)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"encoder = MultiLabelBinarizer()\nmulti_label_enc = encoder.fit_transform(target)\nmulti_label_enc_df = pd.DataFrame(multi_label_enc, columns=encoder.classes_)\nmulti_label_enc_df.shape","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Define and train model\nA DNN running on a CPU with dropout regularization and batch normalization layers is used here.","metadata":{"execution":{"iopub.status.busy":"2023-08-21T15:38:12.259386Z","iopub.execute_input":"2023-08-21T15:38:12.259833Z","iopub.status.idle":"2023-08-21T15:38:12.265508Z","shell.execute_reply.started":"2023-08-21T15:38:12.259773Z","shell.execute_reply":"2023-08-21T15:38:12.263903Z"}}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nimport tensorflow as tf\nfrom tensorflow.keras import layers, callbacks\nfrom tensorflow.keras.models import Sequential\n\ntrain_embeds_data, validation_embeds_data, train_go, validation_go = train_test_split(\n    train_embeds_df, multi_label_enc_df, test_size=0.3, random_state=0\n)\n\nprint(f\"The shape of the training data: {train_embeds_data.shape}\")\nprint(f\"The shape of the validation data: {validation_embeds_data.shape}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"input_shape = [train_embeds_data.shape[1]]\noutput_units = train_go.shape[1]\ndropout_rate = 0.3\nbatch_size = 2048\nepochs = 100\n\nmodel = Sequential(\n    [\n        layers.BatchNormalization(input_shape=input_shape),\n        layers.Dense(units=512, activation='relu'),\n        layers.Dropout(rate=dropout_rate),\n        layers.BatchNormalization(),\n        layers.Dense(units=384, activation='relu'),\n        layers.Dropout(rate=dropout_rate),\n        layers.BatchNormalization(),\n        layers.Dense(units=output_units,activation='sigmoid')\n        \n    ]\n)\n\nlearning_rate = 0.001\n\n\nmodel.compile(\n    optimizer = tf.keras.optimizers.Adam(learning_rate = learning_rate),\n    loss = 'binary_crossentropy',\n    metrics = ['binary_accuracy', tf.keras.metrics.AUC()]\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"early_stopping = callbacks.EarlyStopping(\n    min_delta = 0.001,\n    patience = 10,\n    restore_best_weights = True\n)\n\nhistory = model.fit(\n    train_embeds_data, train_go,\n    validation_data = (validation_embeds_data, validation_go),\n    batch_size = batch_size,\n    epochs = epochs,\n    callbacks = [early_stopping]\n)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"history_df = pd.DataFrame(history.history)\nhistory_df.loc[:, ['loss', 'val_loss']].plot(title='Cross entropy')\nhistory_df.loc[:, ['binary_accuracy', 'val_binary_accuracy']].plot(title='Accuracy')\nhistory_df.loc[:, list(history_df.filter(regex='auc').columns)].plot(title='AUC')\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Predictions on test dataset","metadata":{}},{"cell_type":"code","source":"test_embeds_np = np.load(\"/kaggle/input/t5embeds/test_embeds.npy\")\ntest_embeds_df = pd.DataFrame(test_embeds_np, columns = ['Dim' + str(i) for i in np.arange(1, test_embeds_np.shape[1]+1)])\nprint(test_embeds_df.shape)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_ids_file = \"/kaggle/input/t5embeds/test_ids.npy\"\ntest_ids = np.load(test_ids_file)\ntest_ids_df = pd.DataFrame(test_ids, columns=[\"EntryID\"])\nprint(f\"The shape of test ids: {test_ids.shape}\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"predictions = model.predict(test_embeds_df)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Create a submissions file\nThe competition has a limit of a maximum of 1500 predictions for each protein ID, only predictions that have a probability >= 0.1 or 10% are included in the output file.\n","metadata":{}},{"cell_type":"code","source":"import csv\nimport progressbar\n\npred_line_counter = 0\n\noutfile = \"/kaggle/working/submission.tsv\"\nif os.path.exists(outfile):\n    os.remove(outfile)\n    \nbar = progressbar.ProgressBar(maxval=predictions.shape[0], \\\n    widgets=[progressbar.Bar('=', 'Writing to submission.tsv[', ']'), ' ', progressbar.Percentage()])\n    \nwith open(outfile, 'a', newline='') as csv_file:\n    write = csv.writer(csv_file, delimiter=\"\\t\")\n    for pred in predictions:\n        proteinID = test_ids[pred_line_counter]\n        sorted_preds = sorted(pred, reverse=True)\n        \n        #only write GO values if predicted probability >= 0.1 (10%).\n        relevant_pred = [[test_ids[pred_line_counter], train_go.columns[i], f'{val:.03f}'] for i in np.arange(0, len(sorted_preds)) if (val:= sorted_preds[i]) >= 0.1 and val !=0]\n        write.writerows(relevant_pred)\n        \n        pred_line_counter+=1\n        bar.update(pred_line_counter)\n        \n\nbar.finish()\ncsv_file.close()    \n    \n    ","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_embeds_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"res_look_up = pd.read_csv(\"/kaggle/working/submission.tsv\", sep=\"\\t\", header=None)\nres_look_up.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Notes\nThe model outputs similar probabilities for GO terms for almost all proteins. It is very unlikely it is a good predictor for practical purposes. It needs to be improved. Besides, it could have been overfit on the training data.","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}