How to Build an RNA-seq Pipeline in R for Your Own Data

By the end of this guide you will have a reproducible R project that turns raw FASTQ files into four artifacts: a gene-level count matrix with effective-length offsets, a quality-control report that tells you which samples to trust, a table of genes whose expression changes across your timepoints with shrunken effect sizes and adjusted p-values, and a set of pathway-level and cell-composition summaries that are far more stable than single genes. You need the compressed FASTQ files and their checksums, a sample sheet, roughly 60 GB of free disk and 32 GB of RAM for a human transcriptome index, R 4.4 or later with a current Bioconductor release, and a command line for the quantification step. The whole run takes a couple of hours of compute for a dozen samples, most of it in quantification, and considerably longer in the thinking that follows.

A note on framing before the code. RNA sequencing measures the relative abundance of transcripts in whatever tissue was sampled, usually whole blood or PBMCs for a personal profile. It tells you what a set of cells was transcribing at the moment the tube was drawn, which makes it the most state-dependent modality in a molecular profile and the one most sensitive to sampling conditions. Nothing you compute here is a diagnosis. If a result points at something clinically meaningful, that conversation belongs with a physician who can order a validated assay.

1. Organize inputs and write the sample sheet first

Start with the metadata, because every modeling decision downstream is a column in this table. Create a samples.csv with one row per library and columns for sample ID, subject, collection date and time of day, tube type (PAXgene, Tempus, EDTA), library prep kit, whether globin depletion was performed, ribodepletion versus polyA selection, sequencing run, lane, and read length. Time of day matters more than people expect, since a large block of the blood transcriptome is under circadian control, and a 7 a.m. draw is not comparable to a 4 p.m. draw without accounting for it.

Verify file integrity before you spend compute on corrupted input:

md5sum -c checksums.md5
zcat S01_R1.fastq.gz | head -n 4
seqkit stats -a *.fastq.gz > fastq_stats.tsv

Then run FastQC and aggregate with MultiQC. What you are looking for is adapter content above a few percent, per-base quality collapsing after cycle 100, and duplication curves that suggest a low-complexity library. Best-practice surveys of RNA-seq analysis are consistent that library quality and design dominate the outcome far more than the choice of downstream statistical package.¹

2. Quantify with Salmon and a decoy-aware index

We quantify outside R. Alignment-free or selective-alignment quantifiers are faster, require less memory than a genome aligner, and give transcript-level abundances with effective-length correction that R can consume directly. We use Salmon with a decoy-aware index built from the genome, which prevents reads from unannotated genomic regions being forced onto transcripts they do not belong to.

# Build the gentrome: transcriptome + genome as decoys
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) \
  | cut -d " " -f 1 | sed 's/>//g' > decoys.txt
cat gencode.v46.transcripts.fa.gz GRCh38.primary_assembly.genome.fa.gz > gentrome.fa.gz

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_idx_k31 \
  -k 31 -p 16 --gencode

Use -k 31 for 75 bp reads or longer and drop to -k 23 if your reads are shorter than about 50 bp. Then quantify each sample:

salmon quant -i salmon_idx_k31 -l A \
  -1 S01_R1.fastq.gz -2 S01_R2.fastq.gz \
  --gcBias --seqBias --posBias \
  --numGibbsSamples 20 --thinningFactor 16 \
  -p 16 --validateMappings -o quants/S01

-l A lets Salmon infer strandedness, which you should then confirm in lib_format_counts.json rather than trusting silently. --gcBias and --seqBias correct fragment-level biases that otherwise show up as apparent expression differences between batches. The Gibbs samples give you posterior uncertainty per transcript, which matters if you later work at transcript rather than gene level. Record the mapping rate from logs/salmon_quant.log: for polyA human libraries we expect 75 to 90 percent, and anything below 60 percent usually means contamination, wrong index, or degraded RNA.

3. Import into R with tximport

tximport converts Salmon’s transcript-level output into a gene-level matrix plus the average transcript-length offsets that make gene-level counts comparable across samples. Building tx2gene from the same GTF you used to build the index is not optional, because version suffixes and annotation releases must match.

library(tximport); library(txdbmaker); library(readr)

txdb  <- makeTxDbFromGFF("gencode.v46.annotation.gtf.gz")
k     <- keys(txdb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(txdb, k, "GENEID", "TXNAME")

samples <- read_csv("samples.csv")
files   <- file.path("quants", samples$sample_id, "quant.sf")
names(files) <- samples$sample_id
stopifnot(all(file.exists(files)))

txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                ignoreTxVersion = FALSE)

Leave countsFromAbundance at its default ("no") when the destination is DESeq2, which will use the length matrix as an offset. If you plan to use limma-voom instead, re-import with countsFromAbundance = "lengthScaledTPM" and feed the counts directly, since limma has no offset slot in the same sense.

4. Build the DESeqDataSet and filter

DESeq2 is our default for a personal dataset with a modest number of samples. It models counts with a negative binomial distribution and shares dispersion information across genes, which is what makes inference possible when you have five or ten libraries rather than fifty. limma-voom is the better choice above roughly twenty samples per group because it is faster and its linear-model machinery handles complex designs more gracefully, but with few replicates DESeq2’s shrinkage is more forgiving.¹

library(DESeq2)

coldata <- as.data.frame(samples)
rownames(coldata) <- coldata$sample_id
coldata$timepoint <- factor(coldata$timepoint)
coldata$batch     <- factor(coldata$seq_run)

dds <- DESeqDataSetFromTximport(txi, colData = coldata,
                                design = ~ batch + timepoint)

# Filter: keep genes with >=10 counts in at least as many samples
# as the smallest group
smallest <- min(table(coldata$timepoint))
keep <- rowSums(counts(dds) >= 10) >= smallest
dds  <- dds[keep, ]
nrow(dds)   # expect ~14,000-18,000 for whole blood polyA

Filtering to 14,000 or so expressed genes out of 62,000 annotated features is normal and it improves power by reducing the multiple-testing burden. If you retain only 9,000, your library is probably shallow or globin-dominated. Keep in mind that the counts are compositional: they describe proportions of a fixed sequencing budget, so an apparent decrease in one gene can reflect a genuine increase elsewhere. The size-factor normalization DESeq2 applies assumes most genes are unchanged, and that assumption deserves a second look whenever one transcript family takes over a library.²

5. Quality control on the variance-stabilized matrix

Never do exploratory work on raw or log-CPM counts with pseudocounts, because low-count genes will dominate the distance structure. Use the variance-stabilizing transform, which flattens the mean-variance relationship.

vsd <- vst(dds, blind = FALSE)

plotPCA(vsd, intgroup = c("timepoint", "batch"))

# Sample-sample distances
sampleDists <- dist(t(assay(vsd)))
pheatmap::pheatmap(as.matrix(sampleDists),
                   clustering_distance_rows = sampleDists,
                   clustering_distance_cols = sampleDists)

# Globin and mitochondrial fraction per sample
globin <- c("HBB","HBA1","HBA2","HBD","HBG1","HBG2")
frac <- colSums(counts(dds, normalized = FALSE)[
          rownames(dds) %in% gene_ids_for(globin), ]) /
        colSums(counts(dds, normalized = FALSE))
round(frac, 3)

Read PC1 carefully. In whole blood from a single person sampled over months, PC1 is very often globin fraction, neutrophil proportion, or sequencing batch rather than anything you care about. If batch separates cleanly, use limma::removeBatchEffect on the VST matrix for plotting only, and keep batch in the design matrix for the actual model. Correcting the matrix and then testing on the corrected values understates standard errors.

For latent factors you did not record, run sva::svaseq on the normalized counts with your known design as the full model and an intercept-only reduced model, then add the resulting surrogate variables as covariates. Two surrogate variables is usually enough for a dozen samples, and adding more will consume the residual degrees of freedom you need.

6. Choose a model that matches how you sampled

This is the step where personal data differs most from a textbook experiment. With one subject, you have no between-person variance and every comparison is within-person across time or condition, which is statistically favorable because genetic background and most stable environmental factors are held constant.

For two discrete states, the design in step 4 is sufficient and you extract a contrast. For a dense time course, fit a smooth term and use a likelihood ratio test that asks whether any temporal structure exists:

library(splines)
design(dds) <- ~ batch + ns(day, df = 3)
dds <- DESeq(dds, test = "LRT", reduced = ~ batch)
res <- results(dds)
sum(res$padj < 0.05, na.rm = TRUE)

Set df from the number of timepoints, not from hope: with six timepoints, df = 3 is already generous. Time-course transcriptomics gains most of its interpretive value when trajectories are clustered and then mapped onto regulators, and tools built for that purpose group genes by temporal pattern before asking which transcription factors explain each cluster.³ If you collected several samples per day, treat day as a random effect using limma’s duplicateCorrelation rather than pretending the replicates are independent.

7. Extract results with shrunken effect sizes

Raw log2 fold changes for low-count genes are noisy and will mislead any ranking you build on them. Shrink them.

resultsNames(dds)
res <- lfcShrink(dds, coef = "timepoint_T2_vs_T1", type = "apeglm")

res_df <- as.data.frame(res) |>
  tibble::rownames_to_column("gene_id") |>
  dplyr::arrange(padj)

# Map to symbols
library(AnnotationDbi); library(org.Hs.eg.db)
res_df$symbol <- mapIds(org.Hs.eg.db,
  keys = sub("\\..*", "", res_df$gene_id),
  column = "SYMBOL", keytype = "ENSEMBL", multiVals = "first")

write_csv(res_df, "results/de_T2_vs_T1.csv")

Use type = "apeglm" for contrasts expressible as a coefficient and type = "ashr" when you need an arbitrary contrast. Interpret a gene only if it survives the adjusted p-value threshold and has a shrunken absolute log2 fold change above roughly 0.5, meaning a 1.4-fold change. Single-gene results from a handful of libraries are hypotheses, not findings.

8. Move to pathways, where the signal is more stable

Individual genes are noisy across draws. Coordinated sets of genes are not, which is why we spend most interpretive effort at the pathway level. Rank by the Wald or LRT statistic rather than fold change, because the statistic already incorporates precision.

library(fgsea); library(msigdbr)

h <- msigdbr(species = "Homo sapiens", collection = "H")
pathways <- split(h$gene_symbol, h$gs_name)

stats <- res_df$stat
names(stats) <- res_df$symbol
stats <- stats[!is.na(names(stats)) & !is.na(stats)]
stats <- stats[!duplicated(names(stats))]

fg <- fgseaMultilevel(pathways, stats, minSize = 15, maxSize = 500)
fg[order(fg$padj), c("pathway","NES","padj","size")][1:20, ]

For per-sample pathway scores that you can plot as a trajectory over months, use GSVA::gsva with the ssgseaParam or gsvaParam constructor on the VST matrix. Those scores behave like continuous biomarkers and are the natural unit for joining RNA-seq to proteomic and metabolic measurements collected on the same day.

9. Deconvolve blood cell composition

A whole-blood expression change is ambiguous between a shift in what cells are doing and a shift in which cells are present. Estimating cell fractions from the same matrix resolves much of that ambiguity.

library(immunedeconv)
tpm <- txi$abundance[keep, ]
rownames(tpm) <- res_symbols   # gene symbols, deduplicated
frac <- deconvolute(tpm, "quantiseq")

Run the deconvolution before you interpret anything immune-related, and then re-fit your differential expression with the dominant cell fractions as covariates to ask whether a signal survives. Neutrophil percentage alone explains a large share of variance in whole-blood transcriptomes, and it moves with acute illness, exercise, and stress.

10. Join RNA to the rest of your profile

Expression becomes far more interpretable next to your genome and your proteome. If you carry a variant of interest, expression of the affected gene and its neighbors tells you whether the variant has a transcriptional consequence in the tissue you sampled, which is the logic behind expression and protein quantitative trait loci studies that connect a variant to a transcript and then to a circulating protein.⁴ Practically, join your DESeq2 result table to your annotated variant table on gene ID and inspect allele-specific expression at heterozygous sites using phASER or WASP-corrected alignments. Statistical frameworks for integrating modalities measured on the same samples are an active area, and the general lesson is to integrate at the level of latent factors or pathways rather than trying to align raw features.⁵

Finally, make the whole thing reproducible. Pin package versions with renv::snapshot(), wrap the steps in a targets pipeline so that changing the sample sheet re-runs only what depends on it, and save sessionInfo() alongside every result table. A year from now, when you add the next timepoint, you will want to re-run the identical code rather than reconstruct it.

Common problems

Low mapping rate, under 60 percent, almost always traces to one of four causes: an index built from a different annotation release than your tx2gene, ribosomal RNA carryover from failed depletion, genomic DNA contamination, or bacterial or adapter content. Check lib_format_counts.json and run a quick screen with fastq_screen against human, rRNA, and common contaminants before rebuilding anything.

Globin transcripts dominating whole blood. Without globin depletion, hemoglobin transcripts can consume a large fraction of a whole-blood library, leaving too few reads for everything else. You can filter globin genes before computing size factors, but this only partially recovers sensitivity, and the real fix is globin depletion at library prep. Record depletion status in your sample sheet so you never compare a depleted library to an undepleted one without a covariate.

Batch confounded with time. If every T1 sample was sequenced in run A and every T2 in run B, no statistical method can separate them, and DESeq will either error on a rank-deficient model matrix or give you results that mean nothing. Randomize samples across runs, or retain aliquots and sequence a bridge sample in every run.

Too few replicates. Two samples per condition gives you very little power, and reported fold changes will be inflated for the genes that reach significance. Treat such results as a ranked hypothesis list, prioritize pathways over genes, and resist the temptation to read a single dramatic gene as meaningful.¹

Degraded RNA. A low RNA integrity number produces a 3’ bias that mimics differential expression of long transcripts. Check the gene-body coverage plot in your MultiQC report, and if integrity varies across your timepoints, include it as a covariate rather than ignoring it.

Over-normalizing. If one condition genuinely changes total RNA output per cell, median-ratio size factors will absorb that change and hide it. Spike-in controls added at a fixed amount per cell are the only clean solution, and in their absence you should state plainly that your measurements are relative abundances.²

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. Ana Conesa, Pedro Madrigal, Sonia Tarazona, et al. A survey of best practices for RNA-seq data analysis. Genome Biology, 2016. https://doi.org/10.1186/s13059-016-0881-8 ↩ ↩² ↩³

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

3. Ashley Mae Conard, Nathaniel Goodman, Yanhui Hu, et al. TIMEOR: a web-based tool to uncover temporal regulatory mechanisms from multi-omics data. Nucleic Acids Research, 2021. https://doi.org/10.1093/nar/gkab384 ↩

4. Junxiang Lian, Xinjian Zhang, Shanwei Shi, et al. Multi-omics Mendelian randomization integrating RNA-seq, eQTL and pQTL data revealed CPXM1 as a potential drug target for osteoporosis. Hereditas, 2025. https://doi.org/10.1186/s41065-025-00562-w ↩

5. M. Colomé-Tatché, F.J. Theis. Statistical single cell multi-omics integration. Current Opinion in Systems Biology, 2018. https://doi.org/10.1016/j.coisb.2018.01.003 ↩