Notebooks/A1-Preparation reference genome.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "toc": true
   },
   "source": [
    "<h1>Table of Contents<span class=\"tocSkip\"></span></h1>\n",
    "<div class=\"toc\" style=\"margin-top: 1em;\"><ul class=\"toc-item\"><li><span><a href=\"#Search-for-a-reference-genome\" data-toc-modified-id=\"Search-for-a-reference-genome-1\">Search for a reference genome</a></span><ul class=\"toc-item\"><li><span><a href=\"#Homo-sapiens's-reference-genome-sequence\" data-toc-modified-id=\"Homo-sapiens's-reference-genome-sequence-1.1\"><em>Homo sapiens</em>'s reference genome sequence</a></span></li><li><span><a href=\"#Mus-musculus's-reference-genome-sequence\" data-toc-modified-id=\"Mus-musculus's-reference-genome-sequence-1.2\"><em>Mus musculus</em>'s reference genome sequence</a></span></li><li><span><a href=\"#Saccharomyces-cerevisiae's-reference-genome-sequence\" data-toc-modified-id=\"Saccharomyces-cerevisiae's-reference-genome-sequence-1.3\"><em>Saccharomyces cerevisiae</em>'s reference genome sequence</a></span></li></ul></li><li><span><a href=\"#Download-from-the-NCBI\" data-toc-modified-id=\"Download-from-the-NCBI-2\">Download from the NCBI</a></span><ul class=\"toc-item\"><li><span><a href=\"#List-of-chromosome/contigs\" data-toc-modified-id=\"List-of-chromosome/contigs-2.1\">List of chromosome/contigs</a></span></li><li><span><a href=\"#Sequences-of-each-chromosome/contigs\" data-toc-modified-id=\"Sequences-of-each-chromosome/contigs-2.2\">Sequences of each chromosome/contigs</a></span></li><li><span><a href=\"#Concatenate-all-contigs/chromosomes-into-single-files\" data-toc-modified-id=\"Concatenate-all-contigs/chromosomes-into-single-files-2.3\">Concatenate all contigs/chromosomes into single files</a></span></li></ul></li><li><span><a href=\"#Creation-of-an-index-file-for-GEM-mapper\" data-toc-modified-id=\"Creation-of-an-index-file-for-GEM-mapper-3\">Creation of an index file for GEM mapper</a></span><ul class=\"toc-item\"><li><span><a href=\"#Extra:-Compute-mappability-for-oneD-normalization\" data-toc-modified-id=\"Extra:-Compute-mappability-for-oneD-normalization-3.1\">Extra: Compute mappability for oneD normalization</a></span></li></ul></li></ul></div>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Search for a reference genome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We would need two reference genomes. One as a FASTA file with each chromosome, and one that we will use exclusively for the mapping that would contain all contigs.\n",
    "\n",
    "*The use of contigs in the reference genome increases the mapping specificity.*"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *Homo sapiens*'s reference genome sequence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Human (https://www.ncbi.nlm.nih.gov/genome/?term=Homo+sapiens%5Borgn%5D):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "format": "row"
   },
   "outputs": [],
   "source": [
    "species  = 'Homo_sapiens'\n",
    "taxid    = '9606'\n",
    "assembly = 'GRCh38.p12'\n",
    "genbank  = 'GCF_000001405.38'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *Mus musculus*'s reference genome sequence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Mouse (https://www.ncbi.nlm.nih.gov/genome?term=mus%20musculus):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "format": "row"
   },
   "outputs": [],
   "source": [
    "species  = 'Mus_musculus'\n",
    "taxid    = '10090'\n",
    "assembly = 'GRCm38.p6'\n",
    "genbank  = 'GCF_000001635.26'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *Saccharomyces cerevisiae*'s reference genome sequence"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Yeast (https://www.ncbi.nlm.nih.gov/genome/?term=Saccharomyces+cerevisiae):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "format": "row"
   },
   "outputs": [],
   "source": [
    "species = 'Saccharomyces_cerevisiae'\n",
    "taxid   = '4932'\n",
    "assembly = 'R64'\n",
    "genbank = 'GCA_000146045.2'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Download from the NCBI"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  List of chromosome/contigs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "sumurl = ('ftp://ftp.ncbi.nlm.nih.gov/genomes/all/{0}/{1}/{2}/{3}/{4}_{5}/'\n",
    "          '{4}_{5}_assembly_report.txt').format(genbank[:3], genbank[4:7], genbank[7:10], \n",
    "                                                genbank[10:13], genbank, assembly)\n",
    "\n",
    "crmurl = ('https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi'\n",
    "          '?db=nuccore&id=%s&rettype=fasta&retmode=text')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCA/000/146/045/GCA_000146045.2_R64/GCA_000146045.2_R64_assembly_report.txt\n"
     ]
    }
   ],
   "source": [
    "print sumurl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "! wget -q $sumurl -O chromosome_list.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Assembly name:  R64\r",
      "\r\n",
      "# Organism name:  Saccharomyces cerevisiae S288C (baker's yeast)\r",
      "\r\n",
      "# Infraspecific name:  strain=S288C\r",
      "\r\n",
      "# Taxid:          559292\r",
      "\r\n",
      "# BioProject:     PRJNA43747\r",
      "\r\n",
      "# Submitter:      Saccharomyces Genome Database\r",
      "\r\n",
      "# Date:           2011-4-18\r",
      "\r\n",
      "# Synonyms:       sacCer3\t\r",
      "\r\n",
      "# Assembly type:  haploid\r",
      "\r\n",
      "# Release type:   major\r",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! head chromosome_list.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Sequences of each chromosome/contigs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "format": "row"
   },
   "outputs": [],
   "source": [
    "dirname = 'genome'\n",
    "! mkdir -p {dirname}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For each contig/chromosome download the corresponding FASTA file from NCBI"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "format": "row",
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "contig = []\n",
    "for line in open('chromosome_list.txt'):\n",
    "    if line.startswith('#'):\n",
    "        continue\n",
    "    seq_name, seq_role, assigned_molecule, _, genbank, _, refseq, _ = line.split(None, 7)\n",
    "    if seq_role == 'assembled-molecule':\n",
    "        name = 'chr%s.fasta' % assigned_molecule\n",
    "    else:\n",
    "        name = 'chr%s_%s.fasta' % (assigned_molecule, seq_name.replace('/', '-'))\n",
    "    contig.append(name)\n",
    "\n",
    "    outfile = os.path.join(dirname, name)\n",
    "    if os.path.exists(outfile) and os.path.getsize(outfile) > 10:\n",
    "        continue\n",
    "    error_code = os.system('wget \"%s\" --no-check-certificate -O %s' % (crmurl % (genbank), outfile))\n",
    "    if error_code:\n",
    "        error_code = os.system('wget \"%s\" --no-check-certificate -O %s' % (crmurl % (refseq), outfile))\n",
    "    if error_code:\n",
    "        print genbank"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Concatenate all contigs/chromosomes into single files"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "def write_to_fasta(line):\n",
    "    contig_file.write(line)\n",
    "\n",
    "def write_to_fastas(line):\n",
    "    contig_file.write(line)\n",
    "    simple_file.write(line)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "os.system('mkdir -p {}/{}-{}'.format(dirname, species, assembly))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "contig_file = open('{0}/{1}-{2}/{1}-{2}_contigs.fa'.format(dirname, species, assembly),'w')\n",
    "simple_file = open('{0}/{1}-{2}/{1}-{2}.fa'.format(dirname, species, assembly),'w')\n",
    "\n",
    "for molecule in contig:\n",
    "    fh = open('{0}/{1}'.format(dirname, molecule))\n",
    "    oline = '>%s\\n' % (molecule.replace('.fasta', ''))\n",
    "    _ = fh.next()\n",
    "    # if molecule is an assembled chromosome we write to both files, otherwise only to the *_contigs one\n",
    "    write = write_to_fasta if '_' in molecule else write_to_fastas\n",
    "    for line in fh:\n",
    "        write(oline)\n",
    "        oline = line\n",
    "    # last line usually empty...\n",
    "    if line.strip():\n",
    "        write(line)\n",
    "contig_file.close()\n",
    "simple_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove all the other files (with single chromosome/contig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "! rm -f {dirname}/*.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Creation of an index file for GEM mapper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to GEM-indexer build 1.335 (beta) - (2012-10-29 22:24:35 GMT)\n",
      " (c) 2008-2012 Paolo Ribeca <paolo.ribeca@gmail.com>\n",
      " (c) 2010-2012 Santiago Marco Sola <santiagomsola@gmail.com>\n",
      " (c) 2010-2012 Leonor Frias Moya <leonor.frias@gmail.com>\n",
      "For the terms of use, run the program with the option --show-license.\n",
      "************************************************************************\n",
      "* WARNING: this is a beta version, provided for testing purposes only; *\n",
      "*          check for updates at <http://gemlibrary.sourceforge.net>.   *\n",
      "************************************************************************\n",
      "Creating sequence and location files... done.\n",
      "Computing DNA BWT (likely to take long)... done.\n",
      "Generating index (likely to take long)... done.\n",
      "Cleaning up... done.\n"
     ]
    }
   ],
   "source": [
    "! gem-indexer -t 8 -i {dirname}/{species}-{assembly}/{species}-{assembly}_contigs.fa -o {dirname}/{species}-{assembly}/{species}_contigs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index file will be: __`{dirname}/{species}-{assembly}/{species}_contigs.gem`__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extra: Compute mappability for oneD normalization"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this case we can use the FASTA of the genome whithout contigs and follow these step:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to GEM-indexer build 1.335 (beta) - (2012-10-29 22:24:35 GMT)\n",
      " (c) 2008-2012 Paolo Ribeca <paolo.ribeca@gmail.com>\n",
      " (c) 2010-2012 Santiago Marco Sola <santiagomsola@gmail.com>\n",
      " (c) 2010-2012 Leonor Frias Moya <leonor.frias@gmail.com>\n",
      "For the terms of use, run the program with the option --show-license.\n",
      "************************************************************************\n",
      "* WARNING: this is a beta version, provided for testing purposes only; *\n",
      "*          check for updates at <http://gemlibrary.sourceforge.net>.   *\n",
      "************************************************************************\n",
      "Creating sequence and location files... done.\n",
      "Computing DNA BWT (likely to take long)... done.\n",
      "Generating index (likely to take long)... done.\n",
      "Cleaning up... done.\n",
      "Welcome to GEM-mappability build 1.315 (beta) - (2013/03/29 02:59:40 GMT)\n",
      " (c) 2008-2013 Paolo Ribeca <paolo.ribeca@gmail.com>\n",
      " (c) 2010-2013 Santiago Marco Sola <santiagomsola@gmail.com>\n",
      " (c) 2010-2013 Leonor Frias Moya <leonor.frias@gmail.com>\n",
      "For the terms of use, run the program with the option --show-license.\n",
      "************************************************************************\n",
      "* WARNING: this is a beta version, provided for testing purposes only; *\n",
      "*          check for updates at <http://gemlibrary.sourceforge.net>.   *\n",
      "************************************************************************\n",
      "Fri Sep 14 12:18:47 2018 -- Loading index (likely to take long)... done.\n",
      "Fri Sep 14 12:18:48 2018 -- Starting (24314194 positions to go)...\n",
      "\u001b[KFri Sep 14 12:19:09 2018 -- Pos=10.6% Done=3273532(13.5%) Uniq=72.2%%%ep 14 12:18:49 2018 -- Pos=0.494% Done=107182(0.441%) Uniq=91.7% Sep 14 12:18:49 2018 -- Pos=0.123% Done=138169(0.568%) Uniq=84.5%Sep 14 12:18:51 2018 -- Pos=0.987% Done=399459(1.64%) Uniq=45.2%Fri Sep 14 12:18:51 2018 -- Pos=0.905% Done=419459(1.73%) Uniq=45.5%Fri Sep 14 12:18:51 2018 -- Pos=1.27% Done=475921(1.96%) Uniq=46.1%Fri Sep 14 12:18:53 2018 -- Pos=1.85% Done=877271(3.61%) Uniq=49.7% Sep 14 12:18:54 2018 -- Pos=2.67% Done=1047658(4.31%) Uniq=51.5%Fri Sep 14 12:18:55 2018 -- Pos=2.71% Done=1079382(4.44%) Uniq=52.7% Sep 14 12:18:55 2018 -- Pos=2.92% Done=1196101(4.92%) Uniq=49.4%Fri Sep 14 12:18:55 2018 -- Pos=3.17% Done=1276448(5.25%) Uniq=50.4%Fri Sep 14 12:18:56 2018 -- Pos=3.25% Done=1305979(5.37%) Uniq=51.1%Fri Sep 14 12:18:56 2018 -- Pos=3.37% Done=1346192(5.54%) Uniq=52.5%Fri Sep 14 12:18:57 2018 -- Pos=4.03% Done=1487152(6.12%) Uniq=56.8% Sep 14 12:18:57 2018 -- Pos=4.36% Done=1587153(6.53%) Uniq=59.5%Fri Sep 14 12:18:58 2018 -- Pos=4.65% Done=1627153(6.69%) Uniq=60.4% 14 12:18:58 2018 -- Pos=4.81% Done=1687368(6.94%) Uniq=61.8%Fri Sep 14 12:18:59 2018 -- Pos=5.43% Done=1797368(7.39%) Uniq=64.1% Sep 14 12:18:59 2018 -- Pos=5.31% Done=1807368(7.43%) Uniq=64.2%Fri Sep 14 12:18:59 2018 -- Pos=5.8% Done=1897438(7.8%) Uniq=65.7% Sep 14 12:18:59 2018 -- Pos=5.76% Done=1917438(7.89%) Uniq=66% 14 12:19:02 2018 -- Pos=6.91% Done=2211118(9.09%) Uniq=69.3%Fri Sep 14 12:19:03 2018 -- Pos=7.57% Done=2362253(9.72%) Uniq=71.2%Fri Sep 14 12:19:03 2018 -- Pos=7.94% Done=2474182(10.2%) Uniq=71.3%Fri Sep 14 12:19:04 2018 -- Pos=8.06% Done=2513793(10.3%) Uniq=71.5%Fri Sep 14 12:19:04 2018 -- Pos=8.1% Done=2523793(10.4%) Uniq=71.7%Sep 14 12:19:06 2018 -- Pos=9.34% Done=2867903(11.8%) Uniq=72.3%Fri Sep 14 12:19:07 2018 -- Pos=9.75% Done=2957903(12.2%) Uniq=73%Fri Sep 14 12:19:07 2018 -- Pos=9.71% Done=3072329(12.6%) Uniq=71%Fri Sep 14 12:19:08 2018 -- Pos=9.91% Done=3182823(13.1%) Uniq=71.6%Fri Sep 14 12:19:08 2018 -- Pos=10.4% Done=3232823(13.3%) Uniq=72%"
     ]
    },
    {
     "data": {
      "text/html": [
       "<b>limit_output extension: Maximum message size of 20000 exceeded with 20038 characters</b>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "! gem-indexer -i {dirname}/{species}-{assembly}/{species}-{assembly}.fa \\\n",
    "   -o {dirname}/{species}-{assembly}/{species}-{assembly} -t 8\n",
    "\n",
    "! gem-mappability -I {dirname}/{species}-{assembly}/{species}-{assembly}.gem -l 50 \\\n",
    "   -o {dirname}/{species}-{assembly}/{species}-{assembly}.50mer -T 8\n",
    "\n",
    "! gem-2-wig -I {dirname}/{species}-{assembly}/{species}-{assembly}.gem \\\n",
    "   -i {dirname}/{species}-{assembly}/{species}-{assembly}.50mer.mappability \\\n",
    "   -o {dirname}/{species}-{assembly}/{species}-{assembly}.50mer\n",
    "\n",
    "! wigToBigWig {dirname}/{species}-{assembly}/{species}-{assembly}.50mer.wig \\\n",
    "   {dirname}/{species}-{assembly}/{species}-{assembly}.50mer.sizes Saccharomyces_cerevisiae.50mer.bw\n",
    "\n",
    "! bigWigToBedGraph {dirname}/{species}-{assembly}/{species}-{assembly}.50mer.bw  \\\n",
    "   {dirname}/{species}-{assembly}/{species}-{assembly}.50mer.bedGraph"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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": {
   "nav_menu": {},
   "number_sections": false,
   "sideBar": true,
   "skip_h1_title": false,
   "toc_cell": true,
   "toc_position": {},
   "toc_section_display": "block",
   "toc_window_display": true
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}