{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Calculating Annotation Coverage\n",
    "This section shows how to calculate annotation coverage as described here:    \n",
    "\n",
    "    Annotation coverage of Gene Ontology (GO) terms to individual \n",
    "    gene products is high for human or model organisms: \n",
    "      * 87% of ~20k human protein-coding genes have GO annotations\n",
    "      * 76% of ~14k fly protein-coding genes have GO annotations\n",
    "    (Apr 27, 2016)\n",
    "\n",
    "## 1. Download associations\n",
    "NCBI's gene2go file contains annotations of GO terms to Entrez GeneIDs for over 35 different species. We are interested in human and fly which have the taxids 9606 and 7227 repectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  EXISTS: gene2go\n"
     ]
    }
   ],
   "source": [
    "# wget ftp://ftp.ncbi.nlm.nih.gov/gene/DATA/gene2go.gz\n",
    "from goatools.base import download_ncbi_associations\n",
    "gene2go = download_ncbi_associations()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Read associations\n",
    "\n",
    "### 2a. You can read the associations one species at a time..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:07.954549 331,423 annotations, 20,689 genes, 18,627 GOs, 1 taxids READ: gene2go \n",
      "HMS:0:00:04.919504 104,435 annotations, 13,813 genes,  8,683 GOs, 1 taxids READ: gene2go \n"
     ]
    }
   ],
   "source": [
    "from goatools.anno.genetogo_reader import Gene2GoReader\n",
    "\n",
    "geneid2gos_human = Gene2GoReader(gene2go, taxids=[9606])\n",
    "geneid2gos_fly = Gene2GoReader(gene2go, taxids=[7227])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Or you can read 'gene2go' once and load all species..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:12.061128 435,858 annotations, 34,502 genes, 20,059 GOs, 2 taxids READ: gene2go \n"
     ]
    }
   ],
   "source": [
    "from collections import defaultdict, namedtuple\n",
    "\n",
    "geneid2gos_all = Gene2GoReader(gene2go, taxids=[9606, 7227])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Import protein-coding information for human and fly\n",
    "In this example, the background is all human and fly protein-codinge genes.    \n",
    "\n",
    "Follow the instructions in the `background_genes_ncbi` notebook to download a set of background population genes from NCBI."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "19,658 human genes\n",
      "13,968 fly genes\n"
     ]
    }
   ],
   "source": [
    "from genes_ncbi_9606_proteincoding import GENEID2NT as GeneID2nt_human\n",
    "from genes_ncbi_7227_proteincoding import GENEID2NT as GeneID2nt_fly\n",
    "lst = [\n",
    "    (9606, GeneID2nt_human),\n",
    "    (7227, GeneID2nt_fly)\n",
    "]\n",
    "print('{N:,} human genes'.format(N=len(GeneID2nt_human)))\n",
    "print('{N:,} fly genes'.format(N=len(GeneID2nt_fly)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Calculate Gene Ontology coverage\n",
    "Store GO coverage information for *human* and *fly* in the list, **cov_data**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "cov_data = []\n",
    "NtCov = namedtuple(\"NtCov\", \"taxid num_GOs num_covgenes coverage num_allgenes\")\n",
    "for taxid, pcGeneID2nt in lst:\n",
    "    # Get GeneID2GOs association for current species\n",
    "    geneid2gos = geneid2gos_all.get_id2gos_nss(taxid=taxid)\n",
    "    # Restrict GeneID2GOs to only protein-coding genes for this report\n",
    "    pcgene_w_gos = set(geneid2gos.keys()).intersection(set(pcGeneID2nt.keys()))\n",
    "    num_pcgene_w_gos = len(pcgene_w_gos)\n",
    "    num_pc_genes = len(pcGeneID2nt)\n",
    "    # Number of GO terms annotated to protein-coding genes\n",
    "    gos_pcgenes = set()\n",
    "    for geneid in pcgene_w_gos:\n",
    "        gos_pcgenes |= geneid2gos[geneid]\n",
    "    # Print report data\n",
    "    cov_data.append(NtCov(\n",
    "        taxid = taxid,\n",
    "        num_GOs = len(gos_pcgenes),\n",
    "        num_covgenes = num_pcgene_w_gos,\n",
    "        coverage = 100.0*num_pcgene_w_gos/num_pc_genes,\n",
    "        num_allgenes = num_pc_genes))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5 Report Gene Ontology coverage for human and fly\n",
    "Print the *human* and *fly* GO coverage information that is stored in the list, **cov_data**."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " taxid    GOs GeneIDs Coverage\n",
      "------ ------ ------- ----------------------\n",
      "  9606 18,470  18,725  95% GO coverage of 19,658 protein-coding genes\n",
      "  7227  8,587  10,952  78% GO coverage of 13,968 protein-coding genes\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "\n",
    "print(\" taxid    GOs GeneIDs Coverage\")\n",
    "print(\"------ ------ ------- ----------------------\")\n",
    "fmtstr = \"{TAXID:>6} {N:>6,} {M:>7,}  {COV:2.0f}% GO coverage of {TOT:,} protein-coding genes\"\n",
    "for nt in cov_data:\n",
    "    print(fmtstr.format(\n",
    "        TAXID = nt.taxid,\n",
    "        N = nt.num_GOs,\n",
    "        M = nt.num_covgenes,\n",
    "        COV = nt.coverage,\n",
    "        TOT = nt.num_allgenes))"
   ]
  },
  {
   "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": 1
}