Finding Differentially Expressed Genes in R From Your Own RNA-seq Data

Oak
rna-seq deseq2 edger r differential-expression bioinformatics
Two glowing reed stands of different colors face each other across a dark channel while glowing threads rise and sort by brightness.

By the end of this guide you will have a table of genes ranked by evidence for a difference in expression between two groups or two conditions in your own data, with shrunken effect sizes, adjusted p-values, and gene symbols attached, plus the diagnostic plots that tell you whether to believe any of it. You need: FASTQ files or an existing count matrix from an RNA-seq experiment, R 4.3 or later, Bioconductor 3.18 or later, and roughly 16 GB of RAM (quantification of a 40-million-read human sample with salmon needs about 8 GB for the index alone). You also need a sample sheet: one row per library, with columns for every variable you might want to model, including the ones you would rather ignore, such as sequencing run, RNA extraction date, and time of blood draw.

The core idea is simple. A gene is called differentially expressed when the number of reads assigned to it, after correcting for how deeply each library was sequenced, differs between conditions by more than the noise in your data can plausibly explain. Everything hard about the analysis lives in the phrase “the noise in your data,” which is why most of this guide is about the model rather than the test.

1. Get counts you trust, and pin the annotation

Start from FASTQ files if you have them. We prefer selective alignment with salmon over a full genome alignment for differential expression, because it is roughly ten times faster, it handles multi-mapping reads across transcript isoforms probabilistically, and it produces effective length estimates that let you convert between counts and TPM correctly.

salmon index -t gencode.v44.transcripts.fa.gz -i salmon_idx_v44 \
  --gencode -k 31 -p 16 -d decoys.txt

salmon quant -i salmon_idx_v44 -l A \
  -1 sample01_R1.fastq.gz -2 sample01_R2.fastq.gz \
  --validateMappings --seqBias --gcBias --posBias \
  --numBootstraps 30 -p 16 -o quants/sample01

Use a decoy-aware index built from the genome, otherwise reads from unannotated or intronic regions get forced onto transcripts and inflate counts for whatever is nearest. Turn on --gcBias unless you have a specific reason not to: GC bias varies between library preparation batches and is one of the few technical artifacts that mimics real biological signal. The --numBootstraps flag costs little and gives you per-transcript inferential replicates if you later want to do transcript-level work with swish or sleuth.

Record the exact annotation version in your analysis script and never mix versions across samples. This matters more than it sounds. A recent reanalysis of Alzheimer’s brain RNA-seq found that changing the reference genome and annotation build altered not only which genes were called differentially expressed but the direction of change for some of them 1. The same effect applies to non-coding features, which are annotated far less stably than protein-coding genes across releases, so lncRNA and circRNA results are especially sensitive to the build you chose 2.

Import into R with tximport, which sums transcript-level estimates to genes while carrying the effective-length correction into an offset that DESeq2 will use.

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

2. Build the object and write down the design before you look at results

Load the counts and sample sheet into a DESeqDataSet. The design formula is the analysis. Write it first, from what you know about the experiment, and resist the urge to revise it after seeing which version gives more hits.

library(DESeq2)
stopifnot(identical(colnames(txi$counts), samples$id))
samples$condition <- factor(samples$condition, levels = c("control", "treated"))
samples$batch     <- factor(samples$batch)

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

Order matters for readability, not for the fit: put the variable of interest last so that results() defaults to the right contrast. Set the reference level explicitly with factor(levels = ...), because R orders levels alphabetically and “control” versus “case” will silently give you a comparison in the direction you did not intend.

If your samples are repeated measures from one person or a handful of people, the design changes. With paired samples (same subject, two timepoints), add subject as a blocking factor: ~ subject + timepoint. With many timepoints per subject and a continuous covariate, you will do better with limma-voom and duplicateCorrelation(), which models within-subject correlation as a random effect rather than burning a coefficient per subject. edgeR’s generalized linear model framework handles the same class of designs through a design matrix you construct yourself, which is the route we take when the experiment involves nested or crossed factors that a simple formula cannot express 3.

3. Filter low-count genes and check the library sizes

Genes with near-zero counts across all samples contribute nothing but multiple-testing burden and unstable dispersion estimates. Filter on a count-per-million threshold scaled to the smallest group, which is what edgeR’s filterByExpr does; it is more principled than a flat “at least 10 reads in at least 3 samples” rule when library sizes vary.

library(edgeR)
keep <- filterByExpr(counts(dds), group = dds$condition)
dds  <- dds[keep, ]
nrow(dds)   # expect ~14,000-18,000 for bulk human whole blood or tissue

If that number comes back near 30,000 you have probably kept a large tail of noise; if it comes back near 9,000 your sequencing depth or RNA quality is low. Check library sizes with colSums(counts(dds)). Anything below about 10 million assigned reads for a standard bulk differential expression question is thin, and a sample sitting at a third of the median depth should be flagged, not silently normalized away.

4. Look at the structure of the data before testing anything

Transform to a variance-stabilized scale and plot. This is the step people skip and then regret.

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = c("condition", "batch"))

sampleDists <- dist(t(assay(vsd)))
pheatmap::pheatmap(as.matrix(sampleDists))

You are asking one question: what separates the samples on PC1 and PC2? If it is batch, extraction date, or RNA integrity number rather than your condition, the design formula must account for it. If nothing separates the samples, either the effect is small (common and fine) or the labels are wrong (check the sample sheet against a marker gene, for example XIST and RPS4Y1 to confirm sex).

For whole-blood RNA, cell composition dominates. Neutrophil fraction alone can drive tens of percent of the variance, and a condition that shifts the differential count moves thousands of genes without any change in per-cell regulation. Estimate proportions with a deconvolution method and include the dominant one or two as covariates, or accept that your result describes composition rather than transcriptional state and say so in writing.

5. Fit, test, and shrink

The fit itself is one line. DESeq2 estimates a size factor per sample by the median-of-ratios method, fits a gene-wise negative binomial dispersion, shrinks those dispersions toward a fitted trend, and performs a Wald test on the coefficient of interest.

dds <- DESeq(dds)
resultsNames(dds)
res <- results(dds, contrast = c("condition", "treated", "control"),
               alpha = 0.05, independentFiltering = TRUE)
summary(res)

Raw log2 fold changes for low-count genes are wild: a gene going from 2 reads to 12 reads reports a log2 fold change near 2.6 with an enormous standard error. Shrink them so the ranking is usable.

library(apeglm)
res_lfc <- lfcShrink(dds, coef = "condition_treated_vs_control", type = "apeglm")
plotMA(res_lfc, ylim = c(-4, 4))

If you care about a minimum effect size, test against it directly instead of filtering the output table by fold change afterward. The lfcThreshold argument tests the null hypothesis that the absolute log2 fold change is below your threshold, which is a real hypothesis test, whereas post-hoc filtering on fold change breaks the false discovery rate control you just paid for.

res_thr <- results(dds, contrast = c("condition","treated","control"),
                   lfcThreshold = 0.585, altHypothesis = "greaterAbs")

When would we use edgeR or limma-voom instead? edgeR’s quasi-likelihood F-test (glmQLFit then glmQLFTest) is our default when sample sizes are small, because it accounts for uncertainty in the dispersion estimate and tends to hold its nominal error rate better in that regime 3. limma-voom is our default above roughly 15 samples per group, because it is fast, and because once you need random effects or complicated contrasts the linear-model machinery is easier to reason about. For a two-group comparison with 4 to 10 samples per group, all three will agree on the top few hundred genes, and disagreement between them is itself a signal that the result is fragile.

6. Attach gene symbols and read the table properly

Gene identifiers from GENCODE carry version suffixes that will not match annotation databases, so strip them first.

library(AnnotationDbi); library(org.Hs.eg.db)
res_df <- as.data.frame(res_lfc)
res_df$ensembl <- sub("\\..*$", "", rownames(res_df))
res_df$symbol <- mapIds(org.Hs.eg.db, keys = res_df$ensembl,
                        column = "SYMBOL", keytype = "ENSEMBL",
                        multiVals = "first")
res_df <- res_df[order(res_df$padj), ]
write.csv(res_df, "de_results.csv")

Expect 10 to 20 percent of Ensembl IDs to map to no symbol, mostly novel transcripts, pseudogenes, and lncRNAs. That is normal and not a bug in your code. For the top hits, GeneCards is the fastest way to see what a gene is and what is already known about it, since it aggregates across more than a hundred source databases into one record per gene 4.

Read padj, not pvalue. A gene with p = 0.001 in a table of 16,000 tests is unremarkable on its own. Note also that padj will be NA for genes filtered out by independent filtering or flagged as count outliers by Cook’s distance, which is intentional: those genes were removed from the multiple-testing pool.

7. Move from a gene list to a hypothesis

A list of 400 significant genes is not an answer. Rank-based enrichment is more informative than a hypergeometric test on the significant subset, because it uses the whole ranking and does not depend on an arbitrary cutoff.

library(fgsea); library(msigdbr)
stats <- res_df$log2FoldChange / res_df$lfcSE
names(stats) <- res_df$symbol
stats <- sort(stats[!is.na(names(stats))], decreasing = TRUE)

paths <- msigdbr(species = "Homo sapiens", category = "H")
paths <- split(paths$gene_symbol, paths$gs_name)
fg <- fgsea(paths, stats, minSize = 15, maxSize = 500, nperm = 0)
head(fg[order(fg$padj), c("pathway","NES","padj","size")])

Rank by the Wald statistic or the t-statistic rather than by fold change alone, so that precision is part of the ordering. The output is a set of hypotheses, not conclusions. Published differential expression studies routinely treat their gene lists as the starting point for targeted follow-up with independent assays such as qPCR or protein measurement, and that is the right posture for personal data too 5. If any of this touches a health decision, the interpretation belongs with a clinician who can see your full picture; RNA expression from a single tissue at a single time is not a diagnostic.

8. Make the analysis reproducible

Record the software environment and the input hashes in the output, so that a result you look at in a year can be traced.

writeLines(capture.output(sessionInfo()), "sessionInfo.txt")

Pin your R and Bioconductor versions with renv, and run quantification inside a container with the annotation FASTA baked in. The reason is the one from step 1: annotation and genome builds shift, and a difference between two of your own analyses can come entirely from a reference update rather than from anything about you 1.

Common problems

Batch confounded with condition. If every treated sample was sequenced on run A and every control on run B, no statistical method can separate them. DESeq() will either drop the term or return a rank-deficient model matrix error. The fix is at the bench, not in R. When batch is partially confounded, you can include it, but the effective sample size drops and the standard errors widen accordingly.

Unmodeled hidden structure. When PC1 tracks nothing in your sample sheet, estimate surrogate variables with sva or RUVSeq (using empirically derived negative control genes) and add them to the design. Do this before testing, not after seeing that your favorite gene missed significance.

Too few replicates. With one sample per condition there is no way to estimate dispersion from the data. edgeR lets you set a biological coefficient of variation by hand (0.4 for human subjects, 0.1 for cell lines) and test against it, which produces a ranking but not calibrated p-values. Treat the output as exploratory. Longitudinal sampling of a single person is a better answer than a single pair: with 6 to 10 timepoints you can estimate within-person variability directly and ask whether a given timepoint departs from your own baseline.

Length-dependent bias in enrichment. Longer genes accumulate more reads and are therefore more likely to reach significance at a given fold change. If your enriched pathways are all made of very long genes, use goseq or an equivalent method that models the length bias.

Treating a fold change as an effect on you. A log2 fold change of 1.0 in whole blood means the average transcript abundance doubled across a mixture of cell types whose proportions may also have changed. Before interpreting it, check whether cell composition explains it, whether the gene is expressed at all in the relevant tissue, and whether the same direction appears at the protein level, which frequently it does not.

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. Anina N. Lund, Ryan C. Thompson, Ying-Chich Wang, et al. Updates to the reference genome alter the detection and direction of genes differentially expressed in Alzheimer’s disease. Genome Biology, 2026. https://doi.org/10.1186/s13059-026-04213-9 ↩ ↩2

  2. Min Li, Jian Wang, Dewu Liu, et al. High‑throughput sequencing reveals differentially expressed lncRNAs and circRNAs, and their associated functional network, in human hypertrophic scars. Molecular Medicine Reports, 2018. https://doi.org/10.3892/mmr.2018.9557 ↩

  3. Yunshun Chen, Aaron T. L. Lun, Gordon K. Smyth. Differential Expression Analysis of Complex RNA-seq Experiments Using edgeR. Statistical Analysis of Next Generation Sequencing Data, 2014. https://doi.org/10.1007/978-3-319-07212-8_3 ↩ ↩2

  4. Gil Stelzer, Naomi Rosen, Inbar Plaschkes, et al. The GeneCards Suite: From Gene Data Mining to Disease Genome Sequence Analyses. Current Protocols in Bioinformatics, 2016. https://doi.org/10.1002/cpbi.5 ↩

  5. Fang Miao, Xixi Li, Chenglin Wang, et al. Bioinformatics analysis of differentially expressed genes in diabetic foot ulcer and preliminary experimental verification. Annals of Translational Medicine, 2023. https://doi.org/10.21037/atm-22-6437 ↩