notebooks/goea_nbt3102_group_results.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Group Significant GO terms by Frequently Seen Words\n",
    "\n",
    "We use data from a 2014 Nature paper:    \n",
    "[Computational analysis of cell-to-cell heterogeneity\n",
    "in single-cell RNA-sequencing data reveals hidden \n",
    "subpopulations of cells\n",
    "](http://www.nature.com/nbt/journal/v33/n2/full/nbt.3102.html#methods)\n",
    "    \n",
    "This iPython notebook demonstrates one approach to explore Gene Ontology Enrichment Analysis (GOEA) results:\n",
    "1. Create sub-plots containing significant GO terms which share a common word, like **RNA**.\n",
    "2. Create detailed reports showing all significant GO terms and all study gene symbols for the common word.\n",
    "\n",
    "\n",
    "## Summary\n",
    "The code in this notebook generates the data used in this statement in the GOATOOLS manuscript:\n",
    "\n",
    "    We observed:\n",
    "        93 genes associated with RNA,\n",
    "        47 genes associated with translation,\n",
    "        70 genes associated with mitochondrial or mitochondrian, and\n",
    "        37 genes associated with ribosomal, as reported by GOATOOLS.\n",
    "        \n",
    "        \n",
    "## Details\n",
    "Details summarized here are also found in the file, \n",
    "[nbt3102_GO_word_genes.txt](./doc/nbt3102_GO_word_genes.txt), \n",
    "which is generated by this iPython notebook.\n",
    "    \n",
    "* **RNA**: 93 study genes, 6 GOs:    \n",
    "  0) BP GO:0006364 rRNA processing (8 genes)    \n",
    "  1) MF GO:0003723 RNA binding (32 genes)    \n",
    "  2) MF GO:0003729 mRNA binding (11 genes)    \n",
    "  3) MF GO:0008097 5S rRNA binding (4 genes)    \n",
    "  4) MF GO:0019843 rRNA binding (6 genes)    \n",
    "  5) MF GO:0044822 poly(A) RNA binding (86 genes)    \n",
    "    \n",
    "* **translation**: 47 study genes, 5 GOs:         \n",
    "0) BP GO:0006412 translation (41 genes)    \n",
    "1) BP GO:0006414 translational elongation (7 genes)    \n",
    "2) BP GO:0006417 regulation of translation (9 genes)    \n",
    "3) MF GO:0003746 translation elongation factor activity (5 genes)    \n",
    "4) MF GO:0031369 translation initiation factor binding (4 genes)    \n",
    "\n",
    "* **mitochond**: 70 study genes, 8 GOs:     \n",
    "0) BP GO:0051881 regulation of mitochondrial membrane potential (7 genes)    \n",
    "1) CC GO:0000275 mitochondrial proton-transporting ATP synthase complex, catalytic core F(1) (3 genes)    \n",
    "2) CC GO:0005739 mitochondrion (68 genes)    \n",
    "3) CC GO:0005743 mitochondrial inner membrane (28 genes)    \n",
    "4) CC GO:0005747 mitochondrial respiratory chain complex I (5 genes)    \n",
    "5) CC GO:0005753 mitochondrial proton-transporting ATP synthase complex (4 genes)    \n",
    "6) CC GO:0005758 mitochondrial intermembrane space (7 genes)    \n",
    "7) CC GO:0031966 mitochondrial membrane (6 gene0) CC GO:0005925 focal adhesion (53 genes)    \n",
    "\n",
    "* **ribosomal**: 37 study genes, 6 GOs:     \n",
    "0) BP GO:0000028 ribosomal small subunit assembly (9 genes)    \n",
    "1) BP GO:0042274 ribosomal small subunit biogenesis (6 genes)    \n",
    "2) CC GO:0015934 large ribosomal subunit (4 genes)    \n",
    "3) CC GO:0015935 small ribosomal subunit (13 genes)    \n",
    "4) CC GO:0022625 cytosolic large ribosomal subunit (16 genes)    \n",
    "5) CC GO:0022627 cytosolic small ribosomal subunit (19 genes)    \n",
    "\n",
    "Also seen, but not reported in the manuscript:    \n",
    "* **ribosome**: 41 study genes, 2 GOs:    \n",
    "0) CC GO:0005840 ribosome (35 genes)    \n",
    "1) MF GO:0003735 structural constituent of ribosome (38 genes)     \n",
    "\n",
    "* **adhesion**: 53 study genes, 1 GOs:       \n",
    "0) CC GO:0005925 focal adhesion (53 genes)       \n",
    "\n",
    "* **endoplasmic**: 49 study genes, 3 GOs:    \n",
    "0) CC GO:0005783 endoplasmic reticulum (48 genes)    \n",
    "1) CC GO:0005790 smooth endoplasmic reticulum (5 genes)    \n",
    "2) CC GO:0070971 endoplasmic reticulum exit site (4 genes)    \n",
    "\n",
    "* **nucleotide**: 46 study genes, 1 GOs:    \n",
    "0) MF GO:0000166 nucleotide binding (46 genes)        \n",
    "\n",
    "* **apoptotic**: 42 study genes, 2 GOs:    \n",
    "0) BP GO:0006915 apoptotic process (26 genes)    \n",
    "1) BP GO:0043066 negative regulation of apoptotic process (28 genes)        \n",
    "\n",
    "\n",
    "## Methodology    \n",
    "For this exploration, we choose specific sets of GO terms for each plot based on frequently seen words in the GO term name. Examples of GO term names include \"*rRNA processing*\", \"*poly(A) RNA binding*\", and \"*5S rRNA binding*\".  The common word for these GO terms is \"*RNA*\". \n",
    "\n",
    "Steps:\n",
    "1. Run a Gene Ontology Enrichment Analysis.\n",
    "2. Count all words in the significant GO term names.\n",
    "3. Inspect word-count list from step 2.\n",
    "4. Create curated list of words based on frequently seen GO term words.\n",
    "5. Get significant GO terms which contain the words of interest.\n",
    "6. Plot GO terms seen for each word of interest.\n",
    "7. Print a report with full details\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Run GOEA. Save results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  EXISTS: go-basic.obo\n",
      "go-basic.obo: fmt(1.2) rel(2019-04-17) 47,398 GO Terms\n",
      "  EXISTS: gene2go\n",
      "HMS:0:00:06.668335 364,039 annotations READ: gene2go \n",
      "1 taxids stored: 10090\n",
      "\n",
      "Load BP Gene Ontology Analysis ...\n",
      "fisher module not installed.  Falling back on scipy.stats.fisher_exact\n",
      " 59% 16,747 of 28,212 population items found in association\n",
      "\n",
      "Load CC Gene Ontology Analysis ...\n",
      "fisher module not installed.  Falling back on scipy.stats.fisher_exact\n",
      " 65% 18,276 of 28,212 population items found in association\n",
      "\n",
      "Load MF Gene Ontology Analysis ...\n",
      "fisher module not installed.  Falling back on scipy.stats.fisher_exact\n",
      " 58% 16,418 of 28,212 population items found in association\n",
      "\n",
      "Run BP Gene Ontology Analysis: current study set of 400 IDs ...\n",
      " 93%    357 of    382 study items found in association\n",
      " 96%    382 of    400 study items found in population(28212)\n",
      "Calculating 12,189 uncorrected p-values using fisher_scipy_stats\n",
      "  12,189 GO terms are associated with 16,747 of 28,212 population items\n",
      "   2,068 GO terms are associated with    357 of    400 study items\n",
      "  METHOD fdr_bh:\n",
      "      70 GO terms found significant (< 0.05=alpha) ( 68 enriched +   2 purified): statsmodels fdr_bh\n",
      "     230 study items associated with significant GO IDs (enriched)\n",
      "       4 study items associated with significant GO IDs (purified)\n",
      "\n",
      "Run CC Gene Ontology Analysis: current study set of 400 IDs ...\n",
      " 98%    376 of    382 study items found in association\n",
      " 96%    382 of    400 study items found in population(28212)\n",
      "Calculating 1,724 uncorrected p-values using fisher_scipy_stats\n",
      "   1,724 GO terms are associated with 18,276 of 28,212 population items\n",
      "     445 GO terms are associated with    376 of    400 study items\n",
      "  METHOD fdr_bh:\n",
      "      92 GO terms found significant (< 0.05=alpha) ( 92 enriched +   0 purified): statsmodels fdr_bh\n",
      "     373 study items associated with significant GO IDs (enriched)\n",
      "       0 study items associated with significant GO IDs (purified)\n",
      "\n",
      "Run MF Gene Ontology Analysis: current study set of 400 IDs ...\n",
      " 88%    338 of    382 study items found in association\n",
      " 96%    382 of    400 study items found in population(28212)\n",
      "Calculating 4,128 uncorrected p-values using fisher_scipy_stats\n",
      "   4,128 GO terms are associated with 16,418 of 28,212 population items\n",
      "     581 GO terms are associated with    338 of    400 study items\n",
      "  METHOD fdr_bh:\n",
      "      56 GO terms found significant (< 0.05=alpha) ( 54 enriched +   2 purified): statsmodels fdr_bh\n",
      "     273 study items associated with significant GO IDs (enriched)\n",
      "       0 study items associated with significant GO IDs (purified)\n"
     ]
    }
   ],
   "source": [
    "%run goea_nbt3102_fncs.ipynb\n",
    "goeaobj = get_goeaobj_nbt3102('fdr_bh')\n",
    "# Read Nature data from Excel file (~400 study genes)\n",
    "studygeneid2symbol = read_data_nbt3102()\n",
    "# Run Gene Ontology Enrichment Analysis using Benjamini/Hochberg FDR correction\n",
    "geneids_study = studygeneid2symbol.keys()\n",
    "goea_results_all = goeaobj.run_study(geneids_study)\n",
    "goea_results_sig = [r for r in goea_results_all if r.p_fdr_bh < 0.05]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Count all words in the significant GO term names."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2a. Get list of significant GO term names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "218\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "go_names = [r.name for r in goea_results_sig]\n",
    "print(len(go_names)) # Includes ONLY signficant results"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Get word count in significant GO term names"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections as cx\n",
    "word2cnt = cx.Counter([word for name in go_names for word in name.split()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Inspect word-count list generated at step 2.\n",
    "Words like \"mitochondrial\" can be interesting. Some words will not be interesting, such as \"of\"."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[('binding', 40), ('of', 33), ('regulation', 25), ('protein', 23), ('cell', 20), ('complex', 20), ('activity', 18), ('positive', 15), ('process', 12), ('actin', 12)]\n"
     ]
    }
   ],
   "source": [
    "# Print 10 most common words found in significant GO term names\n",
    "print(word2cnt.most_common(10))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Create curated list of words based on frequently seen GO term words."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "freq_seen = ['RNA', 'translation', 'mitochond', 'ribosomal', 'ribosome',\n",
    "             'adhesion', 'endoplasmic', 'nucleotide', 'apoptotic']\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. For each word of interest, create a list of significant GOs whose name contains the word.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Collect significant GOs for words in freq_seen (unordered)\n",
    "word2siggos = cx.defaultdict(set)\n",
    "# Loop through manually curated words of interest\n",
    "for word in freq_seen:\n",
    "    # Check each significant GOEA result for the word of interest\n",
    "    for rec in goea_results_sig:\n",
    "        if word in rec.name:\n",
    "            word2siggos[word].add(rec.GO)\n",
    "# Sort word2gos to have the same order as words in freq_seen\n",
    "word2siggos = cx.OrderedDict([(w, word2siggos[w]) for w in freq_seen])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Plot GO terms seen for each word of interest."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6a. Create a convenient goid-to-goobject dictionary"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "18041\n"
     ]
    }
   ],
   "source": [
    "goid2goobj_all = {nt.GO:nt.goterm for nt in goea_results_all}\n",
    "print(len(goid2goobj_all))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6b. Create plots formed by a shared word in the significant GO term's name"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    6 usr  30 GOs  WROTE: nbt3102_word_RNA.png\n",
      "    6 usr  67 GOs  WROTE: nbt3102_word_translation.png\n",
      "   11 usr  70 GOs  WROTE: nbt3102_word_mitochond.png\n",
      "    6 usr  28 GOs  WROTE: nbt3102_word_ribosomal.png\n",
      "    4 usr  19 GOs  WROTE: nbt3102_word_ribosome.png\n",
      "    1 usr   7 GOs  WROTE: nbt3102_word_adhesion.png\n",
      "    4 usr  34 GOs  WROTE: nbt3102_word_endoplasmic.png\n",
      "    1 usr   7 GOs  WROTE: nbt3102_word_nucleotide.png\n",
      "    2 usr  16 GOs  WROTE: nbt3102_word_apoptotic.png\n"
     ]
    }
   ],
   "source": [
    "# Plot set of GOs for each frequently seen word\n",
    "from goatools.godag_plot import plot_goid2goobj\n",
    "for word, gos in word2siggos.items():\n",
    "    goid2goobj = {go:goid2goobj_all[go] for go in gos}\n",
    "    plot_goid2goobj(\n",
    "        \"nbt3102_word_{WORD}.png\".format(WORD=word),\n",
    "        goid2goobj, # source GOs to plot and their GOTerm object\n",
    "        study_items=15, # Max number of gene symbols to print in each GO term\n",
    "        id2symbol=studygeneid2symbol, # Contains GeneID-to-Symbol from Step 1\n",
    "        goea_results=goea_results_all, # pvals used for GO Term coloring\n",
    "        dpi=150)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 6c. Example plot for \"apoptotic\"\n",
    "\n",
    "**Colors**:\n",
    "Please note that to have colors related to GOEA significance, you must provide the GOEA results, as shown here with the \"goea_results=goea_results_all\" argument. \n",
    "1. **Levels of Statistical Significance**:\n",
    "    1. **light red** => *extremely significant* fdr_bh values (p<0.005)\n",
    "    2. **orange** => *very significant* fdr_bh values (p<0.01)\n",
    "    2. **yellow** => *significant* fdr_bh values (p<0.05)\n",
    "    3. **grey** => study terms which are ***not*** *statistically significant* (p>0.05)\n",
    "2. **High-level GO terms**:\n",
    "    1. **Cyan** => Level-01 GO terms    \n",
    "\n",
    "*Please note* that the variable, *goea_results_all*, contains gene ids and fdr_bh alpha values for **all** study GO terms, significant or not. If the argument had only included the significant results, \"goea_results=goea_results_sig\", the currently colored grey GO terms would be white and would not have study genes annotated inside.\n",
    "\n",
    "**Gene Symbol Names**    \n",
    "Please notice that the study gene symbol names are written in thier associated GO term box.  Symbol names and not gene count nor gene ids are used because of the argument, \"id2symbol=studygeneid2symbol\", to the function, \"plot_goid2goobj\".\n",
    "\n",
    "\n",
    "\n",
    "![apoptotic subplot](images/nbt3102_word_apoptotic.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Print a report with full details\n",
    "### 7a. Create detailed report"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "  WROTE: nbt3102_GO_word_genes.txt\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fout = \"nbt3102_GO_word_genes.txt\"\n",
    "go2res = {nt.GO:nt for nt in goea_results_all}\n",
    "with open(fout, \"w\") as prt:\n",
    "    prt.write(\"\"\"This file is generated by test_nbt3102.py and is intended to confirm\n",
    "this statement in the GOATOOLS manuscript:\n",
    "\n",
    "        We observed:\n",
    "            93 genes associated with RNA,\n",
    "            47 genes associated with translation,\n",
    "            70 genes associated with mitochondrial or mitochondrian, and\n",
    "            37 genes associated with ribosomal, as reported by GOATOOLS.\n",
    "\n",
    "\"\"\")\n",
    "    for word, gos in word2siggos.items():\n",
    "        # Sort first by BP, MF, CC. Sort second by GO id.\n",
    "        gos = sorted(gos, key=lambda go: [go2res[go].NS, go])\n",
    "        genes = set()\n",
    "        for go in gos:\n",
    "            genes |= go2res[go].study_items\n",
    "        genes = sorted([studygeneid2symbol[g] for g in genes])\n",
    "        prt.write(\"\\n{WD}: {N} study genes, {M} GOs\\n\".format(WD=word, N=len(genes), M=len(gos)))\n",
    "        prt.write(\"{WD} GOs: {GOs}\\n\".format(WD=word, GOs=\", \".join(gos)))\n",
    "        for i, go in enumerate(gos):\n",
    "            res = go2res[go]\n",
    "            prt.write(\"{I}) {NS} {GO} {NAME} ({N} genes)\\n\".format(\n",
    "                I=i, NS=res.NS, GO=go, NAME=res.name, N=res.study_count))\n",
    "        prt.write(\"{N} study genes:\\n\".format(N=len(genes)))\n",
    "        N = 10 # Number of genes per line\n",
    "        mult = [genes[i:i+N] for i in range(0, len(genes), N)]\n",
    "        prt.write(\"  {}\\n\".format(\"\\n  \".join([\", \".join(str(g) for g in sl) for sl in mult])))\n",
    "    print(\"  WROTE: {F}\\n\".format(F=fout))\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7b Snippet of Report\n",
    "Shows the word, RNA, and all significant GO terms which contain RNA and all study genes associated with the RNA GO terms.\n",
    "\n",
    "```\n",
    "...\n",
    "\n",
    "RNA: 93 study genes, 6 GOs\n",
    "RNA GOs: GO:0006364, GO:0003723, GO:0003729, GO:0008097, GO:0019843, GO:0044822\n",
    "0) BP GO:0006364 rRNA processing (8 genes)\n",
    "1) MF GO:0003723 RNA binding (32 genes)\n",
    "2) MF GO:0003729 mRNA binding (11 genes)\n",
    "3) MF GO:0008097 5S rRNA binding (4 genes)\n",
    "4) MF GO:0019843 rRNA binding (6 genes)\n",
    "5) MF GO:0044822 poly(A) RNA binding (86 genes)\n",
    "93 study genes:\n",
    "  Ahnak, Aldoa, Anxa2, Arf1, Atp5a1, Auh, Btf3, Calr, Canx, Coro1a\n",
    "  Csde1, Edf1, Eef1a1, Eef2, Eif3h, Eif4e2, Eif5a, Eno1, Fcf1, Fdps\n",
    "  Gdi2, Gltscr2, Gm5506, Gnb2l1, H2-D1, H2-K1, H2-Q10, H2-Q7, Hdlbp, Hnrnpl\n",
    "  Hsp90ab1, Hspa8, Immt, Lgals1, Mdh2, Mrpl45, Ndufv3, Npm1, P4hb, Pabpc1\n",
    "  Park7, Pdia3, Pfn1, Pkm, Poldip3, Psma6, R3hdm1, Rbm39, Rbm5, Rpl13\n",
    "  Rpl13a, Rpl14, Rpl18, Rpl18a, Rpl23, Rpl24, Rpl34, Rpl4, Rpl6, Rpl7\n",
    "  Rpl7a, Rpl8, Rplp0, Rps10, Rps11, Rps14, Rps15, Rps15a, Rps17, Rps19\n",
    "  Rps2, Rps20, Rps24, Rps25, Rps26, Rps3, Rps4x, Rps5, Rps7, Rps9\n",
    "  Rpsa, Samhd1, Slc25a5, Son, Sptbn1, Srp72, Srpk1, Sub1, Suclg1, Sumo1\n",
    "  Syf2, Ywhae, Zc3hav1\n",
    "  \n",
    "...\n",
    "```\n"
   ]
  },
  {
   "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.6.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}