{"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":"### CAFA Submission -  Segregating the training and evaluation processes for each subontology\n1. Divide the training dataset into three subsets, corresponding to each Gene Ontology (GO) subontology.\n2. Further segregate each subontology training dataset into 2 components: Training, Validation.\n3. Identify the top GO Terms from each subontology for use in classification training.\n4. Develop various training models and evaluate their effectiveness for each GO term prediction using an AUC-ROC and F1/threshold combination.\n5. Submit the Go-Term predictions from it's top-performing model.","metadata":{}},{"cell_type":"code","source":"#supress warning\nimport warnings\nwarnings.filterwarnings('ignore')\n\n!pip install obonet -q\nimport obonet\nimport pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\n\nimport tensorflow as tf\nfrom sklearn.metrics import f1_score, precision_score, recall_score, confusion_matrix, classification_report\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import roc_auc_score, roc_curve, auc","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:40:37.518717Z","iopub.execute_input":"2023-06-20T13:40:37.519522Z","iopub.status.idle":"2023-06-20T13:41:05.252448Z","shell.execute_reply.started":"2023-06-20T13:40:37.519484Z","shell.execute_reply":"2023-06-20T13:41:05.251016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#read train terms (EntryID, term, aspect)\ntrain_terms = pd.read_csv(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv\",sep=\"\\t\")\n\n#get list of EntryIDs in each aspect\ntrain_entries = train_terms[\"EntryID\"].unique()\nbpo_entries = train_terms[train_terms[\"aspect\"]==\"BPO\"][\"EntryID\"].unique()\nmfo_entries = train_terms[train_terms[\"aspect\"]==\"MFO\"][\"EntryID\"].unique()\ncco_entries = train_terms[train_terms[\"aspect\"]==\"CCO\"][\"EntryID\"].unique()\n\n#get list of EntryIDs without BPO terms, MFO terms, CCO terms\nno_bpo_entries = set(train_entries) - set(bpo_entries)\nno_mfo_entries = set(train_entries) - set(mfo_entries)\nno_cco_entries = set(train_entries) - set(cco_entries)\n\n#show number of bpo_entries, mfo_entries, cco_entries, train_entries in bar chart\nvalues = [len(bpo_entries), len(mfo_entries), len(cco_entries), len(train_entries), len(no_bpo_entries), len(no_mfo_entries), len(no_cco_entries)]\nweights = ['BPO Entries', 'MFO Entries', 'CCO Entries', 'Train Entries', 'wihtout-BPO', 'without-MFO', 'wihtout-CCO']\n\nplt.figure(figsize=(12,6))\nplt.bar(weights, values, color=['blue', 'blue', 'blue', 'green', 'red', 'red', 'red'])\nplt.xlabel('Subontology')\nplt.ylabel('Number of Unique Entries')\nplt.title('Number of Entries in Each Subontology')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:41:08.991664Z","iopub.execute_input":"2023-06-20T13:41:08.992563Z","iopub.status.idle":"2023-06-20T13:41:18.646629Z","shell.execute_reply.started":"2023-06-20T13:41:08.992525Z","shell.execute_reply":"2023-06-20T13:41:18.64484Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*It appears that the number of protein IDs associated with each subontology annotation is roughly equivalent. Let's further analyze the total annotations within each subontology.*","metadata":{}},{"cell_type":"code","source":"#get total entries for each subontology\nbpo_train_terms = train_terms[train_terms[\"aspect\"]==\"BPO\"]\ncco_train_terms = train_terms[train_terms[\"aspect\"]==\"CCO\"]\nmfo_train_terms = train_terms[train_terms[\"aspect\"]==\"MFO\"]\n\nvalues = [len(bpo_train_terms), len(mfo_train_terms), len(cco_train_terms), len(train_terms)]\nweights = ['BPO terms', 'MFO terms', 'CCO terms', 'Train Entries']\n\nplt.figure(figsize=(8,6))\nplt.bar(weights, values, color=['blue', 'blue', 'blue', 'green'])\nplt.xlabel('Subontology')\nplt.ylabel('Number of Training Entries')\nplt.title('Number of Entries in Each Subontology')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:41:24.302241Z","iopub.execute_input":"2023-06-20T13:41:24.302691Z","iopub.status.idle":"2023-06-20T13:41:27.710211Z","shell.execute_reply.started":"2023-06-20T13:41:24.302659Z","shell.execute_reply":"2023-06-20T13:41:27.708845Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"While the distribution of protein IDs annotations across subontologies appears to be fairly even, the Biological Process Ontology (BPO) seems to have a higher number of annotations. Let's examine these subontology numbers further in the Gene Ontology (GO) graph.","metadata":{}},{"cell_type":"code","source":"# read obo file provided by CAFA in training data\nobo_file = '/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo'\nobo_graph = obonet.read_obo(obo_file)\n#check number of nodes in each subontology\nname_space = ['biological_process', 'molecular_function', 'cellular_component']\nnumber_of_nodes = []\nunknown_count = 0\nfor names in name_space:\n    name_space_count = 0\n    for node, data in obo_graph.nodes(data=True):\n        if data['namespace'] == names:\n            name_space_count += 1\n        #if data['namespace'] not in name_space, then increase others count\n        elif data['namespace'] not in name_space:\n            unknown_count += 1\n    number_of_nodes.append(name_space_count)\n    print(\"The graph has {} nodes with a namespace of '{}'.\".format(name_space_count, names))\n\nprint(\"The graph has {} nodes with a unknown namespaces.\".format(unknown_count))\n#plt.figure(figsize=(12,6))\nplt.bar(name_space, number_of_nodes, color=['blue', 'green', 'red'])\nplt.xlabel('')\nplt.ylabel('Number of Nodes')\nplt.title('Number of Nodes in Each Subontology')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:41:31.542108Z","iopub.execute_input":"2023-06-20T13:41:31.542527Z","iopub.status.idle":"2023-06-20T13:41:51.869096Z","shell.execute_reply.started":"2023-06-20T13:41:31.542496Z","shell.execute_reply":"2023-06-20T13:41:51.867954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*This observation is somewhat consistent with the previous plot. It seems the Gene Ontology graph inherently possesses a larger number of BPO nodes.*","metadata":{}},{"cell_type":"code","source":"#To understand the distribution of GO Terms by each aspect, let's plot it:\ndef plot_go_terms_dist(go_terms, aspect, plot_terms = 100):\n    plot_df = go_terms['term'].value_counts().iloc[:plot_terms]\n    figure, axis = plt.subplots(1, 1, figsize=(24, 3))\n    bp = sns.barplot(ax=axis, x=np.array(plot_df.index), y=plot_df.values)\n    bp.set_xticklabels(bp.get_xticklabels(), rotation=90, size = 6)\n    title = aspect + 'Top ' + str(plot_terms) +  ' frequent GO term IDs'\n    axis.set_title(title)\n    bp.set_xlabel(\"GO term IDs\", fontsize = 12)\n    bp.set_ylabel(\"Count\", fontsize = 12)\n    plt.show()\n\n\nplot_go_terms_dist(bpo_train_terms, 'BPO ', 500)\nplot_go_terms_dist(mfo_train_terms, 'MFO ')\nplot_go_terms_dist(cco_train_terms, 'CCO ')","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:42:55.082539Z","iopub.execute_input":"2023-06-20T13:42:55.082939Z","iopub.status.idle":"2023-06-20T13:43:04.240101Z","shell.execute_reply.started":"2023-06-20T13:42:55.082909Z","shell.execute_reply":"2023-06-20T13:43:04.238593Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*If we visually compare the top (first) GO Term from each aspect in the current plot with the previous one (representing the number of entries by subontology), it seems they align approximately, it's likely that we will not have a '0' target classification for the top GO term in each subontology dataset.*","metadata":{}},{"cell_type":"code","source":"# Help from https://www.kaggle.com/code/gusthema/cafa-5-protein-function-with-tensorflow\n# Set the limit for labels to be considered for each aspect\n# Note: num_of_labels is for each aspect, not for the whole dataset. Total number of labels will be 3x num_of_labels\nnum_of_labels = 300\n\n#generate binary labels for each aspect\ndef generate_binary_labels(subont_df, num_of_labels):\n    # Take value counts in descending order and fetch top GO terms\n    labels = subont_df['term'].value_counts().index[:num_of_labels].tolist()\n    # Fetch the data for the relevant labels only\n    subont_df_top = subont_df.loc[train_terms['term'].isin(labels)]\n\n    # Create an empty dataframe of required size for storing the labels,\n    # i.e, sub ontology entries x num_of_labels\n    train_labels = np.zeros((len(subont_df['EntryID'].unique()) ,num_of_labels))\n    # Convert from numpy to pandas series for better handling\n    series_train_protein_ids = pd.Series(subont_df['EntryID'].unique())\n\n    # Loop through each label\n    for i in range(num_of_labels):\n        n_train_terms = subont_df_top[subont_df_top['term'] ==  labels[i]]\n        # Fetch all the unique EntryId aka proteins related to the current label(GO term ID)\n        label_related_proteins = n_train_terms['EntryID'].unique()\n        \n        # In the series_train_protein_ids pandas series, if a protein is related\n        # to the current label, then mark it as 1, else 0.\n        # Replace the ith column of train_Y with with that pandas series.\n        train_labels[:,i] =  series_train_protein_ids.isin(label_related_proteins).astype(float)\n        \n        # Convert train_labels numpy into pandas dataframe\n    labels_df = pd.DataFrame(data = train_labels, columns = labels)\n    # Add the protein ids as a column\n    labels_df['EntryID'] = series_train_protein_ids\n    return labels_df\n\ncco_binary_terms = generate_binary_labels(cco_train_terms, num_of_labels)\nprint(cco_binary_terms.shape)\nprint(cco_binary_terms.head(1))\nbpo_binary_terms = generate_binary_labels(bpo_train_terms, num_of_labels)\nprint(bpo_binary_terms.shape)\nprint(bpo_binary_terms.head(1))\nmfo_binary_terms = generate_binary_labels(mfo_train_terms, num_of_labels)\nprint(mfo_binary_terms.shape)\nprint(mfo_binary_terms.head(1))\n#check if binary dataframes have same number of entries\nprint('\\n Do the three binary dataframes contain all the required entries?')\nprint(len(bpo_binary_terms) == len(bpo_entries), len(cco_binary_terms) == len(cco_entries), len(mfo_binary_terms) == len(mfo_entries))","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:43:08.964016Z","iopub.execute_input":"2023-06-20T13:43:08.964472Z","iopub.status.idle":"2023-06-20T13:46:49.733043Z","shell.execute_reply.started":"2023-06-20T13:43:08.964437Z","shell.execute_reply":"2023-06-20T13:46:49.731766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#load priotein embeddings \n# embeddings from https://www.kaggle.com/datasets/sergeifironov/t5embeds\ntrain_protein_ids = np.load('/kaggle/input/t5embeds/train_ids.npy')\nprint('pritein ids: ', train_protein_ids.shape)\ntrain_embeddings = np.load('/kaggle/input/t5embeds/train_embeds.npy')\nprint('train_embeddings type: ', type(train_embeddings), 'train_embeddings shape: ', train_embeddings.shape)\n\n# Now lets convert embeddings numpy array(train_embeddings) into pandas dataframe.\ncolumn_num = train_embeddings.shape[1]\ntrain_df = pd.DataFrame(train_embeddings, columns = [\"Column_\" + str(i) for i in range(1, column_num+1)])\n# Add the protein id as a column to help prepare embeddings for each subontology later\ntrain_df['EntryID'] = train_protein_ids\nprint(train_df.shape)\ntrain_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:47:17.502245Z","iopub.execute_input":"2023-06-20T13:47:17.502693Z","iopub.status.idle":"2023-06-20T13:47:33.727087Z","shell.execute_reply.started":"2023-06-20T13:47:17.502656Z","shell.execute_reply":"2023-06-20T13:47:33.725949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# append embeddings to GO Term binary classifications by Entry ID to create seperate embedding sets for each subontology training\nbpo_seq_binary_terms = bpo_binary_terms.merge(train_df, on='EntryID', how='left')\ncco_seq_binary_terms = cco_binary_terms.merge(train_df, on='EntryID', how='left')\nmfo_seq_binary_terms = mfo_binary_terms.merge(train_df, on='EntryID', how='left')","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:47:44.358287Z","iopub.execute_input":"2023-06-20T13:47:44.359074Z","iopub.status.idle":"2023-06-20T13:47:55.260323Z","shell.execute_reply.started":"2023-06-20T13:47:44.359024Z","shell.execute_reply":"2023-06-20T13:47:55.259156Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#split data for train, valuation\ndef split_data(data):\n    train, val = train_test_split(data, test_size=0.3, random_state=42)\n    #get embedding columns as X and binary GO Term columns as Y\n    train_X = train.iloc[:,-1024:].reset_index(drop=True)\n    train_Y = train.iloc[:, :num_of_labels].reset_index(drop=True)\n    val_X = val.iloc[:,-1024:].reset_index(drop=True)\n    val_Y = val.iloc[:, :num_of_labels].reset_index(drop=True)\n\n    return train, val, train_X, train_Y, val_X, val_Y\n\nbpo_train, bpo_val, bpo_train_X, bpo_train_Y, bpo_val_X, bpo_val_Y = split_data(bpo_seq_binary_terms)\nprint(bpo_train.shape, bpo_val.shape, bpo_train_X.shape, bpo_train_Y.shape, bpo_val_X.shape, bpo_val_Y.shape)\n\ncco_train, cco_val, cco_train_X, cco_train_Y, cco_val_X, cco_val_Y = split_data(cco_seq_binary_terms)\nprint(cco_train.shape, cco_val.shape, cco_train_X.shape, cco_train_Y.shape, cco_val_X.shape, cco_val_Y.shape)\n\nmfo_train, mfo_val, mfo_train_X, mfo_train_Y, mfo_val_X, mfo_val_Y = split_data(mfo_seq_binary_terms)\nprint(mfo_train.shape, mfo_val.shape, mfo_train_X.shape, mfo_train_Y.shape, mfo_val_X.shape, mfo_val_Y.shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:47:58.292089Z","iopub.execute_input":"2023-06-20T13:47:58.29252Z","iopub.status.idle":"2023-06-20T13:48:05.052183Z","shell.execute_reply.started":"2023-06-20T13:47:58.292487Z","shell.execute_reply":"2023-06-20T13:48:05.05087Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#check protein sequences lengths\nfrom Bio import SeqIO\nfasta_file = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\nfasta = SeqIO.to_dict(SeqIO.parse(fasta_file, \"fasta\"))\n#convert fasta seq to dataframe \nfasta_list = [{atr: getattr(record, atr) for atr in ['seq']} for record in fasta.values()]\nfasta_df = pd.DataFrame(fasta_list)\n#find minium and maximum length of sequences\nfasta_df['seq_length'] = fasta_df['seq'].str.len()\nprint('Minimum length of Protein: ', fasta_df['seq_length'].min(), '\\nMaximum length of Protein: ', fasta_df['seq_length'].max())\n#get count of proteins by range of seq length\nfasta_df.groupby(pd.cut(fasta_df['seq_length'], np.arange(0, 40000, 5000))).count()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:49:16.160373Z","iopub.execute_input":"2023-06-20T13:49:16.160956Z","iopub.status.idle":"2023-06-20T13:49:21.471226Z","shell.execute_reply.started":"2023-06-20T13:49:16.16091Z","shell.execute_reply":"2023-06-20T13:49:21.469365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#release memory\nimport gc\n\ndel fasta, fasta_df, train_df, train_embeddings, obo_graph\ndel train_entries, bpo_entries, cco_entries, mfo_entries, no_bpo_entries, no_mfo_entries, no_cco_entries\ndel bpo_seq_binary_terms, cco_seq_binary_terms, mfo_seq_binary_terms\ndel cco_binary_terms, bpo_binary_terms, mfo_binary_terms, bpo_train, cco_train, mfo_train\n\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:49:40.302727Z","iopub.execute_input":"2023-06-20T13:49:40.303186Z","iopub.status.idle":"2023-06-20T13:49:40.891353Z","shell.execute_reply.started":"2023-06-20T13:49:40.303149Z","shell.execute_reply":"2023-06-20T13:49:40.890222Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color='red'>*Please note, the upcoming section will execute six training models for each subontology. This process may take approximately 45 minutes to complete.*\n    <font color='green'>*Should Jupyter session encounter any issues requiring a rerun, we can bypass the next two cells and run the cell that loads and runs the model.*","metadata":{}},{"cell_type":"code","source":"#train multiple models and choose the best one for each Go Term\n#Most of the protein sequence lenghts range from 3 to 5k, we will use models from large to small size neural networks\ndef train_predict_model(train_X, train_Y, val_X, subontology):\n\n    INPUT_SHAPE = [train_X.shape[1]]\n    BATCH_SIZE = train_X.shape[0] // 5\n    epochs = 10\n    \n    # Model 1\n    model1 = tf.keras.Sequential([\n        tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),    \n        tf.keras.layers.Dense(units=2048, activation='relu'),\n        tf.keras.layers.Dense(units=2048, activation='relu'),\n        tf.keras.layers.Dense(units=1024, activation='relu'),\n        tf.keras.layers.Dense(units=num_of_labels,activation='sigmoid')\n    ])\n\n    model1.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    history1 = model1.fit(train_X, train_Y, batch_size=BATCH_SIZE, epochs=epochs, verbose = 0)\n        \n    ## Model2\n    \n    \n    # Model 3\n    model3 = tf.keras.Sequential([\n        tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),    \n        tf.keras.layers.Dense(units=1024, activation='relu'),\n        tf.keras.layers.Dense(units=512, activation='relu'),\n        tf.keras.layers.Dense(units=256, activation='relu'),\n        tf.keras.layers.Dense(units=num_of_labels,activation='sigmoid')\n    ])\n\n    model3.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    history3 = model3.fit(train_X, train_Y, batch_size=BATCH_SIZE, epochs=epochs, verbose = 0)\n    \n     # Model 4\n    model4 = tf.keras.Sequential([\n        tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),    \n        tf.keras.layers.Dense(units=512, activation='relu'),\n        tf.keras.layers.Dense(units=256, activation='relu'),\n        tf.keras.layers.Dense(units=128, activation='relu'),\n        tf.keras.layers.Dense(units=num_of_labels,activation='sigmoid')\n    ])\n\n    model4.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    history4 = model4.fit(train_X, train_Y, batch_size=BATCH_SIZE, epochs=epochs, verbose = 0)\n   \n    # Model 5\n    model5 = tf.keras.Sequential([\n        tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),    \n        tf.keras.layers.Dense(units=256, activation='relu'),\n        tf.keras.layers.Dense(units=128, activation='relu'),\n        tf.keras.layers.Dense(units=num_of_labels,activation='sigmoid')\n    ])\n\n    model5.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    history5 = model5.fit(train_X, train_Y, batch_size=BATCH_SIZE, epochs=epochs, verbose = 0)\n    \n    # Model 6\n    model6 = tf.keras.Sequential([\n        tf.keras.layers.BatchNormalization(input_shape=INPUT_SHAPE),    \n        tf.keras.layers.Dense(units=64, activation='relu'),\n        tf.keras.layers.Dense(units=32, activation='relu'),\n        tf.keras.layers.Dense(units=num_of_labels,activation='sigmoid')\n    ])\n\n    model6.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    history6 = model6.fit(train_X, train_Y, batch_size=BATCH_SIZE, epochs=epochs, verbose = 0)\n    \n\n    #predict for val\n    val_pred1 = model1.predict(val_X)\n    #val_pred2 = model2.predict(val_X)\n    val_pred3 = model3.predict(val_X)\n    val_pred4 = model4.predict(val_X)\n    val_pred5 = model5.predict(val_X)\n    val_pred6 = model6.predict(val_X)\n    \n    return (model1, model3, model4, model5, model6, \n            history1, history3, history4, history5, history6,\n            val_pred1, val_pred3, val_pred4, val_pred5, val_pred6)\n\n#train and predict for each subontology\n(bpo_model1, bpo_model3, bpo_model4, bpo_model5, bpo_model6,\n    bpo_history1, bpo_history3, bpo_history4, bpo_history5, bpo_history6,\n    bpo_val_pred1, bpo_val_pred3, bpo_val_pred4, bpo_val_pred5, bpo_val_pred6) = train_predict_model(bpo_train_X, bpo_train_Y, bpo_val_X, 'BPO')\n(cco_model1, cco_model3, cco_model4, cco_model5, cco_model6,\n    cco_history1, cco_history3, cco_history4, cco_history5, cco_history6,\n    cco_val_pred1, cco_val_pred3, cco_val_pred4, cco_val_pred5, cco_val_pred6) = train_predict_model(cco_train_X, cco_train_Y, cco_val_X, 'CCO')\n(mfo_model1, mfo_model3, mfo_model4, mfo_model5, mfo_model6,\n    mfo_history1, mfo_history3, mfo_history4, mfo_history5, mfo_history6,\n    mfo_val_pred1, mfo_val_pred3, mfo_val_pred4, mfo_val_pred5, mfo_val_pred6) = train_predict_model(mfo_train_X, mfo_train_Y, mfo_val_X, 'MFO')\n","metadata":{"execution":{"iopub.status.busy":"2023-06-20T13:51:30.488677Z","iopub.execute_input":"2023-06-20T13:51:30.48917Z","iopub.status.idle":"2023-06-20T14:20:28.049887Z","shell.execute_reply.started":"2023-06-20T13:51:30.489124Z","shell.execute_reply":"2023-06-20T14:20:28.048197Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plot loss and accuracy \ndef plot_loss_accuracy(histories, aspect): \n    plt.figure(figsize=(20,10))\n    plt.subplot(1, 2, 1)\n    for model_num, history in enumerate(histories, start=1):\n        plt.plot(history.history['loss'], label=('Model ' + str(model_num)))\n        \n    plt.title('Loss for ' + aspect)\n    plt.xlabel('Epoch')\n    plt.ylabel('Loss')\n    plt.legend()\n\n    plt.subplot(1, 2, 2)\n    for model_num, history in enumerate(histories, start=1):\n        plt.plot(history.history['binary_accuracy'], label=('Model ' + str(model_num)))\n        \n    plt.title('Accuracy for ' + aspect)\n    plt.xlabel('Epoch')\n    plt.ylabel('Accuracy')\n    plt.show()\n\nplot_loss_accuracy([bpo_history1, bpo_history3, bpo_history4, bpo_history5, bpo_history6], 'BPO')\nplot_loss_accuracy([cco_history1, cco_history3, cco_history4, cco_history5, cco_history6], 'CCO')\nplot_loss_accuracy([mfo_history1, mfo_history3, mfo_history4, mfo_history5, mfo_history6], 'MFO')    ","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:32:01.749841Z","iopub.execute_input":"2023-06-20T14:32:01.75061Z","iopub.status.idle":"2023-06-20T14:32:03.867229Z","shell.execute_reply.started":"2023-06-20T14:32:01.750558Z","shell.execute_reply":"2023-06-20T14:32:03.86591Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*Indeed, it's evident that the smaller neural network model yields a higher loss and correspondingly lower accuracy. To gain further insight, let's compute the Area Under the Curve (AUC) scores for the evaluation dataset and analyze how these models perform for each individual GO Term.*","metadata":{}},{"cell_type":"code","source":"#Save all 18 models to disk\nfrom tensorflow.keras.models import load_model\nbpo_model1.save('CAFAmodels/bpo_model1.h5')\n#bpo_model2.save('CAFAmodels/bpo_model2.h5')\nbpo_model3.save('CAFAmodels/bpo_model3.h5')\nbpo_model4.save('CAFAmodels/bpo_model4.h5')\nbpo_model5.save('CAFAmodels/bpo_model5.h5')\nbpo_model6.save('CAFAmodels/bpo_model6.h5')\n\ncco_model1.save('CAFAmodels/cco_model1.h5')\n#cco_model2.save('CAFAmodels/cco_model2.h5')\ncco_model3.save('CAFAmodels/cco_model3.h5')\ncco_model4.save('CAFAmodels/cco_model4.h5')\ncco_model5.save('CAFAmodels/cco_model5.h5')\ncco_model6.save('CAFAmodels/cco_model6.h5')\n\nmfo_model1.save('CAFAmodels/mfo_model1.h5')\n#mfo_model2.save('CAFAmodels/mfo_model2.h5')\nmfo_model3.save('CAFAmodels/mfo_model3.h5')\nmfo_model4.save('CAFAmodels/mfo_model4.h5')\nmfo_model5.save('CAFAmodels/mfo_model5.h5')\nmfo_model6.save('CAFAmodels/mfo_model6.h5')","metadata":{"execution":{"iopub.status.busy":"2023-06-19T17:33:14.522834Z","iopub.execute_input":"2023-06-19T17:33:14.523248Z","iopub.status.idle":"2023-06-19T17:33:16.892511Z","shell.execute_reply.started":"2023-06-19T17:33:14.523209Z","shell.execute_reply":"2023-06-19T17:33:16.891529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# If need to skip above 2 cells, run this to load the models and proceed\n#cafamodels saved during version 3\n'''\nfrom tensorflow.keras.models import load_model\nbpo_model1 = load_model('/kaggle/input/CAFAmodels/bpo_model1.h5')\nbpo_model2 = load_model('/kaggle/input/CAFAmodels/bpo_model2.h5')\nbpo_model3 = load_model('/kaggle/input/CAFAmodels/bpo_model3.h5')\nbpo_model4 = load_model('/kaggle/input/CAFAmodels/bpo_model4.h5')\nbpo_model5 = load_model('/kaggle/input/CAFAmodels/bpo_model5.h5')\nbpo_model6 = load_model('/kaggle/input/CAFAmodels/bpo_model6.h5')\n\ncco_model1 = load_model('/kaggle/input/CAFAmodels/cco_model1.h5')\ncco_model2 = load_model('/kaggle/input/CAFAmodels/cco_model2.h5')\ncco_model3 = load_model('/kaggle/input/CAFAmodels/cco_model3.h5')\ncco_model4 = load_model('/kaggle/input/CAFAmodels/cco_model4.h5')\ncco_model5 = load_model('/kaggle/input/CAFAmodels/cco_model5.h5')\ncco_model6 = load_model('/kaggle/input/CAFAmodels/cco_model6.h5')\n\nmfo_model1 = load_model('/kaggle/input/CAFAmodels/mfo_model1.h5')\nmfo_model2 = load_model('/kaggle/input/CAFAmodels/mfo_model2.h5')\nmfo_model3 = load_model('/kaggle/input/CAFAmodels/mfo_model3.h5')\nmfo_model4 = load_model('/kaggle/input/CAFAmodels/mfo_model4.h5')\nmfo_model5 = load_model('/kaggle/input/CAFAmodels/mfo_model5.h5')\nmfo_model6 = load_model('/kaggle/input/CAFAmodels/mfo_model6.h5')\n\n\nbpo_val_pred1 = bpo_model1.predict(bpo_val_X)\nbpo_val_pred2 = bpo_model2.predict(bpo_val_X)\nbpo_val_pred3 = bpo_model3.predict(bpo_val_X)\nbpo_val_pred4 = bpo_model4.predict(bpo_val_X)\nbpo_val_pred5 = bpo_model5.predict(bpo_val_X)\nbpo_val_pred6 = bpo_model6.predict(bpo_val_X)\n\ncco_val_pred1 = cco_model1.predict(cco_val_X)\ncco_val_pred2 = cco_model2.predict(cco_val_X)\ncco_val_pred3 = cco_model3.predict(cco_val_X)\ncco_val_pred4 = cco_model4.predict(cco_val_X)\ncco_val_pred5 = cco_model5.predict(cco_val_X)\ncco_val_pred6 = cco_model6.predict(cco_val_X)\n\nmfo_val_pred1 = mfo_model1.predict(mfo_val_X)\nmfo_val_pred2 = mfo_model2.predict(mfo_val_X)\nmfo_val_pred3 = mfo_model3.predict(mfo_val_X)\nmfo_val_pred4 = mfo_model4.predict(mfo_val_X)\nmfo_val_pred5 = mfo_model5.predict(mfo_val_X)\nmfo_val_pred6 = mfo_model6.predict(mfo_val_X)\n'''","metadata":{"execution":{"iopub.status.busy":"2023-06-18T17:10:34.038665Z","iopub.execute_input":"2023-06-18T17:10:34.040065Z","iopub.status.idle":"2023-06-18T17:12:14.140222Z","shell.execute_reply.started":"2023-06-18T17:10:34.039996Z","shell.execute_reply":"2023-06-18T17:12:14.138895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def auc_roc(predictions, val_df, aspect):\n    labels = val_df.columns\n    for model_num, prediction in enumerate(predictions, start=1):\n        auc_dict = {'aspect': aspect, 'model': model_num}\n        prediction_df = pd.DataFrame(data = prediction, columns = [val_df.columns]).reset_index(drop=True)\n        #roc_auc_score fails for top label as it only has one class (1) - so we skip it\n        for label in labels[1:]:  \n            auc_dict[label] = roc_auc_score(val_df[label], prediction_df[label])\n        try:\n            model_auc_df = model_auc_df.append(pd.DataFrame(auc_dict, index=[model_num-1]), ignore_index=True)\n        except NameError:\n            model_auc_df = pd.DataFrame(auc_dict, index=[0])\n    return model_auc_df\n\nbpo_auc_df = auc_roc([bpo_val_pred1, bpo_val_pred3, bpo_val_pred4, bpo_val_pred5, bpo_val_pred6], bpo_val_Y, 'BPO')\ncco_auc_df = auc_roc([cco_val_pred1, cco_val_pred3, cco_val_pred4, cco_val_pred5, cco_val_pred6], cco_val_Y, 'CCO')\nmfo_auc_df = auc_roc([mfo_val_pred1, mfo_val_pred3, mfo_val_pred4, mfo_val_pred5, mfo_val_pred6], mfo_val_Y, 'MFO')\n\n# write bpo_model_auc_df, cco_model_auc_df, mfo_model_auc_df to excel in sheets\n#with pd.ExcelWriter('model_auc_roc v1.xlsx') as writer:  \n#    bpo_auc_df.to_excel(writer, sheet_name='BPO')\n#    cco_auc_df.to_excel(writer, sheet_name='CCO')\n#    mfo_auc_df.to_excel(writer, sheet_name='MFO')\n","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:32:13.271157Z","iopub.execute_input":"2023-06-20T14:32:13.272002Z","iopub.status.idle":"2023-06-20T14:33:04.792739Z","shell.execute_reply.started":"2023-06-20T14:32:13.271963Z","shell.execute_reply":"2023-06-20T14:33:04.791686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def best_model(model_auc_df):\n    #get the best model name for each label (skip the top GO Term as no 0s in valuation data)\n    best_model_df = model_auc_df.iloc[:, 2:].idxmax()\n    best_model_df = pd.DataFrame({'GO_Term':best_model_df.index, 'best_model':best_model_df.values})\n    best_model_df['best_model'] = best_model_df['best_model'].apply(lambda x: model_auc_df['model'][x])\n    best_model_df['aspect'] = model_auc_df['aspect'][0]\n    return best_model_df\n\nbpo_best_model_df = best_model(bpo_auc_df)\ncco_best_model_df = best_model(cco_auc_df)\nmfo_best_model_df = best_model(mfo_auc_df)\n\n#combine all best_model_df into one\nbest_model_df = pd.concat([bpo_best_model_df, cco_best_model_df, mfo_best_model_df], ignore_index=True)\nprint(best_model_df.shape)\n#group by best_model, aspect and count the number of GO Terms for each model\nprint(best_model_df.groupby(['best_model', 'aspect']).count())","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:33:11.561828Z","iopub.execute_input":"2023-06-20T14:33:11.562404Z","iopub.status.idle":"2023-06-20T14:33:11.607154Z","shell.execute_reply.started":"2023-06-20T14:33:11.56236Z","shell.execute_reply":"2023-06-20T14:33:11.605242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Manual Evaluation\nTo ensure the accuracy of our AUC score rankings, it would be beneficial to manually check a few sample by calculating precision, recall, and F1 scores and compare if manual analysis yields to same conclusion on best model for that sample.","metadata":{}},{"cell_type":"code","source":"#take random sample of 2 GO Terms from best_model_df for each best_model to manually analyse the predictions\nsample = 2 #change number as required\nsample_df = best_model_df.groupby('best_model').apply(lambda x: x.sample(n=min(len(x), sample))).reset_index(drop=True)\n#insert top 2 GO Terms from bpo_val_Y, cco_val_Y, mfo_val_Y into sample_df\nsample_df = sample_df.append(pd.DataFrame({'GO_Term':bpo_val_Y.columns[0:2].tolist(), 'best_model':0, 'aspect':'BPO'}), ignore_index=True)\nsample_df = sample_df.append(pd.DataFrame({'GO_Term':cco_val_Y.columns[0:2].tolist(), 'best_model':0, 'aspect':'CCO'}), ignore_index=True)\nsample_df = sample_df.append(pd.DataFrame({'GO_Term':mfo_val_Y.columns[0:2].tolist(), 'best_model':0, 'aspect':'MFO'}), ignore_index=True)\n#run unique on sample_df to remove duplicates\nsample_df = sample_df.drop_duplicates(subset=['GO_Term', 'best_model', 'aspect'], keep='first').reset_index(drop=True)\n#assign best_model value from best_model_df to sample_df for the 2nd GO Terms (we do not have best model rank for the top GO Term)\nsample_df.loc[sample_df['GO_Term'].isin(bpo_val_Y.columns[1:2].tolist()), 'best_model'] = bpo_best_model_df.loc[bpo_best_model_df['GO_Term'].isin(bpo_val_Y.columns[1:2].tolist()), 'best_model'].tolist()\nsample_df.loc[sample_df['GO_Term'].isin(cco_val_Y.columns[1:2].tolist()), 'best_model'] = cco_best_model_df.loc[cco_best_model_df['GO_Term'].isin(cco_val_Y.columns[1:2].tolist()), 'best_model'].tolist()\nsample_df.loc[sample_df['GO_Term'].isin(mfo_val_Y.columns[1:2].tolist()), 'best_model'] = mfo_best_model_df.loc[mfo_best_model_df['GO_Term'].isin(mfo_val_Y.columns[1:2].tolist()), 'best_model'].tolist()\n\nprint(sample_df)\nGO_Terms_sample = sample_df['GO_Term'].tolist()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:33:18.666631Z","iopub.execute_input":"2023-06-20T14:33:18.667104Z","iopub.status.idle":"2023-06-20T14:33:18.707858Z","shell.execute_reply.started":"2023-06-20T14:33:18.667067Z","shell.execute_reply":"2023-06-20T14:33:18.706448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color='red'>Please note, the execution of the following cell may take a few minutes.</font>\n\nThe cell will generate a score chart for each sample GO Term x 5 models. To limit the number of plots, we've utilized a  (number_of_plots) to display the plot only for the top GO Term for each aspect. This variable can be modified if we need to view all the plots. Regardless, it generates scores for all the sample GO Terms across all models.\n\nMy objective is to identify a point where the highest F1 score aligns with precision, recall, and F1 scores that are in close proximity, thereby ensuring a balanced model performance. The standard deviation of F1, precision (P), and recall (R) is utilized to pinpoint the optimal threshold. Following this, the difference (f1 - StdDev) and the threshold values are assessed to determine which model yields the most satisfactory results.","metadata":{}},{"cell_type":"code","source":"#run model_rank to analyse predictions for required GO Terms\n#the function assumes prediction data sets are passed in an order of model numbers\ndef model_rank(predictions, Y, aspect):\n    labels = Y.columns\n    plot_data = []\n    #iterate for GO_Terms_sample\n    for model_num, prediction in enumerate(predictions, start=1):\n        predictions_df = pd.DataFrame(data = prediction, columns = Y.columns).reset_index(drop=True)\n        predictions_df = predictions_df.add_prefix('pred_')\n        pred_labels = predictions_df.columns\n        predictions_df = Y.join(predictions_df)\n        \n        #for Go_Term in GO_Terms_sample:\n        for Go_Term in labels:  #run this to override AUC score by manual scores\n            if Go_Term in Y.columns:\n                pred_GoTerm = 'pred_' + Go_Term\n                #join labels and predictions_df by index\n                thrs = np.arange(0.2, 0.9, 0.01)\n                #create empty dataframe to store precision, recall and f1 score\n                score_df = pd.DataFrame(columns = ['aspect', 'model_number', 'thr', 'precision', 'recall', 'f1', 'positives', 'negatives', 'true_positives', 'true_negatives', 'false_positives', 'false_negatives'])\n                for thr in thrs:\n                    thr_column = 'thr_' + str(thr) + '_result'\n                    predictions_df[thr_column] = np.where(predictions_df[pred_GoTerm] >= thr, 1., 0.)\n                    #calculate precision, recall and f1 score\n                    precision = precision_score(predictions_df[Go_Term], predictions_df[thr_column])\n                    recall = recall_score(predictions_df[Go_Term], predictions_df[thr_column])\n                    f1 = f1_score(predictions_df[Go_Term], predictions_df[thr_column])\n                    positives = predictions_df[thr_column].sum()\n                    negatives = predictions_df.shape[0] - positives\n                    true_positives = ((predictions_df[Go_Term] == 1) & (predictions_df[thr_column] == 1)).sum()\n                    true_negatives = ((predictions_df[Go_Term] == 0) & (predictions_df[thr_column] == 0)).sum()\n                    false_positives = ((predictions_df[Go_Term] == 0) & (predictions_df[thr_column] == 1)).sum()\n                    false_negatives = ((predictions_df[Go_Term] == 1) & (predictions_df[thr_column] == 0)).sum()\n                    #calculate standard deviation to find best score where precision, recall and f1 score are close to each other\n                    std_dev = np.std(np.array([precision, recall, f1]))\n                    #find the threshold where f1 score is high and standard deviation is low using the difference between f1 score and standard deviation\n                    diff_f1_sd = f1 - std_dev\n                    #concat thr, precision, recall and f1 score to score_df\n                    score_df = pd.concat([score_df, pd.DataFrame([[aspect, model_num, Go_Term, thr, precision, recall, f1, std_dev, diff_f1_sd, positives, negatives, true_positives, true_negatives, false_positives, false_negatives]], \n                                                        columns = ['aspect', 'model_number', 'GO_Term', 'thr', 'precision', 'recall', 'f1', 'std_dev', 'diff_f1_sd', 'positives', 'negatives', 'true_positives', 'true_negatives', 'false_positives', 'false_negatives'])])\n\n                score_df = score_df.reset_index(drop=True)        \n                #identify best score row\n                max_score_row = score_df.loc[score_df['diff_f1_sd'].idxmax()]\n                # check if there are multiple rows with the same max score\n                if (score_df['diff_f1_sd'] == max_score_row['diff_f1_sd']).sum() > 1:\n                    # filter dataframe to rows with max score only\n                    df_max_score = score_df[score_df['diff_f1_sd'] == max_score_row['diff_f1_sd']]\n                    # get row with max 'thr' within rows having max 'score'\n                    max_thr_idx = df_max_score['thr'].idxmax()\n                    best_model_score = score_df.loc[max_thr_idx]\n                else:\n                    best_model_score = max_score_row\n\n                #add best_model_score to model_score_df\n                try:\n                    model_score_df = model_score_df.append(best_model_score.to_frame().T)\n                except NameError:\n                    model_score_df = best_model_score.to_frame().T\n\n                plot_data.append((score_df, aspect, Go_Term, model_num))\n\n    # Sort data by Go_Term and model_num\n    Go_Term_Order = {item: index for index, item in enumerate(labels)}\n    plot_data_sorted = sorted(plot_data, key=lambda x: (Go_Term_Order.get(x[2], len(labels)), x[3]))\n\n    #plot precision, recall and f1 score against threshold\n    #plots displayed =  number of GO Terms x 6 models \n    number_of_plots = 5 # Restricting the display to include only one GO_Term plot per aspect\n    for data in plot_data_sorted[:number_of_plots]:\n        score_df, aspect, Go_Term, model_num = data\n        fig = plt.figure(figsize=(18, 3))\n        plt.plot(score_df['thr'], score_df['precision'], label = 'Precision')\n        plt.plot(score_df['thr'], score_df['recall'], label = 'Recall')\n        plt.plot(score_df['thr'], score_df['f1'], label = 'F1 Score')\n        plt.plot(score_df['thr'], score_df['std_dev'], label = 'std_dev', linestyle='dotted')\n        plt.plot(score_df['thr'], score_df['diff_f1_sd'], label = 'diff_f1_sd', linestyle='--')\n        plt.xticks(score_df['thr'])\n        plt.xlabel('Threshold')\n        plt.ylabel('Score')\n        if model_num == 1:\n            plt.title(aspect + ' Model-' + str(model_num) + ' Scores against Threshold for ' + Go_Term)\n        else:\n            plt.title(aspect + ' Model-' + str(model_num+1) + ' Scores against Threshold for ' + Go_Term)\n        plt.legend()\n        plt.show()\n\n    #sort model_score_df by GO_Term in the order of it's occurance rank\n    model_score_df['GO_Term_Order'] = model_score_df['GO_Term'].map(Go_Term_Order)\n    model_score_df = model_score_df.sort_values(by=['GO_Term_Order', 'model_number']).reset_index(drop=True)\n    model_score_df = model_score_df.drop(columns=['GO_Term_Order'])\n    return model_score_df\n\nbpo_model_score_df = model_rank([bpo_val_pred1, bpo_val_pred3, bpo_val_pred4, bpo_val_pred5, bpo_val_pred6], bpo_val_Y, 'BPO')\ncco_model_score_df = model_rank([cco_val_pred1, cco_val_pred3, cco_val_pred4, cco_val_pred5, cco_val_pred6], cco_val_Y, 'CCO')\nmfo_model_score_df = model_rank([mfo_val_pred1, mfo_val_pred3, mfo_val_pred4, mfo_val_pred5, cco_val_pred6], mfo_val_Y, 'MFO')\nthr_score_df = pd.concat([bpo_model_score_df, cco_model_score_df, mfo_model_score_df]).reset_index(drop=True)\n\n#write scores to excel\nthr_score_df.to_excel('model_thr_analysis.xlsx')\n","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:34:09.051228Z","iopub.execute_input":"2023-06-20T14:34:09.051953Z","iopub.status.idle":"2023-06-20T14:38:54.634881Z","shell.execute_reply.started":"2023-06-20T14:34:09.051904Z","shell.execute_reply":"2023-06-20T14:38:54.633407Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"*Upon reviewing the score sheet, the manual score analysis suggests alternative options for the best model for a few GO Terms. However, these scores are closely comparable to those chosen by the AUC method. As such, I have decided to not make any changes.*","metadata":{}},{"cell_type":"code","source":"'''\n#if this cell has to rerun, the best_model function cell has to rerun prior to this.\n\n#Assign best model for top GO Term in each aspect as there was no AUC rank for them, \n#model analysis data shows model [6] is best for all 3 aspects (Ignoring model 2 overfitting)\n#add row GO Term = 'GO:0008150' to best_model_df with best_model = 1\nbest_model_df = pd.concat([best_model_df, pd.DataFrame([['BPO', 5, 'GO:0008150']], columns = ['aspect', 'best_model', 'GO_Term'])])\n#add row GO Term = 'GO:0005575' to best_model_df with best_model = 1\nbest_model_df = pd.concat([best_model_df, pd.DataFrame([['CCO', 5, 'GO:0005575']], columns = ['aspect', 'best_model', 'GO_Term'])])\n#add row GO Term = 'GO:0003674' to best_model_df with best_model = 1\nbest_model_df = pd.concat([best_model_df, pd.DataFrame([['MFO', 5, 'GO:0003674']], columns = ['aspect', 'best_model', 'GO_Term'])])\n'''\n#Override best model df\n#get best thr row for each GO Term\nthr_score_df['thr'] = thr_score_df['thr'].astype(float)\nbest_thr_score_df = thr_score_df.loc[thr_score_df.groupby(['GO_Term'])['thr'].idxmax()].reset_index(drop=True)\nbest_model_df = best_thr_score_df[['aspect', 'GO_Term', 'model_number']].copy()\nbest_model_df = best_model_df.rename(columns={'model_number': 'best_model'})\n\n#check the count by aspect and best_model\nprint(best_model_df.groupby(['best_model', 'aspect']).count())\nprint(best_model_df.shape, ' The number of GO Terms should be (number_of_labels x 3)')","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:40:05.108872Z","iopub.execute_input":"2023-06-20T14:40:05.109858Z","iopub.status.idle":"2023-06-20T14:40:05.133253Z","shell.execute_reply.started":"2023-06-20T14:40:05.109818Z","shell.execute_reply":"2023-06-20T14:40:05.131731Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Release memory\ndel train_terms, bpo_train_terms, cco_val, bpo_val, mfo_val, cco_train_terms, mfo_train_terms, train_protein_ids\ndel bpo_train_X, bpo_val_X, bpo_val_Y, bpo_val_pred1, bpo_val_pred3, bpo_val_pred4, bpo_val_pred5, bpo_val_pred6\ndel cco_train_X, cco_val_X, cco_val_Y, cco_val_pred1, cco_val_pred3, cco_val_pred4, cco_val_pred5, cco_val_pred6\ndel mfo_train_X, mfo_val_X, mfo_val_Y, mfo_val_pred1, mfo_val_pred3, mfo_val_pred4, mfo_val_pred5, mfo_val_pred6\ndel thr_score_df, best_thr_score_df\n\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:40:20.552798Z","iopub.execute_input":"2023-06-20T14:40:20.553231Z","iopub.status.idle":"2023-06-20T14:40:21.153823Z","shell.execute_reply.started":"2023-06-20T14:40:20.553197Z","shell.execute_reply":"2023-06-20T14:40:21.152534Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Ready to run required best_models for each aspect and prepare submission","metadata":{}},{"cell_type":"code","source":"#load test embeddings of protein sequences\ntest_embeddings = np.load('/kaggle/input/t5embeds/test_embeds.npy')\n# Convert test_embeddings to dataframe\ncolumn_num = test_embeddings.shape[1]\ntest_df = pd.DataFrame(test_embeddings, columns = [\"Column_\" + str(i) for i in range(1, column_num+1)])\nprint(test_df.shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:40:29.533963Z","iopub.execute_input":"2023-06-20T14:40:29.534395Z","iopub.status.idle":"2023-06-20T14:40:42.567787Z","shell.execute_reply.started":"2023-06-20T14:40:29.534359Z","shell.execute_reply":"2023-06-20T14:40:42.56655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"<font color='red'>*Note: Running multiple predictions below for the entire test set and we select columns as required from the prediction result*\n    This may take 10-15 minutes","metadata":{}},{"cell_type":"code","source":"#Run predictions as required by best_model conclusion\n# 3 aspects model dictionary {'BPO' : {1: bpo_model_1, 2: bpo_model_2}}\nmodel_dict = {'BPO' : {1: bpo_model1,  2: bpo_model3, 3: bpo_model4, 4: bpo_model5, 5: bpo_model6},\n              'CCO' : {1: cco_model1,  2: cco_model3, 3: cco_model4, 4: cco_model5, 5: cco_model6},\n              'MFO' : {1: mfo_model1,  2: mfo_model3, 3: mfo_model4, 4: mfo_model5, 5: mfo_model6}}\n\nprediction_columns = {'BPO' : bpo_train_Y.columns, 'CCO' : cco_train_Y.columns, 'MFO' : mfo_train_Y.columns}\n\n#Get 'aspect' and 'best_model' columns from best_model_df, unique and sort them\naspect_model = best_model_df[['aspect', 'best_model']].drop_duplicates().sort_values(by=['aspect', 'best_model']).reset_index(drop=True)\n#iterate through best_model_df and predict the labels for test data\npredictions = []\npredictions_df = []\nfor index, row in aspect_model.iterrows():\n    print('Predicting for aspect: ', row['aspect'], ' and best_model: ', row['best_model'])\n    predictions.append(model_dict[row['aspect']][row['best_model']].predict(test_df))\n    predictions_df.append(pd.DataFrame(predictions[index], columns = prediction_columns[row['aspect']]))\n\nprint('Predictions_df shape: ', predictions_df[0].shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:42:00.203925Z","iopub.execute_input":"2023-06-20T14:42:00.204696Z","iopub.status.idle":"2023-06-20T14:49:04.568Z","shell.execute_reply.started":"2023-06-20T14:42:00.204659Z","shell.execute_reply":"2023-06-20T14:49:04.566769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"del predictions, test_embeddings, test_df\ndel cco_train_Y, bpo_train_Y, mfo_train_Y\ngc.collect()","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:49:13.949479Z","iopub.execute_input":"2023-06-20T14:49:13.950378Z","iopub.status.idle":"2023-06-20T14:49:15.74262Z","shell.execute_reply.started":"2023-06-20T14:49:13.950329Z","shell.execute_reply":"2023-06-20T14:49:15.741203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#select columns from predictions_df as per best_model_df GO_Terms where aspect = row['aspect'] and best_model = row['best_model']\nfor index, row in aspect_model.iterrows():\n    predictions_df[index] = predictions_df[index][best_model_df.loc[(best_model_df['aspect'] == row['aspect']) & (best_model_df['best_model'] == row['best_model']), 'GO_Term']]\n\n#concatenate all predictions_df\nfinal_predictions_df = pd.concat(predictions_df, axis=1)\nprint(final_predictions_df.shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:49:52.444345Z","iopub.execute_input":"2023-06-20T14:49:52.445157Z","iopub.status.idle":"2023-06-20T14:49:53.666689Z","shell.execute_reply.started":"2023-06-20T14:49:52.445105Z","shell.execute_reply":"2023-06-20T14:49:53.665201Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#generate and save submission to tsv file\n#Function Help from https://www.kaggle.com/code/gusthema/cafa-5-protein-function-with-tensorflow\ndf_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] * num_of_labels * 3\n\nlabels = list(final_predictions_df.columns)\n\ndf_submission['Protein Id'] = l\ndf_submission['GO Term Id'] = labels * final_predictions_df.shape[0]\ndf_submission['Prediction'] = final_predictions_df.values.ravel()\nprint(df_submission.head())\nprint(df_submission.shape, ' - Must have 141865 * ',  num_of_labels, ' * 3 = ', 141865 * num_of_labels * 3,  ' records')\ndf_submission.to_csv(\"submission.tsv\",header=False, index=False, sep=\"\\t\")","metadata":{"execution":{"iopub.status.busy":"2023-06-20T14:58:01.45091Z","iopub.execute_input":"2023-06-20T14:58:01.451582Z","iopub.status.idle":"2023-06-20T15:08:30.016648Z","shell.execute_reply.started":"2023-06-20T14:58:01.45153Z","shell.execute_reply":"2023-06-20T15:08:30.015309Z"},"trusted":true},"execution_count":null,"outputs":[]}]}