{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3e44c454-8135-48bb-a5e0-6c35f493fe7d",
   "metadata": {},
   "source": [
    "### Tutorial A1: Sequence Manipulation\n",
    "\n",
    "Being able to manipulate sequences is central to analyzing what machine learning models have learned because they allow you to easily ask hypotheticals. Many common downstream analyses, such as <i>in silico marginalizations</i> and <i>feature attributions</i> rely on sequence manipulations to handle motif insertion or construct background sequences, respectively. But, these forms of manipulations are also invaluable tools for dissecting individual regions and doing sequence design. For instance, given a locus of interest, one can shuffle regions (such as known motifs) and see what impact the shuffling has on the prediction. Alternatively, one can iterative insert motifs into a sequence to see whether a predictive model yields a desired result.\n",
    "\n",
    "Given the foundational nature of these operations, they are central to tangermeme and are implemented in `tangermeme.ersatz`. Why `ersatz`?\n",
    "\n",
    "<img src=\"imgs/ersatz.png\" width=500>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7e79b62e-e8a2-4f8a-a06e-4dfa519eefa0",
   "metadata": {},
   "source": [
    "Accordingly, \"ersatz sequences\" in this context are those that are not natural sequences -- either because they are alterations of a native sequence, or they are fully synthetic.\n",
    "\n",
    "#### Substitutions\n",
    "\n",
    "The simplest operation is the substitution, where some positions are switched out for others. Consider the following where an AP-1 motif (https://jaspar.elixir.no/matrix/MA1144.1/) is inserted into uniformly random nucleotides. Colloquially, \"substitutions\" are sometimes called \"insertions\", and while there are similarities between the two they are not exactly the same thing. Let's begin by randomly generating one sequence. Following the PyTorch format, the one-hot encoded sequences should have the shape `(batch_size, alphabet_size, sequence_length)`. Generation can be made deterministic using the `random_state` parameter.\n",
    "\n",
    "As a note, all of the operations in `tangermeme.ersatz` happen on one-hot encoded tensors but, for simplicity of seeing the effect of each operation in this tutorial, we will convert these tensors back to characters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "19cd9bf8-ba98-49fb-9487-619b20810a03",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATCATTTTCTCGATGAAAGC'"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.utils import random_one_hot\n",
    "from tangermeme.utils import characters\n",
    "\n",
    "X = random_one_hot((1, 4, 20), random_state=0)\n",
    "characters(X[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08c4b640-b6ff-40fc-8fb4-6d4d22e15874",
   "metadata": {},
   "source": [
    "Now, let's put the motif in. We can either pass in a string, in which case the same motif is added in to each sequence in the batch at the same position, or we can pass in a one-hot encoding. If the one-hot encoding has shape `(1, alphabet_size, motif_size)`, the same motif is added to each sequence similarly to if a string were used, but if the one-hot encoding has the shape `(batch_size, alphabet_size, motif_size)` then the i-th sequence in the batch has the i-th one-hot encoding substituted in."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "6dd5abe9-d2a4-4311-9864-9b437986f31d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATCATTTTGACTCAGAAAGC'"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import substitute\n",
    "\n",
    "X_ap1 = substitute(X, \"TGACTCA\")\n",
    "characters(X_ap1[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a988f7a0-9b79-4ed8-b643-d2edcc978a66",
   "metadata": {},
   "source": [
    "By default, the substitution will happen in the middle of the sequence. If you'd like to control where it happens you can pass in a parameter `start` with the index to start the substitution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "8272f60d-5ce1-4f63-bdc9-2fd1d5c362ee",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATTGACTCATCGATGAAAGC'"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_ap1 = substitute(X, \"TGACTCA\", start=2)\n",
    "characters(X_ap1[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b51ad7cf-3ab3-4c54-b9c6-c7dabd49b41c",
   "metadata": {},
   "source": [
    "If we have a one-hot encoding rather than a string sequence, we can pass that in instead without needing to first convert it back to characters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "7234ccc5-991a-422f-bf81-b8663451645b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATCATTTTGACTCAGAAAGC'"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.utils import one_hot_encode\n",
    "\n",
    "motif = one_hot_encode(\"TGACTCA\").unsqueeze(0)\n",
    "\n",
    "X_ap1 = substitute(X, motif)\n",
    "characters(X_ap1[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "30726f16-01d0-4d05-a56e-b51f4bd28651",
   "metadata": {},
   "source": [
    "#### Multisubstitutions\n",
    "\n",
    "Sometimes one would like to make multiple substitutions in the same sequence given some spacing between substitutions. Although this could be achieved by calling `substitute` multiple times, we can provide a convenient wrapper for this with the `multisubstitue` function. This function has a very similar signature to `substitute` except that it takes in a list of motifs and the spacing between them. We can try it out first given no spacing between the two motifs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "8c1095eb-3323-4b29-9b18-a334ac345ea3",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATTGACTCATGACTCAAAGC'"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import multisubstitute\n",
    "\n",
    "X_ap12 = multisubstitute(X, [\"TGACTCA\", \"TGACTCA\"], 0, start=2)\n",
    "characters(X_ap12[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ff7ff6ca-5516-4857-9d57-f023836a18af",
   "metadata": {},
   "source": [
    "Now, let's add a little bit of spacing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "2f953c85-9f7d-4e90-b82a-b26bd2c44ed2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATTGACTCATCTGACTCAGC'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_ap12 = multisubstitute(X, [\"TGACTCA\", \"TGACTCA\"], 2, start=2)\n",
    "characters(X_ap12[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8dfb4382-54bc-4853-a493-725d75d98306",
   "metadata": {},
   "source": [
    "Finally, if we have more than two motifs we can optionally provide spacing values between each set of motifs. Note that if we keep the spacing value as an integer but provide more than two motifs that the same spacing is used between all motif pairs."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "9c6c62ae-b670-4ff2-a688-6eb836af6c4e",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATTGATCACTTGATTGACGC'"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_ap12 = multisubstitute(X, [\"TGA\", \"TCA\", \"TGA\", \"TGAC\"], [0, 2, 1], start=2)\n",
    "characters(X_ap12[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "77aee1eb-7e88-48c6-adf8-a217034830d5",
   "metadata": {},
   "source": [
    "#### Insertions\n",
    "\n",
    "Related to substitutions are insertions. As mentioned before, sometimes people say \"insertions\" when what they mean are \"substitutions\", but insertions involve adding the new sequence without modifying or deleting any of the existing sequence. Essentially, the returned sequence will be longer than the original sequence because it now additionally contains the sequence being added. In contrast, substitutions preserve the length of the sequence because characters are explicitly being changed from those in the original sequence to those in the new sequence. Let's see that in action."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "31393491-0aff-4895-988b-7660c1e42228",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATCATTTTCTTGACTCACGATGAAAGC'"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import insert\n",
    "\n",
    "X_ap1 = insert(X, \"TGACTCA\")\n",
    "characters(X_ap1[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "a9a31ff3-9703-42b8-85f7-a55a4ef919ef",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 7, 27)"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape[-1], len(\"TGACTCA\"), len(characters(X_ap1[0]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9ed664ac-1d32-48e0-815c-bffd61d2cd5d",
   "metadata": {},
   "source": [
    "#### Deletion\n",
    "\n",
    "In direct contrast to insertions are deletions: insertions involve adding new sequence to an existing sequence and deletions involve deleting existing sequence and not replacing it with anything. Accordingly, one only passes in sequences and the coordinates to delete instead of passing in any form of new sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "22a4c5be-6003-4480-8bac-56f86b804c51",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('ATCATTTTCTCGATGAAAGC', 'TTTCTCGATGAAAGC')"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import delete\n",
    "\n",
    "X_del = delete(X, start=0, end=5)\n",
    "characters(X[0]), characters(X_del[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ed04384-94e5-4bf7-923c-de9309d89ee4",
   "metadata": {},
   "source": [
    "Next, let's try deleting a portion from the middle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "bdca6e63-5e82-49af-893e-ae818a5bfbf9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'ATCATTTTCTAAAGC'"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_del = delete(X, start=10, end=15)\n",
    "characters(X_del[0])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "450120b2-4100-4346-8953-400cf13f0958",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(20, 15)"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X.shape[-1], X_del.shape[-1]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3f2f885d-c517-45f6-819b-1cecac49ea4c",
   "metadata": {},
   "source": [
    "#### Randomize\n",
    "\n",
    "Deleting positions is one way that we can remove information from a sequence. However, this can pose some issues -- both practically, in terms of needing a sequence of the same length when using machine learning models, and conceptually, in that removing a motif entirely isn't really a biologically plausible alternative to observed sequence. An alternative approach to deleting positions is to replace those positions with randomly generated characters. This would keep the sequence the same length but remove the motif.\n",
    "\n",
    "Here, we replace the first five positions with randomly generated characters and keep the remaining characters the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9dc5f6fc-173f-4795-8618-6030afb3feea",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('ATCATTTTCTCGATGAAAGC', 'GGGGCTTTCTCGATGAAAGC')"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import randomize\n",
    "\n",
    "X_rand = randomize(X, start=0, end=5, random_state=0)\n",
    "characters(X[0]), characters(X_rand[0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6916e27-2d09-4b03-b61f-eb4b3083e222",
   "metadata": {},
   "source": [
    "Frequently, when using randomly generated sequences, one wishes to generate many randomizations so that one can average over the randomness induces by the sequences. To make this easy, `tangermeme` allows you to pass in a parameter `n` specifying the number of randomizations to perform, and returns a tensor with one more dimension than the original tensor whenever randomness is involved. Specifically, the returned tensor will have shape `(n_orig_sequences, n_shuffles, alphabet_size, seq_len)` so that you can shuffle each of many sequence many times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "f236aa06-0b65-4021-8519-1c6faa75ce5c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GGGGCTTTCTCGATGAAAGC\n",
      "GCTTCTTTCTCGATGAAAGC\n",
      "TGGTATTTCTCGATGAAAGC\n",
      "AATTTTTTCTCGATGAAAGC\n",
      "TTCTATTTCTCGATGAAAGC\n",
      "GATGCTTTCTCGATGAAAGC\n",
      "CTCGATTTCTCGATGAAAGC\n",
      "GGGTGTTTCTCGATGAAAGC\n",
      "CCGAGTTTCTCGATGAAAGC\n",
      "GAACCTTTCTCGATGAAAGC\n"
     ]
    }
   ],
   "source": [
    "X_rand = randomize(X, start=0, end=5, n=10, random_state=0)\n",
    "\n",
    "for i in range(10):\n",
    "    print(characters(X_rand[0, i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5040478-8a96-4f6d-a2e1-5f5d594d1761",
   "metadata": {},
   "source": [
    "#### Shuffle\n",
    "\n",
    "A problem with independently generating sequence at each position is that the sampled sequences might have unrealistic compositions. For instance, when you use uniformly randomly generated sequences, the GC content is fairly high compared to naturally occuring sequences. When trying to create backgrounds for specific loci, you might prefer to instead shuffle the positions to ensure that the composition of these backgrounds are the same, while any types of motif are disrupted. \n",
    "\n",
    "We can use the `shuffle` function to completely shuffle a batch of sequences. Each returned sequence will have the same composition as the original sequence that was shuffled, but a different ordering of the elements. In the genomics setting, this means that the same number of each nucleotide will be present but motifs will likely be removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "b124d45b-66b0-4b32-9335-53e0819f0ea9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('ATCATTTTCTCGATGAAAGC', 'GTCCCATTTCTGTTAGAAAA')"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import shuffle\n",
    "\n",
    "X_shuf = shuffle(X, random_state=0)\n",
    "characters(X[0]), characters(X_shuf[0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "73926f4a-a437-42db-a7f5-0fcf022a90d4",
   "metadata": {},
   "source": [
    "Similar to the `randomize` function, if we want to shuffle a sequence many times, we can use the parameter `n`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "727ebcea-dee9-488e-954f-212460a7876a",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GTCCCATTTCTGTTAGAAAA\n",
      "GTGACACACAATACTTTGTT\n",
      "ATATGCCTAATCAGTTATGC\n",
      "GATCAATAGAGCTATCCTTT\n",
      "TTCCTGCAAATTGTCTAAGA\n",
      "ATCCATTCGGGACTTTATAA\n",
      "GGTTTACTAACGCCAATTAT\n",
      "TGTAAGTACCCTATTCTGAA\n",
      "TTTTTATACAGCAGAGACTC\n",
      "CAAACTTTCGCTTAGTATAG\n"
     ]
    }
   ],
   "source": [
    "X_shuf = shuffle(X, n=10, random_state=0)\n",
    "\n",
    "for i in range(10):\n",
    "    print(characters(X_shuf[0, i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "c1d1fa2c-4c7c-4b13-abe3-49cff75eec6e",
   "metadata": {},
   "source": [
    "Furthermore, we can restrict our shuffling to only a portion of the sequence. This can be valuable if you want to knock out a portion of the sequence, such as a known motif or broader regulatory region. All you need to do is specify the start (inclusive) and end (not inclusive) positions that the shuffling should occur."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2aa873ad-80bf-4fdc-a268-9a376c88c343",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('ATCATTTTCTCGATGAAAGC', 'ATCATCTTTGGATCTAAAGC')"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_shuf = shuffle(X, start=5, end=15, random_state=3)\n",
    "characters(X[0]), characters(X_shuf[0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c096324-1742-4982-b7fe-7cc8356f7f63",
   "metadata": {},
   "source": [
    "#### Dinucleotide Shuffle\n",
    "\n",
    "In the genomics setting the `CG` dinucleotide plays an outsized role compared to other dinucleotides and so is significantly underrepresented in the genome. Because normal shuffling will disrupt dinucleotide content, and hence change the proportion of CGs in the sequence, sometimes one wants to use a shuffling strategy that explicitly preserves dinucleotide content.\n",
    "\n",
    "In `tangermeme`, the `dinucleotide_shuffle` function operates similarly to the `shuffle` function. For instance, we can shuffle entire sequences:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "d915d52b-a5d0-490c-bff5-d4295bf4116d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('ATCATTTTCTCGATGAAAGC', 'ATCATTTCTCGATTGAAAGC')"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import dinucleotide_shuffle\n",
    "\n",
    "X_shuf = dinucleotide_shuffle(X, random_state=0)\n",
    "characters(X[0]), characters(X_shuf[0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc74b745-4312-4b57-90e4-0fa2aec2e059",
   "metadata": {},
   "source": [
    "We can generate many shuffles using the `n` parameter."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "0c176901-5fe8-4193-b73b-33478681d0dd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ATCATTTCTCGATTGAAAGC\n",
      "ATCTCAAATTTTCGATGAGC\n",
      "ATTTTCTCAATCGAATGAGC\n",
      "AAATTTTCTCATCGATGAGC\n",
      "ATCTTCATCGATTTGAAAGC\n",
      "ATCAATTCTTTCGAATGAGC\n",
      "AATTCTCAATTCGATTGAGC\n",
      "AAATCTCATTTCGATTGAGC\n",
      "ATCATTTCTTCGAAATGAGC\n",
      "ATCTTTCAATTCGAATGAGC\n"
     ]
    }
   ],
   "source": [
    "X_shuf = dinucleotide_shuffle(X, n=10, random_state=0)\n",
    "\n",
    "for i in range(10):\n",
    "    print(characters(X_shuf[0, i]))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2880509-302e-4ff0-8aec-921e959d1c4f",
   "metadata": {},
   "source": [
    "And we can shuffle only a portion of the sequence."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "746cfb32-32e5-498e-aa4a-a95de10f39ec",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "('ATCATTTTCTCGATGAAAGC', 'ATCATTCTTTCGATGAAAGC')"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X_shuf = dinucleotide_shuffle(X, start=5, end=15, random_state=3)\n",
    "characters(X[0]), characters(X_shuf[0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4e79e546-b6cb-412f-ac5a-66dbbfef48a8",
   "metadata": {},
   "source": [
    "As a note, the strategy for doing dinucleotide shuffling that is implemented will always keep the first and last nucleotides the same. Depending on the sequence composition and the length of the region being shuffled, it can be impossible to produce new dinucleotide shuffled sequences. Passing in `verbose` will raise a warning when at least one position (other than the first and last positions) are always the same character. Regardless of the value of `verbose`, an error will be raised when all returned sequences are identical."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "7855ad87-0202-427a-a808-472d0de72cd5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Warning: At least one position in dinucleotide shuffle is identical across all positions.\n"
     ]
    },
    {
     "ename": "ValueError",
     "evalue": "All dinucleotide shuffles yield identical sequences, potentially due to a lack of diversity in sequence.",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mValueError\u001b[0m                                Traceback (most recent call last)",
      "Cell \u001b[0;32mIn[21], line 1\u001b[0m\n\u001b[0;32m----> 1\u001b[0m X_shuf \u001b[38;5;241m=\u001b[39m \u001b[43mdinucleotide_shuffle\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstart\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m10\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mend\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m15\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mrandom_state\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m0\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mn\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;241;43m100\u001b[39;49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mverbose\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[38;5;28;43;01mTrue\u001b[39;49;00m\u001b[43m)\u001b[49m\n\u001b[1;32m      2\u001b[0m characters(X[\u001b[38;5;241m0\u001b[39m]), characters(X_shuf[\u001b[38;5;241m0\u001b[39m, \u001b[38;5;241m0\u001b[39m])\n",
      "File \u001b[0;32m/shared/aemurphy/tangermeme/tangermeme/ersatz.py:612\u001b[0m, in \u001b[0;36mdinucleotide_shuffle\u001b[0;34m(X, start, end, n, random_state, verbose)\u001b[0m\n\u001b[1;32m    610\u001b[0m X_shufs \u001b[38;5;241m=\u001b[39m []\n\u001b[1;32m    611\u001b[0m \u001b[38;5;28;01mfor\u001b[39;00m i \u001b[38;5;129;01min\u001b[39;00m \u001b[38;5;28mrange\u001b[39m(X\u001b[38;5;241m.\u001b[39mshape[\u001b[38;5;241m0\u001b[39m]):\n\u001b[0;32m--> 612\u001b[0m \tinsert_ \u001b[38;5;241m=\u001b[39m \u001b[43m_dinucleotide_shuffle\u001b[49m\u001b[43m(\u001b[49m\u001b[43mX\u001b[49m\u001b[43m[\u001b[49m\u001b[43mi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43m:\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mstart\u001b[49m\u001b[43m:\u001b[49m\u001b[43mend\u001b[49m\u001b[43m]\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mn_shuffles\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mn\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\n\u001b[1;32m    613\u001b[0m \u001b[43m\t\t\u001b[49m\u001b[43mrandom_state\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mrandom_state\u001b[49m\u001b[38;5;241;43m+\u001b[39;49m\u001b[43mi\u001b[49m\u001b[43m,\u001b[49m\u001b[43m \u001b[49m\u001b[43mverbose\u001b[49m\u001b[38;5;241;43m=\u001b[39;49m\u001b[43mverbose\u001b[49m\u001b[43m)\u001b[49m\n\u001b[1;32m    615\u001b[0m \tX_shuf \u001b[38;5;241m=\u001b[39m torch\u001b[38;5;241m.\u001b[39mclone(X[i:i\u001b[38;5;241m+\u001b[39m\u001b[38;5;241m1\u001b[39m])\u001b[38;5;241m.\u001b[39mrepeat(n, \u001b[38;5;241m1\u001b[39m, \u001b[38;5;241m1\u001b[39m)\n\u001b[1;32m    616\u001b[0m \tX_shuf[:, :, start:end] \u001b[38;5;241m=\u001b[39m insert_\n",
      "File \u001b[0;32m/shared/aemurphy/tangermeme/tangermeme/ersatz.py:552\u001b[0m, in \u001b[0;36m_dinucleotide_shuffle\u001b[0;34m(X, n_shuffles, random_state, verbose)\u001b[0m\n\u001b[1;32m    549\u001b[0m \t\t\u001b[38;5;28mprint\u001b[39m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mWarning: At least one position in dinucleotide shuffle \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m\n\u001b[1;32m    550\u001b[0m \t\t\t\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mis identical across all positions.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m    551\u001b[0m \u001b[38;5;28;01mif\u001b[39;00m conserved\u001b[38;5;241m.\u001b[39mmax(dim\u001b[38;5;241m=\u001b[39m\u001b[38;5;241m0\u001b[39m)\u001b[38;5;241m.\u001b[39mvalues\u001b[38;5;241m.\u001b[39mmin() \u001b[38;5;241m==\u001b[39m n_shuffles \u001b[38;5;129;01mand\u001b[39;00m n_shuffles \u001b[38;5;241m>\u001b[39m \u001b[38;5;241m1\u001b[39m:\n\u001b[0;32m--> 552\u001b[0m \t\u001b[38;5;28;01mraise\u001b[39;00m \u001b[38;5;167;01mValueError\u001b[39;00m(\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124mAll dinucleotide shuffles yield identical \u001b[39m\u001b[38;5;124m\"\u001b[39m \u001b[38;5;241m+\u001b[39m\n\u001b[1;32m    553\u001b[0m \t\t\u001b[38;5;124m\"\u001b[39m\u001b[38;5;124msequences, potentially due to a lack of diversity in sequence.\u001b[39m\u001b[38;5;124m\"\u001b[39m)\n\u001b[1;32m    555\u001b[0m \u001b[38;5;28;01mreturn\u001b[39;00m shuffled_sequences\n",
      "\u001b[0;31mValueError\u001b[0m: All dinucleotide shuffles yield identical sequences, potentially due to a lack of diversity in sequence."
     ]
    }
   ],
   "source": [
    "X_shuf = dinucleotide_shuffle(X, start=10, end=15, random_state=0, n=100, verbose=True)\n",
    "characters(X[0]), characters(X_shuf[0, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e40ff365",
   "metadata": {},
   "source": [
    "#### Reverse Complement\n",
    "\n",
    "A common approach to training/testing genomic deep learning models is to use the reverse complement of a DNA sequence.\n",
    "\n",
    "In `tangermeme` this can be done with character represnetations of DNA or with one-hot encodings using `reverse_complement`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "c4ce4974",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'CGAT'"
      ]
     },
     "execution_count": 22,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from tangermeme.ersatz import reverse_complement\n",
    "\n",
    "reverse_complement(\"ATCG\",ohe=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f68410a6",
   "metadata": {},
   "source": [
    "In the same manner as `one_hot_encode`, `reverse_complement` is designed to be used on one sequence (Tensor shape (alphabet length,sequence length))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "10d8f65d",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "torch.Size([4, 7])\n",
      "torch.Size([4, 7])\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'TGAGTCA'"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "motif = one_hot_encode(\"TGACTCA\")\n",
    "print(motif.shape)\n",
    "rev_comp_motif = reverse_complement(motif)\n",
    "print(rev_comp_motif.shape)\n",
    "characters(rev_comp_motif)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b405a6b7",
   "metadata": {},
   "source": [
    "This function can also handle missing values in DNA sequences, represented by 'N' or by `[0,0,0,0]`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "f3cdf0bd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'NCGAT'"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "reverse_complement(\"ATCGN\",ohe=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "822ef359",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'NTGAGTCA'"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "motif = one_hot_encode(\"TGACTCAN\")\n",
    "characters(reverse_complement(motif),allow_N=True)"
   ]
  }
 ],
 "metadata": {
  "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.12.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}