This is the answered version of the GO Enrichment Analysis practical. Every ❓ question is followed by a worked model answer (green box). All the interactive steps still work, so you can keep clustering and re-running the analysis while you read, and the self-test quiz at the end is unchanged.
This is the same dataset you clustered earlier: 500 genes measured in Arabidopsis
plants that were treated in three ways and sampled at three time points. Clustering put
genes that behave in the same way into the same group.
That is a useful start, but look at what a cluster actually gives you: a list of codes like
AT1G29930.
A list of names does not tell you what those genes do, and that is what you need if
you want to say anything about how the plant defends itself. Closing that gap is what the rest
of this practical is about.
Distance:
Linkage:
Each row is one gene, each column one sample. Blue means the gene is expressed below its own
average, red above it. The coloured bar on the left shows which cluster each gene ended up in,
labelled C1, C2 and so on. The cluster you have selected is outlined in white.
low → high expression
Pick the cluster you want to work with
Everything in steps 2 to 4 follows this choice, so you can come back and change it at any time.
AVR (avirulent bacteria) VIR (virulent bacteria) MOCK (no bacteria)
❓ Question
Which cluster shows an expression pattern that looks promising if you want to study
resistance to this biotic stress? Say what it is about the pattern that convinces you.
✅ Answer
The promising cluster. At the default settings (k = 4, correlation, average) that is C2, the small cluster of 18 genes. What makes it convincing is not that it goes up, but when and where it goes up: it rises in AVR at 6 and 12 hours, much less in VIR, and not at all in MOCK. Because the MOCK plants went through the same handling but received no bacteria, a pattern that is absent in MOCK cannot be explained by the procedure. And because it is strongest in AVR, the treatment in which the plant successfully defends itself, it is a candidate for the resistance response specifically. Compare that with C1, which is much larger and rises sharply at 1 hour, but rises just as much in MOCK, so it cannot be about the bacteria. If you changed the settings your cluster numbers will differ; the pattern to look for is the one that separates infected plants from MOCK.
2 Looking up what a gene does
Biologists have spent decades writing down what individual genes do. The problem is that
everybody writes it down differently: "helps against bacteria", "bacterial defence",
"responds to Pseudomonas". A computer cannot tell that these mean the same thing.
The Gene Ontology (GO) solves this. It is a fixed dictionary of biological
processes that everyone agrees to use, so a result from one lab can be compared with a result
from another. Every entry is called a GO term and has an id and a name, such
as GO:0009617 response to bacterium.
Two things about GO matter for what follows. First, the terms are ordered from
general to specific: response to stress sits above response to
bacterium, which sits above still more specific terms. If a gene is labelled with a
specific term, it automatically counts for all the more general terms above it as well.
Second, a gene usually carries several terms at once, because a gene can be
involved in more than one process.
Look one gene up
The terms of this gene, drawn from general (top) to specific (bottom).
Green terms are the ones the gene is directly
labelled with; grey terms are the more general ones it inherits automatically.
❓ Questions
What functional annotations does your gene have? Is each GO term equally specific?
What is a disadvantage of looking only at the annotations of a single gene?
Do you think all the GO terms you find for your gene make sense?
💡 Hint if every gene you try looks convincing
Try a few more genes from the same cluster, including one or two that you pick at random
rather than choosing yourself. Also look at how many terms each gene has, and check whether
the terms were measured in the lab or predicted by a computer.
✅ Answer
What your gene does, and specificity. The exact terms depend on the gene you picked, but the shape of the answer is always the same: a handful of terms, and they are not equally specific. You will see very broad ones such as response to stress, cellular process or metabolic process alongside narrow ones such as photosynthesis, light harvesting or response to high light intensity. The graph shows this directly: general terms sit at the top, specific ones at the bottom, and a gene inherits every term above the one it is labelled with. This matters because a broad term is close to useless on its own. Saying a gene is involved in "a metabolic process" barely narrows anything down, since a large share of all genes carry that label.
The disadvantage of one gene at a time. Several, and they compound. It does not scale: nobody is going to look up 214 genes by hand. It is biased by you: because you choose which gene to inspect, you can unconsciously pick the one that fits the story you already had, and in a list of 200 genes there is almost always one that fits any story. It tells you nothing about the group: one photosynthesis gene in a cluster is not evidence that the cluster is about photosynthesis. And the annotations themselves are incomplete and uneven, 109 of these 500 genes have no biological-process annotation at all, and heavily studied genes carry far more terms than obscure ones simply because more people have looked at them.
Do all the terms make sense? Usually not, and that is the point. You will find terms such as response to cold, response to light intensity or response to salt on genes in a bacterial infection experiment. That is not an error in the database. A GO term records what a gene has been observed to do in some experiment somewhere, not what it is doing in yours. Many stress responses share the same machinery, so the same gene turns up in cold, light and pathogen studies alike. Also check the "how we know" column: some annotations were measured in the lab, others were predicted by a computer from similarity to another gene, and those two deserve different amounts of trust.
3 Enrichment: doing this for the whole cluster at once
Looking up genes one at a time does not scale, and it is easy to fool yourself: if you go
hunting through 200 genes you will always find one that fits the story you had in mind.
GO enrichment does the job for the whole cluster at once, and it does it by
counting.
The idea is simple. Take one GO term. Count how many genes in your cluster carry it. Then
count how many genes in the whole dataset carry it. If the cluster has clearly more than its
fair share, that term is enriched, and the process it describes is probably
what this group of genes is doing. Repeat for every term.
How do we decide what counts as "clearly more"?
The page works out how likely that count would be if you had simply grabbed the same number
of genes at random out of the dataset. That likelihood is the p-value: a
small p-value means random picking would almost never give you this many. The calculation is
a standard one called Fisher's exact test, and it is exactly what the
goatools package in the notebook does for you.
Why the numbers below are "corrected"
There is a catch. We do not test one term, we test about 1,300 of them. If you run 1,300
tests, some will look impressive purely by chance, in the same way that rolling a die often
enough will eventually give you five sixes in a row. The correction adjusts every p-value to
account for how many tests were run. Always read the corrected column.
You can switch the colouring below to see how much difference this makes.
Enrichment results for your cluster
Colour by:
Show:
p < 0.005 p < 0.01 p < 0.05 p > 0.05, not significant
The enriched terms in their general-to-specific order, coloured by p-value. Click a term to
see which of your cluster's genes carry it.
❓ Questions
Do the enriched GO terms of the cluster make sense, looking at its expression profile in
step 1?
What is the goal of GO enrichment analysis? Try to say it in one sentence.
Find a cluster that returns no enriched terms at all. Is that a failure of the method,
or is it telling you something?
How would you use the enriched GO terms of your cluster to inform a follow-up
experiment? Be concrete: what would you measure, on which plants, and at what moment?
✅ Answer
Do the terms match the profile? At the default settings they line up well, which is the reassuring case. C1 rises at 1 hour and falls afterwards, and comes back enriched for photosynthesis, photosynthesis light reaction and response to high light intensity: processes that track light and time of day rather than infection. C2 rises in AVR at 6 and 12 hours and comes back enriched for response to hypoxia and regulation of programmed cell death. Programmed cell death is exactly what a plant does when it recognises a pathogen, it kills its own cells around the infection site to stop the bacteria spreading, so the profile and the terms tell the same story. C4 rises late in every treatment and is enriched for translation and protein metabolic process, consistent with cells stepping up protein production. Where the profile and the terms agree you have something worth testing. Where they disagree, be suspicious of the clustering, of the annotations, or of both.
The goal, in one sentence. To find out which biological processes are over-represented in a group of genes, so that a list of gene names turns into a statement about biology that you can act on.
A cluster with nothing enriched. At the default settings C3 is one. This is not a failure of the method, and it is not a bug. It means these genes move together but share no known process, and there are several honest reasons for that: the cluster may be too small for any count to reach significance, most of its genes may be unannotated, the genes may be co-regulated for a reason nobody has written down yet, or the cluster may be a technical artefact rather than biology, C3 is driven almost entirely by a single sample being unusual. An empty result is genuinely useful information: it is a warning not to build a story on that cluster, and it is exactly the kind of result you would never notice if you had gone looking through the gene list by hand until you found something that fitted.
Designing the follow-up. The useful move is to stop measuring RNA and start measuring the process the terms point at. If your cluster is enriched for programmed cell death and hypoxia in the AVR treatment, then measure cell death itself, for example by staining leaves or measuring ion leakage, in AVR, VIR and MOCK plants at 6 and 12 hours. Take two or three of the genes carrying those terms, obtain knockout or overexpression lines for them, repeat the infection, and see whether the plant's resistance actually changes, that is the step that turns a correlation into a cause. Sample more densely between 1 and 6 hours, since that is the window in which the response appears to start and you currently have no measurements there. And since the goal is resilient crops rather than Arabidopsis, look up the orthologues of those genes in a crop species and check whether they behave the same way, which is what the next practical does. The general pattern: the enriched term tells you which assay to run, the genes carrying it tell you which plant lines to build, and the expression profile tells you when to sample.
4 Before you trust the result
You have just sketched a follow-up experiment. Before anyone spends a season on it, there is
one check worth doing, and this step is about that check.
An enrichment result is a lead, not a conclusion. It tells you that genes already known to be
involved in some process changed together. It does not tell you that the process caused
anything. And, more awkwardly, it does not tell you that your treatment caused the
change: if something else happened to the plants at the same time as the infection, the
analysis will report that just as happily, and it will look just as convincing.
The genes that change at 1 hour
In the original study the plants were sampled 1, 6 and 12 hours after treatment. The button
below resets the clustering to its default settings and selects the cluster whose genes
respond most strongly at the 1 hour time point, so you can run the analysis on it.
AVR VIR MOCK
❓ Question
Can you use this methodology to explain what seems to cause the gene expression changes
in the plant at 1 hour after infection? Compare the MOCK line with the two infected lines
before you answer, and look at which processes come out enriched.
✅ Answer
What is happening at 1 hour. Look at cluster C1 carefully. It shoots up at 1 hour in AVR, in VIR and in MOCK. The MOCK plants received buffer with no bacteria in it at all, so whatever is driving this cannot be the pathogen. It has to be the part of the procedure that all three groups share: the plants were handled, infiltrated with liquid, wounded slightly, and sampled at one particular moment of the day. The enriched terms confirm it, they are photosynthesis, light reaction, response to high light intensity and response to cold, all things that change when a plant is moved and disturbed, and none of them anything to do with bacteria.
Why that matters more than it looks. This is the largest cluster in the dataset, so a large part of the strongest signal in the whole experiment is the experiment happening to the plant rather than the infection. GO enrichment did not explain the infection here; it exposed a confounder. Two practical consequences follow. First, always compare infected against MOCK at the same time point, never against an earlier time point, because time itself changes the plant. Second, 1 hour is the worst place in this dataset to go looking for defence genes, since the handling response drowns them out. If this feels familiar, it is the same lesson as the brightness confounder in the deep-learning practical: the method faithfully found the strongest pattern in the data, and the strongest pattern was not the one anybody was interested in.
If you have time left
Go back to step 1, change the number of clusters or the linkage, and run the enrichment on
other clusters. Do the same biological processes keep turning up, or does each set of
settings tell a different story? That is worth knowing before you trust any single result.