How to Plot Your Own RNA-Seq Data: PCA, Volcano, and Trajectory Graphs

Oak
rna-seq pca visualization deseq2 python
A laboratory instrument holding glowing sample cartridges projects a floating cloud of light points that separate into two clusters.

By the end of this guide you will have a gene-by-sample count matrix, a variance-stabilized expression matrix, and five plots that answer distinct questions: a principal component analysis (PCA) scatter showing how your samples relate to each other, a sample-distance heatmap that confirms or contradicts the PCA, an MA plot that shows whether your differential expression results are well-behaved, a volcano plot that ranks genes by effect and confidence, and per-gene trajectory plots for longitudinal data. You need FASTQ files from your sequencing provider (paired-end, typically 2 × 100 bp or 2 × 150 bp), roughly 16 GB of RAM and 100 GB of disk, R 4.3+ with DESeq2, and Python 3.10+ with pandas, scikit-learn, and matplotlib. If your provider already gave you a counts matrix, skip to step 3, but read step 2 anyway, because most bad RNA-seq plots are bad because nobody looked at the QC first.

1. Quantify the FASTQs into a counts matrix

We use Salmon for quantification rather than aligning to the genome with STAR. Salmon maps reads to a transcriptome index with selective alignment, runs about ten times faster, needs a fraction of the memory (roughly 4 GB versus the ~32 GB STAR wants for a human genome index), and gives you transcript-level estimates with bootstrap replicates for free. The tradeoff is that you get no genomic alignments, so you cannot inspect coverage in IGV, call variants from the RNA, or detect novel splice junctions. If any of those matter to you, run STAR as well and keep the BAMs.

Build a decoy-aware index against GENCODE. The decoy sequences (the whole genome) prevent reads from genomic or intronic origin from being force-fit to transcripts, which is the single most common source of inflated counts for lowly expressed genes.

# GENCODE v44 transcripts + genome as decoy
grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v44.transcripts.fa GRCh38.primary_assembly.genome.fa > gentrome.fa
salmon index -t gentrome.fa -d decoys.txt -p 12 -i salmon_idx_v44 -k 31 --gencode

Then quantify each sample. -l A lets Salmon infer library type, --gcBias and --seqBias correct for fragment GC content and random-hexamer priming bias, and --numBootstraps 30 gives you inferential replicates you will want later if you ever look at transcript-level results.

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

Expect a mapping rate of 75–90% for polyA-selected human libraries. Anything under 60% means contamination, a wrong index, or heavy rRNA carryover, and you should resolve it before plotting anything. Collapse transcripts to genes in R with tximport, using countsFromAbundance = "lengthScaledTPM" so that changes in transcript-length distribution between samples do not masquerade as changes in gene expression.

library(tximport); library(readr)
tx2gene <- read_csv("tx2gene_gencode_v44.csv")   # tx_id, gene_id
files <- file.path("quant", samples, "quant.sf"); names(files) <- samples
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

2. Check QC before you plot anything

Three numbers per sample decide whether the downstream plots mean anything. First, library size: the column sums of the count matrix. Below about 10 million assigned reads, low-expression genes become noise-dominated and PCA starts separating samples by depth rather than biology. Second, the number of detected genes, meaning genes with at least one count. A whole-blood polyA library typically detects 13,000–16,000 genes. A sample detecting 8,000 is degraded or shallow. Third, the fraction of counts taken by the top handful of genes. In whole blood collected without globin depletion, hemoglobin transcripts (HBB, HBA1, HBA2) routinely consume more than half the library, which pushes everything else down and creates a composition artifact that PCA will find immediately.

cts <- txi$counts
libsize   <- colSums(cts)
detected  <- colSums(cts > 0)
top_frac  <- apply(cts, 2, function(x) sum(sort(x, decreasing = TRUE)[1:10]) / sum(x))
data.frame(libsize, detected, top_frac)

A reference-free sanity check is also useful when you suspect contamination rather than degradation. Counting k-mers directly in the reads, with a tool such as MerCat2, characterizes library composition without committing to a reference annotation, which is how you catch a bacterial or adapter-dominated library that a transcriptome-only pipeline quietly maps at low rate 1. If globin dominates, the right move is not statistical correction but removing those genes from the matrix before variance stabilization, and noting that your effective library depth was lower than the raw number suggests.

3. Build the matrix you will plot

Raw counts are wrong for every plot in this guide. Counts have a mean-variance relationship: high-expression genes have larger absolute variance, so a PCA on raw or naively log-transformed counts is dominated by the most-expressed genes, and a heatmap of raw counts shows library size, not biology. Use DESeq2’s variance stabilizing transformation, which produces values on roughly the log2 scale with variance approximately constant across the expression range.

library(DESeq2)
coldata <- read.csv("coldata.csv", row.names = 1)  # timepoint, batch, condition
dds <- DESeqDataSetFromTximport(txi, colData = coldata, design = ~ batch + condition)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]      # keep genes with 10+ counts in 3+ samples
vsd <- vst(dds, blind = FALSE)                      # blind=FALSE: use design to estimate dispersion
write.csv(assay(vsd), "vst_matrix.csv")

Use vst() rather than rlog() unless you have fewer than about 30 samples and highly variable sequencing depth, in which case rlog() is slightly better behaved and the runtime penalty does not matter. Set blind = FALSE when you already know your design, and blind = TRUE only when you want a transformation that is genuinely independent of the groups, for example when the plot is meant to be evidence about whether the groups separate at all.

To answer a question people ask constantly: you can compute a PCA in Excel, and you should not. Excel has no variance stabilization, no way to restrict to high-variance genes reproducibly, and a well-documented history of silently converting gene symbols like SEPT9 and MARCH1 into dates. Do this in Python or R.

4. Make the PCA plot

PCA is the first plot to make and the one that most often changes what you do next. It reduces a 15,000-dimensional expression vector per sample into two or three coordinates that capture the largest sources of variance, and it tells you whether the dominant structure in your data is the thing you care about or something else. Restrict the input to the most variable genes: DESeq2’s built-in plotPCA uses the top 500 by variance, and 500 to 2,000 is a reasonable range. Including all genes dilutes the signal with thousands of flat, low-information rows.

import pandas as pd, numpy as np
import matplotlib.pyplot as plt
from sklearn.decomposition import PCA

vst = pd.read_csv("vst_matrix.csv", index_col=0)      # genes x samples
meta = pd.read_csv("coldata.csv", index_col=0)

ntop = 1000
v = vst.var(axis=1).sort_values(ascending=False)
X = vst.loc[v.index[:ntop]].T.values                   # samples x genes
X = X - X.mean(axis=0)                                 # center; do NOT scale to unit variance

p = PCA(n_components=5)
S = p.fit_transform(X)
ev = p.explained_variance_ratio_ * 100

fig, ax = plt.subplots(figsize=(5,5))
for grp, idx in meta.groupby("condition").groups.items():
    pos = [vst.columns.get_loc(s) for s in idx]
    ax.scatter(S[pos,0], S[pos,1], label=grp, s=60)
for i, s in enumerate(vst.columns):
    ax.annotate(s, (S[i,0], S[i,1]), fontsize=7)
ax.set_xlabel(f"PC1 ({ev[0]:.1f}%)"); ax.set_ylabel(f"PC2 ({ev[1]:.1f}%)")
ax.legend(); ax.set_aspect("equal"); fig.tight_layout()
fig.savefig("pca.png", dpi=200)

Two details matter. Center the columns but do not scale each gene to unit variance: after VST the variances are already comparable, and rescaling amplifies the noisiest low-expression genes. Set the aspect ratio to equal, because unequal axes make a PC1 that explains 60% of variance look the same width as a PC2 that explains 5%, which misleads the eye about how separated the groups really are.

Then read the plot against the scree values. Print ev and check whether PC1 explains 50% or 12%. Color the same scatter by library size, by batch, and by RNA integrity number if you have it. If PC1 correlates with any of those at r > 0.8, your leading component is technical. For longitudinal data from one person, the informative version of this plot colors points by time and connects them in order: a sample trajectory that traces a smooth path through PC space is the usual signature of a real ordered biological process, which is how ordination has been used to read out staged transcriptional programs in multi-omics time courses 2. Correlate PC scores with candidate covariates rather than eyeballing them:

for pc in range(3):
    print(pc+1, np.corrcoef(S[:,pc], meta["libsize"].values)[0,1])

If batch dominates and you cannot redesign the experiment, you can produce a corrected matrix for visualization with limma::removeBatchEffect(assay(vsd), batch = vsd$batch), or use ComBat_seq on raw counts before the DESeq2 model. Use the corrected matrix for plots only. The differential expression model should include batch as a term instead, which handles the same confounding without distorting the data.

For more than a handful of sample groups, PCA is not the last word. Joint ordination methods designed for multiple data types keep the compositional structure of each modality straight while finding shared axes, which matters once you are plotting RNA alongside proteomics or microbiome counts in the same space 3.

5. Confirm the PCA with a sample-distance heatmap

PCA compresses distances into two dimensions and can hide structure in PC3 and beyond. A heatmap of pairwise Euclidean distances on the VST matrix is the cheapest cross-check, and it uses all the retained variance rather than two components.

library(pheatmap)
sampleDists <- dist(t(assay(vsd)))
pheatmap(as.matrix(sampleDists),
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         annotation_col = as.data.frame(colData(vsd)[, c("condition","batch")]))

If the heatmap’s dendrogram groups samples the same way the PCA does, you can trust the picture. If the two disagree, the structure lives in higher components, and you should plot PC2 against PC3 before drawing any conclusion. Spearman correlation between samples is a useful companion view: technical replicates from the same library should sit above 0.97, and biological replicates from the same person at different times usually land between 0.90 and 0.97 in blood.

6. Fit the model, then make MA and volcano plots

Differential expression is a model fit, and the MA plot is how you check the fit before you interpret any gene. It plots log2 fold change against mean normalized expression, so a healthy result looks like a symmetric cloud centered on zero that narrows as expression increases. Fold changes for low-expression genes are noisy by nature, which is why you shrink them.

dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "post", "pre"), alpha = 0.05)
resLFC <- lfcShrink(dds, coef = "condition_post_vs_pre", type = "apeglm")
plotMA(resLFC, ylim = c(-4, 4))
summary(res)

A trumpet shape opening to the left is expected. A cloud that is visibly off-center means the normalization assumption (most genes unchanged) is violated, which happens with globin-dominated blood or with a global transcriptional shift. A fan that never narrows at high expression means your replicates are more variable than the model expects, and you should revisit the design.

The volcano plot then shows effect size against evidence: log2 fold change on x, −log10 adjusted p-value on y. Use the shrunken fold changes on the x-axis so that the far-left and far-right points are genes with real effects rather than genes with three counts.

import numpy as np, pandas as pd, matplotlib.pyplot as plt
r = pd.read_csv("resLFC.csv", index_col=0).dropna(subset=["padj"])
sig = (r.padj < 0.05) & (r.log2FoldChange.abs() > 1)
plt.figure(figsize=(5,5))
plt.scatter(r.log2FoldChange, -np.log10(r.padj), s=6, c="lightgray")
plt.scatter(r.log2FoldChange[sig], -np.log10(r.padj[sig]), s=8, c="crimson")
for g in r[sig].sort_values("padj").head(15).index:
    plt.annotate(g, (r.log2FoldChange[g], -np.log10(r.padj[g])), fontsize=7)
plt.axhline(-np.log10(0.05), ls="--", lw=0.7); plt.axvline(1, ls="--", lw=0.7); plt.axvline(-1, ls="--", lw=0.7)
plt.xlabel("log2 fold change (apeglm shrunk)"); plt.ylabel("-log10 adjusted p")
plt.tight_layout(); plt.savefig("volcano.png", dpi=200)

The |log2FC| > 1 threshold is a convention, not a statistical statement. If you want a threshold with inferential meaning, test against it directly with results(dds, lfcThreshold = 1), which changes the null hypothesis rather than filtering after the fact.

7. Plot longitudinal trajectories for a single person

If your data are repeated samples from one individual rather than a case-control comparison, the volcano plot is the wrong primary output, because there are no groups to contrast. What you want is each gene’s trajectory over time against that person’s own baseline, plus an estimate of within-person variability so you know what a meaningful deviation looks like.

Take at least three or four baseline timepoints under comparable conditions (same time of day, same fasting state, same collection tube), compute the per-gene mean and standard deviation on the VST scale, and express later samples as z-scores against that personal reference. Genes with within-person SD above roughly 0.5 on the VST scale are too variable in that person to interpret from a single draw.

base = vst[baseline_samples]
mu, sd = base.mean(axis=1), base.std(axis=1)
z = vst.sub(mu, axis=0).div(sd.clip(lower=0.1), axis=0)

genes = ["IFI27", "ISG15", "MKI67"]
fig, ax = plt.subplots(figsize=(6,4))
for g in genes:
    ax.plot(meta["day"], z.loc[g, meta.index], marker="o", label=g)
ax.axhline(0, color="k", lw=0.6); ax.set_ylabel("z vs personal baseline"); ax.legend()

For modeling rather than plotting, fit time as a spline in the DESeq2 design (~ ns(day, df = 3)) and use test = "LRT" against a reduced model to find genes with any time-dependent pattern. Interpretation of these plots for health purposes belongs with a clinician. A transcriptional shift you can see in a graph is a measurement, and mapping it onto anything about your health requires a physician who can put it next to your history and a clinically validated test.

8. Go beyond single-modality plots

Once you have RNA plus proteins, metabolites, or genotype, single-omic PCA stops being sufficient, because the interesting structure lies in the relationships between modalities. Two families of methods are worth knowing. Graph-based integration represents genes and samples as nodes with edges encoding regulatory or correlational relationships, and learns embeddings across modalities. The design choices in how you build that graph, which nodes and which edges, drive performance more than the choice of network architecture 4. Attention-based graph models extended to regulatory network inference use directed prior knowledge to orient edges, producing graphs you can read as hypotheses about regulation rather than as undirected correlation 5.

The second family is interpretable visualization on top of learned representations. Methods that read out neural network training dynamics, rather than only the final embedding, recover structure in single-cell and spatial data that a static UMAP obscures 6. Integrated AI/ML pipelines for multi-omics biomarker discovery increasingly ship three-dimensional interactive visualizations of the combined feature space as their primary output rather than a single static scatter 7. Before reaching for any of them, get the single-modality plots right, because a PCA driven by library size will not be fixed by a graph neural network.

Common problems

Samples separate by batch on PC1 and by nothing else. Confirm with the correlation check in step 4, add batch to the design formula for testing, and use removeBatchEffect only for the figure. If batch is perfectly confounded with condition, no method recovers the signal, and the experiment needs to be redesigned.

The volcano plot has thousands of significant genes with tiny fold changes. This usually means very high power (many samples) plus a real but small global effect, or an unmodeled covariate. Add lfcThreshold to the results call and look at the MA plot for asymmetry.

The PCA looks like one tight cluster plus one distant outlier. Check that sample’s library size and detected-gene count first. If both are normal, check the top-expressed genes for globin or rRNA dominance. Removing an outlier is defensible only with a technical reason you can state.

Gene symbols in your matrix have turned into dates. The file passed through Excel. Re-export from R or Python, and use Ensembl gene IDs as the index with symbols as a separate column.

Fold changes disagree between Salmon and a count-based pipeline like featureCounts. This is expected for genes with multiple isoforms of different length, and it is the reason for countsFromAbundance = "lengthScaledTPM". Check the transcript-level results for that gene before assuming either pipeline is wrong.

PCA looks fine but replicates correlate at 0.85. Suspect RNA degradation. Look at the 3′ bias in your provider’s QC report, and if DV200 or RIN was low, treat low-expression genes with particular skepticism.

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. Jose L Figueroa, Andrew Redinbo, Ajay Panyala, et al. MerCat2: a versatile k -mer counter and diversity estimator for database-independent property analysis obtained from omics data. Bioinformatics Advances, 2024. https://doi.org/10.1093/bioadv/vbae061 ↩

  2. Steffen Israel, Mathias Ernst, Olympia E. Psathaki, et al. An integrated genome-wide multi-omics analysis of gene expression dynamics in the preimplantation mouse embryo. Scientific Reports, 2019. https://doi.org/10.1038/s41598-019-49817-3 ↩

  3. Bianca Cordazzo Vargas, Cameron Martino, Amanda Hazel Dilmore, et al. Joint-RPCA: domain-aware multi-omics integration for systems microbiology. Molecular Systems Biology, 2026. https://doi.org/10.1038/s44320-026-00236-3 ↩

  4. Muhtasim Noor Alif, Khandakar Tanvir Ahmed, Sudipto Baul, et al. Graph designs for deep learning–based multi-omics integration. Briefings in Bioinformatics, 2026. https://doi.org/10.1093/bib/bbag410 ↩

  5. Noor Jamal Alkhateeb, Mamoun Awad. MultiCausGRN: directed prior-guided graph attention model for multi-omics gene regulatory network inference. Frontiers in Bioinformatics, 2026. https://doi.org/10.3389/fbinf.2026.1883130 ↩

  6. Jonathan Karin, Reshef Mintz, Barak Raveh, et al. Interpreting single-cell and spatial omics data using deep neural network training dynamics. Nature Computational Science, 2024. https://doi.org/10.1038/s43588-024-00721-5 ↩

  7. Rishabh Narayanan, Elizabeth Peker, William DeGroat, et al. 3D IntelliGenes: AI/ML application using multi-omics data for biomarker discovery and disease prediction with multi-dimensional visualization. BMC Medical Research Methodology, 2025. https://doi.org/10.1186/s12874-025-02649-4 ↩