Working Through an RNA-Seq Dataset End to End

By the end of this guide you will have a gene-level counts matrix and a TPM matrix for a set of samples, a QC report that tells you whether the libraries are usable, a variance-stabilized expression matrix suitable for clustering and PCA, and a ranked table of genes with log fold changes and adjusted p-values. You need a machine with at least 16 GB of RAM (32 GB is more comfortable), roughly 100 GB of free disk for a small human dataset, a conda or mamba environment, and R 4.3 or later with Bioconductor. Everything below assumes bulk, poly-A selected, paired-end human RNA-seq, which is what you get from a standard blood or tissue transcriptome. Single-cell data is a different pipeline with different failure modes, and we flag where that matters.

1. Get a real example dataset

The fastest path to a usable example is the European Nucleotide Archive (ENA), which serves FASTQ files over plain HTTP and publishes MD5 checksums, so you avoid the SRA toolkit entirely. Pick a study accession and query the portal API for the file manifest:

curl -s "https://www.ebi.ac.uk/ena/portal/api/filereport?accession=PRJNA562734&result=read_run&fields=run_accession,sample_title,fastq_ftp,fastq_md5,read_count&format=tsv" > runs.tsv
head -3 runs.tsv

Each row gives you semicolon-separated URLs for the R1 and R2 files. Download two or three samples first rather than the whole study, because a single human RNA-seq run at 30 million read pairs is typically 2 to 4 GB compressed per direction.

For practice that resembles what you would do with your own data, we prefer a dataset where RNA-seq sits alongside other molecular layers rather than a standalone expression study. The frontotemporal dementia multi-omics resource from Menden and colleagues pairs bulk RNA-seq with methylation and proteomics on the same individuals, which lets you check whether a signal you find in transcript abundance shows up at the protein level.¹ Cross-layer agreement is the cheapest available control against a purely technical artifact.

If you want a synthetic dataset where you know the ground truth, generate one with polyester in R and simulate a two-fold change in 300 genes. Running your pipeline on simulated data once is the only way to find out whether your annotation, strandedness setting, and design matrix are wired together correctly.

2. Read the files before you run anything

A FASTQ file is four lines per read: an identifier beginning with @, the base calls, a + separator, and per-base quality scores encoded as ASCII characters offset by 33 (Phred+33). Look at the first record directly, because it tells you the read length, the instrument, and whether the file has already been trimmed:

zcat sample_R1.fastq.gz | head -4
zcat sample_R1.fastq.gz | awk 'NR%4==2 {print length($0)}' | head -1000 | sort -n | uniq -c

If read lengths are uniform at 100 or 150, the file is untrimmed. Ragged lengths mean someone already ran an adapter trimmer, and you should not run another one blindly. Count reads with echo $(( $(zcat sample_R1.fastq.gz | wc -l) / 4 )) and confirm that R1 and R2 have identical counts, in the same order. Mismatched pairs are the single most common reason a quantifier silently produces garbage.

Downstream you will meet two more formats. A salmon quant.sf file is a five-column TSV with transcript Name, Length, EffectiveLength, TPM, and NumReads, one row per transcript in the index. A counts matrix is genes by samples with integer values, and TPM (transcripts per million) is the same data normalized for transcript length and sequencing depth so that columns sum to one million. Counts are what statistical models consume. TPMs are what you look at when you want to ask how much of a given gene is present in one sample.

3. Run QC and read it properly

Run FastQC on every file and aggregate with MultiQC, which collapses dozens of reports into one HTML page:

mamba create -n rnaseq -c bioconda -c conda-forge fastqc multiqc salmon=1.10 fastp
conda activate rnaseq
fastqc -t 8 -o qc/ fastq/*.fastq.gz
multiqc -o qc/ qc/

Four panels carry most of the information. Per-base sequence quality should stay above Q28 for the bulk of the read, with the usual decline in the last 10 bases. Adapter content rising past 5 to 10 percent toward the 3’ end means insert sizes are shorter than the read length and you should trim. Sequence duplication levels above roughly 60 percent in a library with fewer than 20 million reads usually indicate low input RNA and PCR over-amplification. Per-sequence GC content that is bimodal rather than a single peak near 45 to 50 percent suggests contamination or a large rRNA fraction.

If trimming is needed, fastp is fast and writes its own report:

fastp -i fastq/S1_R1.fastq.gz -I fastq/S1_R2.fastq.gz \
      -o trim/S1_R1.fastq.gz -O trim/S1_R2.fastq.gz \
      --detect_adapter_for_pe --length_required 36 \
      --json qc/S1.fastp.json --html qc/S1.fastp.html --thread 8

One blood-specific check deserves attention. Whole blood RNA is dominated by hemoglobin transcripts, and without globin depletion HBB, HBA1, and HBA2 can consume more than half of all reads, which guts effective depth for everything else. Look at those three genes as a fraction of total TPM as soon as you have quantification, and treat a library with 60 percent globin as having roughly a third of its nominal read count.

4. Quantify with salmon against a decoy-aware index

We recommend salmon in selective-alignment mode over a full genome aligner for personal transcriptome work. It runs in 10 to 20 minutes per sample on 16 threads, uses under 16 GB of RAM, and produces transcript-level estimates with bootstrap replicates that propagate quantification uncertainty into downstream models. The tradeoff is that you get no BAM file, so you cannot inspect read pileups in a genome browser or call variants from the RNA. If you want that, run STAR as well, and be ready for a 32 GB memory footprint for a human index. Loading aligned reads into a browser next to the annotation is worth doing at least once, because seeing where reads fall across exons builds intuition that a counts table cannot.²

Build the index with genome decoys, which prevents reads from unannotated genomic regions being force-assigned to transcripts:

wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_45/gencode.v45.transcripts.fa.gz
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_45/GRCh38.primary_assembly.genome.fa.gz

grep "^>" <(gunzip -c GRCh38.primary_assembly.genome.fa.gz) | cut -d " " -f1 | 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 16 -i idx_gencode_v45 --gencode

The --gencode flag strips the pipe-delimited GENCODE headers down to the bare ENST identifier, which saves you a parsing step later. Index construction takes about 30 to 45 minutes and produces roughly 15 GB on disk.

Then quantify each sample:

for s in S1 S2 S3; do
  salmon quant -i idx_gencode_v45 -l A \
    -1 trim/${s}_R1.fastq.gz -2 trim/${s}_R2.fastq.gz \
    -p 16 --gcBias --seqBias --numBootstraps 30 \
    -o quants/${s}
done

-l A lets salmon infer library type from the data. Check what it inferred in quants/S1/lib_format_counts.json: ISR is the standard dUTP stranded protocol, IU is unstranded. If different samples in one batch infer different library types, something is wrong with the metadata. --gcBias and --seqBias correct fragment-level biases and cost little. --numBootstraps 30 gives you inferential replicates for tools that can use them.

The mapping rate in quants/S1/logs/salmon_quant.log is your first hard quality gate. For poly-A selected human libraries against a decoy-aware GENCODE index, expect 75 to 90 percent. Below 60 percent, stop and diagnose rather than proceeding, because a low rate almost always means the wrong species, residual adapters, heavy rRNA or genomic DNA contamination, or degraded RNA.

5. Import into R and build the expression object

Transcript-level estimates are summarized to genes with tximport, which also passes along average transcript lengths so that DESeq2 can correct for length differences between samples.

library(tximport); library(DESeq2); library(GenomicFeatures)

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

files <- file.path("quants", c("S1","S2","S3"), "quant.sf")
names(files) <- c("S1","S2","S3")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene,
                ignoreTxVersion = FALSE)

coldata <- data.frame(row.names = names(files),
                      condition = factor(c("ctrl","ctrl","treat")))
dds <- DESeqDataSetFromTximport(txi, coldata, ~ condition)

keep <- rowSums(counts(dds) >= 10) >= 2
dds  <- dds[keep, ]
vsd  <- vst(dds, blind = FALSE)
plotPCA(vsd, intgroup = "condition")

The version-suffix question causes more failures than anything else here. GENCODE transcript IDs carry a version (ENST00000456328.2), and if your tx2gene table and your quant.sf disagree on whether that suffix is present, tximport drops every transcript and reports zero overlap. Set ignoreTxVersion = TRUE if the two sources differ.

Before interpreting anything, look at the PCA. The first two principal components of the variance-stabilized matrix should separate samples by the biological factor you care about. If they separate by sequencing date, RNA extraction batch, or library prep kit instead, that variable is your dominant signal and must enter the design. Filtering to genes with at least 10 counts in at least as many samples as your smallest group typically takes a human dataset from about 62,000 annotated genes to 14,000 to 18,000 measurable ones.

6. Test for differential expression

With replicates in each group, run the standard negative-binomial test and shrink the fold changes so that lowly expressed genes with wild ratios do not top your list:

dds <- DESeq(dds)
res <- lfcShrink(dds, coef = "condition_treat_vs_ctrl", type = "apeglm")
res <- res[order(res$padj), ]
sum(res$padj < 0.05 & abs(res$log2FoldChange) > 0.5, na.rm = TRUE)

Three replicates per group at 25 to 30 million read pairs is the usual floor for detecting two-fold changes in moderately expressed genes. Two per group will find only very large effects. One per group supports no formal test at all.

The single-person case needs a different framing. If you have your own transcriptome at several timepoints and no control group, you cannot do a between-group comparison, and forcing one produces confident-looking nonsense. What you can do is treat each gene’s trajectory across your own timepoints as the unit of analysis: variance-stabilize the full matrix, compute each gene’s z-score against your own prior timepoints, and require a consistent direction across at least three consecutive samples before you take a change seriously. Interpretation is more reliable at the pathway level than the single-gene level, so run fgsea on Hallmark or Reactome sets ranked by within-person change rather than staring at individual genes. Placing your own values inside a larger public cohort is also legitimate and is roughly what integrative meta-analyses do when they pool many transcriptomic datasets to find signals no single study resolves.³

Note the boundary plainly: none of this is a clinical result. Expression changes have no diagnostic meaning outside a validated assay, and anything that looks concerning belongs with a physician who can order a test that was designed to answer the question.

7. Sanity-check against biology you already know

Every pipeline needs a positive control that does not depend on the pipeline being right. Three cheap ones for human data:

tpm <- txi$abundance
tpm[c("ENSG00000229807.12",  # XIST
      "ENSG00000129824.16",  # RPS4Y1
      "ENSG00000244734.4"),] # HBB

XIST should be high and Y-linked genes near zero in samples from female donors, and the reverse in male donors. Disagreement with recorded metadata means a sample swap, and sample swaps are common. In whole blood, HBB should be present at high TPM and the neutrophil and lymphocyte marker balance should track a contemporaneous complete blood count if you have one. When the same individuals have proteomic or methylation data, checking whether a transcript-level change moves in the same direction in an orthogonal layer is the strongest confirmation available short of a replication cohort.¹⁴

Common problems

Low mapping rate with clean FastQC output usually means an index and library mismatch: a total-RNA library quantified against a poly-A-oriented expectation, or reads from a different genome build. Check lib_format_counts.json and the top unmapped sequences by BLASTing a handful of reads.

Zero genes after tximport is nearly always the transcript ID version mismatch described above. Compare head(rownames(txi$counts)) against head(tx2gene$TXNAME) before assuming anything more exotic.

A PCA where one sample sits far from all others is typically a degradation artifact. Low-RIN RNA produces 3’ coverage bias, which inflates short transcripts and depletes long ones. You can detect it by plotting median TPM against annotated transcript length: a strong negative slope in one sample and not others is diagnostic. That sample should be dropped rather than modeled around.

Batch confounded with condition cannot be fixed computationally. If all treated samples were prepped on Tuesday and all controls on Thursday, the design is unidentifiable and no amount of removeBatchEffect recovers the truth.

Applying this workflow to single-cell data will not work. Droplet data needs UMI deduplication, empty-droplet filtering, and per-cell normalization, followed by clustering and cell-type annotation, which is a research problem of its own with dozens of competing methods.⁵ Pseudotemporal ordering of cells along a developmental trajectory is a further layer again, built on the assumption that a snapshot population contains cells at many stages of one process.⁶ Use a tool built for it rather than adapting a bulk pipeline.

Finally, memory. makeTxDbFromGFF on a full GENCODE GTF takes several minutes and a few GB. If you are running on a laptop, build the tx2gene table once, save it with saveRDS, and never parse the GTF again.

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. Kevin Menden, Margherita Francescatto, Tenzin Nyima, et al. A multi-omics dataset for the analysis of frontotemporal dementia genetic subtypes. Scientific Data, 2023. https://doi.org/10.1038/s41597-023-02598-x ↩ ↩²

2. Ann E. Loraine, Ivory Clabaugh Blakley, Sridharan Jagadeesan, et al. Analysis and Visualization of RNA-Seq Expression Data Using RStudio, Bioconductor, and Integrated Genome Browser. Methods in Molecular Biology, 2015. https://doi.org/10.1007/978-1-4939-2444-8_24 ↩

3. Maxim N. Shokhirev, Adiv A. Johnson. An integrative machine-learning meta-analysis of high-throughput omics data identifies age-specific hallmarks of Alzheimer’s disease. Ageing Research Reviews, 2022. https://doi.org/10.1016/j.arr.2022.101721 ↩

4. Upasna Srivastava, Swarna Kanchan, Minu Kesheri, et al. Integrative omics approaches for identification of biomarkers. Integrative Omics, 2024. https://doi.org/10.1016/b978-0-443-16092-9.00010-2 ↩

5. Changde Cheng, Wenan Chen, Hongjian Jin, et al. A Review of Single-Cell RNA-Seq Annotation, Integration, and Cell–Cell Communication. Cells, 2023. https://doi.org/10.3390/cells12151970 ↩

6. Jaehoon Shin, Daniel A. Berg, Yunhua Zhu, et al. Single-Cell RNA-Seq with Waterfall Reveals Molecular Cascades underlying Adult Neurogenesis. Cell Stem Cell, 2015. https://doi.org/10.1016/j.stem.2015.07.013 ↩