Experiment 06TranscriptomicsDifferential Expression
Next-Generation RNA-Seq Differential Expression Analysis: Quality Trimming (FastQC), Pseudoalignment & DESeq2 Volcano Plot Modeling
Duration: 6 Hours • Command-Line & R Bioconductor Pipeline
🎯 Experimental Aim & Learning Objectives
To evaluate raw Illumina RNA-Seq FASTQ read quality using FastQC and Trimmomatic; quantify transcript abundance without full alignment using Salmon/Kallisto pseudoalignment; model count data dispersion in R using DESeq2 negative binomial generalized linear models; normalize for library depth via median-of-ratios; apply the Benjamini-Hochberg false discovery rate (FDR) adjustment; and visualize differential expression via interactive volcano plots and heatmaps.
🔬 Theoretical Foundation & Statistical Models
RNA-Seq provides quantitative counts of transcript abundance across biological conditions. Due to biological variability across replicates, count data exhibits overdispersion ($text{Variance} > text{Mean}$), violating Poisson assumptions:
- Phred Quality Score ($Q$): Logarithmic measure of base-calling error probability:
$$Q = -10 log_{10}(P_{text{error}}) implies Q30 equiv 99.9% text{ base accuracy (1 in 1000 error)}.$$ - Transcript Pseudoalignment (Salmon / Kallisto): Matches reads directly to transcript De Bruijn graph target k-mers, skipping full base-by-base dynamic programming alignment. Outputs estimated counts and Transcripts Per Million (TPM):
$$text{TPM}_i = frac{frac{C_i}{L_i}}{sum_j frac{C_j}{L_j}} times 10^6$$ - DESeq2 Negative Binomial GLM: The raw read count $K_{ij}$ for gene $i$ in sample $j$ is modeled as:
$$K_{ij} sim text{NB}(mu_{ij}, alpha_i), quad operatorname{Var}(K_{ij}) = mu_{ij} + alpha_i mu_{ij}^2$$
where $alpha_i$ is the gene-specific dispersion parameter. - Size Factor Normalization (Median-of-Ratios): Corrects for sequencing depth and library composition:
$$hat{s}_j = operatorname{median}_{i} frac{K_{ij}}{left(prod_{v=1}^m K_{iv}right)^{1/m}}$$ - Wald Test with FDR Correction: Tests hypotheses: $H_0: log_2text{FC} = 0$. $P$-values are adjusted using the Benjamini-Hochberg procedure to control False Discovery Rate ($p_{text{adj}} < 0.05$).
Figure 1.6: Complete RNA-Seq computational pipeline from FastQC quality assessment through Salmon/Kallisto pseudoalignment to DESeq2 negative binomial dispersion modeling and volcano plot visualization.
🧪 Computational Tools & Libraries
- Command-Line Tools: FastQC 0.12, Trimmomatic 0.39, Salmon 1.10 (pseudoaligner).
- R Bioconductor Stack: R 4.3+,
DESeq2,tximport,EnhancedVolcano,pheatmap,ggplot2. - Benchmark RNA-Seq Dataset: Paired-end Illumina FASTQ reads from human lung epithelial cells: Control ($n=3$) vs TGF-$beta$ treated ($n=3$, epithelial-mesenchymal transition model).
📋 Step-by-Step Bioinformatic Protocol
- Read Quality Inspection: Run FastQC on all FASTQ files:
fastqc sample_*.fastq.gz -o fastqc_out/. Inspect per-base sequence quality curves, adapter contamination, and GC content. - Transcript Quantification with Salmon: Generate target transcript index:
salmon index -t gencode.v44.transcripts.fa.gz -i gencode_index -k 31
Quantify each biological replicate:
salmon quant -i gencode_index -l A -1 rep1_R1.fq.gz -2 rep1_R2.fq.gz --validateMappings -o quant_rep1/ - DESeq2 R Pipeline: Import quantifications using
tximport:
txi <- tximport(files, type="salmon", tx2gene=tx2gene)dds <- DESeqDataSetFromTximport(txi, colData=metadata, design= ~ condition)dds <- DESeq(dds)res <- results(dds, contrast=c("condition", "treated", "control"), alpha=0.05) - Dispersion & Volcano Plot: Plot dispersion estimates. Generate volcano plot:
EnhancedVolcano(res, lab=rownames(res), x='log2FoldChange', y='padj', pCutoff=0.05, FCcutoff=1.5)
📊 Representative Differential Expression Output
| Gene Symbol | Base Mean | $log_2$ Fold Change | lfcSE | Wald Stat | $p$-value | Adjusted $p$-value ($p_{text{adj}}$) |
|---|---|---|---|---|---|---|
| CDH2 (N-Cadherin) | 2480.5 | +3.42 | 0.21 | 16.28 | $1.4 times 10^{-59}$ | $8.2 times 10^{-56}$ (Up) |
| VIM (Vimentin) | 8950.2 | +2.85 | 0.19 | 15.00 | $7.1 times 10^{-51}$ | $2.4 times 10^{-47}$ (Up) |
| CDH1 (E-Cadherin) | 4120.8 | −3.18 | 0.24 | -13.25 | $4.5 times 10^{-40}$ | $1.1 times 10^{-36}$ (Down) |
| ACTB (Beta-actin) | 35120.0 | +0.04 | 0.08 | 0.50 | 0.617 | 0.742 (NS Housekeeping) |
Biological Assessment: Pronounced upregulation of mesenchymal markers (N-cadherin, Vimentin) and downregulation of epithelial E-cadherin confirms canonical Epithelial-to-Mesenchymal Transition (EMT) induction by TGF-$beta$ ($p_{text{adj}} < 10^{-35}$).
{reading_nav}
