{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.10.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":41875,"databundleVersionId":5521661,"sourceType":"competition"},{"sourceId":7624238,"sourceType":"datasetVersion","datasetId":4441442}],"dockerImageVersionId":30664,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Goals\n\n* To investigate and understand the distribution of amino acid sequences (n-grams) within a dataset of protein sequences. This includes identifying prevalent patterns and understanding how these patterns distribute across different lengths of n-grams.\n* To analyze how these n-gram patterns fit within known statistical distributions, such as Benford's Law, the Pareto Principle, and Zipf's Law. This aim seeks to apply mathematical models to biological data to uncover any inherent order or regularity in the way protein sequences are structured.\n* To determine the n-gram sizes that reveal the most informative and statistically significant patterns. This involves comparing the adherence of different n-gram sizes to the expected distributions and identifying the \"sweet spot\" where patterns become most meaningful or pronounced.\n* Through the lens of n-gram analysis, to gain insights into the complexity and variability of protein sequences. This could help in understanding the level of sequence conservation, the prevalence of certain motifs, and the diversity of amino acid combinations.\n* To refine and validate n-gram analysis as a methodological approach in bioinformatics for studying protein sequences. This includes assessing the utility and limitations of this approach for uncovering meaningful biological patterns.","metadata":{}},{"cell_type":"markdown","source":"# Conclusions\n\n**Benford's Law Approximation:**\n\n* The observed approximation of the Benford's Law curve with increasing n-gram size up to 5-grams suggests that the distribution of leading digits in frequency counts becomes more 'natural' or 'expected' as you look at longer sequences. This could imply that the larger motifs captured by tetra- and pentagrams are more representative of the underlying biological processes that govern protein sequences.\n* The decline in approximation quality beyond 5-grams may indicate that as sequences become longer, the specific frequency of occurrence of each distinct n-gram becomes more uniform or less patterned in a way that diverges from Benford's expected distribution. This could be due to a greater variety of longer motifs that are less universally conserved or functionally significant across different proteins.\n* The lack of fit for 1-grams is not surprising, as individual amino acids may not be expected to follow Benford's Law, which is more commonly associated with sets of numbers that span multiple orders of magnitude.\n\n**Pareto Distribution:**\n\n* The fluctuations in the intersection point between the data curve and the Pareto baseline for n-grams up to 6-grams suggest that certain lengths of motifs are more critical or prevalent in the data set. Specifically, the intersection at X=36 for pentagrams indicates a significant concentration of occurrences among a relatively small proportion of these 5-grams, aligning well with the Pareto principle.\n* The sharp increase in the intersection point for 7-grams and higher implies that these longer n-grams are more evenly distributed in frequency. This could suggest that the specific combinations of seven or more amino acids are less conserved or less essential for the function of the proteins in your data set.\n* The intersection at X=50 for 1-grams again points to a more uniform distribution of individual amino acid occurrences compared to longer motifs.\n\n**Zipf's Law:**\n\n* The monotonic improvement in the approximation of Zipf's Law with increasing n-gram size could indicate that the frequencies of n-grams become more predictable and structured as you consider longer sequences. This might suggest that longer motifs have a more consistent and significant role in the structure and function of proteins, while individual amino acids or shorter motifs are more variable in their distribution.","metadata":{}},{"cell_type":"markdown","source":"# Further research\n\n* **Biological Significance**: The observed patterns suggest that motifs of length four and five (tetragrams and pentagrams) are particularly significant in the protein sequences analyzed. These may represent functional domains or structural elements that are conserved across different proteins.\n\n* **Sequence Diversity**: There is greater diversity and less predictability in the occurrence of longer sequences beyond pentagrams, as indicated by the deviation from Benford's and Pareto distributions.\n\n* **Functional Prediction**: The results could be used to prioritize certain n-grams for functional prediction, with tetragrams and pentagrams being particularly promising targets for further investigation.\n\n* **Data Characterization**: The character of the data suggests a complex landscape where certain motifs are highly significant and potentially functionally relevant, while there is also a rich diversity of less common sequences.\n\n* **Modeling Implications**: For predictive modeling, focusing on tetragrams and pentagrams may yield the most informative features, whereas models that include longer n-grams may need to contend with a larger number of less frequent, more variable features.","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nfrom collections import defaultdict\nfrom Bio import SeqIO\nimport random\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2024-03-05T07:48:23.868387Z","iopub.execute_input":"2024-03-05T07:48:23.868704Z","iopub.status.idle":"2024-03-05T07:48:25.162789Z","shell.execute_reply.started":"2024-03-05T07:48:23.868682Z","shell.execute_reply":"2024-03-05T07:48:25.161258Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Parse the sequences for analysis\nfn = '/kaggle/input/cafa-5-protein-function-prediction/Train/train_sequences.fasta'\nsequences = [str(seq_record.seq).upper() for seq_record in SeqIO.parse(fn, \"fasta\")]","metadata":{"execution":{"iopub.status.busy":"2024-03-05T07:48:27.931755Z","iopub.execute_input":"2024-03-05T07:48:27.932137Z","iopub.status.idle":"2024-03-05T07:48:30.768363Z","shell.execute_reply.started":"2024-03-05T07:48:27.932109Z","shell.execute_reply":"2024-03-05T07:48:30.767118Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to count n-grams\ndef count_ngrams(sequences, n):\n    ngram_counts = defaultdict(int)\n    for seq in sequences:\n        for i in range(len(seq) - n + 1):\n            ngram = seq[i:i+n]\n            ngram_counts[ngram] += 1\n    return ngram_counts","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:47:29.446802Z","iopub.execute_input":"2024-03-04T10:47:29.447239Z","iopub.status.idle":"2024-03-04T10:47:29.460514Z","shell.execute_reply.started":"2024-03-04T10:47:29.44718Z","shell.execute_reply":"2024-03-04T10:47:29.459051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Function to plot n-gram frequencies\ndef plot_ngram_frequencies(ngram_counts, title):\n    sorted_ngrams = dict(sorted(ngram_counts.items(), key=lambda item: item[1], reverse=True))\n    ngrams = list(sorted_ngrams.keys())\n    counts = list(sorted_ngrams.values())\n\n    plt.figure(figsize=(20, 10))\n    plt.bar(ngrams[:100], counts[:100])  # Displaying top 100 for readability\n    plt.xlabel(f'{title}')\n    plt.ylabel('Frequency')\n    plt.title(f'{title} Frequency Distribution in Protein Sequences')\n    plt.xticks(rotation=90)\n    plt.show() ","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:47:45.255416Z","iopub.execute_input":"2024-03-04T10:47:45.255836Z","iopub.status.idle":"2024-03-04T10:47:45.265196Z","shell.execute_reply.started":"2024-03-04T10:47:45.2558Z","shell.execute_reply":"2024-03-04T10:47:45.263784Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def analyze_ngram_benford(ngram_counts, title):\n    \"\"\"\n    Analyzes and plots the leading digit distribution of n-gram frequencies\n    against the expected distribution under Benford's Law.\n    \n    Parameters:\n    - ngram_counts: A dictionary with n-grams as keys and their frequencies as values.\n    - title: A string for the plot title to indicate the n-gram size being analyzed.\n    \"\"\"\n    # Calculate leading digits\n    leading_digits = [int(str(count)[0]) for count in ngram_counts.values()]\n    leading_digit_distribution = defaultdict(int)\n    for digit in leading_digits:\n        leading_digit_distribution[digit] += 1\n    \n    # Normalize the distribution to percentages\n    total_counts = sum(leading_digit_distribution.values())\n    for digit in leading_digit_distribution:\n        leading_digit_distribution[digit] = (leading_digit_distribution[digit] / total_counts) * 100\n    \n    # Expected distribution under Benford's Law\n    benford_distribution = {digit: np.log10(1 + 1/digit) * 100 for digit in range(1, 10)}\n    \n    # Plotting\n    plt.figure(figsize=(10, 6))\n    plt.bar(leading_digit_distribution.keys(), leading_digit_distribution.values(), label='Observed', width=0.4, align='center')\n    plt.plot(benford_distribution.keys(), benford_distribution.values(), color='red', label='Benford\\'s Law', linestyle='-', marker='o')\n    plt.xlabel('Leading Digit')\n    plt.ylabel('Percentage')\n    plt.xticks(range(1, 10))\n    plt.legend()\n    plt.title(f'Benford\\'s Law Comparison for {title}')\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:47:48.886625Z","iopub.execute_input":"2024-03-04T10:47:48.88701Z","iopub.status.idle":"2024-03-04T10:47:48.899878Z","shell.execute_reply.started":"2024-03-04T10:47:48.886979Z","shell.execute_reply":"2024-03-04T10:47:48.898811Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_pareto_distribution(ngram_counts, title):\n    # Convert ngram_counts values to a sorted list in descending order\n    frequencies = sorted(ngram_counts.values(), reverse=True)\n    total = sum(frequencies)\n    cumulative_sum = np.cumsum(frequencies)\n    \n    # Convert cumulative sums to percentages of the total sum\n    cumulative_percentage = 100 * cumulative_sum / total\n    \n    # Calculate the percentage of n-grams for the x-axis\n    ngrams_percentage = 100 * np.arange(1, len(frequencies) + 1) / len(frequencies)\n    \n    # Plotting\n    plt.figure(figsize=(10, 6))\n    plt.plot(ngrams_percentage, cumulative_percentage, marker='o', linestyle='-', color='b')\n    plt.xlabel(f'Percentage of {title}')\n    plt.ylabel('Cumulative Percentage of Frequencies')\n    plt.title(f'Pareto Distribution of {title} Frequencies')\n    plt.grid(True)\n    \n    # Adding a line at 80% for reference\n    plt.axhline(80, color='red', linestyle='--')\n    plt.text(100, 80, '80%', color = 'red')\n    \n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:47:52.355845Z","iopub.execute_input":"2024-03-04T10:47:52.356336Z","iopub.status.idle":"2024-03-04T10:47:52.366957Z","shell.execute_reply.started":"2024-03-04T10:47:52.356302Z","shell.execute_reply":"2024-03-04T10:47:52.365479Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def test_zipfs_law(ngram_counts, title):\n    # Sort n-grams by frequency in descending order\n    sorted_ngrams = sorted(ngram_counts.items(), key=lambda item: item[1], reverse=True)\n    frequencies = [count for ngram, count in sorted_ngrams]\n    \n    # Generate ranks for n-grams\n    ranks = range(1, len(frequencies) + 1)\n    \n    # Plotting on a log-log scale\n    plt.figure(figsize=(10, 6))\n    plt.loglog(ranks, frequencies, marker=\"o\", linestyle='None')\n    \n    # Plot a reference line for comparison\n    max_freq = max(frequencies)\n    plt.loglog(ranks, [max_freq/rank for rank in ranks], label=\"Reference Line for Zipf's Law\", linestyle='--')\n    \n    plt.xlabel(f'Rank of {title}')\n    plt.ylabel(f'Frequency of {title}')\n    plt.title(f'Testing Zipf\\'s Law for {title} Frequencies')\n    plt.legend()\n    plt.grid(True)\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:47:55.58646Z","iopub.execute_input":"2024-03-04T10:47:55.5869Z","iopub.status.idle":"2024-03-04T10:47:55.597271Z","shell.execute_reply.started":"2024-03-04T10:47:55.586862Z","shell.execute_reply":"2024-03-04T10:47:55.595924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Analyze and visualize for 2-grams, 3-grams, 4-grams\nfor n in [2, 3, 4]:\n    ngram_counts = count_ngrams(sequences, n)\n    plot_ngram_frequencies(ngram_counts, f'{n}-gram')\n    analyze_ngram_benford(ngram_counts, f'{n}-gram')\n    plot_pareto_distribution(ngram_counts, f'{n}-gram')\n    test_zipfs_law(ngram_counts, f'{n}-gram')","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:47:59.774728Z","iopub.execute_input":"2024-03-04T10:47:59.775137Z","iopub.status.idle":"2024-03-04T10:50:31.774542Z","shell.execute_reply.started":"2024-03-04T10:47:59.775106Z","shell.execute_reply":"2024-03-04T10:50:31.773305Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Analyze and visualize for 5-grams, 6-grams, 7-grams\nfor n in [5, 6, 7]:\n    ngram_counts = count_ngrams(sequences, n)\n    plot_ngram_frequencies(ngram_counts, f'{n}-gram')\n    analyze_ngram_benford(ngram_counts, f'{n}-gram')\n    plot_pareto_distribution(ngram_counts, f'{n}-gram')\n    test_zipfs_law(ngram_counts, f'{n}-gram')","metadata":{"execution":{"iopub.status.busy":"2024-03-04T10:52:05.802594Z","iopub.execute_input":"2024-03-04T10:52:05.803701Z","iopub.status.idle":"2024-03-04T11:06:51.779965Z","shell.execute_reply.started":"2024-03-04T10:52:05.803659Z","shell.execute_reply":"2024-03-04T11:06:51.778605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Analyze and visualize for 1-grams, 11-grams\nfor n in [1, 11]:\n    ngram_counts = count_ngrams(sequences, n)\n    plot_ngram_frequencies(ngram_counts, f'{n}-gram')\n    analyze_ngram_benford(ngram_counts, f'{n}-gram')\n    plot_pareto_distribution(ngram_counts, f'{n}-gram')\n    test_zipfs_law(ngram_counts, f'{n}-gram')","metadata":{"execution":{"iopub.status.busy":"2024-03-04T11:08:52.31569Z","iopub.execute_input":"2024-03-04T11:08:52.316128Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Investigation of suspicious seqs as 'QQQQQQQQQQQ' and 'XXXXXXXXXXXX' and so on","metadata":{}},{"cell_type":"code","source":"# Identify unusually short or long sequences\n\naverage_length = np.mean([len(seq) for seq in sequences])\nstd_dev = np.std([len(seq) for seq in sequences])\n\n# Defining thresholds for unusual lengths (e.g., beyond 2 standard deviations)\nlower_threshold = average_length - 2 * std_dev\nupper_threshold = average_length + 2 * std_dev\n\n# Identifying sequences that are unusually short or long\nunusual_lengths = [seq for seq in sequences if len(seq) < lower_threshold or len(seq) > upper_threshold]\n\n# Visualization\nlengths = [len(seq) for seq in sequences]\nplt.hist(lengths, bins=150, color='skyblue')\nplt.axvline(lower_threshold, color='red', linestyle='dashed', linewidth=2)\nplt.axvline(upper_threshold, color='red', linestyle='dashed', linewidth=2)\nplt.title('Distribution of Sequence Lengths')\nplt.xlabel('Sequence Length')\nplt.ylabel('Frequency')\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2024-03-05T07:49:42.469353Z","iopub.execute_input":"2024-03-05T07:49:42.469684Z","iopub.status.idle":"2024-03-05T07:49:43.267688Z","shell.execute_reply.started":"2024-03-05T07:49:42.469659Z","shell.execute_reply":"2024-03-05T07:49:43.266363Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get 5 random samples from unusual lengths\nrandom_unusual_lengths = random.sample(unusual_lengths, min(5, len(unusual_lengths)))\n\nprint(\"Random samples from unusually short or long sequences:\")\nfor sample in random_unusual_lengths:\n    print(sample)","metadata":{"execution":{"iopub.status.busy":"2024-03-05T08:02:51.209531Z","iopub.execute_input":"2024-03-05T08:02:51.209924Z","iopub.status.idle":"2024-03-05T08:02:51.217161Z","shell.execute_reply.started":"2024-03-05T08:02:51.209887Z","shell.execute_reply":"2024-03-05T08:02:51.215605Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Identify sequences with unusually long repetitive amino acids\n\ndef has_long_repetition(sequence, threshold):\n    \"\"\"Checks if sequence has a long repetition of any amino acid.\"\"\"\n    for aa in set(sequence):\n        if sequence.find(aa * threshold) != -1:\n            return True\n    return False\n\n# Applying the function to filter sequences\nrepetitive_sequences = [seq for seq in sequences if has_long_repetition(seq, 11)]\n\n# Visualization of the count of sequences with long repetitions\nplt.figure(figsize=(6, 4))\nplt.bar(['Regular', 'Repetitive'], [len(sequences) - len(repetitive_sequences), len(repetitive_sequences)])\nplt.title('Comparison of Sequences with Long Repetitions')\nplt.ylabel('Number of Sequences')\nplt.show()\nprint(len(repetitive_sequences))","metadata":{"execution":{"iopub.status.busy":"2024-03-05T07:59:17.379022Z","iopub.execute_input":"2024-03-05T07:59:17.379351Z","iopub.status.idle":"2024-03-05T07:59:19.83074Z","shell.execute_reply.started":"2024-03-05T07:59:17.379324Z","shell.execute_reply":"2024-03-05T07:59:19.829647Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Get 5 random samples from repetitive sequences\nrandom_repetitive_sequences = random.sample(repetitive_sequences, min(5, len(repetitive_sequences)))\n\nprint(\"\\nRandom samples from sequences with unusually long repetitive amino acids:\")\nfor sample in random_repetitive_sequences:\n    print(sample)","metadata":{"execution":{"iopub.status.busy":"2024-03-05T08:03:34.966427Z","iopub.execute_input":"2024-03-05T08:03:34.966723Z","iopub.status.idle":"2024-03-05T08:03:34.972456Z","shell.execute_reply.started":"2024-03-05T08:03:34.966701Z","shell.execute_reply":"2024-03-05T08:03:34.971531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Identify sequences containing the letter 'X'\n\nsequences_with_x = [seq for seq in sequences if 'X' in seq]\n\n# Visualization of the proportion of sequences containing 'X'\nplt.figure(figsize=(6, 4))\nplt.bar(['Without X', 'With X'], [len(sequences) - len(sequences_with_x), len(sequences_with_x)])\nplt.title('Sequences Containing the Letter X')\nplt.ylabel('Number of Sequences')\nplt.show()\nprint(len(sequences_with_x))","metadata":{"execution":{"iopub.status.busy":"2024-03-05T08:09:58.751664Z","iopub.execute_input":"2024-03-05T08:09:58.751996Z","iopub.status.idle":"2024-03-05T08:09:58.935496Z","shell.execute_reply.started":"2024-03-05T08:09:58.751971Z","shell.execute_reply":"2024-03-05T08:09:58.934244Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Filter the already obtained random_sequences_with_x for those containing 'XX'\nsequences_with_xx_from_random = [seq for seq in sequences_with_x if 'XX' in seq]\n\n# Ensure there are enough sequences for sampling\nif len(sequences_with_xx_from_random) >= 5:\n    random_samples_with_xx = random.sample(sequences_with_xx_from_random, 5)\nelse:\n    random_samples_with_xx = sequences_with_xx_from_random\n\nprint(\"Random samples from sequences (already identified with 'X') containing 'XX' at least once:\")\nfor sample in random_samples_with_xx:\n    print(sample)","metadata":{"execution":{"iopub.status.busy":"2024-03-05T08:11:47.133303Z","iopub.execute_input":"2024-03-05T08:11:47.133631Z","iopub.status.idle":"2024-03-05T08:11:47.141798Z","shell.execute_reply.started":"2024-03-05T08:11:47.133606Z","shell.execute_reply":"2024-03-05T08:11:47.140269Z"},"trusted":true},"execution_count":null,"outputs":[]}]}