Information Content for Custom GO Annotations (id2go format)

This notebook demonstrates how to calculate the information content (IC) of GO terms for any species using a custom annotation file in id2go format, i.e. a tab-separated file where:

• column 1 = gene identifier

• column 2 = semicolon-separated list of GO terms

This is useful when you want to:

• Analyse non-model organisms not covered by standard GOA databases

• Compare annotation quality produced by different tools

• Identify which GO terms are shared between two annotation sets and how specific they are

Steps:

1. Load the GO ontology

2. Create (or load) custom annotations in id2go format

3. Build a TermCounts object

4. Retrieve per-term information content (IC)

5. Compare IC across two different custom annotation sets

6. Calculate semantic similarity using custom annotations

1. Load the GO ontology

Download go-basic.obo if it is not already present.

from goatools.base import get_godag

godag = get_godag("go-basic.obo", optional_attrs=["relationship"])
print(f"{len(godag):,} GO terms loaded")

2. Create custom annotations in id2go format

2a. Define annotations in memory

You can build a plain Python dict mapping gene IDs to sets of GO IDs. This is the id2go format in memory and is accepted directly by TermCounts.

Here we simulate two annotation tools that annotated the same gene set differently:

# --- Annotation set produced by Tool 1 (uses more specific terms) ---
annots_tool1 = {
    "GeneA": {"GO:0048870",  # cell motility (BP)
              "GO:0006928"},  # movement of cell or subcellular component (BP)
    "GeneB": {"GO:0048870",  # cell motility (BP)
              "GO:0007049"},  # cell cycle (BP)
    "GeneC": {"GO:0007049",  # cell cycle (BP)
              "GO:0008283"},  # cell population proliferation (BP)
    "GeneD": {"GO:0008283",  # cell population proliferation (BP)
              "GO:0006928"},  # movement of cell or subcellular component (BP)
}

# --- Annotation set produced by Tool 2 (uses broader terms) ---
annots_tool2 = {
    "GeneA": {"GO:0008150"},  # biological_process (root BP – very broad)
    "GeneB": {"GO:0008150",   # biological_process (root)
              "GO:0007049"},  # cell cycle (BP)
    "GeneC": {"GO:0007049",   # cell cycle (BP)
              "GO:0008283"},  # cell population proliferation (BP)
    "GeneD": {"GO:0008150"},  # biological_process (root)
}

print("Tool 1 annotations:", {g: sorted(gos) for g, gos in annots_tool1.items()})
print()
print("Tool 2 annotations:", {g: sorted(gos) for g, gos in annots_tool2.items()})

2b. Write and read an id2go annotation file

If your annotations already live in a file, skip to 2c. IdToGosReader.wr_id2gos() writes the in-memory dict to the id2go format:

GeneA   GO:0006928;GO:0048870
GeneB   GO:0007049;GO:0048870
...

from goatools.anno.idtogos_reader import IdToGosReader

# Write Tool 1 annotations to file
fout_tool1 = "custom_tool1.txt"
IdToGosReader.wr_id2gos(fout_tool1, annots_tool1)

fout_tool2 = "custom_tool2.txt"
IdToGosReader.wr_id2gos(fout_tool2, annots_tool2)

2c. Read an id2go annotation file

IdToGosReader reads an id2go file and returns namedtuples that can be converted to a gene -> set(GO IDs) dict via get_id2gos().

Setting godag=godag ensures that:

• obsolete GO IDs are flagged

• namespace information (BP / MF / CC) is attached to every annotation

reader1 = IdToGosReader(fout_tool1, godag=godag)
id2gos_tool1 = reader1.get_id2gos()   # dict[gene_id -> frozenset(GO_IDs)]

reader2 = IdToGosReader(fout_tool2, godag=godag)
id2gos_tool2 = reader2.get_id2gos()

print("Genes loaded from Tool 1 file:", sorted(id2gos_tool1.keys()))
print("Genes loaded from Tool 2 file:", sorted(id2gos_tool2.keys()))

3. Build TermCounts objects

TermCounts propagates annotations up the ontology hierarchy (every gene annotated to a child term is also counted for all ancestor terms) and pre-computes term frequencies and information contents.

from goatools.semantic import TermCounts

# You can pass either the id2gos dict read from file (id2gos_tool1)
# or the original in-memory dict (annots_tool1) – both work identically.
tcnt1 = TermCounts(godag, id2gos_tool1)
tcnt2 = TermCounts(godag, id2gos_tool2)

print("Tool 1 - aspect counts:", tcnt1.aspect_counts)
print("Tool 2 - aspect counts:", tcnt2.aspect_counts)

4. Retrieve information content per GO term

get_info_content(go_id, termcounts) returns IC = −log(frequency) where frequency is the proportion of genes (in this annotation set) annotated to the term or any of its descendants.

• High IC → rare / specific term

• IC = 0 → universal root term (every gene is annotated to it)

from goatools.semantic import get_info_content

# Inspect a representative set of GO IDs
go_ids_of_interest = [
    "GO:0008150",  # biological_process (root)
    "GO:0007049",  # cell cycle
    "GO:0048870",  # cell motility
    "GO:0008283",  # cell population proliferation
    "GO:0006928",  # movement of cell or subcellular component
]

print(f"{'GO ID':<14} {'GO name':<45} {'IC tool1':>8}  {'IC tool2':>8}")
print("-" * 80)
for go_id in go_ids_of_interest:
    if go_id not in godag:
        continue
    name = godag[go_id].name
    ic1 = get_info_content(go_id, tcnt1)
    ic2 = get_info_content(go_id, tcnt2)
    print(f"{go_id:<14} {name:<45} {ic1:8.4f}  {ic2:8.4f}")

5. Compare IC between two annotation sets

IC values reflect how informative each GO term is in a given annotation set. The same term can have different IC in different annotation sets:

• Lower IC in Tool 2 for many terms → Tool 2 tends to use broad/root terms that are common across many genes (low information).

• Higher IC in Tool 1 for the same terms → Tool 1 uses more specific annotations (higher information).

# Collect all GO IDs mentioned by either tool (excl. root term)
all_goids = (
    set.union(*id2gos_tool1.values()) | set.union(*id2gos_tool2.values())
) - {"GO:0008150"}  # skip root

print(f"{'GO ID':<14} {'GO name':<45} {'IC tool1':>8}  {'IC tool2':>8}  {'Δ IC':>8}")
print("-" * 88)
for go_id in sorted(all_goids):
    if go_id not in godag:
        continue
    name = godag[go_id].name
    ic1 = get_info_content(go_id, tcnt1)
    ic2 = get_info_content(go_id, tcnt2)
    delta = ic1 - ic2
    print(f"{go_id:<14} {name:<45} {ic1:8.4f}  {ic2:8.4f}  {delta:+8.4f}")

print()
print("Positive Δ IC → tool1 treats that term as MORE informative (more specific)")
print("Negative Δ IC → tool2 treats that term as MORE informative")

5b. GO terms unique to each annotation set

Terms present in one annotation set but not the other can also be compared.

goids_only_tool1 = set.union(*id2gos_tool1.values()) - set.union(*id2gos_tool2.values())
goids_only_tool2 = set.union(*id2gos_tool2.values()) - set.union(*id2gos_tool1.values())

print("GO terms unique to Tool 1:")
for go_id in sorted(goids_only_tool1):
    ic1 = get_info_content(go_id, tcnt1)
    print(f"  {go_id}  {godag[go_id].name:<45}  IC={ic1:.4f}")

print()
print("GO terms unique to Tool 2:")
for go_id in sorted(goids_only_tool2):
    ic2 = get_info_content(go_id, tcnt2)
    print(f"  {go_id}  {godag[go_id].name:<45}  IC={ic2:.4f}")

6. Semantic similarity using custom annotations

Once you have a TermCounts object you can compute pairwise semantic similarity between any pair of GO terms using the annotation set of your choice. Three measures are available:

Measure	Function	Description
Resnik	resnik_sim	IC of the most specific common ancestor (MSCA)
Lin	lin_sim	Normalised Resnik: 2·IC(MSCA) / (IC(t1) + IC(t2))
Schlicker	schlicker_sim	Lin weighted by (1 − freq(MSCA))

from goatools.semantic import resnik_sim, lin_sim, schlicker_sim

# Compare two GO terms using Tool 1's annotation context
go_a = "GO:0048870"  # cell motility
go_b = "GO:0007049"  # cell cycle

print(f"Comparing {go_a} ({godag[go_a].name})")
print(f"      vs  {go_b} ({godag[go_b].name})")
print()

r1 = resnik_sim(go_a, go_b, godag, tcnt1)
l1 = lin_sim(go_a, go_b, godag, tcnt1)
s1 = schlicker_sim(go_a, go_b, godag, tcnt1)

r2 = resnik_sim(go_a, go_b, godag, tcnt2)
l2 = lin_sim(go_a, go_b, godag, tcnt2)
s2 = schlicker_sim(go_a, go_b, godag, tcnt2)

print(f"{'Measure':<12} {'Tool 1':>10}  {'Tool 2':>10}")
print("-" * 36)
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}")
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}")
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}")

7. Clean up temporary annotation files

import os
for f in (fout_tool1, fout_tool2):
    if os.path.exists(f):
        os.remove(f)
        print(f"Removed {f}")

---

Summary

Step	Code
Load ontology	get_godag("go-basic.obo")
In-memory annotations	Plain dict[gene_id -> set(GO_IDs)]
Write id2go file	IdToGosReader.wr_id2gos(filename, id2gos)
Read id2go file	IdToGosReader(filename, godag=godag).get_id2gos()
Compute term counts	TermCounts(godag, id2gos)
Information content	get_info_content(go_id, termcounts)
Resnik similarity	resnik_sim(go1, go2, godag, termcounts)
Lin similarity	lin_sim(go1, go2, godag, termcounts)
Schlicker similarity	schlicker_sim(go1, go2, godag, termcounts)

Copyright (C) 2024, DV Klopfenstein et al. All rights reserved.