{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Computing basic semantic similarities between GO terms"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Adapted from book chapter written by _Alex Warwick Vesztrocy and Christophe Dessimoz_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this section we look at how to compute semantic similarity between GO terms. First we need to write a function that calculates the minimum number of branches connecting two GO terms."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "go-basic.obo: fmt(1.2) rel(2019-04-17) 47,398 GO Terms\n"
     ]
    }
   ],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "from goatools.obo_parser import GODag\n",
    "godag = GODag(\"go-basic.obo\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GO:0048364\tlevel-04\tdepth-04\troot development [biological_process]\n",
      "GO:0032501\tlevel-01\tdepth-01\tmulticellular organismal process [biological_process]\n"
     ]
    }
   ],
   "source": [
    "go_id3 = 'GO:0048364'\n",
    "go_id4 = 'GO:0032501'\n",
    "print(godag[go_id3])\n",
    "print(godag[go_id4])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's get all the annotations from arabidopsis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "HMS:0:00:08.821464 249,902 annotations READ: /mnt/c/Users/note2/Data/git/tmp/goatools/notebooks/tair.gaf \n",
      "19891 IDs in association branch, BP\n"
     ]
    }
   ],
   "source": [
    "# from goatools.associations import read_gaf\n",
    "# associations = read_gaf(\"tair.gaf\")\n",
    "\n",
    "import os\n",
    "from goatools.associations import dnld_assc\n",
    "fin_gaf = os.path.join(os.getcwd(), \"tair.gaf\")\n",
    "associations = dnld_assc(fin_gaf, godag)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "GO:0008150\n"
     ]
    }
   ],
   "source": [
    "# Find deepest common ancestor\n",
    "from goatools.semantic import deepest_common_ancestor\n",
    "go_root = deepest_common_ancestor([go_id3, go_id4], godag)\n",
    "print(go_root)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Plot the two terms of interest and highlight their deepest common ancestor\n",
    "\n",
    "\n",
    "|color |color  | GO Term    | Description\n",
    "|------|-------|------------|------------------------\n",
    "|blue  |#d5ffff| GO:0008150 | deepest common ancestor\n",
    "|green |#d1ffbd| GO:0048364 | User GO Term\n",
    "|green |#d1ffbd| GO:0032501 | User GO Term\n",
    "\n",
    "```\n",
    "$ scripts/go_plot.py GO:0008150#d5ffff GO:0048364#d1ffbd GO:0032501#d1ffdb -o aaa_lin.png --gaf=tair.gaf\n",
    "\n",
    "go-basic.obo: fmt(1.2) rel(2019-02-07) 47,387 GO Terms\n",
    "  READ      236,943 associations: tair.gaf\n",
    "#d5ffff GO:0008150 # BP 29699  3.30 L00 D00       biological_process\n",
    "#f1fbfd GO:0032502 # BP  3220  5.02 L01 D01 A     developmental process\n",
    "#d1ffdb GO:0032501 # BP  1003  5.48 L01 D01 B     multicellular organismal process\n",
    "        GO:0048856 # BP  1040  5.46 L02 D02 A     anatomical structure development\n",
    "        GO:0099402 # BP    17  6.90 L03 D03 A     plant organ development\n",
    "#d1ffbd GO:0048364 # BP     4  7.56 L04 D04 A     root development\n",
    "```\n",
    "\n",
    "<img src=\"images/nb_lin_go.png\" width=\"600pt\">\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can calculate the semantic distance and semantic similarity, as so:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "The semantic similarity between terms GO:0048364 and GO:0032501 is 0.2.\n"
     ]
    }
   ],
   "source": [
    "from goatools.semantic import semantic_similarity\n",
    "\n",
    "sim = semantic_similarity(go_id3, go_id4, godag)\n",
    "print('The semantic similarity between terms {} and {} is {}.'.format(go_id3, go_id4, sim))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can calculate the information content of the single term, <code>GO:0048364</code>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Information content (GO:0048364) = 7.653820637909629\n"
     ]
    }
   ],
   "source": [
    "from goatools.semantic import TermCounts, get_info_content\n",
    "\n",
    "# First get the counts of each GO term.\n",
    "termcounts = TermCounts(godag, associations)\n",
    "\n",
    "# Calculate the information content\n",
    "go_id = \"GO:0048364\"\n",
    "infocontent = get_info_content(go_id, termcounts)\n",
    "print('Information content ({}) = {}'.format(go_id, infocontent))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resnik's similarity measure is defined as the information content of the most informative common ancestor. That is, the most specific common parent-term in the GO. Then we can calculate this as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Resnik similarity score (GO:0048364, GO:0032501) = 3.2732508872692763\n"
     ]
    }
   ],
   "source": [
    "from goatools.semantic import resnik_sim\n",
    "\n",
    "sim_r = resnik_sim(go_id3, go_id4, godag, termcounts)\n",
    "print('Resnik similarity score ({}, {}) = {}'.format(go_id3, go_id4, sim_r))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lin's similarity measure is defined as:\n",
    "$$ \\textrm{sim}_{\\textrm{Lin}}(t_{1}, t_{2}) = \\frac{2*\\textrm{sim}_{\\textrm{Resnik}}(t_1, t_2)}{IC(t_1) + IC(t_2)} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then we can calculate this as"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Lin similarity score (GO:0048364, GO:0032501) = 0.49819994053182504\n"
     ]
    }
   ],
   "source": [
    "from goatools.semantic import lin_sim\n",
    "\n",
    "sim_l = lin_sim(go_id3, go_id4, godag, termcounts)\n",
    "print('Lin similarity score ({}, {}) = {}'.format(go_id3, go_id4, sim_l))"
   ]
  }
 ],
 "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
}