Best Practices for RNA-seq Data Analysis, Applied to Your Own Sample

The short answer is that a defensible RNA-seq pipeline in 2026 has four decisions in it, and the rest is plumbing. You choose a library prep that matches the RNA you care about, and you sequence deep enough for the genes you care about. You then quantify against a transcriptome with a selective aligner such as Salmon or a splice-aware aligner such as STAR. Finally, you normalize before you compare anything to anything. The canonical review of these steps is still Conesa and colleagues, and its structure remains the right mental model. That structure runs from design through quality control, alignment and quantification to differential expression, then everything downstream.¹ What that review does not cover is the case most readers here have, which is a small number of samples from one person over time. A design like that changes the statistics more than it changes the tooling.

Library prep decisions constrain everything downstream

RNA-seq is not one assay. Poly(A) selection captures polyadenylated mRNA and gives you the cleanest gene-level quantification per dollar. Ribosomal depletion (Ribo-Zero style) keeps non-polyadenylated transcripts, including many long non-coding RNAs and unspliced precursors, at the cost of spending reads on intronic and residual rRNA sequence. If you are profiling blood and want a stable, comparable gene-level baseline, poly(A) is the better default.

Whole blood adds one more problem. Globin transcripts from reticulocytes can consume a large share of your reads, which effectively lowers the depth available for everything else. Either use a globin depletion step or collect PBMCs instead. If you do neither, expect to discard the sample when you find HBB and HBA2 dominating the count matrix.

Strandedness matters and is frequently misrecorded. Most current kits are reverse-stranded (dUTP protocols), which is -s 2 in featureCounts and ISR in Salmon’s notation. Do not trust the vendor sheet. Run salmon quant -l A on a subset of reads and read back the inferred library type from lib_format_counts.json. Getting this wrong silently halves your counts for overlapping genes.

Depth, replication, and read length

The depth you need depends entirely on what you want to detect. A few million mapped reads per sample is enough to quantify moderately and highly expressed genes, while lowly expressed transcripts and isoform-level work need considerably more.¹ For a personal gene-level baseline from blood, we target 30 to 50 million read pairs. For isoform-level questions, alternative splicing, or allele-specific expression, plan on 100 million or more.

Read length and pairing follow the same logic. Single-end 50 bp reads are sufficient for gene-level counting, and the review is explicit that cheaper short single-end reads are normally enough for that purpose.¹ Paired-end reads of 100 to 150 bp are what you want for isoform assignment, fusion detection, and any variant calling from RNA. Paired-end also gives you fragment length distributions, which help the effective length correction inside Salmon and kallisto.

Replication is where individual profiling gets uncomfortable. The standard recommendation is at least three biological replicates per condition, because fewer leaves you unable to estimate within-condition variance.¹ With one person and one draw you have zero replicates and no variance estimate. The fix is time: repeated draws under matched conditions become your replicate structure.

Quality control before you count anything

Run FastQC on the raw FASTQ files, then aggregate with MultiQC so you look at one report instead of forty. The metrics that change decisions are per-base quality decay, adapter content, duplication rate, and the GC distribution. A bimodal GC curve usually means contamination or rRNA carryover.

For trimming, we use fastp rather than Trimmomatic. The main reason is that it does adapter detection, quality filtering, and reporting in one pass:

fastp -i R1.fq.gz -I R2.fq.gz -o R1.trim.fq.gz -O R2.trim.fq.gz \
  --detect_adapter_for_pe --trim_poly_g --cut_tail --cut_tail_mean_quality 20 \
  --length_required 36 --json fastp.json --html fastp.html

Be conservative. Aggressive quality trimming removes real signal and biases quantification more than the low-quality bases it removes. Do not deduplicate RNA-seq reads without unique molecular identifiers. Highly expressed transcripts generate genuine duplicate fragments, and removing them distorts the count distribution.

Post-alignment QC is where you catch the failures that FastQC cannot see. The numbers to record for every sample are the uniquely mapped percentage from STAR’s Log.final.out, the fraction of reads in exons, the rRNA fraction, and the 5’ to 3’ coverage uniformity. Degraded RNA shows up as a strong 3’ bias. Check the RIN from the lab if you have it, and treat anything below about 7 as suspect for isoform work.

Selective quantification versus genome alignment

For gene and transcript expression, we default to Salmon against a GENCODE transcriptome with a decoy-aware index built from the genome. It is fast, it models fragment GC and positional bias, and it gives you transcript-level estimates with uncertainty:

salmon index -t gentrome.fa.gz -d decoys.txt -i salmon_index --gencode -k 31
salmon quant -i salmon_index -l A -1 R1.trim.fq.gz -2 R2.trim.fq.gz \
  --validateMappings --seqBias --gcBias --numGibbsSamples 20 \
  -p 16 -o quant/sample01

Then import to gene level with tximport, using the countsFromAbundance = "lengthScaledTPM" option so that downstream count-based models see length-corrected values.

Align to the genome with STAR when you need the BAM for something else: novel splice junctions, fusion calling with Arriba or STAR-Fusion, variant calling, or visual inspection in IGV. Use --quantMode TranscriptomeSAM GeneCounts and a genome index built with --sjdbOverhang set to read length minus one. The tradeoff is real. STAR needs roughly 32 GB of RAM for a human index and produces BAM files that are ten to twenty times larger than a quant.sf file. Convert finished BAMs to CRAM with samtools view -C -T ref.fa if you are keeping them.

Pin your annotation version and never change it mid-project. Gene identifiers move between GENCODE releases, and a silent annotation swap will make two of your own timepoints look differentially expressed.

Normalization and what TPM cannot do

TPM and FPKM correct for library size and transcript length, which makes them reasonable for comparing genes within one sample. They are poor for comparing one gene across samples, because a few dominant transcripts shift the denominator for everything else. This composition problem is exactly what the count-based normalizations are built to handle.

Use raw or length-scaled counts with DESeq2’s median-of-ratios normalization, or edgeR’s TMM. Both assume most genes are unchanged between samples, which is usually true for repeated blood draws from one person. If you suspect global transcriptional shifts, you need spike-in controls to normalize against, and without them the assumption fails quietly.

Differential expression itself is a negative binomial problem, and DESeq2 or edgeR handles it with shrinkage of dispersion estimates across genes.¹ limma-voom is the better choice once you have more than roughly a dozen samples, because it is faster and handles complex designs cleanly.

Analyzing repeated samples from one person

With one individual, the comparison is you against yourself, and your enemy is batch effect. A single sequencing run is a single batch. Any two samples processed months apart differ by library prep, flow cell, and time of day of the draw. Standardize what you can. Use the same tube type, the same time of morning, the same fasting state, and the same lab.

For analysis, include the batch as a covariate in the design formula rather than correcting counts in place. Something like ~ batch + timepoint in DESeq2 preserves the correct degrees of freedom. Use limma::removeBatchEffect only for visualization, never as input to a differential test.

Expect a sizeable fraction of variance to come from cell composition rather than regulation. Whole blood expression tracks neutrophil and lymphocyte proportions strongly, so a shift in your neutrophil fraction will look like hundreds of differentially expressed genes. Deconvolve with CIBERSORTx or a reference signature matrix, then include the estimated fractions as covariates. Pair the RNA with a CBC from the same draw when you can.

Finally, set an interpretation floor. Gene-level fold changes below about 1.3 in a single pair of samples are not distinguishable from technical variation without replicates. Look at pathway-level and module-level aggregates instead, which average over gene-level noise. Any finding that touches your health belongs in front of a clinician who can order a confirmatory clinical-grade test. Research RNA-seq is not a diagnostic.

Questions people also ask

Should I use Salmon or STAR? Use Salmon for expression quantification, because it is faster, needs less memory, and models fragment-level bias. Use STAR when you need genomic coordinates for junctions, fusions, variants, or browser inspection. Running both on the same sample is cheap and worth doing once.

How many reads do I need? For gene-level quantification of moderately to highly expressed genes, a few million mapped reads is enough, while lowly expressed genes and isoform-level analysis require substantially deeper sequencing.¹ We target 30 to 50 million read pairs for a general baseline.

Do I need biological replicates if I only have my own data? Yes, and repeated draws over time are your replicates. Three or more samples per condition is the usual minimum for estimating variance.¹ A single timepoint supports description. It does not support inference.

Is TPM fine for comparing across samples? Not for per-gene comparisons. Use DESeq2 median-of-ratios or edgeR TMM normalized counts for cross-sample work, and reserve TPM for within-sample comparisons between genes.

Should I deduplicate reads? Not unless your library has UMIs. Duplicate fragments from highly expressed transcripts are real signal, and removing them biases quantification.

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 ↩ ↩² ↩³ ↩⁴ ↩⁵ ↩⁶ ↩⁷