{ "cells": [ { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [] }, "source": [ "# sequences and biopython\n", "## 09/21/2023\n", "\n", "print view
\n", "notebook" ] }, { "cell_type": "code", "execution_count": 30, "metadata": { "slideshow": { "slide_type": "-" }, "tags": [] }, "outputs": [ { "data": { "text/html": [ "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%html\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Sequence data\n", "\n", "[http://www.ncbi.nlm.nih.gov](http://www.ncbi.nlm.nih.gov)\n", "\n", "[Example](http://www.ncbi.nlm.nih.gov/gene/5216)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# FASTA\n", "\n", "First line is description of sequence and starts with `>`\n", "\n", "All lines up to the next `>` are part of the same sequence\n", "\n", "Usually less than 80 characters per line\n", "\n", "```\n", ">gi|568815581:c4949086-4945650 Homo sapiens chromosome 17, GRCh38.p2 Primary Assembly\n", "CCCGCAGGGTCCACACGGGTCGGGCCGGGCGCGCTCCCGTGCAGCCGGCTCCGGCCCCGACCGCCCCATG\n", "CACTCCCGGCCCCGGCGCAGGCGCAGGCGCGGGCACACGCGCCGCCGCCCGCCGGTCCTTCCCTTCGGCG\n", "GAGGTGGGGGAAGGAGGAGTCATCCCGTTTAACCCTGGGCTCCCCGAACTCTCCTTAATTTGCTAAATTT\n", "GCAGCTTGCTAATTCCTCCTGCTTTCTCCTTCCTTCCTTCTTCTGGCTCACTCCCTGCCCCGATACCAAA\n", "GTCTGGTTTATATTCAGTGCAAATTGGAGCAAACCCTACCCTTCACCTCTCTCCCGCCACCCCCCATCCT\n", "TCTGCATTGCTTTCCATCGAACTCTGCAAATTTTGCAATAGGGGGAGGGATTTTTAAAATTGCATTTGCA\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Genbank\n", "\n", "Annotated format. Starts with `LOCUS` field. Can have several other annotation (e.g. `KEYWORDS`, `SOURCE`, `REFERENCE`, `FEATURES`).\n", "\n", "`ORIGIN` record indicates start of sequence\n", "\n", "'`\\\\`' indicates the end of sequence\n", "\n", "```\n", "LOCUS CAG28598 140 aa linear PRI 16-OCT-2008\n", "DEFINITION PFN1, partial [Homo sapiens].\n", "ACCESSION CAG28598\n", "VERSION CAG28598.1 GI:47115277\n", "DBSOURCE embl accession CR407670.1\n", "KEYWORDS .\n", "SOURCE Homo sapiens (human)\n", " ORGANISM Homo sapiens\n", " Eukaryota; Metazoa; Chordata; Craniata; Vertebrata; Euteleostomi;\n", " Mammalia; Eutheria; Euarchontoglires; Primates; Haplorrhini;\n", " Catarrhini; Hominidae; Homo.\n", "ORIGIN \n", " 1 magwnayidn lmadgtcqda aivgykdsps vwaavpgktf vnitpaevgv lvgkdrssfy\n", " 61 vngltlggqk csvirdsllq dgefsmdlrt kstggaptfn vtvtktdktl vllmgkegvh\n", " 121 gglinkkcye mashlrrsqy\n", "//\n", "```" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Biopython\n", "\n", "*Biopython features include parsers for various Bioinformatics file formats (BLAST, Clustalw, FASTA, Genbank,...), access to online services (NCBI, Expasy,...), interfaces to common and not-so-common programs (Clustalw, DSSP, MSMS...), a standard sequence class, various clustering modules, a KD tree data structure etc. and even documentation.*\n", "\n", "Other modules that might be of interest:\n", "\n", " * Pycogent: http://pycogent.org/\n", " * bx-python: http://bitbucket.org/james_taylor/bx-python/wiki/Home\n", " * DendroPy: http://packages.python.org/DendroPy/\n", " * Pygr: http://code.google.com/p/pygr/\n", " * bioservices: https://bioservices.readthedocs.io/en/master/\n", " \n", "Biopython is **not** for performing sequencing itself (see: https://crc.pitt.edu/training/fall-2021-next-generation-sequencing-workshops)." ] }, { "cell_type": "markdown", "metadata": { "collapsed": true, "jupyter": { "outputs_hidden": true }, "slideshow": { "slide_type": "slide" } }, "source": [ "# Sequence Objects" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('GATTACA')" ] }, "execution_count": 31, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from Bio.Seq import Seq # the submodule structure of biopython is a little awkward\n", "\n", "s = Seq(\"GATTACA\")\n", "s" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Sequences act a lot like strings, but have additional methods.\n", "\n", "Methods shared with `str`: `count`, `endswith`, `find`, `lower`, `lstrip`, `rfind`, `split`, `startswith`, `strip`,`upper`\n", "\n", "`Seq` methods:`back_transcribe`, `complement`, `reverse_complement`, `tomutable`, `tostring`, `transcribe`, `translate`, `ungap`" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Accessing Seq data\n", "\n", "Sequences act like strings (indexed from 0)" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'G'" ] }, "execution_count": 32, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s[0]" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('TT')" ] }, "execution_count": 33, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s[2:4] #returns sequence" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('gattaca')" ] }, "execution_count": 34, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s.lower()" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('GATTACAGATTACA')" ] }, "execution_count": 35, "metadata": {}, "output_type": "execute_result" } ], "source": [ "s + s" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# The Central Dogma" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```\n", "DNA coding strand (aka Crick strand, strand +1)\t \n", "5’\tATGGCCATTGTAATGGGCCGCTGAAAGGGTGCCCGATAG\t3’\n", " \t|||||||||||||||||||||||||||||||||||||||\t \n", "3’\tTACCGGTAACATTACCCGGCGACTTTCCCACGGGCTATC\t5’\n", " \tDNA template strand (aka Watson strand, strand −1)\t \n", " |\t \n", " Transcription\t \n", " ↓\t \n", " \n", "5’\tAUGGCCAUUGUAAUGGGCCGCUGAAAGGGUGCCCGAUAG\t3’\n", " \tSingle stranded messenger RNA\t \n", " |\t \n", " Translation\t \n", " ↓\t \n", " MAIVMGR*KGAR*\n", " amino acid sequence (w/stop codons)\n", " ```" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "(Seq('CTAATGTCTAATGTCTAATGT'), Seq('TGTAATCTGTAATCTGTAATC'))" ] }, "execution_count": 36, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dna = Seq('GATTACAGATTACAGATTACA')\n", "dna.complement(),dna.reverse_complement()" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# The Central Dogma" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('GATTACAGATTACAGATTACA')" ] }, "execution_count": 37, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dna" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('GAUUACAGAUUACAGAUUACA', RNAAlphabet())" ] }, "execution_count": 38, "metadata": {}, "output_type": "execute_result" } ], "source": [ "rna = dna.transcribe()\n", "rna" ] }, { "cell_type": "code", "execution_count": 39, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('DYRLQIT', ExtendedIUPACProtein())" ] }, "execution_count": 39, "metadata": {}, "output_type": "execute_result" } ], "source": [ "protein = rna.translate()\n", "protein" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "Seq('DYRLQIT', ExtendedIUPACProtein())" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "dna.translate() #unlike cells, don't actually need rna" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Codon Tables\n", "\n", "By default the standard translation table is used, but others can be provided to the translate method." ] }, { "cell_type": "code", "execution_count": 41, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['Alternative Flatworm Mitochondrial', 'Alternative Yeast Nuclear', 'Archaeal', 'Ascidian Mitochondrial', 'Bacterial', 'Balanophoraceae Plastid', 'Blastocrithidia Nuclear', 'Blepharisma Macronuclear', 'Candidate Division SR1', 'Chlorophycean Mitochondrial', 'Ciliate Nuclear', 'Coelenterate Mitochondrial', 'Condylostoma Nuclear', 'Dasycladacean Nuclear', 'Echinoderm Mitochondrial', 'Euplotid Nuclear', 'Flatworm Mitochondrial', 'Gracilibacteria', 'Hexamita Nuclear', 'Invertebrate Mitochondrial', 'Karyorelict Nuclear', 'Mesodinium Nuclear', 'Mold Mitochondrial', 'Mycoplasma', 'Pachysolen tannophilus Nuclear', 'Peritrich Nuclear', 'Plant Plastid', 'Protozoan Mitochondrial', 'Pterobranchia Mitochondrial', 'SGC0', 'SGC1', 'SGC2', 'SGC3', 'SGC4', 'SGC5', 'SGC8', 'SGC9', 'Scenedesmus obliquus Mitochondrial', 'Spiroplasma', 'Standard', 'Thraustochytrium Mitochondrial', 'Trematode Mitochondrial', 'Vertebrate Mitochondrial', 'Yeast Mitochondrial']\n" ] } ], "source": [ "from Bio.Data import CodonTable\n", "print(sorted(CodonTable.unambiguous_dna_by_name.keys()))" ] }, { "cell_type": "code", "execution_count": 42, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Table 1 Standard, SGC0\n", "\n", " | T | C | A | G |\n", "--+---------+---------+---------+---------+--\n", "T | TTT F | TCT S | TAT Y | TGT C | T\n", "T | TTC F | TCC S | TAC Y | TGC C | C\n", "T | TTA L | TCA S | TAA Stop| TGA Stop| A\n", "T | TTG L(s)| TCG S | TAG Stop| TGG W | G\n", "--+---------+---------+---------+---------+--\n", "C | CTT L | CCT P | CAT H | CGT R | T\n", "C | CTC L | CCC P | CAC H | CGC R | C\n", "C | CTA L | CCA P | CAA Q | CGA R | A\n", "C | CTG L(s)| CCG P | CAG Q | CGG R | G\n", "--+---------+---------+---------+---------+--\n", "A | ATT I | ACT T | AAT N | AGT S | T\n", "A | ATC I | ACC T | AAC N | AGC S | C\n", "A | ATA I | ACA T | AAA K | AGA R | A\n", "A | ATG M(s)| ACG T | AAG K | AGG R | G\n", "--+---------+---------+---------+---------+--\n", "G | GTT V | GCT A | GAT D | GGT G | T\n", "G | GTC V | GCC A | GAC D | GGC G | C\n", "G | GTA V | GCA A | GAA E | GGA G | A\n", "G | GTG V | GCG A | GAG E | GGG G | G\n", "--+---------+---------+---------+---------+--\n" ] } ], "source": [ "print(CodonTable.unambiguous_dna_by_name['Standard'])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "slideshow": { "slide_type": "slide" }, "tags": [] }, "outputs": [ { "data": { "text/html": [ "
\n", "\n" ], "text/plain": [ "" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "%%html\n", "
\n", "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# SeqRecord\n", "\n", "Sequence data is read/written as `SeqRecord` objects.\n", "\n", "These objects store additional information about the sequence (name, id, description, features)\n", "\n", "`SeqIO` reads sequence records: \n", " * must specify format\n", " * `read` to read a file with a single record\n", " * `parse` to iterate over file with mulitple records" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "SeqRecord(seq=Seq('GATGGGATTGGGGTTTTCCCCTCCCATGTGCTCAAGACTGGCGCTAAAAGTTTT...GTG', IUPACAmbiguousDNA()), id='NG_017013.2', name='NG_017013', description='Homo sapiens tumor protein p53 (TP53), RefSeqGene (LRG_321) on chromosome 17', dbxrefs=[])" ] }, "execution_count": 44, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from Bio.SeqRecord import SeqRecord\n", "from Bio import SeqIO\n", "seq = SeqIO.read('../files/p53.gb','genbank')\n", "seq" ] }, { "cell_type": "code", "execution_count": 45, "metadata": {}, "outputs": [], "source": [ "seqs = []\n", "# https://MSCBIO2025.github.io/files/hydra.fasta\n", "for s in SeqIO.parse('../files/hydra.fasta','fasta'):\n", " seqs.append(s)" ] }, { "cell_type": "code", "execution_count": 46, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "40" ] }, "execution_count": 46, "metadata": {}, "output_type": "execute_result" } ], "source": [ "len(seqs)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Fetching sequences from the Internet\n", "\n", "Biopython provides and interface to the [NCBI \"Entrez\" search engine](http://www.ncbi.nlm.nih.gov/sites/gquery)\n", "\n", "The results of internet queries are returned as file-like objects." ] }, { "cell_type": "code", "execution_count": 47, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'DbList': ['pubmed', 'protein', 'nuccore', 'ipg', 'nucleotide', 'structure', 'genome', 'annotinfo', 'assembly', 'bioproject', 'biosample', 'blastdbinfo', 'books', 'cdd', 'clinvar', 'gap', 'gapplus', 'grasp', 'dbvar', 'gene', 'gds', 'geoprofiles', 'homologene', 'medgen', 'mesh', 'nlmcatalog', 'omim', 'orgtrack', 'pmc', 'popset', 'proteinclusters', 'pcassay', 'protfam', 'pccompound', 'pcsubstance', 'seqannot', 'snp', 'sra', 'taxonomy', 'biocollections', 'gtr']}" ] }, "execution_count": 47, "metadata": {}, "output_type": "execute_result" } ], "source": [ "from Bio import Entrez\n", "Entrez.email = \"jpbarton@pitt.edu\" # biopython forces you to provide your email\n", "res = Entrez.read(Entrez.einfo()) # the names of all available databases\n", "res" ] }, { "cell_type": "code", "execution_count": 48, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "['annotinfo', 'assembly', 'biocollections', 'bioproject', 'biosample', 'blastdbinfo', 'books', 'cdd', 'clinvar', 'dbvar', 'gap', 'gapplus', 'gds', 'gene', 'genome', 'geoprofiles', 'grasp', 'gtr', 'homologene', 'ipg', 'medgen', 'mesh', 'nlmcatalog', 'nuccore', 'nucleotide', 'omim', 'orgtrack', 'pcassay', 'pccompound', 'pcsubstance', 'pmc', 'popset', 'protein', 'proteinclusters', 'protfam', 'pubmed', 'seqannot', 'snp', 'sra', 'structure', 'taxonomy']\n" ] } ], "source": [ "print(sorted(res['DbList']))" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# ESearch\n", "\n", "You can search any database for a given term and it will return the ids of all the relevant records." ] }, { "cell_type": "code", "execution_count": 49, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "{'Count': '16953', 'RetMax': '20', 'RetStart': '0', 'IdList': ['2579134893', '2577097206', '1675046266', '1621312268', '1621310586', '1621309831', '1565520846', '1565520832', '1519314995', '1519314531', '1519245996', '1519245289', '1492459888', '2544609391', '2544608770', '2544608358', '2544608205', '1732746165', '2578505094', '2578491630'], 'TranslationSet': [], 'TranslationStack': [{'Term': 'tp53[All Fields]', 'Field': 'All Fields', 'Count': '16953', 'Explode': 'N'}, 'GROUP'], 'QueryTranslation': 'tp53[All Fields]'}\n" ] } ], "source": [ "result = Entrez.esearch(db='nucleotide', term='tp53') # the result is a file-like object of the raw xml data\n", "records = Entrez.read(result) # put into a more palatable form (dictionary)\n", "print(records)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There were 2923 7616 8238 10725 11850 12751 13952 16943 hits, but by default only 20 are returned. We can change this ([and other parameters](http://www.ncbi.nlm.nih.gov/books/NBK25499/#chapter4.ESearch)) by changing our search terms." ] }, { "cell_type": "code", "execution_count": 50, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "data": { "text/plain": [ "{'Count': '16953', 'RetMax': '50', 'RetStart': '0', 'IdList': ['2579134893', '2577097206', '1675046266', '1621312268', '1621310586', '1621309831', '1565520846', '1565520832', '1519314995', '1519314531', '1519245996', '1519245289', '1492459888', '2544609391', '2544608770', '2544608358', '2544608205', '1732746165', '2578505094', '2578491630', '2544608591', '1677501386', '1677500090', '1653961881', '349585153', '50344732', '2577326767', '2577319026', '2577314085', '2577314083', '2577314081', '2577144630', '2577144628', '2577144626', '2577144624', '2577144622', '2577144620', '2577144618', '2577112865', '2577109742', '2577244898', '2577244896', '2577244894', '2577238979', '2577226876', '2577220539', '2577220537', '2577179306', '2577179304', '2577176327'], 'TranslationSet': [], 'TranslationStack': [{'Term': 'tp53[All Fields]', 'Field': 'All Fields', 'Count': '16953', 'Explode': 'N'}, 'GROUP'], 'QueryTranslation': 'tp53[All Fields]'}" ] }, "execution_count": 50, "metadata": {}, "output_type": "execute_result" } ], "source": [ "records = Entrez.read(Entrez.esearch(db='nucleotide', term='tp53', retmax=50))\n", "records" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# EFetch\n", "\n", "To get the full record for a given id, use `efetch`. \n", "\n", "Must provide `rettype` ([available options](http://www.ncbi.nlm.nih.gov/books/NBK25499/table/chapter4.T._valid_values_of__retmode_and/?report=objectonly) include fasta and gb)\n", "\n", "`retmode` can be text or xml." ] }, { "cell_type": "code", "execution_count": 51, "metadata": {}, "outputs": [], "source": [ "#fetch the genbank file for the first id from our search\n", "result = Entrez.efetch(db=\"nucleotide\",id=records['IdList'][0],rettype=\"gb\",retmode='text')\n", "#parse into a seqrecord\n", "p53 = SeqIO.read(result,'gb')" ] }, { "cell_type": "code", "execution_count": 52, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "<_io.TextIOWrapper encoding='latin-1'>" ] }, "execution_count": 52, "metadata": {}, "output_type": "execute_result" } ], "source": [ "result" ] }, { "cell_type": "code", "execution_count": 53, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "SeqRecord(seq=Seq('TCAGGGTCCAGGTGATTCTGAACGAGCTGAAAATCGAAGATAAACAACACTTTC...AAA', IUPACAmbiguousDNA()), id='XM_059624316.1', name='XM_059624316', description='PREDICTED: Neocloeon triangulifer TP53-regulated inhibitor of apoptosis 1-like (LOC132199519), mRNA', dbxrefs=['BioProject:PRJNA1017488'])" ] }, "execution_count": 53, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p53" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Features\n", "\n", "Genbank files are typically annotated with features, which refer to portions of the full sequence\n", "\n", "`SeqRecord` objects track these features and you can extract the corresponding subsequence\n", "\n", "*CDS* - coding sequence" ] }, { "cell_type": "code", "execution_count": 54, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "[SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(712), strand=1), type='source'),\n", " SeqFeature(FeatureLocation(ExactPosition(0), ExactPosition(712), strand=1), type='gene'),\n", " SeqFeature(FeatureLocation(ExactPosition(197), ExactPosition(437), strand=1), type='CDS'),\n", " SeqFeature(FeatureLocation(ExactPosition(711), ExactPosition(712), strand=1), type='polyA_site')]" ] }, "execution_count": 54, "metadata": {}, "output_type": "execute_result" } ], "source": [ "p53.features" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Extracting subsequences" ] }, { "cell_type": "code", "execution_count": 55, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "type: CDS\n", "location: [197:437](+)\n", "qualifiers:\n", " Key: codon_start, Value: ['1']\n", " Key: db_xref, Value: ['GeneID:132199519']\n", " Key: gene, Value: ['LOC132199519']\n", " Key: product, Value: ['TP53-regulated inhibitor of apoptosis 1-like']\n", " Key: protein_id, Value: ['XP_059480299.1']\n", " Key: translation, Value: ['MNSVGEECNDMKKTYDSCFNAWFSHKFLRGQHDDQMCADLFLKYQQCVKNAMKQQNISIEDVQHFQLGTEEEKKAPPKK']\n", "\n" ] } ], "source": [ "cdsfeature = p53.features[2]\n", "print(cdsfeature)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "The subsequence the feature refers to can be extracted from the original full sequence using the feature." ] }, { "cell_type": "code", "execution_count": 56, "metadata": { "slideshow": { "slide_type": "-" } }, "outputs": [ { "data": { "text/plain": [ "SeqRecord(seq=Seq('ATGAACAGTGTCGGCGAGGAATGCAACGACATGAAAAAAACTTATGATTCCTGC...TAG', IUPACAmbiguousDNA()), id='XM_059624316.1', name='XM_059624316', description='PREDICTED: Neocloeon triangulifer TP53-regulated inhibitor of apoptosis 1-like (LOC132199519), mRNA', dbxrefs=[])" ] }, "execution_count": 56, "metadata": {}, "output_type": "execute_result" } ], "source": [ "coding = cdsfeature.extract(p53) #pass the full record (p53) to the feature\n", "coding" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# BLAST!\n", "\n", "Biopython let's you use [NCBI's BLAST](http://blast.ncbi.nlm.nih.gov/Blast.cgi) to search for similar sequences with `qblast` which has three required arguments:\n", " * **program**: blastn, blastp, blastx, tblastn, tblastx\n", " * **database**: see website\n", " * **sequence**: a sequence object" ] }, { "cell_type": "code", "execution_count": 57, "metadata": {}, "outputs": [], "source": [ "from Bio.Blast import NCBIWWW\n", "result = NCBIWWW.qblast(\"blastn\",\"nt\",coding.seq,hitlist_size=250)\n", "#result is a file-like object with xml in it - it can take a while to get results" ] }, { "cell_type": "code", "execution_count": 58, "metadata": {}, "outputs": [], "source": [ "from Bio.Blast import NCBIXML #for parsing xmls\n", "blast_records = NCBIXML.read(result)" ] }, { "cell_type": "code", "execution_count": 59, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "147 147\n" ] } ], "source": [ "print(len(blast_records.alignments),len(blast_records.descriptions))" ] }, { "cell_type": "code", "execution_count": 60, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1\n" ] } ], "source": [ "alignment = blast_records.alignments[0]\n", "print(len(alignment.hsps))" ] }, { "cell_type": "code", "execution_count": 61, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "****Alignment****\n", "sequence: gi|2489302823|emb|OX463782.1| Siphlonurus alternatus genome assembly, chromosome: 5\n", "length: 40728708\n", "e value: 2.75518e-13\n", "ATGAACAGTGTCGGCGAGGAATGCAACGACATGAAAAAAACTTATGATTCCTGCTTCAACGCATGGTTCAGCCAC...\n", "||||| ||||| || || || || | || ||||| || ||||||||||| || | ||||| ||||| ...\n", "ATGAATAGTGTTGGAGAAGAGTGTACTGATATGAAGAAGCAATATGATTCCTGTTTTCATCTTTGGTTTAGCCAA...\n" ] } ], "source": [ "hsp = alignment.hsps[0] # high scoring segment pairs\n", "print('****Alignment****')\n", "print('sequence:', alignment.title)\n", "print('length:', alignment.length)\n", "print('e value:', hsp.expect)\n", "print(hsp.query[0:75] + '...') # what we searched with\n", "print(hsp.match[0:75] + '...')\n", "print(hsp.sbjct[0:75] + '...') # what we matched to" ] }, { "cell_type": "code", "execution_count": 62, "metadata": { "slideshow": { "slide_type": "slide" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "****Alignment****\n", "sequence: gi|2362810685|emb|OX371006.1| Yponomeuta malinellus genome assembly, chromosome: 12\n", "length: 19760461\n", "e value: 2.94431\n", "AGAATGCAATGAAGCAGCAAAACAT...\n", "|||||||||||||||||||||||||...\n", "AGAATGCAATGAAGCAGCAAAACAT...\n" ] } ], "source": [ "alignment = blast_records.alignments[-1] # get last alignment\n", "hsp = alignment.hsps[0]\n", "print('****Alignment****')\n", "print('sequence:', alignment.title)\n", "print('length:', alignment.length)\n", "print('e value:', hsp.expect)\n", "print(hsp.query[0:75] + '...') # what we searched with\n", "print(hsp.match[0:75] + '...')\n", "print(hsp.sbjct[0:75] + '...') # what we matched to" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Creating alignments\n", "\n", "Biopython uses `clustal` to perform multiple alignments\n", "\n", "The interface is very simple - you construct the `clustal` commandline and must read the result files" ] }, { "cell_type": "code", "execution_count": 63, "metadata": {}, "outputs": [], "source": [ "from Bio.Align.Applications import ClustalwCommandline\n", "cline = ClustalwCommandline('clustalw', infile='../files/hydra.fasta',outfile='alignment.aln')" ] }, { "cell_type": "code", "execution_count": 64, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "clustalw -infile=../files/hydra.fasta -outfile=alignment.aln\n" ] } ], "source": [ "print(cline)" ] }, { "cell_type": "code", "execution_count": 65, "metadata": {}, "outputs": [], "source": [ "output = cline() #this calls the above, may take a while # If it doesn't work, sudo apt install clustalw" ] }, { "cell_type": "code", "execution_count": 66, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "CLUSTAL 2.1 multiple sequence alignment\r\n", "\r\n", "\r\n", "gi|164609123|gb|ABY62783.1|/1- ------------MLNIAILCLSCYINYVSCLSLTSPAIIEEVVGRSVTIT\r\n", "model1_F.pdb ------------------------------LSLTSPAIIEEVVGRSVTIT\r\n", "gi|225423262|gb|ACN91137.1|/1- ------------MLNIAILSLSCYINYVSCLSLTSPAIIEEVVGRSVTIT\r\n", "gi|302171750|gb|ADK97776.1|/1- ------------MLNIAILSLSCYINYVSCLSLTSPAIIEEVVGRSVTIT\r\n", "gi|302171774|gb|ADK97788.1|/1- ------------MLNIAILCLSCYINYVSCLSLTSSAIIEEVVNRSVTIT\r\n", "gi|225423244|gb|ACN91128.1|/1- ------------MLNIAILCLSCYINYVSCLSLTSSAIIEEVVNRSVTIT\r\n", "gi|225423272|gb|ACN91142.1|/1- ------------MLNIAILLLSCYINYVSSLSLTSPAIIEEVVNRSVTIT\r\n" ] } ], "source": [ "!head alignment.aln" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Alignments\n", "\n", "`AlignIO` is used to read alignment files (must provide format)" ] }, { "cell_type": "code", "execution_count": 67, "metadata": {}, "outputs": [], "source": [ "from Bio import AlignIO\n", "align = AlignIO.read('alignment.aln','clustal')" ] }, { "cell_type": "code", "execution_count": 68, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "< instance (40 records of length 691, SingleLetterAlphabet()) at 1068300a0>" ] }, "execution_count": 68, "metadata": {}, "output_type": "execute_result" } ], "source": [ "align" ] }, { "cell_type": "code", "execution_count": 69, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SingleLetterAlphabet() alignment with 40 rows and 691 columns\n", "------------MLNIAILCLSCYINYVSCLSLTSPAIIEEVVG...QKK gi|164609123|gb|ABY62783.1|/1-\n", "------------------------------LSLTSPAIIEEVVG...--- model1_F.pdb\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIEEVVG...QKK gi|225423262|gb|ACN91137.1|/1-\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIEEVVG...QKK gi|302171750|gb|ADK97776.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSLTSSAIIEEVVN...QKK gi|302171774|gb|ADK97788.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSLTSSAIIEEVVN...QKK gi|225423244|gb|ACN91128.1|/1-\n", "------------MLNIAILLLSCYINYVSSLSLTSPAIIEEVVN...QKK gi|225423272|gb|ACN91142.1|/1-\n", "------------MLNIAILCLSYYVNYVSCLSLTSPAIIEEVVG...QKK gi|302171742|gb|ADK97772.1|/1-\n", "------------MLNIAILCLSYYVNYVSCLSLTSPAIIEEVVG...QKK gi|225423264|gb|ACN91138.1|/1-\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIEEVVG...QKK gi|302171780|gb|ADK97791.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSLTSPAIIEEVVG...QKK gi|302171756|gb|ADK97779.1|/1-\n", "------------MLNVVILCLSCYINYVSCLSLTSPAIIEEVVG...QKK gi|225423288|gb|ACN91150.1|/1-\n", "MKLVIVIPNDENLLNIAILCLSCYLNYVSCLSLTSPAIIEKVVG...QKK gi|225423270|gb|ACN91141.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSVTSPPIIEEVVG...QKK gi|302171766|gb|ADK97784.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSVTSPPIIEEVVG...QKK gi|225423256|gb|ACN91134.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSVTSPPIIEEVVG...QKK gi|302171748|gb|ADK97775.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSLTSSAIIEEVVG...QKK gi|302171772|gb|ADK97787.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSVTSPAIIEEVVG...QKK gi|225423294|gb|ACN91153.1|/1-\n", "...\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIEKVVG...QKK gi|302171764|gb|ADK97783.1|/1-\n" ] } ], "source": [ "print(align)" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Slicing Alignments\n", "\n", "Alignments are sliced just like `numpy` arrays" ] }, { "cell_type": "code", "execution_count": 70, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "SeqRecord(seq=Seq('------------MLNIAILCLSCYINYVSCLSLTSPAIIEEVVGRSVTITYVTD...QKK', SingleLetterAlphabet()), id='gi|164609123|gb|ABY62783.1|/1-', name='', description='gi|164609123|gb|ABY62783.1|/1-', dbxrefs=[])" ] }, "execution_count": 70, "metadata": {}, "output_type": "execute_result" } ], "source": [ "align[0] #first row" ] }, { "cell_type": "code", "execution_count": 71, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'------------M-----------------M---------'" ] }, "execution_count": 71, "metadata": {}, "output_type": "execute_result" } ], "source": [ "align[:,0] #first column" ] }, { "cell_type": "code", "execution_count": 72, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "SingleLetterAlphabet() alignment with 40 rows and 40 columns\n", "------------MLNIAILCLSCYINYVSCLSLTSPAIIE gi|164609123|gb|ABY62783.1|/1-\n", "------------------------------LSLTSPAIIE model1_F.pdb\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIE gi|225423262|gb|ACN91137.1|/1-\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIE gi|302171750|gb|ADK97776.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSLTSSAIIE gi|302171774|gb|ADK97788.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSLTSSAIIE gi|225423244|gb|ACN91128.1|/1-\n", "------------MLNIAILLLSCYINYVSSLSLTSPAIIE gi|225423272|gb|ACN91142.1|/1-\n", "------------MLNIAILCLSYYVNYVSCLSLTSPAIIE gi|302171742|gb|ADK97772.1|/1-\n", "------------MLNIAILCLSYYVNYVSCLSLTSPAIIE gi|225423264|gb|ACN91138.1|/1-\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIE gi|302171780|gb|ADK97791.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSLTSPAIIE gi|302171756|gb|ADK97779.1|/1-\n", "------------MLNVVILCLSCYINYVSCLSLTSPAIIE gi|225423288|gb|ACN91150.1|/1-\n", "MKLVIVIPNDENLLNIAILCLSCYLNYVSCLSLTSPAIIE gi|225423270|gb|ACN91141.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSVTSPPIIE gi|302171766|gb|ADK97784.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSVTSPPIIE gi|225423256|gb|ACN91134.1|/1-\n", "------------MLNIAILCLSCYINYVSCVSVTSPPIIE gi|302171748|gb|ADK97775.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSLTSSAIIE gi|302171772|gb|ADK97787.1|/1-\n", "------------MLNIAILCLSCYINYVSCLSVTSPAIIE gi|225423294|gb|ACN91153.1|/1-\n", "...\n", "------------MLNIAILSLSCYINYVSCLSLTSPAIIE gi|302171764|gb|ADK97783.1|/1-\n" ] } ], "source": [ "print(align[:,0:40])" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# And now for a brief foray into marine microbiology...\n", "\n", "https://MSCBIO2025.github.io/files/hydra.fasta" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "slide" } }, "source": [ "# Project\n", "\n", "We have a gene (alr2), but what part of the gene is responsible for allorecognition?\n", "\n", "Given 179 sequences, how might we find out?\n", "\n", "https://MSCBIO2025.github.io/files/hydra179.aln" ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "fragment" } }, "source": [ "Find the part of the sequence that changes the most.\n", " * Count number of distinct residues at each position. Plot.\n", " * Count number of distinct *subsequences* of length 10 at each position. Plot." ] }, { "cell_type": "markdown", "metadata": { "slideshow": { "slide_type": "subslide" } }, "source": [ "" ] }, { "cell_type": "code", "execution_count": 73, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "--2023-09-21 12:27:54-- https://mscbio2025.github.io/files/hydra179.aln\n", "Resolving mscbio2025.github.io (mscbio2025.github.io)... 185.199.110.153, 185.199.109.153, 185.199.108.153, ...\n", "Connecting to mscbio2025.github.io (mscbio2025.github.io)|185.199.110.153|:443... connected.\n", "HTTP request sent, awaiting response... 301 Moved Permanently\n", "Location: https://mscbio2025.net/files/hydra179.aln [following]\n", "--2023-09-21 12:27:54-- https://mscbio2025.net/files/hydra179.aln\n", "Resolving mscbio2025.net (mscbio2025.net)... 185.199.111.153, 185.199.108.153, 185.199.109.153, ...\n", "Connecting to mscbio2025.net (mscbio2025.net)|185.199.111.153|:443... connected.\n", "HTTP request sent, awaiting response... 200 OK\n", "Length: 218936 (214K) [application/octet-stream]\n", "Saving to: ‘hydra179.aln’\n", "\n", "hydra179.aln 100%[===================>] 213.80K --.-KB/s in 0.01s \n", "\n", "2023-09-21 12:27:55 (16.1 MB/s) - ‘hydra179.aln’ saved [218936/218936]\n", "\n" ] } ], "source": [ "!wget https://MSCBIO2025.github.io/files/hydra179.aln" ] }, { "cell_type": "code", "execution_count": 74, "metadata": { "slideshow": { "slide_type": "notes" } }, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXAAAAD4CAYAAAD1jb0+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8/fFQqAAAACXBIWXMAAAsTAAALEwEAmpwYAAARJ0lEQVR4nO3db4xldX3H8fenLIhFy4JMNxvWdDASDA8qkAmFaEwL1SgY4QEhENNuGppNWttobGKXNmli0gfQB/5p0mg3ot0H/oGiFgLxD10xTZtmdVZAgZWy0CUuYdnRiv8e1KLfPrhnYZy9O3Pn/pv7m32/ksk953fOved779757LnfOefcVBWSpPb82kYXIEkajgEuSY0ywCWpUQa4JDXKAJekRm2Z5sbOO++8mp+fn+YmJal5Bw4c+H5Vza0cn2qAz8/Ps7i4OM1NSlLzkjzTb9wWiiQ1ygCXpEYZ4JLUKANckhplgEtSowxwSWrUQAGeZGuSu5N8N8nBJFcmOTfJA0me7G7PmXSxkqSXDboH/lHgy1X1BuCNwEFgN7Cvqi4E9nXzkqQpWTPAk5wNvAW4A6Cqfl5VLwDXAXu71fYC10+mRElSP4PsgV8ALAGfSvJQkk8kOQvYVlXPdescBbb1u3OSXUkWkywuLS2Np2pJ0kABvgW4DPhYVV0K/IwV7ZLqfa1P36/2qao9VbVQVQtzcyecyi9JGtIgAX4EOFJV+7v5u+kF+vNJtgN0t8cmU6IkqZ81A7yqjgLfS3JRN3Q18DhwL7CzG9sJ3DORCiVJfQ16NcI/Bz6d5AzgaeCP6IX/XUluAZ4BbpxMiZKkfgYK8Kp6GFjos+jqsVYjSRqYZ2JKUqMMcElqlAEuSY0ywCWpUQb4CvO77z+ltiupXQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNMsAlqVEGuCQ1ygCfIZ6NKWk9DHBJapQBLkmNMsAlqVEGuCQ1ygCXpEYZ4JLUKAN8nSZxqJ+HD0oahgEuSY0ywCWpUQa4JDVqyyArJTkM/AT4BfBiVS0kORe4E5gHDgM3VtUPJ1OmJGml9eyB/15VXVJVC938bmBfVV0I7OvmJUlTMkoL5Tpgbze9F7h+5GokSQMbNMAL+GqSA0l2dWPbquq5bvoosK3fHZPsSrKYZHFpaWnEcjeWh/tJmiUD9cCBN1fVs0l+E3ggyXeXL6yqSlL97lhVe4A9AAsLC33XkSSt30B74FX1bHd7DPgicDnwfJLtAN3tsUkVKUk60ZoBnuSsJK8+Pg28DXgUuBfY2a22E7hnUkVKkk40SAtlG/DFJMfX/0xVfTnJN4G7ktwCPAPcOLkyNy/76pKGtWaAV9XTwBv7jP8AuHoSRUmS1uaZmJLUKAN8xthSkTQoA1ySGmWAS1KjDHBJapQBvoz9Z0ktMcAlqVEGuCQ1ygDvw1aKpBYY4JLUKANckhplgHfW0zaZ332/bRZJG84Al6RGGeCS1CgDXJIaZYBLUqMMcElqlAEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGmWAS1KjDHBJapQBLkmNGjjAk5yW5KEk93XzFyTZn+RQkjuTnDG5MiVJK61nD/y9wMFl87cDH66q1wM/BG4ZZ2GSpNUNFOBJdgDXAp/o5gNcBdzdrbIXuH4C9UmSTmLQPfCPAB8AftnNvwZ4oape7OaPAOf3u2OSXUkWkywuLS2NUqskaZk1AzzJO4FjVXVgmA1U1Z6qWqiqhbm5uWEeQpLUx5YB1nkT8K4k1wBnAr8BfBTYmmRLtxe+A3h2cmVKklZacw+8qm6tqh1VNQ/cBHytqt4NPAjc0K22E7hnYlVKkk4wynHgfwm8P8khej3xO8ZTkiRpEIO0UF5SVV8Hvt5NPw1cPv6SJEmD8ExMSWqUAS5JjTLAJalRBrgkNcoAH8H87vubelxJm4sBLkmNMsAlqVEGuCQ1ygA/CfvQkmadAS5JjTLAJalRBvgqhm2j2H6RNA0GuCQ1ygCXpEYZ4AMatC1i+0TStBjgktQoA1ySGmWAS1KjDPABDNPXthcuadIMcElqlAEuSY0ywMfItomkaTLAJalRBrgkNcoAl6RGrRngSc5M8o0kjyR5LMkHu/ELkuxPcijJnUnOmHy5kqTjBtkD/1/gqqp6I3AJ8PYkVwC3Ax+uqtcDPwRumViVkqQTrBng1fPTbvb07qeAq4C7u/G9wPWTKFCS1N9APfAkpyV5GDgGPAA8BbxQVS92qxwBzj/JfXclWUyyuLS0NIaSp8tDAyXNqoECvKp+UVWXADuAy4E3DLqBqtpTVQtVtTA3NzdclZKkE6zrKJSqegF4ELgS2JpkS7doB/DseEuTJK1mkKNQ5pJs7aZfCbwVOEgvyG/oVtsJ3DOhGiVJfWxZexW2A3uTnEYv8O+qqvuSPA58LsnfAg8Bd0ywzomyzy2pRWsGeFV9G7i0z/jT9PrhkqQN4JmYktQoA1ySGmWAS1KjDHBJapQBLkmNMsBnlIc2SlqLAS5JjTLAJalRBrgkNcoAl6RGGeCS1CgDXJIaZYBP2GqHA3qooKRRGOCS1CgDXJIaZYBP0PEWia0SSZNggEtSowxwSWrUIN+JqTXYIpG0EdwDl6RGGeCS1CgDXJIaZYCPyP63pI1igEtSowxwSWqUAS5JjVozwJO8NsmDSR5P8liS93bj5yZ5IMmT3e05ky9XknTcIHvgLwJ/UVUXA1cA70lyMbAb2FdVFwL7unlJ0pSsGeBV9VxVfaub/glwEDgfuA7Y2622F7h+QjVKkvpYVw88yTxwKbAf2FZVz3WLjgLbTnKfXUkWkywuLS2NUqumxEMjpTYMHOBJXgV8HnhfVf14+bKqKqD63a+q9lTVQlUtzM3NjVSsJOllAwV4ktPphfenq+oL3fDzSbZ3y7cDxyZToiSpn0GOQglwB3Cwqj60bNG9wM5ueidwz/jL0zgM0xLZrG2U+d33v/Qzqcdfz7g0ikEuJ/sm4A+A7yR5uBv7K+A24K4ktwDPADdOpEJJUl9rBnhV/TuQkyy+erzlSJIG5ZmYktSoTR3gk+x1ngoGfe1aeI1bqFFar00d4JK0mRngktSoUyLAV/v4fKp8tD7+PFc7zG3WW06j1rbRz22jt6/N55QIcEnajAxwSWqUAS5JjdoUAb68t3iyXu9G93enue3N0Gtd+e+4/HY9/5arrTvu1+lktUqTsikCXJJORQa4JDVq0wR4vzaKftUwrYTVWlLr2eYw7YS1Dn0c9P4b4WTvxxYO11Q7Nk2AS9KpxgCXpEYZ4BtkEh+hh21RDLreakf2jPNs12EuomVLQqciA1ySGmWAS1KjDHBJatQg34m5aW2mvumwh9yNc/353fdz+LZrV1135TqjGmcffOX9x1nnRlvP6z7ufyNNjnvgktQoA1ySGtV0gK88o62F73Cc1AWUTrbsZG2MSWnxSxfWep1GPSN0PXWMup3lZ3oO83jTeq6zovXn2XSAS9KpzACXpEYZ4JLUqDUDPMknkxxL8uiysXOTPJDkye72nMmW+atG7RFulHH1Fwfpza61jWn1czdi2+sx6Os0iddzGn+fmNTlBjbTVRVbfg6D7IH/E/D2FWO7gX1VdSGwr5uXJE3RmgFeVf8G/M+K4euAvd30XuD68ZYlSVrLsD3wbVX1XDd9FNh2shWT7EqymGRxaWlpyM2dqMXD1darlTbEMFqtf1qHgfb7DtBJtyxG/dKO1rRa93Ij/xGzqgqoVZbvqaqFqlqYm5sbdXOSpM6wAf58ku0A3e2x8ZUkSRrEsAF+L7Czm94J3DOeciRJgxrkMMLPAv8JXJTkSJJbgNuAtyZ5Evj9bn7iNsthS+qvlX/bYXvFgz7eJLXyGk/aZnkd1rycbFXdfJJFV4+5FknSOngmpiQ1qpkA3ywfeUbVwutwKnwp8biv2reeszLXOz7sepO6/8kec1pXzpz2FTonqZkAlyT9KgNckhplgEtSowzwhrTap1tNy89plmof9cupR+2fj+sKof2+SWhar/O4Dw+dBgNckhplgEtSo9Y8kUcbb373/Ry+7dqNLmNiWvioOms26sufp/VYJ2tnjPJ7MMqXO8/q75974JLUKANckhplgEtSowzwRtgn1rgNcnXPYfrG41pvtfuu/Mai1R53PVcxHcfv2TSvmmqAS1KjDHBJalR6X2k5HQsLC7W4uDjUfW0hSBvj8G3Xruv37/ghd5vpd3bQwwj7PedxHIKY5EBVLawcdw9ckhplgEtSo2yhSNIAlrdCVubRWm2mUdsotlAkaZMxwCWpUQa4JDXKqxFK0gBW63Fv1N/o3AOXpEYZ4JLUKANckho1UoAneXuSJ5IcSrJ7XEVJktY2dIAnOQ34B+AdwMXAzUkuHldhkqTVjbIHfjlwqKqerqqfA58DrhtPWZKktYxyGOH5wPeWzR8BfmflSkl2Abu62Z8meWLI7Z0HfH/I+06btU5OS/Va62Q0V2tuH/lxfqvf4MSPA6+qPcCeUR8nyWK/awHMImudnJbqtdbJsNaXjdJCeRZ47bL5Hd2YJGkKRgnwbwIXJrkgyRnATcC94ylLkrSWoVsoVfVikj8DvgKcBnyyqh4bW2UnGrkNM0XWOjkt1Wutk2GtnaleD1ySND6eiSlJjTLAJalRTQT4rJ2yn+STSY4leXTZ2LlJHkjyZHd7TjeeJH/f1f7tJJdNudbXJnkwyeNJHkvy3lmtN8mZSb6R5JGu1g924xck2d/VdGf3R3OSvKKbP9Qtn59WrctqPi3JQ0num+VakxxO8p0kDydZ7MZm7j3QbX9rkruTfDfJwSRXzmKtSS7qXs/jPz9O8r6p1lpVM/1D7w+kTwGvA84AHgEu3uCa3gJcBjy6bOzvgN3d9G7g9m76GuBLQIArgP1TrnU7cFk3/Wrgv+hd+mDm6u22+apu+nRgf1fDXcBN3fjHgT/ppv8U+Hg3fRNw5wa8F94PfAa4r5ufyVqBw8B5K8Zm7j3QbX8v8Mfd9BnA1lmtdVnNpwFH6Z1wM7Vap/5Eh3hhrgS+smz+VuDWGahrfkWAPwFs76a3A0900/8I3NxvvQ2q+x7grbNeL/DrwLfond37fWDLyvcDvSOgruymt3TrZYo17gD2AVcB93W/mLNaa78An7n3AHA28N8rX5tZrHVFfW8D/mPatbbQQul3yv75G1TLarZV1XPd9FFgWzc9M/V3H9svpbdnO5P1di2Jh4FjwAP0Pn29UFUv9qnnpVq75T8CXjOtWoGPAB8AftnNv4bZrbWAryY5kN7lLWA23wMXAEvAp7rW1CeSnDWjtS53E/DZbnpqtbYQ4M2p3n+vM3V8ZpJXAZ8H3ldVP16+bJbqrapfVNUl9PZuLwfesLEV9ZfkncCxqjqw0bUM6M1VdRm9q4e+J8lbli+coffAFnrtyY9V1aXAz+i1IV4yQ7UC0P2d413AP69cNulaWwjwVk7Zfz7JdoDu9lg3vuH1JzmdXnh/uqq+0A3PbL0AVfUC8CC9NsTWJMdPOltez0u1dsvPBn4wpRLfBLwryWF6V+K8CvjojNZKVT3b3R4DvkjvP8dZfA8cAY5U1f5u/m56gT6LtR73DuBbVfV8Nz+1WlsI8FZO2b8X2NlN76TXaz4+/ofdX6CvAH607OPVxCUJcAdwsKo+NMv1JplLsrWbfiW9Xv1BekF+w0lqPf4cbgC+1u3xTFxV3VpVO6pqnt578mtV9e5ZrDXJWUlefXyaXr/2UWbwPVBVR4HvJbmoG7oaeHwWa13mZl5unxyvaTq1TrvZP+QfCK6hd/TEU8Bfz0A9nwWeA/6P3h7DLfT6mfuAJ4F/Bc7t1g29L754CvgOsDDlWt9M7yPct4GHu59rZrFe4LeBh7paHwX+pht/HfAN4BC9j6mv6MbP7OYPdctft0Hvh9/l5aNQZq7WrqZHup/Hjv8OzeJ7oNv+JcBi9z74F+CcGa71LHqfpM5eNja1Wj2VXpIa1UILRZLUhwEuSY0ywCWpUQa4JDXKAJekRhngktQoA1ySGvX/5ps9HOI3NtwAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "import matplotlib.pyplot as plt\n", "from Bio import AlignIO\n", "\n", "def countatpos(seqs,pos):\n", " '''count the number of unique subsequences of length 10 at pos in seqs'''\n", " vals = set()\n", " for s in seqs:\n", " vals.add(str(s[pos:pos+10].seq))\n", " return len(vals)\n", "\n", "seqs = AlignIO.read(\"hydra179.aln\",'clustal');\n", "\n", "xaxis = list()\n", "yaxis = list()\n", "for i in range(seqs.get_alignment_length()):\n", " xaxis.append(i)\n", " yaxis.append(countatpos(seqs,i))\n", "\n", "plt.bar(xaxis,yaxis,width=1);" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "celltoolbar": "Slideshow", "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.6" } }, "nbformat": 4, "nbformat_minor": 4 }