How to Normalize Your Own RNA-seq Data

By the end of this guide you will have three artifacts from your own RNA-seq data: a gene-level count matrix that has not been silently transformed by anything upstream, a set of per-sample scaling factors you computed yourself and can inspect, and two derived matrices for downstream work, log-CPM for ranking and testing, and a variance-stabilized matrix for clustering and any distance-based method. You will also have diagnostic plots that tell you whether the normalization assumption was violated. You need the following.

• a gene- or transcript-level quantification per sample (Salmon quant.sf, kallisto abundance.h5, or a featureCounts matrix)
• the matching annotation used for quantification
• R 4.3 or later with four packages: edgeR, DESeq2, tximport and limma
• roughly 8 GB of RAM for a few dozen human samples

Everything below is measurement and interpretation. Deciding what an expression change means for your health is a conversation with a clinician. It is not a normalization step.

1. Start from counts and do not use FPKM or TPM

The single most consequential choice happens before you normalize anything: what quantity you carry into the analysis. RPKM and FPKM divide by the library’s total mapped reads, which makes them within-sample quantities that are not comparable across samples. If one sample is dominated by a handful of very high-expression transcripts, the denominator absorbs that dominance and every other gene in that sample is pushed down. TPM at least has the property that its values sum to a constant across samples, which makes it a proportion, but that also means it inherits every problem of proportions.

Use estimated counts. With Salmon, quantify with bias correction on, because effective-length estimation and fragment GC bias interact with the normalization you do later:

salmon quant -i gencode_v44_index -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  -p 16 -o quant/sample01

Then aggregate to genes with tximport. The countsFromAbundance argument is where the within-sample length correction happens, and lengthScaledTPM is what we use: it produces gene-level counts scaled so that library size differences are preserved while transcript-length and effective-length differences across samples are corrected.

library(tximport); library(readr)
tx2gene <- read_tsv("gencode_v44_tx2gene.tsv", col_names = c("tx","gene"))
files <- file.path("quant", samples$id, "quant.sf"); names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

If you instead have a featureCounts matrix, you have raw counts with no effective-length correction, which is acceptable for gene-level differential testing but will bias any comparison you make between different genes in the same sample. Keep that limitation in mind and do not treat a featureCounts matrix as a set of expression levels you can rank across genes.

2. Filter low-count genes before you compute anything

Filtering comes first because normalization factors estimated from a matrix full of near-zero genes are noisier and, in the case of median-of-ratios, can be undefined. In a typical human polyA library against GENCODE, more than 40% of annotated genes carry almost no reads, and those rows contribute nothing but variance.

library(edgeR)
y <- DGEList(counts = txi$counts, samples = samples)
keep <- filterByExpr(y, group = samples$condition)   # min.count = 10 by default
y <- y[keep, , keep.lib.sizes = FALSE]
sum(keep)   # expect ~14,000-18,000 for human whole blood, more for tissue

filterByExpr works on CPM rather than raw counts, choosing a CPM cutoff equivalent to about 10 reads in the smallest library and requiring that threshold to be met in as many samples as the smallest group. Setting keep.lib.sizes = FALSE recomputes library sizes from the retained genes, which is what you want before estimating scaling factors.

3. Know which of the two problems each method solves

Normalization in RNA-seq is two separate corrections that are often described as one. The first is within-sample: reads accumulate in proportion to transcript length, so a 10 kb transcript and a 1 kb transcript at the same molar concentration produce very different counts. That is what effective-length division in TPM, and lengthScaledTPM in tximport, addresses.

The second is between-sample, and it has two components. Sequencing depth is the easy one: a library with 40 million reads yields roughly twice the counts of one with 20 million at the same underlying concentration. Composition is the hard one. Because sequencing measures relative abundance under a fixed read budget, a large increase in a few transcripts mathematically reduces the count share of everything else, even though nothing changed in their absolute concentration. This is why total-count scaling is not sufficient and why methods like TMM and median-of-ratios exist: they estimate a scaling factor from a central subset of genes rather than from the sum. The choice among between-sample methods is best understood as a choice among assumptions about what is invariant across your samples, and those assumptions are explicit and checkable ¹.

4. Compute scaling factors with TMM or median-of-ratios

Both methods rest on the same assumption. Most genes are not differentially expressed between any two samples, so an outlier-resistant central tendency of the per-gene ratio between two samples estimates the technical scaling difference. The methods differ in how they find that center.

TMM (trimmed mean of M-values) picks a reference sample whose upper-quartile count is closest to the mean upper-quartile across samples. For every other sample it computes, per gene, the average log expression A and the log ratio M relative to the reference. It then trims the most extreme 30% of M values and 5% of A values, and takes an inverse-variance-weighted mean of what remains. That weighted mean, exponentiated, is the normalization factor.

y <- calcNormFactors(y, method = "TMM")   # logratioTrim = 0.3, sumTrim = 0.05
y$samples$norm.factors

Read those factors. They are multiplicative adjustments to the effective library size, centered so their product is one. Values inside roughly 0.85 to 1.2 are unremarkable. A factor of 0.6 or 1.7 says that sample’s composition differs sharply from the reference, and you should find out why before proceeding.

DESeq2’s median-of-ratios takes a per-gene geometric mean across all samples as a synthetic reference, divides each sample’s count by it, and takes the median of those ratios as the size factor. Genes with a zero in any sample are dropped from the calculation, which is why step 2 matters.

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ condition)
dds <- dds[rowSums(counts(dds)) > 0, ]
dds <- estimateSizeFactors(dds)
sizeFactors(dds)
norm_counts <- counts(dds, normalized = TRUE)

We use TMM when the downstream work is in limma/voom or when we want log-CPM values to carry around, and median-of-ratios when the downstream work is a DESeq2 negative-binomial test. The practical difference between them is small on well-behaved data, and comparative evaluations have repeatedly placed both in the top tier while showing that total-count and quantile normalization behave worse when library composition varies ². If your samples include a condition where a large fraction of the transcriptome genuinely shifts in one direction, use estimateSizeFactors(dds, type = "poscounts") or switch to a control-gene-based approach in step 7, because the “most genes unchanged” assumption is the thing that failed.

5. Produce the matrices you will use

Different downstream tasks want different scales, and the common mistake is to use one matrix for everything. Three outputs cover nearly all cases.

For anything involving fold changes, ranking, or linear modeling, use log-CPM with a prior count that damps the variance of low-count genes:

logcpm <- cpm(y, log = TRUE, prior.count = 2, normalized.lib.sizes = TRUE)

For differential testing with limma, do not feed it log-CPM directly. voom estimates the mean-variance relationship and turns it into precision weights, which is the point of the method:

design <- model.matrix(~ condition, data = y$samples)
v <- voom(y, design, plot = TRUE)
fit <- eBayes(lmFit(v, design))

For clustering, PCA, correlation heatmaps, and any distance metric, use a variance-stabilizing transform so that high-expression genes do not dominate the distance purely by having larger variance on the log scale:

vsd <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = "batch")

Set blind = FALSE when you have a design you trust. The blind = TRUE setting is for the case where you want the transform to know nothing about your groups, which is more conservative and slightly less efficient. The rlog transform is an alternative that behaves better on very small sample counts but scales poorly past a few dozen samples.

On the two scaling questions people often ask directly: normalizing a sample to 100% is exactly what CPM and TPM do, since both convert counts to proportions with a different constant, and once you have done that your data are compositional and standard correlation and distance measures on them are distorted. Min-max scaling a gene to a 0-to-1 range, (x - min(x)) / (max(x) - min(x)), is a presentation choice for heatmaps and should never precede a statistical test, because it makes every gene’s range identical and discards the magnitude information the test needs. If you want to work in a compositionally coherent space, use a log-ratio transform such as centered log-ratio instead of raw proportions ³.

6. Check that the normalization held

Normalization is an estimate, and estimates need diagnostics. Two plots catch most failures. The relative log expression (RLE) plot shows, per sample, the distribution of each gene’s log count minus that gene’s median across samples. After correct normalization every box sits on zero with similar spread. A sample whose box is shifted or visibly wider is telling you that a single scaling factor did not capture what is different about it.

library(EDASeq)
plotRLE(cpm(y, normalized.lib.sizes = TRUE), outline = FALSE, las = 2)

The second is a per-sample MA plot against the reference or against the row-wise mean. A trend in M as a function of A means the discrepancy is expression-level dependent, which no single multiplicative factor can fix and which usually points to a length, GC, or degradation difference. If you want to compare several normalization choices quantitatively rather than by eye, published protocols score methods on how well they preserve known sample relationships and how they affect downstream results, and running that kind of comparison on your own data is inexpensive ⁴. Tools built for cross-omics comparison do the same thing over a grid of methods and report the spread of within-group variance under each ⁵.

Also expect the choice to propagate. In a benchmark of RNA-seq normalization methods used to build genome-scale metabolic models, different normalization and thresholding combinations produced noticeably different reconstructed networks from the same underlying counts, which is a concrete demonstration that this step is not cosmetic ⁶.

7. Handle the longitudinal single-person case

If your samples are you, sampled repeatedly, the structure of the problem changes. There is no biological replication at any timepoint, so what you are estimating is a personal trajectory against a personal baseline, and every non-biological source of variation between draws now competes directly with the signal you care about. Three things matter more than the choice between TMM and median-of-ratios.

First, identify and model batch explicitly. Record the following.

• collection time of day
• tube type
• time from draw to stabilization
• RNA integrity number
• library prep date
• flow cell

Put the ones you cannot randomize into the design matrix. Use limma::removeBatchEffect only to produce a matrix for plotting, never as input to a test, since the test needs to account for the degrees of freedom you spent.

Second, consider control-gene normalization when the composition assumption is in doubt. In whole blood, globin transcripts can consume a very large share of reads, and residual differences in globin or rRNA depletion efficiency between draws will move your scaling factors for reasons that have nothing to do with biology. RUVSeq estimates unwanted factors from either spike-in controls or from a set of samples you designate as replicates:

library(RUVSeq)
set <- betweenLaneNormalization(newSeqExpressionSet(counts), which = "upper")
set <- RUVg(set, cIdx = ercc_gene_ids, k = 2)   # ERCC spike-ins as controls

Choose k by looking at how many factors are needed before the RLE plot flattens, and be aware that increasing k will eventually start removing signal you wanted.

Third, keep the raw counts and the exact normalization code with the data. Six months later you will add three new samples, and if you re-estimate TMM factors on the expanded set every previously computed value shifts slightly. That is correct behavior and not a bug. It does mean the normalized matrix is a function of the cohort and must be regenerated. New samples cannot be appended to the old matrix.

8. Single-cell data follows different rules

If your data are single-cell rather than bulk, the between-sample logic above does not transfer directly, because each cell has a tiny library, a majority of zeros, and no length bias in 10x-style 3’ protocols. The standard starting point is per-cell scaling to a fixed total followed by log1p, which in Seurat is NormalizeData(obj, normalization.method = "LogNormalize", scale.factor = 1e4). That assumes total counts per cell reflect only technical capture efficiency, which is false when cell types genuinely differ in total RNA content. The scran::computeSumFactors function addresses this by pooling similar cells before deconvolving per-cell factors, and it is what we would use when cell-type composition is heterogeneous. The sctransform method fits a regularized negative-binomial model per gene instead and returns Pearson residuals, which behave well for clustering and poorly as inputs to methods expecting counts. Do not carry residuals into a differential expression test that assumes a count distribution.

Common problems

A few failures recur often enough to name. Normalization factors far from one across most samples usually means a library composition problem rather than a depth problem: check for globin, rRNA, and mitochondrial fraction per sample before touching the method. Counts that will not converge in estimateSizeFactors, or a size factor of NA, almost always means too many genes with zeros. Filter first or use type = "poscounts".

Beautiful clustering that tracks library prep date rather than condition means batch dominates, and no scaling factor will fix a shift that differs by gene. Model it. A strong trend on the MA plot after normalization points to a within-sample bias that must be corrected upstream, so re-quantify with --gcBias and --seqBias enabled rather than attempting to patch it downstream.

Finally, resist comparing normalized values across projects. Log-CPM from your samples and log-CPM from a public cohort are not on a shared scale unless you normalize the combined count matrix together, and even then annotation version, library chemistry, and read length will leave residual differences. When you need a cross-study comparison, requantify the public raw reads against your annotation and normalize jointly. Any clinical reading of the result belongs with a physician who can see the rest of your picture.

Footnotes

1. Ciaran Evans, Johanna Hardin, Daniel M Stoebel. Selecting between-sample RNA-Seq normalization methods from the perspective of their assumptions. Briefings in Bioinformatics, 2017. https://doi.org/10.1093/bib/bbx008 ↩

2. M.-A. Dillies, A. Rau, J. Aubert, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Briefings in Bioinformatics, 2012. https://doi.org/10.1093/bib/bbs046 ↩

3. Thomas P Quinn, Ionas Erb, Greg Gloor, et al. A field guide for the compositional analysis of any-omics data. GigaScience, 2019. https://doi.org/10.1093/gigascience/giz107 ↩

4. Zachary B. Abrams, Travis S. Johnson, Kun Huang, et al. A protocol to evaluate RNA sequencing normalization methods. BMC Bioinformatics, 2019. https://doi.org/10.1186/s12859-019-3247-x ↩

5. Aakash Chawade, Erik Alexandersson, Fredrik Levander. Normalyzer: A Tool for Rapid Evaluation of Normalization Methods for Omics Data Sets. Journal of Proteome Research, 2014. https://doi.org/10.1021/pr401264n ↩

6. Hatice Büşra Lüleci, Dilara Uzuner, Müberra Fatma Cesur, et al. A benchmark of RNA-seq data normalization methods for transcriptome mapping on human genome-scale metabolic networks. npj Systems Biology and Applications, 2024. https://doi.org/10.1038/s41540-024-00448-z ↩