{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run a Gene Ontology Enrichment Analysis (GOEA)\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",
    "Note: you must have the Python packages, **xlrd**, **pandas** and **openpyxl**, installed to run this example. \n",
    "\n",
    "Note: To create plots, you must have:\n",
    "  * Python packages: **pyparsing**, **pydot**\n",
    "  * [Graphviz](http://www.graphviz.org/) loaded and your PATH environmental variable pointing to the Graphviz bin directory."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Download Ontologies and Associations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1a. 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": [
    "### 1b. Download Associations, if necessary\n",
    "The NCBI gene2go file contains numerous species. We will select mouse shortly."
   ]
  },
  {
   "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",
    "fin_gene2go = download_ncbi_associations()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Load Ontologies, Associations and Background gene set "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2a. Load Ontologies"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "go-basic.obo: fmt(1.2) rel(2020-01-01) 47,337 GO Terms\n"
     ]
    }
   ],
   "source": [
    "from goatools.obo_parser import GODag\n",
    "\n",
    "obodag = GODag(\"go-basic.obo\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2b. Load Associations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:07.052179 367,335 annotations, 24,267 genes, 18,190 GOs, 1 taxids READ: gene2go \n",
      "BP 17,859 annotated mouse genes\n",
      "MF 16,721 annotated mouse genes\n",
      "CC 18,824 annotated mouse genes\n"
     ]
    }
   ],
   "source": [
    "from __future__ import print_function\n",
    "from goatools.anno.genetogo_reader import Gene2GoReader\n",
    "\n",
    "# Read NCBI's gene2go. Store annotations in a list of namedtuples\n",
    "objanno = Gene2GoReader(fin_gene2go, taxids=[10090])\n",
    "\n",
    "# Get namespace2association where:\n",
    "#    namespace is:\n",
    "#        BP: biological_process               \n",
    "#        MF: molecular_function\n",
    "#        CC: cellular_component\n",
    "#    assocation is a dict:\n",
    "#        key: NCBI GeneID\n",
    "#        value: A set of GO IDs associated with that gene\n",
    "ns2assoc = objanno.get_ns2assc()\n",
    "\n",
    "for nspc, id2gos in ns2assoc.items():\n",
    "    print(\"{NS} {N:,} annotated mouse genes\".format(NS=nspc, N=len(id2gos)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 2c. Load Background gene set\n",
    "In this example, the background is all mouse protein-codinge genes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "26033\n"
     ]
    }
   ],
   "source": [
    "from goatools.test_data.genes_NCBI_10090_ProteinCoding import GENEID2NT as GeneID2nt_mus\n",
    "print(len(GeneID2nt_mus))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Initialize a GOEA object\n",
    "The GOEA object holds the Ontologies, Associations, and background.    \n",
    "Numerous studies can then be run withough needing to re-load the above items.    \n",
    "In this case, we only run one GOEA.    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Load BP Gene Ontology Analysis ...\n",
      "fisher module not installed.  Falling back on scipy.stats.fisher_exact\n",
      " 65% 16,811 of 26,033 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",
      " 70% 18,198 of 26,033 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",
      " 63% 16,331 of 26,033 population items found in association\n"
     ]
    }
   ],
   "source": [
    "from goatools.goea.go_enrichment_ns import GOEnrichmentStudyNS\n",
    "\n",
    "goeaobj = GOEnrichmentStudyNS(\n",
    "        GeneID2nt_mus.keys(), # List of mouse protein-coding genes\n",
    "        ns2assoc, # geneid/GO associations\n",
    "        obodag, # Ontologies\n",
    "        propagate_counts = False,\n",
    "        alpha = 0.05, # default significance cut-off\n",
    "        methods = ['fdr_bh']) # defult multipletest correction method\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Read study genes\n",
    "~400 genes from the Nature paper supplemental table 4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "400 genes READ: /mnt/c/Users/note2/Data/git/goatools/goatools/test_data/nbt_3102/nbt.3102-S4_GeneIDs.xlsx\n"
     ]
    }
   ],
   "source": [
    "import os\n",
    "import goatools\n",
    "# Get xlsx filename where data is stored\n",
    "ROOT = os.path.dirname(goatools.__file__) # goatool root directory\n",
    "din_xlsx = os.path.join(ROOT, \"test_data\", \"nbt_3102\", \"nbt.3102-S4_GeneIDs.xlsx\")\n",
    "# Read data\n",
    "if os.path.isfile(din_xlsx):  \n",
    "    import pandas as pd\n",
    "    df = pd.read_excel(\n",
    "        din_xlsx,\n",
    "        header=None, \n",
    "        names=[\"symbol\", \"geneid\", \"pval\"], \n",
    "        dtype={\"symbol\": str, \"geneid\": int, \"pval\": float},\n",
    "        index_col=1,\n",
    "    )  # requires openpyxl\n",
    "    geneid2symbol = df[\"symbol\"].to_dict()\n",
    "    print('{N} genes READ: {XLSX}'.format(N=len(geneid2symbol), XLSX=din_xlsx))\n",
    "else:\n",
    "    raise RuntimeError('FILE NOT FOUND: {XLSX}'.format(XLSX=din_xlsx))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5. Run Gene Ontology Enrichment Analysis (GOEA)\n",
    "You may choose to keep all results or just the significant results. In this example, we choose to keep only the significant results."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Run BP Gene Ontology Analysis: current study set of 400 IDs ... 94%    357 of    381 study items found in association\n",
      " 95%    381 of    400 study items found in population(26033)\n",
      "Calculating 12,254 uncorrected p-values using fisher_scipy_stats\n",
      "  12,254 GO terms are associated with 16,811 of 26,033 population items\n",
      "   2,085 GO terms are associated with    357 of    400 study items\n",
      "  METHOD fdr_bh:\n",
      "      55 GO terms found significant (< 0.05=alpha) ( 53 enriched +   2 purified): statsmodels fdr_bh\n",
      "     209 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 ... 99%    376 of    381 study items found in association\n",
      " 95%    381 of    400 study items found in population(26033)\n",
      "Calculating 1,724 uncorrected p-values using fisher_scipy_stats\n",
      "   1,724 GO terms are associated with 18,198 of 26,033 population items\n",
      "     449 GO terms are associated with    376 of    400 study items\n",
      "  METHOD fdr_bh:\n",
      "      83 GO terms found significant (< 0.05=alpha) ( 83 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 ... 89%    339 of    381 study items found in association\n",
      " 95%    381 of    400 study items found in population(26033)\n",
      "Calculating 4,145 uncorrected p-values using fisher_scipy_stats\n",
      "   4,145 GO terms are associated with 16,331 of 26,033 population items\n",
      "     580 GO terms are associated with    339 of    400 study items\n",
      "  METHOD fdr_bh:\n",
      "      54 GO terms found significant (< 0.05=alpha) ( 52 enriched +   2 purified): statsmodels fdr_bh\n",
      "     276 study items associated with significant GO IDs (enriched)\n",
      "       0 study items associated with significant GO IDs (purified)\n"
     ]
    }
   ],
   "source": [
    "# 'p_' means \"pvalue\". 'fdr_bh' is the multipletest method we are currently using.\n",
    "geneids_study = geneid2symbol.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": [
    "### 5a. Quietly Run Gene Ontology Enrichment Analysis (GOEA)\n",
    "GOEAs can be run quietly using `prt=None`:\n",
    "```\n",
    "goea_results = goeaobj.run_study(geneids_study, prt=None)\n",
    "```\n",
    "#### No output is printed if `prt=None`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "goea_quiet_all = goeaobj.run_study(geneids_study, prt=None)\n",
    "goea_quiet_sig = [r for r in goea_results_all if r.p_fdr_bh < 0.05]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Print customized results summaries\n",
    "##### Example 1: Significant v All GOEA results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "192 of 18,123 results were significant\n"
     ]
    }
   ],
   "source": [
    "print('{N} of {M:,} results were significant'.format(\n",
    "    N=len(goea_quiet_sig),\n",
    "    M=len(goea_quiet_all)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example 2: Enriched v Purified GOEA results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Significant results: 188 enriched, 4 purified\n"
     ]
    }
   ],
   "source": [
    "print('Significant results: {E} enriched, {P} purified'.format(\n",
    "    E=sum(1 for r in goea_quiet_sig if r.enrichment=='e'),\n",
    "    P=sum(1 for r in goea_quiet_sig if r.enrichment=='p')))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example 3: Significant GOEA results by namespace"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Significant results[192] = 55 BP + 54 MF + 83 CC\n"
     ]
    }
   ],
   "source": [
    "import collections as cx\n",
    "ctr = cx.Counter([r.NS for r in goea_quiet_sig])\n",
    "print('Significant results[{TOTAL}] = {BP} BP + {MF} MF + {CC} CC'.format(\n",
    "    TOTAL=len(goea_quiet_sig),\n",
    "    BP=ctr['BP'],  # biological_process\n",
    "    MF=ctr['MF'],  # molecular_function\n",
    "    CC=ctr['CC'])) # cellular_component"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6. Write results to an Excel file and to a text file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    192 items WROTE: nbt3102.xlsx\n",
      "    192 GOEA results for   375 study items. WROTE: nbt3102.txt\n"
     ]
    }
   ],
   "source": [
    "goeaobj.wr_xlsx(\"nbt3102.xlsx\", goea_results_sig)\n",
    "goeaobj.wr_txt(\"nbt3102.txt\", goea_results_sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 7. Plot all significant GO terms\n",
    "Plotting all significant GO terms produces a messy spaghetti plot. Such a plot can be useful sometimes because you can open it and zoom and scroll around. But sometimes it is just too messy to be of use.\n",
    "\n",
    "The **\"{NS}\"** in **\"nbt3102_{NS}.png\"** indicates that you will see three plots, one for \"biological_process\"(BP), \"molecular_function\"(MF), and \"cellular_component\"(CC)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "   55 usr 444 GOs  WROTE: nbt3102_BP.png\n",
      "   83 usr 149 GOs  WROTE: nbt3102_CC.png\n",
      "   54 usr 156 GOs  WROTE: nbt3102_MF.png\n"
     ]
    }
   ],
   "source": [
    "from goatools.godag_plot import plot_gos, plot_results, plot_goid2goobj\n",
    "\n",
    "plot_results(\"nbt3102_{NS}.png\", goea_results_sig)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7a. These plots are likely to messy\n",
    "The *Cellular Component* plot is the smallest plot...\n",
    "![BIG CC PLOT](images/nbt3102_CC.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7b. So make a smaller sub-plot\n",
    "This plot contains GOEA results:\n",
    "  * GO terms colored by P-value:\n",
    "    * pval < 0.005 (light red)\n",
    "    * pval < 0.01 (light orange)\n",
    "    * pval < 0.05 (yellow)\n",
    "    * pval > 0.05 (grey) Study terms that are not statistically significant\n",
    "  * GO terms with study gene counts printed. e.g., \"32 genes\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    5 usr  13 GOs  WROTE: nbt3102_MF_RNA_genecnt.png\n"
     ]
    }
   ],
   "source": [
    "# Plot subset starting from these significant GO terms\n",
    "goid_subset = [\n",
    "    'GO:0003723', # MF D04 RNA binding (32 genes)\n",
    "    'GO:0044822', # MF D05 poly(A) RNA binding (86 genes)\n",
    "    'GO:0003729', # MF D06 mRNA binding (11 genes)\n",
    "    'GO:0019843', # MF D05 rRNA binding (6 genes)\n",
    "    'GO:0003746', # MF D06 translation elongation factor activity (5 genes)\n",
    "]\n",
    "plot_gos(\"nbt3102_MF_RNA_genecnt.png\", \n",
    "    goid_subset, # Source GO ids\n",
    "    obodag, \n",
    "    goea_results=goea_results_all) # Use pvals for coloring\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![RNA subplot](images/nbt3102_MF_RNA_genecnt.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 7c. Add study gene Symbols to plot\n",
    "e.g., *11 genes: Calr, Eef1a1, Pabpc1*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    5 usr  13 GOs  WROTE: nbt3102_MF_RNA_Symbols.png\n"
     ]
    }
   ],
   "source": [
    "plot_gos(\"nbt3102_MF_RNA_Symbols.png\", \n",
    "    goid_subset, # Source GO ids\n",
    "    obodag,\n",
    "    goea_results=goea_results_all, # use pvals for coloring\n",
    "    # We can further configure the plot...\n",
    "    id2symbol=geneid2symbol, # Print study gene Symbols, not Entrez GeneIDs\n",
    "    study_items=6, # Only only 6 gene Symbols max on GO terms\n",
    "    items_p_line=3, # Print 3 genes per line\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![RNA subplot](images/nbt3102_MF_RNA_Symbols.png)"
   ]
  },
  {
   "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
}