{ "nbformat": 4, "nbformat_minor": 5, "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "name": "python", "version": "3.9.0" } }, "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Information Content for Custom GO Annotations (id2go format)\n", "\n", "This notebook demonstrates how to calculate the **information content (IC)** of GO\n", "terms for any species using a **custom annotation file in `id2go` format**, i.e. a\n", "tab-separated file where:\n", "\n", "* column 1 = gene identifier\n", "* column 2 = semicolon-separated list of GO terms\n", "\n", "This is useful when you want to:\n", "* Analyse non-model organisms not covered by standard GOA databases\n", "* Compare annotation quality produced by different tools\n", "* Identify which GO terms are shared between two annotation sets and how specific they are\n", "\n", "Steps:\n", " 1. Load the GO ontology\n", " 2. Create (or load) custom annotations in id2go format\n", " 3. Build a `TermCounts` object\n", " 4. Retrieve per-term information content (IC)\n", " 5. Compare IC across two different custom annotation sets\n", " 6. Calculate semantic similarity using custom annotations" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 1. Load the GO ontology\n", "\n", "Download `go-basic.obo` if it is not already present." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from goatools.base import get_godag\n", "\n", "godag = get_godag(\"go-basic.obo\", optional_attrs=[\"relationship\"])\n", "print(f\"{len(godag):,} GO terms loaded\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 2. Create custom annotations in id2go format\n", "\n", "### 2a. Define annotations in memory\n", "\n", "You can build a plain Python `dict` mapping gene IDs to sets of GO IDs. This is\n", "the **id2go** format in memory and is accepted directly by `TermCounts`.\n", "\n", "Here we simulate two annotation tools that annotated the same gene set differently:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# --- Annotation set produced by Tool 1 (uses more specific terms) ---\n", "annots_tool1 = {\n", " \"GeneA\": {\"GO:0048870\", # cell motility (BP)\n", " \"GO:0006928\"}, # movement of cell or subcellular component (BP)\n", " \"GeneB\": {\"GO:0048870\", # cell motility (BP)\n", " \"GO:0007049\"}, # cell cycle (BP)\n", " \"GeneC\": {\"GO:0007049\", # cell cycle (BP)\n", " \"GO:0008283\"}, # cell population proliferation (BP)\n", " \"GeneD\": {\"GO:0008283\", # cell population proliferation (BP)\n", " \"GO:0006928\"}, # movement of cell or subcellular component (BP)\n", "}\n", "\n", "# --- Annotation set produced by Tool 2 (uses broader terms) ---\n", "annots_tool2 = {\n", " \"GeneA\": {\"GO:0008150\"}, # biological_process (root BP \u2013 very broad)\n", " \"GeneB\": {\"GO:0008150\", # biological_process (root)\n", " \"GO:0007049\"}, # cell cycle (BP)\n", " \"GeneC\": {\"GO:0007049\", # cell cycle (BP)\n", " \"GO:0008283\"}, # cell population proliferation (BP)\n", " \"GeneD\": {\"GO:0008150\"}, # biological_process (root)\n", "}\n", "\n", "print(\"Tool 1 annotations:\", {g: sorted(gos) for g, gos in annots_tool1.items()})\n", "print()\n", "print(\"Tool 2 annotations:\", {g: sorted(gos) for g, gos in annots_tool2.items()})" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2b. Write and read an id2go annotation file\n", "\n", "If your annotations already live in a file, skip to **2c**.\n", "`IdToGosReader.wr_id2gos()` writes the in-memory dict to the id2go format:\n", "\n", "```\n", "GeneA GO:0006928;GO:0048870\n", "GeneB GO:0007049;GO:0048870\n", "...\n", "```" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from goatools.anno.idtogos_reader import IdToGosReader\n", "\n", "# Write Tool 1 annotations to file\n", "fout_tool1 = \"custom_tool1.txt\"\n", "IdToGosReader.wr_id2gos(fout_tool1, annots_tool1)\n", "\n", "fout_tool2 = \"custom_tool2.txt\"\n", "IdToGosReader.wr_id2gos(fout_tool2, annots_tool2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 2c. Read an id2go annotation file\n", "\n", "`IdToGosReader` reads an id2go file and returns namedtuples that can be\n", "converted to a `gene -> set(GO IDs)` dict via `get_id2gos()`.\n", "\n", "Setting `godag=godag` ensures that:\n", "* obsolete GO IDs are flagged\n", "* namespace information (BP / MF / CC) is attached to every annotation" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "reader1 = IdToGosReader(fout_tool1, godag=godag)\n", "id2gos_tool1 = reader1.get_id2gos() # dict[gene_id -> frozenset(GO_IDs)]\n", "\n", "reader2 = IdToGosReader(fout_tool2, godag=godag)\n", "id2gos_tool2 = reader2.get_id2gos()\n", "\n", "print(\"Genes loaded from Tool 1 file:\", sorted(id2gos_tool1.keys()))\n", "print(\"Genes loaded from Tool 2 file:\", sorted(id2gos_tool2.keys()))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Build TermCounts objects\n", "\n", "`TermCounts` propagates annotations up the ontology hierarchy (every gene\n", "annotated to a child term is also counted for all ancestor terms) and\n", "pre-computes term frequencies and information contents." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from goatools.semantic import TermCounts\n", "\n", "# You can pass either the id2gos dict read from file (id2gos_tool1)\n", "# or the original in-memory dict (annots_tool1) \u2013 both work identically.\n", "tcnt1 = TermCounts(godag, id2gos_tool1)\n", "tcnt2 = TermCounts(godag, id2gos_tool2)\n", "\n", "print(\"Tool 1 - aspect counts:\", tcnt1.aspect_counts)\n", "print(\"Tool 2 - aspect counts:\", tcnt2.aspect_counts)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Retrieve information content per GO term\n", "\n", "`get_info_content(go_id, termcounts)` returns **IC = \u2212log(frequency)** where\n", "*frequency* is the proportion of genes (in this annotation set) annotated to\n", "the term or any of its descendants.\n", "\n", "* **High IC** \u2192 rare / specific term\n", "* **IC = 0** \u2192 universal root term (every gene is annotated to it)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from goatools.semantic import get_info_content\n", "\n", "# Inspect a representative set of GO IDs\n", "go_ids_of_interest = [\n", " \"GO:0008150\", # biological_process (root)\n", " \"GO:0007049\", # cell cycle\n", " \"GO:0048870\", # cell motility\n", " \"GO:0008283\", # cell population proliferation\n", " \"GO:0006928\", # movement of cell or subcellular component\n", "]\n", "\n", "print(f\"{'GO ID':<14} {'GO name':<45} {'IC tool1':>8} {'IC tool2':>8}\")\n", "print(\"-\" * 80)\n", "for go_id in go_ids_of_interest:\n", " if go_id not in godag:\n", " continue\n", " name = godag[go_id].name\n", " ic1 = get_info_content(go_id, tcnt1)\n", " ic2 = get_info_content(go_id, tcnt2)\n", " print(f\"{go_id:<14} {name:<45} {ic1:8.4f} {ic2:8.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 5. Compare IC between two annotation sets\n", "\n", "IC values reflect **how informative each GO term is in a given annotation set**.\n", "The same term can have different IC in different annotation sets:\n", "\n", "* **Lower IC** in Tool 2 for many terms \u2192 Tool 2 tends to use broad/root terms\n", " that are common across many genes (low information).\n", "* **Higher IC** in Tool 1 for the same terms \u2192 Tool 1 uses more specific\n", " annotations (higher information)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# Collect all GO IDs mentioned by either tool (excl. root term)\n", "all_goids = (\n", " set.union(*id2gos_tool1.values()) | set.union(*id2gos_tool2.values())\n", ") - {\"GO:0008150\"} # skip root\n", "\n", "print(f\"{'GO ID':<14} {'GO name':<45} {'IC tool1':>8} {'IC tool2':>8} {'\u0394 IC':>8}\")\n", "print(\"-\" * 88)\n", "for go_id in sorted(all_goids):\n", " if go_id not in godag:\n", " continue\n", " name = godag[go_id].name\n", " ic1 = get_info_content(go_id, tcnt1)\n", " ic2 = get_info_content(go_id, tcnt2)\n", " delta = ic1 - ic2\n", " print(f\"{go_id:<14} {name:<45} {ic1:8.4f} {ic2:8.4f} {delta:+8.4f}\")\n", "\n", "print()\n", "print(\"Positive \u0394 IC \u2192 tool1 treats that term as MORE informative (more specific)\")\n", "print(\"Negative \u0394 IC \u2192 tool2 treats that term as MORE informative\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### 5b. GO terms unique to each annotation set\n", "\n", "Terms present in one annotation set but not the other can also be compared." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "goids_only_tool1 = set.union(*id2gos_tool1.values()) - set.union(*id2gos_tool2.values())\n", "goids_only_tool2 = set.union(*id2gos_tool2.values()) - set.union(*id2gos_tool1.values())\n", "\n", "print(\"GO terms unique to Tool 1:\")\n", "for go_id in sorted(goids_only_tool1):\n", " ic1 = get_info_content(go_id, tcnt1)\n", " print(f\" {go_id} {godag[go_id].name:<45} IC={ic1:.4f}\")\n", "\n", "print()\n", "print(\"GO terms unique to Tool 2:\")\n", "for go_id in sorted(goids_only_tool2):\n", " ic2 = get_info_content(go_id, tcnt2)\n", " print(f\" {go_id} {godag[go_id].name:<45} IC={ic2:.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 6. Semantic similarity using custom annotations\n", "\n", "Once you have a `TermCounts` object you can compute pairwise **semantic\n", "similarity** between any pair of GO terms using the annotation set of your\n", "choice. Three measures are available:\n", "\n", "| Measure | Function | Description |\n", "|---------|----------|-------------|\n", "| Resnik | `resnik_sim` | IC of the most specific common ancestor (MSCA) |\n", "| Lin | `lin_sim` | Normalised Resnik: `2\u00b7IC(MSCA) / (IC(t1) + IC(t2))` |\n", "| Schlicker | `schlicker_sim` | Lin weighted by `(1 \u2212 freq(MSCA))` |" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from goatools.semantic import resnik_sim, lin_sim, schlicker_sim\n", "\n", "# Compare two GO terms using Tool 1's annotation context\n", "go_a = \"GO:0048870\" # cell motility\n", "go_b = \"GO:0007049\" # cell cycle\n", "\n", "print(f\"Comparing {go_a} ({godag[go_a].name})\")\n", "print(f\" vs {go_b} ({godag[go_b].name})\")\n", "print()\n", "\n", "r1 = resnik_sim(go_a, go_b, godag, tcnt1)\n", "l1 = lin_sim(go_a, go_b, godag, tcnt1)\n", "s1 = schlicker_sim(go_a, go_b, godag, tcnt1)\n", "\n", "r2 = resnik_sim(go_a, go_b, godag, tcnt2)\n", "l2 = lin_sim(go_a, go_b, godag, tcnt2)\n", "s2 = schlicker_sim(go_a, go_b, godag, tcnt2)\n", "\n", "print(f\"{'Measure':<12} {'Tool 1':>10} {'Tool 2':>10}\")\n", "print(\"-\" * 36)\n", "print(f\"{'Resnik':<12} {r1 if r1 is not None else 'N/A':>10.4f} {r2 if r2 is not None else 'N/A':>10.4f}\")\n", "print(f\"{'Lin':<12} {l1 if l1 is not None else 'N/A':>10.4f} {l2 if l2 is not None else 'N/A':>10.4f}\")\n", "print(f\"{'Schlicker':<12} {s1 if s1 is not None else 'N/A':>10.4f} {s2 if s2 is not None else 'N/A':>10.4f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## 7. Clean up temporary annotation files" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import os\n", "for f in (fout_tool1, fout_tool2):\n", " if os.path.exists(f):\n", " os.remove(f)\n", " print(f\"Removed {f}\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "---\n", "\n", "## Summary\n", "\n", "| Step | Code |\n", "|------|------|\n", "| Load ontology | `get_godag(\"go-basic.obo\")` |\n", "| In-memory annotations | Plain `dict[gene_id -> set(GO_IDs)]` |\n", "| Write id2go file | `IdToGosReader.wr_id2gos(filename, id2gos)` |\n", "| Read id2go file | `IdToGosReader(filename, godag=godag).get_id2gos()` |\n", "| Compute term counts | `TermCounts(godag, id2gos)` |\n", "| Information content | `get_info_content(go_id, termcounts)` |\n", "| Resnik similarity | `resnik_sim(go1, go2, godag, termcounts)` |\n", "| Lin similarity | `lin_sim(go1, go2, godag, termcounts)` |\n", "| Schlicker similarity | `schlicker_sim(go1, go2, godag, termcounts)` |\n", "\n", "Copyright (C) 2024, DV Klopfenstein et al. All rights reserved." ] } ] }