{"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":"# EXPLORATORY DATA ANALYSIS","metadata":{}},{"cell_type":"code","source":"!pip install obonet -q\n!pip install pyvis -q","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2023-07-07T20:12:08.919813Z","iopub.execute_input":"2023-07-07T20:12:08.920698Z","iopub.status.idle":"2023-07-07T20:12:38.903555Z","shell.execute_reply.started":"2023-07-07T20:12:08.920643Z","shell.execute_reply":"2023-07-07T20:12:38.901659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import os\nimport json\nfrom PIL import Image\nfrom typing import Dict\nfrom collections import Counter\n\nimport random\nimport cv2\nimport obonet\nimport networkx\nimport pandas as pd\nimport numpy as np\nimport plotly.express as px\nimport plotly.graph_objects as go\nimport matplotlib.pyplot as plt\nimport matplotlib.patches as mpatch\nfrom Bio import SeqIO\nfrom pyvis.network import Network","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:12:38.908905Z","iopub.execute_input":"2023-07-07T20:12:38.90952Z","iopub.status.idle":"2023-07-07T20:12:40.794359Z","shell.execute_reply.started":"2023-07-07T20:12:38.909458Z","shell.execute_reply":"2023-07-07T20:12:40.793067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### The CAFA dataset contains the following important files:\n- **Go-basic.obo:** GO graph data. Each node of the graph contains info on GO terms and relationships with other GO terms.\n- **Train_sequences.fasta:** The list of proteins with unique ids, some meta info and sequence.\n- **Train_taxonomy.tsv:** It contains the taxonomy ID of proteins\n- **Train_term.tsv:** Contains mapping of the protein ids with the GO terms ids.\n- **IA.txt:** Information Accretion for each term. This is used to weight precision and recall","metadata":{}},{"cell_type":"code","source":"class CFG:\n    train_go_obo_path: str = \"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\"\n    train_seq_fasta_path: str = \"/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta\"\n    train_terms_path: str = \"/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv\"\n    train_taxonomy_path: str = \"/kaggle/input/cafa-5-protein-function-prediction/Train/train_taxonomy.tsv\"\n    train_ia_path: str = \"/kaggle/input/cafa-5-protein-function-prediction/IA.txt\"","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:12:40.797421Z","iopub.execute_input":"2023-07-07T20:12:40.797928Z","iopub.status.idle":"2023-07-07T20:12:40.80511Z","shell.execute_reply.started":"2023-07-07T20:12:40.797879Z","shell.execute_reply":"2023-07-07T20:12:40.802964Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_dag(graph, term, radius=1):\n    # create smaller subgraph\n    # radius - include all neighbors of distance<=radius from n (increse it to add further parent's branches).\n    ng_graph = networkx.ego_graph(graph, term, radius=radius)\n\n    for n in ng_graph.nodes(data=True):\n        # concatenate label of the node with its attribute\n        n[1][\"label\"] = n[0] + \" \" +n[1][\"name\"]\n\n    nt = Network(directed=True, notebook=True, cdn_resources=\"in_line\")\n    nt.from_nx(ng_graph)\n    return nt.show(\"network.html\")","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:12:40.807904Z","iopub.execute_input":"2023-07-07T20:12:40.808285Z","iopub.status.idle":"2023-07-07T20:12:40.821789Z","shell.execute_reply.started":"2023-07-07T20:12:40.808252Z","shell.execute_reply":"2023-07-07T20:12:40.820631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### load the GO terms graph using the obonet tool","metadata":{}},{"cell_type":"code","source":"graph = obonet.read_obo(CFG.train_go_obo_path)","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:12:40.824586Z","iopub.execute_input":"2023-07-07T20:12:40.825728Z","iopub.status.idle":"2023-07-07T20:13:00.669271Z","shell.execute_reply.started":"2023-07-07T20:12:40.825675Z","shell.execute_reply":"2023-07-07T20:13:00.667868Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Number of nodes: {len(graph)}\")","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:13:00.67112Z","iopub.execute_input":"2023-07-07T20:13:00.671606Z","iopub.status.idle":"2023-07-07T20:13:00.678239Z","shell.execute_reply.started":"2023-07-07T20:13:00.671561Z","shell.execute_reply":"2023-07-07T20:13:00.6769Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Number of edges: {graph.number_of_edges()}\")","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:13:00.685036Z","iopub.execute_input":"2023-07-07T20:13:00.685779Z","iopub.status.idle":"2023-07-07T20:13:00.84253Z","shell.execute_reply.started":"2023-07-07T20:13:00.685743Z","shell.execute_reply":"2023-07-07T20:13:00.841306Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each graph node contains info on a GO term, its aspect, short description and relations with other nodes","metadata":{}},{"cell_type":"code","source":"graph.nodes[\"GO:0034655\"]","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:32:19.810761Z","iopub.execute_input":"2023-07-07T20:32:19.811287Z","iopub.status.idle":"2023-07-07T20:32:19.819896Z","shell.execute_reply.started":"2023-07-07T20:32:19.81124Z","shell.execute_reply":"2023-07-07T20:32:19.818544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The immediate neighbours of the GO term:","metadata":{}},{"cell_type":"code","source":"plot_dag(graph, \"GO:0034655\", radius=1)","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:33:57.567774Z","iopub.execute_input":"2023-07-07T20:33:57.56827Z","iopub.status.idle":"2023-07-07T20:33:57.884958Z","shell.execute_reply.started":"2023-07-07T20:33:57.568222Z","shell.execute_reply":"2023-07-07T20:33:57.883753Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now lets look at larger amount of neighbours","metadata":{}},{"cell_type":"code","source":"plot_dag(graph, \"GO:0034655\", radius=1000)","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:35:10.404906Z","iopub.execute_input":"2023-07-07T20:35:10.405333Z","iopub.status.idle":"2023-07-07T20:35:10.693995Z","shell.execute_reply.started":"2023-07-07T20:35:10.405298Z","shell.execute_reply":"2023-07-07T20:35:10.692813Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Next load and analyze the train protein sequences.\n#### This file is 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.","metadata":{}},{"cell_type":"code","source":"sequences = SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\")\nnum_sequences = sum(1 for seq in sequences)\nprint(num_sequences)","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:13:36.697281Z","iopub.execute_input":"2023-07-07T20:13:36.697679Z","iopub.status.idle":"2023-07-07T20:13:39.271896Z","shell.execute_reply.started":"2023-07-07T20:13:36.697648Z","shell.execute_reply":"2023-07-07T20:13:39.270592Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The following code provides distribution of lengths of the protein sequences. Most of the proteins have length of 100 to 1000 amino acids. The peak occours at 300","metadata":{}},{"cell_type":"code","source":"sequences = SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\")\n\n# get the length of each sequence\nlengths = [len(seq) for seq in sequences]\n\nfig = px.histogram(x=lengths, nbins=1000, color_discrete_sequence=['blue'])\nfig.update_layout(\n    title={\n        'text': \"Distribution of protein sequence lengths\",\n        'y':0.95,\n        'x':0.5,\n        'xanchor': 'center',\n        'yanchor': 'top'\n    },\n    xaxis_title=\"Sequence length\", yaxis_title=\"Count\"\n)\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:17:12.429834Z","iopub.execute_input":"2023-07-07T20:17:12.430358Z","iopub.status.idle":"2023-07-07T20:17:14.87331Z","shell.execute_reply.started":"2023-07-07T20:17:12.430314Z","shell.execute_reply":"2023-07-07T20:17:14.872136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### The composition of amino acids among all the proteins is given by the following code","metadata":{}},{"cell_type":"code","source":"records = SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\")\n\naa_list = [aa for record in records for aa in record.seq]\naa_count = Counter(aa_list)\n\nfig = px.bar(\n    x=list(aa_count.values()), y=list(aa_count.keys()),\n    color_discrete_sequence=['goldenrod'],\n    orientation='h', height=700\n)\nfig.update_layout(\n    title={\n        'y':0.95,\n        'x':0.5,\n        'xanchor': 'center',\n        'yanchor': 'top'\n    },\n    xaxis_title=\"frequencies\", yaxis_title=\"Amino Acids\"\n)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:25:02.568271Z","iopub.execute_input":"2023-07-07T20:25:02.568716Z","iopub.status.idle":"2023-07-07T20:25:13.317521Z","shell.execute_reply.started":"2023-07-07T20:25:02.568682Z","shell.execute_reply":"2023-07-07T20:25:13.316372Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### load the train terms file","metadata":{}},{"cell_type":"code","source":"train_terms_df = pd.read_csv(CFG.train_terms_path, sep=\"\\t\")\ntrain_terms_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:25:57.237728Z","iopub.execute_input":"2023-07-07T20:25:57.23816Z","iopub.status.idle":"2023-07-07T20:26:01.284555Z","shell.execute_reply.started":"2023-07-07T20:25:57.238123Z","shell.execute_reply":"2023-07-07T20:26:01.283327Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms_df.describe()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:26:07.553587Z","iopub.execute_input":"2023-07-07T20:26:07.553994Z","iopub.status.idle":"2023-07-07T20:26:12.161752Z","shell.execute_reply.started":"2023-07-07T20:26:07.553962Z","shell.execute_reply":"2023-07-07T20:26:12.160664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### most of the go terms belong to the BPO class:","metadata":{}},{"cell_type":"code","source":"aspect_counts = train_terms_df.aspect.value_counts()\n\nfig = px.pie(values=aspect_counts.values, names=aspect_counts.index)\nfig.update_traces(textposition='inside', textfont_size=14)\nfig.update_layout(\n    title={\n        'text': \"Pie distribution of aspect values\",\n        'y':0.95,\n        'x':0.5,\n        'xanchor': 'center',\n        'yanchor': 'top'\n    },\n    legend_title_text='Aspect:'\n)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:26:12.164258Z","iopub.execute_input":"2023-07-07T20:26:12.164773Z","iopub.status.idle":"2023-07-07T20:26:13.188798Z","shell.execute_reply.started":"2023-07-07T20:26:12.164724Z","shell.execute_reply":"2023-07-07T20:26:13.187599Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_taxonomy_df = pd.read_csv(CFG.train_taxonomy_path, sep=\"\\t\")\ntrain_taxonomy_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:30:59.444292Z","iopub.execute_input":"2023-07-07T20:30:59.444797Z","iopub.status.idle":"2023-07-07T20:30:59.577055Z","shell.execute_reply.started":"2023-07-07T20:30:59.444759Z","shell.execute_reply":"2023-07-07T20:30:59.575975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_taxonomy_df.describe()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:30:59.86389Z","iopub.execute_input":"2023-07-07T20:30:59.864563Z","iopub.status.idle":"2023-07-07T20:30:59.891224Z","shell.execute_reply.started":"2023-07-07T20:30:59.864527Z","shell.execute_reply":"2023-07-07T20:30:59.890301Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"len(train_taxonomy_df)\n# matches with the num of unique enteries in train_terms.fasta","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:31:00.279452Z","iopub.execute_input":"2023-07-07T20:31:00.280042Z","iopub.status.idle":"2023-07-07T20:31:00.286725Z","shell.execute_reply.started":"2023-07-07T20:31:00.280007Z","shell.execute_reply":"2023-07-07T20:31:00.285798Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"merged_df = pd.merge(train_terms_df,train_taxonomy_df,on='EntryID')\nmerged_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:31:03.016625Z","iopub.execute_input":"2023-07-07T20:31:03.017042Z","iopub.status.idle":"2023-07-07T20:31:05.010659Z","shell.execute_reply.started":"2023-07-07T20:31:03.01701Z","shell.execute_reply":"2023-07-07T20:31:05.009284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"limit = 10\n\nwith open(CFG.train_ia_path) as f:\n    ia_weights = [x.replace(\"\\n\", \"\").split(\"\\t\") for x in f.readlines()]\n\nia_weights[:limit]","metadata":{"execution":{"iopub.status.busy":"2023-07-07T20:31:06.096908Z","iopub.execute_input":"2023-07-07T20:31:06.097887Z","iopub.status.idle":"2023-07-07T20:31:06.172498Z","shell.execute_reply.started":"2023-07-07T20:31:06.09785Z","shell.execute_reply":"2023-07-07T20:31:06.171288Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}