{"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":"# we will use bio python to read the sequences\nfrom Bio import SeqIO\n# numpy to work with arrays\nimport numpy as np\n# plotly to plot the graphs\nimport plotly.graph_objects as go\n# Counter to count stuff\nfrom collections import Counter","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:41.067621Z","iopub.execute_input":"2023-05-13T21:20:41.068031Z","iopub.status.idle":"2023-05-13T21:20:41.258911Z","shell.execute_reply.started":"2023-05-13T21:20:41.067998Z","shell.execute_reply":"2023-05-13T21:20:41.257843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-05-13T21:20:41.26048Z","iopub.execute_input":"2023-05-13T21:20:41.260795Z","iopub.status.idle":"2023-05-13T21:20:41.26582Z","shell.execute_reply.started":"2023-05-13T21:20:41.260769Z","shell.execute_reply":"2023-05-13T21:20:41.264788Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# read train and test fasta files\ntrain_sequences = SeqIO.parse(train_fasta, 'fasta')\ntest_sequences = SeqIO.parse(test_fasta, 'fasta')","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:41.266944Z","iopub.execute_input":"2023-05-13T21:20:41.267508Z","iopub.status.idle":"2023-05-13T21:20:41.315342Z","shell.execute_reply.started":"2023-05-13T21:20:41.267478Z","shell.execute_reply":"2023-05-13T21:20:41.314057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# take a look at the first sequence\nprint('First sequence in train fasta file:')\nprint(next(train_sequences))","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:41.317749Z","iopub.execute_input":"2023-05-13T21:20:41.318145Z","iopub.status.idle":"2023-05-13T21:20:41.328015Z","shell.execute_reply.started":"2023-05-13T21:20:41.318112Z","shell.execute_reply":"2023-05-13T21:20:41.326895Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# put ids and sequences in separate numpy arrays\ntrain_ids = np.array([seq.id for seq in SeqIO.parse(train_fasta, 'fasta')], dtype=object)\ntrain_sequences = np.array([seq.seq for seq in SeqIO.parse(train_fasta, 'fasta')], dtype=object)\ntest_ids = np.array([seq.id for seq in SeqIO.parse(test_fasta, 'fasta')], dtype=object)\ntest_sequences = np.array([seq.seq for seq in SeqIO.parse(test_fasta, 'fasta')], dtype=object)","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:41.329597Z","iopub.execute_input":"2023-05-13T21:20:41.329952Z","iopub.status.idle":"2023-05-13T21:20:52.16731Z","shell.execute_reply.started":"2023-05-13T21:20:41.329902Z","shell.execute_reply":"2023-05-13T21:20:52.166117Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# basic info: how many ids and sequences are there in train and test fasta files\n# how many unique ids and sequences are there in train and test fasta files\nprint('Train fasta file:')\nprint('Number of ids: ', len(train_ids))\nprint('Number of sequences: ', len(train_sequences))\nprint('Number of unique ids: ', len(np.unique(train_ids)))\nprint('Number of unique sequences: ', len(np.unique(train_sequences)))\nprint('Test fasta file:')\nprint('Number of ids: ', len(test_ids))\nprint('Number of sequences: ', len(test_sequences))\nprint('Number of unique ids: ', len(np.unique(test_ids)))\nprint('Number of unique sequences: ', len(np.unique(test_sequences)))","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:52.169077Z","iopub.execute_input":"2023-05-13T21:20:52.169415Z","iopub.status.idle":"2023-05-13T21:20:59.284778Z","shell.execute_reply.started":"2023-05-13T21:20:52.169385Z","shell.execute_reply":"2023-05-13T21:20:59.283698Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# it seems that there are ids which have the same sequence","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:59.287433Z","iopub.execute_input":"2023-05-13T21:20:59.287758Z","iopub.status.idle":"2023-05-13T21:20:59.292727Z","shell.execute_reply.started":"2023-05-13T21:20:59.287731Z","shell.execute_reply":"2023-05-13T21:20:59.291229Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# put only unique sequences in a numpy array\nunique_train_sequences = np.unique(train_sequences)\nunique_test_sequences = np.unique(test_sequences)\nunique_train_sequences.shape, unique_test_sequences.shape","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:20:59.294516Z","iopub.execute_input":"2023-05-13T21:20:59.297176Z","iopub.status.idle":"2023-05-13T21:21:05.773184Z","shell.execute_reply.started":"2023-05-13T21:20:59.297104Z","shell.execute_reply":"2023-05-13T21:21:05.771979Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# plot the distribution of sequence lengths in train and test fasta files\nfig = go.Figure()\nfig.add_trace(go.Histogram(x=[len(seq) for seq in train_sequences], name='Train', opacity=0.5))\nfig.add_trace(go.Histogram(x=[len(seq) for seq in test_sequences], name='Test', opacity=0.5))\nfig.update_layout(title='Distribution of sequence lengths in train and test fasta files',\n                    xaxis_title='Sequence length',\n                    yaxis_title='Count',\n                    bargap=0.2,\n                    bargroupgap=0.1)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:21:05.774557Z","iopub.execute_input":"2023-05-13T21:21:05.774876Z","iopub.status.idle":"2023-05-13T21:21:09.326869Z","shell.execute_reply.started":"2023-05-13T21:21:05.774848Z","shell.execute_reply":"2023-05-13T21:21:09.325842Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we see that the distribution rougly agrees for train and test fasta files\n# plot again but cut off at the 95th percentile to see the distribution better\nfig = go.Figure()\nfig.add_trace(go.Histogram(x=[len(seq) for seq in train_sequences], name='Train', opacity=0.5))\nfig.add_trace(go.Histogram(x=[len(seq) for seq in test_sequences], name='Test', opacity=0.5))\nfig.update_layout(title='Distribution of sequence lengths in train and test fasta files',\n                    xaxis_title='Sequence length',\n                    yaxis_title='Count',\n                    bargap=0.2,\n                    bargroupgap=0.1,\n                    xaxis_range=[0, np.percentile([len(seq) for seq in train_sequences], 95)])\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:21:09.328216Z","iopub.execute_input":"2023-05-13T21:21:09.329302Z","iopub.status.idle":"2023-05-13T21:21:12.571823Z","shell.execute_reply.started":"2023-05-13T21:21:09.329267Z","shell.execute_reply":"2023-05-13T21:21:12.570631Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# the sequence arrays contain sequences \n# convert them to strings to count the number of each amino acid in the sequences\ntrain_sequences = np.array([str(seq) for seq in train_sequences], dtype=object)\ntest_sequences = np.array([str(seq) for seq in test_sequences], dtype=object)\n# count the number of each amino acid in the sequences\ntrain_aa_counts = Counter(''.join(train_sequences))\ntest_aa_counts = Counter(''.join(test_sequences))\n# sort the amino acids by their counts\ntrain_aa_counts = {k: v for k, v in sorted(train_aa_counts.items(), key=lambda item: item[1], reverse=True)}\ntest_aa_counts = {k: v for k, v in sorted(test_aa_counts.items(), key=lambda item: item[1], reverse=True)}\n# plot the amino acid log counts\nfig = go.Figure()\nfig.add_trace(go.Bar(x=list(train_aa_counts.keys()), y=np.log(list(train_aa_counts.values())), name='Train', opacity=0.5))\nfig.add_trace(go.Bar(x=list(test_aa_counts.keys()), y=np.log(list(test_aa_counts.values())), name='Test', opacity=0.5))\nfig.update_layout(title='Log counts of amino acids in train and test fasta files',\n                    xaxis_title='Amino acid',\n                    yaxis_title='Log count',\n                    bargap=0.2,\n                    bargroupgap=0.1)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:21:12.573421Z","iopub.execute_input":"2023-05-13T21:21:12.573825Z","iopub.status.idle":"2023-05-13T21:21:21.455171Z","shell.execute_reply.started":"2023-05-13T21:21:12.573792Z","shell.execute_reply":"2023-05-13T21:21:21.454049Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# we see that the distribution rougly agrees for train and test fasta files\n# we have the following values for the amino acids:\n# use counter to get the letters\nprint('Amino acids in train fasta file: ', train_aa_counts.keys())\n# number of amino acids\nprint('Number of amino acids in train fasta file: ', len(train_aa_counts.keys()))","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:21:21.456377Z","iopub.execute_input":"2023-05-13T21:21:21.456704Z","iopub.status.idle":"2023-05-13T21:21:21.463195Z","shell.execute_reply.started":"2023-05-13T21:21:21.456676Z","shell.execute_reply":"2023-05-13T21:21:21.462079Z"},"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-05-13T21:21:21.464572Z","iopub.execute_input":"2023-05-13T21:21:21.465523Z","iopub.status.idle":"2023-05-13T21:21:21.47909Z","shell.execute_reply.started":"2023-05-13T21:21:21.465491Z","shell.execute_reply":"2023-05-13T21:21:21.477937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# redo the plot but with the amino acid names\nfig = go.Figure()\nfig.add_trace(go.Bar(x=[amino_acid_dict[aa] for aa in train_aa_counts.keys()], y=np.log(list(train_aa_counts.values())), name='Train', opacity=0.5))\nfig.add_trace(go.Bar(x=[amino_acid_dict[aa] for aa in test_aa_counts.keys()], y=np.log(list(test_aa_counts.values())), name='Test', opacity=0.5))\nfig.update_layout(title='Log counts of amino acids in train and test fasta files',\n                    xaxis_title='Amino acid',\n                    yaxis_title='Log count',\n                    bargap=0.2,\n                    bargroupgap=0.1)\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2023-05-13T21:21:21.482562Z","iopub.execute_input":"2023-05-13T21:21:21.482985Z","iopub.status.idle":"2023-05-13T21:21:21.513747Z","shell.execute_reply.started":"2023-05-13T21:21:21.482944Z","shell.execute_reply":"2023-05-13T21:21:21.512496Z"},"trusted":true},"execution_count":null,"outputs":[]}]}