How to Analyze Your Own Bulk RNA-seq Data, Start to Finish

Oak
rna-seq transcriptomics bioinformatics salmon deseq2
A long laboratory instrument passing a glowing thread through six glass chambers where it resolves into separate beads of light.

By the end of this guide you will have five artifacts from your raw sequencing files: a per-sample quality report, a gene-level count matrix with transcript-level abundance estimates, a variance-stabilized expression matrix suitable for clustering and visualization, a ranked table of genes that change across your timepoints or conditions, and an estimate of the immune cell composition driving those changes. You need the raw FASTQ files (paired-end, gzipped), roughly 60 GB of free disk per 10 samples, 32 GB of RAM, 8 or more cores, and a working conda or mamba installation. Everything below runs on a laptop with enough disk, though a 16-core machine will take a 30-million-read sample from FASTQ to counts in about four minutes rather than twenty. RNA-seq is a research measurement, not a clinical test; nothing produced by this pipeline is a diagnosis, and anything that looks medically consequential belongs in front of a physician.

1. Write down what was sequenced before you touch a FASTQ

Nearly every downstream decision depends on library chemistry, and guessing wastes hours. Find the library prep protocol and record six things: whether the library was polyA-selected or ribosomal-RNA-depleted, whether globin mRNA was depleted (relevant for any whole-blood sample), read length and whether reads are paired, strandedness, whether unique molecular identifiers (UMIs) were used, and the target depth. For whole blood, a polyA library without globin depletion routinely gives up 50 to 70 percent of reads to HBA1, HBA2, and HBB, which means a nominal 30-million-read library behaves like a 10-million-read one for everything else.

Strandedness is the single most common silent error. A stranded library analyzed as unstranded inflates antisense noise and shifts counts for overlapping gene pairs. You can infer it empirically rather than trust the datasheet, and step 4 does exactly that.

2. Build the environment and a decoy-aware index

We recommend GENCODE as the annotation source because the transcript FASTA and GTF are guaranteed consistent with each other, which avoids a whole class of identifier mismatches later.

mamba create -n rnaseq -c bioconda -c conda-forge \
  salmon=1.10.3 fastp=0.23.4 fastqc multiqc samtools=1.20 \
  bioconductor-tximport bioconductor-deseq2 r-base=4.4
conda activate rnaseq

REF=~/ref/gencode_v47
mkdir -p $REF && cd $REF
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/GRCh38.primary_assembly.genome.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_47/gencode.v47.primary_assembly.annotation.gtf.gz

# decoy set: genome sequence names, so genomic reads don't get forced onto transcripts
grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) | cut -d " " -f1 | sed 's/>//g' > decoys.txt
cat gencode.v47.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_idx_v47 --gencode

The --gencode flag strips the pipe-delimited header fields so transcript IDs come out clean. The decoy genome matters: without it, intronic and intergenic reads (which are abundant in any library made from partially degraded or pre-mRNA-rich material) get assigned to whatever transcript they resemble most, inflating a subset of genes reproducibly enough to look real.

3. Trim and quality-check with fastp

fastp -i sample_R1.fastq.gz -I sample_R2.fastq.gz \
      -o trim_R1.fq.gz -O trim_R2.fq.gz \
      --detect_adapter_for_pe \
      --qualified_quality_phred 15 --unqualified_percent_limit 40 \
      --length_required 36 --trim_poly_g \
      --thread 8 --json sample.fastp.json --html sample.fastp.html

Do not aggressively quality-trim. Pseudoalignment tolerates mismatches well, and hard trimming at Q30 removes real signal while changing effective fragment lengths. The --trim_poly_g flag matters on NovaSeq and NextSeq two-color chemistry, where a dark base is read as G and unclipped polyG tails depress mapping rates by several percent.

Run FastQC and aggregate with multiqc .. Four numbers are worth reading carefully: duplication rate (high duplication with low library complexity suggests over-amplification of a low-input sample), adapter content after trimming (should be near zero), per-sequence GC distribution (a sharp secondary peak often means rRNA or adapter dimer), and overrepresented sequences (in blood, expect hemoglobin). Comprehensive per-sample QC before any statistical modeling is the step that catches most artifacts, and recent quality control frameworks make the same argument for multi-omic data generally: flag samples on multiple independent metrics rather than one summary score.1

4. Quantify with salmon

salmon quant -i $REF/salmon_idx_v47 -l A \
  -1 trim_R1.fq.gz -2 trim_R2.fq.gz \
  -p 12 --seqBias --gcBias --posBias \
  --numGibbsSamples 20 --validateMappings \
  -o quant/sample

-l A asks salmon to infer library type; check quant/sample/lib_format_counts.json afterward and confirm the inferred type (ISR for a typical dUTP stranded library, IU for unstranded) is consistent across every sample. A single sample disagreeing with the rest means a library prep or demultiplexing problem, not a biological one.

--seqBias and --gcBias model random-hexamer priming bias and fragment GC bias respectively, both of which vary between batches and will otherwise masquerade as differential expression. --numGibbsSamples 20 produces posterior inferential replicates, which let downstream tools distinguish genes with genuinely uncertain quantification (short transcripts, paralog families like the HLA and immunoglobulin loci) from confidently measured ones.

Expect a mapping rate of 75 to 90 percent for a good human library. Below 60 percent, stop and diagnose rather than proceed: check the index species, check for rRNA carryover, and check adapter trimming. Use STAR instead of salmon when you need something quantification cannot give you: splice junction counts, allele-specific expression, fusion detection, or variant calls from RNA. Those questions require an actual genomic alignment, and STAR --quantMode GeneCounts TranscriptomeSAM plus salmon quant --alignments gives you both at the cost of roughly 32 GB RAM for the human index.

5. Import to gene level with tximport

library(tximport); library(readr); library(GenomicFeatures)

txdb <- makeTxDbFromGFF("~/ref/gencode_v47/gencode.v47.primary_assembly.annotation.gtf.gz")
k    <- keys(txdb, keytype = "TXNAME")
tx2gene <- AnnotationDbi::select(txdb, k, "GENEID", "TXNAME")

files <- file.path("quant", samples$id, "quant.sf")
names(files) <- samples$id
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                countsFromAbundance = "lengthScaledTPM")

Use lengthScaledTPM when you plan to feed the matrix to limma-voom, and plain counts with the offset (countsFromAbundance = "no") when using DESeq2 via DESeqDataSetFromTximport, which handles the average transcript length offset itself. That offset is what corrects for differential isoform usage changing a gene’s effective length between samples, and skipping it is a real source of false positives in genes with divergent isoform lengths.

6. Filter, normalize, and look at structure before testing anything

Filtering low-count genes is not cosmetic. It improves the dispersion estimates every subsequent test depends on and reduces the multiple-testing burden by half or more. Regulatory-grade omics frameworks specify the filter explicitly as part of the analysis plan rather than tuning it after seeing results, and we recommend the same discipline: decide your threshold before you look at p-values.23

library(DESeq2)
dds <- DESeqDataSetFromTximport(txi, colData = samples, design = ~ subject + timepoint)
keep <- rowSums(counts(dds) >= 10) >= (0.75 * min(table(samples$timepoint)))
dds <- dds[keep, ]

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

Before interpreting the PCA, run two identity checks. First, confirm sex-linked expression matches expectation: XIST high with RPS4Y1 and DDX3Y near zero, or the reverse. A sample that disagrees is mislabeled. Sex is also a genuine biological axis in expression data well beyond the sex chromosomes, and tools now exist specifically to separate sex-dependent molecular effects from the variable you care about, so include sex as a covariate in any multi-subject comparison.4 Second, if you have matched whole-genome data, call variants from the RNA alignments and check genotype concordance. This is the only reliable way to catch a sample swap.

Then look at what PC1 is. In whole blood it is very often neutrophil fraction, hemoglobin content, or RNA degradation, none of which is the signal you are after. If PC1 correlates with sequencing batch or draw date, note it and include the term in the design.

7. Test the design you have, not the one you wish you had

For a longitudinal personal profile, the standard two-group contrast is the wrong model. You have repeated measures on one individual, and the correct question is whether expression follows a coherent trajectory over time. Time-course RNA-seq has its own methodology, and the choice between a discrete-timepoint test and a continuous model changes what you can detect.5

For a small number of discrete timepoints, use the likelihood ratio test, which asks whether time explains variance at all rather than testing each pairwise contrast:

dds <- DESeq(dds, test = "LRT", reduced = ~ subject)
res <- results(dds, alpha = 0.05)
summary(res)

For densely sampled series (six or more timepoints), fit a natural spline so you model shape rather than arbitrary pairwise jumps:

library(splines)
design(dds) <- ~ subject + ns(time_days, df = 3)
dds <- DESeq(dds, test = "LRT", reduced = ~ subject)

Be direct with yourself about what a single-subject design supports. With one person, your replication is technical and temporal, not biological, and a gene that changes between two of your own draws could reflect an infection, a poor night of sleep, the time of day of the blood draw, or what you ate. Treat the output as a hypothesis list to be tested with more timepoints, not as a finding.

8. Interpret at the level of gene sets and cell types

Single-gene results from one individual are fragile. Aggregate statistics over gene sets are far more stable, because coordinated movement of forty genes in a pathway is much harder to produce by chance than movement of one.

library(fgsea); library(msigdbr)

ranks <- res$stat; names(ranks) <- rownames(res)
ranks <- sort(ranks[!is.na(ranks)])
h <- msigdbr(species = "Homo sapiens", collection = "H")
pathways <- split(h$gene_symbol, h$gs_name)

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

Rank by the Wald or LRT statistic, not by fold change alone, since fold change without a variance term puts noisy low-count genes at the extremes.

Then estimate cell composition, because in whole blood most apparent expression changes are composition changes. A three-fold rise in a neutrophil-specific transcript usually means more neutrophils, not more transcription per neutrophil. Deconvolution tools in the immunedeconv R package (CIBERSORTx, quanTIseq, EPIC) take a TPM matrix and return proportion estimates. Once you have them, re-run your differential test including the dominant fractions as covariates and see which results survive. The ones that do are the interesting ones.

If you have proteomic and metabolomic measurements from the same draws, correlating them with transcript levels is informative but should be done with low expectations for agreement at the single-gene level: transcript and protein abundance are coupled loosely, and each omic layer carries distinct noise and distinct blind spots.6 Integration adds power for pathway-level conclusions more than for individual analytes.7

Common problems

Globin dominance in whole blood is the most frequent surprise. If HBB alone takes 30 percent of your reads, your effective depth is a fraction of what you paid for, and low-expression transcripts become unmeasurable. Check quant.sf for HBA1, HBA2, HBB, and if they dominate, request globin-depleted or ribo-depleted library prep for subsequent draws. Removing globin reads computationally after the fact recovers some sensitivity but cannot recover reads that were never generated.

Degraded RNA produces 3-prime coverage bias that mimics differential expression across whole functional categories, since long transcripts lose proportionally more coverage than short ones. The tell is a correlation between log fold change and transcript length. Prevent it at the bench by checking RIN or DV200 before library prep; blood drawn into PAXgene tubes and frozen promptly typically holds RIN above 7.

A mapping rate under 60 percent usually has one of three causes: the wrong reference index, substantial rRNA carryover from an incomplete depletion, or untrimmed adapters. Salmon’s log reports these separately, and running a quick bowtie2 against an rRNA-only index distinguishes the second from the others.

Gene identifier mismatches waste an evening. GENCODE transcript and gene IDs carry version suffixes (ENSG00000141510.17), and most annotation packages do not. Strip them once at import with sub("\\..*", "", ids) and keep the annotation release fixed for the life of the project. Re-quantifying old samples against a new GENCODE release changes counts enough to produce spurious differences between old and new batches.

Batch effects from draw conditions are underappreciated in personal longitudinal data. Circadian time, fasting state, recent exercise, and even posture during the draw move blood transcriptomes measurably. Standardize the draw (same clinic, same hour, fasted) and record the covariates you cannot control, so you can include them in the model. Consistency in collection protocol buys more statistical power than doubling sequencing depth.

Finally, resist reading a single gene as a result. A pathway moving coherently across three consecutive timepoints is signal worth pursuing. One transcript at 2.5-fold with an adjusted p-value of 0.04 in a single individual is a coin flip. If something in the output looks medically important, the next step is a clinician and a validated clinical assay, not a deeper dive into the count matrix.

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. Daihan Ji, Mei Han, Shuting Lu, et al. SingleCellMQC: A comprehensive quality control workflow for single-cell multi-omics. iScience, 2026. https://doi.org/10.1016/j.isci.2026.117398 ↩

  2. M.C. Verheijen, T.W. Gant, W. Tong, et al. R-ODAF: Omics data analysis framework for regulatory application. Toxicology Letters, 2021. https://doi.org/10.1016/s0378-4274(21)00539-7 ↩

  3. Marcha Verheijen, Weida Tong, Leming Shi, et al. Towards the development of an omics data analysis framework. Regulatory Toxicology and Pharmacology, 2020. https://doi.org/10.1016/j.yrtph.2020.104621 ↩

  4. Sophie Le Bars, Mohamed Soudy, Enrico Glaab. XYomics: detecting sex-dependent molecular mechanisms in omics data. Nucleic Acids Research, 2026. https://doi.org/10.1093/nar/gkag759 ↩

  5. Daniel Spies, Constance Ciaudo. Dynamics in Transcriptomics: Advancements in RNA-seq Time Course and Downstream Analysis. Computational and Structural Biotechnology Journal, 2015. https://doi.org/10.1016/j.csbj.2015.08.004 ↩

  6. C. Nelson Hayes, Hikaru Nakahara, Atsushi Ono, et al. From Omics to Multi-Omics: A Review of Advantages and Tradeoffs. Genes, 2024. https://doi.org/10.3390/genes15121551 ↩

  7. Shikhi Baruri, Lalit Batra, Sohome Adhikari, et al. Integrative Epigenomics: Bioinformatics Strategies for Multi-Omics Data Analysis in Health and Disease. Epigenomes, 2026. https://doi.org/10.3390/epigenomes10030053 ↩