How to Read Your Own RNA-seq Expression Levels

By the end of this guide you will have a gene-by-sample matrix of expression estimates from your own RNA sequencing data, in three forms that serve different purposes. Raw counts are for statistical modeling. Transcript-per-million (TPM) values are for asking “how much of this gene is in this sample”, and variance-stabilized values are for clustering and visualization. You will also have a defensible sense of which numbers in that matrix mean something and which are noise. You need: paired-end FASTQ files from a bulk RNA-seq run (a typical whole-blood library at 30–50 million read pairs), roughly 40 GB of disk, 16 GB of RAM, and a machine where you can install command-line tools. If you have only one sample, say one blood draw, you can still compute expression levels. You cannot compute differential expression, and the distinction matters more than most first-time readers expect. Everything below is measurement and interpretation. Nothing here is a diagnosis, and any result you want to act on medically belongs in front of a clinician who can see the rest of your chart.

1. Understand what an “expression level” is before you compute one

An expression level from RNA-seq is a proportion. It is not a concentration. The sequencer samples fragments from a library built out of the RNA in your tube, and the number of reads assigned to a gene reflects that gene’s share of the library, scaled by however many reads you bought. Two facts follow. First, all RNA-seq measurements are compositional: if one transcript family balloons, every other gene’s share falls even if its absolute copy number per cell is unchanged. In whole blood, the hemoglobin transcripts HBB, HBA1 and HBA2 routinely consume 30–70% of reads in an unglobin-depleted library. That compresses everything else. Second, a gene’s read count scales with its length, because longer transcripts yield more fragments at equal molar abundance. Any comparison across genes has to divide that out, which is what TPM does.

It also helps to hold the biology in mind. Gene expression proceeds in an order: chromatin becomes accessible, transcription produces pre-mRNA, splicing and polyadenylation produce mature mRNA, the transcript is exported and translated, and both mRNA and protein are degraded on their own timescales. RNA-seq measures one point in that chain, the steady-state mRNA pool, and it measures it well. It does not measure transcription rate, and it predicts protein abundance only loosely, because translation efficiency and protein turnover vary by orders of magnitude across genes. Cross-omics comparisons in model organisms show that mRNA and protein levels share broad expression strategies but diverge substantially gene by gene.¹ This is the single most common over-read of an expression matrix, so it is worth internalizing early: high mRNA means a lot of that transcript was present in that sample, and nothing more.

2. Check read quality and library composition first

Run FastQC and aggregate with MultiQC. This takes minutes and catches the problems that would otherwise silently corrupt your matrix.

fastqc -t 8 -o qc/ sample_R1.fastq.gz sample_R2.fastq.gz
multiqc qc/ -o qc/report

The things to look at, in order of how often they matter: per-base sequence content at the 5’ end (the first 10–12 bases in a random-primed library will look biased, and that is normal), duplication level, adapter content, and the overrepresented-sequences table. That last table is where you will see globin. If your top overrepresented sequence maps to HBB, you have a whole-blood library without globin depletion, and your effective depth for everything else is a fraction of your nominal depth.

Also note which library chemistry you have, because it changes interpretation. PolyA selection captures polyadenylated mRNA and gives you a clean, mRNA-dominated library. Ribosomal RNA depletion keeps non-polyadenylated species, including many long non-coding RNAs and histone mRNAs, and it retains more intronic signal. Paired data from human T cells show the two protocols agree well for protein-coding genes but diverge for non-coding and unprocessed transcripts, so a gene absent from a polyA library may be well-measured in an rRNA-depleted one.² Know which you bought before you conclude a gene is off.

3. Quantify against a transcriptome with selective alignment

We recommend Salmon with a decoy-aware index and GC bias correction. It quantifies at the transcript level in about ten minutes per sample on a laptop, handles multi-mapping reads by probabilistic assignment rather than discarding them, and produces the effective-length estimates you need for correct TPM. Full genome alignment with STAR is the right choice if you also care about splice junctions, novel transcripts, or variant calling from RNA, but it costs roughly 32 GB of RAM and considerably more time, and for expression levels alone the gain is small.

Build the index once, from GENCODE transcript sequences plus the genome as decoy sequence. Decoys matter: without them, reads from unannotated genomic regions get force-assigned to whatever transcript they resemble most, inflating low-expression genes.

# concatenate transcriptome + genome, mark genome sequences as decoys
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 12 \
  -i salmon_index_gencode45 --gencode -k 31

Then quantify:

salmon quant -i salmon_index_gencode45 -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  -p 12 --validateMappings --seqBias --gcBias \
  --numGibbsSamples 30 \
  -o quant/sample01

A few flags earn their place. Setting -l A lets Salmon infer library type, and you should read the inferred type out of quant/sample01/lib_format_counts.json to confirm it matches what the protocol claims. A stranded library reported as unstranded usually means something went wrong at library prep. Adding --gcBias corrects fragment-level GC dependence, which is the largest systematic bias left in most modern libraries. Setting --numGibbsSamples 30 gives you posterior samples of transcript abundance, which is how you find out whether a transcript-level estimate is identifiable at all. Check the mapping rate in logs/salmon_quant.log: for human blood with a k=31 index, 75–90% is typical, and anything under 60% means contamination, wrong species, or adapter carryover.

4. Collapse to gene level and build the three matrices

Quantify transcripts, then sum to genes. Transcript-level estimates are genuinely uncertain for genes with many similar isoforms, and gene-level inference is better powered and more reproducible for most questions a person is asking about their own biology. Use tximport with countsFromAbundance = "lengthScaledTPM", which produces counts that are consistent with the abundance estimates and already account for changes in average transcript length across samples.

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

files <- file.path("quant", samples$id, "quant.sf")
names(files) <- samples$id

# tx2gene from the GENCODE FASTA headers or a GTF
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM",
                ignoreTxVersion = TRUE)

dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ timepoint)
keep <- rowSums(counts(dds) >= 10) >= 3      # gene seen in >=3 samples
dds <- dds[keep, ]
dds <- DESeq(dds)

tpm  <- txi$abundance[keep, ]                 # cross-gene comparable
vsd  <- assay(vst(dds, blind = FALSE))        # for PCA / clustering / heatmaps
norm <- counts(dds, normalized = TRUE)        # median-of-ratios scaled counts

Keep all three and use each for its job. TPM is the one to read when you want to know how abundant a gene is in a sample, because it divides out both transcript length and library size. It is the wrong input for differential testing, because it discards the count information that the negative binomial model needs to estimate dispersion. The variance-stabilized matrix is what you feed to PCA, correlation heatmaps, and any distance calculation, because untransformed counts let the few thousand highest-expressed genes dominate every distance. Never run statistics on z-scores of expression: a z-score tells you where a sample sits relative to your other samples for that gene, which is useful for a heatmap and meaningless as evidence on its own.

5. Read the numbers

Now the interpretation. In a typical whole-blood polyA library with 30 million assigned reads, you will detect roughly 12,000–15,000 genes at non-trivial levels, and the TPM distribution will span six orders of magnitude. Useful anchors: TPM below about 1 means the gene is at or near the detection floor and you should treat presence/absence as unresolved. TPM of 1–10 is low but real expression. TPM of 10–100 is the middle of the distribution where most expressed protein-coding genes sit. TPM above 1,000 means the gene is a major component of the library, which in blood is globins, ribosomal protein genes, mitochondrial genes, and a handful of immune transcripts like S100A8/A9 and LYZ. These bands are guides to how much statistical weight a number can bear. They are not clinical thresholds.

The constraint on a single sample is that you have no reference distribution. “MYC at 45 TPM” means nothing by itself, because you do not know your own baseline and cross-study reference ranges are confounded by library protocol, globin depletion, time of day, and cell composition. Longitudinal sampling changes the character of the analysis for that reason. With four or more draws from yourself over time, you become your own control. A two-fold change in a gene against your own history is then interpretable in a way that a single absolute value is not. RNA-seq precision studies show that technical replicate reproducibility is high and that most observed variance in well-run experiments is biological rather than technical, which is what makes within-person comparison over time worth doing.³ Depth also matters more for low-expressed genes than high: doubling reads sharpens estimates in the TPM 1–10 range and barely moves the top decile.⁴

One more interpretive trap specific to blood. A large share of the variance in whole-blood expression across draws is cell composition. Per-cell regulation does not account for it. If your neutrophil fraction rose from 55% to 70% between two draws, hundreds of neutrophil-specific genes will look “up” and lymphocyte genes will look “down”, and no gene changed its regulation. Estimate cell fractions from the expression matrix itself using a deconvolution method, or better, get a complete blood count from the same draw and include the fractions as covariates in your design. Interpreting bulk blood expression without cell composition is the most common way to arrive at a confident wrong conclusion.

6. Test for change over time, if you have the samples

With three or more timepoints you can model change properly. Use DESeq2 with a design that encodes the contrast you care about, and shrink the fold changes before you rank anything.

res <- results(dds, contrast = c("timepoint", "T3", "T1"), alpha = 0.05)
res <- lfcShrink(dds, coef = "timepoint_T3_vs_T1", type = "apeglm", res = res)

sum(res$padj < 0.05 & abs(res$log2FoldChange) > 0.585, na.rm = TRUE)  # >1.5x

Two points on reading the output. The padj column is a Benjamini-Hochberg adjusted p-value across roughly 15,000 tests, so a nominal p of 0.01 is not a finding. And lfcShrink exists because raw log fold changes for low-count genes are wildly unstable: a gene going from 2 counts to 8 counts is a 4x change and essentially no evidence. Shrunken estimates pull those toward zero in proportion to their uncertainty, which makes the ranked list usable.

Once you have a ranked list, interpret at the pathway level rather than gene by gene. Single-gene stories from a personal expression matrix are fragile. Coordinated shifts across a pathway are more likely to reflect something real, and this is also where expression data becomes most informative when combined with other layers, since transcript changes constrain but do not determine protein and metabolite changes.⁵ Reviews of multi-omics integration make the same point from the other direction. The value comes from cross-layer consistency. It does not come from any single layer read in isolation.⁶⁷ And because every one of these analyses involves many tests over many genes, the interpretation step needs the same skepticism about multiplicity and model assumptions that the statistics step did.⁸

Common problems

Low mapping rate (under 60%). Check for adapter read-through with FastQC, confirm you used a human index, and check lib_format_counts.json for a strandedness mismatch. A stranded library quantified as unstranded typically loses a large fraction of assignments.

Globin domination. If HBB alone takes 40% of reads, your effective depth for everything else is 60% of nominal, and your low-expression detection suffers accordingly. The fix is at the wet-lab stage (globin depletion during library prep). It is not in software. You can drop globin genes from the matrix before normalizing, which helps the size-factor estimate, but it does not recover the lost reads.

High duplication rate. In RNA-seq, duplicates are expected for highly expressed genes and do not indicate a problem the way they would in whole-genome sequencing. Do not deduplicate bulk RNA-seq without unique molecular identifiers. You will preferentially delete real signal from the most-expressed genes.

3’ bias and RNA degradation. If coverage falls off sharply toward transcript 5’ ends, your RNA integrity number was low. Degradation shortens effective transcript lengths unevenly and systematically distorts TPM for long transcripts. Check the RIN from the lab report before trusting long-gene estimates.

Transcript-level estimates that will not hold still. Look at the Gibbs posterior samples from step 3. If a transcript’s posterior spans an order of magnitude, the reads cannot distinguish that isoform from its siblings. You should not report the isoform in that case. Report the gene. This is common for genes with many overlapping annotated isoforms.

Batch effects across timepoints. If draws were sequenced months apart on different runs, batch will show up as the first principal component. Include batch in the DESeq2 design, and where possible have the lab sequence stored libraries together. Do not correct batch out and then run statistics on the corrected values without accounting for it in the model.

Conclusions that outrun the assay. RNA-seq tells you what transcript pool was present in that tube at that moment. It does not tell you protein abundance, enzyme activity, or clinical risk, and the read-through from mRNA to phenotype is where most over-interpretation happens.⁹ If a result looks medically consequential, take the data and the analysis to a physician rather than acting on it.

Footnotes

1. Jaume Forés-Martos, Anabel Forte, José García-Martínez, et al. A Trans-Omics Comparison Reveals Common Gene Expression Strategies in Four Model Organisms and Exposes Similarities and Differences between Them. Cells, 2021. https://doi.org/10.3390/cells10020334 ↩

2. Li Chen, Ruirui Yang, Tony Kwan, et al. Paired rRNA-depleted and polyA-selected RNA sequencing data and supporting multi-omics data from human T cells. Scientific Data, 2020. https://doi.org/10.1038/s41597-020-00719-4 ↩

3. Paweł P. Łabaj, Germán G. Leparc, Bryan E. Linggi, et al. Characterization and improvement of RNA-Seq precision in quantitative transcript expression profiling. Bioinformatics, 2011. https://doi.org/10.1093/bioinformatics/btr247 ↩

4. Paul A McGettigan. Transcriptomics in the RNA-seq era. Current Opinion in Chemical Biology, 2013. https://doi.org/10.1016/j.cbpa.2012.12.008 ↩

5. Gun Tae Jung, Kwang-Pyo Kim, Kwoneel Kim. How to interpret and integrate multi-omics data at systems level. Animal Cells and Systems, 2020. https://doi.org/10.1080/19768354.2020.1721321 ↩

6. Joerg Martin Buescher, Edward M Driggers. Integration of omics: more than the sum of its parts. Cancer & Metabolism, 2016. https://doi.org/10.1186/s40170-016-0143-y ↩

7. Biswapriya B Misra, Carl Langefeld, Michael Olivier, et al. Integrated omics: tools, advances and future approaches. Journal of Molecular Endocrinology, 2019. https://doi.org/10.1530/jme-18-0055 ↩

8. Ryo Yamada, Daigo Okada, Juan Wang, et al. Interpretation of omics data analyses. Journal of Human Genetics, 2020. https://doi.org/10.1038/s10038-020-0763-5 ↩

9. Xi Qian, Yi Ba, Qianfeng Zhuang, et al. RNA-Seq Technology and Its Application in Fish Transcriptomics. OMICS: A Journal of Integrative Biology, 2014. https://doi.org/10.1089/omi.2013.0110 ↩