omics.co.in logo omics.co.in

4. Bulk RNA-Seq Quantification & Differential Expression

⚡ Pipeline 4 of 16  ·  🧬 RNA & Transcriptomics

End-to-end transcriptomic profiling from raw paired-end sequencing through splice-aware alignment with STAR 2-Pass mode, transcript quantification via Salmon/RSEM, multi-factor statistical modeling with DESeq2, and GSEA pathway discovery.

🔍 Ctrl + K
RNA • Transcriptomics Authority: ENCODE / nf-core/rnaseq / DESeq2
Production Blueprint v4.5

4. Bulk RNA-Seq Quantification & Differential Expression

End-to-end transcriptomic profiling from raw paired-end sequencing through splice-aware alignment with STAR 2-Pass mode, transcript quantification via Salmon/RSEM, multi-factor statistical modeling with DESeq2, and GSEA pathway discovery.

Metric / RequirementSpecificationTarget Configuration
Core Engine / FrameworkSTAR 2-Pass + RSEM + Salmon + DESeq2Nextflow / Snakemake Compatible
Benchmark Runtime~3 – 5 hours (50M PE reads per sample)16 CPU cores / 40 GB RAM
Standard Containerquay.io/biocontainers/star:2.7.11b--h43eeafb_0Docker • Singularity • Conda
Input → DeliverablePaired-end RNA FASTQs → Gene Counts, Normalized FPKM/TPM, DEG Tables, Volcano & GSEA Plots
Workflow Chain: FASTQ → fastp → STAR 2-Pass → RSEM / featureCounts → Salmon (quasi-mapping) → DESeq2 → clusterProfiler
⚡

Step 1: Quality Profiling, Poly-A / Adapter Trimming with fastp

fastp

Removes Illumina TruSeq/Nextera adapter sequences, clips low quality read tails (Phred < 20), trims poly-A tails, and evaluates duplication levels.

fastp -i raw_R1.fq.gz -I raw_R2.fq.gz \
      -o clean_R1.fq.gz -O clean_R2.fq.gz \
      --detect_adapter_for_pe \
      --trim_poly_x --trim_poly_g \
      --length_required 36 --thread 8 \
      --html rna_fastp.html --json rna_fastp.json
Key Flags & Notes: --trim_poly_x trims poly-A and poly-T tails common in eukaryotic mRNA-seq.
⚡

Step 2: Splice-Aware Genomic Alignment with STAR 2-Pass Mode

STAR 2-Pass

Aligns RNA-seq reads across splice junctions. The 2-pass mode discovers novel, non-annotated splice junctions in pass 1 and re-aligns all reads in pass 2 with maximum sensitivity.

STAR --genomeDir /references/star_index_grch38/ \
     --readFilesIn clean_R1.fq.gz clean_R2.fq.gz \
     --readFilesCommand zcat \
     --runThreadN 16 \
     --twopassMode Basic \
     --outSAMtype BAM SortedByCoordinate \
     --quantMode TranscriptomeSAM GeneCounts \
     --outWigType bedGraph \
     --outFileNamePrefix star_out/sample_
Key Flags & Notes: --quantMode TranscriptomeSAM GeneCounts simultaneously produces genomic BAM, transcriptomic BAM for RSEM, and raw exon gene count matrices.
⚡

Step 3: High-Precision Transcript Abundance Quantification with Salmon

Salmon v1.10+

Performs alignment-free selective quasi-mapping directly against the cDNA transcriptome with GC-bias, sequence-specific bias, and positional fragment corrections.

salmon quant -i /references/salmon_index_grch38/ \
             -l A -1 clean_R1.fq.gz -2 clean_R2.fq.gz \
             --gcBias --seqBias \
             --validateMappings \
             -p 16 -o salmon_quant/sample_01
Key Flags & Notes: --gcBias corrects for severe cDNA amplification dropouts in high/low GC transcript segments; -l A auto-detects strandedness.
⚡

Step 4: Differential Gene Expression Modeling with DESeq2 in R

DESeq2 • Bioconductor

Models read counts using Negative Binomial generalized linear models (GLMs), estimates dispersion parameters, and executes Wald significance tests.

library(DESeq2)
library(tximport)

# Ingest Salmon transcript abundances into gene-level counts
tx2gene <- read.csv("tx2gene.csv")
files <- c("sample1" = "salmon_quant/sample_01/quant.sf", "sample2" = "salmon_quant/sample_02/quant.sf")
txi <- tximport(files, type = "salmon", tx2gene = tx2gene)

# Experimental Design Matrix
colData <- data.frame(condition = factor(c("Control", "Treated")), row.names = names(files))
dds <- DESeqDataSetFromTximport(txi, colData = colData, design = ~ condition)

# Filter low count genes and run Wald test
dds <- dds[rowSums(counts(dds)) >= 10, ]
dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "Treated", "Control"), alpha = 0.05)
res_sig <- subset(res, padj < 0.05 & abs(log2FoldChange) >= 1.0)
write.csv(as.data.frame(res_sig), file = "significant_DEGs.csv")
Key Flags & Notes: tximport correctly accounts for effective transcript length changes across conditions; padj < 0.05 enforces Benjamini-Hochberg FDR control.
⚡

Step 5: Pathway Over-Representation Analysis & GSEA Enrichment

clusterProfiler • fgsea

Uncovers enriched biological processes, KEGG signaling cascades, and Reactome pathways among differentially expressed genes.

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

gene_list <- res$log2FoldChange
names(gene_list) <- rownames(res)
gene_list <- sort(gene_list, decreasing = TRUE)

# Gene Set Enrichment Analysis (GSEA) on GO Biological Processes
gsea_results <- gseGO(geneList = gene_list,
                      OrgDb = org.Hs.eg.db,
                      keyType = "SYMBOL",
                      ont = "BP",
                      pvalueCutoff = 0.05,
                      verbose = FALSE)

write.csv(as.data.frame(gsea_results), "gsea_pathways.csv")
Key Flags & Notes: GSEA evaluates all ranked genes without arbitrary log-fold cutoff thresholds.

🚀 Automated Nextflow (nf-core) Production Execution

Execute this complete workflow across any SLURM, AWS Batch, or Google Cloud HPC cluster using the verified community container pipeline:

nextflow run nf-core/rnaseq -r 3.14.0 \
  --input samplesheet.csv \
  --genome GRCh38 \
  --aligner star_salmon \
  --outdir ./results_rnaseq/ \
  -profile docker
← Return to Home Hub