notebooks/Enrichment_analyses_human_phenotype_ontology.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Do enrichment analyses using the Human Phenotype Ontology (HPO)\n",
    "HPO files located in goatools/notebooks/data/hpo"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "fin_study = 'data/hpo/genes.list'         # Study genes\n",
    "fin_pop   = 'data/hpo/gobackground.list'  # Population genes\n",
    "fin_obo   = 'data/hpo/hp.obo'             # DAG containing HPO terms\n",
    "fin_anno  = 'data/hpo/hpo.annotation.tab' # Annotation of genes-to-HPO terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1) Read list of study genes and population genes"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "    2,252 READ: data/hpo/genes.list\n",
      "   14,446 READ: data/hpo/gobackground.list\n"
     ]
    }
   ],
   "source": [
    "from goatools.utils import read_geneset\n",
    "\n",
    "study_ids = read_geneset(fin_study)\n",
    "population_ids = read_geneset(fin_pop)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2) Load the human phenotype ontology DAG"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "data/hpo/hp.obo: fmt(1.2) rel(hp/2021-02-28) 19,498 Terms\n"
     ]
    }
   ],
   "source": [
    "from goatools.obo_parser import GODag\n",
    "\n",
    "godag = GODag(fin_obo)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3) Load the annotations of genes to sets of HPO terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:00.594349 187,934 annotations READ: data/hpo/hpo.annotation.tab \n",
      "4531 IDs in loaded association branch, human_phenotype\n"
     ]
    }
   ],
   "source": [
    "from goatools.anno.idtogos_reader import IdToGosReader\n",
    "\n",
    "annoobj = IdToGosReader(fin_anno, godag=godag)\n",
    "\n",
    "id2gos = annoobj.get_id2gos()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4) Run enrichment analysis on HPO terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "4531 IDs in loaded association branch, human_phenotype\n",
      "\n",
      "Load  Ontology Enrichment Analysis ...\n",
      "Propagating term counts up: is_a\n",
      " 26%  3,823 of 14,446 population items found in association\n"
     ]
    }
   ],
   "source": [
    "from goatools.go_enrichment import GOEnrichmentStudy\n",
    "\n",
    "goeaobj = GOEnrichmentStudy(\n",
    "    population_ids,\n",
    "    annoobj.get_id2gos(),\n",
    "    godag,\n",
    "    methods=['bonferroni', 'fdr_bh'],\n",
    "    pvalcalc='fisher_scipy_stats')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5) Run an enrichment analysis on HPO terms"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Runing  Ontology Analysis: current study set of 2252 IDs.\n",
      " 29%    650 of  2,252 study items found in association\n",
      "100%  2,252 of  2,252 study items found in population(14446)\n",
      "Calculating 9,331 uncorrected p-values using fisher_scipy_stats\n",
      "   9,331 terms are associated with  3,823 of 14,446 population items\n",
      "   5,802 terms are associated with    650 of  2,252 study items\n",
      "  METHOD bonferroni:\n",
      "       2 GO terms found significant (< 0.05=alpha) (  2 enriched +   0 purified): local bonferroni\n",
      "     489 study items associated with significant GO IDs (enriched)\n",
      "       0 study items associated with significant GO IDs (purified)\n",
      "  METHOD fdr_bh:\n",
      "      12 GO terms found significant (< 0.05=alpha) ( 12 enriched +   0 purified): statsmodels fdr_bh\n",
      "     498 study items associated with significant GO IDs (enriched)\n",
      "       0 study items associated with significant GO IDs (purified)\n"
     ]
    }
   ],
   "source": [
    "results = goeaobj.run_study_nts(study_ids)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 6) Print the results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "namespace       term_id  e/p pval_uncorr Benjamimi/Hochberg Bonferroni  study_ratio population_ratio\n",
      "--------------- -------- --- ----------- ------------------ ----------  ----------- ----------------\n",
      "human_phenotype HP:0000007 e    3.40e-09           3.18e-05   3.18e-05     459/2252   2322/14446\n",
      "human_phenotype HP:0000080 e    4.22e-08           1.97e-04   3.94e-04     142/2252    590/14446\n",
      "human_phenotype HP:0000556 e    7.80e-06           2.26e-02   7.28e-02      74/2252    288/14446\n",
      "human_phenotype HP:0000512 e    1.21e-05           2.26e-02   1.13e-01      48/2252    166/14446\n",
      "human_phenotype HP:0000572 e    1.21e-05           2.26e-02   1.13e-01      48/2252    166/14446\n",
      "human_phenotype HP:0012868 e    2.69e-05           3.85e-02   2.51e-01      10/2252     16/14446\n",
      "human_phenotype HP:0003138 e    3.65e-05           3.85e-02   3.40e-01       8/2252     11/14446\n",
      "human_phenotype HP:0031970 e    3.65e-05           3.85e-02   3.40e-01       8/2252     11/14446\n",
      "human_phenotype HP:0000510 e    4.07e-05           3.85e-02   3.80e-01      49/2252    178/14446\n",
      "human_phenotype HP:0030453 e    4.12e-05           3.85e-02   3.85e-01      53/2252    197/14446\n",
      "human_phenotype HP:0032558 e    5.07e-05           3.94e-02   4.73e-01       9/2252     14/14446\n",
      "human_phenotype HP:0032560 e    5.07e-05           3.94e-02   4.73e-01       9/2252     14/14446\n",
      "e: enriched\n",
      "p: purified\n"
     ]
    }
   ],
   "source": [
    "print('namespace       term_id  e/p pval_uncorr Benjamimi/Hochberg Bonferroni  study_ratio population_ratio')\n",
    "print('--------------- -------- --- ----------- ------------------ ----------  ----------- ----------------')\n",
    "pat = '{NS} {GO} {e}    {PVAL:8.2e}           {BH:8.2e}   {BONF:8.2e} {RS:>12} {RP:>12}'\n",
    "for ntd in sorted(results, key=lambda nt: [nt.p_uncorrected, nt.GO]):\n",
    "    if ntd.p_fdr_bh < 0.05:\n",
    "        print(pat.format(\n",
    "            NS=ntd.NS,\n",
    "            GO=ntd.GO,\n",
    "            e=ntd.enrichment,\n",
    "            RS='{}/{}'.format(*ntd.ratio_in_study),\n",
    "            RP='{}/{}'.format(*ntd.ratio_in_pop),\n",
    "            PVAL=ntd.p_uncorrected,\n",
    "            BONF=ntd.p_bonferroni,\n",
    "            BH=ntd.p_fdr_bh))\n",
    "print('e: enriched')\n",
    "print('p: purified')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Copyright (C) 2021-present, DV Klopfenstein and Haibao 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": 2
}