{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Report counts of GO terms at various levels and depths\n",
    "\n",
    "Reports the number of GO terms at each level and depth.   \n",
    "\n",
    "  * **Level** refers to the length of the shortest path from the top.   \n",
    "  * **Depth** refers to the length of the longest path from the top.\n",
    "\n",
    "See the Gene Ontology Consorium's (GOC) advice regarding \n",
    "[levels and depths of a GO term](http://geneontology.org/faq/how-can-i-calculate-level-go-term)\n",
    "    \n",
    "## GO level and depth reporting\n",
    "\n",
    "GO terms reported can be all GO terms in an ontology.     \n",
    "Or subsets of GO terms can be reported.     \n",
    "GO subset examples include all GO terms annotated for a species or all GO terms in a study.\n",
    "\n",
    "Example report on full Ontology:\n",
    "\n",
    "```\n",
    "go-basic.obo: fmt(1.2) rel(2019-01-12) 47,374 GO Terms\n",
    "\n",
    "Summary for all Ontologies:\n",
    "Dep <-Depth Counts->  <-Level Counts->\n",
    "Lev   BP    MF    CC    BP    MF    CC\n",
    "--- ----  ----  ----  ----  ----  ----\n",
    "00     1     1     1     1     1     1\n",
    "01    29    16    21    29    16    21\n",
    "02   264   125   345   421   154   746\n",
    "03  1273   570   494  2205   866  1073\n",
    "04  2376  1516   735  4825  2072  1359\n",
    "05  3692  4801   913  7297  5035   697\n",
    "06  4474  1834   787  7287  1934   230\n",
    "07  4699  1029   600  4696   728    68\n",
    "08  4214   508   254  2018   194    10\n",
    "09  3516   312    51   631    79     1\n",
    "10  2399   153     4   241    13     0\n",
    "11  1511   140     1    38    19     0\n",
    "12   854    42     0     0     0     0\n",
    "13   303    35     0     0     0     0\n",
    "14    66    21     0     0     0     0\n",
    "15    14     7     0     0     0     0\n",
    "16     4     1     0     0     0     0\n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Download Ontologies, if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  EXISTS: go-basic.obo\n"
     ]
    }
   ],
   "source": [
    "# Get http://geneontology.org/ontology/go-basic.obo\n",
    "from goatools.base import download_go_basic_obo\n",
    "obo_fname = download_go_basic_obo()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Download Associations, if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  EXISTS: gene2go\n"
     ]
    }
   ],
   "source": [
    "# Get 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": [
    "## 3. Initialize GODag object"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "go-basic.obo: fmt(1.2) rel(2021-05-01) 47,284 GO Terms\n"
     ]
    }
   ],
   "source": [
    "from goatools.obo_parser import GODag\n",
    "\n",
    "obodag = GODag(\"go-basic.obo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Initialize Reporter class"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from goatools.rpt.rpt_lev_depth import RptLevDepth\n",
    "\n",
    "rptobj = RptLevDepth(obodag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Generate depth/level report for all GO terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Dep <-Depth Counts->  <-Level Counts->\n",
      "Lev   BP    MF    CC    BP    MF    CC\n",
      "--- ----  ----  ----  ----  ----  ----\n",
      "00     1     1     1     1     1     1\n",
      "01    30    30     3    30    30     3\n",
      "02   278   190   788   430   254   788\n",
      "03  1201   593   772  2215   977   942\n",
      "04  2351  1545   933  4841  2050  1098\n",
      "05  3626  4815   675  6859  4984   623\n",
      "06  4557  1922   465  6796  1947   386\n",
      "07  4727  1158   302  4628   704   199\n",
      "08  4285   519   159  1977   183    83\n",
      "09  3433   230    59   614    39    43\n",
      "10  2129    87    18   209     0    10\n",
      "11  1158    72     1    42     0     0\n",
      "12   608     7     0     0     0     0\n",
      "13   203     0     0     0     0     0\n",
      "14    42     0     0     0     0     0\n",
      "15     9     0     0     0     0     0\n",
      "16     4     0     0     0     0     0\n"
     ]
    }
   ],
   "source": [
    "rptobj.write_summary_cnts_all()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Count of GO terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "All terms: 47284, Unique terms: 43987\n"
     ]
    }
   ],
   "source": [
    "all_terms = obodag.values()\n",
    "all_terms_unique = set(all_terms)\n",
    "print(f\"All terms: {len(all_terms)}, Unique terms: {len(all_terms_unique)}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "MF: 11168, CC: 4175, BP: 28641, Total terms: 43987\n"
     ]
    }
   ],
   "source": [
    "MF = len(obodag[\"GO:0003674\"].get_all_children())\n",
    "CC = len(obodag[\"GO:0005575\"].get_all_children())\n",
    "BP = len(obodag[\"GO:0008150\"].get_all_children())\n",
    "total = MF + CC + BP + 3 # Should match \"Unique terms\" above\n",
    "print(f\"MF: {MF}, CC: {CC}, BP: {BP}, Total terms: {total}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GO:0003887\tlevel-04\tdepth-06\tDNA-directed DNA polymerase activity [molecular_function]\n",
      "This term shows up 17 times in GoDag\n"
     ]
    }
   ],
   "source": [
    "# Find some examples of duplicates\n",
    "from collections import Counter\n",
    "\n",
    "counter = Counter(all_terms)\n",
    "(most_common, most_common_count), = counter.most_common(1)\n",
    "print(most_common)\n",
    "print(f\"This term shows up {most_common_count} times in GoDag\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Alternative ids: ['GO:0003887', 'GO:0016449', 'GO:0016450', 'GO:0016451', 'GO:0016000', 'GO:0016452', 'GO:0003889', 'GO:0003890', 'GO:0003891', 'GO:0015999', 'GO:0003888', 'GO:0019984', 'GO:0016448', 'GO:0003895', 'GO:0003893', 'GO:0003894', 'GO:0008723']\n"
     ]
    }
   ],
   "source": [
    "# Find which keys led to GO:0003887\n",
    "go_0003887 = [k for k, v in obodag.items() if v == most_common]\n",
    "print(f\"Alternative ids: {go_0003887}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GO:0003887\tlevel-04\tdepth-06\tDNA-directed DNA polymerase activity [molecular_function]\n"
     ]
    }
   ],
   "source": [
    "print(obodag[\"GO:0003888\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When we do `obodag.values` or `obodag.items`, we end up showing all the alternative ids. In the above example, we see \"GO:0003887\" 17 times, due to its 16 alternative ids. We can also check the online version of [GO:0003887](https://salivaryproteome.nidcr.nih.gov/public/index.php/Special:Ontology_Term/GO:0003887) and confirm these ids."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (C) 2016-2019, 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.8.8"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}