RNA-Seq Data Analysis: Complete Step-by-Step Guide for Beginners

RNA sequencing has become the default readout for transcriptome-scale biology. A modern bulk RNA-Seq experiment can go from wet-bench extraction to a ranked list of differentially expressed genes in under a day of compute — if your pipeline is set up correctly. This guide walks through a reproducible, production-grade RNA-Seq workflow for a human study, using tools that are widely used in publications and industry.

We assume you have a Linux workstation or HPC allocation with at least 32 GB RAM and access to conda or mamba. All commands are runnable; substitute your paths where needed.

The end-to-end pipeline at a glance

A canonical short-read bulk RNA-Seq analysis has six computational stages:

Stage	Purpose	Typical tool
1. Quality control	Detect adapter contamination, base-quality drops, duplication	FastQC, MultiQC
2. Trimming	Remove adapters and low-quality tails	fastp, Trim Galore
3. Alignment	Map reads to a reference genome	STAR, HISAT2
4. Post-alignment QC	Insert-size, strand specificity, gene-body coverage	RSeQC, Qualimap
5. Quantification	Count reads per gene / transcript	featureCounts, Salmon
6. Differential expression	Test for changes across conditions	DESeq2, edgeR, limma-voom

You don’t need to memorize this — bookmark it. Every subsection below maps to one row.

1. Setting up a reproducible environment

Reproducibility begins with your environment. Use conda (or the much faster mamba) to pin exact versions:

mamba create -n rnaseq -c bioconda -c conda-forge \
  fastqc=0.12.1 \
  fastp=0.23.4 \
  star=2.7.11b \
  samtools=1.20 \
  subread=2.0.6 \
  rseqc=5.0.3 \
  multiqc=1.25 \
  r-base=4.4.1 \
  bioconductor-deseq2=1.44.0
mamba activate rnaseq

Freeze the environment before sharing:

mamba env export --no-builds > rnaseq.yml

Tip. In a real study, drop the whole pipeline inside a Snakemake or Nextflow workflow and record the git hash of your config in every report. Ad-hoc bash pipelines rot within months.

2. Quality control — read your FastQC report before you run anything else

Before alignment, run FastQC on every FASTQ file and aggregate with MultiQC:

mkdir -p qc/raw
fastqc -o qc/raw fastq/*.fastq.gz
multiqc qc/raw -o qc/raw/multiqc

What you’re looking for:

• Per-base sequence quality: Q30+ across most of the read length.
• Adapter content: should be near-zero after trimming.
• Duplication levels: RNA-Seq is expected to have some duplication (highly expressed transcripts), but library-prep artifacts spike this.
• Overrepresented sequences: look for rRNA (28S, 18S) — a rRNA-depletion problem, or globin — a blood sample problem.

3. Trimming with fastp

fastp is fast (multi-threaded), applies quality trimming and adapter detection automatically, and produces a JSON report you can pipe into MultiQC:

for r1 in fastq/*_R1.fastq.gz; do
  sample=$(basename "$r1" _R1.fastq.gz)
  r2=${r1/_R1/_R2}
  fastp \
    -i "$r1" -I "$r2" \
    -o trimmed/"${sample}"_R1.fastq.gz \
    -O trimmed/"${sample}"_R2.fastq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 30 \
    --thread 8 \
    --json qc/fastp/"${sample}".json \
    --html qc/fastp/"${sample}".html
done

4. Alignment with STAR

STAR is the workhorse spliced aligner for RNA-Seq. Build the index once:

STAR --runMode genomeGenerate \
  --genomeDir ref/star_index \
  --genomeFastaFiles ref/GRCh38.primary_assembly.genome.fa \
  --sjdbGTFfile ref/gencode.v46.annotation.gtf \
  --sjdbOverhang 100 \
  --runThreadN 16

Then align each sample:

STAR --runMode alignReads \
  --genomeDir ref/star_index \
  --sjdbGTFfile ref/gencode.v46.annotation.gtf \
  --readFilesIn trimmed/"${sample}"_R1.fastq.gz trimmed/"${sample}"_R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMstrandField intronMotif \
  --outFileNamePrefix bam/"${sample}"_ \
  --twopassMode Basic \
  --runThreadN 16
samtools index "bam/${sample}_Aligned.sortedByCoord.out.bam"

Two-pass mode discovers novel splice junctions during pass one and uses them in pass two — a real accuracy improvement for less-annotated species and important for splice-QTL work.

5. Post-alignment QC

Once you have BAMs, run three sanity checks:

# Strand specificity
infer_experiment.py -i bam/sample_Aligned.sortedByCoord.out.bam \
  -r ref/gencode.v46.bed

# Gene body coverage (5'/3' bias)
geneBody_coverage.py -i bam/*.bam -r ref/hg38.HouseKeepingGenes.bed \
  -o qc/geneBody

# Overall stats
samtools flagstat bam/sample_Aligned.sortedByCoord.out.bam

A well-prepared library shows even 5’-to-3’ coverage. A strong 3’ bias means degradation — flag those samples in your DE model.

6. Quantification with featureCounts

featureCounts (from the subread package) is fast and battle-tested for gene-level counts:

featureCounts \
  -a ref/gencode.v46.annotation.gtf \
  -o counts/gene_counts.tsv \
  -p --countReadPairs \
  -s 2 \
  -T 16 \
  bam/*_Aligned.sortedByCoord.out.bam

Key flags:

• -p --countReadPairs — count fragments, not reads, for paired-end data.
• -s 2 — reverse-stranded protocol (typical for Illumina TruSeq stranded); confirm with the RSeQC result above.

7. Differential expression with DESeq2

Load counts into R and build the DE model. A minimal script:

library(DESeq2)
library(tidyverse)

counts <- read.table("counts/gene_counts.tsv", header = TRUE,
                     sep = "\t", skip = 1, row.names = 1)
counts <- counts[, 6:ncol(counts)]  # drop the annotation columns
colnames(counts) <- gsub("_Aligned.sortedByCoord.out.bam", "",
                         basename(colnames(counts)))

coldata <- data.frame(
  sample    = colnames(counts),
  condition = factor(c("ctrl","ctrl","ctrl","trt","trt","trt")),
  row.names = colnames(counts)
)

dds <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
dds <- dds[rowSums(counts(dds)) >= 10, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "trt", "ctrl"),
               alpha = 0.05, lfcThreshold = 0)
summary(res)

For a full DESeq2 walkthrough (LFC shrinkage, MA plots, batch correction), see our dedicated guide: How to Use DESeq2 for Differential Expression Analysis.

8. Reporting — MultiQC ties it all together

Finally, run MultiQC across every step’s output. It will pick up FastQC, fastp, STAR log files, RSeQC, and featureCounts summaries automatically:

multiqc . -o reports/multiqc_final

The resulting HTML is the single source of truth you attach to a supplementary methods file.

Common pitfalls & how to avoid them

• Wrong strand setting silently halves your counts. Always run infer_experiment.py and set -s accordingly in featureCounts.
• Batch effects disguised as biology. If samples were sequenced across lanes, always include a batch term in the DE design.
• Mismatched annotation between STAR’s --sjdbGTFfile and featureCounts’ -a produces phantom “unassigned” reads. Use one GTF everywhere.
• rRNA contamination. If more than ~10 % of reads map to rRNA loci, your poly(A) selection or rRNA depletion failed. Do not paper over it in silico.

Next steps

• Move to Salmon + tximport for transcript-level analysis and 5-10× faster runtimes.
• Add RSEM if you need TPM values for cross-lab comparisons.
• Learn single-cell RNA-Seq — the sample-preparation logic is different but a lot of the downstream statistics carry over. Read our scRNA-Seq workflow with Seurat next.

RNA-Seq is deceptively simple to run and deceptively easy to get wrong. Build your pipeline as if a stranger will have to reproduce your figures in three years — because they will, and that stranger is you.