unwrapping-the-covid-19-genomes.ipynb

{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.6.6","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Unwrapping the genomes\nThis notebook is an attempt to answering some of the questions of bio-informatics in regards to the recent panademic outbreak of Corona virus, from a non-medical background person's perspective. Every analysis presented here is based on the stuff I learnt and understood (or felt that I actually understood) about the genomes from the internet.","metadata":{}},{"cell_type":"markdown","source":"## Table of Contents:\n- [Some Context](#Some-Context)\n- [Comparison of composition of genomes](#Comparison-of-composition-of-genomes)\n    - [Oligonucleotide composition](#Oligonucleotide-composition)\n    - [GC Content](#GC-Content)\n- [Proteins & Amino Acids](#Proteins-Amino-Acids)\n    - [Protein Strands by genome](#Protein-Strands-by-genome)\n    - [Amino Acids Distribution](#Amino-Acids-Distribution)\n    - [Finding the ORFs](#Finding-the-ORFs)\n    - [Summary of Protein Sequences](#summary-of-protein-seq)\n- [Sequence Alignment](#Sequence-Alignement)\n- [Pairwise Sequence Alignment](#Pairwise-Sequence-Alignment)\n- [Genome Sequence Similarity](#Genome-Sequence-Similarity)\n- [Generating a Phylogenetic Tree](#Generating-Phylogenetic-Tree)\n- [Clustering all COVID-19 patients' genomes](#Clustering-patients-genomes)\n- [Conclusions](#Conclusions)","metadata":{}},{"cell_type":"markdown","source":"<a id=\"Some-Context\"></a>\n# Some Context\n**[Ques.] What is Genome?**<br>\n[Genomes](https://en.wikipedia.org/wiki/Genome) are considered to be genetic material of any organism. It consists of [DNA](https://en.wikipedia.org/wiki/DNA)(or [RNA](https://en.wikipedia.org/wiki/RNA)) which in turn can be considered as the blueprints of any organisms' origin.\n<img src=\"https://www.genome.gov/sites/default/files/tg/en/illustration/acgt.jpg\" width=550 height=400>\n\n**[Ques.] How to study Genomes?**<br>\nA genome consists of a sequence of [nucleotides](https://en.wikipedia.org/wiki/Nucleotide) that together make up all the [chromosomes](https://en.wikipedia.org/wiki/Chromosome) of any organism. These nucleotides are the basic building blocks of DNA and RNA. A DNA genome has 4 types of base nucleotides -<br> `Adenine` also represented as `A`,<br>`Thymine` represented as `T`,<br>`Guanine` represented as `G` and<br>`Cytosine` represented as `C`.<br> In RNA genome, the `Thymine` nucleotide is replaced with what is known as `Uracil` represented as `U`.<br> In short, any genome sequence is basically a combination of all these nucleotide in smoe order.\n<img src=\"https://upload.wikimedia.org/wikipedia/commons/3/37/Difference_DNA_RNA-EN.svg\" width=550 height=600>\n\n**[Ques.] Can one find a cure for [COVID-19](https://en.wikipedia.org/wiki/Coronavirus_disease_2019) by just studying the genomes?**<br>\n**Not sure**, as the [molecular biology](https://en.wikipedia.org/wiki/Molecular_biology) has not really been my field of study (*never liked Biology tbh*) but also because its a very vast subject where to confirm on a single hypothesis also takes months/years of evaluation and testing. Having said that, the whole purpose of this analysis is to help the ones who actually understand molecular biology get some data-driven insights to be able to preapare a RNA vaccine at the earliest.\n<img src=\"https://images.theconversation.com/files/324780/original/file-20200402-23151-18gxxtu.png\" width=500 height=500>","metadata":{}},{"cell_type":"markdown","source":"## Loading libraries & genomes","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 in \n\nimport numpy as np\nimport pandas as pd\nimport re\nimport plotly.offline as py\nimport plotly.graph_objects as go\nfrom plotly.subplots import make_subplots\n\nimport os\nos.environ['QT_QPA_PLATFORM']='offscreen' # https://github.com/etetoolkit/ete/issues/296\n\nimport matplotlib.pyplot as plt\n%matplotlib inline\nfrom colorama import Back, Style, Fore\n\ntry:\n    from Bio.Align.Applications import ClustalwCommandline\n    from Bio import Align, pairwise2\nexcept:\n    !pip install biopython\n    from Bio.Align.Applications import ClustalwCommandline\n    from Bio import Align, pairwise2\ntry:\n    from ete3 import Tree, TreeStyle\nexcept:\n    !pip install ete3\n    from ete3 import Tree, TreeStyle","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-11-27T06:28:45.211147Z","iopub.execute_input":"2021-11-27T06:28:45.211466Z","iopub.status.idle":"2021-11-27T06:29:02.099768Z","shell.execute_reply.started":"2021-11-27T06:28:45.211437Z","shell.execute_reply":"2021-11-27T06:29:02.098683Z"},"trusted":true},"execution_count":1,"outputs":[]},{"cell_type":"code","source":"# covid-19 genome\nwith open(\"/kaggle/input/coronavirus-accession-sars-mers-cov2/sars.fasta\", \"r\") as file:\n    sars_genome = file.read().split(\"genome\")[1].replace(\"\\n\", \"\")\n# mers genome\nwith open (\"/kaggle/input/coronavirus-accession-sars-mers-cov2/mers.fasta\", \"r\") as file:\n    mers_genome = file.read().split(\"genome\")[1].replace(\"\\n\", \"\")\n# sars genome\nwith open (\"/kaggle/input/coronavirus-accession-sars-mers-cov2/cov2.fasta\", \"r\") as file:\n    cov_genome = file.read().split(\"genome\")[1].replace(\"\\n\", \"\")\n# civet with sars like cov genome\nwith open(\"/kaggle/input/ebolav-vs-sarscov-vs-mers/Civet-SARS.fasta\", \"r\") as file:\n    civet_sl_cov_genome = file.read().split(\"complete genome\")[-1].replace(\"\\n\",\"\")\n# bat with sars like cov genome\nwith open(\"/kaggle/input/ebolav-vs-sarscov-vs-mers/BAT-SARS.fasta\", \"r\") as file:\n    bat_sl_cov_genome = file.read().split(\"complete genome\")[-1].replace(\"\\n\",\"\")\n# ebola-5 genome\nwith open(\"/kaggle/input/ebolav-vs-sarscov-vs-mers/EBOLAV.fasta\", \"r\") as file:\n    ebola5_genome = file.read().split(\"complete genome\")[-1].replace(\"\\n\",\"\")\n# camel with cov genome\nwith open(\"/kaggle/input/cameluscovfasta/camelus.fasta\", \"r\") as file:\n    camel_cov_genome = file.read().split(\"complete genome\")[-1].replace(\"\\n\", \"\")\n# malaria genome\nwith open(\"/kaggle/input/cameluscovfasta/plasmodium-malariae.fasta\", \"r\") as file:\n    malaria_genome = file.read().split(\"complete sequence\")[-1].replace(\"\\n\", \"\")\n# hiv-2 genome\nwith open(\"/kaggle/input/cameluscovfasta/hiv2.fasta\", \"r\") as file:\n    hiv2_genome = file.read().split(\"complete genome\")[-1].replace(\"\\n\", \"\")\n# hedgehog with cov genome\nwith open(\"/kaggle/input/cameluscovfasta/hedgehog.fasta\", \"r\") as file:\n    hedgehog_cov_genome = file.read().split(\"complete genome\")[-1].replace(\"\\n\", \"\")\n","metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","_kg_hide-input":true,"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:02.101753Z","iopub.execute_input":"2021-11-27T06:29:02.102008Z","iopub.status.idle":"2021-11-27T06:29:02.206618Z","shell.execute_reply.started":"2021-11-27T06:29:02.101977Z","shell.execute_reply":"2021-11-27T06:29:02.205834Z"},"trusted":true},"execution_count":2,"outputs":[]},{"cell_type":"markdown","source":"lets have a look at a snippet of one of the genomes","metadata":{}},{"cell_type":"code","source":"hiv2_genome[:3000]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:02.210307Z","iopub.execute_input":"2021-11-27T06:29:02.210627Z","iopub.status.idle":"2021-11-27T06:29:02.218547Z","shell.execute_reply.started":"2021-11-27T06:29:02.210595Z","shell.execute_reply":"2021-11-27T06:29:02.217615Z"},"trusted":true},"execution_count":3,"outputs":[]},{"cell_type":"markdown","source":"having a huge sequence and with everyting in single color could be difficult to read. So let's color-code these nucleotides in the sequences and have them represented as they are normally written.","metadata":{}},{"cell_type":"code","source":"# some helper functions for representation of sequences\n\n# usually DNA sequences are written in 10-nucleotide fashion\ndef ten_nucleotide_seq(genome):\n    genes = []\n    for ix, char in enumerate(genome):\n        if ix != 0 and ix%10 == 0:\n            genes.append(' ')\n        genes.append(char)\n    return ''.join(genes)\n\n# color code to represent genome sequences\nnu_clr_switcher = {\n    # standard color-codes\n    'A': Back.GREEN,\n    'C': Back.YELLOW,\n    'G': Back.RED,\n    'T': Back.BLUE,\n    ' ': Style.RESET_ALL\n}\nprotein_clr_switcher = {\n    # color-code by proteinfamily's polarity\n    'A': Back.BLUE,\n    'V': Back.BLUE,\n    'I': Back.BLUE,\n    'L': Back.BLUE,\n    'M': Back.BLUE,\n    'F': Back.BLUE,\n    'Y': Back.CYAN,\n    'W': Back.BLUE,\n    'H': Back.CYAN,\n    'R': Back.RED,\n    'K': Back.RED,\n    'N': Back.GREEN,\n    'Q': Back.GREEN,\n    'E': Back.MAGENTA,\n    'D': Back.MAGENTA,\n    'S': Back.GREEN,\n    'T': Back.GREEN,\n    'G': Back.YELLOW,\n    'P': Back.YELLOW,\n    'C': Back.BLUE,\n    ' ': Style.RESET_ALL\n}\ndef seq_repr(genome_str, strand ='dna'):\n    if strand == 'dna':\n        genome_str = ten_nucleotide_seq(genome=genome_str)\n        line_break_cntr = 0\n        for i in range(len(genome_str)):\n            if genome_str[i] == ' ':\n                line_break_cntr += 1\n                if line_break_cntr>0 and line_break_cntr%6==0:\n                    text = \"\\n\"\n                else:\n                    text = nu_clr_switcher[genome_str[i]] + genome_str[i]\n            else:\n                text = nu_clr_switcher[genome_str[i]] + genome_str[i]\n            print(text, end=\"\")\n        Style.RESET_ALL\n        \n    if strand == 'protein':\n        for i in range(len(genome_str)):\n            if genome_str[i] in protein_clr_switcher:\n                if genome_str[i] == 'S' and genome_str[i+1:i+4] == 'TOP':\n                    text = Style.RESET_ALL + 'S'\n                elif genome_str[i] == 'T' and genome_str[i-1] == 'S' and genome_str[i+1:i+3] == 'OP':\n                    text = Style.RESET_ALL + 'T'\n                elif genome_str[i] == 'P' and genome_str[i-3:i] == 'STO':\n                    text = Style.RESET_ALL + 'P'\n                else:\n                    text = protein_clr_switcher[genome_str[i]] + genome_str[i]\n            else:\n                Style.RESET_ALL\n                text = genome_str[i]\n            print(text, end=\"\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:02.220195Z","iopub.execute_input":"2021-11-27T06:29:02.220806Z","iopub.status.idle":"2021-11-27T06:29:02.245351Z","shell.execute_reply.started":"2021-11-27T06:29:02.220761Z","shell.execute_reply":"2021-11-27T06:29:02.243993Z"},"trusted":true},"execution_count":4,"outputs":[]},{"cell_type":"code","source":"print(\"COVID-19 genome: \")\nseq_repr(cov_genome)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:02.248502Z","iopub.execute_input":"2021-11-27T06:29:02.248823Z","iopub.status.idle":"2021-11-27T06:29:05.698020Z","shell.execute_reply.started":"2021-11-27T06:29:02.248791Z","shell.execute_reply":"2021-11-27T06:29:05.697052Z"},"trusted":true},"execution_count":5,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Comparison-of-composition-of-genomes\"></a>\n# Comparison of composition of genomes\nComparing the genomic composition, helps Vaccinologist do [reverse vaccinology](https://en.wikipedia.org/wiki/Reverse_vaccinology) to discover candidate [antigents](https://en.wikipedia.org/wiki/Antigen) for vaccine development by analyzing the genome of a pathogen or a family of pathogens.\n\nFor the comparison of composition of these genomes, I'll be looking at the<br>\n- Oligonucleotide composition<br>\n- GC content\n\n<a id=\"Oligonucleotide-composition\"></a>\n### 1. Oligonucleotide composition:\n[Oligonucleotides](https://en.wikipedia.org/wiki/Oligonucleotide) are DNA or RNA molecules made up of a sequence of nucleotides. The length of the oligonucleotide is represented as \"[k-mer](https://en.wikipedia.org/wiki/K-mer)\" ($k$ being the number of nucleotides in one oligonulceotide). In this analysis, I'll be looking at the **3-mers** (or **trinucleotides**) and **4-mer** (or **tetranucleotides**) compositions against their normalized frquency of occurrences.","metadata":{}},{"cell_type":"code","source":"# tri-nucleotide compositions\ntrimers = [\"AAA\", \"AAC\", \"AAG\", \"AAT\", \"ACA\", \"ACC\", \"ACG\", \"ACT\", \"AGA\", \"AGC\", \"AGG\", \"ATA\", \"ATC\", \"ATG\", \"CAA\", \n           \"CAC\", \"CAG\", \"CCA\",\"CCC\",\"CCG\",\"CGA\",\"CGC\",\"CTA\",\"CTC\",\"GAA\",\"GAC\",\"GCA\",\"GCC\",\"GGA\",\"GTA\",\"TAA\",\"TCA\"]\n\ndef trimer_composition(genome):\n    trimer_dict = dict()\n    for trimer in trimers:\n        trimer_dict[trimer] = genome.count(trimer)\n    return trimer_dict\n\n# tetra-nucleotide compositions\ntetramers = ['AAAA','AAAC','AAAG','AAAT','AACA','AACC','AACG','AACT','AAGA','AAGC','AAGG','AAGT','AATA','AATC',\n             'AATG','AATT','ACAA','ACAC','ACAG','ACAT','ACCA','ACCC','ACCG','ACCT','ACGA','ACGC','ACGG','ACGT',\n             'ACTA','ACTC','ACTG','AGAA','AGAC','AGAG','AGAT','AGCA','AGCC','AGCG','AGCT','AGGA','AGGC','AGGG',\n             'AGTA','AGTC','AGTG','ATAA','ATAC','ATAG','ATAT','ATCA','ATCC','ATCG','ATGA','ATGC','ATGG','ATTA',\n             'ATTC','ATTG','CAAA','CAAC','CAAG','CACA','CACC','CACG','CAGA','CAGC','CAGG','CATA','CATC','CATG',\n             'CCAA','CCAC','CCAG','CCCA','CCCC','CCCG','CCGA','CCGC','CCGG','CCTA','CCTC','CGAA','CGAC','CGAG',\n             'CGCA','CGCC','CGCG','CGGA','CGGC','CGTA','CGTC','CTAA','CTAC','CTAG','CTCA','CTCC','CTGA','CTGC',\n             'CTTA','CTTC','GAAA','GAAC','GACA','GACC','GAGA','GAGC','GATA','GATC','GCAA','GCAC','GCCA','GCCC',\n             'GCGA','GCGC','GCTA','GGAA','GGAC','GGCA','GGCC','GGGA','GGTA','GTAA','GTAC','GTCA','GTGA','GTTA',\n             'TAAA','TACA','TAGA','TATA','TCAA','TCCA','TCGA','TGAA','TGCA','TTAA']\n\ndef tetramer_composition(genome):\n    tetramer_dict = dict()\n    for tetramer in tetramers:\n        if tetramer in genome:\n            tetramer_dict[tetramer] = genome.count(tetramer)\n        else:\n            tetramer_dict[tetramer] = 0\n    return tetramer_dict","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:05.703590Z","iopub.execute_input":"2021-11-27T06:29:05.703912Z","iopub.status.idle":"2021-11-27T06:29:05.721576Z","shell.execute_reply.started":"2021-11-27T06:29:05.703881Z","shell.execute_reply":"2021-11-27T06:29:05.720590Z"},"trusted":true},"execution_count":6,"outputs":[]},{"cell_type":"code","source":"labels = [\"SARS\", \"MERS\", \"COVID-19\", \"Civet_SL_CoV\",\"Bat_SL_CoV\",\"Ebola\",\"Camel_CoV\",\"Malaria\",\"HIV\",\"Hedgehog_CoV\"]\n\ntraces = []\nfor i, genome in enumerate((sars_genome,mers_genome,cov_genome,civet_sl_cov_genome,bat_sl_cov_genome,ebola5_genome,\n                            camel_cov_genome,malaria_genome,hiv2_genome,hedgehog_cov_genome)):\n    composition = trimer_composition(genome)\n    total_composition = sum(composition.values())\n    norm_freq = [count/total_composition for count in composition.values()]\n    traces.append(go.Scatter(x=trimers,y=norm_freq,mode='lines',name=labels[i]))\n    \nlayout=go.Layout(\n    title=\"Trinucleotide Composition of genomes\",\n    xaxis_title=\"Trinucleotides\",\n    yaxis_title=\"Normalized frequencies\",\n    width=800, height=600\n)\n\nfig = go.Figure(data=traces, layout=layout)\npy.iplot(fig, filename=\"3mer-composition-frequency\")\n","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:05.723572Z","iopub.execute_input":"2021-11-27T06:29:05.723967Z","iopub.status.idle":"2021-11-27T06:29:07.799655Z","shell.execute_reply.started":"2021-11-27T06:29:05.723920Z","shell.execute_reply":"2021-11-27T06:29:07.798627Z"},"trusted":true},"execution_count":7,"outputs":[]},{"cell_type":"markdown","source":"Now lets have a look at the tetramer composition of these genomes","metadata":{}},{"cell_type":"code","source":"labels = [\"SARS\", \"MERS\", \"COVID-19\", \"Civet_SL_CoV\",\"Bat_SL_CoV\",\"Ebola\",\"Camel_CoV\",\"Malaria\",\"HIV\",\"Hedgehog_CoV\"]\n\ntraces = []\nfor i, genome in enumerate((sars_genome,mers_genome,cov_genome,civet_sl_cov_genome,bat_sl_cov_genome,ebola5_genome,\n                            camel_cov_genome,malaria_genome,hiv2_genome,hedgehog_cov_genome)):\n    composition = tetramer_composition(genome)\n    total_composition = sum(composition.values())\n    norm_freq = [count/total_composition for count in composition.values()]\n    traces.append(go.Scatter(x=tetramers,y=norm_freq,mode='lines',name=labels[i]))\n    \nlayout=go.Layout(\n    title=\"Tetranucleotides Composition of genomes\",\n    xaxis_title=\"Tetranucleotides\",\n    yaxis_title=\"Normalized frequencies\",\n    width=1000, height=600\n)\n\nfig = go.Figure(data=traces, layout=layout)\npy.iplot(fig, filename=\"4mer-composition-frequency\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:07.801514Z","iopub.execute_input":"2021-11-27T06:29:07.801875Z","iopub.status.idle":"2021-11-27T06:29:08.598351Z","shell.execute_reply.started":"2021-11-27T06:29:07.801832Z","shell.execute_reply":"2021-11-27T06:29:08.597413Z"},"trusted":true},"execution_count":8,"outputs":[]},{"cell_type":"markdown","source":"Well that's alot to put in one plot, but the good thing with **plotly** is that you can always turn off a few legends to compare the values just among the selected few.\n\nFor a moment, lets just look at the tetranucleotide composition of all the \"corona-virus infected species\".","metadata":{}},{"cell_type":"code","source":"labels = [\"SARS\", \"MERS\", \"COVID-19\", \"Civet_SL_CoV\",\"Bat_SL_CoV\",\"Camel_CoV\",\"Hedgehog_CoV\"]\n\ntraces = []\nfor i, genome in enumerate((sars_genome,mers_genome,cov_genome,civet_sl_cov_genome,\n                            bat_sl_cov_genome,camel_cov_genome,hedgehog_cov_genome)):\n    composition = tetramer_composition(genome)\n    total_composition = sum(composition.values())\n    norm_freq = [count/total_composition for count in composition.values()]\n    traces.append(go.Scatter(x=tetramers,y=norm_freq,mode='lines',name=labels[i]))\n    \nlayout=go.Layout(\n    title=\"Tetranucleotides Composition of genomes\",\n    xaxis_title=\"Tetranucleotides\",\n    yaxis_title=\"Normalized frequencies\",\n    width=1000, height=600\n)\n\nfig = go.Figure(data=traces, layout=layout)\npy.iplot(fig, filename=\"4mer-composition-frequency\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:08.600357Z","iopub.execute_input":"2021-11-27T06:29:08.600715Z","iopub.status.idle":"2021-11-27T06:29:09.336400Z","shell.execute_reply.started":"2021-11-27T06:29:08.600675Z","shell.execute_reply":"2021-11-27T06:29:09.335404Z"},"trusted":true},"execution_count":9,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"GC-Content\"></a>\n### 2. GC Content:\n**[GC content](https://en.wikipedia.org/wiki/GC-content)** (or **Guanine-Cytosine content**) is the percentage of nucleotides bases in a DNA/RNA molecule that are either guanine (G) or cytosine (C). GC content is always expressed as a percentage value and is supposed to remaine the same among genomes of same specie.\n\n$$\\frac{G + C}{A + C + G + T} \\times 100\\ \\%$$","metadata":{}},{"cell_type":"code","source":"def gc_content(genome):\n    g_content = genome.count(\"G\")\n    c_content = genome.count(\"C\")\n    a_content = genome.count(\"A\")\n    t_content = genome.count(\"T\")\n    return 100*((g_content+c_content)/(a_content+c_content+g_content+t_content))\n\ngc_content_perc = [gc_content(sars_genome),gc_content(mers_genome),gc_content(cov_genome),gc_content(civet_sl_cov_genome),\n                   gc_content(bat_sl_cov_genome),gc_content(ebola5_genome),gc_content(camel_cov_genome),\n                   gc_content(malaria_genome),gc_content(hiv2_genome), gc_content(hedgehog_cov_genome)]\nviruses = [\"SARS\", \"MERS\", \"COVID-19\", \"Civet_SL_CoV\",\"Bat_SL_CoV\",\"Ebola\",\"Camel_CoV\",\"Malaria\",\"HIV\",\"Hedgehog_CoV\"]\ncolors = [\"lightslategray\",] * 10\ncolors[2] = \"crimson\"\n\nfig = go.Figure(data=[go.Bar(\n    x=viruses,\n    y=gc_content_perc,\n    text=gc_content_perc,\n    marker_color=colors\n)])\nfig.update_traces(texttemplate=\"%{text:.2f}%\", textposition=\"outside\")\nfig.update_layout(\n    title_text=\"GC content distribution among genomes\",\n    xaxis_title=\"Genomes\",\n    yaxis_title=\"GC content (%)\"\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:09.338181Z","iopub.execute_input":"2021-11-27T06:29:09.338511Z","iopub.status.idle":"2021-11-27T06:29:09.531030Z","shell.execute_reply.started":"2021-11-27T06:29:09.338474Z","shell.execute_reply":"2021-11-27T06:29:09.530059Z"},"trusted":true},"execution_count":10,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Proteins-Amino-Acids\"></a>\n# Proteins & Amino Acids\nThe most commonly occuring amino acids are abbreviated by using 20 letters from the english alphabet (all letters except for B,J,O,U,X and Z). The protein strings are constructed from these 20 alphabets as symbols. Henceforth, a Protein Sequence can also be incorporated using the nucleotide sequences.\n\n<img src=\"https://1.bp.blogspot.com/-QpZ1VXFbJhE/V4d5EjT3FYI/AAAAAAAADvk/UGZK71q3KV4AwrOfR_bdEk2FnC93Cb-JwCLcB/s1600/codon.png\" width=550 height=400>\n\nThe DNA codon table on [wiki](https://en.wikipedia.org/wiki/DNA_codon_table) dictates the details regarding the encoding of specific codons into amino acids alphabets.\n\n","metadata":{}},{"cell_type":"code","source":"def gen_protein_seq(genome_str):\n    # codon table from wiki\n    protein = {\"TTT\" : \"F\", \"CTT\" : \"L\", \"ATT\" : \"I\", \"GTT\" : \"V\",\n           \"TTC\" : \"F\", \"CTC\" : \"L\", \"ATC\" : \"I\", \"GTC\" : \"V\",\n           \"TTA\" : \"L\", \"CTA\" : \"L\", \"ATA\" : \"I\", \"GTA\" : \"V\",\n           \"TTG\" : \"L\", \"CTG\" : \"L\", \"ATG\" : \"M\", \"GTG\" : \"V\",\n           \"TCT\" : \"S\", \"CCT\" : \"P\", \"ACT\" : \"T\", \"GCT\" : \"A\",\n           \"TCC\" : \"S\", \"CCC\" : \"P\", \"ACC\" : \"T\", \"GCC\" : \"A\",\n           \"TCA\" : \"S\", \"CCA\" : \"P\", \"ACA\" : \"T\", \"GCA\" : \"A\",\n           \"TCG\" : \"S\", \"CCG\" : \"P\", \"ACG\" : \"T\", \"GCG\" : \"A\",\n           \"TAT\" : \"Y\", \"CAT\" : \"H\", \"AAT\" : \"N\", \"GAT\" : \"D\",\n           \"TAC\" : \"Y\", \"CAC\" : \"H\", \"AAC\" : \"N\", \"GAC\" : \"D\",\n           \"TAA\" : \"STOP\", \"CAA\" : \"Q\", \"AAA\" : \"K\", \"GAA\" : \"E\",\n           \"TAG\" : \"STOP\", \"CAG\" : \"Q\", \"AAG\" : \"K\", \"GAG\" : \"E\",\n           \"TGT\" : \"C\", \"CGT\" : \"R\", \"AGT\" : \"S\", \"GGT\" : \"G\",\n           \"TGC\" : \"C\", \"CGC\" : \"R\", \"AGC\" : \"S\", \"GGC\" : \"G\",\n           \"TGA\" : \"STOP\", \"CGA\" : \"R\", \"AGA\" : \"R\", \"GGA\" : \"G\",\n           \"TGG\" : \"W\", \"CGG\" : \"R\", \"AGG\" : \"R\", \"GGG\" : \"G\" \n           }\n    protein_seq = \"\"\n    \n    # generating the protein seq\n    for i in range(0, len(genome_str)-(3+len(genome_str)%3), 3):\n        protein_seq += protein[genome_str[i:i+3]]\n    return protein_seq","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:09.533048Z","iopub.execute_input":"2021-11-27T06:29:09.533412Z","iopub.status.idle":"2021-11-27T06:29:09.547481Z","shell.execute_reply.started":"2021-11-27T06:29:09.533370Z","shell.execute_reply":"2021-11-27T06:29:09.546383Z"},"trusted":true},"execution_count":11,"outputs":[]},{"cell_type":"markdown","source":"lets have a look at the protein sequences of a few of these genome before we go for a comparative analysis","metadata":{}},{"cell_type":"code","source":"print(\"Protein Sequences in Malaria: \")\nseq_repr(gen_protein_seq(malaria_genome), strand='protein')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:09.549137Z","iopub.execute_input":"2021-11-27T06:29:09.549515Z","iopub.status.idle":"2021-11-27T06:29:09.648539Z","shell.execute_reply.started":"2021-11-27T06:29:09.549455Z","shell.execute_reply":"2021-11-27T06:29:09.647688Z"},"trusted":true},"execution_count":12,"outputs":[]},{"cell_type":"code","source":"print(\"Protein sequences in CoV: \")\nseq_repr(gen_protein_seq(cov_genome), strand='protein')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:09.653158Z","iopub.execute_input":"2021-11-27T06:29:09.653423Z","iopub.status.idle":"2021-11-27T06:29:10.927359Z","shell.execute_reply.started":"2021-11-27T06:29:09.653396Z","shell.execute_reply":"2021-11-27T06:29:10.926246Z"},"trusted":true},"execution_count":13,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Protein-Strands-by-genome\"></a>\n### 1. Protein Strands by genome","metadata":{}},{"cell_type":"code","source":"genomes = {\"SARS\": sars_genome,\"MERS\": mers_genome,\"COVID-19\": cov_genome,\"Civet_SL_CoV\": civet_sl_cov_genome,\n           \"Bat_SL_CoV\": bat_sl_cov_genome,\"Ebola\": ebola5_genome,\"Camel_CoV\": camel_cov_genome,\n           \"Malaria\": malaria_genome,\"HIV\": hiv2_genome,\"Hedgehog_CoV\": hedgehog_cov_genome}\n\ncolors = [\"lightslategray\",] * 10\ncolors[2] = \"crimson\"\n\n# for virus_name, virus_gene in genomes.items():\n#     print(virus_name, len(gen_protein_seq(virus_gene).split('STOP')))\n\nfig = go.Figure(data=[go.Bar(\n    x=list(genomes.keys()),\n    y=[len(gen_protein_seq(genome).split('STOP')) for genome in genomes.values()],\n    text=[len(gen_protein_seq(genome).split('STOP')) for genome in genomes.values()],\n    marker_color=colors\n)])\n\nfig.update_traces(texttemplate=\"%{text:.2s}\", textposition=\"outside\")\nfig.update_layout(\n    title_text=\"Number of protein strands among genomes\",\n    xaxis_title=\"Genomes\",\n    yaxis_title=\"number of protein strands\"\n)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:10.928646Z","iopub.execute_input":"2021-11-27T06:29:10.928866Z","iopub.status.idle":"2021-11-27T06:29:11.214182Z","shell.execute_reply.started":"2021-11-27T06:29:10.928842Z","shell.execute_reply":"2021-11-27T06:29:11.212992Z"},"trusted":true},"execution_count":14,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Amino-Acids-Distribution\"></a>\n### 2. Amino Acids Distribution","metadata":{}},{"cell_type":"code","source":"def amino_acid_composition(protein_seq):\n    amino_acids_dict = {'Y': 0, 'A': 0, 'Q': 0, 'D': 0, 'C': 0, 'G': 0, 'V': 0, 'T': 0, 'E': 0, 'N': 0, \n                       'K': 0, 'R': 0, 'S': 0, 'I': 0, 'H': 0, 'M': 0, 'F': 0, 'L': 0, 'W': 0, 'P': 0}\n    for amino_acid in amino_acids_dict:\n        amino_acids_dict[amino_acid] = protein_seq.count(amino_acid)/len(protein_seq)*100\n    return amino_acids_dict","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:11.215758Z","iopub.execute_input":"2021-11-27T06:29:11.216112Z","iopub.status.idle":"2021-11-27T06:29:11.225301Z","shell.execute_reply.started":"2021-11-27T06:29:11.216053Z","shell.execute_reply":"2021-11-27T06:29:11.224204Z"},"trusted":true},"execution_count":15,"outputs":[]},{"cell_type":"code","source":"fig = go.Figure()\nfor i in range(10):\n    fig.add_trace(go.Bar(\n        x=list(amino_acid_composition(gen_protein_seq(list(genomes.values())[i])).keys()),\n        y=list(amino_acid_composition(gen_protein_seq(list(genomes.values())[i])).values()),\n        name=list(genomes.keys())[i]\n    ))\n\n# Here we modify the tickangle of the xaxis, resulting in rotated labels.\nfig.update_layout(barmode=\"group\",\n                  title=\"Distribution of Amino Acids by genomes\", \n                  xaxis_title=\"Amino Acids\", \n                  yaxis_title=\"Normalized Frequency (%)\")\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:11.226966Z","iopub.execute_input":"2021-11-27T06:29:11.227315Z","iopub.status.idle":"2021-11-27T06:29:11.720585Z","shell.execute_reply.started":"2021-11-27T06:29:11.227276Z","shell.execute_reply":"2021-11-27T06:29:11.719709Z"},"trusted":true},"execution_count":16,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Finding-the-ORFs\"></a>\n### 3. Finding the ORFs (Open Reading Frames)\nThe stretch of codons from start-codon (ATG) to stop-codon (TAA/TAG/TGA) is called [Open Reading Frame (ORF)](https://en.wikipedia.org/wiki/Open_reading_frame)\n","metadata":{}},{"cell_type":"code","source":"def get_orfs(protein_seq):\n    orf_strands = []\n    # generate ORFs\n    for seq in protein_seq.split(\"STOP\"):\n        for nu in range(len(seq)-1, -1, -1):\n            if seq[nu] == \"M\":\n                orf_strands.append(seq[nu:]+\"STOP\")\n    # search in original seq\n    patterns = \"|\".join(orf_strands)\n    res_seq = re.sub(patterns, lambda x: x.group(0).lower(), protein_seq)\n    # highlight the ORFs in the seq\n    for nu in res_seq:\n        if nu.isupper():\n            print(Style.RESET_ALL + nu, end=\"\")\n        else:\n            print(Back.BLUE + nu.upper(), end=\"\")\n    print(Style.RESET_ALL+\"\\n\")","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:11.722694Z","iopub.execute_input":"2021-11-27T06:29:11.722983Z","iopub.status.idle":"2021-11-27T06:29:11.732398Z","shell.execute_reply.started":"2021-11-27T06:29:11.722953Z","shell.execute_reply":"2021-11-27T06:29:11.731640Z"},"trusted":true},"execution_count":17,"outputs":[]},{"cell_type":"markdown","source":"the highlighted sequences represent the ORFs","metadata":{}},{"cell_type":"code","source":"for genome_name, genome in genomes.items():\n    print(\"\\nORFs in {}'s protein sequences:\".format(genome_name))\n    get_orfs(gen_protein_seq(genome))\n    print()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:11.733816Z","iopub.execute_input":"2021-11-27T06:29:11.734106Z","iopub.status.idle":"2021-11-27T06:29:23.738730Z","shell.execute_reply.started":"2021-11-27T06:29:11.734053Z","shell.execute_reply":"2021-11-27T06:29:23.737857Z"},"trusted":true},"execution_count":18,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"summary-of-protein-seq\"></a>\n### 4. Comparative Summary of Protein Sequences","metadata":{}},{"cell_type":"code","source":"genomes = {\"SARS\": sars_genome,\"MERS\": mers_genome,\"COVID-19\": cov_genome,\"Civet_SL_CoV\": civet_sl_cov_genome,\n           \"Bat_SL_CoV\": bat_sl_cov_genome,\"Ebola\": ebola5_genome,\"Camel_CoV\": camel_cov_genome,\n           \"Malaria\": malaria_genome,\"HIV\": hiv2_genome,\"Hedgehog_CoV\": hedgehog_cov_genome}\n\ncolors = [\"lightslategray\",] * 10\ncolors[2] = \"crimson\"\n\nproteins_seq = dict()\nfor genome_name, genome in genomes.items():\n    proteins_seq[genome_name] = gen_protein_seq(genome)\n    \nproteins_seq_lengths = dict()\nfor genome_name, protein_seq in proteins_seq.items():\n    len_of_pro_seq = []\n    for seq in protein_seq.split('STOP'):\n        len_of_pro_seq.append(len(seq)+1) # adding 1 for including STOP codon in total len of protein seq\n    proteins_seq_lengths[genome_name] = len_of_pro_seq\n\nfig = make_subplots(\n    rows=2, cols=2,\n    specs=[[{}, {}],\n           [{\"colspan\": 2}, None]],\n    subplot_titles=(\"Length of Shortest Protein Seq\", \"Length of Longest Protein Seq\", \"Average Length of Protein Seq\"))\n\nfig.add_trace(\n    go.Bar(\n    x=list(proteins_seq_lengths.keys()),\n    y=[min(_len) for _len in proteins_seq_lengths.values()],\n    text=[min(_len) for _len in proteins_seq_lengths.values()],\n    texttemplate=\"%{text:.1s}\", textposition=\"outside\",\n    marker_color=colors),\n    row=1, col=1\n)\n\nfig.add_trace(\n    go.Bar(\n    x=list(proteins_seq_lengths.keys()),\n    y=[max(_len) for _len in proteins_seq_lengths.values()],\n    text=[max(_len) for _len in proteins_seq_lengths.values()],\n    texttemplate=\"%{text:.1s}\", textposition=\"outside\",\n    marker_color=colors),\n    row=1, col=2\n)\n\nfig.add_trace(\n    go.Bar(\n    x=list(proteins_seq_lengths.keys()),\n    y=[sum(lens)/len(lens) for lens in proteins_seq_lengths.values()],\n    text=[sum(lens)/len(lens) for lens in proteins_seq_lengths.values()],\n    texttemplate=\"%{text:.2f}\", textposition=\"outside\",\n    marker_color=colors),\n    row=2, col=1\n)\n\nfig.update_layout(title_text=\"Comparative Summary of Protein Sequences\",\n                  height=960, width=960, showlegend=False)\n\nfig.show()","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:23.740309Z","iopub.execute_input":"2021-11-27T06:29:23.740538Z","iopub.status.idle":"2021-11-27T06:29:24.073597Z","shell.execute_reply.started":"2021-11-27T06:29:23.740513Z","shell.execute_reply":"2021-11-27T06:29:24.072224Z"},"trusted":true},"execution_count":19,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Sequence-Alignement\"></a>\n# Sequence Alignement\nThe sequence alignment is a way of arranging the sequences of DNA/RNA or protein to identify regions of similarity that may be a consequence of functional, structural, or evolutionary relationships between the sequences. These in turn help in identifying the sequence similarities and developing homology models of protein structures. Alignments are often assumed to reflect a degree of evolutionary change between sequences descended from a common ancestor.","metadata":{}},{"cell_type":"code","source":"sars_genome_aligned = ClustalwCommandline('clustalw2', infile=sars_genome)\nmers_genome_aligned = ClustalwCommandline('clustalw2', infile=mers_genome)\ncov_genome_aligned = ClustalwCommandline('clustalw2', infile=cov_genome)\nbat_genome_aligned = ClustalwCommandline('clustalw2', infile=bat_sl_cov_genome)\ncivet_genome_aligned = ClustalwCommandline('clustalw2', infile=civet_sl_cov_genome)\ncamel_genome_aligned = ClustalwCommandline('clustalw2', infile=camel_cov_genome)\nhedgehog_genome_aligned = ClustalwCommandline('clustalw2', infile=hedgehog_cov_genome)\nebola_genome_aligned = ClustalwCommandline('clustalw2', infile=ebola5_genome)\nmalaria_genome_aligned = ClustalwCommandline('clustalw2', infile=malaria_genome)\nhiv_genome_aligned = ClustalwCommandline('clustalw2', infile=hiv2_genome)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.099841Z","iopub.execute_input":"2021-11-27T06:29:24.100253Z","iopub.status.idle":"2021-11-27T06:29:24.117421Z","shell.execute_reply.started":"2021-11-27T06:29:24.100208Z","shell.execute_reply":"2021-11-27T06:29:24.116593Z"},"trusted":true},"execution_count":20,"outputs":[]},{"cell_type":"markdown","source":"as the genomes sequences are too big in length, lets align and have a look at just the first few nucleodtides in the sequences","metadata":{}},{"cell_type":"code","source":"# comparing first 57 nucleotides' of all the genomes after aligning them\nprint(\"SARS:\\t \",end=\" \")\nseq_repr(sars_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nMERS:\\t \",end=\" \")\nseq_repr(mers_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nCOVID:\\t \",end=\" \")\nseq_repr(cov_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nBAT:\\t \",end=\" \")\nseq_repr(bat_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nCIVET:\\t \",end=\" \")\nseq_repr(civet_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nCAMEL:\\t \",end=\" \")\nseq_repr(camel_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nHEDGEHOG:\",end=\" \")\nseq_repr(hedgehog_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nEBOLA:\\t \",end=\" \")\nseq_repr(ebola_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nMALARIA: \",end=\" \")\nseq_repr(malaria_genome_aligned.infile[:57])\nprint(Style.RESET_ALL + \"\\nHIV:\\t \",end=\" \")\nseq_repr(hiv_genome_aligned.infile[:57])","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.118931Z","iopub.execute_input":"2021-11-27T06:29:24.119349Z","iopub.status.idle":"2021-11-27T06:29:24.183188Z","shell.execute_reply.started":"2021-11-27T06:29:24.119306Z","shell.execute_reply":"2021-11-27T06:29:24.182452Z"},"trusted":true},"execution_count":21,"outputs":[]},{"cell_type":"markdown","source":"let us also have a look at the first few **aligned protein sequences** of these genomes","metadata":{}},{"cell_type":"code","source":"sars_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(sars_genome))\nmers_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(mers_genome))\ncov_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(cov_genome))\nbat_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(bat_sl_cov_genome))\ncivet_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(civet_sl_cov_genome))\ncamel_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(camel_cov_genome))\nhedgehog_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(hedgehog_cov_genome))\nebola_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(ebola5_genome))\nmalaria_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(malaria_genome))\nhiv_pro_aligned = ClustalwCommandline('clustalw2', infile=gen_protein_seq(hiv2_genome))\n\n# comparing first 50 nucleotides' of all the protein seq after aligning them\nprint(\"SARS Proteins:\\t  \",end=\" \")\nseq_repr(sars_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nMERS Proteins: \\t  \",end=\" \")\nseq_repr(mers_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nCOVID Proteins:   \",end=\" \")\nseq_repr(cov_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nBAT Proteins: \\t  \",end=\" \")\nseq_repr(bat_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nCIVET Proteins:   \",end=\" \")\nseq_repr(civet_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nCAMEL Proteins:   \",end=\" \")\nseq_repr(camel_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nHEDGEHOG Proteins:\",end=\" \")\nseq_repr(hedgehog_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nEBOLA Proteins:\\t  \",end=\" \")\nseq_repr(ebola_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nMALARIA Proteins: \",end=\" \")\nseq_repr(malaria_pro_aligned.infile[:50], strand='protein')\nprint(Style.RESET_ALL + \"\\nHIV Proteins: \\t  \",end=\" \")\nseq_repr(hiv_pro_aligned.infile[:50], strand='protein')","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.186969Z","iopub.execute_input":"2021-11-27T06:29:24.187294Z","iopub.status.idle":"2021-11-27T06:29:24.300611Z","shell.execute_reply.started":"2021-11-27T06:29:24.187265Z","shell.execute_reply":"2021-11-27T06:29:24.299550Z"},"trusted":true},"execution_count":22,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Pairwise-Sequence-Alignment\"></a>\n# Pairwise Sequence Alignment\nearlier we had aligned every sequence based on the [Clustal](https://en.wikipedia.org/wiki/Clustal) algorithm, now let's align every other genome with the covid19 genomes.","metadata":{}},{"cell_type":"code","source":"print(\"Aligning the first 100 nucleotides of COVID-19 genome with every other genomes:\\n\\n\")\nfor genome_name, genome in genomes.items():\n    if genome_name != 'COVID-19':\n        print(\"** COVID-19 and {} genome **\".format(genome_name))\n        alignments = pairwise2.align.globalxx(cov_genome[:100], genome[:100])\n        print(pairwise2.format_alignment(*alignments[0], full_sequences=True))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.302202Z","iopub.execute_input":"2021-11-27T06:29:24.302466Z","iopub.status.idle":"2021-11-27T06:29:24.527155Z","shell.execute_reply.started":"2021-11-27T06:29:24.302431Z","shell.execute_reply":"2021-11-27T06:29:24.526369Z"},"trusted":true},"execution_count":23,"outputs":[]},{"cell_type":"markdown","source":"Let us also align the Protein Sequences of these genomes. This sort of analysis is often used by Researchers to spot **Similarity in Protein Spikes**\n\n**[NOTE]** in order to avoid the confussion of misinterpretting the encoding of 'STOP' codon with the encoding for S, T and P amino acids, the STOP codons in the following protein seq alignment has been replaces by $'*'$ symbol!","metadata":{}},{"cell_type":"code","source":"print(\"Aligning the first 100 amino acids of COVID-19 genome's protein sequence with every other genomes's protein sequence:\\n\\n\")\n\nfor genome_name, protein_seq in proteins_seq.items():\n    if genome_name != 'COVID-19':\n        print(\"COVID-19 and {} genomes' protein sequences\\n\".format(genome_name))\n        alignments = pairwise2.align.globalxx(gen_protein_seq(cov_genome[:100]).replace('STOP','*'), protein_seq[:100].replace('STOP','*'))\n        print(pairwise2.format_alignment(*alignments[0], full_sequences=True))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.528469Z","iopub.execute_input":"2021-11-27T06:29:24.528714Z","iopub.status.idle":"2021-11-27T06:29:24.692359Z","shell.execute_reply.started":"2021-11-27T06:29:24.528687Z","shell.execute_reply":"2021-11-27T06:29:24.691538Z"},"trusted":true},"execution_count":24,"outputs":[]},{"cell_type":"markdown","source":"the `score` in all the above alignment is the similarity score base on the number of matches in nucleotides and the dashes that showed up upon alignment are what is known as `gaps`. It is these gaps that account for the genetic mutations occuring from the insertions or deletions in the sequence, sometimes also referred to as **indels**","metadata":{}},{"cell_type":"markdown","source":"<a id=\"Genome-Sequence-Similarity\"></a>\n# Genome Sequence Similarity\n**using Clustal algorithm**<br>\n\nfor establishing the similarity among genomes, we'll be changing the scoring parameter<br>\nfor every match there will be a $+1$ point<br>\nfor every mismatch there will be a $-1$ point<br>\nfor every open gap there will be a $-0.5$ points<br>\nfor every extended gap there will be a $-0.1$ point<br>","metadata":{}},{"cell_type":"code","source":"aligner = Align.PairwiseAligner()\naligner.mode = 'global'\naligner.match_score = 1\naligner.mismatch_score = -1\naligner.open_gap_score = -0.5\naligner.extend_gap_score = -0.1","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.693547Z","iopub.execute_input":"2021-11-27T06:29:24.693757Z","iopub.status.idle":"2021-11-27T06:29:24.703260Z","shell.execute_reply.started":"2021-11-27T06:29:24.693733Z","shell.execute_reply":"2021-11-27T06:29:24.702291Z"},"trusted":true},"execution_count":25,"outputs":[]},{"cell_type":"code","source":"print(\"Similarity scores between\\n\")\nfor genome_name, genome in genomes.items():\n    if genome_name != 'COVID-19':\n        score = aligner.score(cov_genome, genome)\n        print(\"COVID-19 & {} genome sequences:\\t {} ({:.2f}%)\".format(genome_name, score,100*(score/len(cov_genome))))","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:29:24.705060Z","iopub.execute_input":"2021-11-27T06:29:24.705466Z","iopub.status.idle":"2021-11-27T06:30:05.031626Z","shell.execute_reply.started":"2021-11-27T06:29:24.705413Z","shell.execute_reply":"2021-11-27T06:30:05.030690Z"},"trusted":true},"execution_count":26,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Generating-Phylogenetic-Tree\"></a>\n# Generating a Phylogenetic Tree\nA [Phylogenetic Tree](https://en.wikipedia.org/wiki/Phylogenetic_tree) is the tree of evolutionary relationship among various species based on their genetic characteristics.\n### UPGMA (Unweighted Pair Group Method with Arithmetic mean) algorithm\nThis is where you can learn more about it [wiki](https://en.wikipedia.org/wiki/UPGMA)<br>\n[NOTE] unlike the above sequence similarity, here only the mismatch cases are scored using the aligner and no scoring is done for gaps and extended gaps.","metadata":{}},{"cell_type":"code","source":"# set aligner to score on mismatches\naligner.match_score = 0\naligner.mismatch_score = 1\naligner.open_gap_score = 0\naligner.extend_gap_score = 0\n\n# code UPGMA by hand\ndef low_value_position(matrix):\n    # cell value with min value\n    minimum_value = float('inf')\n    position_x = -1\n    position_y = -1\n\n    # traverse the entire matrix  to find the smallest value\n    for i in range(len(matrix)):\n        for j in range(len(matrix[i])):\n            if (matrix[i][j] < minimum_value):\n                minimum_value = matrix[i][j]\n                position_x = i\n                position_y = j\n    return position_x, position_y\n\n# Combine two labels into one and erase the position with less value\n# params: List of labels, positions\ndef join_headers(headers, a, b):\n    if b < a:\n        a,b = b,a\n\n    # join labels in the first position\n    headers[a] = \"(\" + headers[a] + \",\" + headers[b] + \")\"\n\n    # delete label from secnd position\n    del headers[b]\n\n# Join two columns of the matrix and save the result (average) in the first column\n# Params: matrix of differences, position of the columns\ndef join_matrix(matrix,a ,b):\n    if b < a:\n        a,b = b,a\n\n    # recalculate row for smaller index\n    row = []\n    for i in range(0,a):\n        row.append((matrix[a][i] + matrix[b][i])/2)\n    matrix[a] = row\n\n    # recalculate entire column\n    for i in range(a+1, b):\n        matrix[i][a] = (matrix[i][a] + matrix[b][i])/2\n        \n    # get the rest of the values in row i\n    for i in range(b+1, len(matrix)):\n        matrix[i][a] = (matrix[i][a] + matrix[i][b])/2\n        # delete the second column\n        del matrix[i][b]\n\n    # delete second row\n    del matrix[b]\n\ndef upgma(matrix, headers):\n\n    iteration = 0\n    while len(headers) > 1:\n#         print('\\niteration: ', iteration)\n#         print('total clusters: ', len(matrix))\n\n        # locate the position in the matrix containing the lowest value\n        x, y = low_value_position(matrix)\n\n        # update the headers\n        join_headers(headers, x, y)\n#         print('family grouping: ', headers)\n\n        # join two columns of the difference matrix in the position coordinates\n        join_matrix(matrix,x,y)\n        iteration += 1\n\n    # final result stored in the first position\n    return headers[0]","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:30:05.033219Z","iopub.execute_input":"2021-11-27T06:30:05.033568Z","iopub.status.idle":"2021-11-27T06:30:05.052678Z","shell.execute_reply.started":"2021-11-27T06:30:05.033532Z","shell.execute_reply":"2021-11-27T06:30:05.051694Z"},"trusted":true},"execution_count":27,"outputs":[]},{"cell_type":"code","source":"headers = [\"COVID19\",\"SARS\",\"MERS\",\"Bat with SARS like CoV\",\"Civet with SARS like CoV\",\"Camel with CoV\",\"Hedgehog with CoV\",\"EBOLA\",\"Malaria\",\"HIV\"]\nviruses = [cov_genome, sars_genome, mers_genome, bat_sl_cov_genome, civet_sl_cov_genome, camel_cov_genome, hedgehog_cov_genome, ebola5_genome, malaria_genome, hiv2_genome]\n\ndifference_matrix = [\n    [],\n    [aligner.score(viruses[1], viruses[0])],# mers-sars\n    [aligner.score(viruses[2], viruses[0]),aligner.score(viruses[2], viruses[1])],# cov-sars, cov-mers\n    [aligner.score(viruses[3], viruses[0]),aligner.score(viruses[3], viruses[1]), aligner.score(viruses[3], viruses[2])],# civet-sars, civet-mers, civet-cov\n    [aligner.score(viruses[4], viruses[0]),aligner.score(viruses[4], viruses[1]), aligner.score(viruses[4], viruses[2]), aligner.score(viruses[4], viruses[3])],# bat-sars, bat-mers, bat-cov, bat-civet\n    [aligner.score(viruses[5], viruses[0]),aligner.score(viruses[5], viruses[1]), aligner.score(viruses[5], viruses[2]), aligner.score(viruses[5], viruses[3]), aligner.score(viruses[5], viruses[4])],\n    [aligner.score(viruses[6], viruses[0]),aligner.score(viruses[6], viruses[1]), aligner.score(viruses[6], viruses[2]), aligner.score(viruses[6], viruses[3]), aligner.score(viruses[6], viruses[4]), aligner.score(viruses[6], viruses[5])],\n    [aligner.score(viruses[7], viruses[0]),aligner.score(viruses[7], viruses[1]), aligner.score(viruses[7], viruses[2]), aligner.score(viruses[7], viruses[3]), aligner.score(viruses[7], viruses[4]), aligner.score(viruses[7], viruses[5]), aligner.score(viruses[7], viruses[6])],\n    [aligner.score(viruses[8], viruses[0]),aligner.score(viruses[8], viruses[1]), aligner.score(viruses[8], viruses[2]), aligner.score(viruses[8], viruses[3]), aligner.score(viruses[8], viruses[4]), aligner.score(viruses[8], viruses[5]), aligner.score(viruses[8], viruses[6]), aligner.score(viruses[8], viruses[7])],\n    [aligner.score(viruses[9], viruses[0]),aligner.score(viruses[9], viruses[1]), aligner.score(viruses[9], viruses[2]), aligner.score(viruses[9], viruses[3]), aligner.score(viruses[9], viruses[4]), aligner.score(viruses[9], viruses[5]), aligner.score(viruses[9], viruses[6]), aligner.score(viruses[9], viruses[7]), aligner.score(viruses[9],viruses[8])]\n]\n\nphylogenetic_tree = upgma(difference_matrix, headers)\nphylogenetic_tree","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:30:05.053937Z","iopub.execute_input":"2021-11-27T06:30:05.054206Z","iopub.status.idle":"2021-11-27T06:32:17.405796Z","shell.execute_reply.started":"2021-11-27T06:30:05.054179Z","shell.execute_reply":"2021-11-27T06:32:17.405083Z"},"trusted":true},"execution_count":28,"outputs":[]},{"cell_type":"code","source":"tree = Tree(phylogenetic_tree + \";\")\n# render tree image\ntree.render(\"%%inline\", w=500)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:32:17.410367Z","iopub.execute_input":"2021-11-27T06:32:17.410627Z","iopub.status.idle":"2021-11-27T06:32:17.812987Z","shell.execute_reply.started":"2021-11-27T06:32:17.410600Z","shell.execute_reply":"2021-11-27T06:32:17.812218Z"},"trusted":true},"execution_count":29,"outputs":[]},{"cell_type":"code","source":"# show the phylogenetic family as circular tree\nts = TreeStyle()\nts.show_leaf_name = True\nts.mode = 'c'\nts.arc_start = 10 # 0 degrees = 3 o'clock\nts.arc_span = 180\n\ntree.render(\"%%inline\", w=800, tree_style=ts)","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-11-27T06:32:17.814569Z","iopub.execute_input":"2021-11-27T06:32:17.814826Z","iopub.status.idle":"2021-11-27T06:32:17.895918Z","shell.execute_reply.started":"2021-11-27T06:32:17.814798Z","shell.execute_reply":"2021-11-27T06:32:17.895243Z"},"trusted":true},"execution_count":30,"outputs":[]},{"cell_type":"markdown","source":"<a id=\"Conclusions\"></a>\n# Conclusions\n\n- Based on the oligonucleotide compositions, all the genomes from coronaviridae family show similar composition (trinucleotide and tetranucleotide composition)\n- Based on GC content, the composition of Hedgehog infected Corona virus is closely similar to COVID-19 genome's GC content.\n- Judging by the number of protein strands, the count in COVID19's genome is very close to those found in Bat infected with SARS like corona virus.\n- If you filter down some of the categories in the plotly-plot of distribution of amino acids, you'll see the percentage of amino acids in COVID-19's genome closely match to those of MERS, Civet infected with SARS like CoV, Camel with alpha CoV and Malaria except for some and low high contents among L (Leucine), P (Proline) and T (Threonine) amino acids.\n- Checking by the generated Phylogenetic tree, the genomes of COVID-19 and malaria virus show the maximum similarity. Perhaps this also explains why there's a lot of buzz on the news and internet about trying [Hydroxycholorquine](https://en.wikipedia.org/wiki/Hydroxychloroquine) and [anti-retroviral drugs](https://www.deccanherald.com/city/top-bengaluru-stories/coronavirus-hiv-anti-retroviral-drugs-to-be-used-for-covid-19-treatment-814534.html) which happen to be the medicines for malaria and HIV viruses, respectively.\n- Of all the genomes of COVID-19 infected patients avaialable on Kaggle, some patients show similar traits when clustered based on their genome sequence similarity.I really don't know how helpful this would be but if something works or had worked for a patient, then same should also work for patients in the same cluster.","metadata":{}},{"cell_type":"markdown","source":"I really hope this turns out to be of any help to any microbiologist or vaccinologist out there.\nIf you liked this kernel or found it useful, please consider upvoting this notebook.<br>Thank you. Stay Home & Stay Safe!","metadata":{}}]}