Notebooks/Methods-3-Hi-C quality check.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=\"#FASTQ-format\" data-toc-modified-id=\"FASTQ-format-1\">FASTQ format</a></span></li><li><span><a href=\"#Quality-check-before-mapping\" data-toc-modified-id=\"Quality-check-before-mapping-2\">Quality check before mapping</a></span><ul class=\"toc-item\"><li><span><a href=\"#Check-restriction-enzyme-used\" data-toc-modified-id=\"Check-restriction-enzyme-used-2.1\">Check restriction-enzyme used</a></span></li><li><span><a href=\"#Plot-PHRED-score-and-ligation/digestion-sites\" data-toc-modified-id=\"Plot-PHRED-score-and-ligation/digestion-sites-2.2\">Plot PHRED score and ligation/digestion sites</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FASTQ format"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The file is organized in 4 lines per read:\n",
    " 1. starting with `@`, the header of the DNA sequence with the read id (plus optional fields)\n",
    " 2. the DNA sequence\n",
    " 3. starting with `+`, the header of the sequence quality (this line could be either a repetition of first line or empty)\n",
    " 4. the sequence quality (it is provided as PHRED score and it is not human readable. Check   https://en.wikipedia.org/wiki/Phred_quality_score for more details)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "@SRR5344921.1\n",
      "ACCACTGATGGCACCACTGTTGTACTTTAGGATCGATCACTCCTTATACGCACATAGGAAACACTGACTAGAGGG\n",
      "+\n",
      "@B@B@>1>=;1@EG0B0;1=E1:>>1111<?=1E9=0><C:<>11=11?1=ECDC1>1=@FD1<:1?C#######\n",
      "@SRR5344921.2\n",
      "GGGCTTTTTTTGGCTGGGAGACTATTTATAACTGCTTCTATTTCTTTAGGGGATATGGGACTGTTTAGAAGGTCA\n",
      "+\n",
      "BCCBBFGGGG@GEGGGGDGEBGEF@<1FGGEFGGECGG@FEGGGEG>@1=09/BFGGEG111C<F1>11BGF1=1\n"
     ]
    }
   ],
   "source": [
    "%%bash \n",
    "\n",
    "dsrc d -s FASTQs/mouse_B_rep1_1.fastq.dsrc | head -n 8"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Count the number of lines in the file (4 times the number of reads)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "400000000\n"
     ]
    }
   ],
   "source": [
    "%%bash\n",
    "\n",
    "dsrc d -s FASTQs/mouse_B_rep1_1.fastq.dsrc | wc -l"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are 400 M lines in the file, which means 100M reads in total."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Quality check before mapping"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Check restriction-enzyme used"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most important to analyze Hi-C dataset is the restriction enzyme used in the experiment. TADbit provides a simple function to check for it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pytadbit.mapping.restriction_enzymes import identify_re"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "- Most probable pattern: GATC, matching enzymes: BssMI,DpnII,Sau3AI,MboI,Kzo9I,Bsp143I,NdeII,BstMBI,BfuCI\n"
     ]
    }
   ],
   "source": [
    "pat, enz, pv = identify_re('FASTQs/mouse_B_rep1_1.fastq.dsrc')\n",
    "\n",
    "print '- Most probable pattern: %s, matching enzymes: %s' % (pat, ','.join(enz))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Plot PHRED score and ligation/digestion sites"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to quickly assess the quality of the HiC experiment (before mapping), and given that we know the restriction enzyme used, we can check the proportion of reads with ligation sites as well as the number of reads starting by a cut-site.\n",
    "\n",
    "These numbers will give us a first hint on the __efficiencies of two critical steps in the HiC experiment, the digestion and the ligation.__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from pytadbit.utils.fastq_utils import quality_plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "r_enz = 'MboI'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "cell = 'B'\n",
    "repl = 'rep1'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The plot on the top represents the typical per nucleotide quality profile of NGS reads, with, in addition, the proportion of `N` found at each position. \n",
    "\n",
    "The second plot, is specific to Hi-C experiments. Given a restriction enzyme the function searches for the presence of ligation sites and of undigested restriction enzyme sites. Depending on the enzyme used the function can differentiate between dangling-ends and undigested sites.\n",
    "\n",
    "From these proportions some quality statistics can be inferred before mapping:\n",
    "- The __PHRED score__ and the __number of unidentified nucleotides__ (Ns) in the read sequence, which are routinely computed to address the quality of high-throughput sequencing experiments.\n",
    "- The numbers of __undigested and unligated RE sites per-nucleotide__ along the read to assess the quality of the Hi-C experiment.\n",
    "- The overall __percentage of digested sites__, which relates directly to the RE efficiency.\n",
    "- The __percentage of non-ligated digested__ (dangling-ends), which relates to the ligation efficiency.\n",
    "- The __percentage of read-ends with a ligation site__, which is negatively correlated with the percentage of dangling-ends."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "quality_plot('FASTQs/mouse_{0}_{1}_1.fastq.dsrc'.format(cell, repl), r_enz=r_enz, nreads=1000000)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "_Note: this plot is __compatible with the use of multiple restriction enzymes__ (which is why the ligation site is labeled as MboI-MboI)_"
   ]
  }
 ],
 "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": {
    "height": "889.5px",
    "left": "0px",
    "right": "1353px",
    "top": "110.5px",
    "width": "202px"
   },
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}