How to Build a Transcriptome Heatmap From Your Own RNA-Seq Data

Oak
rna-seq transcriptome heatmap deseq2 visualization
A tidal flat covered in a grid of flat coral cells glowing at different intensities, reflected in shallow water at dusk.

By the end of this you will have two figures: a sample-by-sample correlation heatmap that tells you whether your data is structured the way you think it is, and a gene-by-sample expression heatmap with annotated columns, clustered rows, and a color scale that does not lie about effect size. You need FASTQs from at least four RNA-seq libraries (more is better, and for a single person a time course of six to twelve draws is where this gets interesting), a reference transcriptome and genome FASTA for decoy-aware indexing, roughly 16 GB of RAM and 8 cores, R 4.3+, and the packages tximport, DESeq2, ComplexHeatmap, circlize, and limma. Everything below assumes Illumina paired-end polyA or globin-depleted whole blood at 30–50M read pairs per sample. A heatmap is a picture of a matrix, so almost all the work is deciding which matrix.

1. Quantify with salmon, decoy-aware

Alignment-free quantification is the right default here. It runs in minutes per sample, gives you transcript-level estimates with bootstraps, and you do not need a BAM unless you plan to look at splicing or variants.

Build the index once:

# gencode v45 human
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) | cut -d " " -f 1 \
  | sed 's/>//g' > decoys.txt
cat gencode.v45.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index -t gentrome.fa.gz -d decoys.txt -p 16 -i salmon_idx_v45 -k 31 --gencode

The decoy genome matters. Without it, reads from unannotated or intronic regions get force-assigned to whatever transcript they partially match, and you get a handful of genes with fake expression that then dominate your variable-gene selection in step 4.

Quantify:

for s in $(cat samples.txt); do
  salmon quant -i salmon_idx_v45 -l A \
    -1 fastq/${s}_R1.fastq.gz -2 fastq/${s}_R2.fastq.gz \
    --gcBias --seqBias --posBias \
    --numBootstraps 30 \
    --validateMappings -p 16 \
    -o quant/${s}
done

Check quant/*/logs/salmon_quant.log for mapping rate. Below 70% on human polyA data means something is wrong: wrong index version, heavy rRNA or globin carryover, or adapter contamination. Do not proceed with a sample under 60% and expect it to cluster with the others. --gcBias costs a few minutes and is worth it if libraries were prepped in different batches.

2. Import to gene level

library(tximport); library(DESeq2)
library(AnnotationDbi); library(org.Hs.eg.db)

samples <- read.csv("samples.csv")         # sample_id, draw_date, batch, timepoint
files <- file.path("quant", samples$sample_id, "quant.sf")
names(files) <- samples$sample_id

tx2gene <- read.delim("tx2gene_gencode_v45.tsv")   # TXNAME, GENEID
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ batch + timepoint)

countsFromAbundance = "lengthScaledTPM" gives you counts that are already corrected for transcript-length shifts across samples, which matters when isoform usage changes. Use plain "no" if you want DESeq2 to handle length via its offset matrix; either is defensible, but pick one and note it.

Filter before anything else:

keep <- rowSums(counts(dds) >= 10) >= max(3, floor(0.25 * ncol(dds)))
dds <- dds[keep, ]
nrow(dds)   # expect ~13,000-16,000 for whole blood, ~15,000-18,000 for a solid tissue

If you are left with 25,000 genes you filtered too loosely and half of them are lncRNAs with counts of 11, 0, 0, 14. Those genes have enormous variance after scaling and they will win the top-variable-gene competition in step 4.

3. Variance-stabilize, do not log-CPM by hand

log2(CPM + 1) inflates the variance of low-count genes, and on a z-scored heatmap that shows up as a band of screaming red-and-blue rows that are entirely counting noise. Use vst:

vsd <- vst(dds, blind = FALSE, nsub = 1000)
mat <- assay(vsd)          # genes x samples, roughly log2 scale, homoscedastic

blind = FALSE tells the dispersion fit to use your design, which is what you want for visualization of a known contrast. Use blind = TRUE only when the heatmap is a QC artifact and you want it uninformed by the design. If you have fewer than about 8 samples, rlog(dds, blind = FALSE) is more stable than vst and slower; at 30 samples vst is fine and takes a second.

If batch is real and orthogonal to your variable of interest, remove it for display only:

library(limma)
mat_adj <- removeBatchEffect(mat, batch = vsd$batch,
                             design = model.matrix(~ vsd$timepoint))

Keep the unadjusted matrix for any statistics. removeBatchEffect is a plotting tool. If you fit models on its output, your standard errors are wrong.

4. Look at sample-sample structure before you look at genes

This is the figure people skip and then regret. Compute the correlation between samples on the top few thousand variable genes and cluster it.

library(ComplexHeatmap); library(circlize)

rv <- rowVars(mat_adj)
top <- order(rv, decreasing = TRUE)[1:2000]
cm <- cor(mat_adj[top, ], method = "spearman")

Heatmap(cm,
        name = "Spearman",
        col = colorRamp2(c(0.80, 0.90, 1.00), c("#2166AC", "#F7F7F7", "#B2182B")),
        clustering_distance_rows = "pearson",
        clustering_method_rows = "ward.D2",
        top_annotation = HeatmapAnnotation(batch = vsd$batch,
                                           timepoint = vsd$timepoint))

For longitudinal draws from one person on the same platform, expect Spearman correlations of 0.95–0.99 between adjacent timepoints. If two samples sit at 0.85 while everything else is at 0.97, you have a technical outlier, not a biological state change. Look at library size, duplication rate, and RNA integrity before you interpret it. The color scale limits here matter: if you let the default stretch from 0 to 1 when all your values are between 0.94 and 0.99, the plot is one flat color.

5. Choose the genes you display, and be honest about how

Three defensible strategies, in the order we reach for them.

Top variable genes. Fast, unbiased, and the right choice when you do not have a hypothesis:

sel <- order(rowVars(mat_adj), decreasing = TRUE)[1:50]

Fifty rows fits on a screen with readable labels. Two thousand rows is a texture, not a figure, and is only useful for showing gross block structure.

A curated panel. If you care about a specific process (interferon response, mitochondrial oxidative phosphorylation, circadian output), pull the gene set from MSigDB Hallmark or Reactome and plot those rows. The heatmap then answers a question instead of describing the data. This is how multi-omics papers use heatmaps: a module of genes defined by an external screen, shown across conditions, with the omics layers stacked so you can see whether the transcript-level pattern is echoed elsewhere.1

Differentially expressed genes. Legitimate, with one rule: do not select genes by a test across your conditions and then present the resulting block structure as evidence that those conditions differ. The blocks are guaranteed by construction. If you select on DE, say so in the caption and let the heatmap show the shape and consistency of the change, not its existence.

res <- results(DESeq(dds), contrast = c("timepoint", "T2", "T0"))
sel <- rownames(res)[which(res$padj < 0.05 & abs(res$log2FoldChange) > 1)]

6. Scale rows, and set the color limits by hand

Genes differ in mean expression by four orders of magnitude. Unscaled, your heatmap is a picture of which genes are abundant, which you already knew. Z-score each row across samples:

z <- t(scale(t(mat_adj[sel, ])))

Now fix the limits. The default behavior of most heatmap functions is to map the extremes of the data to the extremes of the palette, so a single outlier gene compresses everything else to white:

col_fun <- colorRamp2(c(-2, 0, 2), c("#2166AC", "#FFFFFF", "#B2182B"))

Symmetric limits at ±2 SD, a diverging palette, white at zero. Blue-white-red is conventional; if colorblind accessibility matters use colorRamp2(c(-2,0,2), c("#0571B0","#F7F7F7","#E66101")) or a viridis-family sequential map if your data is not naturally centered. Never use a rainbow palette: the perceptual jumps in it create edges that are not in the data.

Two caveats on z-scoring. With fewer than about five samples, the per-gene SD is estimated so poorly that the z-scores are close to arbitrary. And after scaling, a gene that moves from 2.0 to 2.4 on the vst scale (a 1.3x change) looks identical to one that moves from 4 to 12. Put the un-scaled values somewhere, either as a second annotation column showing row means or in the supplementary table.

7. Draw the real heatmap

ha <- HeatmapAnnotation(
  timepoint = vsd$timepoint,
  batch     = vsd$batch,
  hsCRP     = anno_barplot(samples$hscrp_mg_L),
  col = list(timepoint = c(T0="#DDDDDD", T1="#999999", T2="#333333"),
             batch = c(A="#8DD3C7", B="#FB8072"))
)

ht <- Heatmap(
  z,
  name = "z-score",
  col = col_fun,
  top_annotation = ha,
  cluster_rows = TRUE,
  clustering_distance_rows = "pearson",
  clustering_method_rows = "ward.D2",
  cluster_columns = FALSE,            # time is ordered; do not let clustering scramble it
  row_names_gp = gpar(fontsize = 7),
  row_split = 3,
  border = TRUE,
  heatmap_legend_param = list(at = c(-2, -1, 0, 1, 2))
)

pdf("transcriptome_heatmap.pdf", width = 7, height = 9)
draw(ht, merge_legend = TRUE)
dev.off()

The three decisions that carry the most weight:

Column clustering off when columns have an intrinsic order (time, dose, developmental stage). Clustering them hides the ordering and invites you to see groups that are really just noise in the dendrogram. Turn it on when samples are exchangeable and you want to discover grouping.

Correlation distance for rows, not Euclidean. On z-scored rows the two are monotonically related, so it matters less than people claim, but Pearson distance makes the intent explicit: group genes with the same shape over samples.

Ward.D2 linkage for compact, similar-sized clusters. Average linkage gives you one giant cluster and several singletons more often than not. Whatever you choose, note it in the caption, because the dendrogram shape is a function of linkage and readers will over-interpret it. Then use row_split to cut the tree at a fixed number of clusters, and inspect what is in each with row_order(ht).

Export to PDF, not PNG. A heatmap with 50 rows and gene labels needs vector text.

8. Say what the clusters are

A heatmap with three unlabeled row blocks is not a result. Pull the genes from each block and run an over-representation test:

library(clusterProfiler)
ro <- row_order(ht)
gl <- lapply(ro, function(i) rownames(z)[i])
ego <- compareCluster(gl, fun = "enrichGO", OrgDb = org.Hs.eg.db,
                      keyType = "SYMBOL", ont = "BP",
                      pAdjustMethod = "BH", qvalueCutoff = 0.05)

Then annotate the blocks in the figure with the top term per cluster. Expect blood to give you neutrophil degranulation, interferon signaling, and translation as the three loudest modules almost regardless of the question. Also expect that a pathway-level signal in transcript data is a hypothesis, not a finding. Published work in this space pairs the expression signature with an orthogonal measurement before claiming a mechanism.2 For your own data the orthogonal layer is usually proteomics or a targeted assay on the same draw.

If you want external context for what a given expression pattern looks like in other cohorts, public transcriptome resources with harmonized processing let you place your matrix against a reference rather than reading it in isolation.3 Reprocess the reference FASTQs through your own pipeline if you can. Cross-pipeline comparison of quantifications is one of the most common sources of fake batch effects.

None of this is a clinical test. RNA expression in whole blood shifts with sleep, recent exercise, a cold coming on, and time of draw, and no expression pattern here establishes a diagnosis. If something in the data concerns you, bring it and your standard labs to a physician.

Common problems

One gene occupies half the dynamic range. In whole blood without globin depletion, HBB, HBA1, and HBA2 can consume a majority of reads and crush your effective depth for everything else. Use globin-depleted or PAXgene-with-GLOBINclear prep. If you already have the data, drop those genes before variance ranking and note the reduced sensitivity.

The heatmap splits perfectly by batch. Check whether batch is confounded with your variable of interest using table(samples$batch, samples$timepoint). If every T0 draw was sequenced in batch A, no statistical correction saves you. Randomize draws across batches at the point of collection, or process the whole time course together.

Every row looks like the same wave. In blood, most of the top-variable genes track cell composition, not per-cell regulation. Neutrophil fraction alone can move thousands of genes. Get a CBC with differential on the same draw and either include it as an annotation track or regress it out with removeBatchEffect(mat, covariates = cbind(neut_pct, lymph_pct)) for display. Deconvolution (CIBERSORTx, or a signature-matrix method of your choice) is the alternative when you have no CBC.

Rows with near-zero variance after scaling. t(scale(t(x))) produces NaN for any row with zero SD, and Heatmap will error or drop them silently. Filter with z <- z[complete.cases(z), ] and check how many you lost.

Dendrogram changes every time you rerun. You are clustering on a different top-variable-gene set because of a nondeterministic step upstream, or you changed the filter. Set set.seed(), pin your gene selection to a saved vector, and version the sel object with the figure.

Fold changes look smaller than in the DE table. vst shrinks low-count differences toward the mean by design, so the heatmap is a conservative view. That is the correct behavior. If you want the raw effect, plot log2FoldChange from results() as a side annotation next to the z-scored matrix.

Oak builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

  1. Chunze Xie, Jingjing Long, Tiansi Luo, et al. Integrated GWAS Candidate Prioritization, Moso Bamboo Transcriptome Screening and Diverse Bamboo Shoot Multi-Omics Reveal a Multi-Layer Framework for Bamboo Property Formation. Horticulturae, 2026. https://doi.org/10.3390/horticulturae12081030 ↩

  2. Yuhan Su, Haifeng Huang, Hui Qin. Exploring the molecular mechanism of danshensu targeting PANoptosis therapy for stroke based on multi omics and experimental verification. Biochemical and Biophysical Research Communications, 2026. https://doi.org/10.1016/j.bbrc.2026.154441 ↩

  3. Musun Park, Su-Jin Baek, Eun-Hye Seo, et al. KORE-Map 1.2: Korean medicine Omics Resource Extension Map on transcriptome data of Qi-related herbal medicine. Scientific Data, 2026. https://doi.org/10.1038/s41597-026-07937-2 ↩