{"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 5 Protein Function Prediction 😍 EDA\n\nSometimes you see a Kaggle competition and it pulls you in... have a feeling this one will be like that for me. It took me a while to fully understand what's going on, but finally completed the EDA and hopefully it's going to be a useful resource for you. \n\n### Please upvote if you find this helpful ❤️ ","metadata":{}},{"cell_type":"markdown","source":"# Proteins\n\nProteins are complex biological molecules composed of chains of amino acids. They play a vital role in living organisms and can serve a wide range of functions. Alpha-amino acids, a subset of amino acids, serve as the building blocks of proteins and can be thought of as the \"vocabulary\" of the protein language. There are around 20 alpha-amino acids in the genetic code, which are often represented by single letters (e.g. G for glycine and A for alanine).\n\nThe sequence of amino acids in a protein determines its 3D shape, which is crucial for its function. Predicting protein shape from sequence is a challenging and active area of research (see, for example, Alphafold).\n\n![](https://upload.wikimedia.org/wikipedia/commons/2/2d/Protein_CPM_PDB_1uwy.png)\n","metadata":{}},{"cell_type":"markdown","source":"# Protein Function and Gene Ontology\n\nWe are predicting protein function in this competition, so we need a set of classes that denote functions. The organizers have chosen the Gene Ontology (GO) to serve as our protein function class hierarchy. Genes encode proteins, so it makes sense that a gene ontology can also serve as a protein function hierarchy. \n\n## Gene Ontology\n\nImagine you have a big box of different LEGO pieces (genes/proteins) that you can use to build many different things. The Gene Ontology (GO) is like a guide that helps us understand how these LEGO pieces work together to build and make living things, like plants and animals, function properly.\n\nAs we want to properly understand function, we need to look at the genes/proteins from different perspectives. This is why we have the three subdomains of Gene Ontology:\n\n**Molecular Function (MF)**: Some LEGO pieces do a specific job, like a hinge that helps a door open and close. This subdomain describes the specific activities performed by a gene product at the molecular level, such as binding to a particular molecule or catalyzing a specific biochemical reaction. MF terms do not specify where or when these activities occur.\n\n**Biological Process (BP)**: These are like the steps you follow to build something with your LEGO pieces, like building a house or a car. It represents a collection of related events or pathways that a gene or gene product participates in, such as cellular metabolism, immune response, or cell division.\n\n**Cellular Component (CC)**: These are like the different parts of a LEGO model, like the walls, roof, and windows of a house. CC terms refer to cellular structures, such as organelles, protein complexes, or cell membrane, providing context for gene product function.\n\n## GO as a Graph\n\nThe information in the GO is organized using special relationships, like \"is_a,\" \"part_of,\" \"has_part,\" and \"regulates.\" These relationships connect its components in different ways, creating a web-like structure called a Directed Acyclic Graph (DAG). This structure lets us explore how genes and their functions are related, and it's more flexible than a simple hierarchy because each term can have multiple connections to broader parent terms and more specific child terms.\n\nEach gene (protein) is associated with the most specific terms that describe its functions. When a gene is associated with a term, it is also connected to all the parent terms. The structure is always evolving, with new relationships and terms being added, and obsolete ones being removed. However, even if a term is removed, its identifier remains valid but is labeled as obsolete.","metadata":{}},{"cell_type":"markdown","source":"# Dataset","metadata":{}},{"cell_type":"code","source":"import pandas as pd\nfrom pathlib import Path\nimport random\n\nrandom.seed(42)\npath = Path('../input/cafa-5-protein-function-prediction')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:09.644951Z","iopub.execute_input":"2023-06-17T11:34:09.645366Z","iopub.status.idle":"2023-06-17T11:34:09.651696Z","shell.execute_reply.started":"2023-06-17T11:34:09.645326Z","shell.execute_reply":"2023-06-17T11:34:09.650524Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Gene Ontology\nThe ontology data is in the file `go-basic.obo`. This structure is the 2023-01-01 release of the GO graph. This file is in OBO format, for which there exist many parsing libraries. For example, the obonet package is available for Python. The nodes in this graph are indexed by the term name, for example the roots of the three onotlogies are:\n\n```\nsubontology_roots = {'BPO':'GO:0008150',\n                     'CCO':'GO:0005575',\n                     'MFO':'GO:0003674'}\n```\n\n[This notebook](https://github.com/dhimmel/obonet/blob/main/examples/go-obonet.ipynb) is a good tutorial for using `obonet` to work with the ontology data. ","metadata":{}},{"cell_type":"code","source":"!pip install obonet networkx -qq","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:10.387868Z","iopub.execute_input":"2023-06-17T11:34:10.388361Z","iopub.status.idle":"2023-06-17T11:34:22.084421Z","shell.execute_reply.started":"2023-06-17T11:34:10.388315Z","shell.execute_reply":"2023-06-17T11:34:22.082903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import networkx\nimport obonet\n\ngraph = obonet.read_obo(path/'Train/go-basic.obo')\n\n# Number of nodes & edges\nlen(graph), graph.number_of_edges()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:22.087526Z","iopub.execute_input":"2023-06-17T11:34:22.087942Z","iopub.status.idle":"2023-06-17T11:34:30.667915Z","shell.execute_reply.started":"2023-06-17T11:34:22.087896Z","shell.execute_reply":"2023-06-17T11:34:30.666696Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create name mappings\nid_to_name = {id_: data.get('name') for id_, data in graph.nodes(data=True)}\nname_to_id = {data['name']: id_ for id_, data in graph.nodes(data=True) if 'name' in data}","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.669221Z","iopub.execute_input":"2023-06-17T11:34:30.669555Z","iopub.status.idle":"2023-06-17T11:34:30.731032Z","shell.execute_reply.started":"2023-06-17T11:34:30.669514Z","shell.execute_reply":"2023-06-17T11:34:30.72971Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Randomly select a  node\nrandom_node = random.choice(list(graph))\nrandom_node, id_to_name[random_node]","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.733784Z","iopub.execute_input":"2023-06-17T11:34:30.734161Z","iopub.status.idle":"2023-06-17T11:34:30.744759Z","shell.execute_reply.started":"2023-06-17T11:34:30.734125Z","shell.execute_reply":"2023-06-17T11:34:30.743139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find edges to parent terms\nfor child, parent, key in graph.out_edges(random_node, keys=True):\n    print(f'• {id_to_name[child]} ⟶ {key} ⟶ {id_to_name[parent]}')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.746439Z","iopub.execute_input":"2023-06-17T11:34:30.747566Z","iopub.status.idle":"2023-06-17T11:34:30.760531Z","shell.execute_reply.started":"2023-06-17T11:34:30.747503Z","shell.execute_reply":"2023-06-17T11:34:30.759072Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find edges to children terms\nnode = name_to_id['pilus']\nfor parent, child, key in graph.in_edges(random_node, keys=True):\n    print(f'• {id_to_name[child]} ⟵ {key} ⟵ {id_to_name[parent]}')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.761705Z","iopub.execute_input":"2023-06-17T11:34:30.761991Z","iopub.status.idle":"2023-06-17T11:34:30.777095Z","shell.execute_reply.started":"2023-06-17T11:34:30.761963Z","shell.execute_reply":"2023-06-17T11:34:30.775722Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find all superterms\nsorted(id_to_name[superterm] for superterm in networkx.descendants(graph, random_node))","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.778242Z","iopub.execute_input":"2023-06-17T11:34:30.778565Z","iopub.status.idle":"2023-06-17T11:34:30.793441Z","shell.execute_reply.started":"2023-06-17T11:34:30.778534Z","shell.execute_reply":"2023-06-17T11:34:30.791902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find all subterms\nsorted(id_to_name[subterm] for subterm in networkx.ancestors(graph, random_node))","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.795594Z","iopub.execute_input":"2023-06-17T11:34:30.796491Z","iopub.status.idle":"2023-06-17T11:34:30.807468Z","shell.execute_reply.started":"2023-06-17T11:34:30.796422Z","shell.execute_reply":"2023-06-17T11:34:30.805898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Find all paths to the root\npaths = networkx.all_simple_paths(\n    graph,\n    source=random_node,\n    target=name_to_id['biological_process']\n)\nfor pth in paths:\n    print('•', ' ⟶ '.join(id_to_name[node] for node in pth))","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.808782Z","iopub.execute_input":"2023-06-17T11:34:30.809242Z","iopub.status.idle":"2023-06-17T11:34:30.820122Z","shell.execute_reply.started":"2023-06-17T11:34:30.809198Z","shell.execute_reply":"2023-06-17T11:34:30.819095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### What does this mean? \n\n- There's a lot of potential labels (43248 nodes in the ontology graph), but not all of them may be relevant, we'll look more into this later\n- If we're sophisticated, we can leverage the DAG format of the ontologies when making predictions","metadata":{}},{"cell_type":"markdown","source":"## Train set\n\n### Training sequences\n\n`train_sequences.fasta` contains the protein sequences for the training dataset.\n\nThis files are in FASTA format, a standard format for describing protein sequences. The proteins were all retrieved from the UniProt data set curated at the European Bioinformatics Institute.\n\nThe header contains the protein's UniProt accession ID and additional information about the protein. Most protein sequences were extracted from the Swiss-Prot database, but a subset of proteins that are not represented in Swiss-Prot were extracted from the TrEMBL database. In both cases, the sequences come from the 2022_05 release from 14-Dec-2022. \n\nThe `train_sequences.fasta` file will indicate from which database the sequence originate. For example, `sp|P9WHI7|RECN_MYCT` in the FASTA header indicates the protein with UniProt ID `P9WHI7` and gene name `RECN_MYCT` was taken from Swiss-Prot (sp). Any sequences taken from TrEMBL will have tr in the header instead of sp. Swiss-Prot and TrEMBL are both parts of UniProtKB.\n\nThis file contains only sequences for proteins with annotations in the dataset (labeled proteins).","metadata":{}},{"cell_type":"code","source":"!head {path}/'Train/train_sequences.fasta'","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:30.823612Z","iopub.execute_input":"2023-06-17T11:34:30.823927Z","iopub.status.idle":"2023-06-17T11:34:31.104027Z","shell.execute_reply.started":"2023-06-17T11:34:30.823896Z","shell.execute_reply":"2023-06-17T11:34:31.10251Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see the `fasta` file contains protein identifiers, name, some metadata, and the aminoacid sequence. Let's read the protein ids below. ","metadata":{}},{"cell_type":"code","source":"def load_sequences(fasta_file):\n    sequences = {}\n    current_id = None\n    with open(fasta_file, 'rt') as file:\n        for line in file:\n            line = line.strip()\n            if line.startswith('>'):\n                # Extract the protein identifier from the line\n                identifier = line[1:].split()[0]\n                current_id = identifier\n                sequences[current_id] = ''\n            else:\n                sequences[current_id] += line\n\n    # Print number of protein sequences\n    print(str(len(sequences)) + ' ids loaded.')\n    return sequences\n\n\ntrain_file = path/\"Train/train_sequences.fasta\"\ntrain_sequences = load_sequences(train_file)","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:31.105582Z","iopub.execute_input":"2023-06-17T11:34:31.105969Z","iopub.status.idle":"2023-06-17T11:34:32.071689Z","shell.execute_reply.started":"2023-06-17T11:34:31.105926Z","shell.execute_reply":"2023-06-17T11:34:32.069621Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's check how many aminoacids are represented in train sequences.","metadata":{}},{"cell_type":"code","source":"train_alphabet = list(set(''.join(train_sequences.values())))\nprint(f'Number of aminoacids: {len(train_alphabet)}, aminoacids: {\"\".join(sorted(train_alphabet))}')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:16:16.826587Z","iopub.execute_input":"2023-06-17T13:16:16.82699Z","iopub.status.idle":"2023-06-17T13:16:18.10711Z","shell.execute_reply.started":"2023-06-17T13:16:16.826949Z","shell.execute_reply":"2023-06-17T13:16:18.104911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It's important to see how long the sequences are, as certain models have a max cap on the sequence length. We can see both very short (3 residues) and very long (35375 residues sequences) that might require some special treatment. ","metadata":{}},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nsequence_lengths = [len(sequence) for sequence in train_sequences.values()]\n\nprint(f'Max: {max(sequence_lengths)}, min: {min(sequence_lengths)}.')\n\n# Plot histogram\nplt.hist([x for x in sequence_lengths if x < 6000], bins=50)\nplt.xlabel('Sequence Length')\nplt.ylabel('Count')\nplt.title('Histogram of Sequence Lengths (capped at 6.000)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:33.3451Z","iopub.execute_input":"2023-06-17T11:34:33.345432Z","iopub.status.idle":"2023-06-17T11:34:34.081507Z","shell.execute_reply.started":"2023-06-17T11:34:33.345402Z","shell.execute_reply":"2023-06-17T11:34:34.079857Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Taxonomy\n\n`train_taxonomy.tsv` contains the list of proteins and the species to which they belong, represented by a \"taxonomic identifier\" (taxon ID) number. The first column is the protein UniProt accession ID and the second is the taxon ID. ","metadata":{}},{"cell_type":"code","source":"taxonomy = pd.read_csv(path/'Train/train_taxonomy.tsv', sep='\\t')\ntaxonomy.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:34.0842Z","iopub.execute_input":"2023-06-17T11:34:34.085259Z","iopub.status.idle":"2023-06-17T11:34:34.214913Z","shell.execute_reply.started":"2023-06-17T11:34:34.085201Z","shell.execute_reply":"2023-06-17T11:34:34.21342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Group the DataFrame by Taxon ID and count the number of occurrences\ntrain_counts = taxonomy.groupby(\"taxonomyID\").size().reset_index(name=\"Count\")\n\n# Sort the counts in descending order\ntrain_counts = train_counts.sort_values(\"Count\", ascending=False).reset_index(drop=True)\n\ntrain_counts.head(100).Count.plot(kind='line');","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:27:18.186615Z","iopub.execute_input":"2023-06-17T13:27:18.187038Z","iopub.status.idle":"2023-06-17T13:27:18.351051Z","shell.execute_reply.started":"2023-06-17T13:27:18.187002Z","shell.execute_reply":"2023-06-17T13:27:18.349544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see a large number of proteins coming from few species, and a long tail of species with very few proteins. ","metadata":{}},{"cell_type":"markdown","source":"### Labels\n\n`train_terms.tsv` contains the list of annotated terms (ground truth) for the proteins in `train_sequences.fasta`. The first column indicates the protein's UniProt accession ID, the second is the GO term ID, and the third indicates in which ontology the term appears.\n","metadata":{}},{"cell_type":"code","source":"!head {path}/'Train/train_terms.tsv'","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:35:59.320801Z","iopub.execute_input":"2023-06-17T11:35:59.322132Z","iopub.status.idle":"2023-06-17T11:35:59.595853Z","shell.execute_reply.started":"2023-06-17T11:35:59.322061Z","shell.execute_reply":"2023-06-17T11:35:59.594086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_labels = pd.read_csv(path/'Train/train_terms.tsv', sep='\\t')\ntrain_labels.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:34.444715Z","iopub.execute_input":"2023-06-17T11:34:34.445016Z","iopub.status.idle":"2023-06-17T11:34:37.126888Z","shell.execute_reply.started":"2023-06-17T11:34:34.444986Z","shell.execute_reply":"2023-06-17T11:34:37.125547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's explore this in a bit more detail. ","metadata":{}},{"cell_type":"code","source":"len(train_labels), train_labels.term.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:34:37.128192Z","iopub.execute_input":"2023-06-17T11:34:37.128529Z","iopub.status.idle":"2023-06-17T11:34:37.549446Z","shell.execute_reply.started":"2023-06-17T11:34:37.128498Z","shell.execute_reply":"2023-06-17T11:34:37.547881Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate min, max, median, and average number of GO terms per protein id\ngo_terms_per_protein = train_labels.groupby('EntryID')['term'].count()\nmin_go_terms_per_protein = go_terms_per_protein.min()\nmax_go_terms_per_protein = go_terms_per_protein.max()\nmedian_go_terms_per_protein = go_terms_per_protein.median()\naverage_go_terms_per_protein = go_terms_per_protein.mean()\n\n# Print the results\nprint(\"Min GO terms per protein ID:\", min_go_terms_per_protein)\nprint(\"Max GO terms per protein ID:\", max_go_terms_per_protein)\nprint(\"Median GO terms per protein ID:\", median_go_terms_per_protein)\nprint(\"Average GO terms per protein ID:\", average_go_terms_per_protein)","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:55:36.758116Z","iopub.execute_input":"2023-06-17T11:55:36.759072Z","iopub.status.idle":"2023-06-17T11:55:37.498775Z","shell.execute_reply.started":"2023-06-17T11:55:36.75903Z","shell.execute_reply":"2023-06-17T11:55:37.497258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate min, max, median, and average number of unique protein id per GO term\nproteins_per_go_term = train_labels.groupby('term')['EntryID'].nunique()\nmin_proteins_per_go_term = proteins_per_go_term.min()\nmax_proteins_per_go_term = proteins_per_go_term.max()\nmedian_proteins_per_go_term = proteins_per_go_term.median()\naverage_proteins_per_go_term = proteins_per_go_term.mean()\n\nprint(\"Min unique protein ID per GO term:\", min_proteins_per_go_term)\nprint(\"Max unique protein ID per GO term:\", max_proteins_per_go_term)\nprint(\"Median unique protein ID per GO term:\", median_proteins_per_go_term)\nprint(\"Average unique protein ID per GO term:\", average_proteins_per_go_term)","metadata":{"execution":{"iopub.status.busy":"2023-06-17T11:57:48.488837Z","iopub.execute_input":"2023-06-17T11:57:48.489252Z","iopub.status.idle":"2023-06-17T11:57:50.078115Z","shell.execute_reply.started":"2023-06-17T11:57:48.489217Z","shell.execute_reply":"2023-06-17T11:57:50.076185Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate min, max, median, and average number of unique protein ID per GO term for each aspect\nunique_proteins_per_go_term = train_labels.groupby(['aspect', 'term'])['EntryID'].nunique()\nmin_proteins_per_go_term = unique_proteins_per_go_term.groupby('aspect').min()\nmax_proteins_per_go_term = unique_proteins_per_go_term.groupby('aspect').max()\nmedian_proteins_per_go_term = unique_proteins_per_go_term.groupby('aspect').median()\naverage_proteins_per_go_term = unique_proteins_per_go_term.groupby('aspect').mean()\n\nprint(\"Min unique protein ID per GO term per aspect:\")\nfor aspect, value in min_proteins_per_go_term.items():\n    print(aspect, value)\nprint()\nprint(\"Max unique protein ID per GO term per aspect:\")\nfor aspect, value in max_proteins_per_go_term.items():\n    print(aspect, value)\nprint()\nprint(\"Median unique protein ID per GO term per aspect:\")\nfor aspect, value in median_proteins_per_go_term.items():\n    print(aspect, value)\nprint()\nprint(\"Average unique protein ID per GO term per aspect:\")\nfor aspect, value in average_proteins_per_go_term.items():\n    print(aspect, value)\n","metadata":{"execution":{"iopub.status.busy":"2023-06-17T12:03:37.988582Z","iopub.execute_input":"2023-06-17T12:03:37.989601Z","iopub.status.idle":"2023-06-17T12:03:39.976709Z","shell.execute_reply.started":"2023-06-17T12:03:37.989554Z","shell.execute_reply":"2023-06-17T12:03:39.974927Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Count the number of unique terms per aspect (ontology)\nunique_terms_per_aspect = train_labels.groupby('aspect')['term'].nunique()\n\nprint(\"Unique terms per\", unique_terms_per_aspect)","metadata":{"execution":{"iopub.status.busy":"2023-06-17T12:06:08.155581Z","iopub.execute_input":"2023-06-17T12:06:08.155964Z","iopub.status.idle":"2023-06-17T12:06:09.561983Z","shell.execute_reply.started":"2023-06-17T12:06:08.155928Z","shell.execute_reply":"2023-06-17T12:06:09.56039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create histograms of GO terms per aspect\nhistograms = train_labels.groupby('aspect')['term'].value_counts()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Convert the output into a readable chart\nchart_data = []\nfor aspect, histogram in histograms.groupby(level=0):\n    aspect_chart_data = []\n    for nterms in [10, 50, 100, 200, 300, 400, 500]:\n        min_freq = min(histogram.head(nterms))\n        aspect_chart_data.append((nterms, min_freq))\n    chart_data.append((aspect, aspect_chart_data))\n    print(f'{aspect}: top term has {max(histogram)} annotations, 500th top term has {min(histogram.head(500))} annotations.')\n\nprint()\n    \n# Plot the chart\nfig, ax = plt.subplots()\n\n# Iterate over each aspect's data\nfor aspect, data in chart_data:\n    x_values, y_values = zip(*data)\n    ax.plot(x_values, y_values, marker='o', label=aspect)\n\n# Set chart title and labels\nax.set_title(\"Minimum Frequency of Top Terms in Each Aspect at Different Cutoffs\")\nax.set_xlabel(\"Number of Top Terms\")\nax.set_ylabel(\"Minimum Frequency\")\n\n# Add legend\nax.legend()\n\n# Show the chart\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Calculate the number of terms needed for a minimum of 50 annotations per term\nnum_terms_needed = []\nfor aspect, histogram in histograms.groupby(level=0):\n    terms_needed = 0\n    for term, freq in histogram.iteritems():\n        terms_needed += 1\n        if freq < 50:\n            num_terms_needed.append({aspect: terms_needed})\n            break\n            \nprint('This is how many terms we can keep if we set threshold at minimum 50 annotations per term:')\nfor x in num_terms_needed:\n    print(x)","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:01:17.353188Z","iopub.execute_input":"2023-06-17T13:01:17.353787Z","iopub.status.idle":"2023-06-17T13:01:17.379374Z","shell.execute_reply.started":"2023-06-17T13:01:17.353747Z","shell.execute_reply":"2023-06-17T13:01:17.378079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Test\n### Test sequences\n\n`testsuperset.fasta` contains protein sequences on which the participants are asked to submit predictions. The header for each sequence in `testsuperset.fasta` contains the protein's UniProt accession ID and the Taxon ID of the species this protein belongs to. Only a small subset of those sequences will accumulate functional annotations and will constitute the test set. ","metadata":{}},{"cell_type":"code","source":"!head {path}/'Test (Targets)/testsuperset.fasta'","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:01:44.887231Z","iopub.execute_input":"2023-06-17T13:01:44.888423Z","iopub.status.idle":"2023-06-17T13:01:45.252392Z","shell.execute_reply.started":"2023-06-17T13:01:44.888372Z","shell.execute_reply":"2023-06-17T13:01:45.251037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's read this and check similar properties to train set and how it might be different. This time we will use the `BioPython` package. ","metadata":{}},{"cell_type":"code","source":"from Bio import SeqIO\nimport pandas as pd\n\ntest_file = path/'Test (Targets)/testsuperset.fasta'\n\ndata = []\nfor record in SeqIO.parse(test_file, \"fasta\"):\n    protein_id, taxon = record.description.split(\"\\t\")\n    sequence = str(record.seq)\n    data.append((protein_id, sequence, taxon))\n\ntest_sequences = pd.DataFrame(data, columns=[\"Protein ID\", \"Sequence\", \"Taxon\"])\nprint(f'Loaded {len(test_sequences)} test sequences.\\n')\ntest_sequences.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:12:01.515813Z","iopub.execute_input":"2023-06-17T13:12:01.516248Z","iopub.status.idle":"2023-06-17T13:12:03.326309Z","shell.execute_reply.started":"2023-06-17T13:12:01.516208Z","shell.execute_reply":"2023-06-17T13:12:03.324893Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len([x for x in test_sequences[\"Protein ID\"].values if x in train_sequences.keys()])","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:21:15.80965Z","iopub.execute_input":"2023-06-17T13:21:15.810024Z","iopub.status.idle":"2023-06-17T13:21:15.880383Z","shell.execute_reply.started":"2023-06-17T13:21:15.809993Z","shell.execute_reply":"2023-06-17T13:21:15.879687Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A bunch of proteins in the test superset overlaps with train - this is ok, because we're only evaluated on new experimentally verified function annotations. ","metadata":{}},{"cell_type":"code","source":"test_alphabet = list(set(''.join(test_sequences.Sequence.values)))\nprint(f'Number of aminoacids: {len(test_alphabet)}, aminoacids: {\"\".join(sorted(test_alphabet))}')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:24:14.056134Z","iopub.execute_input":"2023-06-17T13:24:14.056513Z","iopub.status.idle":"2023-06-17T13:24:15.128763Z","shell.execute_reply.started":"2023-06-17T13:24:14.056481Z","shell.execute_reply":"2023-06-17T13:24:15.127592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like test set doesn't include any sequences with a rare aminoacid 'O' (pyrrolysine).","metadata":{}},{"cell_type":"code","source":"sequence_lengths = [len(sequence) for sequence in test_sequences.Sequence.values]\n\nprint(f'Max: {max(sequence_lengths)}, min: {min(sequence_lengths)}.')\n\n# Plot histogram\nplt.hist([x for x in sequence_lengths if x < 6000], bins=50)\nplt.xlabel('Sequence Length')\nplt.ylabel('Count')\nplt.title('Histogram of Test Sequence Lengths (capped at 6.000)')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:24:56.170591Z","iopub.execute_input":"2023-06-17T13:24:56.170955Z","iopub.status.idle":"2023-06-17T13:24:56.894968Z","shell.execute_reply.started":"2023-06-17T13:24:56.170926Z","shell.execute_reply":"2023-06-17T13:24:56.894241Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Group the DataFrame by Taxon ID and count the number of occurrences\ntest_counts = test_sequences.groupby(\"Taxon\").size().reset_index(name=\"Count\")\n\n# Sort the counts in descending order\ntest_counts = test_counts.sort_values(\"Count\", ascending=False).reset_index(drop=True)\n\ntest_counts.head(100).Count.plot(kind='line');","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:27:44.326324Z","iopub.execute_input":"2023-06-17T13:27:44.326926Z","iopub.status.idle":"2023-06-17T13:27:44.522727Z","shell.execute_reply.started":"2023-06-17T13:27:44.326871Z","shell.execute_reply":"2023-06-17T13:27:44.522055Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's compare the most popular taxons in train and test: ","metadata":{}},{"cell_type":"code","source":"test_counts.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:30:08.263735Z","iopub.execute_input":"2023-06-17T13:30:08.264117Z","iopub.status.idle":"2023-06-17T13:30:08.278428Z","shell.execute_reply.started":"2023-06-17T13:30:08.264082Z","shell.execute_reply":"2023-06-17T13:30:08.277346Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_counts['Taxon'] = train_counts.taxonomyID.astype('str')\ncompare_counts = pd.merge(test_counts, train_counts[['Taxon', 'Count']], how='left', on='Taxon', suffixes=['_test', '_train'])","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:33:20.462032Z","iopub.execute_input":"2023-06-17T13:33:20.46246Z","iopub.status.idle":"2023-06-17T13:33:20.475938Z","shell.execute_reply.started":"2023-06-17T13:33:20.462426Z","shell.execute_reply":"2023-06-17T13:33:20.474705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"compare_counts['Diff'] = compare_counts['Count_test'] - compare_counts['Count_train']\ncompare_counts = compare_counts.sort_values('Diff', ascending=False)\ncompare_counts.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:35:37.173866Z","iopub.execute_input":"2023-06-17T13:35:37.174293Z","iopub.status.idle":"2023-06-17T13:35:37.189781Z","shell.execute_reply.started":"2023-06-17T13:35:37.174255Z","shell.execute_reply":"2023-06-17T13:35:37.18831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks like there are some species that have much more representation in test superset than train. ","metadata":{}},{"cell_type":"markdown","source":"The file `testsuperset-taxon-list.tsv` is a set of taxon IDs for the proteins in the test superset.","metadata":{}},{"cell_type":"code","source":"!head {path}/'Test (Targets)/testsuperset-taxon-list.tsv'","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:01:46.377259Z","iopub.execute_input":"2023-06-17T13:01:46.37765Z","iopub.status.idle":"2023-06-17T13:01:46.664277Z","shell.execute_reply.started":"2023-06-17T13:01:46.377615Z","shell.execute_reply":"2023-06-17T13:01:46.663322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Submission File\n\nThe list of predictions contains a list of pairs between **protein targets** and **GO terms**, followed by the **probabilistic estimate of the relationship** (one association per line). The target name must correspond to the target ID listed in the test set (in the FASTA header for each sequence). The GO ID must correspond to valid terms in GO's version listed in the Data section---invalid terms are automatically excluded from evaluation. Molecular Function (MF), Biological Process (BP), and Cellular Component (CC) subontologies of GO are to be combined in the prediction files, but they will be evaluated independently and combined at the end as described above. The score must be in the interval (0, 1.000] and contain up to 3 (three) significant figures. A score of 0 is not allowed; that is, the team should simply not list such pairs. In case the predictions in the submitted files are not propagated to the root of ontology, the predictions will be recursively propagated by assigning each parent term a score that is the maximum score among its children's scores. Finally, to limit prediction file sizes, one target cannot be associated with more than 1500 terms for MF, BP, and CC subontologies combined.\n\nFor any protein ID in the test superset, you must list a set of GO terms and assign your estimated probability. If a protein ID is not listed in your submitted file, the organizers will assume that all predictions are 0. The file should not contain a header, columns must be tab or space separated. An example submission file may look as follows:","metadata":{}},{"cell_type":"code","source":"!head {path}/'sample_submission.tsv'","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:37:16.792154Z","iopub.execute_input":"2023-06-17T13:37:16.792579Z","iopub.status.idle":"2023-06-17T13:37:17.093519Z","shell.execute_reply.started":"2023-06-17T13:37:16.792539Z","shell.execute_reply":"2023-06-17T13:37:17.091889Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Naive submission\n\nLet's try the following for the naive submission - we'll take 1500 most frequent GO terms and assign them with their relative frequencies to each test sequence. \n\n","metadata":{}},{"cell_type":"code","source":"from collections import Counter\n\ndef assign_frequencies(train_df, test_sequences, num_terms):\n    # Extract the frequencies of GO terms from the train dataframe\n    term_frequencies = Counter(train_df['term'])\n\n    # Select the top 'num_terms' frequent GO terms\n    most_common_terms = [term for term, _ in term_frequencies.most_common(num_terms)]\n\n    # Assign frequencies to test sequences\n    test_assignments = []\n    for sequence in test_sequences:\n        term_counts = Counter(sequence)\n        frequencies = [term_counts[term] for term in most_common_terms]\n        test_assignments.append(frequencies)\n\n    return test_assignments\n\nnum_terms = 1500\n\nassignments = assign_frequencies(train_labels, test_sequences.Sequence.values, num_terms)","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:53:37.111347Z","iopub.execute_input":"2023-06-17T13:53:37.111679Z","iopub.status.idle":"2023-06-17T13:54:37.530752Z","shell.execute_reply.started":"2023-06-17T13:53:37.111648Z","shell.execute_reply":"2023-06-17T13:54:37.529697Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term_frequencies = Counter(train_labels['term'])","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:55:26.500786Z","iopub.execute_input":"2023-06-17T13:55:26.501196Z","iopub.status.idle":"2023-06-17T13:55:27.485857Z","shell.execute_reply.started":"2023-06-17T13:55:26.501142Z","shell.execute_reply":"2023-06-17T13:55:27.484438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term_frequencies = pd.DataFrame(data=term_frequencies.most_common(num_terms), columns=['term', 'confidence'])","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:57:20.215002Z","iopub.execute_input":"2023-06-17T13:57:20.215394Z","iopub.status.idle":"2023-06-17T13:57:20.232478Z","shell.execute_reply.started":"2023-06-17T13:57:20.215361Z","shell.execute_reply":"2023-06-17T13:57:20.230931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term_frequencies.confidence = term_frequencies.confidence / len(train_sequences)\nterm_frequencies.confidence.min(), term_frequencies.confidence.max()","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:58:18.865542Z","iopub.execute_input":"2023-06-17T13:58:18.86598Z","iopub.status.idle":"2023-06-17T13:58:18.877105Z","shell.execute_reply.started":"2023-06-17T13:58:18.865934Z","shell.execute_reply":"2023-06-17T13:58:18.875795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term_frequencies['ids'] = [test_sequences['Protein ID'].values.tolist() for _ in range(len(term_frequencies))]\n                        ","metadata":{"execution":{"iopub.status.busy":"2023-06-17T14:04:02.812969Z","iopub.execute_input":"2023-06-17T14:04:02.81332Z","iopub.status.idle":"2023-06-17T14:04:11.613543Z","shell.execute_reply.started":"2023-06-17T14:04:02.813289Z","shell.execute_reply":"2023-06-17T14:04:11.611937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"naive_sub = term_frequencies.explode('ids')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T14:04:32.199427Z","iopub.execute_input":"2023-06-17T14:04:32.199788Z","iopub.status.idle":"2023-06-17T14:05:26.795033Z","shell.execute_reply.started":"2023-06-17T14:04:32.199758Z","shell.execute_reply":"2023-06-17T14:05:26.793853Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"naive_sub.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"naive_sub[['ids', 'term', 'confidence']].to_csv('submission.tsv', header=None, index=False, sep='\\t')","metadata":{"execution":{"iopub.status.busy":"2023-06-17T14:06:15.387695Z","iopub.execute_input":"2023-06-17T14:06:15.388065Z","iopub.status.idle":"2023-06-17T14:15:03.643449Z","shell.execute_reply.started":"2023-06-17T14:06:15.388035Z","shell.execute_reply":"2023-06-17T14:15:03.641473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Evaluation\n\nThis competition focuses on predicting Gene Ontology (GO) terms for protein sequences in three different categories: Molecular Function (MF), Biological Process (BP), and Cellular Component (CC). Proteins in the test set, which were initially lacking experimentally determined functions, accumulate these annotations between the submission and evaluation phase. These proteins form three different test sets for evaluation, each corresponding to one of the GO categories. The evaluation metric is the maximum F-measure, a weighted average of precision and recall, calculated for each test set. The final performance is the arithmetic mean of the maximum F-measure of the three categories. Weights are based on the term's occurrence frequency in a large pool of proteins. The root terms have a weight of 0 as they appear in every protein's annotation, while terms deeper in the ontology are assigned larger weights as they appear less frequently and are harder to predict. The submission file should contain pairs of protein targets and GO terms along with their associated probability estimates. The evaluation is carried out on the combined set of no-knowledge and limited-knowledge protein targets.\n\n\n### Information accretion\n\n`IA.txt` contains the information accretion (weights) for each GO term. These weights are used to compute weighted precision and recall, as described in the [Evaluation section](https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/overview/evaluation). Here is a [discussion post by host on IA](https://www.kaggle.com/competitions/cafa-5-protein-function-prediction/discussion/405237).\n\nInformation Accretion (IA) is a concept used in bioinformatics to quantify the amount of novel information a specific node or term brings into an ontology, assuming its parent terms are already annotated. In simpler terms, it's a way to assess how much \"new knowledge\" a term brings to the table if we already know about its parent terms. When dealing with a large ontology like protein functions, each term represents a certain aspect of a protein's function, and the 'parent' terms are more general functions that encompass the more specific 'child' terms. IA thus measures how much unique information is added by knowing about the specific child term, on top of what we already know from the parent terms.\n\nTo calculate IA, we use the observed annotations of terms in a dataset. Specifically, the probability of a term is estimated based on how often it appears in the dataset, as is the probability of its parent terms. The IA for a given term is then computed as the log ratio of the probability of the parent terms to the probability of the term itself. A key point to remember is that if a term hasn't appeared in the training set, we add an artificial count of 1 to prevent a zero probability. Also, the IA of a term could be zero if every protein annotated with the parent terms also has the term annotated, implying no additional information is brought by this term. Lastly, it's worth noting that terms deeper in the ontology usually have a lower IA, as they are more specific and thus less frequently annotated in proteins.","metadata":{}},{"cell_type":"code","source":"!head {path}/'IA.txt'","metadata":{"execution":{"iopub.status.busy":"2023-06-17T13:37:20.406106Z","iopub.execute_input":"2023-06-17T13:37:20.407057Z","iopub.status.idle":"2023-06-17T13:37:20.69454Z","shell.execute_reply.started":"2023-06-17T13:37:20.407013Z","shell.execute_reply":"2023-06-17T13:37:20.692662Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Disclaimer\n\nI'm not an expert in this domain and I'm working with GPT copilot. The notebook might contain some factual mistakes or errors, but this process helps me learn so please take it as-is. ","metadata":{}},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}