{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to download background genes from NCBI\n",
    "\n",
    "## Example\n",
    "### 1) Download mouse (TaxID=10090) protein-coding genes\n",
    "1. **Query [NCBI Gene](https://www.ncbi.nlm.nih.gov/gene):**    \n",
    "   `\"10090\"[Taxonomy ID] AND alive[property] AND genetype protein coding[Properties]`\n",
    "2. **Click \"Send to:\"**\n",
    "3. **Select \"File\"**\n",
    "4. **Select \"Create File\" button**\n",
    "   The default name of the tsv file is `gene_result.txt`\n",
    "\n",
    "Note: To download all mouse DNA items:    \n",
    "`\"10090\"[Taxonomy ID] AND alive[property]`\n",
    "\n",
    "![](images/dnld_mouse_pcd_genes.png)\n",
    "\n",
    "## 2) Convert NCBI Gene tab separated values (tsv) file to a Python module\n",
    "Use the command line or a Python script to convert a NCBI Gene tsv file to a Python module\n",
    "\n",
    "### 2a) Run a script from the command line\n",
    "```\n",
    "$ goatools ncbi_gene_results_to_python gene_result.txt -o genes_ncbi_10090_proteincoding.py\n",
    "      26,386 lines READ:  gene_result.txt\n",
    "      26,376 geneids WROTE: genes_ncbi_10090_proteincoding.py\n",
    "```\n",
    "\n",
    "### 2b) Run a function from inside your Python script"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "      26,386 lines READ:  gene_result.txt\n",
      "      26,376 geneids WROTE: genes_ncbi_10090_proteincoding.py\n"
     ]
    }
   ],
   "source": [
    "from goatools.cli.ncbi_gene_results_to_python import ncbi_tsv_to_py\n",
    "\n",
    "ncbi_tsv = 'gene_result.txt'\n",
    "output_py = 'genes_ncbi_10090_proteincoding.py'\n",
    "ncbi_tsv_to_py(ncbi_tsv, output_py)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Explore NCBI gene data\n",
    "### 3a) Import NCBI data from new NCBI gene Python module"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "from genes_ncbi_10090_proteincoding import GENEID2NT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3b) Examine fields stored in a namedtuple for a gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aliases         ['A1m', 'A2m', 'AI893533', 'MAM']\n",
      "CurrentID       0\n",
      "GeneID          11287\n",
      "OMIM            []\n",
      "Org_name        Mus musculus\n",
      "Status          live\n",
      "Symbol          Pzp\n",
      "chromosome      6\n",
      "description     PZP, alpha-2-macroglobulin like\n",
      "end_position_on_the_genomic_accession 128503683\n",
      "exon_count      36\n",
      "genomic_nucleotide_accession_version NC_000072.7\n",
      "map_location    6 63.02 cM\n",
      "no_hdr0         \n",
      "orientation     minus\n",
      "other_designations pregnancy zone protein|alpha 1 macroglobulin|alpha-2-M|alpha-2-macroglobulin\n",
      "start_position_on_the_genomic_accession 128460530\n",
      "tax_id          10090\n"
     ]
    }
   ],
   "source": [
    "# Get the data for one gene\n",
    "nt_gene = next(iter(sorted(GENEID2NT.values())))\n",
    "\n",
    "# Print the field name and value for all fields for one gene\n",
    "for key, val in sorted(nt_gene._asdict().items()):\n",
    "    print('{:15} {}'.format(key, val))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3c) Get genes which have specific genomic locations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22,216 genes have specific genomic basepair locations\n"
     ]
    }
   ],
   "source": [
    "nts = [nt for nt in GENEID2NT.values() if nt.start_position_on_the_genomic_accession != '']\n",
    "nts = sorted(nts, key=lambda nt: nt.GeneID)\n",
    "print('{N:,} genes have specific genomic basepair locations'.format(N=len(nts)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3d) Print GeneID, Symbol, and description of some genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GeneID Symbol   Description\n",
      "------ -------  --------------------------------------------------------\n",
      " 11287 Pzp      PZP, alpha-2-macroglobulin like\n",
      " 11298 Aanat    arylalkylamine N-acetyltransferase\n",
      " 11302 Aatk     apoptosis-associated tyrosine kinase\n",
      " 11303 Abca1    ATP-binding cassette, sub-family A (ABC1), member 1\n",
      " 11304 Abca4    ATP-binding cassette, sub-family A (ABC1), member 4\n",
      " 11305 Abca2    ATP-binding cassette, sub-family A (ABC1), member 2\n",
      " 11306 Abcb7    ATP-binding cassette, sub-family B (MDR/TAP), member 7\n",
      " 11307 Abcg1    ATP binding cassette subfamily G member 1\n",
      " 11308 Abi1     abl interactor 1\n",
      " 11350 Abl1     c-abl oncogene 1, non-receptor tyrosine kinase\n",
      " 11352 Abl2     v-abl Abelson murine leukemia viral oncogene 2 (arg, Abelson-related gene)\n",
      " 11354 Scgb1b27 secretoglobin, family 1B, member 27\n",
      " 11363 Acadl    acyl-Coenzyme A dehydrogenase, long-chain\n",
      " 11364 Acadm    acyl-Coenzyme A dehydrogenase, medium chain\n",
      " 11370 Acadvl   acyl-Coenzyme A dehydrogenase, very long chain\n",
      " 11409 Acads    acyl-Coenzyme A dehydrogenase, short chain\n",
      " 11416 Slc33a1  solute carrier family 33 (acetyl-CoA transporter), member 1\n",
      " 11418 Asic2    acid-sensing (proton-gated) ion channel 2\n",
      " 11419 Asic1    acid-sensing (proton-gated) ion channel 1\n",
      " 11421 Ace      angiotensin I converting enzyme (peptidyl-dipeptidase A) 1\n"
     ]
    }
   ],
   "source": [
    "print('GeneID Symbol   Description')\n",
    "print('------ -------  --------------------------------------------------------')\n",
    "for nt_gene in nts[:20]:\n",
    "    print('{GeneID:6} {Symbol:8} {description}'.format(**nt_gene._asdict()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3e) Create a symbol2nt dict"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "22,216 gene symbols\n"
     ]
    }
   ],
   "source": [
    "sym2nt = {nt.Symbol:nt for nt in nts}\n",
    "print('{N:,} gene symbols'.format(N=len(sym2nt)))\n",
    "assert len(nts) == len(sym2nt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 3f) Print NCBI information for a specific gene"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Aliases         ['AW208573', 'CD143']\n",
      "CurrentID       0\n",
      "GeneID          11421\n",
      "OMIM            []\n",
      "Org_name        Mus musculus\n",
      "Status          live\n",
      "Symbol          Ace\n",
      "chromosome      11\n",
      "description     angiotensin I converting enzyme (peptidyl-dipeptidase A) 1\n",
      "end_position_on_the_genomic_accession 105880790\n",
      "exon_count      26\n",
      "genomic_nucleotide_accession_version NC_000077.7\n",
      "map_location    11 68.84 cM\n",
      "no_hdr0         \n",
      "orientation     plus\n",
      "other_designations angiotensin-converting enzyme|dipeptidyl carboxypeptidase I|dipeptidyl peptidase|kininase II\n",
      "start_position_on_the_genomic_accession 105858774\n",
      "tax_id          10090\n"
     ]
    }
   ],
   "source": [
    "# Choose a specific gene\n",
    "symbol = 'Ace'\n",
    "\n",
    "# Print NCBI information for the chosen gene\n",
    "for field, value in sorted(sym2nt[symbol]._asdict().items()):\n",
    "    print('{FLD:15} {VAL:}'.format(FLD=field, VAL=value))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (C) 2016-present, DV Klopfenstein, H Tang. All rights reserved."
   ]
  }
 ],
 "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}