What is Dr. Zubair Khalid's research focus?

Dr. Zubair Khalid specializes in molecular virology, mRNA vaccine development, and computational biology, with a focus on avian pathogens like IBDV and Avian Reovirus.

Where is Dr. Zubair Khalid currently working?

Dr. Zubair Khalid is a Postdoctoral Research Associate at the University of Maryland (UMD), specifically within the Department of Animal and Avian Sciences.

Heatmap of Differentially Expressed Genes: A Practical Guide

By Dr. Zubair Khalid, DVM, MS, PhD · 2026-08-11

Introduction to Heatmaps in Differential Expression Analysis

What is a differential expression heatmap?

A heatmap of differentially expressed genes is a two-dimensional graphical representation of a gene expression matrix in which individual expression values are encoded as colors. The rows represent genes, the columns represent biological samples (or conditions), and the color intensity at each intersection reflects the abundance of that gene's transcript in that sample. In the context of differential expression analysis, the heatmap is typically restricted to genes that pass statistical thresholds for differential expression between experimental groups, and the displayed values are usually transformed—often to Z-scores or variance-stabilized counts—so that the visual pattern emphasizes relative changes rather than absolute abundance.

The core utility of this visualization lies in its capacity to compress thousands of data points into a single, human-readable image that reveals structure: which genes behave similarly across conditions, which samples cluster together, and whether the transcriptional program aligns with the experimental design. A well-constructed heatmap is not merely a figure; it is an analytical tool that can expose batch effects, reveal co-regulated gene modules, and generate hypotheses about regulatory mechanisms.

Why heatmaps are essential in omics research

Transcriptomic experiments routinely measure 20,000 or more genes across dozens of samples. Raw tables of fold changes and p-values are unwieldy and fail to convey the gestalt of the biological response. Heatmaps solve this by providing an immediate visual summary of the data's dominant patterns. They are indispensable for quality assessment—a heatmap of the top variable genes will often reveal sample mislabeling, batch structure, or outlier libraries before any formal statistical test is applied.

Beyond quality control, heatmaps serve as the primary vehicle for communicating differential expression results in publications. Reviewers and readers expect to see a heatmap that demonstrates the reproducibility of biological replicates, the separation between treatment groups, and the coherence of gene expression programs. In the context of Biomarker Discovery, heatmaps are frequently used to identify gene signatures that discriminate between disease states or patient subgroups. Similarly, in studies of Target Genes of Nf KB Signalling, a heatmap of NF-κB-responsive genes across time points or perturbations can reveal the kinetics and magnitude of the transcriptional response in a way that no table of statistics can match.

Data Preprocessing for Reliable Heatmaps

The quality of a heatmap is entirely determined by the quality of the data that feeds it. Generating a heatmap from raw count data without proper preprocessing is a common and serious error. The following steps are mandatory before any visualization.

Quality control and filtering

Before normalization, each sample must pass quality control. For RNA-seq data, this includes assessing sequencing depth (total mapped reads), the proportion of reads mapping to exonic regions, the number of genes detected at a minimum threshold (commonly at least 10 counts in at least a defined fraction of samples), and the absence of anomalous GC bias or 3' bias. Tools such as FastQC, MultiQC, and RSeQC provide these metrics. Samples with fewer than 10–20 million mapped reads, or with an unusually low gene detection rate, should be flagged and potentially excluded.

After sample-level QC, gene-level filtering is required. Genes with very low counts across all samples contribute noise and cannot be reliably assessed for differential expression. A standard filter retains genes with a minimum count of 10–15 in at least a specified number of samples (e.g., the smallest group size). For example, in a comparison of three treated versus three control samples, a gene must have at least 10 counts in at least three samples to be retained. This filter removes genes whose expression is essentially zero and reduces the multiple testing burden.

Normalization methods (TPM, RPKM, DESeq2, edgeR)

Normalization is the process of making gene expression values comparable within and between samples. The choice of normalization method depends on the downstream analysis and the biological question.

TPM (Transcripts Per Million) and RPKM/FPKM (Reads/Fragments Per Kilobase of transcript per Million mapped reads) are length-normalized measures that allow comparison of expression levels between genes within a sample. TPM is preferred over RPKM/FPKM because it is calculated after dividing read counts by gene length, then normalizing to the total sum of length-normalized counts, making the sum of all TPM values equal to one million across samples. However, these measures are not appropriate for between-sample differential expression testing because they do not account for differences in library composition—a few highly expressed genes can skew the total count and create false differences.

For differential expression analysis, DESeq2 and edgeR use raw count data and apply their own internal normalization. DESeq2 uses the median-of-ratios method: for each gene, the geometric mean of counts across samples is computed, and the ratio of each sample's count to this geometric mean is calculated. The median of these ratios across all genes is the sample's size factor, which is used to scale the counts. edgeR uses the trimmed mean of M-values (TMM) method, which calculates a weighted mean of log fold changes between samples after trimming the most extreme values. Both methods are robust to compositional differences and are the standard for identifying differentially expressed genes.

For heatmap visualization, the choice of normalization matters. If you are displaying raw expression values, TPM or FPKM is appropriate. However, for differential expression heatmaps, it is almost always better to display a variance-stabilized or regularized log transformation of the DESeq2-normalized counts, as described below.

Variance stabilizing transformation

Raw count data exhibit a mean-variance relationship: genes with higher expression have higher variance. This heteroscedasticity is problematic for clustering and visualization because highly expressed genes will dominate the distance calculations. The variance stabilizing transformation (VST) and the regularized logarithm (rlog) in DESeq2 address this by applying a transformation that stabilizes the variance across the range of expression values.

The VST uses the fitted dispersion-mean relationship to find a transformation that makes the variance approximately constant. The rlog is similar but uses a shrinkage approach that borrows information across genes, making it particularly useful for small sample sizes. After transformation, the data are approximately homoscedastic and can be used directly for clustering and heatmap generation. For edgeR, the cpm (counts per million) values with a log2 transformation and a small offset (e.g., log2(cpm + 1)) are often used, though the voom transformation provides a mean-variance modeling approach that is more principled.

For a heatmap of differentially expressed genes, the standard workflow is: (1) run DESeq2 or edgeR to identify significant genes, (2) extract the VST or rlog-transformed values for those genes, and (3) optionally Z-score normalize each gene across samples so that the mean is 0 and the standard deviation is 1. Z-scoring is critical when genes have vastly different baseline expression levels; without it, highly expressed genes will dominate the color scale and subtle changes in lowly expressed genes will be invisible.

Selecting Differentially Expressed Genes

The genes that enter the heatmap define its biological message. Selecting too many genes produces a noisy, unreadable figure; selecting too few risks missing important biology. The selection process is a balance between statistical rigor and biological relevance.

Statistical cutoffs (p-value, FDR)

The first filter is statistical significance. For each gene, a p-value is computed from the differential expression test (e.g., the Wald test in DESeq2 or the exact test in edgeR). Because thousands of genes are tested simultaneously, the p-value must be adjusted for multiple testing. The Benjamini-Hochberg (BH) procedure controls the false discovery rate (FDR) and is the default in most tools. A common threshold is an adjusted p-value (padj) of less than 0.05, though more stringent thresholds (padj < 0.01) are appropriate when the number of differentially expressed genes is very large and the goal is to focus on the most robust changes.

Effect size (log2 fold change)

Statistical significance alone is insufficient. A gene with a padj of 1e-10 but a log2 fold change of 0.1 is statistically robust but biologically trivial. Conversely, a gene with a log2 fold change of 5 but a padj of 0.1 is likely too noisy to trust. The standard practice is to apply a combined threshold: padj < 0.05 and |log2 fold change| ≥ 1 (equivalent to a 2-fold change). For some applications, a more stringent fold change cutoff (|log2FC| ≥ 2) is used to focus on the most dramatic transcriptional changes.

The choice of fold change cutoff should be guided by the expected biology. For example, transcription factor perturbation experiments often produce large fold changes in direct targets, while subtle metabolic rewiring may involve many genes with modest fold changes. In the latter case, a lower fold change threshold (|log2FC| ≥ 0.5) combined with a stringent FDR may be more appropriate.

Filtering for biological significance

Statistical filters do not guarantee biological relevance. After applying the statistical cutoffs, it is common to further filter genes based on their expression level (e.g., requiring a minimum mean normalized count of 50) to exclude genes whose expression is so low that the fold change is unreliable. Additionally, if the heatmap is intended to tell a specific story—for example, focusing on Target Genes of Nf KB Signalling—the gene list can be intersected with a curated gene set from a database such as MSigDB or KEGG.

A practical approach is to rank the significant genes by a combined score, such as the product of the log2 fold change and the negative log10 of the adjusted p-value, and then take the top N genes (e.g., the top 50–200) for the heatmap. This ensures that the displayed genes are both statistically robust and biologically meaningful.

Heatmap Construction: Algorithms and Tools

Hierarchical clustering

The defining feature of most heatmaps is the dendrogram—the tree structure that groups similar genes and similar samples. Hierarchical clustering is an agglomerative process that begins with each gene (or sample) as its own cluster and iteratively merges the two most similar clusters until a single cluster remains. The result is a binary tree whose leaves are the individual genes or samples.

The clustering is performed on the transformed, Z-scored expression matrix. For genes, the rows are clustered; for samples, the columns are clustered. The order of the leaves in the dendrogram is determined by the clustering algorithm and reflects the similarity structure of the data. The dendrogram is then used to reorder the rows and columns of the heatmap so that similar genes and similar samples are adjacent.

Distance metrics (Euclidean, Pearson correlation)

The distance metric defines what "similar" means. The two most common choices are Euclidean distance and Pearson correlation.

Euclidean distance between two genes is the square root of the sum of squared differences in their expression values across samples. It is sensitive to the magnitude of expression: two genes with identical patterns but different baseline levels will have a large Euclidean distance. This is why Z-scoring is essential before using Euclidean distance—it centers each gene at zero and scales to unit variance, making the distance reflect pattern similarity rather than absolute abundance.

Pearson correlation measures the linear association between two genes' expression profiles. A correlation of 1 indicates perfect positive correlation, -1 indicates perfect negative correlation, and 0 indicates no correlation. For clustering, the distance is often defined as 1 - correlation, so that highly correlated genes have small distances. Pearson correlation is invariant to linear transformations (adding a constant or multiplying by a constant), which means it naturally handles differences in baseline expression without requiring Z-scoring. However, it is sensitive to outliers and can be misleading when the expression profiles are not approximately linear.

In practice, Euclidean distance on Z-scored data and Pearson correlation-based distance often give similar results, but they can differ when genes have different variances or when the data contain outliers. A useful heuristic: if the goal is to identify genes with similar patterns regardless of magnitude, use Pearson correlation; if the goal is to identify genes with similar absolute expression levels, use Euclidean distance on Z-scored data.

Tools: pheatmap, ComplexHeatmap, Morpheus

Several tools are available for heatmap generation, each with its own strengths.

pheatmap (R package) is a simple, reliable choice for standard heatmaps. It supports hierarchical clustering, row and column annotations, and customization of color scales. Its syntax is straightforward, making it ideal for quick exploratory analysis.

ComplexHeatmap (R package) is a more powerful and flexible option. It supports complex annotations, multiple heatmaps arranged in a single figure, and fine-grained control over every graphical element. It is the tool of choice for publication-quality figures that require sample annotations, gene annotations, and significance markers.

Morpheus (web-based, from the Broad Institute) is an interactive tool that allows users to upload a matrix, perform clustering, and explore the heatmap dynamically. It is useful for sharing data with collaborators who may not use R.

For Python users, seaborn's clustermap function provides a convenient interface for clustered heatmaps, and plotly offers interactive options.

The choice of tool is less important than the underlying decisions about normalization, transformation, and clustering. A heatmap generated with any of these tools will be valid if the data preprocessing is correct.

Interpreting Heatmaps: Patterns and Biological Meaning

Reading color scales and dendrograms

The color scale is the most immediately visible element of a heatmap. A typical scale ranges from blue (low expression) through white (mean expression) to red (high expression), with the specific colors chosen for accessibility and contrast. The scale is usually centered at zero for Z-scored data, so blue indicates below-mean expression and red indicates above-mean expression for that gene across the samples.

The dendrogram on the left (for genes) and top (for samples) encodes the similarity structure. The height at which two branches merge indicates their dissimilarity: the higher the merge, the more dissimilar the clusters. A long vertical branch separating two major groups of samples indicates a strong transcriptional difference between those groups, which is often the primary biological finding of the experiment.

When reading a heatmap, the first step is to examine the sample dendrogram. Do the samples cluster by treatment group, time point, or genotype? If biological replicates do not cluster together, this is a red flag for batch effects or sample mislabeling. The second step is to examine the gene dendrogram. Are there clear modules of genes that are coordinately up- or down-regulated across conditions?

Identifying gene modules

A gene module is a cluster of genes that share a similar expression pattern across samples. These modules often represent co-regulated genes—for example, genes controlled by the same transcription factor, genes in the same metabolic pathway, or genes encoding subunits of the same protein complex.

To identify modules, look for blocks of color in the heatmap: a contiguous block of red across a set of samples indicates a group of genes that are coordinately upregulated in those samples. The boundaries of these blocks are defined by the dendrogram branches. A well-defined module will have a clear branch in the gene dendrogram and a consistent color pattern across the samples.

For example, in a time-course experiment of immune activation, one module might contain genes that peak at 2 hours post-stimulation (early response genes such as FOS, JUN, and NFKBIA), while another module contains genes that peak at 24 hours (late response genes such as CXCL10 and IRF7). The heatmap makes these temporal programs immediately apparent.

Linking clusters to functional enrichment

Once gene modules are identified, the next step is to ask what biological functions they represent. This is done through functional enrichment analysis, which tests whether a module is enriched for genes annotated with a particular Gene Ontology (GO) term, KEGG pathway, or other gene set. Tools such as DAVID, Enrichr, or clusterProfiler (R package) perform this analysis. The output is a list of enriched terms with associated p-values, which can be used to assign a biological label to each module.

For example, a module that is enriched for genes in the "inflammatory response" GO term and the "TNF signaling pathway" KEGG pathway is likely to be a pro-inflammatory module. This functional annotation transforms the heatmap from a purely descriptive figure into a hypothesis-generating tool. The connection between heatmap modules and functional enrichment is central to Gene Ontology Pathway Enrichment analyses, and the two approaches are often combined in a single figure.

Enhancing Heatmaps with Annotations

Sample annotations (treatment, time, batch)

A heatmap without annotations is difficult to interpret. Sample annotations are colored bars placed above the heatmap that indicate metadata for each column. Common annotations include treatment group (e.g., drug vs. vehicle), time point, genotype, sex, and batch. These annotations allow the reader to immediately see whether the sample clustering corresponds to the experimental design.

In R, the pheatmap and ComplexHeatmap packages support annotation data frames. For example, in pheatmap, the annotation_col argument takes a data frame where each row corresponds to a sample and each column is a metadata variable. The colors for each annotation level are specified with annotation_colors.

Batch annotations are particularly important. If samples from different sequencing batches are included, annotating the batch can reveal whether the sample clustering is driven by biology or by technical artifacts. If batch effects are present, they should be addressed with methods such as Combat Batch Effect Removal before generating the final heatmap.

Gene annotations (pathways, GO terms)

Gene annotations are colored bars placed to the right of the heatmap that indicate properties of each gene. These can include membership in a pathway, a GO term, a protein domain, or a classification such as "transcription factor" or "cytokine." Gene annotations help the reader identify whether specific functional categories are enriched in particular modules.

For example, if the heatmap includes genes from multiple signaling pathways, a gene annotation bar indicating pathway membership (e.g., "NF-κB", "JAK-STAT", "Wnt") will immediately reveal whether the modules correspond to distinct pathways. This is particularly useful when the heatmap contains hundreds of genes and the dendrogram alone is insufficient to convey the functional structure.

Highlighting key genes

In some cases, the goal is to draw attention to specific genes of interest. This can be achieved by labeling the rows with gene names (rather than only showing the dendrogram), or by adding a side bar that highlights the genes of interest in a distinct color. For example, in a heatmap of genes regulated by a transcription factor, the known direct targets could be highlighted, allowing the reader to see whether they cluster together and whether their expression patterns match the expected kinetics.

In ComplexHeatmap, the rowAnnotation function can be used to add a custom annotation that highlights specific genes. Alternatively, the row labels can be colored or bolded for the genes of interest.

Common Pitfalls and Best Practices

Overinterpretation of clustering

The most common error in heatmap interpretation is overinterpreting the dendrogram. Hierarchical clustering will always produce a tree, even for random data. The fact that two genes are adjacent in the dendrogram does not necessarily mean they are co-regulated; it may simply reflect noise. The robustness of a cluster should be assessed, for example by bootstrapping (resampling the data and checking whether the cluster is reproduced) or by comparing the clustering to an independent dataset.

Similarly, the sample dendrogram should not be overinterpreted. If the sample clustering does not match the experimental groups, this may indicate a batch effect, but it may also indicate that the biological signal is weak relative to the technical noise. Before concluding that the biology is unexpected, check the normalization and transformation steps—a failure to properly normalize can easily produce spurious clustering.

Improper normalization

Generating a heatmap from raw counts or from improperly normalized data is a frequent and serious error. Raw counts are not comparable between samples with different sequencing depths, and RPKM/FPKM values are not appropriate for between-sample comparisons. Using these values directly in a heatmap will produce a figure that primarily reflects sequencing depth and library composition rather than biology.

The correct approach is to use variance-stabilized or regularized log-transformed counts from DESeq2, or TMM-normalized log2-cpm values from edgeR. These transformations account for the mean-variance relationship and make the data suitable for clustering.

Color scale misrepresentation

The choice of color scale can dramatically affect the interpretation of a heatmap. A scale that is too narrow will make noise appear as signal; a scale that is too wide will obscure real differences. For Z-scored data, a common range is -2 to +2 (i.e., two standard deviations from the mean). Values beyond this range are clipped to the extreme colors, which prevents a few extreme genes from dominating the color scale.

It is also important to choose a color scale that is perceptually uniform and accessible to color-blind readers. The "RdBu" (red-blue) or "viridis" color scales are good choices. Avoid rainbow color scales, which are not perceptually uniform and can create false boundaries.

Reproducibility and reporting

A heatmap is only useful if it can be reproduced. The code used to generate the heatmap should be saved and version-controlled, and the exact parameters (normalization method, transformation, distance metric, clustering method, color scale limits) should be reported in the methods section of any publication. The raw data and the processed matrix should be deposited in a public repository.

A common failure is to generate a heatmap interactively and then lose the code. This makes it impossible to regenerate the figure if a reviewer requests a change (e.g., a different color scale or a different set of genes). Best practice is to write a script that generates the heatmap from the processed data, so that the figure can be regenerated with a single command.

Frequently Asked Questions

What is a heatmap of differentially expressed genes?

A heatmap of differentially expressed genes is a visualization in which rows represent genes that are significantly differentially expressed between conditions, columns represent biological samples, and colors represent the expression level (or Z-score) of each gene in each sample. It is used to visualize the overall pattern of gene expression changes, identify co-regulated gene clusters, and assess the reproducibility of biological replicates.

How do I choose genes for a differential expression heatmap?

Genes are selected by applying statistical thresholds: an adjusted p-value (FDR) below a cutoff (typically 0.05) and an absolute log2 fold change above a cutoff (typically 1). After applying these filters, you may further restrict the list to the top N genes ranked by significance, or intersect the list with a curated gene set of interest. The goal is to include enough genes to show the biological pattern but not so many that the heatmap becomes unreadable.

What normalization should I use before making a heatmap?

For RNA-seq data, use a variance-stabilizing transformation (VST) or regularized log (rlog) transformation from DESeq2, or TMM-normalized log2-counts-per-million from edgeR. These transformations account for the mean-variance relationship and make the data suitable for clustering. Do not use raw counts or RPKM/FPKM values for between-sample comparisons. After transformation, Z-score each gene across samples so that the mean is 0 and the standard deviation is 1.

Why do my heatmap clusters not match my expectations?

There are several possible reasons. First, the normalization or transformation may be incorrect, leading to spurious clustering driven by technical artifacts. Second, the biological signal may be weak relative to noise, especially if the sample size is small. Third, the distance metric may be inappropriate—for example, using Euclidean distance without Z-scoring will make the clustering reflect absolute expression levels rather than patterns. Fourth, there may be batch effects that are not accounted for. Check each of these in order.

What is the difference between Euclidean distance and Pearson correlation in clustering?

Euclidean distance measures the absolute difference in expression values between two genes (or samples) across all dimensions. It is sensitive to the magnitude of expression, so Z-scoring is required to make it reflect pattern similarity. Pearson correlation measures the linear association between two expression profiles and is invariant to linear transformations, so it naturally reflects pattern similarity regardless of baseline expression. In practice, Euclidean distance on Z-scored data and 1 - Pearson correlation often give similar results, but they can differ when the data contain outliers or when genes have different variances.

How can I add annotations to my heatmap?

In R, the pheatmap package supports the annotation_col and annotation_row arguments, which take data frames of sample and gene metadata, respectively. The ComplexHeatmap package provides more flexible annotation options, including the ability to add multiple annotation bars, custom colors, and labels. For web-based tools, Morpheus allows annotations to be added directly to the heatmap interface.

What are common mistakes when making a heatmap?

Common mistakes include: using raw or improperly normalized counts; failing to Z-score the data before clustering; using an inappropriate distance metric; choosing a color scale that obscures or exaggerates differences; overinterpreting the dendrogram; ignoring batch effects; and failing to report the exact parameters used to generate the figure. Each of these can lead to a misleading or uninterpretable heatmap.

Key Takeaways

A heatmap of differentially expressed genes is a powerful visualization tool that reveals co-regulated gene modules and sample relationships, but its reliability depends entirely on proper data preprocessing.
Always normalize RNA-seq data with DESeq2's variance-stabilizing transformation or edgeR's TMM method before generating a heatmap, and Z-score each gene across samples to make patterns comparable.
Select genes for the heatmap using combined statistical thresholds (adjusted p-value < 0.05 and |log2 fold change| ≥ 1), and consider further filtering for biological relevance.
Choose the distance metric deliberately: Euclidean distance on Z-scored data for pattern similarity, or Pearson correlation for invariance to baseline expression.
Use annotations for samples (treatment, time, batch) and genes (pathways, GO terms) to make the heatmap interpretable without reference to external tables.
Avoid common pitfalls: do not overinterpret clustering, do not use improper normalization, choose a perceptually uniform color scale, and always save the code for reproducibility.
Link heatmap modules to functional enrichment using tools like clusterProfiler to transform the visualization from a descriptive figure into a hypothesis-generating analysis.