docs/notebooks/Dseq.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Representing sequences in pydna\n",
    "> Visit the full library documentation [here](https://pydna-group.github.io/pydna/)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<a target=\"_blank\" href=\"https://colab.research.google.com/github/pydna-group/pydna/blob/master/docs/notebooks/Dseq.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": [
    "Pydna contains classes to represent double stranded DNA sequences that can:\n",
    "\n",
    "* Be linear\n",
    "* Be circular\n",
    "* Contain overhangs (sticky ends).\n",
    "\n",
    "These sequences can be used to simulate molecular biology methods such as cloning and PCR. The main classes used to represent sequences are `Dseq` and `Dseqrecord`.\n",
    "* `Dseq` represents the sequence only. Think of it as a FASTA file.\n",
    "* `Dseqrecord` can contain sequence features and other info such as publication, authors, etc. Think of it as a Genbank file.\n",
    "\n",
    "> NOTE: The `Dseq` class is a subclass of biopython's `Seq`, whose documentation can be found [here](https://biopython.org/wiki/Seq). `Dseqrecord` is a subclass of biopython's `SeqRecord`, whose documentation can be found [here](https://biopython.org/wiki/SeqRecord)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Dseq Class\n",
    "\n",
    "We can create a `Dseq` object in different ways.\n",
    "\n",
    "For a linear sequence without overhangs, we create a `Dseq` object passing a string with the sequence. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-5)\n",
       "aatat\n",
       "ttata"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from pydna.dseq import Dseq\n",
    "my_seq = Dseq(\"aatat\")\n",
    "my_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the console representation above, there are three lines:\n",
    "1. `Dseq(-5)` indicates that the sequence is linear and has 5 basepairs.\n",
    "2. `aatat`, the top / sense / watson strand, referred from now on as **watson** strand..\n",
    "3. `ttata`, the bottom / anti-sense / crick strand, referred from now on as **crick** strand.\n",
    "\n",
    "Now, let's create a circular sequence:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(o5)\n",
       "aatat\n",
       "ttata"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq = Dseq(\"aatat\", circular=True)\n",
    "my_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note how `o5` indicates that the sequence is circular and has 5 basepairs.\n",
    "\n",
    "One way to represent a linear sequence with overhangs is to instantiate `Dseq` with the following arguments:\n",
    "* The `watson` strand as a string in the 5'-3' direction.\n",
    "* The `crick` strand as a string in the 5'-3' direction.\n",
    "* The 5' overhang `ovhg` (overhang), which can be positive or negative, and represents the number of basepairs that the `watson` strand extends beyond the `crick` strand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-5)\n",
       "actag\n",
       " gatc"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dseq(\"actag\", \"ctag\", -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> Note how the bottom strand is passed in the 5'-3' direction, but it is represented in the 3'-5' direction in the console output.\n",
    "\n",
    "If you omit the `ovhg` argument, pydna will try to find the value that makes the `watson` and `crick` strands complementary."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-5)\n",
       "actag\n",
       " gatc"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dseq(\"actag\", \"ctag\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The best way to get a feeling for the meaning of `ovhg` is to visualise the possible scenarios as such:\n",
    "\n",
    "```\n",
    "dsDNA       overhang\n",
    "\n",
    "  nnn...    2\n",
    "nnnnn...\n",
    "\n",
    "  nnnn...   1\n",
    "nnnnn...\n",
    "\n",
    "nnnnn...    0\n",
    "nnnnn...\n",
    "\n",
    "nnnnn...   -1\n",
    "  nnnn...\n",
    "\n",
    "nnnnn...   -2\n",
    "  nnn...\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of note, the DNA sequence can be passed in both lower case and upper case, and are not restricted to the conventional ATCG nucleotides (E.g ), The class supports the IUPAC ambiguous nucleotide code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-5)\n",
       "Actag\n",
       " gatC"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dseq(\"Actag\", \"Ctag\", -1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to pass the overhangs is to use the `from_full_sequence_and_overhangs` classmethod, which only needs the `watson`/sense strand. This is useful you can only store the entire sequence (e.g. in a FASTA file), or if you want to specify overhangs on both sides of the double stranded DNA when you create the object.\n",
    "\n",
    "Both the `watson_ovhg` and `crick_ovhg` can be passed following the same rules as above. Specifically, the `crick_ovhg` argument is identical to the conventional `ovhg` argument. The `watson_ovhg` argument is the `ovhg` argument applied to the reverse complementary sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-8)\n",
       "aaatta\n",
       "   aattt"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq = Dseq.from_full_sequence_and_overhangs(\"aaattaaa\", crick_ovhg=-3, watson_ovhg=-2)\n",
    "my_seq"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A list of possible scenarios, applying positive and negative `crick_ovhg` and `watson_ovhg` to a `Dseq` object are visualised in the output of the code below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "watson_ovhg is -3, crick_ovhg is -2\n",
      "Dseq(-8)\n",
      "aaatt\n",
      "  taattt\n",
      "\n",
      "watson_ovhg is 3, crick_ovhg is -2\n",
      "Dseq(-8)\n",
      "aaattaaa\n",
      "  taa\n",
      "\n",
      "watson_ovhg is -3, crick_ovhg is 2\n",
      "Dseq(-8)\n",
      "  att\n",
      "tttaattt\n",
      "\n",
      "watson_ovhg is 3, crick_ovhg is 2\n",
      "Dseq(-8)\n",
      "  attaaa\n",
      "tttaa\n",
      "\n"
     ]
    }
   ],
   "source": [
    "for crick_ovhg in [-2, 2]:\n",
    "    for watson_ovhg in [-3, 3]:\n",
    "        print(\"watson_ovhg is \" + str(watson_ovhg) + \", crick_ovhg is \" + str(crick_ovhg))\n",
    "        my_seq = Dseq.from_full_sequence_and_overhangs(\"aaattaaa\", crick_ovhg, watson_ovhg)\n",
    "        print(my_seq.__repr__() + \"\\n\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The drawing below can help visualize the meaning of the overhangs.\n",
    "```\n",
    "  (-3)--(-2)--(-1)--(x)--(x)--(x)--(-1)--(-2)\n",
    "\n",
    "5'( a)--( a)--( a)--(t)--(t)--(a)--( a)--( a)3'\n",
    "3'( a)--( a)--( a)--(t)--(t)--(a)--( a)--( a)5'\n",
    "\n",
    "5'( a)--( a)--( a)--(t)--(t)--(a)--(  )--(  )3'\n",
    "3'(  )--(  )--(  )--(t)--(t)--(a)--( a)--( a)5'\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you would like to check the overhangs for a `Dseq` object, it can be done by calling the methods `five_prime_end` and `three_prime_end` to show the 5' and 3' overhangs, respectively. An example of a `Dseq` object, and examples showing what the print-out of the methods looks like are demonstrated here:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dseq(-7)\n",
      "aatat\n",
      "  atatt\n",
      "(\"5'\", 'aa')\n",
      "(\"5'\", 'tt')\n"
     ]
    }
   ],
   "source": [
    "my_seq = Dseq(\"aatat\", \"ttata\", ovhg=-2)\n",
    "print(my_seq.__repr__())\n",
    "print(my_seq.five_prime_end())\n",
    "print(my_seq.three_prime_end())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you now want to join your sequence's sticky ends to make a circular sequence (i.e Plasmid), you can use the `looped` method. The sticky ends must be compatible to do so.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(o5)\n",
       "aatat\n",
       "ttata"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq = Dseq(\"aatat\", \"ttata\", ovhg=-2)\n",
    "my_seq.looped()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you want to change the circular origin of the sequence/plasmid, this can be easily done using the `shifted` method. This can be done by providing the number of bases between the original origin with the new origin: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(o5)\n",
       "tataa\n",
       "atatt"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq = Dseq(\"aatat\", circular=True)\n",
    "my_seq.shifted(2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## __getitem__, __repr__, and  __str__ methods\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Slicing sequences (`__getitem__`)\n",
    "\n",
    "`__getitem__` is the method that is called when you use the square brackets `[]` after a python object. Below is an example of the builtin python `list`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "using square brackets: [2, 3]\n",
      "is the same as using __getitem__: [2, 3]\n"
     ]
    }
   ],
   "source": [
    "my_list = [1, 2, 3]\n",
    "\n",
    "print('using square brackets:', my_list[1:])\n",
    "print('is the same as using __getitem__:', my_list.__getitem__(slice(1, None)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `__getitem__` method is modified in pydna to deal with `Dseq` objects and returns a slice of the `Dseq` object, defined by the a start value and a stop value, similarly to string indexing. In other words, `__getitem__` indexes `Dseq`. Note that '__getitem__' (and, consequently, `[]`) uses zero-based indexing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-3)\n",
       "tat\n",
       "ata"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq = Dseq(\"aatataa\")\n",
    "my_seq[2:5]\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`__getitem__` respects overhangs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-5)\n",
       "tata\n",
       "atatt"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "my_seq = Dseq.from_full_sequence_and_overhangs(\"aatataa\", crick_ovhg=0, watson_ovhg=-1)\n",
    "my_seq[2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that index zero corresponds to the leftmost base of the sequence, which might not necessarily be on the `watson` strand. Let's create a sequence that has an overhang on the left side."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-10)\n",
       "   acgttcc\n",
       "ttatgcaagg"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sequence_with_overhangs = Dseq.from_full_sequence_and_overhangs(\"aatacgttcc\", crick_ovhg=3, watson_ovhg=0)\n",
    "sequence_with_overhangs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we index starting from `2`, we don't start counting on the watson, but on the crick strand since it is the leftmost one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-8)\n",
       " acgttcc\n",
       "atgcaagg"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sequence_with_overhangs[2:]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Slicing circular sequences\n",
    "When slicing circular `Dseq` objects we get linear sequences."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-4)\n",
       "atct\n",
       "taga"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "circular_seq = Dseq(\"aatctaa\", circular=True)\n",
    "circular_seq[1:5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can slice circular sequences across the origin (where index is zero) if the first index is bigger than the second index. This is demonstrated in the example below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-4)\n",
       "aaaa\n",
       "tttt"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "circular_seq[5:2]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Printing sequences to the console: `__repr__` and `__str__`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`__repr__` and `__str__` are methods present in all python classes that return a string representation of an object. `__str__` is called by the `print` function, and `__repr__` is used by the console or notebook output when the object is not assigned to a variable. Below is an example with a `date` object:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> print statement: 2023-08-15\n",
      "> repr: datetime.date(2023, 8, 15)\n",
      "> repr from class method: datetime.date(2023, 8, 15)\n",
      "\n",
      "> console output:\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "datetime.date(2023, 8, 15)"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import datetime\n",
    "\n",
    "my_date = datetime.date(2023, 8, 15)\n",
    "\n",
    "print('> print statement:', my_date)\n",
    "print('> repr:', repr(my_date))\n",
    "print('> repr from class method:', my_date.__repr__())\n",
    "\n",
    "print()\n",
    "print('> console output:')\n",
    "my_date"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In a similar way, `__repr__` and `__str__` methods are used by pydna to represent sequences as strings for different purposes:\n",
    "\n",
    "* `__repr__` is used to make a figure-like representation that shows both strands and the overhangs.\n",
    "* `__str__` is used to return the entire sequence as a string of characters (from the left-most to the right-most base of both strands), the way we would store it in a FASTA file.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "> figure-like representation:\n",
      " Dseq(-8)\n",
      "aaatta\n",
      "   aattt\n",
      "\n",
      "> string representation:\n",
      " aaattaaa\n"
     ]
    }
   ],
   "source": [
    "my_seq = Dseq.from_full_sequence_and_overhangs(\"aaattaaa\", crick_ovhg=-3, watson_ovhg=-2)\n",
    "print('> figure-like representation:\\n', my_seq.__repr__())\n",
    "print()\n",
    "print('> string representation:\\n', my_seq)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that on the string representation, the bases correspond to the entire sequence provided, even when they are only present on either the `watson` or `crick` strand. In the example above, the last two `aa` bases are missing from the `watson` strand,  and that only the `crick` strand has them."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Edge cases\n",
    "\n",
    "You can create arbitrary double-stranded sequences that are not complementary if you specify both strands and an overhang, but you won't be able to use them for molecular biology simulations. For example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Dseq(-6)\n",
       "  xxxx\n",
       "tata"
      ]
     },
     "execution_count": null,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dseq(\"xxxx\", \"atat\", ovhg=2)"
   ]
  }
 ],
 "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
}