Analyzing Your Own RNA-Seq Data in R

By the end of this guide you will have an R project that reads transcript-level quantifications into gene-level counts, filters them by a defensible rule, checks the samples for technical problems, fits a model appropriate to repeated draws from a single person, and produces a ranked gene list plus pathway enrichment you can inspect. You need: quantification output from a transcript-level quantifier (salmon quant.sf files are assumed here), a sample sheet with draw date, time of day, and anything else you recorded, R 4.4 or later with Bioconductor 3.20, roughly 16 GB of RAM for a few dozen samples, and disk space for the transcriptome index if you are quantifying yourself. Everything below is measurement and interpretation. Nothing here is a clinical result, and any finding you would act on belongs in front of a physician who can order confirmatory testing.

1. Get to transcript quantifications before you open R

R is the wrong tool for read alignment. Do the FASTQ-to-counts step in a workflow manager and treat its output as your starting point. The most reliable path is nf-core/rnaseq, which wraps trimming, quantification, and QC with pinned container versions:

nextflow run nf-core/rnaseq -r 3.14.0 \
  --input samplesheet.csv \
  --outdir results \
  --genome GRCh38 \
  --pseudo_aligner salmon \
  --skip_alignment \
  -profile docker

If you would rather run salmon directly, build a decoy-aware index so that reads from unannotated genomic regions do not get forced onto transcripts:

grep "^>" GRCh38.primary_assembly.genome.fa | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v45.transcripts.fa GRCh38.primary_assembly.genome.fa > gentrome.fa
salmon index -t gentrome.fa -d decoys.txt -p 16 -i salmon_idx_v45 --gencode

salmon quant -i salmon_idx_v45 -l A \
  -1 sample_R1.fastq.gz -2 sample_R2.fastq.gz \
  -p 16 --gcBias --seqBias --validateMappings \
  --numBootstraps 30 -o quants/sample

The three flags that matter are --gcBias and --seqBias, which correct for fragment GC content and random-hexamer priming bias, and --numBootstraps, which gives you per-transcript inferential replicates. Bootstraps cost little and are the only way to know whether a transcript-level change is real or an artifact of multimapping reads being shuffled between isoforms. Record the mapping rate from quants/sample/logs/salmon_quant.log for every sample: below about 70% for a poly-A blood library usually means contamination, adapter carryover, or the wrong index.

One sample-handling note, because it changes what you can do later. If the same blood or tissue aliquot is being used for metabolomics, the extraction order is a real variable. A controlled comparison in mouse liver found that extracting metabolites first, then RNA from the same sample, preserved RNA integrity and sequencing quality well enough to support integrated analysis ¹. That is useful when your material is limited and you want the transcriptome and the metabolome measured on identical biology rather than on two different draws.

2. Set up the project so it still runs in a year

Use renv from the first line. Bioconductor packages change APIs often enough that an un-pinned analysis will not reproduce.

install.packages("renv"); renv::init()
install.packages("BiocManager")
BiocManager::install(c("tximport", "DESeq2", "apeglm", "limma", "edgeR",
                       "AnnotationHub", "fgsea", "vsn", "sva"))
renv::snapshot()

Keep the analysis in numbered scripts (01_import.R, 02_qc.R, 03_de.R) that each write an .rds to a cache/ directory, rather than one long notebook. When you rerun with three new timepoints, you want to invalidate one step, not all of them.

3. Import with tximport and build the count object

tximport aggregates transcript estimates to genes and, importantly, computes an average transcript-length offset per gene per sample. That offset corrects for the case where a gene’s expression looks unchanged but its dominant isoform switched to a shorter one, which would otherwise register as a spurious change in counts.

library(tximport); library(AnnotationHub); library(DESeq2)

meta <- read.csv("samplesheet.csv", stringsAsFactors = FALSE)
meta$draw_date <- as.Date(meta$draw_date)
meta$clock_hour <- as.numeric(meta$draw_time_hhmm) / 100
files <- file.path("quants", meta$sample_id, "quant.sf")
names(files) <- meta$sample_id
stopifnot(all(file.exists(files)))

ah <- AnnotationHub()
edb <- ah[["AH113665"]]   # EnsDb, Homo sapiens 110 — verify the ID for your release
tx2gene <- transcripts(edb, return.type = "data.frame")[, c("tx_id", "gene_id")]

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                ignoreTxVersion = TRUE, countsFromAbundance = "no")
dds <- DESeqDataSetFromTximport(txi, colData = meta, design = ~ 1)

Set ignoreTxVersion = TRUE only if your GTF and your index come from the same Ensembl or GENCODE release with differing version suffixes. If the annotation releases genuinely differ, fix that instead of silencing the mismatch. tximport will warn about transcripts missing from tx2gene; a loss of more than about 1% of transcripts means you have the wrong annotation file.

4. Filter low-count genes with a rule you can state

The default habit of keeping genes with more than ten reads total is too loose and varies with library size. We prefer the relevance filter from the omics data analysis framework for regulatory application (R-ODAF), which keeps a gene only if its counts per million exceed 1 in at least 75% of the samples of at least one experimental group ². It is explicit, scales with sequencing depth, and does not quietly retain genes that are zero everywhere except one sample.

library(edgeR)
cpm_mat <- cpm(counts(dds))
groups  <- split(colnames(dds), dds$condition)   # or a single group for an n=1 series
keep <- Reduce(`|`, lapply(groups, function(s)
          rowMeans(cpm_mat[, s, drop = FALSE] > 1) >= 0.75))
dds <- dds[keep, ]
nrow(dds)

For whole blood on a 30–40 million read library this typically leaves 12,000 to 15,000 genes. The same framework also specifies post-hoc checks on the resulting gene lists, including comparing the observed number of significant genes against what the false discovery rate alone would produce, which is a sensible guard when you are working with few samples ².

Before moving on, look at globin. In non-depleted whole blood, HBB, HBA1, HBA2, and HBD can consume a large share of reads, and that share varies between draws. Quantify it (colSums(counts(dds)[globin_ids, ]) / colSums(counts(dds))), and if it swings widely, remove those genes and recompute size factors so that a globin excursion does not shift the normalization for every other gene.

5. Transform, then look at the samples before modeling anything

Variance-stabilizing transformation gives you a matrix on a roughly homoscedastic scale, which is what principal components and clustering need.

vsd <- vst(dds, blind = TRUE)
plotPCA(vsd, intgroup = c("draw_date", "clock_hour"))

library(pheatmap)
d <- dist(t(assay(vsd)))
pheatmap(as.matrix(d), clustering_distance_rows = d, clustering_distance_cols = d)

For a personal longitudinal series, PC1 is frequently not biology. It is usually library prep batch, RNA integrity, or time of day. Cortisol-driven leukocyte trafficking makes the 7 a.m. blood transcriptome measurably different from the 4 p.m. one in the same person, so record draw time and carry it as a covariate. Also plot meanSdPlot(assay(vsd)) from the vsn package: a strongly rising trend at low means says your filter was too permissive.

Outlier handling deserves a rule set in advance. We flag any sample whose median Spearman correlation to all others falls more than three median absolute deviations below the group median, then inspect its salmon mapping rate and RIN before deciding. Dropping samples after seeing the differential expression result is how you get findings that do not replicate.

6. Choose a model that matches an n=1 design

This is where guides written for two-group lab experiments stop being useful. With one person, there is no between-subject replication, and the question changes from “which genes differ between groups” to “which genes move systematically over time, or with a covariate I measured.”

If you have enough timepoints (we would want at least eight), model time as a smooth term and let the residual degrees of freedom serve as your replication:

library(splines)
dds$days <- as.numeric(dds$draw_date - min(dds$draw_date))
design(dds) <- ~ clock_hour + ns(days, df = 3)
dds <- DESeq(dds, test = "LRT", reduced = ~ clock_hour)
res_time <- results(dds)

The likelihood ratio test here asks whether the spline in time explains variance beyond the time-of-day covariate. With df = 3 and, say, twelve samples you retain seven residual degrees of freedom, which is thin but workable because DESeq2 shrinks dispersion estimates toward a fitted trend across genes.

If instead you have a defined before/after contrast with repeated draws in each state, use limma-voom with duplicateCorrelation so that draws from the same period are not treated as independent:

library(limma)
dge <- DGEList(counts(dds)); dge <- calcNormFactors(dge)
mm  <- model.matrix(~ clock_hour + period, data = as.data.frame(colData(dds)))
v   <- voom(dge, mm, plot = TRUE)
corfit <- duplicateCorrelation(v, mm, block = dds$draw_week)
v   <- voom(dge, mm, block = dds$draw_week, correlation = corfit$consensus)
fit <- eBayes(lmFit(v, mm, block = dds$draw_week, correlation = corfit$consensus))
topTable(fit, coef = "periodafter", n = 30)

Two practical warnings. First, with fewer than four samples per state, neither method produces a fold change you should trust for a single gene. Read the results as a ranked list for enrichment, not as per-gene claims. Second, if the covariate you care about is perfectly confounded with sequencing batch, no model can separate them. Randomize sample-to-batch assignment when you submit, or hold a bridging aliquot across batches.

7. Shrink effect sizes and export a result table

Raw log2 fold changes for low-count genes are enormous and meaningless. Shrink them:

resLFC <- lfcShrink(dds, coef = "periodafter_vs_before", type = "apeglm")
out <- as.data.frame(resLFC)
out$symbol <- mapIds(edb, rownames(out), "SYMBOL", "GENEID")
write.csv(out[order(out$padj), ], "results/de_period.csv")

Report the shrunken log fold change and the adjusted p-value together, and keep the unshrunken table as well so the ranking used downstream is reproducible. For ranking genes into enrichment, we use the signed test statistic (stat from DESeq2 or t from limma) rather than fold change, because it incorporates the precision of each estimate.

8. Ask what the blood cell composition did

A large fraction of apparent transcriptome change in whole blood is cell composition rather than per-cell regulation. A neutrophil shift of a few percent moves hundreds of genes. Deconvolution packages such as granulator or immunedeconv estimate proportions from your bulk expression using a reference signature matrix, and those proportions can go into the model as covariates.

Treat the estimates as approximate. Benchmarking of cell type annotation methods on single-cell data has shown that accuracy varies substantially with the method and the reference used, and bulk deconvolution inherits the same dependency on reference choice ³. If a result survives adjustment for estimated composition, it is more interesting than one that does not. If you have a same-day complete blood count with differential, use those measured proportions instead of inferred ones.

9. Move from genes to pathways

Single genes from a small design are fragile. Pathway-level signal is more stable because it aggregates across correlated genes. Run a ranked enrichment rather than an over-representation test on a thresholded list, since thresholding throws away information you paid for.

library(fgsea)
ranks <- setNames(out$stat, out$symbol)
ranks <- sort(ranks[!is.na(ranks) & !is.na(names(ranks))], decreasing = TRUE)
gs <- gmtPathways("c2.cp.reactome.v2024.1.Hs.symbols.gmt")
fg <- fgseaMultilevel(gs, ranks, minSize = 15, maxSize = 500)
head(fg[order(fg$padj), c("pathway","NES","padj","size")], 20)

Gene set size bounds matter: sets below 15 genes are noisy and sets above 500 are too diffuse to interpret. For organizing results afterward, the g:Profiler, GSEA, Cytoscape, and EnrichmentMap protocol is the clearest published route from a ranked list to a network of related gene sets, which collapses the hundred redundant Reactome terms you will get into a handful of themes ⁴.

When you have several contrasts at once, or transcriptome and proteome measured on the same samples, mitch is the tool we reach for. It performs multi-contrast enrichment by ranking genes in each dimension and testing gene sets for coordinated movement across dimensions, so you can distinguish sets that move together in RNA and protein from sets that move in only one ⁵. That distinction is often more informative than either assay alone, because concordance argues against an assay-specific artifact. Integrative multi-omics designs of this kind have identified pathway-level programs that were not apparent from transcriptomics by itself ⁶.

Common problems

Low mapping rate with high duplication usually means a low-input library that was over-amplified. Check the salmon log alongside the duplication estimate from fastqc, and do not try to rescue it computationally.

A PCA where samples separate perfectly by sequencing run tells you batch is confounded. sva::ComBat_seq can adjust counts, but if batch and your variable of interest are collinear, the adjustment will remove your signal along with the batch effect. Check table(batch, condition) before you run anything.

Thousands of significant genes from a tiny design is a red flag, not a success. Compare your count of significant genes against a permuted-label run of the same model. The R-ODAF post-hoc checks formalize this comparison ².

tximport failing with “transcripts missing from tx2gene” almost always means mismatched annotation versions between the index and the EnsDb. Rebuild the index from the same release rather than dropping the unmatched transcripts.

Genes that look dramatic but are all mitochondrial or all ribosomal protein genes point to RNA degradation or a library prep difference, not biology. Compute the mitochondrial read fraction per sample and put it in the model if it varies.

Finally, remember what a single person’s transcriptome can and cannot tell you. RNA-seq measures abundance across the transcriptome at high dynamic range and can resolve isoform usage and allele-specific expression when depth allows ⁷. It does not tell you why a transcript moved, and it is not a diagnostic test. Interesting findings are hypotheses for a next measurement or for a conversation with a clinician who can order the appropriate assay.

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. Zachary B Madaj, Michael S. Dahabieh, Vijayvardhan Kamalumpundi, et al. Prior metabolite extraction fully preserves RNAseq quality and enables integrative multi-‘omics analysis of the liver metabolic response to viral infection. RNA Biology, 2023. https://doi.org/10.1080/15476286.2023.2204586 ↩

2. Yuko Ogata. Applying the omics data analysis framework for regulatory application (R-ODAF) filtering criteria: A case study with 3-methylcholanthrene. Regulatory Toxicology and Pharmacology, 2026. https://doi.org/10.1016/j.yrtph.2026.106175 ↩ ↩² ↩³

3. Qianhui Huang, Yu Liu, Yuheng Du, et al. Evaluation of Cell Type Annotation R Packages on Single-Cell RNA-Seq Data. Genomics, Proteomics & Bioinformatics, 2020. https://doi.org/10.1016/j.gpb.2020.07.004 ↩

4. Jüri Reimand, Ruth Isserlin, Veronique Voisin, et al. Pathway enrichment analysis and visualization of omics data using g:Profiler, GSEA, Cytoscape and EnrichmentMap. Nature Protocols, 2019. https://doi.org/10.1038/s41596-018-0103-9 ↩

5. Antony Kaspi, Mark Ziemann. mitch: multi-contrast pathway enrichment for multi-omics and single-cell profiling data. BMC Genomics, 2020. https://doi.org/10.1186/s12864-020-06856-9 ↩

6. Vincent Gureghian, Hailee Herbst, Ines Kozar, et al. A multi-omics integrative approach unravels novel genes and pathways associated with senescence escape after targeted therapy in NRAS mutant melanoma. Cancer Gene Therapy, 2023. https://doi.org/10.1038/s41417-023-00640-z ↩

7. 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 ↩