Notebooks/Methods-4-Mapping.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\"><ul class=\"toc-item\"><li><span><a href=\"#Iterative-vs-fragment-based-mapping\" data-toc-modified-id=\"Iterative-vs-fragment-based-mapping-1\">Iterative vs fragment-based mapping</a></span><ul class=\"toc-item\"><li><span><a href=\"#Advantages-of-iterative-mapping\" data-toc-modified-id=\"Advantages-of-iterative-mapping-1.1\">Advantages of iterative mapping</a></span></li><li><span><a href=\"#Advantages-of-fragment-based-mapping\" data-toc-modified-id=\"Advantages-of-fragment-based-mapping-1.2\">Advantages of fragment-based mapping</a></span></li></ul></li><li><span><a href=\"#Mapping\" data-toc-modified-id=\"Mapping-2\">Mapping</a></span><ul class=\"toc-item\"><li><span><a href=\"#Iterative-mapping\" data-toc-modified-id=\"Iterative-mapping-2.1\">Iterative mapping</a></span></li><li><span><a href=\"#Fragment-based-mapping\" data-toc-modified-id=\"Fragment-based-mapping-2.2\">Fragment-based mapping</a></span></li></ul></li><li><span><a href=\"#References\" data-toc-modified-id=\"References-3\">References</a></span></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Iterative vs fragment-based mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Iterative mapping first proposed by <a name=\"ref-1\"/>[(Imakaev et al., 2012)](#cite-Imakaev2012a), allows to map usually a high number of reads. However other methodologies, less \"brute-force\" can be used to take into account the chimeric nature of the Hi-C reads.\n",
    "\n",
    "A simple alternative is to allow split mapping.\n",
    "\n",
    "Another way consists in _pre-truncating_ <a name=\"ref-1\"/>[(Ay and Noble, 2015)](#cite-Ay2015a) reads that contain a ligation site and map only the longest part of the read <a name=\"ref-2\"/>[(Wingett et al., 2015)](#cite-Wingett2015).\n",
    "\n",
    "Finally, an intermediate approach, _fragment-based_, consists in mapping full length reads first, and than splitting unmapped reads at the ligation sites <a name=\"ref-1\"/>[(Serra et al. 2017)](#cite-Serra2017)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![mapping strategies](images/mapping.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advantages of iterative mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- It's the only solution when no restriction enzyme has been used (i.e. micro-C)\n",
    "- Can be faster when few windows (2 or 3) are used"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Advantages of fragment-based mapping "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Generally faster\n",
    "- Safer: mapped reads are generally larger than 25-30 nm (the largest window used in iterative mapping). Less reads are mapped, but the difference is usually canceled or reversed when looking for \"valid-pairs\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "_Note:_ We use __GEM2__ <a name=\"ref-1\"/>[(Marco-Sola et al. 2012)](#cite-Marco-Sola2012), performance are very similar to Bowtie2, and in some cases slightly better. \n",
    "\n",
    "_For now TADbit is only compatible with GEM2._"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pytadbit.mapping.full_mapper import full_mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The full mapping function can be used to perform either iterative or fragment-based mapping, or a combination of both.\n",
    "\n",
    "__It is important to note__ that although the default mapping parameters used by TADbit are relatively strict, __a non negligible proportion of the reads will be mis-mapped__, and this applies __at each iteration of the mapping.__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Iterative mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here an example of use as iterative mapping:\n",
    "(Estimated time 15h with 8 cores)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = 'mouse_B'  # or mouse_PSC\n",
    "rep = 'rep1'  # or rep2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "_Note: the execution of this notebook should be repeated for each of the 4 replicates_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir -p results/iterativ/$cell\\_$rep\n",
    "! mkdir -p results/iterativ/$cell\\_$rep/01_mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preparing FASTQ file\n",
      "  - conversion to MAP format\n",
      "  - trimming reads 1-25\n",
      "Mapping reads in window 1-25...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_PSC_rep1/01_mapping/mapped_mouse_PSC_rep1_r1_tmp/mouse_PSC_rep1_1.fastq_khgGNE\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_PSC_rep1/01_mapping/mapped_mouse_PSC_rep1_r1_tmp/mouse_PSC_rep1_1.fastq_khgGNE -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_PSC_rep1/01_mapping/mapped_mouse_PSC_rep1_r1_tmp/mouse_PSC_rep1_1.fastq_khgGNE_full_1-25\n",
      "Parsing result...\n"
     ]
    }
   ],
   "source": [
    "# for the first side of the reads\n",
    "full_mapping(gem_index_path='genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem',\n",
    "             out_map_dir='results/iterativ/{0}_{1}/01_mapping/mapped_{0}_{1}_r1/'.format(cell, rep),\n",
    "             fastq_path='FASTQs/%s_%s_1.fastq.dsrc' % (cell,rep),\n",
    "             frag_map=False,  clean=True, nthreads=8,\n",
    "             windows=((1,25),(1,35),(1,45),(1,55),(1,65),(1,75)),\n",
    "             temp_dir='results/iterativ/{0}_{1}/01_mapping/mapped_{0}_{1}_r1_tmp/'.format(cell, rep))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And for the second side of the read-end:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preparing FASTQ file\n",
      "  - conversion to MAP format\n",
      "  - trimming reads 1-25\n",
      "Mapping reads in window 1-25...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_rErDYK\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_rErDYK -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_rErDYK_full_1-25\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_rErDYK\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_rErDYK_full_1-25.map\n",
      "Preparing MAP file\n",
      "  - trimming reads 1-35\n",
      "   x removing original input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_rErDYK_filt_1-25.map\n",
      "Mapping reads in window 1-35...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_14i9tR\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_14i9tR -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_14i9tR_full_1-35\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_14i9tR\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_14i9tR_full_1-35.map\n",
      "Preparing MAP file\n",
      "  - trimming reads 1-45\n",
      "   x removing original input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_14i9tR_filt_1-35.map\n",
      "Mapping reads in window 1-45...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_xiEdeg\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_xiEdeg -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_xiEdeg_full_1-45\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_xiEdeg\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_xiEdeg_full_1-45.map\n",
      "Preparing MAP file\n",
      "  - trimming reads 1-55\n",
      "   x removing original input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_xiEdeg_filt_1-45.map\n",
      "Mapping reads in window 1-55...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_EYDFyI\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_EYDFyI -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_EYDFyI_full_1-55\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_EYDFyI\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_EYDFyI_full_1-55.map\n",
      "Preparing MAP file\n",
      "  - trimming reads 1-65\n",
      "   x removing original input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_EYDFyI_filt_1-55.map\n",
      "Mapping reads in window 1-65...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_pe6yQw\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_pe6yQw -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_pe6yQw_full_1-65\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_pe6yQw\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_pe6yQw_full_1-65.map\n",
      "Preparing MAP file\n",
      "  - trimming reads 1-75\n",
      "   x removing original input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_pe6yQw_filt_1-65.map\n",
      "Mapping reads in window 1-75...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_MkfcZ5\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_MkfcZ5 -o /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_MkfcZ5_full_1-75\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_MkfcZ5\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_MkfcZ5_full_1-75.map\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-25.map',\n",
       " 'results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-35.map',\n",
       " 'results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-45.map',\n",
       " 'results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-55.map',\n",
       " 'results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-65.map',\n",
       " 'results/iterativ/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-75.map']"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# for the second side of the reads\n",
    "full_mapping(gem_index_path='genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem',\n",
    "             out_map_dir='results/iterativ/{0}_{1}/01_mapping/mapped_{0}_{1}_r2/'.format(cell, rep),\n",
    "             fastq_path='FASTQs/%s_%s_2.fastq.dsrc' % (cell,rep),\n",
    "             frag_map=False,  clean=True, nthreads=8,\n",
    "             windows=((1,25),(1,35),(1,45),(1,55),(1,65),(1,75)), \n",
    "             temp_dir='results/iterativ/{0}_{1}/01_mapping/mapped_{0}_{1}_r2_tmp/'.format(cell, rep))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fragment-based mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The fragment-based mapping strategy works in 2 steps: \n",
    "1. The read-ends are mapped entirely, assuming that no ligation occurred in them. \n",
    "2. For unmapped read-ends, the function searches for a ligation site (e.g. in the case of MboI this would correspond to `GATCGATC` and in the case of HindIII to `AAGCTAGCTT`). The read-end is split accordingly replacing the ligation site by two RE sites:\n",
    "\n",
    "```\n",
    "                        read-end-part-one---AAGCTAGCTT----read-end-part-two\n",
    "     will be split in:\n",
    "                        read-end-part-one---AAGCTT\n",
    "     and:\n",
    "                                                AAGCTT----read-end-part-two\n",
    "```\n",
    "\n",
    "_Note: __if no ligation site is found__, step two will be repeated using digested RE site as split point (`AAGCT` in the case of HindIII). This is done in order to be protected against sequencing errors. When this path is followed the digested RE site is removed, but not replaced. __If digested RE sites are not found either, the read will be classified as unmapped__._\n",
    "\n",
    "_Note: __both mapping strategies can be combined__, for example defining the windows as previously (iterative mapping), but also give a RE name_ `r_enz=MboI` _and setting_ `frag_map=True` _like this if a read has not been mapped in any window, TADbit will also try to apply the fragment-based strategy._"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "! mkdir -p results/fragment/$cell\\_$rep\n",
    "! mkdir -p results/fragment/$cell\\_$rep/01_mapping"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preparing FASTQ file\n",
      "  - conversion to MAP format\n",
      "Mapping reads in window 1-end...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN -o /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN_full_1-end\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN_full_1-end.map\n",
      "  - splitting into restriction enzyme (RE) fragments using ligation sites\n",
      "  - ligation sites are replaced by RE sites to match the reference genome\n",
      "    * enzymes: MboI & MboI, ligation site: GATCGATC, RE site: GATC & GATC\n",
      "Preparing MAP file\n",
      "   x removing pre-GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN_filt_1-end.map\n",
      "Mapping fragments of remaining reads...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_lYQC9s\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_lYQC9s -o /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_lYQC9s_frag_1-end\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_lYQC9s\n",
      "   x removing failed to map /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_Q50vCN_fail.map\n",
      "   x removing tmp mapped /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1_tmp/mouse_B_rep1_1.fastq_lYQC9s_frag_1-end.map\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1/mouse_B_rep1_1.fastq_full_1-end.map',\n",
       " 'results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r1/mouse_B_rep1_1.fastq_frag_1-end.map']"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# for the first side of the reads \n",
    "full_mapping(gem_index_path='genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem',\n",
    "             out_map_dir='results/fragment/{0}_{1}/01_mapping/mapped_{0}_{1}_r1/'.format(cell, rep),\n",
    "             fastq_path='FASTQs/%s_%s_1.fastq.dsrc' % (cell, rep),\n",
    "             r_enz='MboI', frag_map=True, clean=True, nthreads=8, \n",
    "             temp_dir='results/fragment/{0}_{1}/01_mapping/mapped_{0}_{1}_r1_tmp/'.format(cell, rep))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Preparing FASTQ file\n",
      "  - conversion to MAP format\n",
      "Mapping reads in window 1-end...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj -o /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj_full_1-end\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj\n",
      "   x removing map /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj_full_1-end.map\n",
      "  - splitting into restriction enzyme (RE) fragments using ligation sites\n",
      "  - ligation sites are replaced by RE sites to match the reference genome\n",
      "    * enzymes: MboI & MboI, ligation site: GATCGATC, RE site: GATC & GATC\n",
      "Preparing MAP file\n",
      "   x removing pre-GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj_filt_1-end.map\n",
      "Mapping fragments of remaining reads...\n",
      "TO GEM /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_V5kpMT\n",
      "/home/dcastillo/miniconda2/bin/gem-mapper -I /scratch/workspace/MethodsMolBiol/Notebooks/genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem -q offset-33 -m 0.04 -s 0 --allow-incomplete-strata 0.00 --granularity 10000 --max-decoded-matches 1 --min-decoded-strata 0 --min-insert-size 0 --max-insert-size 0 --min-matched-bases 0.8 --gem-quality-threshold 26 --max-big-indel-length 15 --mismatch-alphabet ACGT -E 0.30 --max-extendable-matches 20 --max-extensions-per-match 1 -e 0.04 -T 8 -i /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_V5kpMT -o /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_V5kpMT_frag_1-end\n",
      "Parsing result...\n",
      "   x removing GEM input /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_V5kpMT\n",
      "   x removing failed to map /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_aAwRcj_fail.map\n",
      "   x removing tmp mapped /scratch/workspace/MethodsMolBiol/Notebooks/results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2_tmp/mouse_B_rep1_2.fastq_V5kpMT_frag_1-end.map\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "['results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_full_1-end.map',\n",
       " 'results/fragment/mouse_B_rep1/01_mapping/mapped_mouse_B_rep1_r2/mouse_B_rep1_2.fastq_frag_1-end.map']"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# for the second side of the reads\n",
    "full_mapping(gem_index_path='genome/Mus_musculus-GRCm38.p6/Mus_musculus-GRCm38.p6_contigs.gem',\n",
    "             out_map_dir='results/fragment/{0}_{1}/01_mapping/mapped_{0}_{1}_r2/'.format(cell, rep),\n",
    "             fastq_path='FASTQs/%s_%s_2.fastq.dsrc' % (cell, rep),\n",
    "             r_enz='MboI', frag_map=True, clean=True, nthreads=8, \n",
    "             temp_dir='results/fragment/{0}_{1}/01_mapping/mapped_{0}_{1}_r2_tmp/'.format(cell, rep))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "<!--bibtex\n",
    "@article{Imakaev2012a,\n",
    "abstract = {Extracting biologically meaningful information from chromosomal interactions obtained with genome-wide chromosome conformation capture (3C) analyses requires the elimination of systematic biases. We present a computational pipeline that integrates a strategy to map sequencing reads with a data-driven method for iterative correction of biases, yielding genome-wide maps of relative contact probabilities. We validate this ICE (iterative correction and eigenvector decomposition) technique on published data obtained by the high-throughput 3C method Hi-C, and we demonstrate that eigenvector decomposition of the obtained maps provides insights into local chromatin states, global patterns of chromosomal interactions, and the conserved organization of human and mouse chromosomes.},\n",
    "author = {Imakaev, Maxim V and Fudenberg, Geoffrey and McCord, Rachel Patton and Naumova, Natalia and Goloborodko, Anton and Lajoie, Bryan R and Dekker, Job and Mirny, Leonid A},\n",
    "doi = {10.1038/nmeth.2148},\n",
    "file = {:home/fransua/.local/share/data/Mendeley Ltd./Mendeley Desktop/Downloaded/Imakaev et al. - 2012 - Iterative correction of Hi-C data reveals hallmarks of chromosome organization.pdf:pdf},\n",
    "issn = {1548-7105},\n",
    "journal = {Nature methods},\n",
    "keywords = {Hi-C},\n",
    "mendeley-groups = {stats/Hi-C,Research articles},\n",
    "mendeley-tags = {Hi-C},\n",
    "month = {oct},\n",
    "number = {10},\n",
    "pages = {999--1003},\n",
    "pmid = {22941365},\n",
    "title = {{Iterative correction of Hi-C data reveals hallmarks of chromosome organization.}},\n",
    "url = {http://www.ncbi.nlm.nih.gov/pubmed/22941365},\n",
    "volume = {9},\n",
    "year = {2012}\n",
    "}\n",
    "@article{Ay2015a,\n",
    "author = {Ay, Ferhat and Noble, William Stafford},\n",
    "doi = {10.1186/s13059-015-0745-7},\n",
    "file = {:home/fransua/.local/share/data/Mendeley Ltd./Mendeley Desktop/Downloaded/Ay, Noble - 2015 - Analysis methods for studying the 3D architecture of the genome.pdf:pdf},\n",
    "issn = {1474-760X},\n",
    "journal = {Genome Biology},\n",
    "keywords = {Chromatin conformation capture,Genome architecture,Hi-C,Three-dimensional genome,Three-dimensional modeling,chromatin,conformation capture,genome architecture,in many other applications,ranging from genome assem-,review,three-dimensional genome,three-dimensional modeling},\n",
    "mendeley-groups = {Research articles},\n",
    "mendeley-tags = {Hi-C,review},\n",
    "number = {1},\n",
    "pages = {183},\n",
    "publisher = {Genome Biology},\n",
    "title = {{Analysis methods for studying the 3D architecture of the genome}},\n",
    "url = {http://genomebiology.com/2015/16/1/183},\n",
    "volume = {16},\n",
    "year = {2015}\n",
    "}\n",
    "@article{Wingett2015,\n",
    "abstract = {HiCUP is a pipeline for processing sequence data generated by Hi-C and Capture Hi-C (CHi-C) experiments, which are techniques used to investigate three-dimensional genomic organisation. The pipeline maps data to a specified reference genome and removes artefacts that would otherwise hinder subsequent analysis. HiCUP also produces an easy-to-interpret yet detailed quality control (QC) report that assists in refining experimental protocols for future studies. The software is freely available and has already been used for processing Hi-C and CHi-C data in several recently published peer-reviewed studies.},\n",
    "author = {Wingett, Steven and Ewels, Philip and Furlan-Magaril, Mayra and Nagano, Takashi and Schoenfelder, Stefan and Fraser, Peter and Andrews, Simon},\n",
    "doi = {10.12688/f1000research.7334.1},\n",
    "file = {:home/fransua/Downloads/f1000research-4-7903.pdf:pdf},\n",
    "issn = {2046-1402},\n",
    "journal = {F1000Research},\n",
    "mendeley-groups = {Computer programs/Hi-C/Hi-C processing},\n",
    "pages = {1310},\n",
    "pmid = {26835000},\n",
    "title = {{HiCUP: pipeline for mapping and processing Hi-C data.}},\n",
    "url = {http://f1000research.com/articles/4-1310/v1},\n",
    "volume = {4},\n",
    "year = {2015}\n",
    "}\n",
    "@article{Serra2016,\n",
    "abstract = {The sequence of a genome is insufficient to understand all genomic processes carried out in the cell nucleus. To achieve this, the knowledge of its three- dimensional architecture is necessary. Advances in genomic technologies and the development of new analytical methods, such as Chromosome Conformation Capture (3C) and its derivatives, now permit to investigate the spatial organization of genomes. However, inferring structures from raw contact data is a tedious process for shortage of available tools. Here we present TADbit, a computational framework to analyze and model the chromatin fiber in three dimensions. To illustrate the use of TADbit, we automatically modeled 50 genomic domains from the fly genome revealing differential structural features of the previously defined chromatin colors, establishing a link between the conformation of the genome and the local chromatin composition. More generally, TADbit allows to obtain three-dimensional models ready for visualization from 3C-based experiments and to characterize their relation to gene expression and epigenetic states. TADbit is open-source and available for download from http://www.3DGenomes.org.},\n",
    "author = {Serra, Fran{\\c{c}}ois and Ba{\\`{u}}, Davide and Filion, Guillaume and Marti-Renom, Marc A.},\n",
    "doi = {10.1101/036764},\n",
    "file = {:home/fransua/.local/share/data/Mendeley Ltd./Mendeley Desktop/Downloaded/Serra et al. - 2016 - Structural features of the fly chromatin colors revealed by automatic three-dimensional modeling.pdf:pdf},\n",
    "journal = {bioRxiv},\n",
    "keywords = {3d,Hi-C,capture,genome architecture,genome reconstruction chromosome conformation,modeling,optimization,resampling,restraint-based,three-dimensional genome reconstruction},\n",
    "mendeley-groups = {Research articles,projects/WIREs{\\_}review/modeling{\\_}tools,Computer programs/Hi-C/Hi-C modeling},\n",
    "mendeley-tags = {Hi-C,modeling,optimization,resampling},\n",
    "pages = {1--29},\n",
    "title = {{Structural features of the fly chromatin colors revealed by automatic three-dimensional modeling.}},\n",
    "url = {http://biorxiv.org/content/early/2016/01/15/036764},\n",
    "year = {2016}\n",
    "}\n",
    "@misc{Marco-Sola2012,\n",
    "abstract = {Because of ever-increasing throughput requirements of sequencing data, most existing short-read aligners have been designed to focus on speed at the expense of accuracy. The Genome Multitool (GEM) mapper can leverage string matching by filtration to search the alignment space more efficiently, simultaneously delivering precision (performing fully tunable exhaustive searches that return all existing matches, including gapped ones) and speed (being several times faster than comparable state-of-the-art tools).},\n",
    "author = {Marco-Sola, Santiago and Sammeth, Michael and Guig{\\'{o}}, Roderic and Ribeca, Paolo},\n",
    "booktitle = {Nature Methods},\n",
    "doi = {10.1038/nmeth.2221},\n",
    "isbn = {1548-7105 (Electronic)$\\backslash$r1548-7091 (Linking)},\n",
    "issn = {1548-7091},\n",
    "mendeley-groups = {Research articles},\n",
    "pmid = {23103880},\n",
    "title = {{The GEM mapper: fast, accurate and versatile alignment by filtration}},\n",
    "year = {2012}\n",
    "}\n",
    "\n",
    "-->\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# References\n",
    "\n",
    "<a name=\"cite-Imakaev2012a\"/><sup>[^](#ref-1) </sup>Imakaev, Maxim V and Fudenberg, Geoffrey and McCord, Rachel Patton and Naumova, Natalia and Goloborodko, Anton and Lajoie, Bryan R and Dekker, Job and Mirny, Leonid A. 2012. _Iterative correction of Hi-C data reveals hallmarks of chromosome organization._. [URL](http://www.ncbi.nlm.nih.gov/pubmed/22941365)\n",
    "\n",
    "<a name=\"cite-Ay2015a\"/><sup>[^](#ref-2) </sup>Ay, Ferhat and Noble, William Stafford. 2015. _Analysis methods for studying the 3D architecture of the genome_. [URL](http://genomebiology.com/2015/16/1/183)\n",
    "\n",
    "<a name=\"cite-Wingett2015\"/><sup>[^](#ref-3) </sup>Wingett, Steven and Ewels, Philip and Furlan-Magaril, Mayra and Nagano, Takashi and Schoenfelder, Stefan and Fraser, Peter and Andrews, Simon. 2015. _HiCUP: pipeline for mapping and processing Hi-C data._. [URL](http://f1000research.com/articles/4-1310/v1)\n",
    "\n",
    "<a name=\"cite-Serra2017\"/><sup>[^](#ref-4) </sup>François Serra, Davide Baù, Mike Goodstadt, David Castillo,  Guillaume Filion and Marc A. Marti-Renom. 2017. _Automatic analysis and 3D-modelling of Hi-C data using TADbit reveals structural features of the fly chromatin colors_. [URL](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005665)\n",
    "\n",
    "<a name=\"cite-Marco-Sola2012\"/><sup>[^](#ref-5) </sup>Marco-Sola, Santiago and Sammeth, Michael and Guigó, Roderic and Ribeca, Paolo. 2012. _The GEM mapper: fast, accurate and versatile alignment by filtration_.\n",
    "\n"
   ]
  }
 ],
 "metadata": {
  "hide_input": false,
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.15"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}