How to Visualize Your Own RNA-Seq Data

Oak
rna-seq visualization deseq2 r bioinformatics
A feathered reptile on a specimen stand, its wing feathers glowing in ordered bands of violet, amber and teal against black.

By the end of this guide you will have a small set of plots built from your own transcriptome: a principal components plot and sample-distance heatmap that tell you whether your timepoints cluster by biology or by batch, per-gene trajectory plots over time, an MA plot and volcano plot from a model fit across timepoints, a heatmap of the genes that move most, and a pathway-level summary. You need per-sample transcript quantifications (Salmon or kallisto output, one quant.sf or abundance.h5 per sample), a sample metadata table with collection dates and any covariates you tracked, R 4.3 or later with tximport, DESeq2, ggplot2, ComplexHeatmap, and fgsea, and ideally the aligned BAMs plus IGV for spot checks. If you only have a counts matrix from a provider, most of this still works, but you lose the effective-length correction that makes cross-gene comparison meaningful. RNA-seq measures relative transcript abundance across the whole transcriptome in one assay, which is what makes these plots possible; a qPCR panel gives you absolute-ish quantification of a handful of targets instead, and the two disagree in predictable ways at low expression 1.

1. Assemble the quantification files and a metadata table

Start by getting every sample into the same reference and the same tool version. Mixing a GENCODE v39 index with a v44 index across timepoints will introduce transcript-level differences that look like biology. If you are quantifying yourself, run Salmon with automatic library type detection, GC bias correction, and bootstraps so you have a per-sample uncertainty estimate:

salmon index -t gencode.v44.transcripts.fa.gz -i salmon_v44_idx \
  --gencode -k 31 -p 16

for s in T01 T02 T03 T04 T05 T06; do
  salmon quant -i salmon_v44_idx -l A \
    -1 fastq/${s}_R1.fastq.gz -2 fastq/${s}_R2.fastq.gz \
    --validateMappings --gcBias --seqBias --numBootstraps 30 \
    -p 16 -o quants/${s}
done

Keep the logs/salmon_quant.log from each run. The mapping rate line in that file is the first number to look at: for poly-A selected human blood RNA, 70 to 85 percent is typical, and anything under about 60 percent usually means contamination, adapter carryover, or the wrong index. Write a metadata CSV alongside it with one row per sample and columns for sample ID, collection date, time since a reference date in days, RNA integrity number if your lab reported it, library prep batch, and sequencing run. Batch and run are the two covariates you will most often need to regress out before plotting.

2. Import to gene level and choose a transformation

Transcript-level estimates are noisy for individual isoforms, so summarize to genes for most visualization. tximport does the summarization and, importantly, carries the average transcript length per gene per sample into an offset that DESeq2 uses, which corrects for isoform switching changing a gene’s effective length between timepoints.

library(tximport); library(DESeq2)
meta  <- read.csv("metadata.csv", stringsAsFactors = FALSE)
files <- file.path("quants", meta$sample, "quant.sf"); names(files) <- meta$sample
tx2gene <- read.delim("tx2gene_v44.tsv")          # tx_id, gene_id, gene_symbol
txi <- tximport(files, type = "salmon", tx2gene = tx2gene[,1:2],
                countsFromAbundance = "no")

dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ batch + day)
keep <- rowSums(counts(dds) >= 10) >= 3           # gene seen in >= 3 samples
dds  <- dds[keep, ]
vsd  <- vst(dds, blind = FALSE)                   # blind=FALSE: design-aware

Use variance-stabilized values (vst) or regularized log (rlog) for every distance-based or clustering plot, never raw counts and never plain log2 TPM with a pseudocount. The reason is heteroskedasticity: in count data the variance grows with the mean, so log-transformed low-count genes become the noisiest rows in your matrix and dominate any Euclidean distance or PCA. With fewer than about 30 samples rlog is the better-behaved choice and vst is faster; at a personal-profile scale of six to twenty timepoints either is fine.

For plots where you want an interpretable unit on the y-axis, use TPM (transcripts per million) rather than RPKM or FPKM. TPM divides by transcript length first and then normalizes so every sample’s values sum to exactly one million, which means a given TPM has the same meaning in every sample. RPKM normalizes by total reads first, so its per-sample sum depends on the shape of the expression distribution and the same number in two samples does not represent the same molar fraction. Neither corrects for composition effects, so use DESeq2 median-of-ratios size factors when you are comparing across samples formally.

3. Look at global structure before anything else

The first two plots answer a single question: do your samples separate by something you care about, or by how they were processed? Run PCA on the variance-stabilized matrix and a sample-to-sample distance heatmap on the same values.

library(ggplot2); library(ComplexHeatmap); library(circlize)

pcaData <- plotPCA(vsd, intgroup = c("batch", "day"), ntop = 500,
                   returnData = TRUE)
pv <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color = day, shape = batch)) +
  geom_point(size = 4) +
  geom_text(aes(label = name), vjust = -1.2, size = 3, show.legend = FALSE) +
  labs(x = paste0("PC1: ", pv[1], "% variance"),
       y = paste0("PC2: ", pv[2], "% variance")) +
  coord_fixed() + theme_bw()

sampleDists <- dist(t(assay(vsd)))
Heatmap(as.matrix(sampleDists), name = "dist",
        col = colorRamp2(c(0, max(sampleDists)), c("#2166AC", "white")),
        top_annotation = HeatmapAnnotation(batch = vsd$batch))

Note that plotPCA uses only the top 500 most variable genes by default. Rerun with ntop = 2000 and ntop = nrow(vsd); if the structure flips, your PC1 is being driven by a small set of genes and you should find out which ones with prcomp rotations before interpreting anything. In whole blood, PC1 very often tracks the neutrophil-to-lymphocyte composition of that draw rather than any cell-intrinsic change, which is a real biological signal but not the one people usually think they are seeing.

If batch separates your samples cleanly, remove it for display only, and keep the covariate in the model for testing:

mat <- assay(vsd)
mat <- limma::removeBatchEffect(mat, batch = vsd$batch,
                                design = model.matrix(~ vsd$day))

Do not feed batch-corrected values back into DESeq2. Correct once for plotting, and let the design formula handle it for inference.

For single-cell data the analogous first look is a t-SNE or UMAP embedding, and the same caution applies with more force: distances between clusters in those embeddings are not quantitative, and the apparent tightness of a cluster is a function of perplexity or n_neighbors rather than of biology 2.

4. Plot trajectories for individual genes

With a longitudinal profile of one person, the most informative plot is often the least sophisticated one: normalized expression for a single gene across every timepoint, with the points connected in collection order. This is where you see whether something moved.

plotGene <- function(dds, symbol, tx2gene) {
  gid <- unique(tx2gene$gene_id[tx2gene$gene_symbol == symbol])
  d <- plotCounts(dds, gene = gid, intgroup = "day",
                  normalized = TRUE, returnData = TRUE)
  ggplot(d, aes(day, count)) +
    geom_line(color = "grey50") + geom_point(size = 3) +
    scale_y_log10() +
    labs(title = symbol, y = "normalized count (log10)", x = "day") +
    theme_bw()
}
plotGene(dds, "HBB", tx2gene)

Plot on a log scale by default. A gene going from 40 to 120 normalized counts and one going from 4,000 to 12,000 are the same fold change, and a linear axis makes the second look dramatic and hides the first. When you want to compare several genes on one panel, convert each to a z-score across your own timepoints, which is the right within-person unit: it asks how unusual today is relative to your personal distribution, not relative to a population reference.

Treat single-gene trajectories for low-expressed transcripts with suspicion. Benchmarking of RNA-seq workflows against whole-transcriptome RT-qPCR showed that agreement between pipelines and with qPCR degrades substantially for low-abundance genes, so a transcript sitting at 15 counts that appears to double is within the noise of quantification rather than a finding 3.

5. Fit a model across time and draw MA and volcano plots

With one person and one sample per timepoint, there are no replicates in the usual sense, so the design has to borrow information across timepoints. Treating day as a continuous covariate (~ batch + day) tests for monotone drift. Treating a binary state you recorded (before and after a dietary change, say, with several samples in each period) as a factor tests for a shift between periods. Either way, DESeq2 estimates dispersion by pooling across genes, which is what makes inference possible at all here.

dds <- DESeq(dds)
res <- results(dds, name = "day", alpha = 0.05)
res <- lfcShrink(dds, coef = "day", type = "apeglm", res = res)

plotMA(res, ylim = c(-2, 2), alpha = 0.05)

Always apply shrinkage before plotting log fold changes. Without it, the left side of the MA plot fills with low-count genes showing implausible 6-fold swings, and the plot becomes a picture of your sequencing depth. The MA plot is the diagnostic to look at first: a well-behaved fit is a cloud centered on zero across the whole expression range, and a systematic tilt at high expression means normalization has failed, usually because a few very abundant transcripts, globin or immunoglobulin in blood, dominate the library.

The volcano plot is the summary view:

library(ggrepel)
df <- as.data.frame(res); df$symbol <- tx2gene$gene_symbol[match(rownames(df), tx2gene$gene_id)]
df$sig <- !is.na(df$padj) & df$padj < 0.05 & abs(df$log2FoldChange) > 0.5

ggplot(df, aes(log2FoldChange, -log10(pvalue), color = sig)) +
  geom_point(alpha = 0.5, size = 1) +
  scale_color_manual(values = c("grey70", "#B2182B")) +
  geom_vline(xintercept = c(-0.5, 0.5), linetype = 2) +
  ggrepel::geom_text_repel(data = subset(df, sig)[1:25, ],
                           aes(label = symbol), size = 3, max.overlaps = 40) +
  theme_bw()

Plot -log10(pvalue) on the y-axis but color by adjusted p-value. Plotting adjusted values directly produces a flat shelf of ties from the Benjamini-Hochberg step that obscures the shape of the distribution.

6. Heatmap the genes that move, with the rows scaled

A heatmap is only as good as its row selection and its scaling. Select rows by the model result or by variance across timepoints, scale each row to a z-score so that color encodes change rather than absolute abundance, and cluster rows but not columns when columns are ordered by time.

sel <- head(order(res$padj), 60)
m <- t(scale(t(assay(vsd)[sel, order(vsd$day)])))
rownames(m) <- tx2gene$gene_symbol[match(rownames(m), tx2gene$gene_id)]

Heatmap(m, name = "z", cluster_columns = FALSE,
        col = colorRamp2(c(-2, 0, 2), c("#2166AC", "white", "#B2182B")),
        row_names_gp = gpar(fontsize = 8),
        top_annotation = HeatmapAnnotation(day = anno_points(sort(vsd$day))))

Cap the color scale at roughly ±2 or ±2.5 standard deviations. Letting the scale run to the data maximum means a single outlier gene washes the rest of the map to white. If you cluster columns as a sanity check and the timepoints do not group in time order, that is informative, and worth looking into before you interpret any row.

7. Summarize at the pathway level

Individual genes in an n-of-1 series are noisy, and coordinated movement across a gene set is more stable than any single member of it. Use pre-ranked GSEA on the full ranked statistic rather than an overlap test on a thresholded list, because the ranked version does not throw away the genes that sat just outside your cutoff.

library(fgsea); library(msigdbr)
stats <- res$stat; names(stats) <- tx2gene$gene_symbol[match(rownames(res), tx2gene$gene_id)]
stats <- stats[!is.na(stats) & !duplicated(names(stats))]
h <- msigdbr(species = "Homo sapiens", category = "H")
pathways <- split(h$gene_symbol, h$gs_name)

fg <- fgsea(pathways, stats, minSize = 15, maxSize = 500, nPermSimple = 10000)
fg <- fg[order(fg$padj), ]

ggplot(head(fg, 20), aes(NES, reorder(pathway, NES), size = size, color = padj)) +
  geom_point() + geom_vline(xintercept = 0) + theme_bw()

The MSigDB hallmark collection is the right starting point at 50 well-curated sets. Reactome and GO add resolution and a great deal of redundancy, so if you use them, collapse near-duplicate sets with collapsePathways before plotting. A pathway view is also the natural place to join transcript data to the rest of a molecular profile, since proteomic and epigenetic measurements map onto the same gene identifiers, and several integration platforms are built around exactly that join 45.

8. Verify surprising genes in a genome browser

Before you believe any single gene, load the BAM in IGV and look at the coverage. Sort and index first, and set IGV to show junctions.

samtools sort -@ 8 -o T03.sorted.bam T03.bam && samtools index T03.sorted.bam

Three patterns explain most false positives. Coverage piled on one exon with the rest empty usually means a repeat or a mismapped paralog. Coverage spread evenly across an intron means genomic DNA carryover or unspliced pre-mRNA. Coverage restricted to one isoform when the gene has several means the change is a splicing shift rather than a change in total output, which a sashimi plot of junction counts will show directly. Polymorphic gene families deserve extra scepticism: for HLA genes, reads from divergent alleles align poorly to the reference and quantification is biased in ways that persist through the whole pipeline, so apparent expression changes there need personalized reference sequences to be trusted 6. Established best-practice treatments of RNA-seq quantification cover the alignment and reference choices behind these artifacts in more depth 7.

Common problems

Globin transcripts swamping a whole-blood library is the most frequent one. HBB, HBA1, and HBA2 can take up a large share of reads in non-depleted whole blood, which compresses the depth available for everything else and shows up as a compressed dynamic range in your MA plot. If your provider offered globin depletion or PAXgene-based prep, prefer it, and either way exclude the globin genes from the variable-gene selection you feed to PCA.

Interpreting PC1 as a biological phenotype when it is cell composition is the second. Whole blood and PBMC transcriptomes track the fraction of neutrophils, monocytes, and lymphocytes in that draw, and a deconvolution step (CIBERSORTx or similar) run on your TPM matrix gives you estimated fractions you can plot as a column annotation next to any heatmap.

Reading fold changes off unshrunk estimates is third, and it is entirely avoidable with lfcShrink. Related to it is drawing conclusions from genes below roughly 20 to 50 counts, where quantification methods diverge from each other and from qPCR 3.

Finally, a plot that looks abnormal is a measurement, not a finding about your health. Transcript abundance in blood shifts with time of day, fasting state, recent exercise, and any concurrent infection, and none of these plots distinguish among those causes. If something in your data looks clinically relevant, bring the underlying values and the collection context to a physician or a genetic counselor, who can order a validated clinical assay; research-grade RNA-seq is not a diagnostic test.

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. C. A. Fassbinder-Orth. Methods for Quantifying Gene Expression in Ecoimmunology: From qPCR to RNA-Seq. Integrative and Comparative Biology, 2014. https://doi.org/10.1093/icb/icu023 ↩

  2. Yan Wu, Kun Zhang. Tools for the analysis of high-dimensional single-cell RNA sequencing data. Nature Reviews Nephrology, 2020. https://doi.org/10.1038/s41581-020-0262-0 ↩

  3. Celine Everaert, Manuel Luypaert, Jesper L. V. Maag, et al. Benchmarking of RNA-sequencing analysis workflows using whole-transcriptome RT-qPCR expression data. Scientific Reports, 2017. https://doi.org/10.1038/s41598-017-01617-3 ↩ ↩2

  4. Zeyu Zhang, Zuoling Ma, Tian Hua, et al. OmicsCanvas: A multi‐omics platform for integration and visualization of epigenetic regulation. Journal of Integrative Plant Biology, 2026. https://doi.org/10.1111/jipb.70258 ↩

  5. David Márquez-Oller, Andrea Pauli, Jörg Fallmann. BRIDGE: an interactive application for multi-omics data analysis, visualization and integration. Bioinformatics, 2026. https://doi.org/10.1093/bioinformatics/btag558 ↩

  6. Vitor R. C. Aguiar, Erick C. Castelli, Richard M. Single, et al. Comparison between qPCR and RNA-seq reveals challenges of quantifying HLA expression. Immunogenetics, 2023. https://doi.org/10.1007/s00251-023-01296-7 ↩

  7. Michele Araújo Pereira, Eddie Luidy Imada, Rafael Lucas Muniz Guedes. RNA‐seq: Applications and Best Practices. Applications of RNA-Seq and Omics Strategies - From Microorganisms to Human Health, 2017. https://doi.org/10.5772/intechopen.69250 ↩