{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Unpaired assembly challenge"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You will implement software to assemble a genome from synthetic reads.  We supply Python code snippets that you might use or adapt in your solutions, but you don't have to."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 1: Get and parse the reads\n",
    "\n",
    "**10% of the points for this question**\n",
    "\n",
    "Download the reads:\n",
    "\n",
    "http://www.cs.jhu.edu/~langmea/resources/f2020_hw4_reads.fa\n",
    "\n",
    "All the reads come from the same synthetic genome and each is 100 nt long.  For simplicity, these reads don't have any quality values.\n",
    "\n",
    "The following Python code will download the data to a file called `reads.fa` in the current directory.  (Caveat: I don't think the code below works in Python 3.  Sorry about that.  Go here for details on how to fix: http://python-future.org/compatible_idioms.html#urllib-module.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Download the file containing the reads to \"reads.fa\" in current directory\n",
    "! wget http://www.cs.jhu.edu/~langmea/resources/f2020_hw4_reads.fa\n",
    "\n",
    "# Following line is so we can see the first few lines of the reads file\n",
    "# from within IPython -- don't paste this into your Python code\n",
    "! head f2020_hw4_reads.fa"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or you can download them manually from your browser."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 2: Build an overlap graph"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**40% of the points for this question**\n",
    "\n",
    "Goal: Write a file containing each read's *best buddy to the right*.  Let's define that.\n",
    "\n",
    "For each read $A$, find the other read $B$ that has the longest suffix/prefix match with $A$, i.e. a suffix of $A$ matches a prefix of $B$. $B$ is $A$'s *best buddy* to the right. **However**, if there is a tie, or if the longest suffix/prefix match is less than 40 nucleotides long, then $A$ *has no best buddy to the right*.  For each read, your program should output either (a) nothing, if there is no best buddy to the right, or (b) a single, space-separated line with the IDs of $A$ and $B$ and the length of the overlap, like this:\n",
    "\n",
    "    0255/2 2065/1 88\n",
    "\n",
    "This indicates a **88** bp suffix of the read with ID **0255/2** is a prefix of the read with ID **2065/1**.  Because of how we defined best buddy, it also means no other read besides **2065/1** has a prefix of **88**+ bp that is also a suffix of read **0255/2**.  A corrolary of this is that a particular read ID should appear in the first column of your program's output at most once.  Also, since we require the overlap to be at least 40 bases long, no number less than 40 should every appear in the last column.\n",
    "\n",
    "Notes:\n",
    "* You can assume all reads are error-free and from the forward strand.  You do **not** need to consider sequencing errors or reverse complements.\n",
    "* Below is a hint that can make things speedier.\n",
    "* The order of the output lines is not important.\n",
    "\n",
    "[FASTA parsing code]: http://bit.ly/CG_FASTA\n",
    "[overlap detection code]: http://bit.ly/CG_SuffixPrefix_1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint 1**: the following function groups reads such that you can avoid comparing every read to every other read when looking for suffix/prefix matches.  It builds a dictionary called `table` where the keys are k-mers and the values are sets containing the names of all reads containing that k-mer.  Since you are looking for overlaps of length at least 40, you only need to compare reads if they have at least 1 40-mer in common."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_kmer_table(seqs, k):\n",
    "    ''' Given dictionary (e.g. output of parse_fasta) and integer k,\n",
    "        return a dictionary that maps each k-mer to the set of names\n",
    "        of reads containing the k-mer. '''\n",
    "    table = {}  # maps k-mer to set of names of reads containing k-mer\n",
    "    for name, seq in seqs.items():\n",
    "        for i in range(0, len(seq) - k + 1):\n",
    "            kmer = seq[i:i+k]\n",
    "            if kmer not in table:\n",
    "                table[kmer] = set()\n",
    "            table[kmer].add(name)\n",
    "    return table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint 2**: here's a function for finding suffix/prefix matches; we saw this in class:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def suffixPrefixMatch(str1, str2, min_overlap):\n",
    "    ''' Returns length of longest suffix of str1 that is prefix of\n",
    "        str2, as long as that suffix is at least as long as min_overlap. '''\n",
    "    if len(str2) < min_overlap: return 0\n",
    "    str2_prefix = str2[:min_overlap]\n",
    "    str1_pos = -1\n",
    "    while True:\n",
    "        str1_pos = str1.find(str2_prefix, str1_pos + 1)\n",
    "        if str1_pos == -1: return 0\n",
    "        str1_suffix = str1[str1_pos:]\n",
    "        if str2.startswith(str1_suffix): return len(str1_suffix)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Part 3: Build unitigs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**50% of the points for this question**\n",
    "\n",
    "Goal: Write a program that takes the output of the overlap program from part 1 and creates **uni**quely assemblable con**tigs** (unitigs), using the best buddy algorithm described below.\n",
    "\n",
    "We already determined each read's best buddy *to the right*.  I'll abbreviate this as bbr.  We did not attempt to compute each read's best buddy *to the left* (bbl), but we can infer it from the bbrs.  Consider the following output:\n",
    "\n",
    "    A B 60\n",
    "    E A 40\n",
    "    C B 70\n",
    "    D C 40\n",
    "\n",
    "$A$'s bbr is $B$.  But $B$'s bbl is $C$, *not* $A$!  Your program should form unitigs by joining together two reads $X$ and $Y$ if they are *mutual best buddies*.  $X$ and $Y$ are mutual best buddies if $X$'s bbr is $Y$ **and** $Y$'s bbl is $X$, or vice versa.  In this example, we would join $D$, $C$, and $B$ into a single unitig (and in that order), and would join reads $E$ and $A$ into a single unitig (also in that order).\n",
    "\n",
    "Your program's output should consist of several entries like the following, with one entry per unitig:\n",
    "\n",
    "    START UNITIG 1 D\n",
    "      C 40\n",
    "      B 70\n",
    "    END UNITIG 1\n",
    "    START UNITIG 2 E\n",
    "      A 40\n",
    "    END UNITIG 2\n",
    "\n",
    "The first entry represents a unitig with ID **1** consisting of 3 reads.  The first (leftmost) read is **D**.  The second read, **C**, has a **40** nt prefix that is a suffix of the previous read (**D**).  The third (rightmost) read in the contig (**B**) has a **70** bp prefix that is a suffix of the previous read (**C**).\n",
    "\n",
    "Each read should be contained in exactly one unitig.  The order of unitigs in the file is not important, but the unitig IDs should be integers and assigned in ascending order.\n",
    "\n",
    "Note: we will never provide an input that can result in a circular unitig (i.e. one where a chain of mutual best buddies loops back on itself.)\n",
    "\n",
    "**Hint: the correct solution consists of exactly 4 unitigs**."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## OPTIONAL Part 4: Finish the assembly\n",
    "\n",
    "**This part is _optional_.  You can submit your solution if you like.  No extra credit will be awarded.**\n",
    "\n",
    "Goal: Assemble the genome!  Report the sequence of the original genome as a FASTA file.\n",
    "\n",
    "This requires that you compare the unitigs to each other, think about what order they must go in, and then put them together accordingly.  Submit your solution as a single FASTA file containing a single sequence named \"`solution`\".  The FASTA file should be \"wrapped\" so that no line has more than 60 characters.  You can use the following Python code to write out your answer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sys\n",
    "\n",
    "def write_solution(genome, per_line=60, out=sys.stdout):\n",
    "    offset = 0\n",
    "    out.write('>solution\\n')\n",
    "    while offset < len(genome):\n",
    "        nchars = min(len(genome) - offset, per_line)\n",
    "        line = genome[offset:offset+nchars]\n",
    "        offset += nchars\n",
    "        out.write(line + '\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here is an example of what your output should look like.  Note how the sequence is spread over many lines."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import random\n",
    "random.seed(5234)\n",
    "write_solution(''.join([random.choice('ACGT') for _ in range(500)]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Hint 1: the correct genome is 7959 nucleotides long**\n",
    "\n",
    "**Hint 2: to learn how the unitigs should go together, try overlapping them with each other**\n",
    "\n",
    "**Hint 3: a unitig might be used more than once in the solution**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.8.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}