Clustering: Biological Applications

This page provides worked examples of clustering applied to common data types in biology and ecology. The methods are the same as those covered in the main clustering page; what changes is the preprocessing and the choice of distance or similarity measure.

---

Microbiome Sample Clustering

OTU and ASV tables are compositional count data, so the choice of transformation and distance measure matters. Bray-Curtis dissimilarity is a common choice for abundance data. When the compositional nature of the data is a primary concern, CLR transformation followed by Euclidean distance provides an Aitchison-distance approach. In practice, Bray-Curtis and Aitchison approaches answer slightly different questions about how compositional differences should be represented.

Here we use CLR transformation and Ward's method, which is appropriate because Ward's method operates in Euclidean space.

library(vegan)
library(phyloseq)
library(microbiome)
library(dendextend)

# Load example data
data(GlobalPatterns)
gp <- GlobalPatterns

# 1. Filter rare taxa
gp_filtered <- filter_taxa(
  gp,
  function(x) sum(x > 3) > (0.1 * length(x)),
  TRUE
)

# 2. CLR transformation
gp_clr <- microbiome::transform(gp_filtered, "clr")
otu_clr <- as.matrix(otu_table(gp_clr))

# 3. Aitchison distance
# Euclidean distance on CLR-transformed data
dist_aitchison <- dist(t(otu_clr), method = "euclidean")

# 4. Hierarchical clustering
hc <- hclust(dist_aitchison, method = "ward.D2")

# 5. Colour dendrogram by sample type
sample_meta <- as(sample_data(gp_filtered), "data.frame")
dend <- as.dendrogram(hc)

sample_types <- sample_meta[labels(dend), "SampleType"]
sample_levels <- unique(sample_types)

colors <- rainbow(length(sample_levels))
names(colors) <- sample_levels
labels_colors(dend) <- colors[as.character(sample_types)]

plot(dend, main = "Microbiome Sample Clustering")

The important lesson is that the clustering algorithm cannot be separated from the data representation. CLR transformation changes the geometry of the data, so the appropriate distance and clustering method change with it.

---

Gene Expression Clustering

RNA-seq count data contain strong mean-variance relationships, so raw counts should not be clustered directly. Variance-stabilising transformations such as the DESeq2 VST are commonly used before clustering.

Clustering can be applied both to samples, for example to detect batch effects or treatment groups, and to genes, for example to identify co-expression patterns.

library(DESeq2)
library(matrixStats)
library(pheatmap)
library(airway)

data("airway")

dds <- DESeqDataSet(
  airway,
  design = ~ cell + dex
)

# 1. Variance-stabilising transformation
vsd <- vst(dds, blind = FALSE)
mat <- assay(vsd)

# 2. Select the 50 most variable genes
rv <- rowVars(mat)
top_genes <- order(rv, decreasing = TRUE)[1:50]
mat_subset <- mat[top_genes, ]

# 3. Scale across samples
mat_scaled <- t(scale(t(mat_subset)))

# 4. Cluster genes and samples simultaneously
pheatmap(
  mat_scaled,
  clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean",
  clustering_method = "ward.D2",
  show_rownames = FALSE,
  annotation_col = as.data.frame(
    colData(dds)[, c("cell", "dex")]
  ),
  main = "Gene Expression Heatmap"
)

The column annotation helps identify whether samples cluster primarily by treatment (dex) or by cell line (cell). This is often one of the first checks for unexpected structure or batch effects in an RNA-seq dataset.

---

Single-Cell RNA-seq

Single-cell datasets contain thousands of cells and thousands of genes, making clustering in the original gene space impractical. A common workflow first reduces the dimensionality with PCA, then constructs a neighbour graph and performs graph-based clustering.

library(Seurat)

# Assuming a Seurat object 'pbmc' is already loaded

# 1. Normalise and find variable features
pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc)

# 2. Scale and run PCA
pbmc <- ScaleData(pbmc)
pbmc <- RunPCA(pbmc)

# 3. Build neighbour graph and cluster
# Clustering happens in PCA space, not raw gene space
pbmc <- FindNeighbors(pbmc, dims = 1:10)
pbmc <- FindClusters(pbmc, resolution = 0.5)

# 4. UMAP for visualisation
pbmc <- RunUMAP(pbmc, dims = 1:10)
DimPlot(pbmc, reduction = "umap", label = TRUE)

The resolution parameter influences clustering granularity. Lower values generally produce fewer, broader clusters, while higher values generally produce more, finer clusters.

Biological annotation comes after clustering. Marker gene expression is then examined to determine what biological cell types or states the clusters may represent.

---

Ecological Species Data

For ecological species abundance data, Bray-Curtis is a common choice when analysing abundance differences directly. Hellinger transformation provides another useful approach: after Hellinger transformation, Euclidean distance can be used while giving less weight to highly abundant species.

library(vegan)

data(dune)

# 1. Hellinger transformation
dune_hell <- decostand(dune, method = "hellinger")

# 2. Euclidean distance on Hellinger-transformed data
dist_hell <- dist(dune_hell, method = "euclidean")

# 3. Cluster sites
hc <- hclust(dist_hell, method = "ward.D2")

plot(hc, main = "Dune Meadow Sites")

# 4. Relate community structure to an independently defined
#    environmental variable
data(dune.env)

adonis2(
  dist_hell ~ Management,
  data = dune.env,
  permutations = 999
)

The clustering step explores the structure of the community data without using Management. PERMANOVA then asks a different question: does an independently defined environmental variable explain variation in community composition?

This illustrates the broader course workflow: first explore multivariate structure, then ask what might explain it.

---

Quick Reference: Data Type to Workflow

Data type	Transformation	Distance	Clustering / analysis
Continuous environmental variables	scale()	Euclidean	Ward.D2 or k-means
Species counts	Hellinger or none	Euclidean or Bray-Curtis	Ward.D2 or average linkage
OTU/ASV table	CLR	Euclidean (Aitchison)	Ward.D2
RNA-seq counts	VST	Euclidean	Hierarchical clustering
Single-cell RNA-seq	PCA scores	Neighbour graph	Graph-based clustering

The key lesson

There is no universally best clustering workflow. The appropriate transformation, distance measure, and clustering algorithm depend on the type of biological data and the question you are asking.