Expression Clustering — Finding Gene Co-Expression Groups Using Correlation and Ordering

After Completing This Topic

By combining the numpy, correlation, and sorting concepts learned in the textbook, you will be able to implement hierarchical clustering yourself to find co-expression groups in gene expression data across multiple samples. You will gain a code-level understanding of the principles behind scikit-learn or scipy.cluster.

This article is an educational general example. For real-world clustering, use pandas/scanpy/scikit-learn. Here, we focus on precisely understanding the underlying concepts of these tools.

"Out of these 5,000 genes, which ones move together?" — You can't see it with your eyes.

You have RNA-seq data. An expression matrix of 5,000 genes × 30 samples. You ask the following question:

• Which genes are turned on and off together?
• What biological function does this co-expression group correspond to?

It's impossible to compare 30-dimensional vectors of 5,000 genes visually. You need to automatically find and sort the relationships.

Naive Approach: Plot the expression profile for each gene. 5,000 graphs. You still won't see the relationships.

Real Approach: There are two axes:

1. Correlation: Quantify how much the expression profiles of two genes move together. Pearson correlation.
2. Hierarchical Clustering: Merge the most similar pairs into groups. The result is a dendrogram (phylogenetic tree-like).

By combining these two, groups of genes that move together are naturally represented in a hierarchical structure.

From Black Box to Components

Component 1: Expression Matrix and Correlation

import numpy as np


def load_expression_matrix(path: str) -> tuple[np.ndarray, list[str], list[str]]:
    """
    Returns: (matrix, gene_names, sample_names)
    matrix shape: (n_genes, n_samples)
    """
    with open(path) as f:
        header = f.readline().strip().split("\t")
        sample_names = header[1:]
        
        gene_names = []
        rows = []
        for line in f:
            parts = line.strip().split("\t")
            gene_names.append(parts[0])
            rows.append([float(x) for x in parts[1:]])
    
    return np.array(rows), gene_names, sample_names


expr, genes, samples = load_expression_matrix("expression.tsv")
print(f"Genes: {len(genes)}, Samples: {len(samples)}")
print(f"Matrix shape: {expr.shape}")

Now, the correlation matrix. Each element represents the Pearson correlation between two genes.

def compute_correlation_matrix(expr: np.ndarray) -> np.ndarray:
    """
    expr shape: (n_genes, n_samples)
    Returns: (n_genes, n_genes) correlation matrix
    """
    return np.corrcoef(expr)


corr = compute_correlation_matrix(expr)
print(f"corr shape: {corr.shape}")  # (5000, 5000)
print(f"Diagonal: {corr[0, 0]}")  # 1.0 (self-correlation)

np.corrcoef treats each row as a single variable. That is, if expr is of shape (genes, samples), the result is the gene-gene correlation.

Component 2: Distance Matrix

In clustering, we use distance instead of similarity. Correlation is a similarity measure ranging from -1 to +1, so we transform it into a distance using the following:

def correlation_to_distance(corr: np.ndarray) -> np.ndarray:
    """
    distance = 1 - correlation.
    High positive correlation → small distance.
    Negative correlation (move in opposite directions) → large distance.
    """
    return 1.0 - corr


dist = correlation_to_distance(corr)

Note: Other definitions exist. sqrt(2 * (1 - corr)) is a definition consistent with Euclidean distance. Choose based on the purpose.

Component 3: Hierarchical Clustering (Single Linkage)

Agglomerative approach: Initially, each gene is a separate cluster. We repeatedly merge the two closest clusters. This continues until only one cluster remains.

We need to define the distance between two clusters. Single linkage is the simplest definition: the minimum distance between any pair of elements in the two clusters.

def hierarchical_clustering_single(dist: np.ndarray) -> list[tuple[int, int, float]]:
    """
    Returns: A list of merge events: [(cluster_a, cluster_b, merge_distance), ...]
    After each merge, a new cluster ID is assigned as the original count plus the event index.
    """
    n = dist.shape[0]
    active_clusters = {i: [i] for i in range(n)}
    cluster_distances = {(i, j): dist[i, j] for i in range(n) for j in range(i + 1, n)}
    
    events: list[tuple[int, int, float]] = []
    next_id = n
    
    while len(active_clusters) > 1:
        best_pair = min(cluster_distances, key=cluster_distances.get)
        i, j = best_pair
        merge_dist = cluster_distances[best_pair]
        events.append((i, j, merge_dist))
        
        new_cluster = active_clusters[i] + active_clusters[j]
        del active_clusters[i]
        del active_clusters[j]
        active_clusters[next_id] = new_cluster
        
        # Calculate distances from the new cluster to the remaining clusters (single linkage = min)
        new_distances = {}
        for existing_id in active_clusters:
            if existing_id == next_id:
                continue
            existing = active_clusters[existing_id]
            min_d = min(
                dist[a, b] for a in new_cluster for b in existing
            )
            new_distances[(min(existing_id, next_id), max(existing_id, next_id))] = min_d
        
        cluster_distances = {
            k: v for k, v in cluster_distances.items()
            if i not in k and j not in k
        }
        cluster_distances.update(new_distances)
        
        next_id += 1
    
    return events

Time complexity: The naive implementation is O(n³). With 5,000 genes, this is difficult to handle, and in practice, we use O(n² log n) algorithms. This tutorial focuses on conceptual understanding, so we use the naive implementation.

Component 4: Extracting Clusters from the Result

We extract a specific number of clusters from the list of merge events.

def cut_dendrogram(events: list[tuple[int, int, float]], n_leaves: int, num_clusters: int):
    parent = list(range(n_leaves + len(events)))
    
    def find(x):
        while parent[x] != x:
            parent[x] = parent[parent[x]]
            x = parent[x]
        return x
    
    def union(x, y, new_id):
        rx, ry = find(x), find(y)
        parent[rx] = new_id
        parent[ry] = new_id
    
    # Ignore the last (num_clusters - 1) merges and apply the union up to that point
    events_to_apply = events[:-num_clusters + 1] if num_clusters > 1 else events
    
    next_id = n_leaves
    for a, b, _ in events_to_apply:
        union(a, b, next_id)
        next_id += 1
    
    clusters = {}
    for leaf in range(n_leaves):
        root = find(leaf)
        clusters.setdefault(root, []).append(leaf)
    
    return list(clusters.values())


cluster_lists = cut_dendrogram(events, n_leaves=len(genes), num_clusters=10)
for i, members in enumerate(cluster_lists):
    print(f"Cluster {i}: {len(members)} genes")

The Union-Find data structure allows for O(α(n)) ≈ O(1) merging and querying. This is useful for large-scale problems.

Sorting and Sorted Heatmap

The visual result of clustering becomes apparent in the reordered heatmap. The genes within each cluster are sorted so that they become adjacent rows.

def cluster_order(cluster_lists: list[list[int]]) -> list[int]:
    """Sorts the members of each cluster while preserving the original index order within each cluster."""
    result = []
    for members in cluster_lists:
        result.extend(sorted(members))
    return result


new_order = cluster_order(cluster_lists)
reordered_expr = expr[new_order]

When this reordered matrix is plotted as a heatmap, each cluster appears as a distinct block.

import matplotlib.pyplot as plt


def plot_heatmap(matrix: np.ndarray, ax=None) -> None:
    if ax is None:
        _, ax = plt.subplots(figsize=(6, 8))
    im = ax.imshow(matrix, cmap="RdBu_r", aspect="auto")
    plt.colorbar(im, ax=ax)


fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 8))
plot_heatmap(expr, ax1)
ax1.set_title("Original")
plot_heatmap(reordered_expr, ax2)
ax2.set_title("Clustered")
plt.tight_layout()
plt.show()

This before-and-after comparison of the reordering visually demonstrates the power of clustering.

Fading – Three Blanks for You to Fill In

Blank 1: Complete Linkage / Average Linkage

Implement either complete linkage (the maximum distance between all pairs of elements in two clusters) or average linkage (the average distance) instead of single linkage.

def hierarchical_clustering_complete(dist: np.ndarray) -> list[tuple[int, int, float]]:
    # Mostly identical to the single linkage code, replace min with max
    # ...
    # Cluster distance calculation section:
    # TODO: Replace min(dist[a, b] for ...) with max(dist[a, b] for ...)
    pass

Hint: Parameterize only one function. linkage: Callable[[list[float]], float] = min or max or (lambda xs: sum(xs) / len(xs)).

Blank 2: Representative Gene of a Cluster

In each cluster, the gene with the highest average correlation with other members is the representative gene.

def find_hub_genes(
    corr: np.ndarray,
    clusters: list[list[int]],
    gene_names: list[str]
) -> list[tuple[str, float]]:
    """Representative gene and its average correlation for each cluster."""
    hubs = []
    for members in clusters:
        # TODO: Calculate the average correlation of each member with other members
        # Select the gene with the highest value as the hub
        pass
    return hubs

Hint:

best_score = -1
best_gene = None
for m in members:
    score = np.mean([corr[m, other] for other in members if other != m])
    if score > best_score:
        best_score = score
        best_gene = gene_names[m]

Blank 3: Functional Annotation (Integration with External Tool)

Map the list of genes in a cluster to KEGG pathways or GO terms. Here, we load the mapping from a CSV file.

def enrich_clusters(
    clusters: list[list[str]],
    gene_to_pathway_csv: str
) -> dict:
    """
    Top 3 pathways that appear most frequently in each cluster.
    """
    # TODO: Load the CSV (gene, pathway columns)
    # Calculate the frequency of pathways in each cluster
    # Return the top 3 pathways for each cluster
    pass

Hint: from collections import Counter; counter = Counter(); for g in members: counter.update(gene_pathways.get(g, [])).

Reflections – Differences from Real-World Co-Expression Tools

scipy.cluster.hierarchy: Real-world tools use the linkage function from this module. It's an O(n² log n) algorithm implemented in C. This is over 100 times faster than your naive O(n³) implementation.

WGCNA: A standard tool for gene co-expression networks. It uses sophisticated concepts such as signed adjacency matrices, soft thresholding, and module preservation. Your tool represents a minimal framework for this.

scanpy: A standard stack for single-cell RNA-seq. Clustering is typically done using the Leiden algorithm – a graph-based approach, which is a different family of algorithms from the hierarchical clustering you implemented.

Batch effect: Real-world expression data often involves combining data from different experimental batches. Batch effects can contaminate the true co-expression signal. Tools like Combat and Harmony address this issue.

Dynamic tree cutting: Your cut_dendrogram specifies the number of clusters. Real-world WGCNA uses an algorithm to dynamically cut the dendrogram based on its shape (Dynamic Tree Cut).

Extension Project

1. Reimplementation with SciPy: Replace your naive clustering with scipy.cluster.hierarchy.linkage and compare performance.

2. Dendrogram Visualization: Visualize the tree using scipy.cluster.hierarchy.dendrogram. Convert your list of events into the SciPy format.

3. GO Term Enrichment: Perform actual GO term enrichment for each cluster using gseapy.

4. Streamlit Dashboard: Create a dashboard where users can adjust the number of clusters using a slider, and the heatmap and representative genes are updated in real-time.

Component Guide for This Module

• [F] numpy: Manipulation of expression matrices and correlation matrices. Practical application of np.corrcoef.
• [F] Correlation: Pearson correlation is used as a quantitative measure of similarity between two genes.
• [F] Sorting: Heatmap blocks are created by reordering clusters. Results are extracted using Union-Find.
• [W] matplotlib: Heatmap visualization (complete script provided).

[F] = You will implement this yourself / [W] = Provided as complete code.