How to Run an RNA-seq Analysis on Your Own Data, Step by Step

By the end of this guide you will have taken a set of paired-end FASTQ files from a bulk RNA-seq run and produced three things: a per-sample transcript quantification table, a gene-level count matrix, and a ranked list of differentially expressed genes with shrunken effect sizes and adjusted p-values. Along the way you will have the diagnostic artifacts that tell you whether to trust any of it: per-base quality profiles, ribosomal and mitochondrial read fractions, mapping rates, a sample-distance heatmap, and a principal component plot. You need a Linux or macOS machine with at least 16 GB of RAM (32 GB is more comfortable), roughly 100 GB of free disk per dozen samples, conda or a container runtime, and R 4.3 or newer. The commands below assume human data and the GENCODE annotation, but the shape of the workflow is species-agnostic. If you have a single time point from a single person, the differential expression section still applies conceptually but you will be comparing across time or across tissue rather than across a case/control design, and the statistical caveats at the end matter a great deal.

Before any of this, understand what the wet-lab side did, because it constrains everything downstream. RNA was extracted, assessed for integrity (RIN, the RNA integrity number, on a Bioanalyzer or TapeStation), then either poly-A selected or ribo-depleted, fragmented, reverse transcribed to cDNA, adapter-ligated, PCR-amplified, and sequenced. Poly-A selection captures mature mRNA and gives you the cleanest gene expression signal per read. Ribo-depletion keeps non-polyadenylated species including most long non-coding RNAs and pre-mRNA, at the cost of a higher intronic fraction. Ask which one was used, ask whether the library is stranded (almost all modern kits are dUTP-based reverse-stranded), and ask for the RIN values. A RIN below 7 means degraded RNA, which biases coverage toward 3’ ends and will show up later as a coverage skew you cannot fix computationally.

1. Set up the environment and get the annotation right

Reproducibility here comes from pinning versions, not from remembering what you typed. Create an isolated environment:

conda create -n rnaseq -c conda-forge -c bioconda \
  fastqc=0.12.1 multiqc=1.21 fastp=0.23.4 \
  salmon=1.10.3 samtools=1.19 star=2.7.11b
conda activate rnaseq

Annotation choice matters more than most people expect, because gene identifiers, transcript sets, and which loci exist at all differ between GENCODE, Ensembl, and RefSeq. We use GENCODE for human, because the transcript FASTA and GTF are versioned together and the identifiers carry version suffixes that make provenance explicit. Pick one release and use it for every sample in a study.

REL=45
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/gencode.v${REL}.annotation.gtf.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_${REL}/GRCh38.primary_assembly.genome.fa.gz

Record the release number in your analysis notes. Six months from now, when a gene symbol has been renamed and your figure no longer matches a collaborator’s, that number is how you resolve the discrepancy.

2. Quality control the raw reads

Run FastQC on every FASTQ and aggregate with MultiQC. This takes minutes and catches problems that would otherwise surface as inexplicable results three steps later.

fastqc -t 8 -o qc/fastqc fastq/*.fastq.gz
multiqc -o qc qc/fastqc

Read the report in a specific order. First, per-base sequence quality: Illumina reads normally hold a Phred score above 30 for most of the read with a decline in the last 10 to 15 bases. Second, adapter content: a rising curve at the 3’ end means short inserts, which is common when the library was over-fragmented. Third, overrepresented sequences: if the top hit is a poly-G run you are looking at a NovaSeq two-color chemistry artifact where no signal is read as G, and if it is a ribosomal RNA sequence your depletion or selection underperformed. Fourth, per-sequence GC content: a sharp secondary peak usually means contamination or a dominant rRNA population.

The single most useful number at this stage is the duplicate rate, but interpret it carefully. High duplication in RNA-seq is often biological, because highly expressed transcripts genuinely produce identical fragments. Do not deduplicate bulk RNA-seq reads unless your library has unique molecular identifiers.

3. Trim adapters and low-quality tails

We use fastp rather than Trimmomatic because it is faster, detects adapters automatically for paired-end data, handles the poly-G problem directly, and emits a JSON report that MultiQC ingests. Trimmomatic remains fine and is well documented in the older literature, where trimming, quality assessment, and alignment were commonly bundled into a single interface for users who did not want to script them ¹.

for s in $(cat samples.txt); do
  fastp \
    -i fastq/${s}_R1.fastq.gz -I fastq/${s}_R2.fastq.gz \
    -o trim/${s}_R1.fq.gz     -O trim/${s}_R2.fq.gz \
    --detect_adapter_for_pe \
    --trim_poly_g \
    --qualified_quality_phred 15 \
    --unqualified_percent_limit 40 \
    --length_required 36 \
    --thread 8 \
    --json qc/fastp/${s}.json --html qc/fastp/${s}.html
done

Resist aggressive quality trimming. Cutting every base below Q30 shortens reads, reduces mappability, and introduces a length bias that varies by sample. The defaults above remove adapters and genuinely bad tails and leave the rest alone. Aim to retain more than 90 percent of reads after trimming. If a sample drops below 80 percent, open its fastp HTML report and find out why before proceeding.

4. Quantify with salmon

Here is the main decision in the workflow. You can align reads to the genome with a splice-aware aligner (STAR or HISAT2) and then count reads per gene (featureCounts or htseq-count), or you can quantify directly against the transcriptome with a selective-alignment method (salmon or kallisto). Most published pipelines follow the align-then-count path, and pipeline frameworks built around STAR plus featureCounts are well established and reproducible across compute environments ². Automated pipelines that wrap trimming, alignment, quantification, and differential expression into one invocation exist and work well if you want the whole flowchart in a single command ³.

We recommend salmon for standard gene expression work. It is roughly an order of magnitude faster, needs far less memory, handles multi-mapping reads across transcript isoforms probabilistically instead of discarding them, and corrects for GC and positional bias in the fragment model. The tradeoff is that you get no BAM file, so you cannot inspect coverage over a locus, call variants from the RNA, or detect novel splice junctions. If you need any of those, run STAR as well.

Build the index with a decoy set of genomic sequence, which prevents reads from intronic and intergenic regions from being forced onto transcripts:

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_idx --gencode

The index build takes about 30 minutes and 20 GB of RAM for human. Then quantify:

for s in $(cat samples.txt); do
  salmon quant -i salmon_idx -l A \
    -1 trim/${s}_R1.fq.gz -2 trim/${s}_R2.fq.gz \
    --validateMappings --gcBias --seqBias \
    --numBootstraps 30 \
    -p 12 -o quant/${s}
done

-l A lets salmon infer library type; check quant/${s}/lib_format_counts.json afterward and confirm it reports ISR for a standard reverse-stranded paired-end kit. --gcBias and --seqBias cost a little runtime and meaningfully reduce technical variation between samples prepared in different batches. The bootstraps give you inferential replicates, which matter if you later do transcript-level analysis. Each sample takes 5 to 15 minutes on 12 threads for 30 million read pairs.

Check the mapping rate in quant/${s}/logs/salmon_quant.log. For poly-A selected human libraries, expect 75 to 90 percent. Below 60 percent, suspect contamination, rRNA carryover, or the wrong reference species.

5. Import to gene level and build the count matrix

Transcript-level estimates need to be summarized to genes, with an offset that accounts for the fact that average transcript length differs between samples. tximport handles both.

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

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

samples <- read_tsv("samples.tsv")   # columns: sample, condition, batch
files   <- file.path("quant", samples$sample, "quant.sf")
names(files) <- samples$sample
stopifnot(all(file.exists(files)))

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

If tximport warns about missing transcripts, your GTF and transcript FASTA are from different releases. Fix that rather than setting ignoreTxVersion = TRUE to make the warning go away, because mismatched annotations silently drop genes.

6. Exploratory analysis before any statistics

Look at the data before you test anything. This step catches swapped sample labels, failed libraries, and batch effects that would otherwise be attributed to biology.

library(DESeq2); library(pheatmap)

dds <- DESeqDataSetFromTximport(txi, colData = samples,
                                design = ~ batch + condition)
keep <- rowSums(counts(dds) >= 10) >= min(table(samples$condition))
dds  <- dds[keep, ]

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

sampleDists <- dist(t(assay(vsd)))
pheatmap(as.matrix(sampleDists), clustering_distance_rows = sampleDists)

The variance stabilizing transform is the right choice over a simple log because it removes the mean-variance dependence that otherwise makes low-count genes dominate distance calculations. On the PCA, ask whether PC1 separates your conditions or your batches. If batch dominates and batch is confounded with condition, no software will rescue the experiment.

One check that is easy to skip and often informative: confirm genetic sex from expression of XIST and the Y-linked genes RPS4Y1, DDX3Y, and UTY. A mismatch against the recorded metadata means a sample swap. Sex differences in expression are pervasive and extend well beyond the sex chromosomes, and tools now exist specifically to test for sex-dependent effects across omics data types ⁴.

7. Differential expression

With the design set and the exploratory plots clean, the test itself is short.

dds <- DESeq(dds)
resultsNames(dds)

res <- results(dds, contrast = c("condition", "treated", "control"),
               alpha = 0.05)
res_shrunk <- lfcShrink(dds, coef = "condition_treated_vs_control",
                        type = "apeglm")

summary(res)
write.csv(as.data.frame(res_shrunk[order(res_shrunk$padj), ]),
          "de_results.csv")

Two points about interpretation. First, always report the shrunken log fold changes, because raw fold changes for low-count genes are wildly inflated and shrinkage pulls them toward zero in proportion to their uncertainty. Second, filter by adjusted p-value first and fold change second, never the reverse. Benjamini-Hochberg control at 0.05 means roughly 5 percent of your called genes are expected to be false positives, which is a property of the list as a whole.

DESeq2 is our default. edgeR with quasi-likelihood F-tests performs comparably, and limma-voom is faster and better behaved with large sample counts (above roughly 30 per group). The three disagree most at small sample sizes and on genes with extreme dispersion. Comparative reviews of differential expression pipelines find that method choice, normalization, and filtering thresholds all shift the resulting gene lists, which argues for fixing your choices before you look at results rather than after ⁵.

8. Functional interpretation

A ranked gene list is an intermediate, not a result. Move to pathway level with gene set enrichment analysis, which uses the full ranked list rather than an arbitrary significance cutoff.

library(clusterProfiler); library(org.Hs.eg.db)

ranks <- res_shrunk$log2FoldChange
names(ranks) <- sub("\\..*", "", rownames(res_shrunk))
ranks <- sort(ranks[!is.na(ranks)], decreasing = TRUE)

gsea <- gseGO(ranks, OrgDb = org.Hs.eg.db, keyType = "ENSEMBL",
              ont = "BP", minGSSize = 15, maxGSSize = 500,
              pvalueCutoff = 0.05, eps = 0)

The standard toolchain for this (over-representation tests, gene set enrichment, and pathway topology methods) has been stable for over a decade and is well described in the methods literature ⁶. Report the normalized enrichment score, the adjusted p-value, and the size of the gene set. A pathway with eight genes and a spectacular score is usually noise.

Expression alone rarely settles a question. Studies that combine transcriptomic signal with proteomic, methylation, or clinical layers consistently resolve mechanisms that single-layer analysis leaves ambiguous, because transcript abundance and protein abundance correlate imperfectly ⁷. Integration across layers is now the standard approach in the personalized-profiling literature ⁸, and cross-condition integrative designs can surface shared regulatory signatures that no single dataset would reveal ⁹.

Common problems

Low mapping rate with high-quality reads almost always means an annotation or contamination problem. Take 100,000 unmapped reads and BLAST a handful. Bacterial rRNA, adapter dimers, and the wrong species all show up immediately.

A high intronic fraction (above roughly 30 percent for a poly-A library) points to genomic DNA carryover or nuclear pre-mRNA. You can estimate this only if you produced a BAM, which is one argument for running STAR alongside salmon on at least a few samples.

Degraded RNA produces a 3’ coverage bias that mimics differential expression of long transcripts. If RIN values vary across your groups, include RIN as a covariate in the design and check whether transcript length correlates with your fold changes.

Batch effects that are partially confounded with condition can sometimes be modeled with surrogate variable analysis (the sva package) or RUVSeq, which estimate unwanted variation from control genes. These help when the confounding is partial. They cannot separate two variables that move together perfectly.

If you are working from a single person’s longitudinal samples rather than a group comparison, the multiple-testing framework above still holds, but “replicate” means repeated sampling of one individual, and the inference is about that person’s variation over time rather than about a population. The same principle applies at single-cell resolution, where cells from one individual are not independent replicates of a biological condition ¹⁰.

Finally, a note on scope. Nothing in this workflow constitutes a clinical finding. Expression changes are measurements, and mapping them to health status requires clinical context, orthogonal confirmation, and a physician who can see the rest of your record. If a result concerns you, take it to one.

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. M. Lohse, A. M. Bolger, A. Nagel, et al. RobiNA: a user-friendly, integrated software solution for RNA-Seq-based transcriptomics. Nucleic Acids Research, 2012. https://doi.org/10.1093/nar/gks540 ↩

2. Phelelani Mpangase, Jacqueline Frost, Mohammed Tikly, et al. nf-rnaSeqCount: A Nextflow pipeline for obtaining raw read counts from RNA-seq data. South African Computer Journal, 2021. https://doi.org/10.18489/sacj.v33i2.830 ↩

3. Giulio Spinozzi, Valentina Tini, Alessia Adorni, et al. ARPIR: automatic RNA-Seq pipelines with interactive report. BMC Bioinformatics, 2020. https://doi.org/10.1186/s12859-020-03846-2 ↩

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. Diletta Rosati, Maria Palmieri, Giulia Brunelli, et al. Differential gene expression analysis pipelines and bioinformatic tools for the identification of specific biomarkers: A review. Computational and Structural Biotechnology Journal, 2024. https://doi.org/10.1016/j.csbj.2024.02.018 ↩

6. Moritz Kebschull, Melanie Julia Fittler, Ryan T. Demmer, et al. Differential Expression and Functional Analysis of High-Throughput -Omics Data Using Open Source Tools. Methods in Molecular Biology, 2016. https://doi.org/10.1007/978-1-4939-6685-1_19 ↩

7. Xinglin Li, Yiqi Xiong, Jiyin Wang, et al. Integrative multi-omics identifies PANX2 as a ferroptosis suppressor and potential therapeutic target in clear cell renal cell carcinoma. Frontiers in Immunology, 2026. https://doi.org/10.3389/fimmu.2026.1883667 ↩

8. Somayah Albaradei. AI-Driven Approaches to Utilization of Multi-Omics Data for Personalized Diagnosis and Treatment of Cancer: A Comprehensive Review. Computer Modeling in Engineering & Sciences, 2025. https://doi.org/10.32604/cmes.2025.072584 ↩

9. Danni Huang, Junying Wu, Jinhua Chen, et al. Integrated multi-omics analysis identifies candidate eRNA-associated signatures shared between osteoarthritis and type 2 diabetes. Frontiers in Genetics, 2026. https://doi.org/10.3389/fgene.2026.1875546 ↩

10. Xun Zhu, Lana X. Garmire. Data Analysis in Single-Cell RNA-Seq. Single-Cell Omics, 2019. https://doi.org/10.1016/b978-0-12-814919-5.00019-1 ↩