{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# FASTA files EDA:\n\nThe goal of this notebook is to do some preliminary EDA just on the sequences and discuss the FASTA files in which they are stored\n\n## Takeaways:\n1. The train sequences are significantly longer than the test sequences.\n2. There is no significant difference in the amino acid compositions between the train and test sequences\n\n## Next steps: \n1. Explore GO Terms and Taxonomy data\n2. Analyze embeddings \n3. Start modeling strategy\n\nCredit to Simon Veitner for some of the figures: https://www.kaggle.com/code/simonveitner/basic-eda-of-sequences\n\n## What is a FASTA file?\n\nA FASTA file is a simple and commonly used file format in bioinformatics for storing biological sequence data, such as DNA or protein sequences. It's like a text document that organizes genetic information in a specific way.\n\nImagine you have a long, complex code that represents a sequence of DNA or protein molecules. It's like a secret message that only scientists can understand. The FASTA file provides a standardized way to write down and share  a sequence of residues with other researchers.\n\nIn 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\nTo make it more readable, the code is often split into lines of a certain length, called line wrapping. This helps scientists visually process the information and ensures that the file doesn't become too wide or difficult to work with.\n\nSo, in simpler terms, a FASTA file is like a special document that stores genetic information. It's a way for scientists to share and organize sequences of DNA or proteins, making it easier for them to study and understand the building blocks of life.\n\nWe'll use Biopython.SeqIO to work with the FASTA files.","metadata":{"execution":{"iopub.status.busy":"2023-05-23T19:00:50.883531Z","iopub.execute_input":"2023-05-23T19:00:50.883917Z","iopub.status.idle":"2023-05-23T19:00:50.897651Z","shell.execute_reply.started":"2023-05-23T19:00:50.883888Z","shell.execute_reply":"2023-05-23T19:00:50.89591Z"}}},{"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\nimport pandas as pd\nimport matplotlib.pyplot as plt\nfrom scipy import stats\nimport re","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2023-05-26T19:39:37.624611Z","iopub.execute_input":"2023-05-26T19:39:37.625007Z","iopub.status.idle":"2023-05-26T19:39:38.136159Z","shell.execute_reply.started":"2023-05-26T19:39:37.624976Z","shell.execute_reply":"2023-05-26T19:39:38.134939Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Think of Bio.SeqIO as a toolbox that provides various tools to handle biological sequences, such as DNA, RNA, or protein sequences. It helps you read, write, and manipulate these sequences in different file formats.\n\nImagine you have a bookshelf filled with books, each representing a different biological sequence. The Bio.SeqIO toolbox allows you to pick up these books, read the sequences written inside, and perform different operations on them.\n\nFor example, we can use Bio.SeqIO to:\n\n- Read sequences from files in formats like FASTA or GenBank.\n- Write sequences to files in different formats.\n- Extract specific information from sequences, like their length or identification.\n- Manipulate sequences by combining, splitting, or altering them.\n- Perform calculations or analyses on sequences, such as sequence alignment or translation.","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-05-26T19:39:43.623682Z","iopub.execute_input":"2023-05-26T19:39:43.624489Z","iopub.status.idle":"2023-05-26T19:39:43.630003Z","shell.execute_reply.started":"2023-05-26T19:39:43.624449Z","shell.execute_reply":"2023-05-26T19:39:43.62879Z"},"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-26T19:40:45.554005Z","iopub.execute_input":"2023-05-26T19:40:45.554431Z","iopub.status.idle":"2023-05-26T19:40:45.59186Z","shell.execute_reply.started":"2023-05-26T19:40:45.554374Z","shell.execute_reply":"2023-05-26T19:40:45.590555Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**SeqIO.parse(train_fasta, 'fasta')** reads the FASTA file specified by the variable train_fasta and returns an iterator that can be used to iterate over the sequences in the file.","metadata":{}},{"cell_type":"code","source":"help(SeqIO.parse)","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:27:19.384292Z","iopub.execute_input":"2023-05-26T19:27:19.384656Z","iopub.status.idle":"2023-05-26T19:27:19.390792Z","shell.execute_reply.started":"2023-05-26T19:27:19.384629Z","shell.execute_reply":"2023-05-26T19:27:19.389472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We use the **SeqIO.parse()** function to read in the FASTA file.\n\nHere's what the different parts of the above explanation mean:\n\n- parse(handle, format, alphabet=None): This is the function signature that tells us how to use the SeqIO.parse() function. It takes three arguments:\n\n    1. **handle**: It refers to the file or filename that we want to read. It can be a file handle or a string specifying the filename.\n    2. **format**: It specifies the format of the file we're reading, such as \"fasta\" for FASTA format.\n    3. **alphabet**: This argument is no longer used, so we ignore it by setting it to None.\n- The explanation shows an example of how to use SeqIO.parse() to read a FASTA file and loop over the sequence records.\n    1. It imports **SeqIO** from Biopython and specifies the filename and format.\n    2. By using a for loop with **SeqIO.parse()**, it iterates over each record in the file.\n    3. Inside the loop, it demonstrates printing the ID and the length of each sequence.\n- It also provides an example of how to use **SeqIO.parse()** with a string of sequence data instead of a file.\n    1. It shows how to convert the string data into a file-like handle using **StringIO**.\n    2. The example then demonstrates parsing the data and printing the ID and sequence of each record.","metadata":{}},{"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-26T19:27:21.59016Z","iopub.execute_input":"2023-05-26T19:27:21.59051Z","iopub.status.idle":"2023-05-26T19:27:21.600616Z","shell.execute_reply.started":"2023-05-26T19:27:21.590483Z","shell.execute_reply":"2023-05-26T19:27:21.599382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### What is going on in this output?\n\n- The output begins with \"ID: P20536\", which represents the unique identifier of the sequence.\n- The line \"Name: P20536\" suggests that \"P20536\" might be the name associated with the sequence.\n- The \"Description\" field provides additional details about the sequence, including its source (Vaccinia virus), strain (Copenhagen), organism, and other information.\n    - OX is the taxon ID\n    - GN is the gene name\n    - PE is the protein existance level. In this case, \"1\" indicates that experimental evidence exists for the protein encoded by the sequence (Protein existence levels are used to categorize the reliability and confidence of protein predictions).\n    - SV denotes the sequence version. In this casem \"1\" indicates this is the first version of the sequence\n- \"Number of features: 0\" indicates that there are no additional features or annotations associated with this particular sequence.\n- Finally, \"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\n### Lets mess around with these data some more:","metadata":{}},{"cell_type":"code","source":"# retrieve the next sequence record from the train_sequences iterator\nrecord = next(train_sequences)\n# Get all attributes of the record\nrecord_attributes = record.__dict__\n\n# Print all attributes and their values\nfor attribute, value in record_attributes.items():\n    print(attribute, \":\", value)\n    print(type(attribute))","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:27:22.759265Z","iopub.execute_input":"2023-05-26T19:27:22.759622Z","iopub.status.idle":"2023-05-26T19:27:22.76676Z","shell.execute_reply.started":"2023-05-26T19:27:22.759594Z","shell.execute_reply":"2023-05-26T19:27:22.765493Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's examine the output:\n\n- _seq represents the sequence itself, with the attribute name _seq and the sequence string as its value.\n- id represents the ID attribute of the sequence record, with the attribute name id and the ID value.\n- name represents the name attribute of the sequence record, with the attribute name name and the name value. In this case, the ID and name values are the same.\n- description represents the description attribute of the sequence record, with the attribute name description and the description value containing additional information about the sequence. It includes details like the ID, species, gene name, and more.\n- dbxrefs represents the dbxrefs attribute of the sequence record, which refers to any database cross-references associated with the sequence. In this case, the attribute value is an empty list [], indicating no cross-references are present.\n- annotations represents the annotations attribute of the sequence record, which contains any additional annotations or metadata associated with the sequence. In this case, the attribute value is an empty dictionary {}, indicating no annotations are present.\n- _per_letter_annotations line represents the _per_letter_annotations attribute of the sequence record, which allows annotations to be associated with each letter of the sequence individually. In this case, the attribute value is an empty dictionary {}, indicating no per-letter annotations are present.\n- features represents the features attribute of the sequence record, which can store any features or regions of interest present in the sequence. In this case, the attribute value is an empty list [], indicating no features are present.","metadata":{}},{"cell_type":"code","source":"sequence_count = sum(1 for _ in train_sequences)\ntrain_sequences = SeqIO.parse(train_fasta, 'fasta')\nsequence_count","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:27:26.099815Z","iopub.execute_input":"2023-05-26T19:27:26.100209Z","iopub.status.idle":"2023-05-26T19:27:28.189376Z","shell.execute_reply.started":"2023-05-26T19:27:26.100185Z","shell.execute_reply":"2023-05-26T19:27:28.188145Z"},"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-26T19:39:52.197603Z","iopub.execute_input":"2023-05-26T19:39:52.19809Z","iopub.status.idle":"2023-05-26T19:40:02.375203Z","shell.execute_reply.started":"2023-05-26T19:39:52.198053Z","shell.execute_reply":"2023-05-26T19:40:02.374242Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# EDA:\n## Basic Statistics","metadata":{}},{"cell_type":"markdown","source":"### Number of sequences","metadata":{}},{"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-26T19:27:34.508237Z","iopub.execute_input":"2023-05-26T19:27:34.508556Z","iopub.status.idle":"2023-05-26T19:27:38.927049Z","shell.execute_reply.started":"2023-05-26T19:27:34.508529Z","shell.execute_reply":"2023-05-26T19:27:38.925995Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Length Distribution of sequences","metadata":{}},{"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":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-26T19:27:41.253805Z","iopub.execute_input":"2023-05-26T19:27:41.254269Z","iopub.status.idle":"2023-05-26T19:27:43.343077Z","shell.execute_reply.started":"2023-05-26T19:27:41.254233Z","shell.execute_reply":"2023-05-26T19:27:43.342128Z"},"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\n# Recreate the iterators\ntrain_sequences = SeqIO.parse(train_fasta, 'fasta')\ntest_sequences = SeqIO.parse(test_fasta, 'fasta')\n\n# Calculate sequence lengths\ntrain_sequence_lengths = [len(seq) for seq in train_sequences]\ntest_sequence_lengths = [len(seq) for seq in test_sequences]\n\n# Create the histogram figure\nfig = go.Figure()\nfig.add_trace(go.Histogram(x=train_sequence_lengths, name='Train', opacity=0.5))\nfig.add_trace(go.Histogram(x=test_sequence_lengths, name='Test', opacity=0.5))\nfig.update_layout(\n    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(train_sequence_lengths, 95)]\n)\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-26T19:27:48.25977Z","iopub.execute_input":"2023-05-26T19:27:48.260217Z","iopub.status.idle":"2023-05-26T19:27:52.474471Z","shell.execute_reply.started":"2023-05-26T19:27:48.26018Z","shell.execute_reply":"2023-05-26T19:27:52.473098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(train_sequence_lengths).describe()","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:27:53.133573Z","iopub.execute_input":"2023-05-26T19:27:53.133923Z","iopub.status.idle":"2023-05-26T19:27:53.218862Z","shell.execute_reply.started":"2023-05-26T19:27:53.133896Z","shell.execute_reply":"2023-05-26T19:27:53.217351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pd.DataFrame(test_sequence_lengths).describe()","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:27:58.498309Z","iopub.execute_input":"2023-05-26T19:27:58.498803Z","iopub.status.idle":"2023-05-26T19:27:58.545947Z","shell.execute_reply.started":"2023-05-26T19:27:58.498763Z","shell.execute_reply":"2023-05-26T19:27:58.544752Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Lets perform a Kolmogorov-Smirnov test to determine if there are significant differences between the length distributions. \n\nThe KS test is a non-parametric test that does not assume any specific distribution for the data. It is commonly used to compare the distributions of two samples, in this case, the sequence lengths from the train and test files. The KS test calculates the maximum difference between the cumulative distribution functions (CDFs) of the two samples and provides a p-value indicating the likelihood of observing such a difference by chance.","metadata":{}},{"cell_type":"code","source":"# Perform Kolmogorov-Smirnov test\nks_statistic, p_value = stats.ks_2samp(train_sequence_lengths, test_sequence_lengths)\nprint(\"Kolmogorov-Smirnov test:\")\nprint(\"KS statistic:\", ks_statistic)\nprint(\"p-value:\", p_value)","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:28:02.642144Z","iopub.execute_input":"2023-05-26T19:28:02.642511Z","iopub.status.idle":"2023-05-26T19:28:02.752737Z","shell.execute_reply.started":"2023-05-26T19:28:02.642487Z","shell.execute_reply":"2023-05-26T19:28:02.751701Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The Kolmogorov-Smirnov test results indicate the following:\n\nKS statistic: The KS statistic is a measure of the maximum difference between the cumulative distribution functions (CDFs) of the two sample distributions. In this case, **the KS statistic is 0.0597**, which represents the maximum discrepancy between the cumulative distribution of sequence lengths in the train and test files.\n\np-value: The p-value is a measure of the evidence against the null hypothesis. In this case, **the p-value is very close to zero** (approximately 1.28e-220), which suggests extremely strong evidence against the null hypothesis. A small p-value indicates that **the two sample distributions significantly differ from each other in terms of sequence length distribution.**\n\nBased on the results, **there is a statistically significant difference between the distribution of sequence lengths in the train and test files.** The small p-value indicates that the observed difference is unlikely to be due to random chance alone. Therefore, you can conclude that the train and test files have significantly different distributions of sequence lengths.\n\nLets also perform a Mann-Whitney U test: This non-parametric test assesses whether the distributions of two independent samples differ significantly. It is commonly used when comparing two groups of data to determine if one group tends to have higher or lower values than the other.","metadata":{}},{"cell_type":"code","source":"# Perform the Mann-Whitney U test\nstatistic, p_value = stats.mannwhitneyu(train_sequence_lengths, test_sequence_lengths)\nprint(\"Mann-Whitney U test:\")\nprint(\"U statistic:\", statistic)\nprint(\"p-value:\", p_value)","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:28:06.066232Z","iopub.execute_input":"2023-05-26T19:28:06.066587Z","iopub.status.idle":"2023-05-26T19:28:06.340727Z","shell.execute_reply.started":"2023-05-26T19:28:06.066558Z","shell.execute_reply":"2023-05-26T19:28:06.339646Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The Mann-Whitney U test results indicate the following:\n\nU statistic: The U statistic is a measure of the rank-sum of the two samples. In this case, the U statistic is 10900212030.5, which represents the sum of the ranks of the values in the smaller sample (train_lengths) when combined with the ranks of the values in the larger sample (test_lengths).\n\np-value: The p-value is a measure of the evidence against the null hypothesis. In this case, the p-value is approximately 7.60e-301, which is extremely close to zero. This suggests strong evidence against the null hypothesis and indicates a significant difference between the sequence length distributions of the train and test datasets.\n\nThis corroborates the KS test results.","metadata":{}},{"cell_type":"markdown","source":"## Sequence Composition Analysis","metadata":{}},{"cell_type":"markdown","source":"### Frequency of each residue","metadata":{}},{"cell_type":"code","source":"# Parse the train and test FASTA files\ntrain_sequences = SeqIO.parse(train_fasta, \"fasta\")\ntest_sequences = SeqIO.parse(test_fasta, \"fasta\")\n\n# Calculate the frequency of each amino acid in train sequences\ntrain_frequency = Counter()\nfor record in train_sequences:\n    sequence = str(record.seq)\n    train_frequency.update(sequence)\n\n# Calculate the frequency of each amino acid in test sequences\ntest_frequency = Counter()\nfor record in test_sequences:\n    sequence = str(record.seq)\n    test_frequency.update(sequence)\n\n# Print the frequency of each amino acid in train sequences\nprint(\"Frequency of amino acids in train sequences:\")\nfor item in train_frequency.most_common():\n    print(item[0], \":\", item[1])\n\n# Print the frequency of each amino acid in test sequences\nprint(\"\\nFrequency of amino acids in test sequences:\")\nfor item in test_frequency.most_common():\n    print(item[0], \":\", item[1])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-26T19:28:11.850486Z","iopub.execute_input":"2023-05-26T19:28:11.851271Z","iopub.status.idle":"2023-05-26T19:28:21.725585Z","shell.execute_reply.started":"2023-05-26T19:28:11.851235Z","shell.execute_reply":"2023-05-26T19:28:21.723954Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Frequencies of amino acids in train sequences\ntrain_freq = {\n    'L': 7386109, 'S': 6493284, 'A': 5601393, 'E': 5388873, 'G': 5061691,\n    'V': 4931367, 'K': 4594229, 'P': 4382337, 'T': 4378528, 'R': 4283020,\n    'D': 4151113, 'I': 3879567, 'Q': 3563728, 'N': 3400751, 'F': 2930638,\n    'Y': 2239153, 'H': 1914456, 'M': 1776672, 'C': 1504314, 'W': 888449,\n    'X': 2771, 'U': 148, 'O': 4, 'B': 4, 'Z': 4\n}\n\n# Frequencies of amino acids in test sequences\ntest_freq = {\n    'L': 6571441, 'S': 5512318, 'A': 4750271, 'E': 4588833, 'G': 4318376,\n    'V': 4223179, 'K': 4060313, 'T': 3650414, 'P': 3648374, 'R': 3611308,\n    'I': 3501751, 'D': 3501462, 'Q': 2950714, 'N': 2919742, 'F': 2667898,\n    'Y': 1970413, 'H': 1628870, 'M': 1545971, 'C': 1264828, 'W': 798942,\n    'X': 657, 'U': 218, 'Z': 45, 'B': 39\n}\n\n# Determine all unique amino acids\nall_amino_acids = list(set(list(train_freq.keys()) + list(test_freq.keys())))\n\n# Add zero frequencies for missing amino acids\nfor aa in all_amino_acids:\n    if aa not in train_freq:\n        train_freq[aa] = 0\n    if aa not in test_freq:\n        test_freq[aa] = 0\n\n# Calculate relative frequencies by dividing by the total number of sequences\ntrain_total = sum(train_freq.values())\ntest_total = sum(test_freq.values())\n\ntrain_relative_freq = {aa: freq / train_total for aa, freq in train_freq.items()}\ntest_relative_freq = {aa: freq / test_total for aa, freq in test_freq.items()}\n\n# Visualize the distributions using bar plots\nlabels = train_relative_freq.keys()\nx = np.arange(len(labels))\nwidth = 0.35\n\nfig, ax = plt.subplots()\nrects1 = ax.bar(x - width/2, train_relative_freq.values(), width, label='Train')\nrects2 = ax.bar(x + width/2, test_relative_freq.values(), width, label='Test')\n\nax.set_ylabel('Relative Frequency')\nax.set_xlabel('Amino Acid')\nax.set_title('Comparison of Amino Acid Frequencies')\nax.set_xticks(x)\nax.set_xticklabels(labels)\nax.legend()\n\nplt.xticks(rotation=90)\nplt.tight_layout()\nplt.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2023-05-26T19:28:21.727983Z","iopub.execute_input":"2023-05-26T19:28:21.728685Z","iopub.status.idle":"2023-05-26T19:28:22.131549Z","shell.execute_reply.started":"2023-05-26T19:28:21.728656Z","shell.execute_reply":"2023-05-26T19:28:22.130215Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from scipy.stats import chi2_contingency\n# Perform chi-square test\nobserved = np.array(list(train_relative_freq.values()))\nexpected = np.array(list(test_relative_freq.values()))\n\nchi2, p_value, _, _ = chi2_contingency([observed, expected])\n\nprint(\"Chi-square test results:\")\nprint(\"p-value:\", p_value)","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:28:52.25918Z","iopub.execute_input":"2023-05-26T19:28:52.259521Z","iopub.status.idle":"2023-05-26T19:28:52.267365Z","shell.execute_reply.started":"2023-05-26T19:28:52.259497Z","shell.execute_reply":"2023-05-26T19:28:52.266434Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The chi-square test results indicate that there is no significant difference in the amino acid frequencies between the train and test sequences. The p-value of 1.0 suggests that the observed differences in frequencies can be attributed to chance, and there is no evidence to reject the null hypothesis of equal distributions.\n\nIn other words, based on the comparison of amino acid frequencies, there is no notable difference between the train and test sequences in terms of amino acid composition.","metadata":{}},{"cell_type":"markdown","source":"#### Biases and patterns\n\nTo identify biases or patterns in the frequencies of different amino acids, we can compare the relative frequencies of amino acids between the train and test sequences. One approach is to calculate the difference in relative frequencies and analyze which amino acids show the largest deviations.\n","metadata":{}},{"cell_type":"code","source":"# Define the list of amino acids\namino_acids = ['A', 'C', 'D', 'E', 'F', 'G', 'H', 'I', 'K', 'L', 'M', 'N', 'P', 'Q', 'R', 'S', 'T', 'V', 'W', 'Y']\n\n# Calculate the difference in relative frequencies\nfrequency_difference = {\n    amino_acid: abs(train_relative_freq[amino_acid] - test_relative_freq[amino_acid])\n    for amino_acid in amino_acids\n}\n\n# Sort the amino acids based on the difference in frequencies\nsorted_amino_acids = sorted(frequency_difference.keys(), key=lambda x: frequency_difference[x], reverse=True)\n\n# Print the amino acids with the largest deviations\nnum_top_amino_acids = 5  # Number of top amino acids to display\nprint(f\"Top {num_top_amino_acids} amino acids with the largest deviations:\")\nfor amino_acid in sorted_amino_acids[:num_top_amino_acids]:\n    print(f\"Amino Acid: {amino_acid}\")\n    print(f\"Train Frequency: {train_relative_freq[amino_acid]:.4f}\")\n    print(f\"Test Frequency: {test_relative_freq[amino_acid]:.4f}\")\n    print(f\"Difference: {frequency_difference[amino_acid]:.4f}\")\n    print()","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:28:58.429662Z","iopub.execute_input":"2023-05-26T19:28:58.430028Z","iopub.status.idle":"2023-05-26T19:28:58.438345Z","shell.execute_reply.started":"2023-05-26T19:28:58.430001Z","shell.execute_reply":"2023-05-26T19:28:58.436786Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"To compare the observed frequencies of amino acids to the provided background distribution, we can calculate the expected frequencies based on a background distribution and perform a statistical test to assess the deviation.","metadata":{}},{"cell_type":"code","source":"# Observed frequencies of amino acids (based on previous calculations)\nobserved_frequencies = {\n    'A': 0.0938,\n    'Q': 0.0517,\n    'L': 0.0971,\n    'S': 0.0394,\n    'R': 0.0556,\n    'E': 0.0394,\n    'K': 0.0556,\n    'T': 0.0539,\n    'N': 0.0517,\n    'G': 0.0556,\n    'M': 0.0372,\n    'W': 0.0372,\n    'D': 0.0556,\n    'H': 0.0372,\n    'F': 0.0394,\n    'Y': 0.0372,\n    'C': 0.0372,\n    'I': 0.0556,\n    'P': 0.0539,\n    'V': 0.0556,\n}\n\n# Background distribution frequencies (based on https://www.ebi.ac.uk/uniprot/TrEMBLstats)\nbackground_frequencies = {\n    'A': 9.06,\n    'Q': 3.79,\n    'L': 9.87,\n    'S': 6.78,\n    'R': 5.84,\n    'E': 6.24,\n    'K': 4.92,\n    'T': 5.54,\n    'N': 3.79,\n    'G': 7.29,\n    'M': 2.33,\n    'W': 1.30,\n    'D': 5.47,\n    'H': 2.21,\n    'F': 3.89,\n    'Y': 2.88,\n    'C': 1.28,\n    'I': 5.55,\n    'P': 4.97,\n    'V': 6.88,\n}\n\n# Calculate expected frequencies based on the background distribution\ntotal_sequences = len(train_sequence_lengths) + len(train_sequence_lengths)\nexpected_frequencies = {\n    amino_acid: background_frequencies[amino_acid] / 100 * total_sequences\n    for amino_acid in observed_frequencies\n}\n\n# Perform statistical test to assess the deviation\nchi2, p_value, _, _ = chi2_contingency([list(observed_frequencies.values()), list(expected_frequencies.values())])\n\nprint(\"Chi-square test results:\")\nprint(\"p-value:\", p_value)","metadata":{"execution":{"iopub.status.busy":"2023-05-26T19:29:06.423794Z","iopub.execute_input":"2023-05-26T19:29:06.42418Z","iopub.status.idle":"2023-05-26T19:29:06.433802Z","shell.execute_reply.started":"2023-05-26T19:29:06.42415Z","shell.execute_reply":"2023-05-26T19:29:06.432869Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In this code, the observed frequencies are compared to the expected frequencies based on the provided background distribution using the chi-square test. The resulting p-value of 1.0 indicates that there is no significant deviation between the observed frequencies of amino acids and the expected frequencies based on the provided background distribution. This suggests that the observed frequencies are consistent with the background distribution, and there are no notable biases or patterns in the composition of amino acids compared to the expected distribution.\n\nA p-value of 1.0 suggests that the observed frequencies are likely due to random chance and not indicative of a significant difference from the expected distribution. Therefore, based on this analysis, it appears that the frequencies of amino acids in our sequences are in line with the background distribution.","metadata":{}}]}