RNA-seq Analysis

intermediate

RNA-seq Analysis Pipeline

RNA-seq from FASTQ to Pathways

🧪

RNA Extraction & QC

RIN ≥7 required; NanoDrop A260/A280 ≥1.8; Bioanalyzer/TapeStation

📚

Library Prep

Poly-A selection (mRNA) or ribo-depletion (all RNA); dUTP strand-specific

⭐

STAR Alignment

Splice-aware; 20-30M reads/sample; --quantMode GeneCounts

📊

Count Matrix

featureCounts or HTSeq; rows=genes, cols=samples; raw integer counts

📈

DESeq2 Analysis

Normalisation → model fitting → Wald test; padj<0.05, |log2FC|>1

🗺

Pathway Analysis

GO enrichment (clusterProfiler); GSEA; KEGG; STRING network

RNA-seq Workflow

RNA Extraction

QC (RIN ≥7)

Library Prep (poly-A or ribo-depletion)

Sequencing

FASTQ QC

Trimming

Alignment (STAR)

Count Matrix (featureCounts)

Normalisation

Differential Expression (DESeq2)

Pathway Analysis

Poly-A selection: enriches mRNA; removes rRNA (>95% of total RNA); standard for transcriptomics
Ribosomal depletion: keeps non-polyadenylated transcripts (lncRNA, histone mRNAs); better for degraded samples
Strand specificity: dUTP method (most common); always check stranded-ness before analysis
Sequencing depth: 20–30M reads for differential expression; 50M+ for transcript discovery
Replicates: minimum 3 biological replicates per condition; 6+ for clinical studies

STAR Alignment

code

# Build genome index (once)
STAR --runMode genomeGenerate \
  --genomeDir STAR_index/ \
  --genomeFastaFiles hg38.fa \
  --sjdbGTFfile gencode.v45.annotation.gtf \
  --runThreadN 16

# Align reads
STAR --runThreadN 16 \
  --genomeDir STAR_index/ \
  --readFilesIn R1.fastq.gz R2.fastq.gz \
  --readFilesCommand zcat \
  --outSAMtype BAM SortedByCoordinate \
  --outSAMattributes NH HI AS NM MD \
  --outFileNamePrefix sample_ \
  --quantMode GeneCounts \
  --outFilterMismatchNmax 2

DESeq2 Differential Expression

code

library(DESeq2)
library(ggplot2)

# Load count matrix
counts <- read.table("counts.txt", header=TRUE, row.names=1)
metadata <- data.frame(
  condition = c("control","control","control","treated","treated","treated"),
  row.names = colnames(counts)
)

# Create DESeq object
dds <- DESeqDataSetFromMatrix(
  countData = counts,
  colData = metadata,
  design = ~ condition
)

# Pre-filter low-count genes
dds <- dds[rowSums(counts(dds)) >= 10, ]

# Run analysis
dds <- DESeq(dds)
res <- results(dds, contrast=c("condition","treated","control"))
res <- lfcShrink(dds, coef="condition_treated_vs_control", type="apeglm")

# Filter significant DEGs
sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
write.csv(as.data.frame(sig), "DEGs.csv")

# Volcano plot
ggplot(as.data.frame(res), aes(log2FoldChange, -log10(padj))) +
  geom_point(aes(colour = padj < 0.05 & abs(log2FoldChange) > 1)) +
  theme_bw()

Pathway Analysis

Gene Ontology (GO): Biological Process, Molecular Function, Cellular Component enrichment
KEGG pathways: disease and metabolic pathway mapping; KEGGREST R package
GSEA (Gene Set Enrichment Analysis): uses ranked gene list; more powerful than over-representation
clusterProfiler (R): comprehensive enrichment analysis; enrichGO(), enrichKEGG(), gseGO()
MSigDB: curated gene sets; hallmark (50 pathways), C2 (chemical/genetic), C5 (ontology)
STRING: protein-protein interaction network; identifies hub genes
Upstream regulator analysis (IPA): identifies transcription factors driving observed changes

ACMG Variant Interpretation Single-cell Genomics