{"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":"#import libraiares:\nimport tensorflow as tf\nimport pandas as pd\nfrom collections import Counter\nimport numpy as np\n\n# we will use bio python to read the sequences\nfrom Bio import SeqIO\n\nimport seaborn as sns\nimport matplotlib.pyplot as plt","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:30.653264Z","iopub.execute_input":"2023-06-23T10:57:30.653912Z","iopub.status.idle":"2023-06-23T10:57:30.660217Z","shell.execute_reply.started":"2023-06-23T10:57:30.653882Z","shell.execute_reply":"2023-06-23T10:57:30.658895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Problem definition:\n1. The goal of this competition is to predict the function of a set of proteins. \n2. The accurate assignment of biological function to the protein is key to understanding life at the molecular level. However, assigning function to any specific protein can be made difficult due to the multiple functions many proteins have, along with their ability to interact with multiple partners. \n3. We are provided with a set of protein sequences on which the participants are asked to predict Gene Ontology (GO) terms in each of the three subontologies: Molecular Function (MF), Biological Process (BP), and Cellular Component (CC). This set of sequences is referred to as test superset.\n\n# Data:\n\n1. train_sequences.fasta - amino acid sequences for proteins in training set\n2. testsuperset.fasta - amino acid sequences for proteins on which the predictions should be made\n3. train_terms.tsv - the training set of proteins and corresponding annotated GO terms\n4. train_taxonomy.tsv - taxon ID for proteins in training set\n5. go-basic.obo - ontology graph structure\n6. testsuperset-taxon-list.tsv - taxon ID for proteins in test superset (Note: you may need to use encoding=\"ISO-8859-1\" to read this file in pandas)\n7. IA.txt - Information Accretion for each term. This is used to weight precision and recall (see Evaluation)\n8. sample_submission.csv - a sample submission file in the correct format\n\n","metadata":{}},{"cell_type":"markdown","source":"# 1. train_sequences.fasta\n\n1. train_sequences.fasta contains the protein sequences for the training dataset.\n2. This 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.\n3. In a FASTA file, each sequence starts with a unique identifier, which is like a name or label for that particular sequence. It helps scientists keep track of different sequences they are working with. After the identifier, you have the actual genetic code represented by letters: A, C, G, and T for DNA sequences, or letters representing different amino acids for protein sequences.\n\nWe'll use Biopython.SeqIO to work with the FASTA files. then we will convert it into pandas dataframe \"df_sequences\".\n\n# df_sequence dataframe:\n It has 4 columns and 142246 instances.\n1. ID: 142246 unique ids.\n2. Sequence- 138924 protein sequences.\n3. description- 142246\n\nA look at a single sequence \"ID: P20536\":\n\n    1. \"ID: P20536\"- unique identifier of the sequence.\n    2. \"Name: P20536\" - \"P20536\" might be the name associated with the sequence.\n    3. \"Description\"-  field provides additional details about the sequence, including its source (Vaccinia virus), strain (Copenhagen), organism, and other information.\n          a. OX == taxon ID\n          b. GN == gene name\n          c. PE == the protein existance level\n          d. SV == sequence version. In this casem \"1\" indicates this is the first  version of the sequence\n          e. \"Number of features: 0\" indicates that there are no additional features or annotations associated with this particular sequence.\n          f. \"Seq('MNSVTVSHAPYTITYHDDWEPVMSQLVEFYNEVASWLLRDETSPIPDKFFIQLK...FIY')\" represents the actual sequence data. The sequence is displayed using the \"Seq()\" object, with the specific sequence letters shown as a string.\n          \n3. sequence_length: this is a derived feature to get an insight into distribution of sequences length.\n   1. average length is 554 amino acids\n   2. min length is 3 amino acids \n   3. max length is 35375 amino acids","metadata":{}},{"cell_type":"code","source":"# path to the train and test fasta files\ntrain_fasta = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\ntest_fasta = '/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset.fasta'","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:30.661329Z","iopub.execute_input":"2023-06-23T10:57:30.661713Z","iopub.status.idle":"2023-06-23T10:57:30.673374Z","shell.execute_reply.started":"2023-06-23T10:57:30.661683Z","shell.execute_reply":"2023-06-23T10:57:30.672215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# read train and test fasta files\ntrain_sequences = SeqIO.parse(train_fasta, 'fasta')\n#test_sequences = SeqIO.parse(test_fasta, 'fasta')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:30.676211Z","iopub.execute_input":"2023-06-23T10:57:30.676728Z","iopub.status.idle":"2023-06-23T10:57:30.714922Z","shell.execute_reply.started":"2023-06-23T10:57:30.676685Z","shell.execute_reply":"2023-06-23T10:57:30.713692Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sequences_dict = SeqIO.to_dict(train_sequences)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:30.716831Z","iopub.execute_input":"2023-06-23T10:57:30.717371Z","iopub.status.idle":"2023-06-23T10:57:34.895536Z","shell.execute_reply.started":"2023-06-23T10:57:30.717324Z","shell.execute_reply":"2023-06-23T10:57:34.894254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# take a look at the first three instances:\nfirst_three_keys = list(sequences_dict.keys())[:3]\nfirst_three_instances = {key: sequences_dict[key] for key in first_three_keys}\n\n# Printing the first three instances\nfor seq_id, seq_record in first_three_instances.items():\n    print(\"ID:\", seq_id)\n    print(\"Sequence:\", seq_record.seq)\n    print('Description:', seq_record.description)\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:34.11593Z","iopub.execute_input":"2023-06-23T07:08:34.116272Z","iopub.status.idle":"2023-06-23T07:08:34.132898Z","shell.execute_reply.started":"2023-06-23T07:08:34.116241Z","shell.execute_reply":"2023-06-23T07:08:34.131514Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create a dataframe of the above fasta files dictionary:\ndata = []\nfor seq_id, seq_record in sequences_dict.items():\n    sequence = str(seq_record.seq)\n    seq_des = str(seq_record.description)\n    data.append({\"ID\": seq_id, \"Sequence\": sequence, \"description\": seq_des})\n\ndf_sequences = pd.DataFrame(data)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:34.134274Z","iopub.execute_input":"2023-06-23T07:08:34.134667Z","iopub.status.idle":"2023-06-23T07:08:34.611973Z","shell.execute_reply.started":"2023-06-23T07:08:34.134622Z","shell.execute_reply":"2023-06-23T07:08:34.611052Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:34.613384Z","iopub.execute_input":"2023-06-23T07:08:34.614193Z","iopub.status.idle":"2023-06-23T07:08:34.640647Z","shell.execute_reply.started":"2023-06-23T07:08:34.614159Z","shell.execute_reply":"2023-06-23T07:08:34.639903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:34.641615Z","iopub.execute_input":"2023-06-23T07:08:34.642618Z","iopub.status.idle":"2023-06-23T07:08:34.648283Z","shell.execute_reply.started":"2023-06-23T07:08:34.642588Z","shell.execute_reply":"2023-06-23T07:08:34.647177Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:34.649744Z","iopub.execute_input":"2023-06-23T07:08:34.650191Z","iopub.status.idle":"2023-06-23T07:08:35.068305Z","shell.execute_reply.started":"2023-06-23T07:08:34.650154Z","shell.execute_reply":"2023-06-23T07:08:35.067263Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create a column with length of sequences:\ndf_sequences['sequence_len'] = df_sequences['Sequence'].apply(len)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:35.072472Z","iopub.execute_input":"2023-06-23T07:08:35.072778Z","iopub.status.idle":"2023-06-23T07:08:35.134471Z","shell.execute_reply.started":"2023-06-23T07:08:35.072753Z","shell.execute_reply":"2023-06-23T07:08:35.133356Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:35.135884Z","iopub.execute_input":"2023-06-23T07:08:35.136243Z","iopub.status.idle":"2023-06-23T07:08:35.147296Z","shell.execute_reply.started":"2023-06-23T07:08:35.136216Z","shell.execute_reply":"2023-06-23T07:08:35.146019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:35.149063Z","iopub.execute_input":"2023-06-23T07:08:35.149514Z","iopub.status.idle":"2023-06-23T07:08:35.178227Z","shell.execute_reply.started":"2023-06-23T07:08:35.149476Z","shell.execute_reply":"2023-06-23T07:08:35.177196Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# distribution of sequence lengths:\n\n# Adjust the dimensions of the plot\nfig, ax = plt.subplots(figsize=(8, 4)) \n\n# Set a specific color palette\n#sns.set_palette('husl')\n\n# Create the histogram plot\nsns.histplot(df_sequences, x = 'sequence_len', kde=False, color = 'blue')\n\n# Add rugplot for outliers\nsns.rugplot(df_sequences, height=0.03, color='red')\n\n# create title\nplt.title('Distribution of length Protein Sequences in Train data set')\n\n# Show the plot\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:35.179657Z","iopub.execute_input":"2023-06-23T07:08:35.180043Z","iopub.status.idle":"2023-06-23T07:08:41.701395Z","shell.execute_reply.started":"2023-06-23T07:08:35.180012Z","shell.execute_reply":"2023-06-23T07:08:41.700371Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Log distribution of length of protein sequences:\nplt.hist(np.log10(df_sequences['sequence_len']), bins=\"fd\")\nplt.title(\"Distribution of lengths of the protein sequences\")\nplt.xlabel(\"log$_1$$_0$(Length of the protein sequence)\")\nplt.ylabel(\"Number of protein sequences\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:41.702978Z","iopub.execute_input":"2023-06-23T07:08:41.703621Z","iopub.status.idle":"2023-06-23T07:08:42.330335Z","shell.execute_reply.started":"2023-06-23T07:08:41.703584Z","shell.execute_reply":"2023-06-23T07:08:42.328275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The length of a protein appears to follow a log-normal distribution. The number of sequences with length more than 5K are nearly negligible although the length goes upto 35K.","metadata":{}},{"cell_type":"markdown","source":"# 2. test fasts files:\nThe test superset is a set of protein sequences on which the participants are asked to predict GO terms.\n\nWe'll use Biopython.SeqIO to get the 'test' FASTA files. Then we will convert it into pandas dataframe \"df_sequences_test\".\n\n# df_sequences_test\nIt has data for 141864 sequences \n1. ID- 141864 unique ids\n2. Sequence- 139155 protien sequences\n3. description\n4. sequence_len- a feature enginered col \n   1. average protein length is about 477 amino acids\n   2. min protein length is about 2 amino acids.\n   3. max protein length is about 35213 amino acids.\n\n\n\n","metadata":{}},{"cell_type":"code","source":"# read test fasta files\ntest_sequences = SeqIO.parse(test_fasta, 'fasta')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:42.332777Z","iopub.execute_input":"2023-06-23T07:08:42.333252Z","iopub.status.idle":"2023-06-23T07:08:42.343016Z","shell.execute_reply.started":"2023-06-23T07:08:42.333211Z","shell.execute_reply":"2023-06-23T07:08:42.34186Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Remove duplicates based on ID\ntest_sequences = {rec.id: rec for rec in test_sequences}.values()\n\n# Convert the sequences to a dictionary\nsequences_dict_test = SeqIO.to_dict(test_sequences)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:42.344391Z","iopub.execute_input":"2023-06-23T07:08:42.344728Z","iopub.status.idle":"2023-06-23T07:08:44.823476Z","shell.execute_reply.started":"2023-06-23T07:08:42.344698Z","shell.execute_reply":"2023-06-23T07:08:44.822483Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create a dataframe of the above fasta files dictionary:\ndata_test = []\nfor seq_id, seq_record in sequences_dict_test.items():\n    sequence = str(seq_record.seq)\n    seq_des = str(seq_record.description)\n    data_test.append({\"ID\": seq_id, \"Sequence\": sequence, \"description\": seq_des})\n\ndf_sequences_test = pd.DataFrame(data_test)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:44.824726Z","iopub.execute_input":"2023-06-23T07:08:44.825131Z","iopub.status.idle":"2023-06-23T07:08:45.266981Z","shell.execute_reply.started":"2023-06-23T07:08:44.825096Z","shell.execute_reply":"2023-06-23T07:08:45.266049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences_test.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:45.268303Z","iopub.execute_input":"2023-06-23T07:08:45.26865Z","iopub.status.idle":"2023-06-23T07:08:45.275465Z","shell.execute_reply.started":"2023-06-23T07:08:45.268621Z","shell.execute_reply":"2023-06-23T07:08:45.274503Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create a column with length of sequences:\ndf_sequences_test['sequence_len'] = df_sequences_test['Sequence'].apply(len)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:45.276725Z","iopub.execute_input":"2023-06-23T07:08:45.277058Z","iopub.status.idle":"2023-06-23T07:08:45.348237Z","shell.execute_reply.started":"2023-06-23T07:08:45.27703Z","shell.execute_reply":"2023-06-23T07:08:45.346936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences_test.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:45.349587Z","iopub.execute_input":"2023-06-23T07:08:45.349946Z","iopub.status.idle":"2023-06-23T07:08:45.360457Z","shell.execute_reply.started":"2023-06-23T07:08:45.349916Z","shell.execute_reply":"2023-06-23T07:08:45.359711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences_test.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:08:45.361412Z","iopub.execute_input":"2023-06-23T07:08:45.361925Z","iopub.status.idle":"2023-06-23T07:08:45.72172Z","shell.execute_reply.started":"2023-06-23T07:08:45.361895Z","shell.execute_reply":"2023-06-23T07:08:45.720544Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create a column with length of sequences:\ndf_sequences_test['sequence_len'] = df_sequences_test['Sequence'].apply(len)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:09:48.268476Z","iopub.execute_input":"2023-06-23T07:09:48.268887Z","iopub.status.idle":"2023-06-23T07:09:48.331675Z","shell.execute_reply.started":"2023-06-23T07:09:48.268851Z","shell.execute_reply":"2023-06-23T07:09:48.330557Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sequences_test.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:09:48.923307Z","iopub.execute_input":"2023-06-23T07:09:48.923991Z","iopub.status.idle":"2023-06-23T07:09:48.945447Z","shell.execute_reply.started":"2023-06-23T07:09:48.923953Z","shell.execute_reply":"2023-06-23T07:09:48.944489Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# distribution of sequence lengths :\n\n# Adjust the dimensions of the plot\nfig, ax = plt.subplots(figsize=(8, 4)) \n\n# Set a specific color palette\n#sns.set_palette('husl')\n\n# Create the histogram plot\nsns.histplot(df_sequences_test, x = 'sequence_len', kde=False, color='blue')\n\n# Add rugplot for outliers\nsns.rugplot(df_sequences_test, height=0.03, color='red')\n\n# create title\nplt.title('Distribution of length of Protein Sequences in Test data set')\n\n# Show the plot\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:09:49.596434Z","iopub.execute_input":"2023-06-23T07:09:49.59678Z","iopub.status.idle":"2023-06-23T07:09:57.02745Z","shell.execute_reply.started":"2023-06-23T07:09:49.596754Z","shell.execute_reply":"2023-06-23T07:09:57.026244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Log distribution of length of protein sequences:\nplt.hist(np.log10(df_sequences_test['sequence_len']), bins=\"fd\")\nplt.title(\"Distribution of lengths of the protein sequences in Test data\")\nplt.xlabel(\"log$_1$$_0$(Length of the protein sequence in TEst data)\")\nplt.ylabel(\"Number of protein sequences in test data\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:09:57.029284Z","iopub.execute_input":"2023-06-23T07:09:57.02965Z","iopub.status.idle":"2023-06-23T07:09:57.639407Z","shell.execute_reply.started":"2023-06-23T07:09:57.029619Z","shell.execute_reply":"2023-06-23T07:09:57.638208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# A comparison of train vs test sequences:\n1. A look at the fasta files show:\n\n   I. Train fasta file:\n\n   1. Number of ids:  142246\n   2. Number of sequences:  142246\n   3. Number of unique ids:  142246\n   4. Number of unique sequences:  138924\n\n   II. Test fasta file:\n\n   1. Number of ids:  141865\n   2. Number of sequences:  141865\n   3. Number of unique ids:  141864\n   4. Number of unique sequences:  139155\n\n2. Both train and test protein sequences show a similar distribtuion of protein sequence lengths.\n\n3. Train and Test sequences show similar amino acid composition.\n\nHere are a few observations that can be made from the obtained frequency values:\n\n1. The most common amino acids in this dataset are leucine (L), serine (S), alanine (A), and glycine (G). These amino acids are known to be abundant in proteins and play important roles in protein structure and function.\n\n2. The least common amino acids in this dataset are cysteine (C), methionine (M), tryptophan (W), and histidine (H). These amino acids are typically less abundant in proteins, but they can be important for specific functions, such as catalysis, metal binding, or protein-protein interactions.\n\n3. The presence of the amino acid selenocysteine (U) in the dataset suggests that some of the proteins may be selenoproteins, which contain selenium in the form of selenocysteine instead of cysteine.\n\n4. The presence of ambiguous amino acids (X, B, Z) and rare amino acids (O, U) in the dataset suggests that some of the sequences may be incomplete or contain errors.","metadata":{}},{"cell_type":"code","source":"# amino acid composition of  train 'sequence':\namino_acid_count_train=df_sequences[\"Sequence\"].str.split(\"\").explode(\"Sequence\").value_counts().drop(\"\")\namino_acid_count_train\n","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:09:57.640721Z","iopub.execute_input":"2023-06-23T07:09:57.641047Z","iopub.status.idle":"2023-06-23T07:10:14.980572Z","shell.execute_reply.started":"2023-06-23T07:09:57.641021Z","shell.execute_reply":"2023-06-23T07:10:14.979533Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# amino acid composition of  test 'sequence':\namino_acid_count_test=df_sequences_test[\"Sequence\"].str.split(\"\").explode(\"Sequence\").value_counts().drop(\"\")\namino_acid_count_test\n","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:10:14.982799Z","iopub.execute_input":"2023-06-23T07:10:14.983178Z","iopub.status.idle":"2023-06-23T07:10:30.439761Z","shell.execute_reply.started":"2023-06-23T07:10:14.983148Z","shell.execute_reply":"2023-06-23T07:10:30.43854Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# create dataframes from the series:\ntrain_aa = pd.DataFrame(amino_acid_count_train)\ntrain_aa = train_aa.reset_index().rename(columns={'index': 'amino acid', 'Sequence': 'train_seq'})\n\ntest_aa = pd.DataFrame(amino_acid_count_test)\ntest_aa = test_aa.reset_index().rename(columns={'index': 'amino acid', 'Sequence': 'test_seq'})\n","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:10:30.441352Z","iopub.execute_input":"2023-06-23T07:10:30.441769Z","iopub.status.idle":"2023-06-23T07:10:30.451126Z","shell.execute_reply.started":"2023-06-23T07:10:30.441728Z","shell.execute_reply":"2023-06-23T07:10:30.44986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# merge the datfrems and then melt it to get amino acid compsition:\nmerged_df = pd.merge(train_aa, test_aa, on='amino acid')\nmelted_df = merged_df.melt('amino acid', var_name='sequence_label', value_name='seq_count')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:10:30.452612Z","iopub.execute_input":"2023-06-23T07:10:30.453052Z","iopub.status.idle":"2023-06-23T07:10:30.477618Z","shell.execute_reply.started":"2023-06-23T07:10:30.453014Z","shell.execute_reply":"2023-06-23T07:10:30.476645Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# dict from the letters to their name\namino_acid_dict = {\n    'A': 'Alanine',\n    'R': 'Arginine',\n    'N': 'Asparagine',\n    'D': 'Aspartic Acid',\n    'C': 'Cysteine',\n    'E': 'Glutamic Acid',\n    'Q': 'Glutamine',\n    'G': 'Glycine',\n    'H': 'Histidine',\n    'I': 'Isoleucine',\n    'L': 'Leucine',\n    'K': 'Lysine',\n    'M': 'Methionine',\n    'F': 'Phenylalanine',\n    'P': 'Proline',\n    'S': 'Serine',\n    'T': 'Threonine',\n    'W': 'Tryptophan',\n    'Y': 'Tyrosine',\n    'V': 'Valine',\n    'X': 'Any/Unknown',\n    'O': 'Pyrrolysine',\n    'U': 'Selenocysteine',\n    'B': 'Asparagine or Aspartic Acid',\n    'Z': 'Glutamine or Glutamic Acid',\n}","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:10:30.478866Z","iopub.execute_input":"2023-06-23T07:10:30.479268Z","iopub.status.idle":"2023-06-23T07:10:30.485592Z","shell.execute_reply.started":"2023-06-23T07:10:30.479241Z","shell.execute_reply":"2023-06-23T07:10:30.484585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# map the name to amino acid dictionary:\nmelted_df['amino_acid_names'] = melted_df['amino acid'].map(amino_acid_dict)\nmelted_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:10:30.48687Z","iopub.execute_input":"2023-06-23T07:10:30.487864Z","iopub.status.idle":"2023-06-23T07:10:30.512768Z","shell.execute_reply.started":"2023-06-23T07:10:30.487811Z","shell.execute_reply":"2023-06-23T07:10:30.51158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Plot the bar graph using Seaborn\nsns.barplot(data=melted_df, y='amino_acid_names', x='seq_count', hue='sequence_label', orient ='h')\n\n# Set labels and title\nplt.title('Comparison of Amio Acid compositions of  Train and Test Protein-Sequences')\n\nplt.xlabel('Amino Acid Names')\nplt.ylabel('amino acid count')\n\n# Show the plot\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:10:30.514152Z","iopub.execute_input":"2023-06-23T07:10:30.514463Z","iopub.status.idle":"2023-06-23T07:10:31.170023Z","shell.execute_reply.started":"2023-06-23T07:10:30.514436Z","shell.execute_reply":"2023-06-23T07:10:31.168982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. train_terms.tsv \n\n1. contains the list of annotated terms (ground truth) for the proteins in train_sequences.fasta. \n2. It has 3 features:\n    1. EntryID: indicates the protein's UniProt accession ID, \n    2. term: the GO term ID.\n    3. aspect: indicates in which ontology the term appears. BPO, CCO, and MFO are abbreviations for different categories of gene ontology terms. \n    \n           BPO: Biological Process Ontology, which describes biological processes, functions, and pathways. \n    \n           CCO: Cellular Component Ontology, which describes the components of a cell or its extracellular environment. \n    \n           MFO: Molecular Function Ontology, which describes the biochemical activities or capabilities of proteins and other molecules. \n    \n    \nThese categories are used in gene ontology to classify genes and gene products based on their biological roles and functions. By using these categories, researchers can better understand the functions and interactions of different genes and gene products within a biological system.","metadata":{}},{"cell_type":"code","source":"#Load the Dataset\ntrain_terms = pd.read_csv(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_terms.tsv\",sep=\"\\t\")\nprint(train_terms.shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:50.294388Z","iopub.execute_input":"2023-06-23T10:57:50.295168Z","iopub.status.idle":"2023-06-23T10:57:54.298023Z","shell.execute_reply.started":"2023-06-23T10:57:50.29513Z","shell.execute_reply":"2023-06-23T10:57:54.296629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:54.3004Z","iopub.execute_input":"2023-06-23T10:57:54.300915Z","iopub.status.idle":"2023-06-23T10:57:54.336363Z","shell.execute_reply.started":"2023-06-23T10:57:54.300868Z","shell.execute_reply":"2023-06-23T10:57:54.334936Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:54.338206Z","iopub.execute_input":"2023-06-23T10:57:54.338711Z","iopub.status.idle":"2023-06-23T10:57:55.888418Z","shell.execute_reply.started":"2023-06-23T10:57:54.338665Z","shell.execute_reply":"2023-06-23T10:57:55.887053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_terms.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:57:55.890633Z","iopub.execute_input":"2023-06-23T10:57:55.890976Z","iopub.status.idle":"2023-06-23T10:58:00.397127Z","shell.execute_reply.started":"2023-06-23T10:57:55.890947Z","shell.execute_reply":"2023-06-23T10:58:00.395814Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#distribution of 'aspect' in train_terms:\nsns.countplot(data = train_terms, y= 'aspect', orient ='h')\nplt.title('Distribution of GO labels ')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:00.398699Z","iopub.execute_input":"2023-06-23T10:58:00.399194Z","iopub.status.idle":"2023-06-23T10:58:06.637551Z","shell.execute_reply.started":"2023-06-23T10:58:00.39915Z","shell.execute_reply":"2023-06-23T10:58:06.636096Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Count the number of unique terms per aspect (ontology)\nunique_terms_per_aspect = train_terms.groupby('aspect')['term'].nunique()\n\nprint(\"Unique terms per\", unique_terms_per_aspect)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:06.639639Z","iopub.execute_input":"2023-06-23T10:58:06.640142Z","iopub.status.idle":"2023-06-23T10:58:08.78708Z","shell.execute_reply.started":"2023-06-23T10:58:06.6401Z","shell.execute_reply":"2023-06-23T10:58:08.785869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colors = [ 'blue','orange','green']\nunique_terms_per_aspect.plot(kind='barh', color = colors)\nplt.title('Distribution of unique terms in labels ')\nplt.xlabel('Count of unique terms')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:08.790156Z","iopub.execute_input":"2023-06-23T10:58:08.791117Z","iopub.status.idle":"2023-06-23T10:58:09.146012Z","shell.execute_reply.started":"2023-06-23T10:58:08.791073Z","shell.execute_reply":"2023-06-23T10:58:09.144757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create df of GO terms per aspect\nhistograms = train_terms.groupby('aspect')['term'].value_counts().reset_index(name='count')\nhistograms","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:09.14727Z","iopub.execute_input":"2023-06-23T10:58:09.147601Z","iopub.status.idle":"2023-06-23T10:58:12.114486Z","shell.execute_reply.started":"2023-06-23T10:58:09.147575Z","shell.execute_reply":"2023-06-23T10:58:12.113136Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#create dfs based on aspect/subontologies:\nhistograms_bpo = histograms[histograms['aspect' ]== 'BPO'].reset_index(drop = True)\nhistograms_cco = histograms[histograms['aspect' ]== 'CCO'].reset_index(drop = True)\nhistograms_mfo = histograms[histograms['aspect' ]== 'MFO'].reset_index(drop = True)\n#rename count in different df:\nhistograms_bpo=histograms_bpo.rename(columns={'count': 'terms_count_bpo'})\nhistograms_cco=histograms_cco.rename(columns={'count': 'terms_count_cco'})\nhistograms_mfo=histograms_mfo.rename(columns={'count': 'terms_count_mfo'})","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:12.119244Z","iopub.execute_input":"2023-06-23T10:58:12.119689Z","iopub.status.idle":"2023-06-23T10:58:12.167469Z","shell.execute_reply.started":"2023-06-23T10:58:12.119655Z","shell.execute_reply":"2023-06-23T10:58:12.166439Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#top twenty term in differnt labels!\nbpo_20_term = histograms_bpo.head(20)\ncco_20_term = histograms_cco.head(20)\nmfo_20_term = histograms_mfo.head(20)\n\n#plot top ten term in all the three labels:\nbpo_20_term.plot(kind='line', x = 'term', color ='blue')\nplt.title('Top 20 Terms in labels ')\nplt.xticks(range(len(bpo_20_term['term'])), bpo_20_term['term'], rotation='vertical')\n\ncco_20_term.plot(kind='line', x = 'term', color='orange')\nplt.xticks(range(len(cco_20_term['term'])), cco_20_term['term'], rotation='vertical')\n\nmfo_20_term.plot(kind='line', x = 'term', color='green')\nplt.xticks(range(len(mfo_20_term['term'])), mfo_20_term['term'], rotation='vertical')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:26.153567Z","iopub.execute_input":"2023-06-23T10:58:26.153973Z","iopub.status.idle":"2023-06-23T10:58:27.399426Z","shell.execute_reply.started":"2023-06-23T10:58:26.153941Z","shell.execute_reply":"2023-06-23T10:58:27.398168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3. testsuperset-taxon-list.tsv\n\nIt contains taxon ID for proteins in test superset. And also gives an insight as to which species each 'taxon id' represent. We haev 90 'taxon ids' represented in this data set.","metadata":{}},{"cell_type":"code","source":"testsupertest_taxon = pd.read_csv('/kaggle/input/cafa-5-protein-function-prediction/Test (Targets)/testsuperset-taxon-list.tsv', sep='\\t', error_bad_lines=False, encoding= 'unicode_escape')\ndisplay(testsupertest_taxon)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.309946Z","iopub.execute_input":"2023-06-23T10:58:13.310447Z","iopub.status.idle":"2023-06-23T10:58:13.334533Z","shell.execute_reply.started":"2023-06-23T10:58:13.310402Z","shell.execute_reply":"2023-06-23T10:58:13.333221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"testsupertest_taxon.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.336287Z","iopub.execute_input":"2023-06-23T10:58:13.338108Z","iopub.status.idle":"2023-06-23T10:58:13.353078Z","shell.execute_reply.started":"2023-06-23T10:58:13.338056Z","shell.execute_reply":"2023-06-23T10:58:13.351674Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4. train_taxonomy.tsv \n1. contains the list of proteins and the species to which they belong, represented by a \"taxonomic identifier\" (taxon ID) number. \n2. There are 3156 unique taxon ids. We can see a large number of proteins coming from few taxon ids(species), and a long tail of species with very few proteins. \n\nDistribution of 'taxon_ids' shows that the prominent taxon of protein sequences are from humans, Arabidopsis(model organism widely used in plant biology research), mouse(lab model for mammals:\n   1. 9606- The taxon ID 9606 corresponds to the scientific name Homo sapiens, which represents the species commonly known as human.\n   2. 3702- The taxon ID 3702 corresponds to the scientific name Arabidopsis thaliana, which represents the species commonly known as thale cress. \n   3. 10090 - The taxon ID 10090 corresponds to the scientific name Mus musculus, which represents the species commonly known as the house mouse.\n   4. 7955 -The taxon ID 7955 corresponds to the scientific name Danio rerio, which represents the species commonly known as zebrafish. \n   5. 7227 - The taxon ID 7227 actually corresponds to the scientific name Drosophila melanogaster, which represents the species commonly known as fruit fly or vinegar fly. \n","metadata":{}},{"cell_type":"code","source":"#Load the Dataset\ntrain_tax = pd.read_csv(\"/kaggle/input/cafa-5-protein-function-prediction/Train/train_taxonomy.tsv\",sep=\"\\t\")\nprint(train_tax.shape)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.355431Z","iopub.execute_input":"2023-06-23T10:58:13.356083Z","iopub.status.idle":"2023-06-23T10:58:13.485299Z","shell.execute_reply.started":"2023-06-23T10:58:13.356037Z","shell.execute_reply":"2023-06-23T10:58:13.484226Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_tax.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.486867Z","iopub.execute_input":"2023-06-23T10:58:13.487345Z","iopub.status.idle":"2023-06-23T10:58:13.501396Z","shell.execute_reply.started":"2023-06-23T10:58:13.487301Z","shell.execute_reply":"2023-06-23T10:58:13.499903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_tax.info()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.502932Z","iopub.execute_input":"2023-06-23T10:58:13.503422Z","iopub.status.idle":"2023-06-23T10:58:13.570232Z","shell.execute_reply.started":"2023-06-23T10:58:13.503377Z","shell.execute_reply":"2023-06-23T10:58:13.568963Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_tax.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.573953Z","iopub.execute_input":"2023-06-23T10:58:13.575344Z","iopub.status.idle":"2023-06-23T10:58:13.637271Z","shell.execute_reply.started":"2023-06-23T10:58:13.575287Z","shell.execute_reply":"2023-06-23T10:58:13.636095Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"top_Tax_id = train_tax.taxonomyID.value_counts().head(30)\ntop_Tax_id","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.639151Z","iopub.execute_input":"2023-06-23T10:58:13.640543Z","iopub.status.idle":"2023-06-23T10:58:13.654289Z","shell.execute_reply.started":"2023-06-23T10:58:13.640492Z","shell.execute_reply":"2023-06-23T10:58:13.652877Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#distribution of 'taxonomyID's Top twenty:\ntop_Tax_id.plot(kind = 'bar')\nplt.xlabel('Taxon IDs')\nplt.ylabel('Count')\nplt.title('Distribution of Taxon IDs ')","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:13.65579Z","iopub.execute_input":"2023-06-23T10:58:13.656212Z","iopub.status.idle":"2023-06-23T10:58:14.244643Z","shell.execute_reply.started":"2023-06-23T10:58:13.656178Z","shell.execute_reply":"2023-06-23T10:58:14.243234Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#least frequent taxon ids\nleast_Tax_id = train_tax.taxonomyID.value_counts().tail(10)\nleast_Tax_id","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:14.246204Z","iopub.execute_input":"2023-06-23T10:58:14.24757Z","iopub.status.idle":"2023-06-23T10:58:14.260727Z","shell.execute_reply.started":"2023-06-23T10:58:14.247524Z","shell.execute_reply":"2023-06-23T10:58:14.259547Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Merge the train labels and the taxonomy data frame.\n\nA look at how proteins are annonated by GO terms shows:\n1. Average number of terms associated with proteins is 37\n2. Min number of Terms associated with a protein is 2 \n3. Max number of Terms associated with a protein is 815.\n\nThe distribution is heavily skewed to the right, and the number of proteins having more than 150-200 terms associated to them are negligible.\n\nDistribution of proteirns in labels:\n\n","metadata":{}},{"cell_type":"code","source":"# Merge the train labels and the taxonomy data frame.\ntraining_df = pd.merge(train_terms,train_tax, on= \"EntryID\")\ntraining_df.head()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:54.988164Z","iopub.execute_input":"2023-06-23T10:58:54.988572Z","iopub.status.idle":"2023-06-23T10:58:56.937171Z","shell.execute_reply.started":"2023-06-23T10:58:54.988543Z","shell.execute_reply":"2023-06-23T10:58:56.935967Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_df.shape","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:56.939289Z","iopub.execute_input":"2023-06-23T10:58:56.940489Z","iopub.status.idle":"2023-06-23T10:58:56.949271Z","shell.execute_reply.started":"2023-06-23T10:58:56.94045Z","shell.execute_reply":"2023-06-23T10:58:56.946898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"training_df.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T10:58:57.272032Z","iopub.execute_input":"2023-06-23T10:58:57.272713Z","iopub.status.idle":"2023-06-23T10:58:58.917784Z","shell.execute_reply.started":"2023-06-23T10:58:57.272677Z","shell.execute_reply":"2023-06-23T10:58:58.916099Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#A look at how proteins are annotated by terms:\nprot_term_df = training_df.groupby('EntryID')['term'].count()\n\nprot_term_df","metadata":{"execution":{"iopub.status.busy":"2023-06-23T11:28:04.528615Z","iopub.execute_input":"2023-06-23T11:28:04.529075Z","iopub.status.idle":"2023-06-23T11:28:06.015085Z","shell.execute_reply.started":"2023-06-23T11:28:04.52904Z","shell.execute_reply":"2023-06-23T11:28:06.013899Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"prot_term_df.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T11:27:52.560341Z","iopub.execute_input":"2023-06-23T11:27:52.560789Z","iopub.status.idle":"2023-06-23T11:27:52.577993Z","shell.execute_reply.started":"2023-06-23T11:27:52.560756Z","shell.execute_reply":"2023-06-23T11:27:52.57673Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#distribution of unique terms over proetin ids:\nun_term_per_prot = training_df.groupby('EntryID')['term'].nunique()\n\nplt.hist(un_term_per_prot,bins=\"fd\")\nplt.title(\"Distribution of Number of Terms Annotating a protein\")\nplt.xlabel(\"# of terms per unit protein\")\nplt.ylabel(\"frequency\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-23T11:33:53.319187Z","iopub.execute_input":"2023-06-23T11:33:53.319656Z","iopub.status.idle":"2023-06-23T11:33:57.458638Z","shell.execute_reply.started":"2023-06-23T11:33:53.319625Z","shell.execute_reply":"2023-06-23T11:33:57.457316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#distribution of log of unique terms over proetin ids:\nun_term_per_prot = training_df.groupby('EntryID')['term'].nunique()\n\nplt.hist(np.log10(un_term_per_prot),bins=\"fd\")\nplt.xlabel(\"log$_1$$_0$(# of proteins per term)\")\nplt.xlabel(\"# of terms per unit protein\")\nplt.ylabel(\"frequency\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2023-06-23T11:36:09.613571Z","iopub.execute_input":"2023-06-23T11:36:09.614001Z","iopub.status.idle":"2023-06-23T11:36:13.082596Z","shell.execute_reply.started":"2023-06-23T11:36:09.613971Z","shell.execute_reply":"2023-06-23T11:36:13.081284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#distribution of proteins in labels 'aspect':\nprot_label_un = training_df.groupby('aspect')['EntryID'].nunique()\n#prot_label_un = prot_label_un.reset_index()\n#prot_label_un.rename(columns={'EntryID': 'protein_count'})","metadata":{"execution":{"iopub.status.busy":"2023-06-23T11:56:30.213946Z","iopub.execute_input":"2023-06-23T11:56:30.215453Z","iopub.status.idle":"2023-06-23T11:56:31.483045Z","shell.execute_reply.started":"2023-06-23T11:56:30.215393Z","shell.execute_reply":"2023-06-23T11:56:31.480766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"colors = [ 'blue','orange','green']\nprot_label_un.plot(kind ='barh', color = colors )\nplt.title('Distribution of Proteins in labels')\nplt.xlabel('unique Protein id count')\n","metadata":{"execution":{"iopub.status.busy":"2023-06-23T11:58:43.016229Z","iopub.execute_input":"2023-06-23T11:58:43.016655Z","iopub.status.idle":"2023-06-23T11:58:43.302905Z","shell.execute_reply.started":"2023-06-23T11:58:43.016624Z","shell.execute_reply":"2023-06-23T11:58:43.302033Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5. go-basic.obo - ontology graph structure\n\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. \nFor example, the obonet package is available for Python. \n\nThe nodes in this graph are indexed by the term name, for example the roots of the three onotlogies are:\n\nsubontology_roots =\n                    {'BPO':'GO:0008150',\n\n                     'CCO':'GO:0005575',\n                   \n                     'MFO':'GO:0003674'}\n                     \n                     \n Further analysis show, \n \n    i) total number of nodes in ontology of proeteins are 43248\n    ii) There are 84805 Edges in the graph\n    iii) The GO ontology is commonly represented as a Directed Acyclic Graph (DAG). In this graph structure, each term or node in the ontology is represented as a vertex, and the relationships between terms are represented as directed edges. The directed edges indicate the \"is-a\" or \"part-of\" relationships between the terms. The DAG structure allows for multiple parents and multiple children for a given term, capturing the hierarchical and overlapping relationships in the ontology.\n\ncredits: https://www.kaggle.com/code/leonidkulyk/eda-cafa5-pfp-interactive-dags-plotly/notebook\n               ","metadata":{}},{"cell_type":"code","source":"!pip install obonet","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:11:30.282194Z","iopub.execute_input":"2023-06-23T07:11:30.282574Z","iopub.status.idle":"2023-06-23T07:11:42.62915Z","shell.execute_reply.started":"2023-06-23T07:11:30.282546Z","shell.execute_reply":"2023-06-23T07:11:42.627762Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install pygraphviz","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:11:42.631849Z","iopub.execute_input":"2023-06-23T07:11:42.632343Z","iopub.status.idle":"2023-06-23T07:11:57.534007Z","shell.execute_reply.started":"2023-06-23T07:11:42.63229Z","shell.execute_reply":"2023-06-23T07:11:57.532391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#libraireis to study GO ontologies:\n# networkx for graph stuff\nimport networkx\n\n# to plot subgraphs\nfrom networkx.drawing.nx_agraph import graphviz_layout\n\n# to plot graphs\nimport matplotlib.pyplot as plt\n\n# use obonet to load the ontology in a networkx graph\nimport obonet","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:11:57.53567Z","iopub.execute_input":"2023-06-23T07:11:57.536066Z","iopub.status.idle":"2023-06-23T07:11:57.800792Z","shell.execute_reply.started":"2023-06-23T07:11:57.536022Z","shell.execute_reply":"2023-06-23T07:11:57.799757Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"go_path = \"/kaggle/input/cafa-5-protein-function-prediction/Train/go-basic.obo\"","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:11:57.802963Z","iopub.execute_input":"2023-06-23T07:11:57.803344Z","iopub.status.idle":"2023-06-23T07:11:57.808589Z","shell.execute_reply.started":"2023-06-23T07:11:57.803312Z","shell.execute_reply":"2023-06-23T07:11:57.807344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"knowledge_graph = obonet.read_obo(go_path)\n# print the number of nodes and edges in the graph\nprint(\"Nodes: {}\".format(len(knowledge_graph.nodes)))\nprint(\"Edges: {}\".format(len(knowledge_graph.edges)))","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:11:57.810111Z","iopub.execute_input":"2023-06-23T07:11:57.810466Z","iopub.status.idle":"2023-06-23T07:12:05.31676Z","shell.execute_reply.started":"2023-06-23T07:11:57.810438Z","shell.execute_reply":"2023-06-23T07:12:05.315909Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# check if the graph is a DAG\nprint(\"Is DAG: {}\".format(networkx.is_directed_acyclic_graph(knowledge_graph)))","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:16.283837Z","iopub.execute_input":"2023-06-19T05:55:16.284358Z","iopub.status.idle":"2023-06-19T05:55:16.787408Z","shell.execute_reply.started":"2023-06-19T05:55:16.284312Z","shell.execute_reply":"2023-06-19T05:55:16.786198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# using networkx we can get the descendants of a node\nterm = \"GO:0003700\"\nprint(\"Descendants of {}: {}\".format(term, len(networkx.descendants(knowledge_graph, term))))","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:16.78919Z","iopub.execute_input":"2023-06-19T05:55:16.789918Z","iopub.status.idle":"2023-06-19T05:55:16.797517Z","shell.execute_reply.started":"2023-06-19T05:55:16.789876Z","shell.execute_reply":"2023-06-19T05:55:16.796267Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# using the above we can get a subgraph for a given term containing the term and all its descendants\ndef get_subgraph(graph, term):\n    descendants = networkx.descendants(graph, term)\n    return networkx.subgraph(graph, [term] + list(descendants))","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:16.80069Z","iopub.execute_input":"2023-06-19T05:55:16.801099Z","iopub.status.idle":"2023-06-19T05:55:16.810839Z","shell.execute_reply.started":"2023-06-19T05:55:16.801065Z","shell.execute_reply":"2023-06-19T05:55:16.809914Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get node atrributes\nknowledge_graph.nodes[term]","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:16.812099Z","iopub.execute_input":"2023-06-19T05:55:16.813185Z","iopub.status.idle":"2023-06-19T05:55:16.828165Z","shell.execute_reply.started":"2023-06-19T05:55:16.813148Z","shell.execute_reply":"2023-06-19T05:55:16.826701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# integrate the name attribute into our plots. (similar other attributes could be added)\n\ndef get_node_labels(graph):\n    labels = {}\n    for node in graph.nodes:\n        id = node\n        namespace = graph.nodes[node].get('name', '')  # Default to empty string if 'name' is not found\n        labels[node] = f'{id}\\n{namespace}'\n    return labels\n\ndef plot_subgraph(graph, term, width=10, height=10):\n    sg = get_subgraph(graph, term)\n    pos = graphviz_layout(sg, prog='dot')\n    labels = get_node_labels(sg)\n    plt.figure(figsize=(width, height))\n    networkx.draw_networkx(sg, pos, labels=labels, with_labels=True, node_color = 'lightblue', node_size=1000)","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:16.831409Z","iopub.execute_input":"2023-06-19T05:55:16.831876Z","iopub.status.idle":"2023-06-19T05:55:16.842994Z","shell.execute_reply.started":"2023-06-19T05:55:16.831842Z","shell.execute_reply":"2023-06-19T05:55:16.841616Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test\nplot_subgraph(knowledge_graph, term, 15, 15)","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:16.844765Z","iopub.execute_input":"2023-06-19T05:55:16.845197Z","iopub.status.idle":"2023-06-19T05:55:17.363755Z","shell.execute_reply.started":"2023-06-19T05:55:16.845164Z","shell.execute_reply":"2023-06-19T05:55:17.362242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# mapping from GO ID to name\nid_to_name = {id_: data.get('name') for id_, data in knowledge_graph.nodes(data=True)}\n# mapping from name to GO ID\nname_to_id = {data['name']: id_ for id_, data in knowledge_graph.nodes(data=True) if 'name' in data}\n","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:17.365244Z","iopub.execute_input":"2023-06-19T05:55:17.365569Z","iopub.status.idle":"2023-06-19T05:55:17.459715Z","shell.execute_reply.started":"2023-06-19T05:55:17.365541Z","shell.execute_reply":"2023-06-19T05:55:17.458242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# test\nid_to_name[\"GO:0006302\"], name_to_id['double-strand break repair']","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:17.461384Z","iopub.execute_input":"2023-06-19T05:55:17.46177Z","iopub.status.idle":"2023-06-19T05:55:17.471821Z","shell.execute_reply.started":"2023-06-19T05:55:17.46174Z","shell.execute_reply":"2023-06-19T05:55:17.470325Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# the name of the roots.\nid_to_name['GO:0008150'], id_to_name['GO:0005575'], id_to_name['GO:0003674']","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:17.475798Z","iopub.execute_input":"2023-06-19T05:55:17.476282Z","iopub.status.idle":"2023-06-19T05:55:17.487081Z","shell.execute_reply.started":"2023-06-19T05:55:17.476218Z","shell.execute_reply":"2023-06-19T05:55:17.485529Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get the corresponding graph\nplot_subgraph(knowledge_graph, 'GO:0003899', 15, 15)","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:17.488877Z","iopub.execute_input":"2023-06-19T05:55:17.48932Z","iopub.status.idle":"2023-06-19T05:55:18.167817Z","shell.execute_reply.started":"2023-06-19T05:55:17.489286Z","shell.execute_reply":"2023-06-19T05:55:18.166442Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we can reverse the graph as well to get the ancestors of a given term\nreversed_graph = knowledge_graph.reverse(copy=True)\n# plot the ancestors of GO:0003899\nplot_subgraph(reversed_graph, 'GO:0003899', 30, 30)","metadata":{"execution":{"iopub.status.busy":"2023-06-19T05:55:18.169406Z","iopub.execute_input":"2023-06-19T05:55:18.169818Z","iopub.status.idle":"2023-06-19T05:55:21.704441Z","shell.execute_reply.started":"2023-06-19T05:55:18.169782Z","shell.execute_reply":"2023-06-19T05:55:21.703003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 6. IA.txt - \nInformation Accretion for each term. It 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 of the competition.","metadata":{}},{"cell_type":"code","source":"limit = 100\nfile_path = \"/kaggle/input/cafa-5-protein-function-prediction/IA.txt\"\nfile = open(file_path, \"r\")\nfile_contents =file.read()\nfile.close()\n\nprint(file_contents[:limit])","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:18.000234Z","iopub.execute_input":"2023-06-23T07:12:18.000632Z","iopub.status.idle":"2023-06-23T07:12:18.021852Z","shell.execute_reply.started":"2023-06-23T07:12:18.0006Z","shell.execute_reply":"2023-06-23T07:12:18.020778Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Read the text file\ndata_ia = pd.read_csv(\"/kaggle/input/cafa-5-protein-function-prediction/IA.txt\", delimiter='\\t', header = 0)  # Replace 'file.txt' with your file name and adjust the delimiter if necessary\n\n# Display the DataFrame\nprint(data_ia)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:18.735487Z","iopub.execute_input":"2023-06-23T07:12:18.736405Z","iopub.status.idle":"2023-06-23T07:12:18.774978Z","shell.execute_reply.started":"2023-06-23T07:12:18.736358Z","shell.execute_reply":"2023-06-23T07:12:18.773661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# clean data by making index as df.loc[1]:\ndata_ia.loc[-1] = data_ia.columns  # Append the column names as new row\ndata_ia.sort_index(inplace=True)\ndata_ia.reset_index(inplace=True, drop = True)\n\n#assign names to columns:\ndata_ia.columns = [\"term\",\"WeightAssigned\"]\n\n# Display the DataFrame\nprint(data_ia)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:19.483532Z","iopub.execute_input":"2023-06-23T07:12:19.484229Z","iopub.status.idle":"2023-06-23T07:12:19.502464Z","shell.execute_reply.started":"2023-06-23T07:12:19.48419Z","shell.execute_reply":"2023-06-23T07:12:19.501374Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_ia.dtypes","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:20.619224Z","iopub.execute_input":"2023-06-23T07:12:20.620215Z","iopub.status.idle":"2023-06-23T07:12:20.628127Z","shell.execute_reply.started":"2023-06-23T07:12:20.620181Z","shell.execute_reply":"2023-06-23T07:12:20.62699Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#assign int to weightedAssigned col:\ndata_ia[\"WeightAssigned\"] = data_ia[\"WeightAssigned\"].astype(np.float64)","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:21.48329Z","iopub.execute_input":"2023-06-23T07:12:21.483898Z","iopub.status.idle":"2023-06-23T07:12:21.490248Z","shell.execute_reply.started":"2023-06-23T07:12:21.483864Z","shell.execute_reply":"2023-06-23T07:12:21.489303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_ia.nunique()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:22.254964Z","iopub.execute_input":"2023-06-23T07:12:22.255846Z","iopub.status.idle":"2023-06-23T07:12:22.278906Z","shell.execute_reply.started":"2023-06-23T07:12:22.255796Z","shell.execute_reply":"2023-06-23T07:12:22.277898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_ia.describe()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:22.983448Z","iopub.execute_input":"2023-06-23T07:12:22.983859Z","iopub.status.idle":"2023-06-23T07:12:22.999189Z","shell.execute_reply.started":"2023-06-23T07:12:22.983811Z","shell.execute_reply":"2023-06-23T07:12:22.998444Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plot of distribution of weights assigned to every GO-term.\nplt.hist((data_ia[\"WeightAssigned\"][data_ia[\"WeightAssigned\"] != 0]),bins=\"fd\")\nplt.title(\"Distribution of Weights Assigned to different GO-terms\")\nplt.xlabel(\"log$_1$$_0$ Weights Assigned\")\nplt.ylabel(\"frequency\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:26.61946Z","iopub.execute_input":"2023-06-23T07:12:26.619885Z","iopub.status.idle":"2023-06-23T07:12:26.975684Z","shell.execute_reply.started":"2023-06-23T07:12:26.619853Z","shell.execute_reply":"2023-06-23T07:12:26.974887Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plot of distribution of log10 weights assigned to every GO-term.\nplt.hist(np.log10(data_ia[\"WeightAssigned\"][data_ia[\"WeightAssigned\"] != 0]),bins=\"fd\")\nplt.title(\"Distribution of Weights Assigned to different GO-terms\")\nplt.xlabel(\"log$_1$$_0$ Weights Assigned\")\nplt.ylabel(\"frequency\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2023-06-23T07:12:27.583454Z","iopub.execute_input":"2023-06-23T07:12:27.583856Z","iopub.status.idle":"2023-06-23T07:12:27.977199Z","shell.execute_reply.started":"2023-06-23T07:12:27.583805Z","shell.execute_reply":"2023-06-23T07:12:27.976019Z"},"trusted":true},"execution_count":null,"outputs":[]}]}