{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# How to get _Descendants Count_ and _Information Content_\n",
    "\n",
    "  1. Download Ontologies, if necessary\n",
    "  2. Choose a set of GO IDs to print\n",
    "  3. GoSubDag contains descendants count\n",
    "  4. Print the **descendants count** (_dcnt_)\n",
    "  5. Compare **information content** (_tinfo_) for human, mouse, and fly\n",
    "  6. Print both **descendants count** and **information content**: _dcnt_ and _tinfo_\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",
      "go-basic.obo: fmt(1.2) rel(2019-05-09) 47,407 GO Terms\n"
     ]
    }
   ],
   "source": [
    "from goatools.base import get_godag\n",
    "godag = get_godag(\"go-basic.obo\")\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Choose a set of GO IDs to print\n",
    "### 2a. Choose a GO IDs related to bacteria"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "24 bacteria GO terms\n"
     ]
    }
   ],
   "source": [
    "# Choose a deep leaf-level GO ID associated with \"bacteria\"\n",
    "DESC = 'bacteria'            # GO Term name contains this\n",
    "NSPC = 'cellular_component'  # Desired namespace\n",
    "\n",
    "# Create a chooser function which returns True or False\n",
    "def chooser(goterm):\n",
    "    \"\"\"Choose a leaf-level GO term based on its name\"\"\"\n",
    "    b_match = DESC in goterm.name\n",
    "    # True if GO term is leaf-level (has no children)\n",
    "    b_leaf = not goterm.children\n",
    "    # True if GO term is in 'cellular_component' namespace (nspc)\n",
    "    b_nspc = goterm.namespace == NSPC\n",
    "    return b_match and b_leaf and b_nspc\n",
    "\n",
    "# Get GO terms with desired name in desired GO DAG branch\n",
    "go_ids_selected = set(o.item_id for o in godag.values() if chooser(o))\n",
    "\n",
    "print('{N} {desc} GO terms'.format(N=len(go_ids_selected), desc=DESC))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Get the GO-DAG subset for your GO IDs"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "INITIALIZING GoSubDag:  24 sources in  60 GOs rcnt(True). 0 alt GO IDs\n",
      "             GoSubDag: namedtuple fields: NS level depth GO alt GO_name dcnt D1 id\n",
      "             GoSubDag: relationships: set()\n"
     ]
    }
   ],
   "source": [
    "from goatools.gosubdag.gosubdag import GoSubDag\n",
    "gosubdag = GoSubDag(go_ids_selected, godag)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2c. Get the deepest GO ID in the GO DAG subset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GO:1990061 bacterial degradosome\n"
     ]
    }
   ],
   "source": [
    "go_id, go_term = max(gosubdag.go2obj.items(), key=lambda t: t[1].depth)\n",
    "\n",
    "# Print GO ID, using print format in gosubdag\n",
    "print(go_id, go_term.name)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2d. Get all parents of the deepest GO ID"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "10 ancestors for GO:1990061 \"bacterial degradosome\"\n"
     ]
    }
   ],
   "source": [
    "go_ids_chosen = go_term.get_all_parents()\n",
    "print('{N} ancestors for {GO} \"{name}\"'.format(\n",
    "    N=len(go_ids_chosen), GO=go_term.item_id, name=go_term.name))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. GoSubDag contains descendants count\n",
    "\n",
    "### gosubdag.go2nt\n",
    "The data member, _**go2nt**_, of the class, _**GoSubDag**_, is a dict where:\n",
    "  * _**go**_ is the GO ID (e.g., GO:1990061)    \n",
    "  * _**nt**_ is the namedtuple\n",
    "\n",
    "\n",
    "### The namedtuple field, _dcnt_, contains the descendants count.\n",
    "\n",
    "Additional namedtuple fields:\n",
    "\n",
    "| field   | description\n",
    "|---------|-----------------\n",
    "| NS      | Namespace: BP, MF, or CC\n",
    "| level   | Minimum path from the top of the branch\n",
    "| depth   | Maximum path from the top of the branch\n",
    "| GO      | GO ID\n",
    "| GO_name | GO Name - A short description\n",
    "| dcnt    | Descendants count"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add the deep GO ID to its list of ancestors for printing\n",
    "go_ids_chosen.add(go_id)\n",
    "nts = [gosubdag.go2nt[go] for go in go_ids_chosen]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Print the descendants count (_dcnt_)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "IDX NS GO ID      Descendants Count Depth Name\n",
      "--- -- ---------- ----------------- ----- --------------------\n",
      " 1) CC GO:0005575        4206        D00  cellular_component\n",
      " 2) CC GO:0044464        3321        D01  cell part\n",
      " 3) CC GO:0032991        2116        D01  protein-containing complex\n",
      " 4) CC GO:0044424        2376        D02  intracellular part\n",
      " 5) CC GO:1902494         546        D02  catalytic complex\n",
      " 6) CC GO:0044444        1271        D03  cytoplasmic part\n",
      " 7) CC GO:1905354           8        D03  exoribonuclease complex\n",
      " 8) CC GO:0044445          87        D04  cytosolic part\n",
      " 9) CC GO:0000178           5        D04  exosome (RNase complex)\n",
      "10) CC GO:0000177           2        D05  cytoplasmic exosome (RNase complex)\n",
      "11) CC GO:1990061           0        D06  bacterial degradosome\n"
     ]
    }
   ],
   "source": [
    "fmt_str = '{I:2}) {NS} {GO:10} {dcnt:11}        D{depth:02}  {GO_name}'\n",
    "\n",
    "# Print selected GO information\n",
    "print('IDX NS GO ID      Descendants Count Depth Name')\n",
    "print('--- -- ---------- ----------------- ----- --------------------')\n",
    "for idx, nt_go in enumerate(sorted(nts, key=lambda nt: nt.depth), 1):\n",
    "    print(fmt_str.format(I=idx, **nt_go._asdict()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$ goatools go_plot GO:1990061\n",
    "<img src=\"images/GO_1990061_bacterial_degradosome.png\" width=\"70%\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Compare information content (_tinfo_) for human, mouse, and fly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Short alias, indicating species and the annotation filename\n",
    "abc2gpad = {\n",
    "    'hsa': 'goa_human.gpad',  # human\n",
    "    'mmu': 'mgi.gpad',        # mouse\n",
    "    'dme': 'fb.gpad'          # fly\n",
    "}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5a. Download annotations, if necessary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "for abc, gpad in abc2gpad.items():\n",
    "    if not os.path.exists(gpad):\n",
    "        fin_gz = '{FILE}.gz'.format(FILE=gpad)\n",
    "        wget_gz = 'current.geneontology.org/annotations/{GZ}'.format(GZ=fin_gz)\n",
    "        !wget $wget_gz\n",
    "        !gunzip $fin_gz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5b. Read annotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:15.276513 479,437 annotations READ: goa_human.gpad \n",
      "HMS:0:00:13.508003 409,743 annotations READ: mgi.gpad \n",
      "HMS:0:00:04.770149 114,704 annotations READ: fb.gpad \n"
     ]
    }
   ],
   "source": [
    "from goatools.anno.factory import get_objanno\n",
    "# Load each annotation into an annotation object\n",
    "abc2objanno = {s:get_objanno(f, 'gpad', godag=godag) for s, f in abc2gpad.items()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5c. Get TermCounts object for the annotation of each species"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from goatools.semantic import TermCounts\n",
    "abc2objtcnt = {s:TermCounts(godag, o.get_id2gos_nss()) for s, o in abc2objanno.items()}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 5d. Print _tinfo_ for GO IDs for human, mouse, and fly"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "                              |<----- tinfo --->|\n",
      "IDX NS GO ID      dcnt Depth   hsa    mmu    dme  Name\n",
      "--- -- ---------- ---- ------ ------ ------ ----- --------------------\n",
      " 1) CC GO:0005575 4206  D00    2.846  2.755  2.411 cellular_component\n",
      " 2) CC GO:0044464 3321  D01    2.937  2.912  2.589 cell part\n",
      " 3) CC GO:0032991 2116  D01    4.059  3.983  3.543 protein-containing complex\n",
      " 4) CC GO:0044424 2376  D02    3.096  3.086  2.734 intracellular part\n",
      " 5) CC GO:1902494  546  D02    5.488  5.428  4.713 catalytic complex\n",
      " 6) CC GO:0044444 1271  D03    3.501  3.583  3.403 cytoplasmic part\n",
      " 7) CC GO:1905354    8  D03    9.469  9.437  8.937 exoribonuclease complex\n",
      " 8) CC GO:0044445   87  D04    7.147  6.904  6.012 cytosolic part\n",
      " 9) CC GO:0000178    5  D04    9.510  9.437  8.937 exosome (RNase complex)\n",
      "10) CC GO:0000177    2  D05   10.049  9.976  9.379 cytoplasmic exosome (RNase complex)\n",
      "11) CC GO:1990061    0  D06    0.000  0.000  0.000 bacterial degradosome\n"
     ]
    }
   ],
   "source": [
    "from goatools.semantic import get_info_content\n",
    "fmt_str = ('{I:2}) {NS} {GO:10} {dcnt:4}  D{depth:02}   '\n",
    "           '{hsa:6.3f} {mmu:6.3f} {dme:6.3f} '\n",
    "           '{GO_name}')\n",
    "\n",
    "# Print selected GO information\n",
    "print('                              |<----- tinfo --->|')\n",
    "print('IDX NS GO ID      dcnt Depth   hsa    mmu    dme  Name')\n",
    "print('--- -- ---------- ---- ------ ------ ------ ----- --------------------')\n",
    "for idx, nt_go in enumerate(sorted(nts, key=lambda nt: nt.depth), 1):\n",
    "    abc2tinfo = {s: get_info_content(nt_go.GO, o) for s, o in abc2objtcnt.items()}\n",
    "    print(fmt_str.format(I=idx, **abc2tinfo, **nt_go._asdict()))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Note: _tinfo_ is 0.000 for 'bacterial degradosome' in human mouse and fly because GO:1990061 is not found in the annotations downloaded May 29, 2019"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Print both descendants count and information content: dcnt and tinfo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GPAD ANNOTATION FILENAME: goa_human.gpad\n"
     ]
    }
   ],
   "source": [
    "gpad_filename = abc2gpad['hsa']\n",
    "print('GPAD ANNOTATION FILENAME: {GPAD}'.format(GPAD=gpad_filename))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:14.915225 479,437 annotations READ: goa_human.gpad \n",
      "INITIALIZING GoSubDag:  24 sources in  60 GOs rcnt(True). 0 alt GO IDs\n",
      "             GoSubDag: namedtuple fields: NS level depth GO alt GO_name dcnt D1 tcnt tfreq tinfo id\n",
      "             GoSubDag: relationships: set()\n"
     ]
    }
   ],
   "source": [
    "# Re-read annotation file for clarity in this tutorial\n",
    "annoobj = get_objanno(gpad_filename, 'gpad', godag=godag)\n",
    "\n",
    "# Get associations (Gene ID - to - set of GO IDs) for all namespaces (BP, MF, and CC)\n",
    "id2goids = annoobj.get_id2gos_nss()\n",
    "\n",
    "# Create TermCounts object\n",
    "tcntobj = TermCounts(godag, id2goids)\n",
    "\n",
    "# Manage a subset of GO IDs using GoSubDag\n",
    "# This time, load the GO term counts\n",
    "gosubdag = GoSubDag(go_ids_selected, godag, tcntobj=tcntobj)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NS GO         dcnt   tinfo depth name\n",
      "-- ---------- ----- ------  ---  -----------------------\n",
      "CC GO:0005575 4206   2.846  D00  cellular_component\n",
      "CC GO:0044464 3321   2.937  D01  cell part\n",
      "CC GO:0032991 2116   4.059  D01  protein-containing complex\n",
      "CC GO:0044424 2376   3.096  D02  intracellular part\n",
      "CC GO:1902494  546   5.488  D02  catalytic complex\n",
      "CC GO:0044444 1271   3.501  D03  cytoplasmic part\n",
      "CC GO:1905354    8   9.469  D03  exoribonuclease complex\n",
      "CC GO:0044445   87   7.147  D04  cytosolic part\n",
      "CC GO:0000178    5   9.510  D04  exosome (RNase complex)\n",
      "CC GO:0000177    2  10.049  D05  cytoplasmic exosome (RNase complex)\n",
      "CC GO:1990061    0   0.000  D06  bacterial degradosome\n"
     ]
    }
   ],
   "source": [
    "# Print both descendants count and information content: dcnt and tinfo\n",
    "\n",
    "fmt_str = '{NS} {GO:10} {dcnt:4}  {tinfo:6.3f}  D{depth:02}  {GO_name}'\n",
    "print('NS GO         dcnt   tinfo depth name')\n",
    "print('-- ---------- ----- ------  ---  -----------------------')\n",
    "nts = gosubdag.prt_goids(go_ids_chosen, prtfmt=fmt_str)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (C) 2016-2019, DV Klopfenstein et al. 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
}