docs/notebooks/primer_design.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Primer design in pydna\n",
    "\n",
    "You can use `pydna` for primer design in different contexts, let's start with some basic primer functionalities."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/pydna-group/pydna/blob/master/docs/notebooks/primer_design.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\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Checking the Tm of a primer\n",
    "\n",
    "Primer design in pydna is very flexible, and supports different methods. For typical use-cases, we recommend using `tm_default`, which uses the method `Bio.SeqUtils.MeltingTemp` from biopython, nearest neighbor thermodynamics values from [SantaLucia & Hicks (2004)](https://pubmed.ncbi.nlm.nih.gov/15139820/) and common values for nucleotide concentration, salt concentration, etc. You can of course change those settings. For a full dive check `src/pydna/tm.py`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "55.047602641480864\n",
      "54.55481807340169\n"
     ]
    }
   ],
   "source": [
    "from pydna.tm import tm_default\n",
    "\n",
    "# The primers from the readme example\n",
    "print(tm_default(\"ATGCAAACAGTAATGATGGA\"))\n",
    "print(tm_default(\"ATTATCTTTTTCAGCAATAGAATCA\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using NEB Tm calculator\n",
    "\n",
    "If you are used to the NEB Tm calculator, and you want to use it programmatically, you can do so as well. The function takes three arguments:\n",
    "\n",
    "- `primer`: The primer sequence.\n",
    "- `conc`: The primer concentration.\n",
    "- `prodcode`: The product code, which you can find on [NEB's website](https://tmapi.neb.com/docs/productcodes).\n",
    "\n",
    "> **NOTE:** When you call the function, it will make a request to the NEB server. This makes it much slower than using the builtin methods. In addition, we cannot guarantee that the NEB server will always be available, nor that the calculations that they use will not change in the future, since the code is not available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "59\n",
      "57\n"
     ]
    }
   ],
   "source": [
    "from pydna.tm import tm_neb\n",
    "\n",
    "print(tm_neb(\"ATGCAAACAGTAATGATGGA\", 0.5, \"q5-0\"))\n",
    "print(tm_neb(\"ATTATCTTTTTCAGCAATAGAATCA\", 0.5, \"q5-0\"))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Designing primers for PCR\n",
    "\n",
    "Let's use pydna to amplify a region from a DNA sequence. You can use the `primer_design` function to design primers for a given target Tm (`target_tm`) and indicate a minimum primer hybridization length in basepairs (`limit`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Forward primer Tm: 59.71997924024873\n",
      "Forward primer sequence: ATGCAAACAGTAATGATGGATGAC\n",
      "\n",
      "Reverse primer Tm: 60.22377911083646\n",
      "Reverse primer sequence: TTATTCAGCAATAGAATCAGTGCTTTG\n"
     ]
    }
   ],
   "source": [
    "from pydna.dseqrecord import Dseqrecord\n",
    "from Bio.SeqFeature import SeqFeature, SimpleLocation\n",
    "from pydna.design import primer_design\n",
    "\n",
    "dna = Dseqrecord(\"ggttcaATGCAAACAGTAATGATGGATGACATTCAAAGCACTGATTCTATTGCTGAATAAcatttacatca\")\n",
    "\n",
    "# Let's add a feature representing the CDS\n",
    "dna.features.append(SeqFeature(SimpleLocation(start=6, end=60), type=\"CDS\"))\n",
    "\n",
    "# To design the primer, we extract the template sequence we want to amplify, and use the `primer_design` method.\n",
    "template = dna.features[0].location.extract(dna)\n",
    "\n",
    "# We get an amplicon object (a subclass of Dseqrecord), that also contains extra info\n",
    "# of where the primers align etc.\n",
    "amplicon = primer_design(template, target_tm=60.0, limit=15)\n",
    "\n",
    "# We extract the primers\n",
    "fwd_primer, rvs_primer = amplicon.primers()\n",
    "\n",
    "# We print the Tms\n",
    "print(\"Forward primer Tm:\", tm_default(fwd_primer.seq))\n",
    "print(\"Forward primer sequence:\", fwd_primer.seq)\n",
    "print()\n",
    "print(\"Reverse primer Tm:\", tm_default(rvs_primer.seq))\n",
    "print(\"Reverse primer sequence:\", rvs_primer.seq)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Special primers\n",
    "\n",
    "We saw an example where we simply want to amplify a region of DNA. But what if we want to design primers for a specific restriction enzyme, or for Gibson Assembly? That's also easy."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Restriction enzyme\n",
    "\n",
    "Simply append the sequence you want at the 5' end of the primers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ttGAATTCATGCAAACAGTAATGATGGATGAC\n",
      "ttGAATTCTTATTCAGCAATAGAATCAGTGCTTTG\n"
     ]
    }
   ],
   "source": [
    "from Bio.Restriction import EcoRI\n",
    "fwd_primer_EcoRI = 'ttGAATTC' + fwd_primer\n",
    "# You can also do it like this!\n",
    "rvs_primer_EcoRI = 'tt' + EcoRI.site + rvs_primer\n",
    "\n",
    "print(fwd_primer_EcoRI.seq)\n",
    "print(rvs_primer_EcoRI.seq)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> **Edge case:** Some recognition sites contain ambiguous bases, for instance `Bst4CI` cuts at site `ACNGT`, where `N` can be any nucleotide. In that case, you can use the dictionary provided by biopython (`from Bio.Data.IUPACData import ambiguous_dna_values`) to produce a concrete DNA sequence you can use in real life."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Gibson Assembly\n",
    "\n",
    "To design primers for Gibson Assembly, you can use the `assembly_fragments` function."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Linear Gibson Assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Primers for fragment 1:\n",
      "ATGCAAACAGTAATGATGGATGAC\n",
      "GAGTGATTATCTTTTTCAGCAATAGAATCAGTGC\n",
      "\n",
      "Primers for fragment 2:\n",
      "ATAATCACTCTAATAATGAATCTAACTTTACTTGGAAA\n",
      "ATGCTTTTCCACTTGTTCACG\n",
      "\n",
      "PCR product 1:\n",
      "ATGCAAACAGTAATGATGGATGACATTCAAAGCACTGATTCTATTGCTGAAAAAGATAATCACTC\n",
      "\n",
      "PCR product 2:\n",
      "ATAATCACTCTAATAATGAATCTAACTTTACTTGGAAAGCGTTTCGTGAACAAGTGGAAAAGCAT\n",
      "\n",
      "Overlap\n",
      "ATGCAAACAGTAATGATGGATGACATTCAAAGCACTGATTCTATTGCTGAAAAAGATAATCACTC\n",
      "                                                       ATAATCACTCTAATAATGAATCTAACTTTACTTGGAAAGCGTTTCGTGAACAAGTGGAAAAGCAT\n"
     ]
    }
   ],
   "source": [
    "from pydna.design import assembly_fragments\n",
    "# Let's imagine we want to join these two sequences together linearly with Gibson Assembly\n",
    "seq1 = Dseqrecord('ATGCAAACAGTAATGATGGATGACATTCAAAGCACTGATTCTATTGCTGAAAAAGATAAT')\n",
    "seq2 = Dseqrecord('CACTCTAATAATGAATCTAACTTTACTTGGAAAGCGTTTCGTGAACAAGTGGAAAAGCAT')\n",
    "\n",
    "# First, we design primers for each fragment, as before:\n",
    "pre_amplicon1 = primer_design(seq1, target_tm=60.0, limit=15)\n",
    "pre_amplicon2 = primer_design(seq2, target_tm=60.0, limit=15)\n",
    "\n",
    "# Then, we use the `assembly_fragments` function to design primers for Gibson Assembly\n",
    "amplicon1, amplicon2 = assembly_fragments([pre_amplicon1, pre_amplicon2], overlap=10)\n",
    "\n",
    "# We print the primers:\n",
    "fwd_1, rvs_1 = amplicon1.primers()\n",
    "fwd_2, rvs_2 = amplicon2.primers()\n",
    "\n",
    "print('Primers for fragment 1:')\n",
    "print(fwd_1.seq)\n",
    "print(rvs_1.seq)\n",
    "print()\n",
    "print('Primers for fragment 2:')\n",
    "print(fwd_2.seq)\n",
    "print(rvs_2.seq)\n",
    "print()\n",
    "\n",
    "# The amplicons contain the PCR products (note the overlap between the two fragments)\n",
    "print('PCR product 1:')\n",
    "print(amplicon1.seq)\n",
    "print()\n",
    "print('PCR product 2:')\n",
    "print(amplicon2.seq)\n",
    "print()\n",
    "\n",
    "print('Overlap')\n",
    "print(amplicon1.seq)\n",
    "print(' '*55,amplicon2.seq, sep='')\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once you have the amplicons, you can use `Assembly` to join them together (see the `Gibson` notebook for more details)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "65bp_PCR_prod|10\n",
      "              \\/\n",
      "              /\\\n",
      "              10|65bp_PCR_prod\n",
      "\n",
      "Dseqrecord(-120)\n",
      "\u001b[48;5;11mATGCAAACAGTAATGATGGATGAC\u001b[0mATTCAAAGCACTGATTCTATTGCTGAAAAAGATAATCACTCTAATAATGAATCTAACTTTACTTGGAAAGCGTTTCGTGAACAAGTGGAAAAGCAT\n",
      "TACGTTTGTCATTACTACCTACTGTAAGTTTCGTGACTAAGATAACGACTTTTTCTATTAGTGAGATTATTACTTAGATTGAAATGAACCTTTCGCAAAGCACTTGTTCACCTTTTCGTA\n"
     ]
    }
   ],
   "source": [
    "from pydna.assembly import Assembly\n",
    "from pydna.common_sub_strings import terminal_overlap\n",
    "\n",
    "\n",
    "assembly = Assembly([amplicon1, amplicon2], limit=10, algorithm=terminal_overlap)\n",
    "product = assembly.assemble_linear()[0]\n",
    "print(product.figure())\n",
    "\n",
    "print()\n",
    "\n",
    "print(Dseqrecord(product).figure())\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Circular Gibson Assembly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Primers for fragment 1:\n",
      "AGCATATGCAAACAGTAATGATGGATGAC\n",
      "GAGTGATTATCTTTTTCAGCAATAGAATCAGTGC\n",
      "\n",
      "Primers for fragment 2:\n",
      "ATAATCACTCTAATAATGAATCTAACTTTACTTGGAAA\n",
      "TGCATATGCTTTTCCACTTGTTCACG\n",
      "\n",
      " -|70bp_PCR_prod|10\n",
      "|                \\/\n",
      "|                /\\\n",
      "|                10|70bp_PCR_prod|10\n",
      "|                                 \\/\n",
      "|                                 /\\\n",
      "|                                 10-\n",
      "|                                    |\n",
      " ------------------------------------\n",
      "\n",
      "Dseqrecord(o120)\n",
      "AGCAT\u001b[48;5;11mATGCAAACAGTAATGATGGATGAC\u001b[0mATTCAAAGCACTGATTCTATTGCTGAAAAAGATAATCACTCTAATAATGAATCTAACTTTACTTGGAAAGCGTTTCGTGAACAAGTGGAAA\n",
      "TCGTATACGTTTGTCATTACTACCTACTGTAAGTTTCGTGACTAAGATAACGACTTTTTCTATTAGTGAGATTATTACTTAGATTGAAATGAACCTTTCGCAAAGCACTTGTTCACCTTT\n"
     ]
    }
   ],
   "source": [
    "# We use the `assembly_fragments` function with `circular=True`\n",
    "amplicon1, amplicon2 = assembly_fragments([pre_amplicon1, pre_amplicon2], overlap=10, circular=True)\n",
    "\n",
    "# We print the primers:\n",
    "fwd_1, rvs_1 = amplicon1.primers()\n",
    "fwd_2, rvs_2 = amplicon2.primers()\n",
    "\n",
    "print('Primers for fragment 1:')\n",
    "print(fwd_1.seq)\n",
    "print(rvs_1.seq)\n",
    "print()\n",
    "print('Primers for fragment 2:')\n",
    "print(fwd_2.seq)\n",
    "print(rvs_2.seq)\n",
    "print()\n",
    "\n",
    "assembly = Assembly([amplicon1, amplicon2], limit=10, algorithm=terminal_overlap)\n",
    "\n",
    "# Here we use assemble_circular!\n",
    "product = assembly.assemble_circular()[0]\n",
    "print(product.figure())\n",
    "\n",
    "print()\n",
    "\n",
    "print(Dseqrecord(product).figure())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Adding spacers / linkers to Gibson Assembly primers\n",
    "\n",
    "In this case, you can pass a list of `Amplicons` and `Dseqrecords` to the `assembly_fragments` function. The `amplicons` will be used as the fragments to assemble, and the `dseqrecords` will be used as spacers between the fragments, as long as they are shorter than the argument `maxlink`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "GCATtttATGCAAACAGTAATGATGGATGACATTCAAAGCACTGATTCTATTGCTGAAAAAGATAATaaaCACTCTAATAATGAATCTAACTTTACTTGGAAAGCGTTTCGTGAACAAGTGGAAAA\n",
      "    ^^^                                                            ^^^\n"
     ]
    }
   ],
   "source": [
    "# We create two spacers as dseqrecords\n",
    "spacer1 = Dseqrecord('aaa')\n",
    "spacer2 = Dseqrecord('ttt')\n",
    "\n",
    "amplicon1, amplicon2 = assembly_fragments([pre_amplicon1, spacer1, pre_amplicon2, spacer2], overlap=10, circular=True)\n",
    "\n",
    "assembly = Assembly([amplicon1, amplicon2], limit=10, algorithm=terminal_overlap)\n",
    "\n",
    "# Here we use assemble_circular!\n",
    "product = assembly.assemble_circular()[0]\n",
    "\n",
    "# See the linkers that have been added\n",
    "print()\n",
    "print(Dseqrecord(product).seq)\n",
    "print(4*' ', '^^^', 60*' ', '^^^', sep='')\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary of assembly_fragments behaviour\n",
    "\n",
    "The behaviour is summarised in the following graphics for linear and circular assembly.\n",
    "\n",
    "![assembly_fragments behaviour](../_static/assembly_fragment_slide_linear.png)\n",
    "![assembly_fragments behaviour](../_static/assembly_fragment_slide_circular.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using the NEB Tm calculator for primer design\n",
    "\n",
    "You can use the NEB Tm calculator for primer design, but it is much slower than the other methods because it relies on making requests to the NEB server.\n",
    "\n",
    "You can use the `estimate_function` argument to pass another function that will be used to estimate a first Tm, then the NEB calculator will be used to refine it.\n",
    "\n",
    "The `estimate_function` and `tm_method` only takes the sequence string as argument, so you need to make an adapter function for the the case of `tm_neb` (see below)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With an estimate function\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "\n",
      "Without an estimate function\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n",
      "Estimating Tm using NEB Tm calculator\n"
     ]
    },
    {
     "data": {
      "text/html": [
       "Amplicon(54)"
      ],
      "text/plain": [
       "Amplicon(54)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "\n",
    "def custom_tm_function(seq):\n",
    "    # We include the print statement to show that is only called a few times\n",
    "    print('Estimating Tm using NEB Tm calculator')\n",
    "    return tm_neb(seq, 0.5, 'q5-0')\n",
    "\n",
    "print(\"With an estimate function\")\n",
    "primer_design(template, tm_func=custom_tm_function, estimate_function=tm_default)\n",
    "print()\n",
    "\n",
    "print(\"Without an estimate function\")\n",
    "primer_design(template, tm_func=custom_tm_function)"
   ]
  }
 ],
 "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.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}