{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run a GOEA. Print study genes as either IDs symbols\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 package, **xlrd**, installed to run this example. "
   ]
  },
  {
   "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"
   ]
  },
  {
   "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",
    "file_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:06.334969 367,335 annotations, 24,267 genes, 18,190 GOs, 1 taxids READ: gene2go \n",
      "CC 18,824 annotated mouse genes\n",
      "BP 17,859 annotated mouse genes\n",
      "MF 16,721 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(file_gene2go, taxids=[10090])\n",
    "\n",
    "# Get associations for each branch of the GO DAG (BP, MF, CC)\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.    \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": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "from genes_ncbi_10090_proteincoding import GENEID2NT as 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, # 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": [
      "READ: /mnt/c/Users/note2/Data/git/goatools/goatools/test_data/nbt_3102/nbt.3102-S4_GeneIDs.xlsx\n"
     ]
    }
   ],
   "source": [
    "# Data will be stored in this variable\n",
    "import os\n",
    "geneid2symbol = {}\n",
    "# Get xlsx filename where data is stored\n",
    "ROOT = os.path.dirname(os.getcwd()) # go up 1 level from current working directory\n",
    "din_xlsx = os.path.join(ROOT, \"goatools/test_data/nbt_3102/nbt.3102-S4_GeneIDs.xlsx\")\n",
    "# Read data\n",
    "if os.path.isfile(din_xlsx):  \n",
    "    import xlrd\n",
    "    book = xlrd.open_workbook(din_xlsx)\n",
    "    pg = book.sheet_by_index(0)\n",
    "    for r in range(pg.nrows):\n",
    "        symbol, geneid, pval = [pg.cell_value(r, c) for c in range(pg.ncols)]\n",
    "        if geneid:\n",
    "            geneid2symbol[int(geneid)] = symbol\n",
    "    print('READ: {XLSX}'.format(XLSX=din_xlsx))\n",
    "else:\n",
    "    raise RuntimeError('CANNOT READ: {XLSX}'.format(XLSX=fin_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": [
    "## 6. Write results to an Excel file and to a text file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    192 items WROTE: nbt3102_symbols.xlsx\n",
      "    192 items WROTE: nbt3102_geneids.xlsx\n"
     ]
    }
   ],
   "source": [
    "goeaobj.wr_xlsx(\"nbt3102_symbols.xlsx\", goea_results_sig, itemid2name=geneid2symbol)\n",
    "goeaobj.wr_xlsx(\"nbt3102_geneids.xlsx\", goea_results_sig)"
   ]
  },
  {
   "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
}