This is the answered version of the Transcriptomics & Clustering practical. Every ❓ question is followed by a worked model answer (green box). All the interactive steps still work, so you can keep re-clustering the data and exploring the network while you read, and the self-test quiz at the end is unchanged.
The researchers took Arabidopsis plants and applied three treatments. Each treatment was sampled at three time points (1, 6 and 12 hours), giving 9 samples in total. For each sample we have the expression of the same 500 genes. The numbers are z-scored per gene: for every gene the nine values were rescaled to have mean 0, so a positive number means "higher than that gene's average" and a negative number means "lower than average". This lets us compare the shape of a gene's response across samples rather than its absolute level.
A small slice of the data: rows are genes, columns are the 9 samples, numbers are z-scored expression.
❓ Questions
Look at the paper. What are the three treatments, and what is the goal of comparing them?
Why do you think the plants were sampled at several time points instead of just one?
💡 Hint on where to find the treatments
In the paper, scroll down to Results and look at the section called Dataset. The very first paragraph there lists what was infiltrated into the leaves for each of the three treatments. While you read it, ask yourself which treatment received no bacteria at all, and what the difference is between the two bacterial strains, because that difference is the reason the experiment was designed with three treatments rather than two.
✅ Answer
The three treatments and why compare them.MOCK is the control: the plant is infiltrated with liquid but receives no bacteria. AVR is infection with an avirulent strain, one the plant recognises and successfully fights off. VIR is infection with a virulent strain, one that evades the plant's defences and causes disease. Comparing either infection against MOCK tells you what the bacteria did rather than what the handling did, and comparing AVR against VIR tells you what a successful defence looks like.
Why several time points. A defence response is a sequence of events, not a state: the plant has to detect the bacteria, pass the signal on, switch on transcription factors and only then make defence proteins, and those stages are hours apart. A single snapshot catches one frame of that with no way of knowing which. Several time points show the order in which genes come on, and protect you from sampling too early or too late.
2 Reading a clustered heatmap
A heatmap shows the whole 500×9 table at once as colour: red means a gene is expressed above its average in that sample, blue means below. On its own that is just a wall of colour. The power comes from clustering: we reorder the rows so that genes with similar response patterns sit next to each other, and reorder the columns so that similar samples sit next to each other. The tree above the columns (a dendrogram) shows which samples were merged together, and how early.
low → high expression (z-score) AVR VIR MOCK
❓ Questions
What do the rows and the columns represent, and what do the colours mean?
Is the avirulent (AVR) sample more similar to the virulent (VIR) sample or to the mock sample? Look at the tree above the heatmap.
Does the treatment or the time point play the bigger role in how the samples group together?
Explain the clustering in your own words: what does it tell you about the relationships between the samples?
💡 Hint on reading the dendrogram
Do not answer this one from the names of the treatments. To work out how far apart two samples are, start at each of them and follow the tree upwards until the two paths meet. The height at which they join is the distance the algorithm found between them, so the lower they meet, the more alike they are. Go through the tree in that order: which two samples get merged together first, which sample joins them next, and which one is left until last? Read that off the tree using only the data. Once you have that order, go back to the biology and ask which explanation fits what you found.
✅ Answer
Rows, columns and colours. Each row is one of the 500 genes, each column is one of the 9 samples (one treatment at one time point), and the colour is that gene's z-scored expression: red means above that gene's own average, blue below, pale about average. Because the scores are computed per gene, colours are only comparable along a row, so a red cell does not mean the gene is highly expressed in absolute terms.
Is AVR closer to VIR or to MOCK? To VIR. The dendrogram shows it twice. At 1 hour, AVR_1 and VIR_1 merge with each other first and MOCK_1 joins only afterwards. Among the later samples, the AVR pair and the VIR pair join each other before the MOCK pair joins them. Once you stop reading the names this makes sense: avirulent and virulent sound like opposites, but from the plant's point of view both mean bacteria are present, and MOCK is the odd one out because nothing is attacking it.
Treatment or time point? Both, at different levels. The very first split separates the three 1-hour samples from all six later ones, so at the coarsest level time dominates. Inside the later group the structure is entirely by treatment, with AVR pairing with AVR and so on. Note that all three 1-hour samples cluster together, MOCK included, so whatever makes them distinctive cannot be the bacteria: it is the handling, wounding and infiltration every plant went through, plus the time of day.
The clustering in your own words. Something like: the algorithm knew nothing about treatments or time points, it only saw 9 columns of numbers, and it still recovered the design of the experiment. That doubles as a sanity check, because if matching conditions had not grouped together you would suspect a mislabelled sample or a batch effect long before interpreting any individual gene.
3 How hierarchical clustering works, and how many clusters?
Hierarchical clustering is unsupervised: nobody tells the algorithm which genes belong together, it discovers structure on its own. It is beautifully simple. Every gene starts as its own tiny cluster. Then, over and over, the algorithm finds the two closest clusters and merges them, until everything is joined into one tree. Two choices decide the outcome:
Distance: how we measure "closeness" between two genes. Correlation distance counts two genes as close if their expression rises and falls together across the samples, regardless of absolute level. Euclidean distance instead compares the raw values point by point.
Linkage: once clusters contain several genes, how we measure the distance between two groups: the average distance between their members, the distance between their closest members (single), or their farthest members (complete).
The figure below draws this process directly onto the heatmap. On the left is the tree; to actually get clusters we cut it at some height (the dashed line), and every branch below the cut becomes one cluster. The coloured bar between the tree and the heatmap shows which cluster each gene row belongs to, so you can see the coloured blocks in the heatmap line up with the clusters. Slide the cut, and switch the distance and linkage, and watch genes move between clusters and the heatmap re-sort itself.
Distance:
Linkage:
low → high expression (z-score) dashed line = the cut
Left: the gene tree, coloured by cluster. Middle bar: each gene's cluster. Right: the heatmap, its rows sorted by the tree so clusters appear as coloured blocks.
Number of genes in each cluster.
Which genes are in cluster:
Notice how the number of clusters, and even which genes end up together, shift as you move the cut and switch the distance or linkage. There is rarely one obviously correct answer, and choosing a clustering is a judgement call that depends on the biological question you are asking.
❓ Questions
How many clusters do you think is most appropriate for this dataset? Why? Is there a single right answer?
Switch between the linkage and distance options. What differences do you observe, and can you explain them?
Which cluster(s) of genes would be most interesting to study further, and how would you decide?
Try to link a cluster back to the underlying biology. Look at the coloured block a cluster forms in the heatmap: in which treatments and at which time points are those genes high, and in which are they low? What function do you think the genes in that cluster might be carrying out for the plant?
✅ Answer
How many clusters, and is there a right answer? No, and that is the point. Four is defensible at the default settings, but so are other choices, and the decision is a judgement rather than a calculation. The real test is whether the groups you get are interpretable and whether your conclusions survive a change of settings: if a story only appears at one particular value of k, it is not a story. Decide what you would pick and why, and bring that to class so we can compare what different people chose.
What changes when you switch the settings.Linkage matters enormously. Average and complete both give sensibly sized groups, while single collapses: at k = 4 it gives one cluster of 497 genes and three single-gene clusters. That is the classic failure called chaining, because single linkage merges two clusters as soon as their closest members are close, so one gene sitting between two groups is enough to fuse them. Distance barely matters here, because the values are already z-scored per gene, which makes Euclidean distance close to a rescaling of correlation distance. On raw expression the two would differ a great deal.
Which cluster to follow up, and how to decide. A cluster of genes that rises in AVR at 6 and 12 hours is the best candidate. The criterion is not that it changes a lot, but that it changes where the biology is interesting: it is high in the treatment where the plant successfully defends itself, much weaker in VIR, and flat in MOCK, so the handling cannot explain it.
Linking a cluster to the biology. Be honest about what a profile can do: it tells you when and in which treatment genes are active, so it lets you propose a function, not establish one. A cluster high at 1 hour in all three treatments cannot be about the bacteria, since MOCK had none, so those genes are more likely responding to handling and wounding or tracking the light. A cluster that rises in AVR but stays flat in MOCK is switched on only when the plant detects an attacker, so defence is the obvious guess. A cluster climbing late in every treatment looks more like general housekeeping. Each of those guesses came from where the block of colour sits rather than from the genes; the gene list under the figure lets you check a handful by hand, and doing it for hundreds at once is what the next practical is for.
If you have extra time, pick a cluster in the list above and click a few of its genes. Each one opens that gene in UniProt, where you can read what the protein is thought to do. Can you link them to a biological process involved in stress response, does the whole cluster tell a consistent story, and does it match what the Howard et al. paper reports?
4 A gene correlation network
Clustering is not the only way to look for structure. We can also build a network: draw one dot (node) per gene, and connect two genes with a line (edge) whenever their expression patterns are highly correlated across the nine samples. Genes that move together end up tied into tightly connected knots. Raise the threshold and only the very strongest relationships survive; lower it and the web fills in.
Colour nodes by:
Each dot is a gene; lines connect strongly correlated genes. Drag nodes around, scroll to zoom, and hover a node to see its name.
node colour & size: low high centrality
❓ Questions
What happens to the network as you change the correlation threshold?
What could it mean, biologically, that two genes are connected by an edge?
Think it through rather than looking it up: would you expect the hierarchical clusters from step 3 and the tightly connected knots here to correspond to each other? Why, or why not?
What can you learn from a network that you cannot learn from clustering, and what can you learn from clustering that you cannot learn from a network?
The colouring highlights "central" genes. If you worked with a plant biologist, which genes would you suggest studying next, and what experiment would you propose?
✅ Answer
What the threshold does. It sets how strong a correlation has to be before two genes are joined. At 0.85 almost everything is connected, roughly 15,000 edges among 489 genes, and the result is a hairball; at 0.95 about 2,900 edges among 378 genes leave distinct knots visible; at 0.99 only about 156 edges among 118 genes survive. The threshold is a choice you make, not something the data tells you, so check that the group you care about holds together across a range of values.
What an edge means biologically. That the two genes rise and fall together across the nine samples, which is evidence they are co-regulated, perhaps by the same signal or as part of the same pathway. It is not evidence that the two proteins touch, that one controls the other, or that either is causal. With only 9 samples and a response as strong as the 1-hour handling effect, plenty of pairs will correlate for reasons that have nothing to do with a shared function. Correlation networks generate hypotheses; they do not test them.
Would you expect a relationship? Yes. Both views are built from the same correlations between the same genes, so a group that really does move together has to show up in both, as a branch of the tree and as a knot in the web. If they came out completely unrelated, the conclusion would not be "interesting biology" but that one of the two is being driven by its settings. Expect the correspondence to be good but not exact, since a cut of the tree and a threshold on the edges are two different ways of drawing a line.
What each view gives you that the other cannot. Only from the network: a gene may belong to no group at all, or sit between two knots and belong to both, which is how you spot something linking two processes; it also gives you a number per gene describing where that one gene sits in relation to all the others, such as its degree (how many genes it is directly tied to, so whether it is a hub) or its betweenness (how often it lies on the path between other genes, so whether it is a bridge). Clustering cannot give you that at all: it tells you which group a gene is in, never where it sits within that group. Sliding the threshold also acts as a robustness check. Only from clustering: a complete partition, every gene in exactly one group, which is what you need to count anything or run an enrichment test on a defined list; a hierarchy rather than one slice; and an ordering of the samples, which is what let you read the design off the heatmap in step 2. Two views of one correlation matrix, so trust the structure that survives in both.
Which genes to study next, and what experiment. Say why, rather than just naming the top of the list: a gene with high degree is correlated with many others so may sit at the heart of a module, while high betweenness puts it on the paths between knots, making it a candidate for connecting two processes. Both are hypotheses about position in a graph, not measured facts. Concretely: take two or three central genes from the knot matching the AVR-specific cluster, check what is known about them in UniProt, obtain knockout and overexpression lines, repeat the infection with AVR, VIR and MOCK, and measure the disease along with the expression of the genes they were correlated with. If knocking the gene out flattens its neighbours, you have evidence it is upstream rather than merely correlated.