docs/notebooks/Dseq_Features.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Working with Features using the Dseqrecord class\n",
    "\n",
    "> Before working with features, check how to import sequences from files in the [Importing_Seqs notebook](./Importing_Seqs.ipynb).\n",
    ">\n",
    "> For full library documentation, visit [here](https://pydna-group.github.io/pydna/).\n",
    "\n",
    "Some sequence file formats (like Genbank) include features, describing key biological properties of sequence regions. In Genbank, features \"include genes, gene products, as well as regions of biological significance reported in the sequence.\" (See [here](https://www.ncbi.nlm.nih.gov/genbank/samplerecord/) for a description of a Genbank file and associated terminologies/annotations) Examples include coding sequences (CDS), introns, promoters, etc.\n",
    "\n",
    "pydna offers many ways to easily view, add, extract, and write features into a Genbank file via the `Dseqrecord` class. After reading a file into a `Dseqrecord` object, we can check out the list of features in the record using the following code. This example uses the sample record [U49845](https://www.ncbi.nlm.nih.gov/genbank/samplerecord/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/pydna-group/pydna/blob/master/docs/notebooks/Dseq_Features.ipynb\">\n",
    "  <img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/>\n",
    "</a>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%capture\n",
    "# Install pydna (only when running on Colab)\n",
    "import sys\n",
    "if 'google.colab' in sys.modules:\n",
    "    %pip install pydna"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: source\n",
      "location: [0:5028](+)\n",
      "qualifiers:\n",
      "    Key: chromosome, Value: ['IX']\n",
      "    Key: db_xref, Value: ['taxon:4932']\n",
      "    Key: mol_type, Value: ['genomic DNA']\n",
      "    Key: organism, Value: ['Saccharomyces cerevisiae']\n",
      "\n",
      "type: mRNA\n",
      "location: [<0:>206](+)\n",
      "qualifiers:\n",
      "    Key: product, Value: ['TCP1-beta']\n",
      "\n",
      "type: CDS\n",
      "location: [<0:206](+)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['3']\n",
      "    Key: product, Value: ['TCP1-beta']\n",
      "    Key: protein_id, Value: ['AAA98665.1']\n",
      "    Key: translation, Value: ['SSIYNGISTSGLDLNNGTIADMRQLGIVESYKLKRAVVSSASEAAEVLLRVDNIIRARPRTANRQHM']\n",
      "\n",
      "type: gene\n",
      "location: [<686:>3158](+)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['AXL2']\n",
      "\n",
      "type: mRNA\n",
      "location: [<686:>3158](+)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['AXL2']\n",
      "    Key: product, Value: ['Axl2p']\n",
      "\n",
      "type: CDS\n",
      "location: [686:3158](+)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['1']\n",
      "    Key: gene, Value: ['AXL2']\n",
      "    Key: note, Value: ['plasma membrane glycoprotein']\n",
      "    Key: product, Value: ['Axl2p']\n",
      "    Key: protein_id, Value: ['AAA98666.1']\n",
      "    Key: translation, Value: ['MTQLQISLLLTATISLLHLVVATPYEAYPIGKQYPPVARVNESFTFQISNDTYKSSVDKTAQITYNCFDLPSWLSFDSSSRTFSGEPSSDLLSDANTTLYFNVILEGTDSADSTSLNNTYQFVVTNRPSISLSSDFNLLALLKNYGYTNGKNALKLDPNEVFNVTFDRSMFTNEESIVSYYGRSQLYNAPLPNWLFFDSGELKFTGTAPVINSAIAPETSYSFVIIATDIEGFSAVEVEFELVIGAHQLTTSIQNSLIINVTDTGNVSYDLPLNYVYLDDDPISSDKLGSINLLDAPDWVALDNATISGSVPDELLGKNSNPANFSVSIYDTYGDVIYFNFEVVSTTDLFAISSLPNINATRGEWFSYYFLPSQFTDYVNTNVSLEFTNSSQDHDWVKFQSSNLTLAGEVPKNFDKLSLGLKANQGSQSQELYFNIIGMDSKITHSNHSANATSTRSSHHSTSTSSYTSSTYTAKISSTSAAATSSAPAALPAANKTSSHNKKAVAIACGVAIPLGVILVALICFLIFWRRRRENPDDENLPHAISGPDLNNPANKPNQENATPLNNPFDDDASSYDDTSIARRLAALNTLKLDNHSATESDISSVDEKRDSLSGMNTYNDQFQSQSKEELLAKPPVQPPESPFFDPQNRSSSVYMDSEPAVNKSWRYTGNLSPVSDIVRDSYGSQKTVDTEKLFDLEAPEKEKRTSRDVTMSSLDPWNSNISPSPVRKSVTPSPYNVTKHRNRHLQNIQDSQSGKNGITPTTMSTSSSDDFVPVKDGENFCWVHSMEPDRRPSKKRLVDFSNKSNVNVGQVKDIHGRIPEML']\n",
      "\n",
      "type: gene\n",
      "location: [<3299:>4037](-)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['REV7']\n",
      "\n",
      "type: mRNA\n",
      "location: [<3299:>4037](-)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['REV7']\n",
      "    Key: product, Value: ['Rev7p']\n",
      "\n",
      "type: CDS\n",
      "location: [3299:4037](-)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['1']\n",
      "    Key: gene, Value: ['REV7']\n",
      "    Key: product, Value: ['Rev7p']\n",
      "    Key: protein_id, Value: ['AAA98667.1']\n",
      "    Key: translation, Value: ['MNRWVEKWLRVYLKCYINLILFYRNVYPPQSFDYTTYQSFNLPQFVPINRHPALIDYIEELILDVLSKLTHVYRFSICIINKKNDLCIEKYVLDFSELQHVDKDDQIITETEVFDEFRSSLNSLIMHLEKLPKVNDDTITFEAVINAIELELGHKLDRNRRVDSLEEKAEIERDSNWVKCQEDENLPDNNGFQPPKIKLTSLVGSDVGPLIIHQFSEKLISGDDKILNGVYSQYEEGESIFGSLF']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "from pydna.dseqrecord import Dseqrecord\n",
    "from pydna.parsers import parse\n",
    "\n",
    "#Import your file into python. \n",
    "file_path = \"./U49845.gb\"\n",
    "records = parse(file_path)\n",
    "sample_record = records[0]\n",
    "\n",
    "# List all features\n",
    "for feature in sample_record.features:\n",
    "    print(feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Additional ways to view and search for particular features are shown at the bottom of the page under \"Other Methods to Viewing Features\""
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding Features and Qualifiers\n",
    "\n",
    "To add new feature to describe a region of interest to a record, for instance a region that you would like to perform a PCR, you need to create a `SeqFeature` (sequence feature). The minimal information required is:\n",
    "* A `FeatureLocation`: position of the feature in the sequence.\n",
    "* The `type` of feature you want to add.\n",
    "\n",
    "\n",
    "🚨🚨 **VERY IMPORTANT** 🚨🚨. Note that `FeatureLocation`s are like python ranges (zero-based open intervals), whereas in GenBank files, locations are one-based closed intervals. For instance, the following code adds a new feature from the 2nd to the 5th nucleotide (`FeatureLocation(3, 15)`), of the `gene` type, but in the GenBank file will be represented as `4..15`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: gene\n",
      "location: [3:15]\n",
      "qualifiers:\n",
      "\n",
      "LOCUS       name                      19 bp    DNA     linear   UNK 01-JAN-1980\n",
      "DEFINITION  description.\n",
      "ACCESSION   id\n",
      "VERSION     id\n",
      "KEYWORDS    .\n",
      "SOURCE      .\n",
      "  ORGANISM  .\n",
      "            .\n",
      "FEATURES             Location/Qualifiers\n",
      "     gene            4..15\n",
      "ORIGIN\n",
      "        1 aaaatgcgta cgtgaacgt\n",
      "//\n"
     ]
    }
   ],
   "source": [
    "from Bio.SeqFeature import FeatureLocation, SeqFeature\n",
    "\n",
    "# Create a dummy record\n",
    "dummy_record = Dseqrecord(\"aaaATGCGTACGTGAacgt\")\n",
    "\n",
    "# Define the locations of a CDS\n",
    "location = FeatureLocation(3, 15)\n",
    "\n",
    "# Create a SeqFeature with the type mRNA\n",
    "my_feature = SeqFeature(location=location, type=\"gene\")\n",
    "\n",
    "# Add my_feature to dummy_record with .append\n",
    "dummy_record.features.append(my_feature)\n",
    "\n",
    "# Confirm that my_feature has been added\n",
    "print(dummy_record.features[-1])\n",
    "\n",
    "# Print the feature in GenBank format (see how the location is `4..15`)\n",
    "print(dummy_record.format(\"genbank\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To give further information about a feature, we can add a qualifier using the `qualifiers` property of `SeqFeature`, which contains a dictionary of qualifiers. For instance, if I would like to note a new feature of type 'domain', between 3-9 bases as my region of interest, I can instantiate the `SeqFeature` class object as such.\n",
    "\n",
    "> Note that a new feature is always added to the last position of the features list."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">> Feature was added:\n",
      "type: domain\n",
      "location: [3:9]\n",
      "qualifiers:\n",
      "    Key: Note, Value: ['Region of interest']\n",
      "\n",
      "\n",
      ">> GenBank format:\n",
      "LOCUS       name                      19 bp    DNA     linear   UNK 01-JAN-1980\n",
      "DEFINITION  description.\n",
      "ACCESSION   id\n",
      "VERSION     id\n",
      "KEYWORDS    .\n",
      "SOURCE      .\n",
      "  ORGANISM  .\n",
      "            .\n",
      "FEATURES             Location/Qualifiers\n",
      "     gene            4..15\n",
      "     domain          4..9\n",
      "                     /Note=\"Region of interest\"\n",
      "ORIGIN\n",
      "        1 aaaatgcgta cgtgaacgt\n",
      "//\n"
     ]
    }
   ],
   "source": [
    "location = FeatureLocation(3, 9)\n",
    "\n",
    "# Create a SeqFeature with a qualifier\n",
    "my_feature2 = SeqFeature(location=location, type=\"domain\", qualifiers={\"Note\": [\"Region of interest\"]})\n",
    "\n",
    "# Add my_feature to my_record with .append\n",
    "dummy_record.features.append(my_feature2)\n",
    "\n",
    "# Confirm that my_feature has been added\n",
    "print('>> Feature was added:')\n",
    "print(dummy_record.features[-1])\n",
    "print()\n",
    "\n",
    "# Print the feature in GenBank format\n",
    "print('>> GenBank format:')\n",
    "print(dummy_record.format(\"genbank\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**🤔 Best practices for qualifiers:**\n",
    "\n",
    "The values in the `qualifiers` dictionary should be lists. The reason for this is that in a GenBank file, a single feature can have multiple values for a single qualifier. Below is a real world of the ase1 CDS example from the _S. pombe_ genome in EMBL format:\n",
    "\n",
    "```\n",
    "FT   CDS             join(1878362..1878785,1878833..1880604)\n",
    "FT                   /colour=2\n",
    "FT                   /primary_name=\"ase1\"\n",
    "FT                   /product=\"antiparallel microtubule cross-linking factor\n",
    "FT                   Ase1\"\n",
    "FT                   /systematic_id=\"SPAPB1A10.09\"\n",
    "FT                   /controlled_curation=\"term=species distribution, conserved\n",
    "FT                   in eukaryotes; date=20081110\"\n",
    "FT                   /controlled_curation=\"term=species distribution, conserved\n",
    "FT                   in metazoa; date=20081110\"\n",
    "FT                   /controlled_curation=\"term=species distribution, conserved\n",
    "FT                   in vertebrates; date=20081110\"\n",
    "FT                   /controlled_curation=\"term=species distribution,\n",
    "FT                   predominantly single copy (one to one); date=20081110\"\n",
    "FT                   /controlled_curation=\"term=species distribution, conserved\n",
    "FT                   in fungi; date=20081110\"\n",
    "FT                   /controlled_curation=\"term=species distribution, conserved\n",
    "FT                   in eukaryotes only; date=20081110\"\n",
    "```\n",
    "\n",
    "Note how there are several `controlled_curation` qualifiers, therefore it makes sense to store them as a list.\n",
    "\n",
    "By default, you can add any type of object in the qualifiers dictionary, and most things will work if you add a string. However, you risk overwriting the existing value for a qualifier, so best practice is:\n",
    "1. Check if the qualifier already exists using `if \"qualifier_name\" in feature.qualifiers`\n",
    "2. If it exists, append to the existing list of values using `feature.qualifiers[\"qualifier_name\"].append(\"new_value\")`\n",
    "3. If it does not exist, add it to the qualifiers dictionary using `feature.qualifiers[\"qualifier_name\"] = [\"new_value\"]`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that `Bio.SeqFeatures` does not automatically assume a sequence strand for the feature. If you would like to refer to a feature on the positive or minus strand, you can add a parameter in `FeatureLocation` specifying `strand=+1` or `strand=-1`. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: domain\n",
      "location: [15:19](-)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['example_domain']\n",
      "\n",
      "LOCUS       name                      19 bp    DNA     linear   UNK 01-JAN-1980\n",
      "DEFINITION  description.\n",
      "ACCESSION   id\n",
      "VERSION     id\n",
      "KEYWORDS    .\n",
      "SOURCE      .\n",
      "  ORGANISM  .\n",
      "            .\n",
      "FEATURES             Location/Qualifiers\n",
      "     gene            4..15\n",
      "     domain          4..9\n",
      "                     /Note=\"Region of interest\"\n",
      "     domain          complement(16..19)\n",
      "                     /gene=\"example_domain\"\n",
      "ORIGIN\n",
      "        1 aaaatgcgta cgtgaacgt\n",
      "//\n"
     ]
    }
   ],
   "source": [
    "#Create a location specifying the minus strand\n",
    "location = FeatureLocation(15, 19, strand=-1)\n",
    "\n",
    "my_feature3 = SeqFeature(location=location, type=\"domain\", qualifiers={\"gene\":[\"example_domain\"]})\n",
    "\n",
    "dummy_record.features.append(my_feature3)\n",
    "\n",
    "print(dummy_record.features[-1])\n",
    "\n",
    "print(dummy_record.format(\"genbank\"))\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding a Feature with Parts\n",
    "\n",
    "To add a feature with parts, like a CDS with introns, we need to use a `CompoundLocation` object when creating a `SeqFeature`.\n",
    "\n",
    "The example code below adds a CDS with two parts, between 3-9bp and 12-15bp, to my features list. In a real-world scenario this would represent a CDS with an intron that skips the `ACG` codon: ATGCGT~~ACG~~TGA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: CDS\n",
      "location: join{[3:9], [12:15]}\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['example_gene']\n",
      "\n",
      "LOCUS       name                      19 bp    DNA     linear   UNK 01-JAN-1980\n",
      "DEFINITION  description.\n",
      "ACCESSION   id\n",
      "VERSION     id\n",
      "KEYWORDS    .\n",
      "SOURCE      .\n",
      "  ORGANISM  .\n",
      "            .\n",
      "FEATURES             Location/Qualifiers\n",
      "     gene            4..15\n",
      "     domain          4..9\n",
      "                     /Note=\"Region of interest\"\n",
      "     domain          complement(16..19)\n",
      "                     /gene=\"example_domain\"\n",
      "     CDS             join(4..9,13..15)\n",
      "                     /gene=\"example_gene\"\n",
      "ORIGIN\n",
      "        1 aaaatgcgta cgtgaacgt\n",
      "//\n"
     ]
    }
   ],
   "source": [
    "from Bio.SeqFeature import CompoundLocation\n",
    "\n",
    "# Define the locations of the CDS\n",
    "locations = [FeatureLocation(3, 9), FeatureLocation(12, 15)]\n",
    "\n",
    "# Create a compound location from these parts\n",
    "compound_location = CompoundLocation(locations)\n",
    "\n",
    "# Create a SeqFeature with this compound location, including type and qualifiers. \n",
    "cds_feature = SeqFeature(location=compound_location, type=\"CDS\", qualifiers={\"gene\": [\"example_gene\"]})\n",
    "\n",
    "# Add the feature to the Dseqrecord\n",
    "dummy_record.features.append(cds_feature)\n",
    "\n",
    "print(dummy_record.features[-1])\n",
    "\n",
    "print(dummy_record.format(\"genbank\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can even extract a protein record as follows (see how the protein sequence is `MR`, skipping the intron):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ID: id\n",
      "Name: name\n",
      "Description: description\n",
      "Number of features: 0\n",
      "/molecule_type=DNA\n",
      "ProteinSeq('MR')\n"
     ]
    }
   ],
   "source": [
    "sub_record = dummy_record.features[-1].extract(dummy_record)\n",
    "\n",
    "print(sub_record.translate())\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Standard Feature Types and Qualifiers\n",
    "\n",
    "`pydna` and `Bio.SeqFeature` suppports all the conventional feature types through the `type` parameters. A non-exhaustive list include gene, CDS, promoter, exon, intron, 5' UTR, 3' UTR, terminator, enhancer, and RBS. You can also define custom features, which could be useful for synthetic biology applications. For instance, you might want to have Bio_brick or spacer features to describe a synthetic standardised plasmid construct.  \n",
    "\n",
    "It is important to note that while `pydna` and `Bio.SeqFeature` does not restrict the feature types you can use, sticking to standard types helps maintain compatibility with other bioinformatics tools and databases. Please refer to the official [GenBank_Feature_Table](https://www.insdc.org/submitting-standards/feature-table/#2), that lists the standard feature types and their associated qualifiers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Further documentation for `SeqFeature`, `CompoundLocation`, and `FeatureLocation` can be found in the `SeqFeature` module [here](https://biopython.org/docs/1.75/api/Bio.SeqFeature.html). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Handling Origin Spanning Features\n",
    "\n",
    "An origin spanning feature is a special type of feature that crosses over a circular sequence's origin. In pydna, such a feature is represented as a feature with parts, joining the part of the sequence immediately before the origin and immediately after the origin. They can be added using `CompoundLocation` as normal. \n",
    "\n",
    "An origin spanning feature, between base 19 to base 6, in a 25bp long circular sequence, is represented like so:   \n",
    "\n",
    "```\n",
    "type: gene \n",
    "location: join{[19:25](+), [0:6](+)}  \n",
    "qualifiers: gene, Value: example_gene  \n",
    "```\n",
    "\n",
    "This feature will be displayed as a single feature in SnapGene viewer and Benchling, since they support this convention."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">> Feature:\n",
      "type: misc\n",
      "location: join{[19:25], [0:6]}\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['example origin spanning gene']\n",
      "\n",
      ">> Feature sequence:\n",
      "ATGCGTACGTGA\n",
      "\n",
      ">> GenBank format:\n",
      "LOCUS       name                      25 bp    DNA     circular UNK 01-JAN-1980\n",
      "DEFINITION  description.\n",
      "ACCESSION   id\n",
      "VERSION     id\n",
      "KEYWORDS    .\n",
      "SOURCE      .\n",
      "  ORGANISM  .\n",
      "            .\n",
      "FEATURES             Location/Qualifiers\n",
      "     misc            join(20..25,1..6)\n",
      "                     /gene=\"example origin spanning gene\"\n",
      "ORIGIN\n",
      "        1 acgtgaaaaa aaaaaaaaaa tgcgt\n",
      "//\n"
     ]
    }
   ],
   "source": [
    "circular_record = Dseqrecord('ACGTGAaaaaaaaaaaaaaATGCGT', circular=True)\n",
    "\n",
    "location = [FeatureLocation(19,25), FeatureLocation(0, 6)]\n",
    "ori_feat_location = CompoundLocation(location)\n",
    "ori_feature = SeqFeature(location=ori_feat_location, type=\"misc\", qualifiers={\"gene\": [\"example origin spanning gene\"]})\n",
    "circular_record.features.append(ori_feature)\n",
    "\n",
    "print('>> Feature:')\n",
    "print(circular_record.features[-1])\n",
    "\n",
    "# Note how the feature sequence is extracted properly across the origin.\n",
    "print('>> Feature sequence:')\n",
    "print(circular_record.features[-1].extract(circular_record).seq)\n",
    "print()\n",
    "\n",
    "print('>> GenBank format:')\n",
    "print(circular_record.format(\"genbank\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Methods to Viewing Features\n",
    "\n",
    "pydna also provides the `list_features` method as a simple way to list all the features in a `Dseqrecord` object. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "+-----+------------------+-----+-------+-------+------+--------+------+\n",
      "| Ft# | Label or Note    | Dir | Sta   | End   |  Len | type   | orf? |\n",
      "+-----+------------------+-----+-------+-------+------+--------+------+\n",
      "|   0 | nd               | --> | 0     | 5028  | 5028 | source |  no  |\n",
      "|   1 | nd               | --> | <0    | >206  |  206 | mRNA   |  no  |\n",
      "|   2 | nd               | --> | <0    | 206   |  206 | CDS    |  no  |\n",
      "|   3 | nd               | --> | <686  | >3158 | 2472 | gene   | yes  |\n",
      "|   4 | nd               | --> | <686  | >3158 | 2472 | mRNA   | yes  |\n",
      "|   5 | N:plasma membran | --> | 686   | 3158  | 2472 | CDS    | yes  |\n",
      "|   6 | nd               | <-- | <3299 | >4037 |  738 | gene   | yes  |\n",
      "|   7 | nd               | <-- | <3299 | >4037 |  738 | mRNA   | yes  |\n",
      "|   8 | nd               | <-- | 3299  | 4037  |  738 | CDS    | yes  |\n",
      "+-----+------------------+-----+-------+-------+------+--------+------+\n"
     ]
    }
   ],
   "source": [
    "print(sample_record.list_features())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This method is convenient for checking-out a brief overview of each feature, without reading through an entire sequence record.\n",
    "\n",
    "Alternatively, we can look for specific features using their qualifiers. For instance:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Getting all CDS features:\n",
      "type: CDS\n",
      "location: [<0:206](+)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['3']\n",
      "    Key: product, Value: ['TCP1-beta']\n",
      "    Key: protein_id, Value: ['AAA98665.1']\n",
      "    Key: translation, Value: ['SSIYNGISTSGLDLNNGTIADMRQLGIVESYKLKRAVVSSASEAAEVLLRVDNIIRARPRTANRQHM']\n",
      "\n",
      "type: CDS\n",
      "location: [686:3158](+)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['1']\n",
      "    Key: gene, Value: ['AXL2']\n",
      "    Key: note, Value: ['plasma membrane glycoprotein']\n",
      "    Key: product, Value: ['Axl2p']\n",
      "    Key: protein_id, Value: ['AAA98666.1']\n",
      "    Key: translation, Value: ['MTQLQISLLLTATISLLHLVVATPYEAYPIGKQYPPVARVNESFTFQISNDTYKSSVDKTAQITYNCFDLPSWLSFDSSSRTFSGEPSSDLLSDANTTLYFNVILEGTDSADSTSLNNTYQFVVTNRPSISLSSDFNLLALLKNYGYTNGKNALKLDPNEVFNVTFDRSMFTNEESIVSYYGRSQLYNAPLPNWLFFDSGELKFTGTAPVINSAIAPETSYSFVIIATDIEGFSAVEVEFELVIGAHQLTTSIQNSLIINVTDTGNVSYDLPLNYVYLDDDPISSDKLGSINLLDAPDWVALDNATISGSVPDELLGKNSNPANFSVSIYDTYGDVIYFNFEVVSTTDLFAISSLPNINATRGEWFSYYFLPSQFTDYVNTNVSLEFTNSSQDHDWVKFQSSNLTLAGEVPKNFDKLSLGLKANQGSQSQELYFNIIGMDSKITHSNHSANATSTRSSHHSTSTSSYTSSTYTAKISSTSAAATSSAPAALPAANKTSSHNKKAVAIACGVAIPLGVILVALICFLIFWRRRRENPDDENLPHAISGPDLNNPANKPNQENATPLNNPFDDDASSYDDTSIARRLAALNTLKLDNHSATESDISSVDEKRDSLSGMNTYNDQFQSQSKEELLAKPPVQPPESPFFDPQNRSSSVYMDSEPAVNKSWRYTGNLSPVSDIVRDSYGSQKTVDTEKLFDLEAPEKEKRTSRDVTMSSLDPWNSNISPSPVRKSVTPSPYNVTKHRNRHLQNIQDSQSGKNGITPTTMSTSSSDDFVPVKDGENFCWVHSMEPDRRPSKKRLVDFSNKSNVNVGQVKDIHGRIPEML']\n",
      "\n",
      "type: CDS\n",
      "location: [3299:4037](-)\n",
      "qualifiers:\n",
      "    Key: codon_start, Value: ['1']\n",
      "    Key: gene, Value: ['REV7']\n",
      "    Key: product, Value: ['Rev7p']\n",
      "    Key: protein_id, Value: ['AAA98667.1']\n",
      "    Key: translation, Value: ['MNRWVEKWLRVYLKCYINLILFYRNVYPPQSFDYTTYQSFNLPQFVPINRHPALIDYIEELILDVLSKLTHVYRFSICIINKKNDLCIEKYVLDFSELQHVDKDDQIITETEVFDEFRSSLNSLIMHLEKLPKVNDDTITFEAVINAIELELGHKLDRNRRVDSLEEKAEIERDSNWVKCQEDENLPDNNGFQPPKIKLTSLVGSDVGPLIIHQFSEKLISGDDKILNGVYSQYEEGESIFGSLF']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Filter based on feature type\n",
    "print('Getting all CDS features:')\n",
    "cds_features = [f for f in sample_record.features if f.type == \"CDS\"]\n",
    "for feature in cds_features:\n",
    "    print(feature)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: gene\n",
      "location: [<3299:>4037](-)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['REV7']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "# Find a particular feature by its qualifier (e.g. gene name)\n",
    "rev7_cds_feature = next(f for f in sample_record.features if \n",
    "     f.type == \"gene\" and\n",
    "     \"gene\" in f.qualifiers and \"REV7\" in f.qualifiers[\"gene\"]\n",
    "     )\n",
    "\n",
    "print(rev7_cds_feature)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you would like to search for another type of features, simply replace the `\"gene\"` with your desired feature type in quotation marks."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Removing Features\n",
    "\n",
    "In pydna, we can search for the feature that we would like to remove using the feature's types or qualififers. For instance, we can modify the features list to exclude all CDS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: source\n",
      "location: [0:5028](+)\n",
      "qualifiers:\n",
      "    Key: chromosome, Value: ['IX']\n",
      "    Key: db_xref, Value: ['taxon:4932']\n",
      "    Key: mol_type, Value: ['genomic DNA']\n",
      "    Key: organism, Value: ['Saccharomyces cerevisiae']\n",
      "\n",
      "type: mRNA\n",
      "location: [<0:>206](+)\n",
      "qualifiers:\n",
      "    Key: product, Value: ['TCP1-beta']\n",
      "\n",
      "type: gene\n",
      "location: [<686:>3158](+)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['AXL2']\n",
      "\n",
      "type: mRNA\n",
      "location: [<686:>3158](+)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['AXL2']\n",
      "    Key: product, Value: ['Axl2p']\n",
      "\n",
      "type: gene\n",
      "location: [<3299:>4037](-)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['REV7']\n",
      "\n",
      "type: mRNA\n",
      "location: [<3299:>4037](-)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['REV7']\n",
      "    Key: product, Value: ['Rev7p']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Remove all CDS type features from my feature list\n",
    "sample_record.features = [f for f in sample_record.features if not (f.type == \"CDS\")]\n",
    "\n",
    "for feature in sample_record.features:\n",
    "    print(feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can also modify the features list to exclude a specific gene:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "type: source\n",
      "location: [0:5028](+)\n",
      "qualifiers:\n",
      "    Key: chromosome, Value: ['IX']\n",
      "    Key: db_xref, Value: ['taxon:4932']\n",
      "    Key: mol_type, Value: ['genomic DNA']\n",
      "    Key: organism, Value: ['Saccharomyces cerevisiae']\n",
      "\n",
      "type: mRNA\n",
      "location: [<0:>206](+)\n",
      "qualifiers:\n",
      "    Key: product, Value: ['TCP1-beta']\n",
      "\n",
      "type: gene\n",
      "location: [<686:>3158](+)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['AXL2']\n",
      "\n",
      "type: mRNA\n",
      "location: [<686:>3158](+)\n",
      "qualifiers:\n",
      "    Key: gene, Value: ['AXL2']\n",
      "    Key: product, Value: ['Axl2p']\n",
      "\n"
     ]
    }
   ],
   "source": [
    "#Exclude REV7 from my feature list\n",
    "sample_record.features = [f for f in sample_record.features if not ('gene' in f.qualifiers and 'REV7' in f.qualifiers['gene'])]\n",
    "\n",
    "for feature in sample_record.features:\n",
    "    print(feature)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": ".venv",
   "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.12.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}