RNA-seq Heatmap: Combining Normalization, Clustering, and Visualization into One

After Completing This Topic

By combining the normalization, hierarchical clustering, and seaborn heatmap techniques learned in the textbook, you can create a tool that takes an RNA-seq counts matrix as input, automatically normalizes it, groups samples and genes with similar expression patterns side-by-side, and generates a heatmap suitable for inclusion in a publication. You will understand that this beautiful image is actually the result of three layers combined.

This article is a general educational example. RNA-seq was chosen as the subject because it is a standard method in current expression studies.

"If I Just Change the Order of the Samples, the Result Looks Completely Different?" — The Pitfalls of Normalization and Ordering

After an RNA-seq experiment, you'll get a counts matrix like this:

gene       tumor1  tumor2  tumor3  normal1  normal2  normal3
BRCA1      1520    1830    1650    340      420      380
TP53       850     920     780     650      710      680
MYC        4200    5100    4700    1200     1350     1280
GAPDH      18000   19500   17800   17200    17900    18300
ACTB       15000   16800   15500   14800    15600    15200
...

What happens if you simply create a heatmap from this table?

• BRCA1 has values in the hundreds to thousands. It will appear as faint blue cells that are hard to see.
• GAPDH and ACTB have values in the tens of thousands. The entire heatmap will be dominated by these two genes, appearing as bright red rows.
• The order of the sample columns was set to the order of the experiment, so tumor and normal samples are mixed, making it difficult to see any patterns.

From the moment you start encountering these problems, you'll realize that you need three different processing steps.

1. Normalization — To make the expression ranges comparable between genes and samples.
2. Ordering (Clustering) — To group similar samples together so that patterns become visible.
3. Color Mapping — To convert values into visually distinguishable colors.

Only when these three layers are well-integrated can you create figures that are suitable for publication. From now on, we will create and assemble these three layers one by one.

Let's See the Final Product First (Run the Black Box First)

Here's how we'll use the tool we're going to create:

plot_rnaseq_heatmap(
    counts_df,
    top_n_variable=50,  # Only use the top 50 genes with the largest variance
    cluster_samples=True,
    cluster_genes=True,
    save_path="heatmap.png",
)

The resulting heatmap conveys the following information:

=== Heatmap Story ===
1. The 3 tumor samples are automatically clustered together on the left.
2. The 3 normal samples are automatically clustered together on the right.
3. A group of genes that are particularly highly expressed in tumors are clustered at the top.
4. A group of genes that are highly expressed in normal samples are clustered at the bottom.
5. Color = z-score (standardized value based on variance between samples).
6. Top dendrogram = sample hierarchy.
7. Left dendrogram = gene hierarchy.

File saved: heatmap.png (300dpi, suitable for inclusion in a paper).

Each time we use the same data, we get the same plot. Reproducibility is captured in a single filename.

What components make up this tool (component breakdown)?

RNA-seq Heatmap Pipeline
   ┌──────────────────────────────────────────────────┐
   │  [Input] Load counts matrix ────── Component: pandas       │  ← Provided as is (tool)
   │              │                                      │
   │              ▼                                      │
   │  [Step 1] Normalization (CPM · z-score)                       │
   │        Component: normalization                            │  ← Created from scratch ★
   │              │                                      │
   │              ▼                                      │
   │  [Step 2] Select genes with large variance                            │
   │        Component: pandas variable calculation                          │  ← Provided as is (tool)
   │              │                                      │
   │              ▼                                      │
   │  [Step 3] Sample and gene hierarchical clustering                    │
   │        Component: hierarchical clustering                   │  ← Created from scratch ★
   │              │                                      │
   │              ▼                                      │
   │  [Output] seaborn clustermap (heatmap + dendrogram)         │  ← Assembled from scratch ★
   └──────────────────────────────────────────────────┘

📌 If you are new to these concepts (links at the top)

• Normalization — the principle of unifying scales
• Hierarchical clustering — the growth rules of dendrograms
• seaborn heatmap — color value mapping

The three new concepts created from scratch are normalization, hierarchical clustering, and heatmap assembly. pandas is provided as a complete tool. Just three — within the limits of cognitive capacity.

Step 1: Creating the Dataset (Provided as Complete)

Let's create a dataset that mimics a real RNA-seq counts matrix. In practice, this data comes out as a CSV file from tools like Salmon, featureCounts, or HTSeq.

import numpy as np
import pandas as pd

rng = np.random.default_rng(42)

def make_rnaseq_counts():
    genes = ["BRCA1", "TP53", "MYC", "KRAS", "EGFR", "PIK3CA", "PTEN", "APC",
             "GAPDH", "ACTB", "B2M", "HPRT1", "TBP"]  # Last 5 are housekeeping genes
    tumor_samples = [f"tumor{i+1}" for i in range(3)]
    normal_samples = [f"normal{i+1}" for i in range(3)]
    samples = tumor_samples + normal_samples

    counts = np.zeros((len(genes), len(samples)), dtype=int)
    for gi, gene in enumerate(genes):
        # Housekeeping genes: similar in both conditions
        base_tumor  = 1000 if gene in {"BRCA1", "MYC", "KRAS"} else 300
        base_normal = 300 if gene in {"BRCA1", "MYC", "KRAS"} else 300
        if gene in {"GAPDH", "ACTB", "B2M", "HPRT1", "TBP"}:
            base_tumor = base_normal = 15000  # Housekeeping
        for si, sample in enumerate(samples):
            base = base_tumor if sample.startswith("tumor") else base_normal
            counts[gi, si] = max(0, int(base * rng.lognormal(0, 0.15)))
    return pd.DataFrame(counts, index=genes, columns=samples)

counts = make_rnaseq_counts()
print(counts)

# Verification
assert counts.shape == (13, 6)
assert (counts >= 0).all().all()
# BRCA1 in tumor should be higher than in normal (by design)
assert counts.loc["BRCA1", "tumor1"] > counts.loc["BRCA1", "normal1"]
# Housekeeping genes should be similar in both conditions
housekeeping_tumor = counts.loc["GAPDH", ["tumor1","tumor2","tumor3"]].mean()
housekeeping_normal = counts.loc["GAPDH", ["normal1","normal2","normal3"]].mean()
assert abs(housekeeping_tumor / housekeeping_normal - 1) < 0.3

Step 2: CPM Normalization ★ (Normalization: Comparison between samples)

✍️ Fill in this section yourself. Component = Normalization. Since each sample has a different total number of reads, absolute counts cannot be compared. We convert them into relative proportions.

The first reason for normalization is the difference in the total number of reads between samples. If one sample has 30 million reads due to good sequencing, and another sample has only 20 million reads, you cannot directly compare the counts of the same gene.

CPM (Counts Per Million): Normalizes by the total number of reads, scaling to one million.

CPM(gene, sample) = counts(gene, sample) / total_counts(sample) × 1,000,000

def to_cpm(counts):
    """Divide by the total counts per sample and scale to one million."""
    library_sizes = counts.sum(axis=0)     # Total counts for each sample
    cpm = counts.div(library_sizes, axis=1) * 1_000_000
    return cpm

cpm = to_cpm(counts)
print(cpm.round(1).head())

# Verification: The sum of CPM for each sample should be exactly 1,000,000
assert np.allclose(cpm.sum(axis=0).values, 1_000_000)
# Comparison: Comparisons that were not possible with raw counts are now possible with CPM
tumor_avg = cpm[["tumor1","tumor2","tumor3"]].mean(axis=1)
normal_avg = cpm[["normal1","normal2","normal3"]].mean(axis=1)
# BRCA1 is indeed more highly expressed in the tumor
assert tumor_avg["BRCA1"] > normal_avg["BRCA1"]

Now it is possible to compare samples. However, comparison between genes is still not possible. The CPM for GAPDH is still 100 times greater than the CPM for BRCA1. The Z-score will solve this.

🔎 Why is CPM so simple (drawer — normalization) It is the simplest form of normalization, "scaling to a standard size." In practice, more sophisticated methods such as TPM, TMM, and the median-of-ratios from DESeq2 are used. The conceptual root is the same — to make it comparable.

🤔 Self-explanatory prompt Why is CPM safer than using raw counts directly? When the total number of reads in two samples is 50 million and 10 million, how can you compare the expression of gene A using raw counts? (Hint: You can't.)

Step 3: Z-score Normalization ★ (Normalization: Comparing Genes)

✍️ Fill-in section. Component = Normalization. Standardize each gene by its mean and standard deviation to eliminate scale differences between genes.

For each gene row:

z = (value - mean) / std

This expresses how many standard deviations away from its mean the expression of each gene is. This makes it possible to compare GAPDH and BRCA1 on the same scale.

def to_zscore(cpm):
    """Z-score per gene. Standardize after log transformation (common practice)."""
    log_cpm = np.log2(cpm + 1)  # +1 prevents log(0)
    means = log_cpm.mean(axis=1)
    stds = log_cpm.std(axis=1, ddof=0)
    # Safely set standard deviation to 0 for genes with zero variance (all values are the same)
    stds = stds.replace(0, 1)
    return log_cpm.sub(means, axis=0).div(stds, axis=0)

zscore = to_zscore(cpm)
print(zscore.round(2))

# Verification: Mean ≈ 0 and standard deviation ≈ 1 for each gene
row_means = zscore.mean(axis=1)
row_stds = zscore.std(axis=1, ddof=0)
assert np.allclose(row_means.values, 0, atol=1e-9)
# Check only genes with non-zero standard deviation (excluding those where standard deviation was set to 0)
active_genes = cpm.std(axis=1) > 0
assert np.allclose(row_stds[active_genes].values, 1, atol=0.05)
# By design, tumor samples for BRCA1 should be positive, and normal samples should be negative
assert (zscore.loc["BRCA1", ["tumor1", "tumor2", "tumor3"]] > 0).all()
assert (zscore.loc["BRCA1", ["normal1", "normal2", "normal3"]] < 0).all()

Now all genes are on the same scale (mean 0, standard deviation 1). The heatmap colors will freely span for each gene. This standardization is the key to determining the information density of the heatmap.

🤔 Self-explanatory prompt Why was log2(cpm + 1) applied before the Z-score? What would happen if you simply applied the Z-score to the raw CPM? (Hint: Expression is naturally on a logarithmic scale, and large values can be outliers that dominate the standard deviation.)

Step 4: Select Only Genes with High Variance (Provided Complete)

To ensure the heatmap contains meaningful information, we must eliminate genes that do not change. Genes that consistently have the same value regardless of the condition only create white lines and noise.

def select_top_variable(zscore, top_n=50):
    """Selects the top 'top_n' genes based on the range of their z-scores."""
    # Here, we use range (max - min) instead of standard deviation – prioritizing genes with strong positive and negative deviations.
    variability = zscore.max(axis=1) - zscore.min(axis=1)
    top_genes = variability.sort_values(ascending=False).head(top_n).index
    return zscore.loc[top_genes]

top = select_top_variable(zscore, top_n=8)
print(top.index.tolist())

# Verification: Only 8 genes remain
assert len(top) == 8
# Housekeeping genes (genes with little variation) should be excluded from the top 8
housekeeping = {"GAPDH", "ACTB", "B2M", "HPRT1", "TBP"}
selected_housekeeping = set(top.index) & housekeeping
# Most housekeeping genes should be excluded (some may be included as noise)
assert len(selected_housekeeping) <= 2

Step 5: Hierarchical Clustering ★ (clustering-hierarchical)

✍️ Fill-in section. Component = Hierarchical Clustering. It calculates the distance between two vectors and sequentially groups them to create a tree.

Idea: Since each sample (or gene) is a vector, we can measure the distance between vectors. We group the two closest ones into a group and treat that group as a single vector, continuing to group them. The result is a dendrogram — a hierarchical tree.

This method is available in a complete form in SciPy.

🔎 What is Hierarchical Clustering (Drawer — clustering-hierarchical)? "Start with each point in its own group → repeatedly merge the two closest groups → until only one remains." Cutting the resulting tree gives you the desired number of clusters. Unlike k-means, you don't need to specify the number of clusters in advance.

from scipy.cluster.hierarchy import linkage, dendrogram

def cluster_order(data, axis="rows", method="average", metric="euclidean"):
    """Returns the index reordered by hierarchical clustering."""
    matrix = data.values if axis == "rows" else data.values.T
    Z = linkage(matrix, method=method, metric=metric)
    # leaves_list is the optimal order of the tree leaves (original)
    from scipy.cluster.hierarchy import leaves_list
    order = leaves_list(Z)
    if axis == "rows":
        return data.index[order]
    return data.columns[order]

# Reorder samples and genes by cluster order
sample_order = cluster_order(top, axis="cols")
gene_order   = cluster_order(top, axis="rows")

reordered = top.loc[gene_order, sample_order]
print("Sample order:", list(sample_order))
print("Gene order:", list(gene_order))

# Verification: The sum of values in the reordered table is the same as the original
assert np.allclose(reordered.values.sum(), top.values.sum())
# By design, 3 tumors and 3 normal samples should be clustered together
sample_names = list(sample_order)
tumor_positions = [i for i, s in enumerate(sample_names) if s.startswith("tumor")]
normal_positions = [i for i, s in enumerate(sample_names) if s.startswith("normal")]
# Tumors should be in consecutive positions
assert max(tumor_positions) - min(tumor_positions) == 2

The 3 tumors are now automatically grouped together next to each other. Even though we didn't explicitly tell the clustering algorithm about the tumor/normal labels, the data itself reveals these groups. This is the value of clustering — discovering structure without labels.

🤔 Self-explanatory prompt I chose method="average" and metric="euclidean". What would happen if I changed it to method="ward"? Or if metric="correlation"? Try changing it and comparing the results.

Step 6: Building the seaborn Clustermap ★

seaborn's clustermap function combines the clustering and heatmap visualization into a single step. Now that we understand the underlying principles of how we assembled these components manually, we can more easily grasp what each parameter of this function does.

import matplotlib
matplotlib.use("Agg")  # For Pyodide/server compatibility
import matplotlib.pyplot as plt
import seaborn as sns

def plot_rnaseq_heatmap(counts_df, top_n=50, save_path=None):
    """Complete pipeline: normalization → top variable selection → clustermap."""
    cpm = to_cpm(counts_df)
    z = to_zscore(cpm)
    top = select_top_variable(z, top_n=top_n)

    g = sns.clustermap(
        top,
        cmap="RdBu_r",           # Red-blue (clear positive/negative)
        center=0,                # Center at 0 (mean)
        vmin=-2, vmax=2,         # Fixed color range
        method="average",
        metric="euclidean",
        figsize=(10, max(5, top_n * 0.15)),
        cbar_kws={"label": "z-score (log2 CPM)"},
    )
    if save_path:
        g.savefig(save_path, dpi=300, bbox_inches="tight")
    return g

# Execution (Colab/Jupyter)
g = plot_rnaseq_heatmap(counts, top_n=8, save_path=None)
plt.close("all")

# Verification: Ensure the returned grid is a seaborn clustermap object
from seaborn.matrix import ClusterGrid
assert isinstance(g, ClusterGrid)
# Verify that the reordered data has the expected shape
assert g.data2d.shape == (8, 6)

This single function assembles all the normalization, clustering, and color mapping that we built in the previous steps. Knowing the principles behind each layer allows you to confidently tune the parameters.

Combining the Pieces: The Complete Dashboard Class

Encapsulate the entire pipeline into a single class.

class RNAseqDashboard:
    def __init__(self, counts_df: pd.DataFrame):
        self.counts = counts_df
        self.cpm = to_cpm(counts_df)
        self.zscore = to_zscore(self.cpm)

    def summary(self) -> pd.DataFrame:
        return pd.DataFrame({
            "n_genes": [len(self.counts)],
            "n_samples": [self.counts.shape[1]],
            "total_reads_mean": [self.counts.sum(axis=0).mean()],
            "top_variable_gene": [
                (self.zscore.max(axis=1) - self.zscore.min(axis=1)).idxmax()
            ],
        })

    def plot(self, top_n=50, save_path=None):
        return plot_rnaseq_heatmap(self.counts, top_n=top_n, save_path=save_path)

dash = RNAseqDashboard(counts)
print(dash.summary())

# Validation
assert dash.summary()["n_genes"].iloc[0] == 13
assert dash.summary()["n_samples"].iloc[0] == 6
# The most variable gene should be one that truly has differences between conditions (BRCA1, MYC, or KRAS)
assert dash.summary()["top_variable_gene"].iloc[0] in {"BRCA1", "MYC", "KRAS"}

This tool is a simplified version of the RNA-seq expression heatmap that you often see in papers today. Practical tools (like pheatmap in R, ComplexHeatmap, etc.) simply add color annotations, custom distance functions, statistical test overlays, and more on top of this.

There Are Other Paths (Multi-Pass Reflection)

• TPM vs. CPM: CPM only normalizes by sample size. TPM (Transcripts Per Million) also normalizes for gene length, making it more accurate for comparisons between genes. When and what: For comparisons between samples, use CPM; for absolute comparisons between genes, use TPM.
• DESeq2 Normalization: If the goal is differential expression analysis, use DESeq2's median-of-ratios or TMM (edgeR) instead of CPM/TPM. These methods better satisfy the assumptions of statistical tests.
• k-means Alternative: Hierarchical clustering has O(n²) memory complexity, which is impractical when dealing with tens of thousands of genes. In such cases, use k-means or UMAP + HDBSCAN instead. The goal remains the same: to discover structure in unlabeled data.
• Choosing a Distance Function: euclidean is sensitive to the absolute magnitude of values. correlation focuses only on the pattern shape (ups and downs). After z-score normalization, the two methods become similar, but knowing the alternatives allows you to choose the best one for the situation.
• Risks with Small Sample Sizes: When the sample size is very small (e.g., 3 + 3), clustering can be easily affected by noise. When including such heatmaps in a paper, be sure to reflect the fact that the sample size is small in the interpretation of the results.

Key takeaway: "Normalization enables comparability, clustering reveals structure, and heatmaps provide visualization." If each of these layers understands its role, the path from the original data to the discovery can be reproduced. What you just created is the embodiment of that path.

Next Steps (Links at the Bottom)

• More details on various normalization methods → Normalization Principles
• Why hierarchical clustering can find structure without labels → Hierarchical Clustering
• Turn the qPCR results from earlier into a heatmap → Automated qPCR Analysis

Try it Yourself (Independent Exercises)

1. TPM Extension: Implement to_tpm(counts, gene_lengths) to normalize using TPM instead of CPM, by additionally taking gene length information (length CSV) as input.
2. Annotation Column: Add a color bar for tumor/normal conditions to the clustermap using the col_colors parameter. (Hint: Create and pass a series of colors for each condition.)
3. Selected Genes Only: Add an option that allows the user to specify a list of genes of interest directly (instead of using top_n).
4. Challenge — Correlation Heatmap: Calculate the gene-gene correlation matrix (Pearson correlation) and create a clustermap. Use it to discover co-expression patterns.

Summary

We've conquered the problem of "creating an informative heatmap from RNA-seq results" by breaking it down into three parts.

• Normalization enabled meaningful comparisons between genes and samples.
• Hierarchical clustering automatically discovered the tumor/normal group structure without labels.
• Seaborn clustermap assembled the results into a visually appealing format suitable for inclusion in a publication.

The magic behind those beautiful heatmaps we see in publications has now been revealed. It's simply the result of three layers each performing its specific role. With the pipeline you've created, you can discover your own stories within your data.

This article is a general educational example. Real-world RNA-seq analysis (using DESeq2, edgeR, Bioconductor, etc.) involves adding statistical tests, batch effect correction, GO enrichment, and more. You can either build these additional components on top of this framework or rely on well-validated libraries.