example_plain.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from scipy.stats import fisher_exact\n",
    "from statsmodels.sandbox.stats.multicomp import multipletests"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "project_directory = '/Users/tstoeger/Dropbox/biomedical_data_science_day/'  # <-- modify"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def complete_path(p):\n",
    "    \"\"\"\n",
    "    Completes the path toward a file.\n",
    "    \n",
    "    Input:\n",
    "        p     str, path relative to project path\n",
    "    \n",
    "    Output:\n",
    "        p     str, complete path\n",
    "    \"\"\"\n",
    "    p = os.path.join(project_directory, p)\n",
    "    return p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def read_go_and_turn_to_annotation(taxon_id, subset):\n",
    "# \"\"\"\n",
    "# Reads Gene Ontology (GO) annotation and formats them\n",
    "# toward format anticipated by our enrichment tool.\n",
    "\n",
    "# Input:\n",
    "#     taxon_id    int, NCBI taxonomy identifier \n",
    "#                      (e.g.: 9606 for human)\n",
    "#     subset      int, 'Component', 'Process', \n",
    "#                     or 'Function', or 'all'\n",
    "\n",
    "# Output:\n",
    "#     gene_annotation  df, with columns\n",
    "#                      'gene_entrez' and 'annotation' where\n",
    "#                      'gene_entrez': int with Entrez Gene Identifier\n",
    "#                      'annotation': str with name of annotation\n",
    "# \"\"\"\n",
    "\n",
    "# load tab-delimited file\n",
    "p = complete_path('ncbi_gene/gene2go.gz')   \n",
    "\n",
    "df = #...... inset\n",
    "\n",
    "# filter according to taxon\n",
    "df = df.loc[df.loc[:, '#tax_id']==taxon_id, :]\n",
    "\n",
    "# filter according to sub-category of GO\n",
    "\n",
    "#...... inset\n",
    "\n",
    "# Exclude annotations for which there is some \n",
    "# negative evidence (e.g.: annotators found that\n",
    "# a plausible annotation actually does not occur)\n",
    "# Note that for a given pair of gene and annotation\n",
    "# there may be positive and negative evidence, but\n",
    "# all pairs should be removed if there is any negative\n",
    "# http://geneontology.org/docs/guide-go-evidence-codes/\n",
    "\n",
    "#...... inset\n",
    "\n",
    "# Get future pairs of  'gene_entrez' and 'annotation',\n",
    "# while removing duplicates; and format them according\n",
    "# to the docstring (help/description ) of this function\n",
    "# 'GeneID' will become 'gene_entrez', and 'GO_term'\n",
    "# will become 'annotation'\n",
    "\n",
    "df = df.loc[:, ['GeneID', 'GO_term']].drop_duplicates()\n",
    "\n",
    "df = df.rename(columns={\n",
    "    'GeneID': 'gene_entrez', \n",
    "    'GO_term': 'annotation'}).reset_index(drop=True)\n",
    "df['gene_entrez'] = df['gene_entrez'].astype(int)\n",
    "\n",
    "#     return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "taxon_id = 9606   # the NCBI taxonomy identifier for Homo sapiens\n",
    "gene_annotation = read_go_and_turn_to_annotation(taxon_id, 'all')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gene_annotation.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def functional_enrichment(\n",
    "#     hits, \n",
    "#     background, \n",
    "#     gene_annotation,\n",
    "#     only_genes_with_annotation=True\n",
    "# ):\n",
    "# \"\"\"\n",
    "# Tool to compute functional enrichments using Fisher's\n",
    "# Exact test.\n",
    "\n",
    "# Input:\n",
    "#     hits         list-like, gene_entrez of hits\n",
    "#     background   list-like, gene_entrez of background (all genes in experiment)\n",
    "#     gene_annotation   df, with columns\n",
    "#                      'gene_entrez' and 'annotation' where\n",
    "#                      'gene_entrez': int with Entre Gene Identifier\n",
    "#                      'annotation': str with name of annotation \n",
    "#     only_genes_with_annotation  optional, default: True\n",
    "#                                 if True, hit and backround\n",
    "#                                 will be filtered for genes\n",
    "#                                 in gene_annotation\n",
    "\n",
    "# Output:\n",
    "#     results  df, with index 'annotation' and columns\n",
    "#                 'fold_enrichment' log2 ratio of enrichment of annotation over backround\n",
    "#                 'pvalue' Fisher's exact test on enrichment of annotation\n",
    "#                 'bonferroni' Bonferroni corrected pvalue (multiple testing)\n",
    "#                 'benjamini_hochberg' Benjamini Hochberg corrected pvalue (multiple testing)\n",
    "# \"\"\"\n",
    "\n",
    "\n",
    "# Create three sets:\n",
    "# - hits: the hits from the study\n",
    "# - background: the user supplied list of considered genes\n",
    "# - not hits: genes in background but not hits\n",
    "\n",
    "#...... inset\n",
    "\n",
    "# Optionally, use filter to only consider genes\n",
    "# that have at least one annotation, e.g.: see\n",
    "# https://www.nature.com/articles/s41598-018-19333-x\n",
    "# on gene annotation bias\n",
    "if only_genes_with_annotation:\n",
    "    unique_annotations = set(gene_annotation['gene_entrez'])\n",
    "\n",
    "    hits = hits.intersection(unique_annotations)\n",
    "    background = background.intersection(unique_annotations)\n",
    "    not_hits = not_hits.intersection(unique_annotations)\n",
    "\n",
    "# Force an error if not all hits are in the background\n",
    "if any([hit not in background for hit in hits]):\n",
    "    raise AssertionError(\n",
    "        'At least one hit is not in the supplied background.'\n",
    "    )\n",
    "\n",
    "# Restrict annotations to those that occur in the background\n",
    "gene_annotation = #...... inset\n",
    "\n",
    "# Restrict annotations to those occuring a certain number\n",
    "# Rationale is to get rid of the bulk of annotations that only affect\n",
    "# a very small number of gens <-- note this is somewhat arbritrary\n",
    "# and subjective cutoff\n",
    "minimally_required_candidates_in_annotation = 2\n",
    "c = gene_annotation['annotation'].value_counts()\n",
    "gene_annotation = gene_annotation[gene_annotation['annotation'].isin(\n",
    "    c[c >= minimally_required_candidates_in_annotation].index)]\n",
    "\n",
    "# Prepare to loop through annotations\n",
    "# First, set 'annotation' to become the index of gene_annotation\n",
    "gene_annotation = #...... inset\n",
    "# Second, get all unique 'annotation'\n",
    "annotations = #...... inset\n",
    "# Third, initialize a dataframe with all unique annotations in\n",
    "# rows, and 'fold_enrichment', 'pvalue', 'bonferroni',\n",
    "# 'benjamini_hochberg' as columns\n",
    "results = #...... inset\n",
    "\n",
    "# Loop through annotations\n",
    "for annotation in annotations:\n",
    "\n",
    "    # Set up a contingency table containing \n",
    "    #   [hit_genes_in_annotation, not_hit_genes_in_annotation], \n",
    "    #   [hit_genes_not_in_annotation, not_hit_genes_not_in_annotation]\n",
    "    # where each of above is total number of such genes\n",
    "    # https://www.pathwaycommons.org/guide/primers/statistics/fishers_exact_test/\n",
    "\n",
    "    genes_with_annotation = #...... inset\n",
    "    genes_without_annotation = #...... inset\n",
    "\n",
    "    hit_genes_in_annotation = #...... inset\n",
    "    not_hit_genes_in_annotation = #...... inset\n",
    "\n",
    "    hit_genes_not_in_annotation = #...... inset\n",
    "    not_hit_genes_not_in_annotation = #...... inset\n",
    "\n",
    "    # Perform Fisher's exact test\n",
    "    _, pvalue = fisher_exact(\n",
    "        [\n",
    "             [hit_genes_in_annotation, not_hit_genes_in_annotation], \n",
    "             [hit_genes_not_in_annotation, not_hit_genes_not_in_annotation]],\n",
    "        alternative='greater'\n",
    "    )\n",
    "\n",
    "    # Count number of annotated genes in hits and backround and compute\n",
    "    # the enrichment among hits\n",
    "    in_hits = #...... inset\n",
    "    in_background = #...... inset\n",
    "\n",
    "    fract_of_hits_w_annot = #...... inset\n",
    "    fract_of_background_w_annot = #...... inset\n",
    "\n",
    "    if fract_of_hits_w_annot>0:\n",
    "        enrichment = np.log2(\n",
    "            fract_of_hits_w_annot / fract_of_background_w_annot)\n",
    "    else:\n",
    "        enrichment = -np.inf\n",
    "\n",
    "    results.loc[annotation, 'fold_enrichment'] = enrichment\n",
    "    results.loc[annotation, 'pvalue'] = pvalue\n",
    "\n",
    "# Add corrections for multiple testing\n",
    "results.loc[:, 'bonferroni'] = #...... inset\n",
    "results.loc[:, 'benjamini_hochberg'] = #...... inset\n",
    "\n",
    "# Sort results by pvalue\n",
    "results = results.sort_values('pvalue')\n",
    "\n",
    "#     return results"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get list of hits and background genes to be considered (all in experiment)\n",
    "# in this case we want to know, which annotations are carried by the \n",
    "# genes that most commonly appear in differential gene expression analysis\n",
    "# Data is from Crow et al. 2019 \"Predictablity of human gene expression\",\n",
    "# PNAS, https://www.pnas.org/content/116/13/6491\n",
    "\n",
    "\n",
    "p = complete_path('candidate_genes/crow_et_al_2019/DE_Prior.txt')\n",
    "study_with_genes = pd.read_csv(p, sep='\\t')  # a table that sorts genes according to occurence\n",
    "background_from_study = study_with_genes['Gene_EntrezID']\n",
    "hits_from_study = study_with_genes.loc[:300, 'Gene_EntrezID']  # lets' treat the 300 most frequent genes as hits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "go_enrichments = functional_enrichment(\n",
    "    hits=hits_from_study, \n",
    "    background=background_from_study, \n",
    "    gene_annotation=gene_annotation,\n",
    "    only_genes_with_annotation=False\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "go_enrichments.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# def read_gwas_and_turn_to_annotation():\n",
    "#     \"\"\"\n",
    "#     Reads EBI GWAS annotation and formats them\n",
    "#     toward format anticipated by our enrichment tool.\n",
    "    \n",
    "#     Output:\n",
    "#         gene_annotation  df, with columns\n",
    "#                          'gene_entrez' and 'annotation' where\n",
    "#                          'gene_entrez': int with Entrez Gene Identifier\n",
    "#                          'annotation': str with name of trait\n",
    "#     \"\"\"\n",
    "    \n",
    "#     p = complete_path('ebi_gwas/full.tsv')\n",
    "#     df_gwas = pd.read_csv(\n",
    "#         p, \n",
    "#         sep='\\t', \n",
    "#         usecols=['DISEASE/TRAIT', 'MAPPED_GENE']\n",
    "#     ).drop_duplicates()\n",
    "\n",
    "#     p = complete_path('ncbi_gene/Homo_sapiens.gene_info.gz')\n",
    "#     df_gene_info = pd.read_csv(p, sep='\\t', usecols=['GeneID', 'Symbol'])\n",
    "\n",
    "#     # combine overly simplistically through merge: unambiguous single SNPs \n",
    "#     # -> only for illustrative purposes of tutorial\n",
    "#     df_annotation = pd.merge(\n",
    "#         df_gene_info.rename(columns={'Symbol': 'MAPPED_GENE'}),\n",
    "#         df_gwas,\n",
    "#         how='inner'   \n",
    "#     ).rename(columns={\n",
    "#         'GeneID': 'gene_entrez',\n",
    "#         'DISEASE/TRAIT': 'annotation'})\n",
    "\n",
    "#     df_annotation = df_annotation[['gene_entrez', 'annotation']].drop_duplicates()\n",
    "\n",
    "#     return df_annotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gene_annotation = read_gwas_and_turn_to_annotation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# gwas_enrichments = functional_enrichment(\n",
    "#     hits=hits_from_study, \n",
    "#     background=background_from_study, \n",
    "#     gene_annotation=gene_annotation,\n",
    "#     only_genes_with_annotation=False\n",
    "# )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "\n",
    "# gwas_enrichments.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}