Notebooks/A1-Preparation_reference_genome.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Search for a reference genome"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## *Homo sapiens*'s reference genome sequence"
   ]
  },
  {
   "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": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "format": "row"
   },
   "outputs": [],
   "source": [
    "species = 'Homo sapiens'\n",
    "taxid   = '9606'\n",
    "genome  = 'GRCh38.p10'\n",
    "refseq, dir1, dir2, dir3 = 'GCF', '000', '001', '405' \n",
    "genbank = 'GCF_000001405.36'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "sumurl = 'ftp://ftp.ncbi.nlm.nih.gov/genomes/all/{0}/{1}/{2}/{3}/{4}_{5}/{4}_{5}_assembly_report.txt'.format(\n",
    "    refseq, dir1, dir2, dir3, genbank, genome)\n",
    "\n",
    "crmurl = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=%s&rettype=fasta&retmode=text'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Download from the NCBI the list of chromosome/contigs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "! wget -q $sumurl -O chromosome_list.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": false,
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "# Assembly name:  GRCh38.p10\r",
      "\r\n",
      "# Description:    Genome Reference Consortium Human Build 38 patch release 10 (GRCh38.p10)\r",
      "\r\n",
      "# Organism name:  Homo sapiens (human)\r",
      "\r\n",
      "# Taxid:          9606\r",
      "\r\n",
      "# BioProject:     PRJNA168\r",
      "\r\n",
      "# Submitter:      Genome Reference Consortium\r",
      "\r\n",
      "# Date:           2017-1-6\r",
      "\r\n",
      "# Assembly type:  haploid-with-alt-loci\r",
      "\r\n",
      "# Release type:   patch\r",
      "\r\n",
      "# Assembly level: Chromosome\r",
      "\r\n"
     ]
    }
   ],
   "source": [
    "! head chromosome_list.txt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "collapsed": true,
    "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": null,
   "metadata": {
    "collapsed": false,
    "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 a single file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "contig_file = open('genome/Homo_sapiens_contigs.fa','w')\n",
    "for molecule in contig:\n",
    "    for line in open('genome/' + molecule):\n",
    "        # replace the header of the sequence in the fasta file\n",
    "        if line == '\\n':\n",
    "            continue\n",
    "        if line.startswith('>'):\n",
    "            line = '>' + molecule[3:].replace('.fasta', '') + '\\n'\n",
    "        contig_file.write(line)\n",
    "contig_file.close()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Remove all the other files (with single chromosome/contig)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "! rm -f genome/*.fasta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Creation of an index file for GEM mapper"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Welcome to GEM-indexer build 1.423 (beta) - (2013/04/01 01:02:13 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",
      "Creating sequence and location files... done.\n",
      "Computing DNA BWT (likely to take long)..."
     ]
    }
   ],
   "source": [
    "! gem-indexer -t 8 -i genome/Homo_sapiens_contigs.fa -o genome/Homo_sapiens_contigs"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The index file will be: __`genome/Homo_sapiens_contigs.gem`__"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*__WARNING__: in more recent versions of GEM the \"-t\" flag should be \"-T\"*"
   ]
  }
 ],
 "metadata": {
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}