{"metadata":{"kernelspec":{"name":"python3","display_name":"Python 3","language":"python"},"language_info":{"name":"python","version":"3.11.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":116062,"databundleVersionId":14084779,"sourceType":"competition"}],"dockerImageVersionId":31153,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# CAFA 6 Protein Function Prediction Challenge\nIn this notebook, we'll...\n- do some light EDA on the data files for [this competition](https://www.kaggle.com/competitions/cafa-6-protein-function-prediction)\n- generate protein descriptors for the FASTA sequences\n- build a simple model for predicting GO terms\n\nLet's do this!","metadata":{}},{"cell_type":"markdown","source":"-------\n-------\n-------\n# Install Necessary Libraries\nThe following libraries are for reading the data files (FASTA files and Gene Ontology graph) and for feature engineering (protein descriptors).","metadata":{}},{"cell_type":"code","source":"!pip install Bio goatools propy3","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:29:38.165741Z","iopub.execute_input":"2025-12-23T21:29:38.166120Z","iopub.status.idle":"2025-12-23T21:30:23.413638Z","shell.execute_reply.started":"2025-12-23T21:29:38.166095Z","shell.execute_reply":"2025-12-23T21:30:23.412359Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"-------\n-------\n-------\n# Imports","metadata":{}},{"cell_type":"code","source":"import os\nimport numpy as np\nimport pandas as pd\nfrom datetime import datetime\nimport time\nimport matplotlib.pyplot as plt\nimport pickle\n\nfrom scipy.sparse import csr_matrix\n\nfrom Bio import SeqIO\nfrom goatools.obo_parser import GODag\nfrom propy import PyPro\nfrom propy.GetProteinFromUniprot import GetProteinSequence","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:30:23.414893Z","iopub.execute_input":"2025-12-23T21:30:23.415197Z","iopub.status.idle":"2025-12-23T21:30:24.700688Z","shell.execute_reply.started":"2025-12-23T21:30:23.415170Z","shell.execute_reply":"2025-12-23T21:30:24.699445Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"-------\n-------\n-------\n# Observe File Contents (Terminal)\nHere, we'll have a look at the first few lines of each file to know what we are dealing with.\n\nFirst, let's see the competition's files.","metadata":{}},{"cell_type":"code","source":"competition_folder = '/kaggle/input/cafa-6-protein-function-prediction/'\nlist_of_filenames = []\nfor dirname, _, filenames in os.walk(competition_folder):\n    for filename in filenames:\n        filepath = os.path.join(dirname, filename)\n        list_of_filenames.append(filepath)\n        print(filepath)","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","execution":{"iopub.status.busy":"2025-12-23T21:31:39.593770Z","iopub.execute_input":"2025-12-23T21:31:39.594506Z","iopub.status.idle":"2025-12-23T21:31:39.619736Z","shell.execute_reply.started":"2025-12-23T21:31:39.594452Z","shell.execute_reply":"2025-12-23T21:31:39.618621Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Now, for the ***train set*** files:\n- FASTA sequences\n- GO terms for each of the sequencecs\n- Taxonomy (i.e., species)","metadata":{}},{"cell_type":"code","source":"# FASTA sequences (train)\ntrain_sequences_file = competition_folder + '/Train/train_sequences.fasta'\n!head {train_sequences_file} -n 10","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:31:42.398675Z","iopub.execute_input":"2025-12-23T21:31:42.399012Z","iopub.status.idle":"2025-12-23T21:31:42.532675Z","shell.execute_reply.started":"2025-12-23T21:31:42.398989Z","shell.execute_reply":"2025-12-23T21:31:42.530859Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# GO terms (train)\ntrain_terms_file = competition_folder + '/Train/train_terms.tsv'\n!head {train_terms_file} -n 10","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:31:43.349460Z","iopub.execute_input":"2025-12-23T21:31:43.349844Z","iopub.status.idle":"2025-12-23T21:31:43.484388Z","shell.execute_reply.started":"2025-12-23T21:31:43.349800Z","shell.execute_reply":"2025-12-23T21:31:43.482628Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# Taxonomy (train)\ntrain_taxonomy_file = competition_folder + '/Train/train_taxonomy.tsv'\n!head {train_taxonomy_file} -n 10","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:31:58.102608Z","iopub.execute_input":"2025-12-23T21:31:58.102969Z","iopub.status.idle":"2025-12-23T21:31:58.240166Z","shell.execute_reply.started":"2025-12-23T21:31:58.102919Z","shell.execute_reply":"2025-12-23T21:31:58.239011Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Now, let's look at the file containing the Gene Ontology (GO) graph.","metadata":{}},{"cell_type":"code","source":"# Gene Ontology\ngo_file = competition_folder + '/Train/go-basic.obo'\n!head {go_file} -n 55","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:02.847189Z","iopub.execute_input":"2025-12-23T21:32:02.847561Z","iopub.status.idle":"2025-12-23T21:32:02.982822Z","shell.execute_reply.started":"2025-12-23T21:32:02.847531Z","shell.execute_reply":"2025-12-23T21:32:02.981321Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Below is the files containing ***test set*** sequences.","metadata":{}},{"cell_type":"code","source":"# FASTA sequences (test)\ntestsuperset_file = competition_folder + '/Test/testsuperset.fasta'\n!head {testsuperset_file} -n 20","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:06.745081Z","iopub.execute_input":"2025-12-23T21:32:06.745441Z","iopub.status.idle":"2025-12-23T21:32:06.880816Z","shell.execute_reply.started":"2025-12-23T21:32:06.745412Z","shell.execute_reply":"2025-12-23T21:32:06.879464Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"The file below shows what each of the species IDs means.","metadata":{}},{"cell_type":"code","source":"# Species lookup\ntestsuperset_taxon_file = competition_folder + '/Test/testsuperset-taxon-list.tsv'\n!head {testsuperset_taxon_file} -n 10","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:07.977736Z","iopub.execute_input":"2025-12-23T21:32:07.978106Z","iopub.status.idle":"2025-12-23T21:32:08.112093Z","shell.execute_reply.started":"2025-12-23T21:32:07.978076Z","shell.execute_reply":"2025-12-23T21:32:08.110731Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"No idea what the file below is for.  \n<left>¯\\\\__ (ツ) __/¯</left>","metadata":{}},{"cell_type":"code","source":"# Information accretion (weights) for each GO term --> used to compute weighted precision and recall during Evaluation\nia_file = competition_folder + '/IA.tsv'\n!head {ia_file} -n 10","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:09.877433Z","iopub.execute_input":"2025-12-23T21:32:09.878513Z","iopub.status.idle":"2025-12-23T21:32:10.011475Z","shell.execute_reply.started":"2025-12-23T21:32:09.878462Z","shell.execute_reply":"2025-12-23T21:32:10.010171Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"And finally, what a sample submission looks like.","metadata":{}},{"cell_type":"code","source":"sample_submission_file = competition_folder + '/sample_submission.tsv'\n!head {sample_submission_file} -n 20","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:11.910430Z","iopub.execute_input":"2025-12-23T21:32:11.910809Z","iopub.status.idle":"2025-12-23T21:32:12.046925Z","shell.execute_reply.started":"2025-12-23T21:32:11.910771Z","shell.execute_reply":"2025-12-23T21:32:12.045678Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"-------\n-------\n-------\n# EDA\nWe shall load the different files, get some stats on the data, and observe some examples.","metadata":{}},{"cell_type":"markdown","source":"## Species Lookup\nSpecies IDs and their corresponding names","metadata":{}},{"cell_type":"code","source":"# species lookup table\n\nspecies_lookup_df = pd.read_csv(testsuperset_taxon_file, sep=\"\\t\")\nspecies_lookup_df.columns = ['species_id', 'species_name']\nspecies_lookup_df","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:14.710761Z","iopub.execute_input":"2025-12-23T21:32:14.711139Z","iopub.status.idle":"2025-12-23T21:32:14.774272Z","shell.execute_reply.started":"2025-12-23T21:32:14.711109Z","shell.execute_reply":"2025-12-23T21:32:14.773221Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"## FASTA Sequences\nProtein sequences in the train and test sets","metadata":{}},{"cell_type":"code","source":"# load the files containing the FASTA sequences\ntrain_fasta_sequences = [{\"id\": record.id, \"seq\": record.seq} for record in SeqIO.parse(train_sequences_file, \"fasta\")]\ntest_fasta_sequences  = [{\"id\": record.id, \"seq\": record.seq} for record in SeqIO.parse(testsuperset_file,    \"fasta\")]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:15.116549Z","iopub.execute_input":"2025-12-23T21:32:15.116876Z","iopub.status.idle":"2025-12-23T21:32:19.510668Z","shell.execute_reply.started":"2025-12-23T21:32:15.116854Z","shell.execute_reply":"2025-12-23T21:32:19.509706Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# number of FASTA sequences in train & test sets\nprint(f\"len( train_fasta_sequences ) = {len(train_fasta_sequences)}\")\nprint(f\"len(  test_fasta_sequences ) = {len(test_fasta_sequences)}\")","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:19.512088Z","iopub.execute_input":"2025-12-23T21:32:19.512461Z","iopub.status.idle":"2025-12-23T21:32:19.517986Z","shell.execute_reply.started":"2025-12-23T21:32:19.512438Z","shell.execute_reply":"2025-12-23T21:32:19.516999Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# look at data sample (first 3 examples of train set)\nfor record in train_fasta_sequences[:3]:\n    print(record['id'])\n    print(record['seq'])\n    print('')","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:19.518882Z","iopub.execute_input":"2025-12-23T21:32:19.519201Z","iopub.status.idle":"2025-12-23T21:32:19.537475Z","shell.execute_reply.started":"2025-12-23T21:32:19.519171Z","shell.execute_reply":"2025-12-23T21:32:19.536472Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# plot distribution of sequence lengths in train set\nlengths = [len(str(entry[\"seq\"])) for entry in train_fasta_sequences]\n\nbin_width = 100\nbins = np.arange(0, max(lengths) + bin_width, bin_width)\nplt.figure()\nplt.hist(lengths, bins = bins)\nplt.xlabel(\"Sequence Length\")\nplt.ylabel(\"Count\")\nplt.title(\"Distribution of FASTA Sequence Lengths (TRAIN)\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:19.539609Z","iopub.execute_input":"2025-12-23T21:32:19.539989Z","iopub.status.idle":"2025-12-23T21:32:20.463607Z","shell.execute_reply.started":"2025-12-23T21:32:19.539937Z","shell.execute_reply":"2025-12-23T21:32:20.462594Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Almost all sequences are below 5K amino acids long. There are a few very long sequences as well (up to 35K amino acids).\n\nLet us have a look at some numeric stats as well as the 90th and 99th percentile.","metadata":{}},{"cell_type":"code","source":"# stats on sequence lengths in train set\nprint(pd.Series(lengths).describe())\nprint('')\np80 = np.percentile(lengths, 80)\nprint(\"80th percentile sequence length:\", p80)\np90 = np.percentile(lengths, 90)\nprint(\"90th percentile sequence length:\", p90)\np99 = np.percentile(lengths, 99)\nprint(\"99th percentile sequence length:\", p99)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:20.464481Z","iopub.execute_input":"2025-12-23T21:32:20.464735Z","iopub.status.idle":"2025-12-23T21:32:20.526126Z","shell.execute_reply.started":"2025-12-23T21:32:20.464713Z","shell.execute_reply":"2025-12-23T21:32:20.525149Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Some of the proteins are really tiny (< 10 amino acids). Guess in that case, they should be called peptides, not proteins. :)\n\nFurthermore, 90% of all proteins are < 1K amino acids long.\n\nNow let's do the same but for the test set!","metadata":{}},{"cell_type":"code","source":"# stats and distribution of sequence lengths in test set\nlengths = [len(str(entry[\"seq\"])) for entry in test_fasta_sequences]\n\nprint(pd.Series(lengths).describe())\nprint('')\np80 = np.percentile(lengths, 80)\nprint(\"80th percentile sequence length:\", p80)\np90 = np.percentile(lengths, 90)\nprint(\"90th percentile sequence length:\", p90)\np99 = np.percentile(lengths, 99)\nprint(\"99th percentile sequence length:\", p99)\n\nbin_width = 100\nbins = np.arange(0, max(lengths) + bin_width, bin_width)\nplt.figure()\nplt.hist(lengths, bins = bins)\nplt.xlabel(\"Sequence Length\")\nplt.ylabel(\"Count\")\nplt.title(\"Distribution of FASTA Sequence Lengths (TEST)\")\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:20.527306Z","iopub.execute_input":"2025-12-23T21:32:20.527721Z","iopub.status.idle":"2025-12-23T21:32:21.900022Z","shell.execute_reply.started":"2025-12-23T21:32:20.527693Z","shell.execute_reply":"2025-12-23T21:32:21.898762Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Stats and distribution of sequence lengths in test set are similar to those of the train set.","metadata":{}},{"cell_type":"markdown","source":"## Train Set GO Terms\nFor each train set sequence, which GO terms are assigned to it?","metadata":{}},{"cell_type":"code","source":"# load train set terms\ntrain_terms_df = pd.read_csv(train_terms_file, sep=\"\\t\")\ntrain_terms_df","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:29.767153Z","iopub.execute_input":"2025-12-23T21:32:29.768095Z","iopub.status.idle":"2025-12-23T21:32:30.124390Z","shell.execute_reply.started":"2025-12-23T21:32:29.768060Z","shell.execute_reply":"2025-12-23T21:32:30.123422Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# number of unique values in each column\nprint(f\"Number of unique entries: {train_terms_df.EntryID.nunique()}\")\nprint(f\"Number of unique terms:   {train_terms_df.term.nunique()}\")\nprint(f\"Number of unique aspects: {train_terms_df.aspect.nunique()}\")","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:30.125786Z","iopub.execute_input":"2025-12-23T21:32:30.126139Z","iopub.status.idle":"2025-12-23T21:32:30.271510Z","shell.execute_reply.started":"2025-12-23T21:32:30.126116Z","shell.execute_reply":"2025-12-23T21:32:30.270580Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"**From the cell above, It seems that what we have here is a *\"multi-label classification problem\"* involving 26,125 labels!**\n\nThat is, for each of the FASTA sequences, we are trying to predict which of the labels apply to it.\n\nNow, let's see how many unique terms there are for each aspect.","metadata":{}},{"cell_type":"code","source":"# total number of terms + number of unique terms per aspect\naspect_terms_df = train_terms_df[['aspect', 'term']].drop_duplicates().reset_index(drop=True)\nprint('Total number of terms:', len(aspect_terms_df))\nprint('')\nprint(aspect_terms_df['aspect'].value_counts())","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:30.273188Z","iopub.execute_input":"2025-12-23T21:32:30.273464Z","iopub.status.idle":"2025-12-23T21:32:30.414138Z","shell.execute_reply.started":"2025-12-23T21:32:30.273442Z","shell.execute_reply":"2025-12-23T21:32:30.413023Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Below is a sanity check to ensure that there are no common terms across aspects.","metadata":{}},{"cell_type":"code","source":"# NO intersections between terms of different aspects (as expected)\np_terms = aspect_terms_df[aspect_terms_df['aspect'] == \"P\"].term\nc_terms = aspect_terms_df[aspect_terms_df['aspect'] == \"C\"].term\nf_terms = aspect_terms_df[aspect_terms_df['aspect'] == \"F\"].term\n\nprint(set(p_terms).intersection(set(c_terms)))  #\nprint(set(p_terms).intersection(set(f_terms)))  # empty sets (i.e., no intersections)\nprint(set(c_terms).intersection(set(f_terms)))  #","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:32:30.672294Z","iopub.execute_input":"2025-12-23T21:32:30.672676Z","iopub.status.idle":"2025-12-23T21:32:30.699650Z","shell.execute_reply.started":"2025-12-23T21:32:30.672650Z","shell.execute_reply":"2025-12-23T21:32:30.698729Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Let's check now some on how many terms are assigned to each of the train sequences.","metadata":{}},{"cell_type":"code","source":"# stats on number of terms per train sequence\ntmp = train_terms_df.value_counts('EntryID')\nprint(tmp.describe())\nprint('')\np80 = np.percentile(tmp.values, 80)\nprint(\"80th percentile sequence length:\", p80)\np90 = np.percentile(tmp.values, 90)\nprint(\"90th percentile sequence length:\", p90)\np99 = np.percentile(tmp.values, 99)\nprint(\"99th percentile sequence length:\", p99)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:30.937252Z","iopub.execute_input":"2025-12-23T21:32:30.937549Z","iopub.status.idle":"2025-12-23T21:32:31.121161Z","shell.execute_reply.started":"2025-12-23T21:32:30.937528Z","shell.execute_reply":"2025-12-23T21:32:31.120210Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Proteins tend to have around 2 to 8 terms assigned to them. Almost all of them have under 50 terms.\n\nBefore we move on, let's have observe how many terms ***of each aspect*** each protein has assigned to it.","metadata":{}},{"cell_type":"code","source":"# counts of each aspect's terms per train sequence\nresult = (\n    train_terms_df.groupby(\"EntryID\")\n      .agg(\n          num_p = (\"aspect\", lambda x: (x == \"P\").sum()),\n          num_c = (\"aspect\", lambda x: (x == \"C\").sum()),\n          num_f = (\"aspect\", lambda x: (x == \"F\").sum())\n      )\n      .reset_index()\n)\n\nresult[\"num_terms\"] = result[\"num_p\"] + result[\"num_c\"] + result[\"num_f\"]\nresult[\"percent_p\"] = result[\"num_p\"] / result[\"num_terms\"] * 100\nresult[\"percent_c\"] = result[\"num_c\"] / result[\"num_terms\"] * 100\nresult[\"percent_f\"] = result[\"num_f\"] / result[\"num_terms\"] * 100\n\nresult","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:31.191538Z","iopub.execute_input":"2025-12-23T21:32:31.191906Z","iopub.status.idle":"2025-12-23T21:32:57.781657Z","shell.execute_reply.started":"2025-12-23T21:32:31.191879Z","shell.execute_reply":"2025-12-23T21:32:57.780516Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Train sequence with < 4 terms:\")\nprint(result[result['num_terms'] < 4].describe())\nprint('')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:57.783142Z","iopub.execute_input":"2025-12-23T21:32:57.783413Z","iopub.status.idle":"2025-12-23T21:32:57.826051Z","shell.execute_reply.started":"2025-12-23T21:32:57.783391Z","shell.execute_reply":"2025-12-23T21:32:57.825029Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Train sequence with 4 to 10 terms:\")\nprint(result[result['num_terms'].between(4, 10)].describe())\nprint('')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:57.826894Z","iopub.execute_input":"2025-12-23T21:32:57.827160Z","iopub.status.idle":"2025-12-23T21:32:57.864815Z","shell.execute_reply.started":"2025-12-23T21:32:57.827140Z","shell.execute_reply":"2025-12-23T21:32:57.863845Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print(\"Train sequence with at > 10 terms:\")\nprint(result[result['num_terms'] > 10].describe())\nprint('')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:32:57.866659Z","iopub.execute_input":"2025-12-23T21:32:57.866987Z","iopub.status.idle":"2025-12-23T21:32:57.896445Z","shell.execute_reply.started":"2025-12-23T21:32:57.866939Z","shell.execute_reply":"2025-12-23T21:32:57.895223Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# box plot of percentages of P, C and F terms per train sequence\n# (we show only sequences with 4 to 10 terms because otherwise the percentages are not very meaningful)\n# (1 to 3 is too few terms, and most proteins have less than 10 terms anyway)\np_terms_counts = result[result['num_terms'].between(4, 10)]['num_p'].value_counts('EntryID')\nc_terms_counts = result[result['num_terms'].between(4, 10)]['num_c'].value_counts('EntryID')\nf_terms_counts = result[result['num_terms'].between(4, 10)]['num_f'].value_counts('EntryID')\nplt.figure()\nplt.boxplot([p_terms_counts, c_terms_counts, f_terms_counts])\nplt.xticks([1, 2, 3], ['P', 'C', 'F'])\nplt.ylabel(\"Number of Terms per Sequence\")\n\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:35:41.470035Z","iopub.execute_input":"2025-12-23T21:35:41.470387Z","iopub.status.idle":"2025-12-23T21:35:41.622697Z","shell.execute_reply.started":"2025-12-23T21:35:41.470363Z","shell.execute_reply":"2025-12-23T21:35:41.621634Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"In conclusion...  \nIn most proteins, percentages of P, C and F terms tend to be similar. However, when there are MANY terms (> 50), percentage of P terms tends to significantly higher than the other two.","metadata":{}},{"cell_type":"markdown","source":"## Train Set Taxonomies (Species)\nBy how many proteins is each species represented in the train set?","metadata":{}},{"cell_type":"code","source":"# species id for each fasta record in the train set\ntrain_species_df = pd.read_csv(train_taxonomy_file, sep=\"\\t\", header=None)\ntrain_species_df.columns = ['fasta_id', 'species_id']\ntrain_species_df","metadata":{"execution":{"iopub.status.busy":"2025-12-23T21:36:22.268409Z","iopub.execute_input":"2025-12-23T21:36:22.268738Z","iopub.status.idle":"2025-12-23T21:36:22.337432Z","shell.execute_reply.started":"2025-12-23T21:36:22.268715Z","shell.execute_reply":"2025-12-23T21:36:22.336505Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_species_df = pd.merge(train_species_df, species_lookup_df, on='species_id', how='inner')\ntrain_species_df[['species_id', 'species_name']].value_counts()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:23.179893Z","iopub.execute_input":"2025-12-23T21:36:23.180256Z","iopub.status.idle":"2025-12-23T21:36:23.223900Z","shell.execute_reply.started":"2025-12-23T21:36:23.180232Z","shell.execute_reply":"2025-12-23T21:36:23.223002Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"Humans are the most represented in the (train) data, followed by mice. Hardly surprising.","metadata":{}},{"cell_type":"markdown","source":"## Gene Ontology DAG\nGene Ontology is represented as a directed acyclic graph showing hierarchy of biological processes (P), cellular components (C) and molecular functions (F).\n\nWe will have a look.","metadata":{}},{"cell_type":"code","source":"ontology = GODag(go_file)\n\n# Example: Look up a GO term\nterm = ontology[\"GO:0006915\"]\nprint('term.name:', term.name)\nprint('term.namespace:', term.namespace)\n# print(term.parents)\n# print(term.children)\n# print(term.level)\n# print(term.depth)\n# print(term.is_obsolete)\n# print(term.alt_ids)\n# print(term.is_a)\n\n# for documentation for the GoTerm class, refer to:\n# https://github.com/tanghaibao/goatools/blob/main/goatools/obo_parser.py#L157","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:23.852612Z","iopub.execute_input":"2025-12-23T21:36:23.852906Z","iopub.status.idle":"2025-12-23T21:36:25.062201Z","shell.execute_reply.started":"2025-12-23T21:36:23.852886Z","shell.execute_reply":"2025-12-23T21:36:25.061185Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# number of terms in the GO graph\nlen(ontology)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:25.064062Z","iopub.execute_input":"2025-12-23T21:36:25.064353Z","iopub.status.idle":"2025-12-23T21:36:25.071614Z","shell.execute_reply.started":"2025-12-23T21:36:25.064332Z","shell.execute_reply":"2025-12-23T21:36:25.070663Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# keep elements of GO graph that also exist in train set\nlist_of_train_terms = list(train_terms_df['term'])\nontology_filtered = {k: ontology[k] for k in list_of_train_terms}\nlen(ontology_filtered)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:25.072846Z","iopub.execute_input":"2025-12-23T21:36:25.073139Z","iopub.status.idle":"2025-12-23T21:36:25.253632Z","shell.execute_reply.started":"2025-12-23T21:36:25.073110Z","shell.execute_reply":"2025-12-23T21:36:25.252518Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# look at an example GO term and its parents\nrandom_go_term = 'GO:0043371'\nrandom_go_term_parents = ontology_filtered[random_go_term].get_all_parents()\n\nprint('GO term:')\nprint(ontology_filtered[random_go_term])\nprint('')\nprint('Parents:')\nfor parent in random_go_term_parents:\n    if parent in ontology_filtered.keys():\n        print(ontology[parent])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:25.255280Z","iopub.execute_input":"2025-12-23T21:36:25.255618Z","iopub.status.idle":"2025-12-23T21:36:25.262848Z","shell.execute_reply.started":"2025-12-23T21:36:25.255598Z","shell.execute_reply":"2025-12-23T21:36:25.261837Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"-------\n-------\n-------\n# Feature Engineering\nWe'll focus only on the train set for now.","metadata":{}},{"cell_type":"markdown","source":"## propy3\nIn this segment, we generate features from protein sequences using the `propy3` library. Extracted features include amino acid composition, dipeptide composition and CTD descriptors.\n\nYou'll find that I skip extracting many of the possible features since I am writing this notebook for demonstration purposes.","metadata":{}},{"cell_type":"code","source":"# function for generating features from a protein sequence\ndef get_protein_descriptors(sequence):\n    DesObj = PyPro.GetProDes(sequence)\n\n    all_descriptors = {}\n    all_descriptors = all_descriptors | DesObj.GetAAComp()     # Amino acid compositon descriptors (20)\n    all_descriptors = all_descriptors | DesObj.GetDPComp()     # Dipeptide composition descriptors (400)\n    all_descriptors = all_descriptors | DesObj.GetCTD()        # Composition Transition Distribution descriptors (147)\n    # all_descriptors = all_descriptors | DesObj.GetGearyAuto()  # Geary autocorrelation descriptors (240)\n    # all_descriptors = all_descriptors | DesObj.GetMoranAuto()  # Moran autocorrelation descriptors (240)\n    # all_descriptors = all_descriptors | DesObj.GetQSO()        # Quasi sequence order descriptors (default is 50)\n    # all_descriptors = all_descriptors | DesObj.GetSOCN()       # Sequence order coupling numbers (default is 45)\n\n    return all_descriptors","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:29.147461Z","iopub.execute_input":"2025-12-23T21:36:29.147834Z","iopub.status.idle":"2025-12-23T21:36:29.153393Z","shell.execute_reply.started":"2025-12-23T21:36:29.147808Z","shell.execute_reply":"2025-12-23T21:36:29.152330Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# example\nget_protein_descriptors(train_fasta_sequences[0]['seq'])","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:29.347200Z","iopub.execute_input":"2025-12-23T21:36:29.347508Z","iopub.status.idle":"2025-12-23T21:36:29.364960Z","shell.execute_reply.started":"2025-12-23T21:36:29.347486Z","shell.execute_reply":"2025-12-23T21:36:29.363479Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_descriptors = {}\nfor i in range(len(train_fasta_sequences)):\n    key_i = train_fasta_sequences[i]['id']\n    seq_i = train_fasta_sequences[i]['seq']\n    train_descriptors[key_i] = get_protein_descriptors(seq_i)\n    if i % 10000 == 0:\n        print(i, '-->', datetime.now())","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:36:41.483025Z","iopub.execute_input":"2025-12-23T21:36:41.483351Z","iopub.status.idle":"2025-12-23T21:43:59.416824Z","shell.execute_reply.started":"2025-12-23T21:36:41.483328Z","shell.execute_reply":"2025-12-23T21:43:59.415587Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# save protein descriptors\nwith open('train_descriptors.pkl', 'wb') as file:\n    pickle.dump(train_descriptors, file)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:43:59.418494Z","iopub.execute_input":"2025-12-23T21:43:59.418801Z","iopub.status.idle":"2025-12-23T21:44:22.411903Z","shell.execute_reply.started":"2025-12-23T21:43:59.418779Z","shell.execute_reply":"2025-12-23T21:44:22.411004Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# test_descriptors = {}\n# for i in range(len(test_fasta_sequences)):\n#     key_i = test_fasta_sequences[i]['id']\n#     seq_i = test_fasta_sequences[i]['seq']\n#     test_descriptors[key_i] = get_protein_descriptors(seq_i)\n#     if i % 10000 == 0:\n#         print(i, '-->', datetime.now())\n\n# test_descriptors","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:44:22.412811Z","iopub.execute_input":"2025-12-23T21:44:22.413071Z","iopub.status.idle":"2025-12-23T21:44:22.417178Z","shell.execute_reply.started":"2025-12-23T21:44:22.413051Z","shell.execute_reply":"2025-12-23T21:44:22.416057Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# # save protein descriptors\n# with open('test_descriptors.pkl', 'wb') as file:\n#     pickle.dump(test_descriptors, file)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:44:22.419014Z","iopub.execute_input":"2025-12-23T21:44:22.419249Z","iopub.status.idle":"2025-12-23T21:44:22.438800Z","shell.execute_reply.started":"2025-12-23T21:44:22.419229Z","shell.execute_reply":"2025-12-23T21:44:22.437674Z"}},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"-------\n-------\n-------\n# Modeling\nBelow, I provide only code for generating the **data** to be used in training the prediction model, not the code for training the model itself. Sadly, I am busy and don't have time to invest in trying to get the model training to work (I've had many failed attempts).","metadata":{}},{"cell_type":"markdown","source":"## Feature Matrix X and Label Matrix Y\nHere, we prepare the required components to train multi-label classification model, namely:\n- the feature matrix `X` representing the proteins in the train set\n- the label matrix `Y` indicating which of the ~26K labels are assigned to which proteins","metadata":{}},{"cell_type":"code","source":"# load pickle file (train set descriptors)\nwith open('train_descriptors.pkl', 'rb') as file:\n    train_descriptors = pickle.load(file)","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:44:22.440037Z","iopub.execute_input":"2025-12-23T21:44:22.440495Z","iopub.status.idle":"2025-12-23T21:44:34.327671Z","shell.execute_reply.started":"2025-12-23T21:44:22.440465Z","shell.execute_reply":"2025-12-23T21:44:34.326713Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# prepare feature matrix X\nX = pd.DataFrame.from_dict(train_descriptors, orient='index')\nX.index = [idx.split('|')[1] for idx in X.index]  # parse fasta ids to match train_terms_df EntryID format\nX","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:44:34.328709Z","iopub.execute_input":"2025-12-23T21:44:34.329087Z","iopub.status.idle":"2025-12-23T21:45:33.608982Z","shell.execute_reply.started":"2025-12-23T21:44:34.328939Z","shell.execute_reply":"2025-12-23T21:45:33.608097Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# sort index\nX.index = X.index.sort_values()\nX","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:45:33.609872Z","iopub.execute_input":"2025-12-23T21:45:33.610164Z","iopub.status.idle":"2025-12-23T21:45:33.720135Z","shell.execute_reply.started":"2025-12-23T21:45:33.610139Z","shell.execute_reply":"2025-12-23T21:45:33.719009Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# prepare label matrix Y\n\ndf = train_terms_df.copy()\n\n# order by EntryID\ndf = df.sort_values(by=['EntryID', 'term']).reset_index(drop=True)\n\n# Encode terms as categories for efficient mapping\ndf['EntryID'] = df['EntryID'].astype('category')\ndf['term']    = df['term'].astype('category')\n\n# Create sparse matrix (to preserve memory)\nrow = df['EntryID'].cat.codes\ncol = df['term'].cat.codes\ndata = [1] * len(df)\nY = csr_matrix((data, (row, col)),\n               shape=(df['EntryID'].nunique(),\n                      df['term'].nunique()))\n\nY","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:45:33.721094Z","iopub.execute_input":"2025-12-23T21:45:33.721435Z","iopub.status.idle":"2025-12-23T21:45:34.457548Z","shell.execute_reply.started":"2025-12-23T21:45:33.721413Z","shell.execute_reply":"2025-12-23T21:45:34.456643Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ensure that the order of EntryIDs in X matches that in Y\nprint(list(df['EntryID'].unique()) == list(X.index))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:45:34.458487Z","iopub.execute_input":"2025-12-23T21:45:34.458743Z","iopub.status.idle":"2025-12-23T21:45:34.489559Z","shell.execute_reply.started":"2025-12-23T21:45:34.458722Z","shell.execute_reply":"2025-12-23T21:45:34.488431Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# get column sums of Y (i.e., number of positive examples per label)\nmy_colsums = []\nfor i in range(Y.shape[1]):\n    col_sum = Y[:, i].toarray().reshape(-1).sum()\n    my_colsums.append(col_sum)\n\n# statistics of column sums\npd.Series(np.array(my_colsums)).describe()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2025-12-23T21:45:34.492148Z","iopub.execute_input":"2025-12-23T21:45:34.492628Z","iopub.status.idle":"2025-12-23T21:48:00.711582Z","shell.execute_reply.started":"2025-12-23T21:45:34.492604Z","shell.execute_reply":"2025-12-23T21:48:00.710651Z"}},"outputs":[],"execution_count":null}]}