{"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":"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-05-22T12:48:27.379358Z","iopub.execute_input":"2023-05-22T12:48:27.379853Z","iopub.status.idle":"2023-05-22T12:48:54.936949Z","shell.execute_reply.started":"2023-05-22T12:48:27.379815Z","shell.execute_reply":"2023-05-22T12:48:54.935462Z"},"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-05-22T12:48:54.940712Z","iopub.execute_input":"2023-05-22T12:48:54.941123Z","iopub.status.idle":"2023-05-22T12:48:54.949611Z","shell.execute_reply.started":"2023-05-22T12:48:54.941087Z","shell.execute_reply":"2023-05-22T12:48:54.948383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-05-22T12:48:54.951192Z","iopub.execute_input":"2023-05-22T12:48:54.95164Z","iopub.status.idle":"2023-05-22T12:48:54.963384Z","shell.execute_reply.started":"2023-05-22T12:48:54.951606Z","shell.execute_reply":"2023-05-22T12:48:54.962049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class color:\n   PURPLE = '\\033[95m'\n   CYAN = '\\033[96m'\n   DARKCYAN = '\\033[36m'\n   BLUE = '\\033[94m'\n   GREEN = '\\033[92m'\n   YELLOW = '\\033[93m'\n   RED = '\\033[91m'\n   BOLD = '\\033[1m'\n   UNDERLINE = '\\033[4m'\n   END = '\\033[0m'","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:48:54.964756Z","iopub.execute_input":"2023-05-22T12:48:54.965159Z","iopub.status.idle":"2023-05-22T12:48:54.976715Z","shell.execute_reply.started":"2023-05-22T12:48:54.965124Z","shell.execute_reply":"2023-05-22T12:48:54.975675Z"},"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-05-22T14:12:04.839635Z","iopub.execute_input":"2023-05-22T14:12:04.84007Z","iopub.status.idle":"2023-05-22T14:12:04.848365Z","shell.execute_reply.started":"2023-05-22T14:12:04.840035Z","shell.execute_reply":"2023-05-22T14:12:04.84704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"subontology_roots = {'BPO':'GO:0008150',\n                     'CCO':'GO:0005575',\n                     'MFO':'GO:0003674'}","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:48:54.99147Z","iopub.execute_input":"2023-05-22T12:48:54.992174Z","iopub.status.idle":"2023-05-22T12:48:55.006384Z","shell.execute_reply.started":"2023-05-22T12:48:54.992138Z","shell.execute_reply":"2023-05-22T12:48:55.005224Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%%time\ngraph = obonet.read_obo(CFG.train_go_obo_path)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:48:55.007592Z","iopub.execute_input":"2023-05-22T12:48:55.008667Z","iopub.status.idle":"2023-05-22T12:49:15.107168Z","shell.execute_reply.started":"2023-05-22T12:48:55.00863Z","shell.execute_reply":"2023-05-22T12:49:15.105796Z"},"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-05-22T12:49:15.108832Z","iopub.execute_input":"2023-05-22T12:49:15.109278Z","iopub.status.idle":"2023-05-22T12:49:15.266555Z","shell.execute_reply.started":"2023-05-22T12:49:15.109237Z","shell.execute_reply":"2023-05-22T12:49:15.265374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Number of nodes: {len(graph)}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:15.26827Z","iopub.execute_input":"2023-05-22T12:49:15.268647Z","iopub.status.idle":"2023-05-22T12:49:15.274695Z","shell.execute_reply.started":"2023-05-22T12:49:15.268615Z","shell.execute_reply":"2023-05-22T12:49:15.273463Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term = \"GO:0000006\"","metadata":{"execution":{"iopub.status.busy":"2023-05-22T14:11:45.676527Z","iopub.execute_input":"2023-05-22T14:11:45.67694Z","iopub.status.idle":"2023-05-22T14:11:45.68215Z","shell.execute_reply.started":"2023-05-22T14:11:45.67691Z","shell.execute_reply":"2023-05-22T14:11:45.680783Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph.nodes[term]","metadata":{"execution":{"iopub.status.busy":"2023-05-22T14:11:46.702672Z","iopub.execute_input":"2023-05-22T14:11:46.703394Z","iopub.status.idle":"2023-05-22T14:11:46.711584Z","shell.execute_reply.started":"2023-05-22T14:11:46.703344Z","shell.execute_reply":"2023-05-22T14:11:46.710611Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_dag(graph, term, radius=10)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T14:12:35.940124Z","iopub.execute_input":"2023-05-22T14:12:35.940667Z","iopub.status.idle":"2023-05-22T14:12:36.225636Z","shell.execute_reply.started":"2023-05-22T14:12:35.940533Z","shell.execute_reply":"2023-05-22T14:12:36.22438Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_dag(graph, term, radius=1000)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:15.595439Z","iopub.execute_input":"2023-05-22T12:49:15.596281Z","iopub.status.idle":"2023-05-22T12:49:15.880599Z","shell.execute_reply.started":"2023-05-22T12:49:15.596245Z","shell.execute_reply":"2023-05-22T12:49:15.879514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"term = \"GO:0044270\"","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:15.886442Z","iopub.execute_input":"2023-05-22T12:49:15.886803Z","iopub.status.idle":"2023-05-22T12:49:15.891856Z","shell.execute_reply.started":"2023-05-22T12:49:15.886774Z","shell.execute_reply":"2023-05-22T12:49:15.890647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"graph.nodes[term]","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:15.893423Z","iopub.execute_input":"2023-05-22T12:49:15.89425Z","iopub.status.idle":"2023-05-22T12:49:15.905208Z","shell.execute_reply.started":"2023-05-22T12:49:15.894214Z","shell.execute_reply":"2023-05-22T12:49:15.904012Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_dag(graph, term, radius=1)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:15.906787Z","iopub.execute_input":"2023-05-22T12:49:15.907157Z","iopub.status.idle":"2023-05-22T12:49:16.194775Z","shell.execute_reply.started":"2023-05-22T12:49:15.907093Z","shell.execute_reply":"2023-05-22T12:49:16.193714Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_dag(graph, term, radius=1000)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:16.196308Z","iopub.execute_input":"2023-05-22T12:49:16.196675Z","iopub.status.idle":"2023-05-22T12:49:16.48074Z","shell.execute_reply.started":"2023-05-22T12:49:16.196646Z","shell.execute_reply":"2023-05-22T12:49:16.479632Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Sequence example:\\n\\n\", next(iter(SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\"))))","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:16.482244Z","iopub.execute_input":"2023-05-22T12:49:16.482612Z","iopub.status.idle":"2023-05-22T12:49:16.490943Z","shell.execute_reply.started":"2023-05-22T12:49:16.48258Z","shell.execute_reply":"2023-05-22T12:49:16.489423Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sequences = SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\")\nnum_sequences = sum(1 for seq in sequences)\n\nprint(\"Number of sequences:\", num_sequences)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:16.492584Z","iopub.execute_input":"2023-05-22T12:49:16.49294Z","iopub.status.idle":"2023-05-22T12:49:18.650643Z","shell.execute_reply.started":"2023-05-22T12:49:16.49291Z","shell.execute_reply":"2023-05-22T12:49:18.649702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sequences = SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\")\navg_length = sum(len(seq) for seq in sequences)/num_sequences\navg_length","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:18.653538Z","iopub.execute_input":"2023-05-22T12:49:18.654764Z","iopub.status.idle":"2023-05-22T12:49:20.885194Z","shell.execute_reply.started":"2023-05-22T12:49:18.654724Z","shell.execute_reply":"2023-05-22T12:49:20.884199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sequences = SeqIO.parse(CFG.train_seq_fasta_path, \"fasta\")\nlengths = [len(seq) for seq in sequences]\nmedian = np.sort(lengths)\nmedian[int(num_sequences/2)]","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:20.886802Z","iopub.execute_input":"2023-05-22T12:49:20.888633Z","iopub.status.idle":"2023-05-22T12:49:23.16754Z","shell.execute_reply.started":"2023-05-22T12:49:20.888575Z","shell.execute_reply":"2023-05-22T12:49:23.166356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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=['goldenrod'],histfunc = \"min\")\nfig.update_layout(\n    title={\n        'text': \"Distribution of protein sequence lengths\",\n        'y':0.95,\n        'x':1.0,\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-05-22T12:49:23.169095Z","iopub.execute_input":"2023-05-22T12:49:23.169543Z","iopub.status.idle":"2023-05-22T12:49:25.658642Z","shell.execute_reply.started":"2023-05-22T12:49:23.169509Z","shell.execute_reply":"2023-05-22T12:49:25.657486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"np.percentile(lengths, 90)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:25.660057Z","iopub.execute_input":"2023-05-22T12:49:25.660798Z","iopub.status.idle":"2023-05-22T12:49:25.698211Z","shell.execute_reply.started":"2023-05-22T12:49:25.660757Z","shell.execute_reply":"2023-05-22T12:49:25.696668Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"taxon_df = pd.read_csv(CFG.train_taxonomy_path,'\\t')\ntaxon_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:25.700111Z","iopub.execute_input":"2023-05-22T12:49:25.700866Z","iopub.status.idle":"2023-05-22T12:49:25.814657Z","shell.execute_reply.started":"2023-05-22T12:49:25.700817Z","shell.execute_reply":"2023-05-22T12:49:25.813413Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"terms_df = pd.read_csv(CFG.train_terms_path,'\\t')\nterms_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:25.81659Z","iopub.execute_input":"2023-05-22T12:49:25.817394Z","iopub.status.idle":"2023-05-22T12:49:28.792181Z","shell.execute_reply.started":"2023-05-22T12:49:25.817347Z","shell.execute_reply":"2023-05-22T12:49:28.791081Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_GO_terms = len(pd.unique(terms_df['term']))\nnum_protein_terms = len(pd.unique(terms_df['EntryID']))\nprint(f\"Number of Unique Proteins: {num_protein_terms}\")\nprint(f\"Number of Unique GO terms: {num_GO_terms}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:28.794126Z","iopub.execute_input":"2023-05-22T12:49:28.79506Z","iopub.status.idle":"2023-05-22T12:49:29.833158Z","shell.execute_reply.started":"2023-05-22T12:49:28.794997Z","shell.execute_reply":"2023-05-22T12:49:29.83186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Number of terms per protein\ndstats_prot = terms_df.groupby('EntryID').agg(term_per_prot=('term', 'nunique')).reset_index()\n\n# Number of proteins per GO term\ndstats_term = terms_df.groupby('term').agg(prot_per_term=('EntryID', 'nunique')).reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:29.83491Z","iopub.execute_input":"2023-05-22T12:49:29.835461Z","iopub.status.idle":"2023-05-22T12:49:34.867049Z","shell.execute_reply.started":"2023-05-22T12:49:29.835418Z","shell.execute_reply":"2023-05-22T12:49:34.865912Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dstats_term.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:34.869397Z","iopub.execute_input":"2023-05-22T12:49:34.870397Z","iopub.status.idle":"2023-05-22T12:49:34.88061Z","shell.execute_reply.started":"2023-05-22T12:49:34.87036Z","shell.execute_reply":"2023-05-22T12:49:34.879412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"dstats_prot.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:34.882351Z","iopub.execute_input":"2023-05-22T12:49:34.882815Z","iopub.status.idle":"2023-05-22T12:49:34.898476Z","shell.execute_reply.started":"2023-05-22T12:49:34.882773Z","shell.execute_reply":"2023-05-22T12:49:34.897297Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colname = 'term_per_prot'\nout = dstats_prot[colname].agg(['min',\n                                lambda x: x.quantile(0.1),\n                                lambda x: x.quantile(0.25),\n                                'median',\n                                'mean',\n                                lambda x: x.quantile(0.75),\n                                lambda x: x.quantile(0.90),\n                                'max']).reset_index()\n\nout.columns = ['stat', colname]\nout[colname] = out[colname].round(0)\nout = out.drop_duplicates(subset='stat')  # Remove duplicates in 'stat' column\nout['stat'] = pd.Categorical(out['stat'], categories=out['stat'], ordered=True)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:34.900269Z","iopub.execute_input":"2023-05-22T12:49:34.900894Z","iopub.status.idle":"2023-05-22T12:49:34.953675Z","shell.execute_reply.started":"2023-05-22T12:49:34.900849Z","shell.execute_reply":"2023-05-22T12:49:34.952532Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"out.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:34.955166Z","iopub.execute_input":"2023-05-22T12:49:34.955537Z","iopub.status.idle":"2023-05-22T12:49:34.969104Z","shell.execute_reply.started":"2023-05-22T12:49:34.955504Z","shell.execute_reply":"2023-05-22T12:49:34.967858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.histogram(x=dstats_prot['term_per_prot'], nbins=200, color_discrete_sequence=['red'])\nfig.update_layout(\n    title={\n        'text': \"Distribution of GO term per protein\",\n        #'y':0.95,\n        #'x':0.5,\n        #'xanchor': 'center',\n        #'yanchor': 'top'\n    },\n    xaxis_title=\"Number of terms per protein\", yaxis_title=\"Count\"\n)\n\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:34.97059Z","iopub.execute_input":"2023-05-22T12:49:34.970935Z","iopub.status.idle":"2023-05-22T12:49:35.071048Z","shell.execute_reply.started":"2023-05-22T12:49:34.970906Z","shell.execute_reply":"2023-05-22T12:49:35.06989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colname = 'prot_per_term'\nout2 = dstats_term[colname].agg(['min',\n                                lambda x: x.quantile(0.1),\n                                lambda x: x.quantile(0.25),\n                                'median',\n                                'mean',\n                                lambda x: x.quantile(0.75),\n                                lambda x: x.quantile(0.90),\n                                'max']).reset_index()\n\nout2.columns = ['stat', colname]\nout2[colname] = out2[colname].round(0)\nout2 = out2.drop_duplicates(subset='stat')  # Remove duplicates in 'stat' column\nout2['stat'] = pd.Categorical(out2['stat'], categories=out2['stat'], ordered=True)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:35.072713Z","iopub.execute_input":"2023-05-22T12:49:35.073542Z","iopub.status.idle":"2023-05-22T12:49:35.105414Z","shell.execute_reply.started":"2023-05-22T12:49:35.073489Z","shell.execute_reply":"2023-05-22T12:49:35.103768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"out2.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:35.106949Z","iopub.execute_input":"2023-05-22T12:49:35.107309Z","iopub.status.idle":"2023-05-22T12:49:35.121208Z","shell.execute_reply.started":"2023-05-22T12:49:35.107276Z","shell.execute_reply":"2023-05-22T12:49:35.120355Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colname = 'prot_per_term'\nout2 = dstats_term[colname].agg(['min',\n                                lambda x: x.quantile(0.1),\n                                lambda x: x.quantile(0.25),\n                                'median',\n                                'mean',\n                                lambda x: x.quantile(0.75),\n                                lambda x: x.quantile(0.90),\n                                'max']).reset_index()\n\nout2.columns = ['stat', colname]\nout2[colname] = out2[colname].round(0)\nout2 = out2.drop_duplicates(subset='stat')  # Remove duplicates in 'stat' column\nout2['stat'] = pd.Categorical(out2['stat'], categories=out2['stat'], ordered=True)\nout2.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:35.122527Z","iopub.execute_input":"2023-05-22T12:49:35.122847Z","iopub.status.idle":"2023-05-22T12:49:35.157237Z","shell.execute_reply.started":"2023-05-22T12:49:35.122819Z","shell.execute_reply":"2023-05-22T12:49:35.155915Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = px.histogram(x=np.log10(dstats_term['prot_per_term']), nbins=50, color_discrete_sequence=['green'])\nfig.update_layout(\n    title={\n        'text': \"Distribution of protein per GO term\",\n        #'y':0.95,\n        #'x':0.5,\n        #'xanchor': 'center',\n        #'yanchor': 'top'\n    },\n    xaxis_title=\"Number of proteins per GO term\", yaxis_title=\"Count\"\n)\nfig.update_layout(xaxis_range=[0,10])\nfig.update_layout(yaxis_range=[0,5000])\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T16:11:27.127988Z","iopub.execute_input":"2023-05-22T16:11:27.128558Z","iopub.status.idle":"2023-05-22T16:11:27.233947Z","shell.execute_reply.started":"2023-05-22T16:11:27.1285Z","shell.execute_reply":"2023-05-22T16:11:27.232818Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"BPO = np.where(terms_df['aspect']=='BPO')\nnp.size(BPO)\nBPO\n#print(f\"Number of BPO terms: {(num_BPO)}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:35.253883Z","iopub.execute_input":"2023-05-22T12:49:35.254395Z","iopub.status.idle":"2023-05-22T12:49:36.25652Z","shell.execute_reply.started":"2023-05-22T12:49:35.254355Z","shell.execute_reply":"2023-05-22T12:49:36.255572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bpo_df = terms_df.loc[terms_df['aspect']=='BPO']\nlen(bpo_df)\nbpo_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:36.257858Z","iopub.execute_input":"2023-05-22T12:49:36.258697Z","iopub.status.idle":"2023-05-22T12:49:37.437968Z","shell.execute_reply.started":"2023-05-22T12:49:36.258659Z","shell.execute_reply":"2023-05-22T12:49:37.436746Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_bpo_proteins = len(pd.unique(bpo_df['EntryID']))\nnum_bpo_terms = len(pd.unique(bpo_df['term']))\nprint(f\"Number of Unique BPO Proteins: {num_bpo_proteins}\") \nprint(f\"Number of Unique BPO Terms: {num_bpo_terms}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:37.439605Z","iopub.execute_input":"2023-05-22T12:49:37.44006Z","iopub.status.idle":"2023-05-22T12:49:38.135355Z","shell.execute_reply.started":"2023-05-22T12:49:37.440018Z","shell.execute_reply":"2023-05-22T12:49:38.134051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cco_df = terms_df.loc[terms_df['aspect']=='CCO']\nlen(cco_df)\ncco_df.head()\nnum_cco_proteins = len(pd.unique(cco_df['EntryID']))\nnum_cco_terms = len(pd.unique(cco_df['term']))\nprint(f\"Number of Unique CCO Proteins: {num_cco_proteins}\") \nprint(f\"Number of Unique CCO Terms: {num_cco_terms}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:38.136666Z","iopub.execute_input":"2023-05-22T12:49:38.137064Z","iopub.status.idle":"2023-05-22T12:49:39.375938Z","shell.execute_reply.started":"2023-05-22T12:49:38.137029Z","shell.execute_reply":"2023-05-22T12:49:39.374779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mfo_df = terms_df.loc[terms_df['aspect']=='MFO']\nlen(mfo_df)\nmfo_df.head()\nnum_mfo_proteins = len(pd.unique(mfo_df['EntryID']))\nnum_mfo_terms = len(pd.unique(mfo_df['term']))\nprint(f\"Number of Unique MFO Proteins: {num_mfo_proteins}\") \nprint(f\"Number of Unique MFO Terms: {num_mfo_terms}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:39.384669Z","iopub.execute_input":"2023-05-22T12:49:39.385066Z","iopub.status.idle":"2023-05-22T12:49:40.488024Z","shell.execute_reply.started":"2023-05-22T12:49:39.385034Z","shell.execute_reply":"2023-05-22T12:49:40.486777Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"terms_df","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:49:40.489345Z","iopub.execute_input":"2023-05-22T12:49:40.489676Z","iopub.status.idle":"2023-05-22T12:49:40.504488Z","shell.execute_reply.started":"2023-05-22T12:49:40.489648Z","shell.execute_reply":"2023-05-22T12:49:40.503256Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bpo_prot = bpo_df.groupby('EntryID').agg(term_per_prot=('term', 'nunique')).reset_index()\nbpo_term = bpo_df.groupby('term').agg(prot_per_term=('EntryID', 'nunique')).reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:49.773594Z","iopub.execute_input":"2023-05-22T18:08:49.77446Z","iopub.status.idle":"2023-05-22T18:08:53.256122Z","shell.execute_reply.started":"2023-05-22T18:08:49.774413Z","shell.execute_reply":"2023-05-22T18:08:53.255079Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bpo_term","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:29:21.371048Z","iopub.execute_input":"2023-05-22T18:29:21.371539Z","iopub.status.idle":"2023-05-22T18:29:21.389711Z","shell.execute_reply.started":"2023-05-22T18:29:21.371507Z","shell.execute_reply":"2023-05-22T18:29:21.388335Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bpo_prot","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:53.258615Z","iopub.execute_input":"2023-05-22T18:08:53.259464Z","iopub.status.idle":"2023-05-22T18:08:53.274707Z","shell.execute_reply.started":"2023-05-22T18:08:53.259414Z","shell.execute_reply":"2023-05-22T18:08:53.273342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mfo_prot = mfo_df.groupby('EntryID').agg(term_per_prot=('term', 'nunique')).reset_index()\nmfo_term = mfo_df.groupby('term').agg(prot_per_term=('EntryID', 'nunique')).reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:53.276176Z","iopub.execute_input":"2023-05-22T18:08:53.276684Z","iopub.status.idle":"2023-05-22T18:08:53.846052Z","shell.execute_reply.started":"2023-05-22T18:08:53.276636Z","shell.execute_reply":"2023-05-22T18:08:53.844812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mfo_prot","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:55.398105Z","iopub.execute_input":"2023-05-22T18:08:55.398761Z","iopub.status.idle":"2023-05-22T18:08:55.411261Z","shell.execute_reply.started":"2023-05-22T18:08:55.398725Z","shell.execute_reply":"2023-05-22T18:08:55.410404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cco_prot = cco_df.groupby('EntryID').agg(term_per_prot=('term', 'nunique')).reset_index()\ncco_term = cco_df.groupby('term').agg(prot_per_term=('EntryID', 'nunique')).reset_index()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:56.089356Z","iopub.execute_input":"2023-05-22T18:08:56.089815Z","iopub.status.idle":"2023-05-22T18:08:57.074815Z","shell.execute_reply.started":"2023-05-22T18:08:56.08978Z","shell.execute_reply":"2023-05-22T18:08:57.073654Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colname = 'prot_per_term'\nout3 = bpo_term[colname].agg(['min',\n                                lambda x: x.quantile(0.1),\n                                lambda x: x.quantile(0.25),\n                                'median',\n                                'mean',\n                                lambda x: x.quantile(0.75),\n                                lambda x: x.quantile(0.90),\n                                'max']).reset_index()\n\nout3.columns = ['stat', colname]\nout3[colname] = out3[colname].round(0)\nout3 = out3.drop_duplicates(subset='stat')  # Remove duplicates in 'stat' column\nout3['stat'] = pd.Categorical(out3['stat'], categories=out3['stat'], ordered=True)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:57.077251Z","iopub.execute_input":"2023-05-22T18:08:57.077678Z","iopub.status.idle":"2023-05-22T18:08:57.101693Z","shell.execute_reply.started":"2023-05-22T18:08:57.077643Z","shell.execute_reply":"2023-05-22T18:08:57.100705Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"out3","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:57.9827Z","iopub.execute_input":"2023-05-22T18:08:57.983766Z","iopub.status.idle":"2023-05-22T18:08:57.995798Z","shell.execute_reply.started":"2023-05-22T18:08:57.983718Z","shell.execute_reply":"2023-05-22T18:08:57.994302Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"out4 = bpo_term['prot_per_term'].agg(lambda x: x.quantile(0.75))\nout4","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:58.327636Z","iopub.execute_input":"2023-05-22T18:08:58.328031Z","iopub.status.idle":"2023-05-22T18:08:58.341012Z","shell.execute_reply.started":"2023-05-22T18:08:58.327999Z","shell.execute_reply":"2023-05-22T18:08:58.339513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(out3.loc[1])","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:58.84802Z","iopub.execute_input":"2023-05-22T18:08:58.848793Z","iopub.status.idle":"2023-05-22T18:08:58.857919Z","shell.execute_reply.started":"2023-05-22T18:08:58.848744Z","shell.execute_reply":"2023-05-22T18:08:58.856415Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def statcalc(df, colname):\n    out = df[colname].agg(['min',\n                                lambda x: x.quantile(0.1),\n                                lambda x: x.quantile(0.25),\n                                'median',\n                                'mean',\n                                lambda x: x.quantile(0.75),\n                                lambda x: x.quantile(0.90),\n                                'max']).reset_index()\n\n    out.columns = ['stat', colname]\n    out[colname] = out[colname].round(0)\n    out = out.drop_duplicates(subset='stat')  # Remove duplicates in 'stat' column\n    out['stat'] = pd.Categorical(out['stat'], categories=out['stat'], ordered=True)\n    return out","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:08:59.809351Z","iopub.execute_input":"2023-05-22T18:08:59.809755Z","iopub.status.idle":"2023-05-22T18:08:59.818595Z","shell.execute_reply.started":"2023-05-22T18:08:59.809722Z","shell.execute_reply":"2023-05-22T18:08:59.817377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bpo_prot","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:13:04.427378Z","iopub.execute_input":"2023-05-22T18:13:04.428562Z","iopub.status.idle":"2023-05-22T18:13:04.444706Z","shell.execute_reply.started":"2023-05-22T18:13:04.42852Z","shell.execute_reply":"2023-05-22T18:13:04.442366Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nimport pandas as pd\n\n# Sample data\ndata = pd.DataFrame({\n    'GO Term': ['GO1', 'GO2', 'GO3', 'GO4', 'GO1', 'GO2', 'GO3', 'GO4', 'GO1', 'GO2', 'GO3', 'GO4'],\n    'Category': ['BP', 'BP', 'BP', 'BP', 'MP', 'MP', 'MP', 'MP', 'CC', 'CC', 'CC', 'CC'],\n    'Protein Count': [10, 20, 15, 5, 12, 8, 6, 10, 30, 25, 20, 15]\n})\n\n# Create violin plot\nfig = px.violin(bpo_term, x='term', y='prot_per_term', color='term', box=True, points='all')\n\n# Update plot layout\nfig.update_layout(\n    title='Proteins per GO Term',\n    xaxis_title='Category',\n    yaxis_title='Number of Proteins'\n)\n\n# Show the plot\nfig.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:30:25.719372Z","iopub.execute_input":"2023-05-22T18:30:25.719795Z","iopub.status.idle":"2023-05-22T18:32:17.794571Z","shell.execute_reply.started":"2023-05-22T18:30:25.719764Z","shell.execute_reply":"2023-05-22T18:32:17.793529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\nfig, (ax1, ax2) = plt.subplots(nrows=1, ncols=2, figsize=(9, 4), sharey=True)\n\n# Plotting violin plot for BPO\nax1.set_title('BPO')\nax1.set_ylabel('Number of Proteins')\nax1.violinplot(bpo_prot)\n\n# Plotting violin plot for MFO\nax2.set_title('MFO')\nax2.violinplot(mfo_prot)\n\n# Set labels for x-axis ticks\nax1.set_xticks([1])\nax1.set_xticklabels(['BPO'])\nax2.set_xticks([1])\nax2.set_xticklabels(['MFO'])\n\n# Add a common y-axis label\nfig.text(0.04, 0.5, 'Number of Proteins', va='center', rotation='vertical')\n\n# Show the plot\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:05:01.135853Z","iopub.execute_input":"2023-05-22T18:05:01.136339Z","iopub.status.idle":"2023-05-22T18:05:01.733153Z","shell.execute_reply.started":"2023-05-22T18:05:01.136291Z","shell.execute_reply":"2023-05-22T18:05:01.731428Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig = plt.figure()\n  \n# Create an axes instance\nax = fig.gca()\nax.violinplot(out3,\n           showmeans=True,\n           showmedians=True)\nax.set_title('violin plot')","metadata":{"execution":{"iopub.status.busy":"2023-05-22T18:03:38.634662Z","iopub.execute_input":"2023-05-22T18:03:38.63507Z","iopub.status.idle":"2023-05-22T18:03:38.973003Z","shell.execute_reply.started":"2023-05-22T18:03:38.63504Z","shell.execute_reply":"2023-05-22T18:03:38.971313Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"terms_df\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:51:07.687956Z","iopub.execute_input":"2023-05-22T12:51:07.688415Z","iopub.status.idle":"2023-05-22T12:51:07.70363Z","shell.execute_reply.started":"2023-05-22T12:51:07.688379Z","shell.execute_reply":"2023-05-22T12:51:07.702378Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"mfo_df\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:51:13.020753Z","iopub.execute_input":"2023-05-22T12:51:13.021429Z","iopub.status.idle":"2023-05-22T12:51:13.036879Z","shell.execute_reply.started":"2023-05-22T12:51:13.021394Z","shell.execute_reply":"2023-05-22T12:51:13.03528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cco_df","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:51:20.402032Z","iopub.execute_input":"2023-05-22T12:51:20.402477Z","iopub.status.idle":"2023-05-22T12:51:20.416351Z","shell.execute_reply.started":"2023-05-22T12:51:20.402441Z","shell.execute_reply":"2023-05-22T12:51:20.415504Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Merge df1 and df2 based on common column(s)\nmerged_df = pd.merge(bpo_df, mfo_df, on='EntryID', how='inner')\n\n# Merge merged_df and df3 based on common column(s)\nfinal_df = pd.merge(merged_df, cco_df, on='EntryID', how='inner')\n\n# The resulting final_df will contain the overlapping data from all three DataFrames\n\n# Find the intersection of the common column(s) in all three DataFrames\ncommon_values = set(bpo_df['EntryID']).intersection(mfo_df['EntryID'], cco_df['EntryID'])\n\n# Filter the original DataFrames based on the common values\n#overlapping_df1 = bpo_df[bpo_df['aspect'].isin(common_values)]\n#overlapping_df2 = mfo_df[mfo_df['aspect'].isin(common_values)]\n#overlapping_df3 = cco_df[cco_df['aspect'].isin(common_values)]\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:51:26.287156Z","iopub.execute_input":"2023-05-22T12:51:26.287582Z","iopub.status.idle":"2023-05-22T12:51:33.296847Z","shell.execute_reply.started":"2023-05-22T12:51:26.28755Z","shell.execute_reply":"2023-05-22T12:51:33.295111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Assuming your data is stored in a DataFrame called 'df' and the column of interest is 'column_name'\nrepetitions = terms_df['EntryID'].duplicated(keep='first')\n\n# Filter the DataFrame to show only the rows with repeated values\nrepeated_rows = terms_df[repetitions]\n\n# Display the repeated rows\nprint(repeated_rows)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T12:56:13.892756Z","iopub.execute_input":"2023-05-22T12:56:13.893173Z","iopub.status.idle":"2023-05-22T12:56:14.461654Z","shell.execute_reply.started":"2023-05-22T12:56:13.893134Z","shell.execute_reply":"2023-05-22T12:56:14.46046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_counts = terms_df.groupby(['EntryID', 'aspect']).size().reset_index(name='Count')\ndf_counts.head(10)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T13:03:59.187968Z","iopub.execute_input":"2023-05-22T13:03:59.188439Z","iopub.status.idle":"2023-05-22T13:04:00.831008Z","shell.execute_reply.started":"2023-05-22T13:03:59.188403Z","shell.execute_reply":"2023-05-22T13:04:00.8299Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_counts = terms_df.pivot_table(index='EntryID', columns='aspect', aggfunc='size', fill_value=0)\ndf_counts = df_counts.reset_index()\n\n# Display the resulting DataFrame\nprint(df_counts)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T13:13:21.658063Z","iopub.execute_input":"2023-05-22T13:13:21.658476Z","iopub.status.idle":"2023-05-22T13:13:23.394759Z","shell.execute_reply.started":"2023-05-22T13:13:21.658443Z","shell.execute_reply":"2023-05-22T13:13:23.39341Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_ontology = []\nfor id in range(len(df_counts)):\n    num_ontology.append(3 - (not int(df_counts['BPO'][id])) - (not int(df_counts['CCO'][id])) - (not int(df_counts['MFO'][id])))\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T13:25:16.302306Z","iopub.execute_input":"2023-05-22T13:25:16.302779Z","iopub.status.idle":"2023-05-22T13:25:21.582505Z","shell.execute_reply.started":"2023-05-22T13:25:16.302745Z","shell.execute_reply":"2023-05-22T13:25:21.581342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"value_counts = pd.Series(num_ontology).value_counts().reset_index()\n\n# Create a DataFrame to store the unique values and their counts\nmeow_df = pd.DataFrame({'Value': value_counts['index'], 'Count': value_counts[0]})\n\n# Display the DataFrame\nprint(meow_df)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T13:32:26.768739Z","iopub.execute_input":"2023-05-22T13:32:26.769165Z","iopub.status.idle":"2023-05-22T13:32:26.851502Z","shell.execute_reply.started":"2023-05-22T13:32:26.769128Z","shell.execute_reply":"2023-05-22T13:32:26.850141Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.pie(meow_df['Count'], labels=meow_df['Value'], autopct='%1.1f%%')\n\n# Add a title\nplt.title('Pie Chart')\n\n# Display the chart\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T13:34:50.624825Z","iopub.execute_input":"2023-05-22T13:34:50.625286Z","iopub.status.idle":"2023-05-22T13:34:50.771216Z","shell.execute_reply.started":"2023-05-22T13:34:50.625249Z","shell.execute_reply":"2023-05-22T13:34:50.769569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install goatools","metadata":{"execution":{"iopub.status.busy":"2023-05-22T13:56:43.136337Z","iopub.execute_input":"2023-05-22T13:56:43.136733Z","iopub.status.idle":"2023-05-22T13:57:37.599286Z","shell.execute_reply.started":"2023-05-22T13:56:43.136704Z","shell.execute_reply":"2023-05-22T13:57:37.597911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from goatools.obo_parser import GODag\n\n# Load the Gene Ontology OBO file\ngo_obo_file = CFG.train_go_obo_path\ngo = GODag(go_obo_file)\n\n# Specify the GO terms for BPO, MFO, and CCO\nbpo_term = 'GO:0008150'  # Biological Process root term\nmfo_term = 'GO:0003674'  # Molecular Function root term\ncco_term = 'GO:0005575'  # Cellular Component root term\n\n# Define the maximum number of nodes per GO level\nmax_nodes_per_level = 3\n\n# Function to recursively build the reduced GO-DAG\ndef build_reduced_dag(go_dag, go_term, max_nodes_per_level):\n    go_node = go_dag[go_term]\n    \n    # Determine the node category based on children\n    if go_node.children and all(child.level == go_node.level + 1 for child in go_node.children):\n        node_category = 'RN'  # Regular node\n    elif go_node.children and any(child.level > go_node.level + 1 for child in go_node.children):\n        node_category = 'JN'  # Jump node\n    else:\n        node_category = 'LN'  # Leaf node\n    \n    # Build the reduced DAG for the current GO term\n    reduced_dag = [(go_node.id, go_node.name, node_category)]\n    children = sorted(go_node.children, key=lambda x: x.name)\n    for child in children[:max_nodes_per_level]:\n        reduced_dag.extend(build_reduced_dag(go_dag, child.id, max_nodes_per_level))\n    \n    return reduced_dag\n\n# Build the reduced DAGs for BPO, MFO, and CCO\nreduced_dag_bpo = build_reduced_dag(go, bpo_term, max_nodes_per_level)\nreduced_dag_mfo = build_reduced_dag(go, mfo_term, max_nodes_per_level)\nreduced_dag_cco = build_reduced_dag(go, cco_term, max_nodes_per_level)\n\n# Print the reduced DAGs\nprint(\"BPO Reduced DAG:\")\nfor go_id, go_name, node_category in reduced_dag_bpo:\n    print(f\"{go_id} - {go_name} ({node_category})\")\nprint()\nprint(\"MFO Reduced DAG:\")\nfor go_id, go_name, node_category in reduced_dag_mfo:\n    print(f\"{go_id} - {go_name} ({node_category})\")\nprint()\nprint(\"CCO Reduced DAG:\")\nfor go_id, go_name, node_category in reduced_dag_cco:\n    print(f\"{go_id} - {go_name} ({node_category})\")\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T14:53:58.00523Z","iopub.execute_input":"2023-05-22T14:53:58.005764Z","iopub.status.idle":"2023-05-22T14:53:59.777961Z","shell.execute_reply.started":"2023-05-22T14:53:58.005729Z","shell.execute_reply":"2023-05-22T14:53:59.776728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import networkx as nx\nimport matplotlib.pyplot as plt\nfrom goatools.obo_parser import GODag\n\n# Load the Gene Ontology OBO file\ngo_obo_file = CFG.train_go_obo_path\ngo = GODag(go_obo_file)\n\n# Specify the GO terms for BPO, MFO, and CCO\nbpo_term = 'GO:0008150'  # Biological Process root term\nmfo_term = 'GO:0003674'  # Molecular Function root term\ncco_term = 'GO:0005575'  # Cellular Component root term\n\n# Define the maximum number of nodes per GO level\nmax_nodes_per_level = 3\n\n# Function to recursively build the reduced GO-DAG as a graph\ndef build_reduced_dag_graph(go_dag, go_term, max_nodes_per_level):\n    go_node = go_dag[go_term]\n    \n    # Determine the node category based on children\n    if go_node.children and all(child.level == go_node.level + 1 for child in go_node.children):\n        node_category = 'RN'  # Regular node\n    elif go_node.children and any(child.level > go_node.level + 1 for child in go_node.children):\n        node_category = 'JN'  # Jump node\n    else:\n        node_category = 'LN'  # Leaf node\n    \n    # Build the reduced DAG graph for the current GO term\n    graph = nx.DiGraph()\n    graph.add_node(go_node.id, name=go_node.name, category=node_category)\n    children = sorted(go_node.children, key=lambda x: x.name)\n    for child in children[:max_nodes_per_level]:\n        child_graph = build_reduced_dag_graph(go_dag, child.id, max_nodes_per_level)\n        graph = nx.compose(graph, child_graph)\n        graph.add_edge(go_node.id, child.id)\n    \n    return graph\n\n# Build the reduced DAG graphs for BPO, MFO, and CCO\ngraph_bpo = build_reduced_dag_graph(go, bpo_term, max_nodes_per_level)\ngraph_mfo = build_reduced_dag_graph(go, mfo_term, max_nodes_per_level)\ngraph_cco = build_reduced_dag_graph(go, cco_term, max_nodes_per_level)\n\n# Set node positions using a hierarchical layout\npos_bpo = nx.kamada_kawai_layout(graph_bpo)\npos_mfo = nx.kamada_kawai_layout(graph_mfo)\npos_cco = nx.kamada_kawai_layout(graph_cco)\n\n# Plot the reduced DAG graphs\nplt.figure(figsize=(30, 10))\n\nplt.subplot(131)\nnx.draw(graph_bpo, pos=pos_bpo, with_labels=True, node_size=500, node_color='lightblue', edge_color='gray')\nplt.title('BPO Reduced DAG')\n\nplt.subplot(132)\nnx.draw(graph_mfo, pos=pos_mfo, with_labels=True, node_size=500, node_color='lightgreen', edge_color='gray')\nplt.title('MFO Reduced DAG')\n\nplt.subplot(133)\nnx.draw(graph_cco, pos=pos_cco, with_labels=True, node_size=500, node_color='lightyellow', edge_color='gray')\nplt.title('CCO Reduced DAG')\n\nplt.tight_layout()\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T17:50:24.960865Z","iopub.execute_input":"2023-05-22T17:50:24.96203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"levels_bpo = nx.shortest_path_length(graph_bpo, bpo_term)\nlevels_mfo = nx.shortest_path_length(graph_mfo, mfo_term)\nlevels_cco = nx.shortest_path_length(graph_cco, cco_term)\n\n# Print the levels for each graph\nprint(\"Levels for BPO:\")\nfor node, level in levels_bpo.items():\n    print(f\"{node}: Level {level}\")\n\nprint(\"\\nLevels for MFO:\")\nfor node, level in levels_mfo.items():\n    print(f\"{node}: Level {level}\")\n\nprint(\"\\nLevels for CCO:\")\nfor node, level in levels_cco.items():\n    print(f\"{node}: Level {level}\")","metadata":{"execution":{"iopub.status.busy":"2023-05-22T15:05:10.341754Z","iopub.execute_input":"2023-05-22T15:05:10.342757Z","iopub.status.idle":"2023-05-22T15:05:10.355935Z","shell.execute_reply.started":"2023-05-22T15:05:10.3427Z","shell.execute_reply":"2023-05-22T15:05:10.354155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_bpo = pd.DataFrame.from_dict(levels_bpo, orient='index', columns=['Level'])\ndf_mfo = pd.DataFrame.from_dict(levels_mfo, orient='index', columns=['Level'])\ndf_cco = pd.DataFrame.from_dict(levels_cco, orient='index', columns=['Level'])\n\n# Perform descriptive statistics on the levels\nstats_bpo = df_bpo.describe()\nstats_mfo = df_mfo.describe()\nstats_cco = df_cco.describe()\n\n# Create a summary DataFrame\ndf_summary = pd.DataFrame({'BPO': [stats_bpo.loc['count', 'Level'], stats_bpo.loc['mean', 'Level'],\n                                   stats_bpo.loc['min', 'Level'], stats_bpo.loc['max', 'Level']],\n                           'MFO': [stats_mfo.loc['count', 'Level'], stats_mfo.loc['mean', 'Level'],\n                                   stats_mfo.loc['min', 'Level'], stats_mfo.loc['max', 'Level']],\n                           'CCO': [stats_cco.loc['count', 'Level'], stats_cco.loc['mean', 'Level'],\n                                   stats_cco.loc['min', 'Level'], stats_cco.loc['max', 'Level']]},\n                          index=['Protein Count', 'Mean Level', 'Min Level', 'Max Level'])\n\n# Print the summary DataFrame\nprint(df_summary)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T15:09:07.799926Z","iopub.execute_input":"2023-05-22T15:09:07.800403Z","iopub.status.idle":"2023-05-22T15:09:07.830624Z","shell.execute_reply.started":"2023-05-22T15:09:07.800372Z","shell.execute_reply":"2023-05-22T15:09:07.829277Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"levels_bpo = list(levels_bpo.values())\nlevels_mfo = list(levels_mfo.values())\nlevels_cco = list(levels_cco.values())\n\n# Plot the histograms\nplt.figure(figsize=(10, 6))\n\nplt.subplot(131)\nplt.hist(levels_bpo, bins=range(min(levels_bpo), max(levels_bpo)+2), edgecolor='black', alpha=0.75)\nplt.xlabel('BPO DAG Levels')\nplt.ylabel('Number of Proteins')\nplt.title('Distribution of Proteins vs BPO DAG Levels')\n\nplt.subplot(132)\nplt.hist(levels_mfo, bins=range(min(levels_mfo), max(levels_mfo)+2), edgecolor='black', alpha=0.75)\nplt.xlabel('MFO DAG Levels')\nplt.ylabel('Number of Proteins')\nplt.title('Distribution of Proteins vs MFO DAG Levels')\n\nplt.subplot(133)\nplt.hist(levels_cco, bins=range(min(levels_cco), max(levels_cco)+2), edgecolor='black', alpha=0.75)\nplt.xlabel('CCO DAG Levels')\nplt.ylabel('Number of Proteins')\nplt.title('Distribution of Proteins vs CCO DAG Levels')\n\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T15:13:13.913262Z","iopub.execute_input":"2023-05-22T15:13:13.91376Z","iopub.status.idle":"2023-05-22T15:13:14.719859Z","shell.execute_reply.started":"2023-05-22T15:13:13.913728Z","shell.execute_reply":"2023-05-22T15:13:14.718437Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"go_obo_file = CFG.train_go_obo_path\ngo = GODag(go_obo_file)\ndef calculate_term_percentage(go_dag, ontology):\n    term_levels = {}\n    for term in go_dag.values():\n        if term.namespace == ontology and not term.children: # Consider only leaf nodes\n            level = term.level\n            term_levels[level] = term_levels.get(level, 0) + 1\n    total_terms = sum(term_levels.values())\n    term_percentages = {level: count / total_terms * 100 for level, count in term_levels.items()}\n    return term_percentages","metadata":{"execution":{"iopub.status.busy":"2023-05-22T15:57:05.835206Z","iopub.execute_input":"2023-05-22T15:57:05.835739Z","iopub.status.idle":"2023-05-22T15:57:08.888594Z","shell.execute_reply.started":"2023-05-22T15:57:05.835707Z","shell.execute_reply":"2023-05-22T15:57:08.887412Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"ontologies = ['molecular_function', 'biological_process', 'cellular_component']\nterm_percentages = {}\nfor ontology in ontologies:\n    term_percentages[ontology] = calculate_term_percentage(go, ontology)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T15:57:12.192255Z","iopub.execute_input":"2023-05-22T15:57:12.192922Z","iopub.status.idle":"2023-05-22T15:57:12.261087Z","shell.execute_reply.started":"2023-05-22T15:57:12.192886Z","shell.execute_reply":"2023-05-22T15:57:12.259982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(12, 4))\n\nfor i, ontology in enumerate(ontologies):\n    levels = list(term_percentages[ontology].keys())\n    percentages = list(term_percentages[ontology].values())\n    plt.subplot(1, 3, i+1)\n    plt.bar(levels, percentages, color='skyblue')\n    plt.xlabel('GO Term Level')\n    plt.ylabel('Percentage')\n    plt.title(f'Percentage of {ontology} GO Terms from Each Level in Train Dataset')\n    plt.xticks(levels)\n    plt.ylim(0, max(percentages) + 5)\n    plt.grid(axis='y', linestyle='--')\n    plt.tight_layout()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-22T15:58:33.062151Z","iopub.execute_input":"2023-05-22T15:58:33.062592Z","iopub.status.idle":"2023-05-22T15:58:34.076626Z","shell.execute_reply.started":"2023-05-22T15:58:33.062556Z","shell.execute_reply":"2023-05-22T15:58:34.075657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from goatools.obo_parser import GODag\nimport pandas as pd\n\n# Load the GO OBO file\ngo_obo_file = CFG.train_go_obo_path # Replace with your GO OBO file path\ngo = GODag(go_obo_file)\n\n# Read the protein annotation file into a pandas DataFrame\nprotein_annotation_file = CFG.train_terms_path  # Replace with your protein annotation file path\ndf_annotations = pd.read_csv(protein_annotation_file, sep='\\t', header=None, names=['EntryID', 'GO Term', 'Aspect'])\n\n# Group the annotations by aspect (BPO, MFO, CCO)\ngrouped_annotations = df_annotations.groupby('Aspect')\n\n# Define a function to count the frequency of GO terms\ndef count_go_terms(go_terms):\n    go_term_counts = {}\n    for go_term in go_terms:\n        if go_term in go_term_counts:\n            go_term_counts[go_term] += 1\n        else:\n            go_term_counts[go_term] = 1\n    return go_term_counts\n\n# Iterate over the aspects (BPO, MFO, CCO)\nfor aspect, annotations in grouped_annotations:\n    print(f\"Aspect: {aspect}\")\n    go_terms = annotations['GO Term'].tolist()\n    go_term_counts = count_go_terms(go_terms)\n    \n    # Find the most frequent GO term\n    most_frequent_go_term = max(go_term_counts, key=go_term_counts.get)\n    most_frequent_count = go_term_counts[most_frequent_go_term]\n    print(f\"Most Frequent GO Term: {most_frequent_go_term}, Count: {most_frequent_count}\")\n    \n    # Find the least frequent GO term\n    least_frequent_go_term = min(go_term_counts, key=go_term_counts.get)\n    least_frequent_count = go_term_counts[least_frequent_go_term]\n    print(f\"Least Frequent GO Term: {least_frequent_go_term}, Count: {least_frequent_count}\")\n    print()\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T16:17:51.228531Z","iopub.execute_input":"2023-05-22T16:17:51.228989Z","iopub.status.idle":"2023-05-22T16:18:00.286413Z","shell.execute_reply.started":"2023-05-22T16:17:51.228953Z","shell.execute_reply":"2023-05-22T16:18:00.284485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from goatools.obo_parser import GODag\nimport pandas as pd\n\n# Load the GO OBO file\ngo_obo_file = CFG.train_go_obo_path  # Replace with your GO OBO file path\ngo = GODag(go_obo_file)\n\n# Read the protein annotation file into a pandas DataFrame\nprotein_annotation_file = CFG.train_terms_path  # Replace with your protein annotation file path\ndf_annotations = pd.read_csv(protein_annotation_file, sep='\\t', header=None, names=['EntryID', 'GO Term', 'Aspect'])\n\n# Count the number of proteins associated with each GO term\ngo_term_counts = df_annotations['GO Term'].value_counts()\n\n# Calculate the total number of proteins\ntotal_proteins = len(df_annotations['EntryID'].unique())\n\n# Filter the GO terms that are associated with at least half of the proteins\nfiltered_go_terms_half = go_term_counts[go_term_counts >= total_proteins / 2]\n\n# Filter the GO terms that are associated with at least one third of the proteins\nfiltered_go_terms_third = go_term_counts[go_term_counts >= total_proteins / 3]\n\n# Retrieve the GO term information for at least half of the proteins\ngo_terms_half = []\nfor go_term in filtered_go_terms_half.index:\n    go_term_info = go[go_term]\n    go_terms_half.append({\n        'GO Term': go_term,\n        'Name': go_term_info.name,\n        'Namespace': go_term_info.namespace,\n        'Count': filtered_go_terms_half[go_term]\n    })\n\n# Retrieve the GO term information for at least one third of the proteins\ngo_terms_third = []\nfor go_term in filtered_go_terms_third.index:\n    go_term_info = go[go_term]\n    go_terms_third.append({\n        'GO Term': go_term,\n        'Name': go_term_info.name,\n        'Namespace': go_term_info.namespace,\n        'Count': filtered_go_terms_third[go_term]\n    })\n\n# Create pandas DataFrames for the GO terms\ndf_go_terms_half = pd.DataFrame(go_terms_half)\ndf_go_terms_third = pd.DataFrame(go_terms_third)\n\n# Print the GO terms associated with at least half of the proteins\nprint(\"GO Terms Associated with at least Half of the Proteins:\")\nprint(df_go_terms_half)\n\n# Print the GO terms associated with at least one third of the proteins\nprint(\"GO Terms Associated with at least One Third of the Proteins:\")\nprint(df_go_terms_third)","metadata":{"execution":{"iopub.status.busy":"2023-05-22T16:26:56.6731Z","iopub.execute_input":"2023-05-22T16:26:56.673608Z","iopub.status.idle":"2023-05-22T16:27:04.213406Z","shell.execute_reply.started":"2023-05-22T16:26:56.673572Z","shell.execute_reply":"2023-05-22T16:27:04.21209Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load the GO OBO file\ngo_obo_file = CFG.train_go_obo_path  # Replace with your GO OBO file path\ngo = GODag(go_obo_file)\n\n# Read the protein annotation file into a pandas DataFrame\nprotein_annotation_file = CFG.train_terms_path  # Replace with your protein annotation file path\ndf_annotations = pd.read_csv(protein_annotation_file, sep='\\t', header=None, names=['EntryID', 'term', 'Aspect'])\n\n# Count the number of proteins associated with each GO term\ngo_term_counts = df_annotations['term'].value_counts()\n\n# Filter the GO terms associated with 100 or fewer proteins\nfiltered_go_terms_100 = go_term_counts[go_term_counts <= 100]\n\n# Filter the GO terms associated with 1 or fewer proteins\nfiltered_go_terms_1 = go_term_counts[go_term_counts <= 1]\n\n# Filter the GO terms associated with 10 or fewer proteins\nfiltered_go_terms_10 = go_term_counts[go_term_counts <= 10]\n\n# Retrieve the GO term information for 100 or fewer proteins\ngo_terms_100 = []\nfor go_term in filtered_go_terms_100.index:\n    go_term_info = go[go_term]\n    go_terms_100.append({\n        'GO Term': go_term,\n        'Name': go_term_info.name,\n        'Namespace': go_term_info.namespace,\n        'Count': filtered_go_terms_100[go_term]\n    })\n\n# Retrieve the GO term information for 1 or fewer proteins\ngo_terms_1 = []\nfor go_term in filtered_go_terms_1.index:\n    go_term_info = go[go_term]\n    go_terms_1.append({\n        'GO Term': go_term,\n        'Name': go_term_info.name,\n        'Namespace': go_term_info.namespace,\n        'Count': filtered_go_terms_1[go_term]\n    })\n\n# Retrieve the GO term information for 10 or fewer proteins\ngo_terms_10 = []\nfor go_term in filtered_go_terms_10.index:\n    go_term_info = go[go_term]\n    go_terms_10.append({\n        'GO Term': go_term,\n        'Name': go_term_info.name,\n        'Namespace': go_term_info.namespace,\n        'Count': filtered_go_terms_10[go_term]\n    })\n\n# Create pandas DataFrames for the GO terms\ndf_go_terms_100 = pd.DataFrame(go_terms_100)\ndf_go_terms_1 = pd.DataFrame(go_terms_1)\ndf_go_terms_10 = pd.DataFrame(go_terms_10)\n\n# Print the GO terms associated with 100 or fewer proteins\nprint(\"GO Terms Associated with 100 or Fewer Proteins:\")\nprint(df_go_terms_100)\n\n# Print the GO terms associated with 1 or fewer proteins\nprint(\"GO Terms Associated with 1 or Fewer Proteins:\")\nprint(df_go_terms_1)\n\n# Print the GO terms associated with 10 or fewer proteins\nprint(\"GO Terms Associated with 10 or Fewer Proteins:\")\nprint(df_go_terms_10)\n","metadata":{"execution":{"iopub.status.busy":"2023-05-22T17:09:42.582807Z","iopub.execute_input":"2023-05-22T17:09:42.583256Z","iopub.status.idle":"2023-05-22T17:09:49.802311Z","shell.execute_reply.started":"2023-05-22T17:09:42.583221Z","shell.execute_reply":"2023-05-22T17:09:49.800711Z"},"trusted":true},"execution_count":null,"outputs":[]}]}