🧠 RNA-Seq Fundamentals

Q: What is RNA-Seq?

RNA-Seq (RNA sequencing) is a technique that uses next-generation sequencing to measure the presence and quantity of RNA in a biological sample. It reveals gene expression levels, alternative splicing, and novel transcripts, enabling transcriptome-wide analysis.

Q: What is the difference between alignment-based and alignment-free quantification?

Alignment-based methods (STAR, HISAT2 + featureCounts) align reads to a reference genome before counting. Alignment-free methods (Salmon, kallisto) use k-mer matching to a transcriptome, which is much faster and memory-efficient, with comparable accuracy for differential expression.

Q: Which differential expression tool should I use?

Use DESeq2 for most bulk RNA-Seq analyses (robust, well-documented), edgeR for small sample sizes or complex designs, and limma-voom for large datasets or when speed is critical. All three are R/Bioconductor packages and produce similar results for typical experiments.

Q: How many biological replicates do I need for RNA-Seq?

A minimum of 3 biological replicates per condition is strongly recommended. More replicates increase statistical power and reduce false positives. For cell lines, 3-4 is common; for human samples, 10+ is ideal. Technical replicates are not a substitute for biological replicates.

Q: What is the difference between bulk and single-cell RNA-Seq?

Bulk RNA-Seq measures average expression across all cells in a sample. Single-cell RNA-Seq (scRNA-Seq) measures expression in individual cells, revealing cell-type heterogeneity, rare cell populations, and developmental trajectories. scRNA-Seq requires specialized tools like Seurat, Scanpy, and Cell Ranger.

🎯 What is RNA-Seq?

RNA-Seq (RNA sequencing) is a technique that uses next-generation sequencing to measure the presence and quantity of RNA in a biological sample. It reveals gene expression levels, alternative splicing, and novel transcripts, enabling transcriptome-wide analysis. It's a cornerstone of functional genomics, disease research, and drug discovery.

Why RNA-Seq Matters

Expression: Gene levels. Splicing: Isoforms. Fusion: Cancer. Variants: eQTL. Discovery: Novel transcripts. Clinical: Diagnostics.

📌 Gene Expression📌 Splicing📌 Fusion Genes📌 Differential Expression

# RNA-Seq Analysis — the big picture
# From raw reads to biological insight.

# CONCEPTUAL OVERVIEW
print("""
=== RNA-Seq Analysis ===

What it is:
  - Sequence RNA (mRNA, total RNA)
  - Quantify gene expression
  - Detect splicing, fusions
  - Compare conditions

Why it's powerful:
  - Whole transcriptome
  - Quantitative
  - Sensitive
  - Discovery-driven
  - Clinically relevant

Types of RNA-Seq:
  - Bulk RNA-Seq: Average expression
  - Single-cell: Per-cell expression
  - Spatial: Location + expression
  - Long-read: Full isoforms
  - Targeted: Specific genes

Workflow:
  1. Experimental design
  2. Library prep
  3. Sequencing
  4. QC (FastQC)
  5. Trimming (fastp)
  6. Alignment (STAR, HISAT2)
  7. Quantification (featureCounts, Salmon)
  8. Differential expression (DESeq2)
  9. Pathway analysis (GSEA)
  10. Visualization

Applications:
  - Disease vs normal
  - Drug response
  - Development
  - Cancer subtypes
  - Biomarker discovery
  - Single-cell atlases
""")

# WHY IT MATTERS
print("=== Why It Matters ===")
print("""
Biological questions:
  - Which genes are differentially expressed?
  - What pathways are affected?
  - What isoforms are present?
  - Are there fusion genes?
  - What cell types are present?

RNA-Seq answers:
  - Expression matrix
  - Fold changes
  - P-values
  - Pathway enrichment
  - Cell clusters
""")

# HISTORICAL CONTEXT
print("=== Historical Context ===")
print("""
1990s: SAGE, MPSS
2000s: Microarrays
2006: First RNA-Seq (Solexa)
2008: Mortazavi et al. (mouse)
2010: TopHat, Cufflinks
2012: STAR, HTSeq
2013: DESeq2, edgeR
2014: Salmon, kallisto
2015: Drop-seq, 10x Genomics
2020s: Single-cell atlases, spatial
""")

# CHALLENGES
print("=== Challenges ===")
print("""
1. EXPERIMENTAL DESIGN
   - Replicates
   - Batch effects
   - Confounding

2. LIBRARY PREP
   - rRNA depletion
   - Strandedness
   - Fragment size

3. COMPUTATIONAL
   - Alignment (splice-aware)
   - Quantification
   - Normalization

4. STATISTICS
   - Multiple testing
   - Low counts
   - Overdispersion

5. INTERPRETATION
   - Biological meaning
   - Pathways
   - Validation
""")

🎯 Why RNA-Seq Matters

🎯 The Value of RNA-Seq

RNA-Seq is the bridge between genotype and phenotype. It reveals which genes are active and how they change.

Biological Insights

Expression: Levels. Splicing: Isoforms. Fusion: Cancer. eQTL: Variants. Discovery: Novel. Clinical: Diagnostics.

🎯 Expression🎯 Splicing🎯 Fusion🎯 Clinical

# Why RNA-Seq Matters
print("=== 1. Gene Expression ===")
print("""
Expression:
  - Which genes are on/off
  - How much
  - Compare conditions

Examples:
  - Disease vs normal
  - Drug treated vs control
  - Time course
  - Tissue-specific

Output:
  - Count matrix
  - TPM/FPKM
  - Fold change
  - P-value
""")

print("\n=== 2. Differential Expression ===")
print("""
Differential expression:
  - Genes changing between conditions
  - Statistical testing
  - Multiple testing correction

Tools:
  - DESeq2
  - edgeR
  - limma-voom
  - sleuth

Output:
  - Log2 fold change
  - Adjusted p-value
  - Ranking
""")

print("\n=== 3. Alternative Splicing ===")
print("""
Splicing:
  - Exon skipping
  - Intron retention
  - Alternative 5'/3' splice sites
  - Mutually exclusive exons

Tools:
  - rMATS
  - SUPPA2
  - MAJIQ
  - LeafCutter

Output:
  - PSI (percent spliced in)
  - Isoform ratios
""")

print("\n=== 4. Fusion Genes ===")
print("""
Fusion:
  - Two genes joined
  - Cancer driver
  - Target

Examples:
  - BCR-ABL (CML)
  - EML4-ALK (lung)
  - TMPRSS2-ERG (prostate)

Tools:
  - STAR-Fusion
  - Arriba
  - SOAPfuse
  - FusionCatcher
""")

print("\n=== 5. eQTL Analysis ===")
print("""
eQTL:
  - Genetic variant → expression
  - Regulatory variation
  - Disease mechanism

Workflow:
  1. Genotype
  2. Expression (RNA-Seq)
  3. Association
  4. Fine mapping

Tools:
  - Matrix eQTL
  - FastQTL
  - QTLtools
""")

print("\n=== 6. Single-Cell ===")
print("""
Single-cell:
  - Per-cell expression
  - Heterogeneity
  - Rare cells
  - Trajectories

Tools:
  - Seurat (R)
  - Scanpy (Python)
  - Cell Ranger
  - Monocle

Output:
  - Cell clusters
  - Marker genes
  - UMAP/tSNE
""")

print("\n=== 7. Clinical ===")
print("""
Clinical:
  - Cancer diagnostics
  - Fusion detection
  - Expression signatures
  - Drug response

Examples:
  - Oncotype DX
  - MammaPrint
  - FoundationOne
  - Caris

Tools:
  - STAR-Fusion
  - Arriba
  - DESeq2
""")

print("\n=== 8. Drug Discovery ===")
print("""
Drug discovery:
  - Mechanism of action
  - Target identification
  - Biomarkers
  - Resistance

Workflow:
  1. Treat cells
  2. RNA-Seq
  3. DE analysis
  4. Pathway
  5. Validation
""")

⚙️ Workflow Overview

🎯 From Reads to Insight

The RNA-Seq workflow has well-defined steps: QC, alignment/quantification, DE analysis, and interpretation.

Workflow

QC: FastQC. Trim: fastp. Align: STAR. Quantify: featureCounts. DE: DESeq2. Pathway: GSEA. Visualize: ggplot2.

⚙️ QC⚙️ Align⚙️ Quantify⚙️ DE

# RNA-Seq Workflow
print("=== Complete Pipeline ===")
print("""
1. EXPERIMENTAL DESIGN
   - Biological replicates (n ≥ 3)
   - Randomization
   - Batch control
   - Power analysis

2. LIBRARY PREP
   - mRNA selection (poly-A)
   - rRNA depletion
   - Stranded or unstranded
   - Paired-end or single-end

3. SEQUENCING
   - Illumina (short reads)
   - 20-30M reads per sample
   - 2×100 bp or 2×150 bp

4. QUALITY CONTROL
   - FastQC
   - MultiQC (aggregate)
   - Check: per-base quality, GC, adapters

5. TRIMMING
   - fastp / Trimmomatic
   - Remove adapters
   - Quality filtering
   - Length filtering

6. ALIGNMENT
   - STAR (fast, splice-aware)
   - HISAT2 (memory-efficient)
   - Salmon (alignment-free)

7. QUANTIFICATION
   - featureCounts (alignment-based)
   - Salmon / kallisto (alignment-free)
   - HTSeq (Python)

8. DIFFERENTIAL EXPRESSION
   - DESeq2 (R)
   - edgeR (R)
   - limma-voom (R)

9. PATHWAY ANALYSIS
   - GSEA
   - clusterProfiler
   - Reactome
   - GO enrichment

10. VISUALIZATION
    - PCA, heatmap
    - Volcano, MA
    - Boxplots, violin
""")

print("\n=== Two Main Approaches ===")
print("""
A. Alignment-based:
   Reads → STAR → BAM → featureCounts → counts → DESeq2
   - Accurate
   - Slower
   - More disk space
   - Splice-aware

B. Alignment-free:
   Reads → Salmon → TPM → DESeq2 (with tximport)
   - Fast
   - Memory-efficient
   - Good accuracy
   - No BAM
""")

print("\n=== Key Decisions ===")
print("""
1. Alignment vs alignment-free?
   → Salmon for speed
   → STAR for accuracy + novel junctions

2. Paired vs single?
   → Paired-end (better)

3. Stranded vs unstranded?
   → Stranded (better)

4. Poly-A vs rRNA depletion?
   → Poly-A for mRNA
   → rRNA depletion for total RNA

5. Replicates?
   → n ≥ 3 per condition
""")

🧹 QC & Trimming

🎯 Clean Before You Analyze

Quality control and trimming are critical for accurate expression estimates. Garbage in, garbage out.

QC

FastQC: Raw reads. MultiQC: Aggregate. fastp: Trim. Adapters: Remove. Quality: Filter. rRNA: Check.

🧹 FastQC🧹 MultiQC🧹 fastp🧹 rRNA

# QC & Trimming
print("=== FastQC ===")
print("""
FastQC:
  - Per-base quality
  - GC content
  - Adapter contamination
  - Duplication
  - Overrepresented sequences

Command:
  fastqc reads_1.fastq.gz reads_2.fastq.gz -o qc/

Output:
  reads_1_fastqc.html
  reads_1_fastqc.zip
""")

print("\n=== MultiQC ===")
print("""
MultiQC:
  - Aggregate all QC
  - One report
  - Compare samples

Command:
  multiqc qc/ -o multiqc_report/

Output:
  multiqc_report.html
""")

print("\n=== fastp ===")
print("""
fastp:
  - Fast
  - All-in-one
  - Adapter detection
  - Quality filtering
  - Deduplication

Command:
  fastp \
    -i reads_1.fastq.gz \
    -I reads_2.fastq.gz \
    -o clean_1.fastq.gz \
    -O clean_2.fastq.gz \
    --detect_adapter_for_pe \
    --qualified_quality_phred 20 \
    --length_required 36 \
    --json report.json \
    --html report.html

Key options:
  -q 20 (quality)
  -l 36 (min length)
  -u 40 (max uncalled)
  --cut_front, --cut_tail
""")

print("\n=== Trimmomatic ===")
print("""
Trimmomatic:
  - Java
  - Flexible
  - Popular

Command:
  trimmomatic PE \
    reads_1.fastq.gz reads_2.fastq.gz \
    clean_1.fq.gz unpaired_1.fq.gz \
    clean_2.fq.gz unpaired_2.fq.gz \
    ILLUMINACLIP:adapters.fa:2:30:10 \
    LEADING:3 TRAILING:3 \
    SLIDINGWINDOW:4:20 MINLEN:36
""")

print("\n=== rRNA Contamination ===")
print("""
rRNA:
  - Major contaminant
  - Use rRNA depletion
  - Check with FastQ Screen

Tools:
  - FastQ Screen
  - SortMeRNA
  - Bowtie2 (to rRNA)

Check:
  - < 5% rRNA (good)
  - > 20% rRNA (bad)
""")

print("\n=== Best Practices ===")
print("""
1. Always run FastQC
2. Aggregate with MultiQC
3. Trim adapters
4. Filter low quality
5. Check rRNA
6. Check duplication
7. Remove overrepresented
8. Document everything
""")

🗺️ Alignment (STAR, HISAT2)

🎯 Mapping Reads to the Genome

Splice-aware alignment is essential for RNA-Seq. STAR and HISAT2 are the standard tools.

Alignment

STAR: Fast, accurate. HISAT2: Memory-efficient. Splice: Aware. Index: Build once. BAM: Output. Sort: samtools.

🗺️ STAR🗺️ HISAT2🗺️ Splice-aware🗺️ BAM

# Alignment
print("=== STAR ===")
print("""
STAR (Spliced Transcripts Alignment to a Reference):
  - Fast
  - Accurate
  - Splice-aware
  - Standard

Build index:
  STAR --runMode genomeGenerate \
    --genomeDir star_index \
    --genomeFastaFiles genome.fa \
    --sjdbGTFfile annotation.gtf \
    --sjdbOverhang 99 \
    --runThreadN 8

Align:
  STAR --genomeDir star_index \
    --readFilesIn clean_1.fq.gz clean_2.fq.gz \
    --readFilesCommand zcat \
    --runThreadN 8 \
    --outFileNamePrefix sample1_ \
    --outSAMtype BAM SortedByCoordinate \
    --outSAMattributes NH HI AS NM MD \
    --quantMode GeneCounts

Output:
  sample1_Aligned.sortedByCoord.out.bam
  sample1_ReadsPerGene.out.tab
""")

print("\n=== HISAT2 ===")
print("""
HISAT2:
  - Memory-efficient
  - Fast
  - Splice-aware
  - Good for large genomes

Build index:
  hisat2-build genome.fa hisat2_index

Align:
  hisat2 -p 8 \
    -x hisat2_index \
    -1 clean_1.fq.gz \
    -2 clean_2.fq.gz \
    --rna-strandness RF \
    -S sample1.sam

Convert to BAM:
  samtools sort -@ 8 -o sample1.bam sample1.sam
  samtools index sample1.bam
""")

print("\n=== STAR vs HISAT2 ===")
print(f"{'Feature':<25} {'STAR':<20} {'HISAT2':<20}")
print("-" * 65)
print(f"{'Speed':<25} {'Very fast':<20} {'Fast':<20}")
print(f"{'Memory':<25} {'High (30 GB)':<20} {'Low (5 GB)':<20}")
print(f"{'Accuracy':<25} {'Excellent':<20} {'Very good':<20}")
print(f"{'Novel junctions':<25} {'Yes':<20} {'Yes':<20}")
print(f"{'Output':<25} {'BAM':<20} {'SAM':<20}")
print(f"{'Best for':<25} {'Large genomes':<20} {'Limited RAM':<20}")

print("\n=== Alignment Stats ===")
print("""
Check:
  - Alignment rate > 70%
  - Unique alignment > 60%
  - rRNA < 5%
  - Exonic > 80%

samtools flagstat sample1.bam

Output:
  total: 20,000,000
  mapped: 18,000,000 (90%)
  properly paired: 17,500,000 (87.5%)
  duplicates: 2,000,000 (10%)
""")

print("\n=== Best Practices ===")
print("""
1. Use splice-aware aligner
2. Use stranded protocol
3. Build index with GTF
4. Sort BAM
5. Index BAM
6. Mark duplicates (optional)
7. Check alignment stats
""")

🔢 Quantification

🎯 Counting Reads per Gene

Quantification converts aligned reads into expression counts. Methods: featureCounts, HTSeq, Salmon, kallisto.

Quantification

featureCounts: Fast, alignment-based. HTSeq: Python, slow. Salmon: Alignment-free. kallisto: Fast. TPM: Normalized. Counts: For DE.

🔢 featureCounts🔢 Salmon🔢 TPM🔢 Counts

# Quantification
print("=== featureCounts ===")
print("""
featureCounts (Subread):
  - Fast
  - Alignment-based
  - Multi-threaded
  - Standard

Command:
  featureCounts -T 8 \
    -a annotation.gtf \
    -o counts.txt \
    -t exon \
    -g gene_id \
    -s 2 \
    sample1.bam sample2.bam sample3.bam

Options:
  -t exon (feature type)
  -g gene_id (group)
  -s 2 (stranded: 0=no, 1=yes, 2=reverse)
  -p (paired-end)
  -M (multi-mapping)
""")

print("\n=== HTSeq ===")
print("""
HTSeq:
  - Python
  - Slow
  - Accurate

Command:
  htseq-count \
    -f bam \
    -r pos \
    -s reverse \
    -t exon \
    -i gene_id \
    sample1.bam annotation.gtf > sample1.counts

Options:
  -s yes/no/reverse
  -m union/intersection-strict
""")

print("\n=== Salmon ===")
print("""
Salmon:
  - Alignment-free
  - Fast
  - Accurate
  - TPM output

Build index:
  salmon index \
    -t transcripts.fa \
    -i salmon_index \
    -k 31

Quantify:
  salmon quant \
    -i salmon_index \
    -l A \
    -1 clean_1.fq.gz \
    -2 clean_2.fq.gz \
    -p 8 \
    --validateMappings \
    -o sample1_quant

Output:
  quant.sf (transcript TPM)
  aux_info/
""")

print("\n=== kallisto ===")
print("""
kallisto:
  - Ultra-fast
  - Alignment-free
  - TPM output

Build index:
  kallisto index -i kallisto_index transcripts.fa

Quantify:
  kallisto quant \
    -i kallisto_index \
    -o sample1_quant \
    -t 8 \
    clean_1.fq.gz clean_2.fq.gz

Output:
  abundance.tsv
  abundance.h5
""")

print("\n=== Normalization ===")
print("""
Raw counts:
  - Not comparable across samples
  - Need normalization

Methods:
  - TPM (transcripts per million)
  - FPKM (fragments per kb per million)
  - CPM (counts per million)
  - TMM (trimmed mean of M-values)
  - RPKM (reads per kb per million)

For DE:
  - DESeq2: median of ratios
  - edgeR: TMM
  - limma-voom: TMM
""")

print("\n=== TPM vs FPKM ===")
print(f"{'Feature':<20} {'TPM':<20} {'FPKM':<20}")
print("-" * 60)
print(f"{'Order':<20} {'Transcript':<20} {'Fragment':<20}")
print(f"{'Normalization':<20} {'Length + depth':<20} {'Length + depth':<20}")
print(f"{'Comparable':<20} {'Across samples':<20} {'Within sample':<20}")
print(f"{'Sum':<20} {'1M':<20} {'Variable':<20}")
print(f"{'Recommended':<20} {'Yes':<20} {'Deprecated':<20}")

print("\n=== Best Practices ===")
print("""
1. Use featureCounts (alignment)
2. Or Salmon (alignment-free)
3. Check strandedness
4. Use GTF annotation
5. Use gene_id
6. Check mapping rate
7. Normalize for DE
""")

📈 Differential Expression

🎯 Finding DE Genes

Differential expression (DE) identifies genes that change between conditions using statistical models that account for count data.

DE

DESeq2: Negative binomial. edgeR: Negative binomial. limma: Linear model. Fold change: Effect size. P-value: Significance. FDR: Multiple testing.

📈 DESeq2📈 edgeR📈 Fold change📈 FDR

# Differential Expression
print("=== DESeq2 ===")
print("""
DESeq2:
  - Negative binomial
  - Median of ratios
  - Shrinkage
  - Gold standard

R code:
  library(DESeq2)
  
  # Load counts
  counts <- read.table("counts.txt", header=TRUE, row.names=1)
  
  # Metadata
  coldata <- data.frame(
    condition = factor(c("control","control","control",
                         "treated","treated","treated"))
  )
  
  # DESeqDataSet
  dds <- DESeqDataSetFromMatrix(
    countData = counts,
    colData = coldata,
    design = ~ condition
  )
  
  # Run
  dds <- DESeq(dds)
  res <- results(dds, contrast=c("condition","treated","control"))
  
  # Order
  res <- res[order(res$padj),]
  
  # Significant
  sig <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
""")

print("\n=== edgeR ===")
print("""
edgeR:
  - Negative binomial
  - TMM normalization
  - Flexible

R code:
  library(edgeR)
  
  # DGEList
  dge <- DGEList(counts=counts, group=condition)
  
  # Filter
  keep <- filterByExpr(dge)
  dge <- dge[keep,,keep.lib.sizes=FALSE]
  
  # Normalize
  dge <- calcNormFactors(dge)
  
  # Design
  design <- model.matrix(~ condition)
  
  # Estimate
  dge <- estimateDisp(dge, design)
  fit <- glmQLFit(dge, design)
  
  # Test
  qlf <- glmQLFTest(fit, coef=2)
  topTags(qlf)
""")

print("\n=== limma-voom ===")
print("""
limma-voom:
  - Linear model
  - Fast
  - Large datasets

R code:
  library(limma)
  library(edgeR)
  
  # DGEList
  dge <- DGEList(counts=counts, group=condition)
  dge <- calcNormFactors(dge)
  
  # Design
  design <- model.matrix(~ condition)
  
  # voom
  v <- voom(dge, design)
  
  # Fit
  fit <- lmFit(v, design)
  fit <- eBayes(fit)
  
  # Results
  topTable(fit, coef=2, adjust="BH")
""")

print("\n=== Results Interpretation ===")
print("""
Output columns:
  - baseMean: Mean normalized counts
  - log2FoldChange: Effect size
  - lfcSE: Standard error
  - stat: Test statistic
  - pvalue: Raw p-value
  - padj: Adjusted p-value (FDR)

Thresholds:
  - padj < 0.05
  - |log2FC| > 1 (2-fold)
  - Or |log2FC| > 0.58 (1.5-fold)

Volcano plot:
  - X: log2FC
  - Y: -log10(padj)
""")

print("\n=== Batch Effects ===")
print("""
Batch effects:
  - Technical variation
  - Confounding

Detection:
  - PCA
  - Hierarchical clustering

Correction:
  - Include in design: ~ batch + condition
  - ComBat-seq (sva)
  - Remove batch effect
""")

print("\n=== Best Practices ===")
print("""
1. n ≥ 3 biological replicates
2. Filter low-count genes
3. Normalize (TMM or median ratios)
4. Use appropriate design
5. Correct for batch
6. Multiple testing correction
7. Effect size threshold
8. Validate with qPCR
""")

🛤️ Pathway Analysis

🎯 From Genes to Biology

Pathway analysis interprets DE gene lists by identifying enriched biological processes, pathways, and functions.

Pathway

GO: Gene ontology. KEGG: Pathways. Reactome: Reactions. GSEA: Ranked. ORA: Overlap. Enrichment: Significance.

🛤️ GO🛤️ KEGG🛤️ GSEA🛤️ Reactome

# Pathway Analysis
print("=== ORA (Over-Representation Analysis) ===")
print("""
ORA:
  - Input: DE gene list
  - Test: Overlap with pathway
  - Hypergeometric test

Tools:
  - clusterProfiler (R)
  - Enrichr (web)
  - DAVID (web)
  - WebGestalt

R code (clusterProfiler):
  library(clusterProfiler)
  
  # GO enrichment
  ego <- enrichGO(
    gene = de_genes,
    OrgDb = org.Hs.eg.db,
    ont = "BP",
    pAdjustMethod = "BH",
    pvalueCutoff = 0.05
  )
  
  # KEGG
  ekk <- enrichKEGG(
    gene = de_genes,
    organism = "hsa",
    pvalueCutoff = 0.05
  )
""")

print("\n=== GSEA (Gene Set Enrichment Analysis) ===")
print("""
GSEA:
  - Input: All genes ranked
  - No cutoff needed
  - More sensitive

Ranking:
  - log2FC
  - -log10(p) * sign(log2FC)
  - stat

R code:
  library(clusterProfiler)
  
  # Rank genes
  gene_list <- res$log2FoldChange
  names(gene_list) <- rownames(res)
  gene_list <- sort(gene_list, decreasing=TRUE)
  
  # GSEA
  gsea <- gseGO(
    geneList = gene_list,
    OrgDb = org.Hs.eg.db,
    ont = "BP",
    pvalueCutoff = 0.05
  )
""")

print("\n=== GO Terms ===")
print("""
GO:
  - BP: Biological Process
  - MF: Molecular Function
  - CC: Cellular Component

Example:
  BP: "immune response"
  MF: "kinase activity"
  CC: "nucleus"
""")

print("\n=== KEGG Pathways ===")
print("""
KEGG:
  - Metabolic pathways
  - Signaling
  - Disease

Example:
  - hsa04110: Cell cycle
  - hsa04151: PI3K-Akt
  - hsa04010: MAPK
""")

print("\n=== Reactome ===")
print("""
Reactome:
  - Curated pathways
  - Detailed
  - Hierarchical

Tools:
  - ReactomePA (R)
  - Reactome web
""")

print("\n=== Visualization ===")
print("""
Plots:
  - Dotplot
  - Barplot
  - Enrichment map
  - Network
  - Heatmap

R code:
  dotplot(ego, showCategory=20)
  barplot(ego, showCategory=20)
  emapplot(ego)
  cnetplot(ego)
""")

print("\n=== Best Practices ===")
print("""
1. Use all genes for GSEA
2. Use DE list for ORA
3. Multiple databases
4. Check enrichment score
5. Visualize networks
6. Validate key pathways
7. Report parameters
""")

🔬 Single-Cell RNA-Seq

🎯 Per-Cell Expression

Single-cell RNA-Seq (scRNA-Seq) measures expression in individual cells, revealing heterogeneity, cell types, and trajectories.

scRNA-Seq

10x: Droplet. Smart-seq: Plate. Seurat: R. Scanpy: Python. Clustering: Cell types. UMAP: Visualization.

🔬 10x🔬 Seurat🔬 Scanpy🔬 UMAP

# Single-Cell RNA-Seq
print("=== Platforms ===")
print("""
1. 10x Genomics Chromium
   - Droplet-based
   - 3' or 5'
   - Thousands of cells
   - Standard

2. Smart-seq2/3
   - Plate-based
   - Full-length
   - Hundreds of cells
   - Higher sensitivity

3. Drop-seq
   - Droplet
   - Academic

4. Seq-Well
   - Portable
""")

print("\n=== Workflow ===")
print("""
1. Cell isolation
2. Library prep
3. Sequencing
4. Cell Ranger (10x)
5. QC (Seurat/Scanpy)
6. Normalization
7. Clustering
8. Marker genes
9. Annotation
10. Trajectory
""")

print("\n=== Cell Ranger ===")
print("""
Cell Ranger:
  - 10x official
  - Fast
  - Standard

Command:
  cellranger count \
    --id=sample1 \
    --transcriptome=refdata-gex-GRCh38 \
    --fastqs=fastq_path \
    --sample=sample1 \
    --localcores=8 \
    --localmem=64

Output:
  filtered_feature_bc_matrix/
  raw_feature_bc_matrix/
  analysis/
""")

print("\n=== Seurat (R) ===")
print("""
Seurat:
  - R
  - Comprehensive
  - Standard

R code:
  library(Seurat)
  
  # Load
  data <- Read10X("filtered_feature_bc_matrix/")
  seurat <- CreateSeuratObject(data, project="sample1")
  
  # QC
  seurat[["percent.mt"]] <- PercentageFeatureSet(seurat, pattern="^MT-")
  seurat <- subset(seurat, subset = nFeature_RNA > 200 & 
                   nFeature_RNA < 2500 & percent.mt < 5)
  
  # Normalize
  seurat <- NormalizeData(seurat)
  seurat <- FindVariableFeatures(seurat)
  seurat <- ScaleData(seurat)
  
  # PCA
  seurat <- RunPCA(seurat)
  
  # Clustering
  seurat <- FindNeighbors(seurat, dims=1:10)
  seurat <- FindClusters(seurat, resolution=0.5)
  
  # UMAP
  seurat <- RunUMAP(seurat, dims=1:10)
  
  # Markers
  markers <- FindAllMarkers(seurat, only.pos=TRUE)
""")

print("\n=== Scanpy (Python) ===")
print("""
Scanpy:
  - Python
  - Fast
  - Scalable

Python code:
  import scanpy as sc
  
  # Load
  adata = sc.read_10x_mtx("filtered_feature_bc_matrix/")
  
  # QC
  adata.var_names_make_unique()
  sc.pp.filter_cells(adata, min_genes=200)
  sc.pp.filter_genes(adata, min_cells=3)
  
  # Normalize
  sc.pp.normalize_total(adata, target_sum=1e4)
  sc.pp.log1p(adata)
  
  # HVG
  sc.pp.highly_variable_genes(adata, n_top_genes=2000)
  
  # PCA
  sc.tl.pca(adata)
  
  # Neighbors
  sc.pp.neighbors(adata, n_pcs=10)
  
  # Clustering
  sc.tl.leiden(adata, resolution=0.5)
  
  # UMAP
  sc.tl.umap(adata)
  
  # Markers
  sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
  
  # Plot
  sc.pl.umap(adata, color=['leiden'])
""")

print("\n=== QC Metrics ===")
print("""
Per cell:
  - nFeature_RNA (genes)
  - nCount_RNA (UMIs)
  - percent.mt (mito %)
  - percent.ribo (ribo %)

Thresholds:
  - nFeature: 200-2500
  - nCount: 500-10000
  - percent.mt: < 5-20%

Doublets:
  - DoubletFinder
  - Scrublet
  - Solo
""")

print("\n=== Clustering ===")
print("""
Steps:
  1. Normalize
  2. HVG
  3. Scale
  4. PCA
  5. Neighbors
  6. Cluster (Leiden/Louvain)
  7. UMAP/tSNE

Resolution:
  - 0.3: Fewer clusters
  - 0.5: Standard
  - 1.0: More clusters
""")

print("\n=== Annotation ===")
print("""
Manual:
  - Marker genes
  - Known markers

Automatic:
  - SingleR
  - CellTypist
  - scType

Markers:
  - CD3D: T cells
  - CD79A: B cells
  - CD14: Monocytes
  - NKG7: NK cells
""")

print("\n=== Trajectory ===")
print("""
Trajectory:
  - Developmental
  - Differentiation
  - Disease progression

Tools:
  - Monocle3
  - Slingshot
  - PAGA
  - scVelo (RNA velocity)
""")

🐍 Python/R Integration

🎯 RNA-Seq from Code

You can run RNA-Seq pipelines from Python and R. Here's how.

Python/R

subprocess: Run CLI. pysam: BAM. pandas: Counts. DESeq2: R. Scanpy: Python. Seurat: R. Plot: Visualize.

🐍 subprocess🐍 pysam🐍 pandas🐍 Scanpy

# Python Integration
print("=== Run STAR ===")
print("""
import subprocess

def run_star(reads_1, reads_2, genome_dir, output_prefix,
             threads=8):
    \"\"\"Run STAR alignment.\"\"\"
    cmd = [
        'STAR',
        '--genomeDir', genome_dir,
        '--readFilesIn', reads_1, reads_2,
        '--readFilesCommand', 'zcat',
        '--runThreadN', str(threads),
        '--outFileNamePrefix', output_prefix,
        '--outSAMtype', 'BAM', 'SortedByCoordinate',
        '--quantMode', 'GeneCounts'
    ]
    
    result = subprocess.run(cmd, capture_output=True, text=True)
    
    if result.returncode != 0:
        raise RuntimeError(f"STAR failed: {result.stderr}")
    
    return f"{output_prefix}Aligned.sortedByCoord.out.bam"

# Use
# bam = run_star("r1.fq.gz", "r2.fq.gz", "star_index", "sample1_")
""")

print("\n=== Run featureCounts ===")
print("""
def run_featurecounts(bam_files, gtf, output,
                     threads=8, stranded=2):
    \"\"\"Run featureCounts.\"\"\"
    cmd = [
        'featureCounts',
        '-T', str(threads),
        '-a', gtf,
        '-o', output,
        '-t', 'exon',
        '-g', 'gene_id',
        '-s', str(stranded),
        '-p'
    ] + bam_files
    
    result = subprocess.run(cmd, capture_output=True, text=True)
    
    if result.returncode != 0:
        raise RuntimeError(f"featureCounts failed: {result.stderr}")
    
    return output

# Use
# counts = run_featurecounts(["s1.bam","s2.bam"], "genes.gtf", "counts.txt")
""")

print("\n=== Run Salmon ===")
print("""
def run_salmon(reads_1, reads_2, index, output_dir,
              threads=8):
    \"\"\"Run Salmon quantification.\"\"\"
    cmd = [
        'salmon', 'quant',
        '-i', index,
        '-l', 'A',
        '-1', reads_1,
        '-2', reads_2,
        '-p', str(threads),
        '--validateMappings',
        '-o', output_dir
    ]
    
    result = subprocess.run(cmd, capture_output=True, text=True)
    
    if result.returncode != 0:
        raise RuntimeError(f"Salmon failed: {result.stderr}")
    
    return f"{output_dir}/quant.sf"

# Use
# quant = run_salmon("r1.fq.gz", "r2.fq.gz", "salmon_index", "sample1_quant")
""")

print("\n=== Parse Counts ===")
print("""
import pandas as pd

def parse_counts(counts_file):
    \"\"\"Parse featureCounts output.\"\"\"
    # Skip first line
    counts = pd.read_csv(counts_file, sep='\\t', skiprows=1)
    
    # Set index
    counts = counts.set_index('Geneid')
    
    # Remove metadata columns
    counts = counts.drop(['Chr', 'Start', 'End', 'Strand', 'Length'],
                         axis=1)
    
    # Rename samples
    counts.columns = [c.replace('.bam','') for c in counts.columns]
    
    print(f"Genes: {len(counts)}")
    print(f"Samples: {len(counts.columns)}")
    
    return counts

# Use
# counts = parse_counts("counts.txt")
""")

print("\n=== DESeq2 from Python ===")
print("""
import subprocess
import pandas as pd

def run_deseq2(counts_file, metadata_file, output_file):
    \"\"\"Run DESeq2 via Rscript.\"\"\"
    r_script = f'''
    library(DESeq2)
    
    counts <- read.table("{counts_file}", header=TRUE, row.names=1)
    coldata <- read.table("{metadata_file}", header=TRUE, row.names=1)
    
    dds <- DESeqDataSetFromMatrix(
        countData = counts,
        colData = coldata,
        design = ~ condition
    )
    
    dds <- DESeq(dds)
    res <- results(dds)
    
    write.csv(as.data.frame(res), "{output_file}")
    '''
    
    with open("temp_deseq2.R", "w") as f:
        f.write(r_script)
    
    result = subprocess.run(['Rscript', 'temp_deseq2.R'],
                           capture_output=True, text=True)
    
    if result.returncode != 0:
        raise RuntimeError(f"DESeq2 failed: {result.stderr}")
    
    return pd.read_csv(output_file, index_col=0)

# Use
# res = run_deseq2("counts.txt", "metadata.txt", "deseq2_results.csv")
""")

print("\n=== Scanpy (Python) ===")
print("""
import scanpy as sc
import matplotlib.pyplot as plt

def analyze_scrna(data_dir):
    \"\"\"Analyze single-cell RNA-Seq.\"\"\"
    # Load
    adata = sc.read_10x_mtx(data_dir)
    adata.var_names_make_unique()
    
    # QC
    sc.pp.filter_cells(adata, min_genes=200)
    sc.pp.filter_genes(adata, min_cells=3)
    
    # Normalize
    sc.pp.normalize_total(adata, target_sum=1e4)
    sc.pp.log1p(adata)
    
    # HVG
    sc.pp.highly_variable_genes(adata, n_top_genes=2000)
    adata = adata[:, adata.var.highly_variable]
    
    # PCA
    sc.tl.pca(adata, n_comps=50)
    
    # Neighbors
    sc.pp.neighbors(adata, n_pcs=20)
    
    # Clustering
    sc.tl.leiden(adata, resolution=0.5)
    
    # UMAP
    sc.tl.umap(adata)
    
    # Markers
    sc.tl.rank_genes_groups(adata, 'leiden', method='wilcoxon')
    
    # Plot
    sc.pl.umap(adata, color=['leiden'], show=False)
    plt.savefig("umap.png")
    
    return adata

# Use
# adata = analyze_scrna("filtered_feature_bc_matrix/")
""")

print("\n=== Complete Pipeline ===")
print("""
import subprocess
import pandas as pd
import os

def rnaseq_pipeline(reads_1, reads_2, genome_dir, gtf,
                   output_dir="rnaseq_output", threads=8):
    \"\"\"Complete RNA-Seq pipeline.\"\"\"
    os.makedirs(output_dir, exist_ok=True)
    
    # 1. Trim
    print("Trimming...")
    subprocess.run([
        'fastp', '-i', reads_1, '-I', reads_2,
        '-o', f'{output_dir}/clean_1.fq.gz',
        '-O', f'{output_dir}/clean_2.fq.gz',
        '--detect_adapter_for_pe'
    ], check=True)
    
    # 2. Align
    print("Aligning...")
    run_star(f'{output_dir}/clean_1.fq.gz',
             f'{output_dir}/clean_2.fq.gz',
             genome_dir, f'{output_dir}/sample_',
             threads)
    
    # 3. Count
    print("Counting...")
    run_featurecounts(
        [f'{output_dir}/sample_Aligned.sortedByCoord.out.bam'],
        gtf, f'{output_dir}/counts.txt', threads
    )
    
    # 4. Parse
    counts = parse_counts(f'{output_dir}/counts.txt')
    
    print("Pipeline complete!")
    return counts

# Use
# counts = rnaseq_pipeline("r1.fq.gz", "r2.fq.gz", "star_index", "genes.gtf")
""")

🎨 Visualization

🎯 Seeing the Data

Visualization is critical for RNA-Seq interpretation: PCA, heatmaps, volcano plots, and more.

Visualization

PCA: Sample clustering. Heatmap: Expression. Volcano: DE. MA: Fold change. Boxplot: Distribution. UMAP: Single-cell.

🎨 PCA🎨 Heatmap🎨 Volcano🎨 UMAP

# Visualization
print("=== PCA ===")
print("""
PCA:
  - Sample clustering
  - Quality control
  - Batch detection

R code (DESeq2):
  vsd <- vst(dds, blind=FALSE)
  plotPCA(vsd, intgroup="condition")

Python:
  from sklearn.decomposition import PCA
  pca = PCA(n_components=2)
  pca_result = pca.fit_transform(log_counts.T)
  plt.scatter(pca_result[:,0], pca_result[:,1])
""")

print("\n=== Heatmap ===")
print("""
Heatmap:
  - Expression patterns
  - Clustering
  - Top genes

R code (pheatmap):
  library(pheatmap)
  
  top_genes <- head(order(res$padj), 50)
  pheatmap(assay(vsd)[top_genes,],
           cluster_rows=TRUE,
           cluster_cols=TRUE,
           annotation_col=coldata)

Python:
  import seaborn as sns
  sns.clustermap(log_counts.loc[top_genes])
""")

print("\n=== Volcano Plot ===")
print("""
Volcano:
  - X: log2FC
  - Y: -log10(padj)
  - Color: significance

R code (ggplot2):
  library(ggplot2)
  
  res$sig <- ifelse(res$padj < 0.05 & abs(res$log2FoldChange) > 1,
                    "Significant", "Not")
  
  ggplot(res, aes(x=log2FoldChange, y=-log10(padj), color=sig)) +
    geom_point() +
    theme_minimal()

Python:
  import matplotlib.pyplot as plt
  plt.scatter(res['log2FoldChange'], -np.log10(res['padj']))
""")

print("\n=== MA Plot ===")
print("""
MA:
  - X: Mean expression
  - Y: log2FC
  - Check bias

R code:
  plotMA(res, ylim=c(-5,5))
""")

print("\n=== Boxplot / Violin ===")
print("""
Boxplot:
  - Expression distribution
  - Per gene

R code:
  ggplot(data, aes(x=condition, y=expression, fill=condition)) +
    geom_boxplot() +
    theme_minimal()

Violin:
  ggplot(data, aes(x=condition, y=expression, fill=condition)) +
    geom_violin() +
    theme_minimal()
""")

print("\n=== UMAP/tSNE ===")
print("""
UMAP:
  - Single-cell
  - Clusters
  - Cell types

R code (Seurat):
  DimPlot(seurat, reduction="umap", group.by="seurat_clusters")

Python (Scanpy):
  sc.pl.umap(adata, color=['leiden'])
""")

print("\n=== Pathway Plots ===")
print("""
Dotplot:
  dotplot(ego, showCategory=20)

Barplot:
  barplot(ego, showCategory=20)

Network:
  emapplot(ego)
  cnetplot(ego)
""")

print("\n=== Best Practices ===")
print("""
1. Label axes clearly
2. Use colorblind-friendly
3. Include sample sizes
4. Show statistics
5. High resolution (300 dpi)
6. Vector format (PDF/SVG)
7. Consistent colors
8. Legend
""")

⚖️ Tool Comparison

🎯 Which Tool?

Choosing the right RNA-Seq tools depends on your goals, data size, and computational resources.

Comparison

STAR: Fast, accurate. HISAT2: Memory-efficient. Salmon: Alignment-free. DESeq2: Gold standard. edgeR: Flexible. limma: Fast.

⚖️ STAR⚖️ Salmon⚖️ DESeq2⚖️ edgeR

Tool	Type	Speed	Accuracy	Best For
STAR	Aligner	Fast	Excellent	General purpose
HISAT2	Aligner	Fast	Very good	Limited memory
Salmon	Quantifier	Very fast	Excellent	Fast quantification
kallisto	Quantifier	Very fast	Excellent	Fast quantification
featureCounts	Counter	Fast	Excellent	Alignment-based
DESeq2	DE	Medium	Excellent	Gold standard
edgeR	DE	Fast	Excellent	Small samples
limma-voom	DE	Very fast	Excellent	Large datasets
Seurat	scRNA	Medium	Excellent	Single-cell (R)
Scanpy	scRNA	Fast	Excellent	Single-cell (Python)

# Decision Guide
print("=== Choose Aligner ===")
print("""
STAR:
  ✓ Best accuracy
  ✓ Novel junctions
  ✓ Large genomes
  ✗ High memory

HISAT2:
  ✓ Low memory
  ✓ Fast
  ✓ Good accuracy
  ✗ Slightly less accurate

Salmon/kallisto:
  ✓ Fastest
  ✓ No BAM
  ✓ Good accuracy
  ✗ No novel junctions

Recommendation:
  → STAR for most
  → HISAT2 for limited RAM
  → Salmon for speed
""")

print("\n=== Choose DE Tool ===")
print("""
DESeq2:
  ✓ Gold standard
  ✓ Shrinkage
  ✓ Well-documented
  ✗ Slower

edgeR:
  ✓ Fast
  ✓ Flexible
  ✓ Small samples
  ✗ Less shrinkage

limma-voom:
  ✓ Fastest
  ✓ Large datasets
  ✓ Linear model
  ✗ Less robust for small n

Recommendation:
  → DESeq2 for most
  → edgeR for small n
  → limma for large n
""")

print("\n=== Speed Comparison ===")
print(f"{'Tool':<20} {'10M reads':<15} {'50M reads':<15} {'100M reads':<15}")
print("-" * 65)
print(f"{'STAR':<20} {'2 min':<15} {'10 min':<15} {'20 min':<15}")
print(f"{'HISAT2':<20} {'5 min':<15} {'25 min':<15} {'50 min':<15}")
print(f"{'Salmon':<20} {'1 min':<15} {'5 min':<15} {'10 min':<15}")
print(f"{'kallisto':<20} {'30 sec':<15} {'2 min':<15} {'5 min':<15}")
print(f"{'featureCounts':<20} {'30 sec':<15} {'2 min':<15} {'5 min':<15}")

print("\n=== When to Use Which ===")
print("""
USE STAR:
  ✓ General purpose
  ✓ Novel junctions
  ✓ Fusion detection

USE HISAT2:
  ✓ Limited memory
  ✓ Large genomes

USE SALMON:
  ✓ Fast quantification
  ✓ No alignment
  ✓ Many samples

USE DESeq2:
  ✓ Gold standard
  ✓ Small samples
  ✓ Shrinkage

USE edgeR:
  ✓ Small samples
  ✓ Complex designs

USE limma-voom:
  ✓ Large datasets
  ✓ Fast
""")

💼 Use Cases

🎯 Where RNA-Seq Is Used

RNA-Seq is used in virtually every area of biology: disease, development, cancer, drug discovery, and more.

Use Cases

Disease: vs normal. Cancer: Subtypes. Drug: Response. Development: Time course. Single-cell: Heterogeneity. Clinical: Diagnostics.

💼 Disease💼 Cancer💼 Drug💼 Development

# Use Cases
print("=== 1. Disease vs Normal ===")
print("""
Use case:
  - Identify disease genes
  - Understand mechanism
  - Biomarkers

Workflow:
  1. Collect samples
  2. RNA-Seq
  3. DE analysis
  4. Pathway
  5. Validate

Examples:
  - Cancer vs normal
  - Alzheimer's
  - Diabetes
  - Autoimmune

Tools:
  - DESeq2
  - clusterProfiler
""")

print("\n=== 2. Cancer Subtypes ===")
print("""
Use case:
  - Classify tumors
  - Prognosis
  - Treatment

Workflow:
  1. Tumor samples
  2. RNA-Seq
  3. Clustering
  4. Subtypes
  5. Survival

Examples:
  - Breast cancer
  - Leukemia
  - Glioblastoma

Tools:
  - DESeq2
  - ConsensusClusterPlus
""")

print("\n=== 3. Drug Response ===")
print("""
Use case:
  - Mechanism of action
  - Resistance
  - Biomarkers

Workflow:
  1. Treat cells
  2. RNA-Seq
  3. DE
  4. Pathway
  5. Validate

Examples:
  - Kinase inhibitors
  - Chemotherapy
  - Immunotherapy

Tools:
  - DESeq2
  - GSEA
""")

print("\n=== 4. Development ===")
print("""
Use case:
  - Time course
  - Differentiation
  - Organogenesis

Workflow:
  1. Time points
  2. RNA-Seq
  3. DE over time
  4. Trajectory

Examples:
  - Embryonic
  - Stem cells
  - Organoids

Tools:
  - DESeq2
  - Monocle
""")

print("\n=== 5. Single-Cell ===")
print("""
Use case:
  - Heterogeneity
  - Cell types
  - Rare cells
  - Trajectories

Workflow:
  1. Cells
  2. 10x
  3. Cell Ranger
  4. Seurat/Scanpy
  5. Clusters
  6. Markers

Examples:
  - Immune atlas
  - Brain atlas
  - Tumor microenvironment

Tools:
  - Seurat
  - Scanpy
""")

print("\n=== 6. Clinical Diagnostics ===")
print("""
Use case:
  - Fusion detection
  - Expression signatures
  - Diagnosis

Workflow:
  1. Patient sample
  2. RNA-Seq
  3. Fusion calling
  4. Report

Examples:
  - Leukemia
  - Sarcoma
  - Lung cancer

Tools:
  - STAR-Fusion
  - Arriba
""")

print("\n=== 7. Agriculture ===")
print("""
Use case:
  - Crop traits
  - Stress response
  - Breeding

Workflow:
  1. Conditions
  2. RNA-Seq
  3. DE
  4. Pathways

Examples:
  - Drought
  - Disease
  - Yield

Tools:
  - DESeq2
  - clusterProfiler
""")

print("\n=== 8. Microbiology ===")
print("""
Use case:
  - Gene regulation
  - Pathogenesis
  - Antibiotic

Workflow:
  1. Conditions
  2. RNA-Seq
  3. DE
  4. Operons

Examples:
  - E. coli
  - S. aureus
  - M. tuberculosis

Tools:
  - DESeq2
  - Rockhopper
""")

print("\n=== Comparison ===")
print(f"{'Use Case':<25} {'Tool':<20} {'Output':<20}")
print("-" * 65)
print(f"{'Disease':<25} {'DESeq2':<20} {'DE genes':<20}")
print(f"{'Cancer':<25} {'DESeq2 + clustering':<20} {'Subtypes':<20}")
print(f"{'Drug':<25} {'GSEA':<20} {'Pathways':<20}")
print(f"{'Development':<25} {'Monocle':<20} {'Trajectory':<20}")
print(f"{'Single-cell':<25} {'Seurat':<20} {'Clusters':<20}")
print(f"{'Clinical':<25} {'STAR-Fusion':<20} {'Fusions':<20}")
print(f"{'Agriculture':<25} {'DESeq2':<20} {'DE genes':<20}")
print(f"{'Microbiology':<25} {'Rockhopper':<20} {'Operons':<20}")

❓ FAQ

🎯 Frequently Asked Questions

Answers to common questions about RNA-Seq analysis.

What is RNA-Seq?

What is the difference between alignment-based and alignment-free quantification?

Alignment-based methods (STAR, HISAT2 + featureCounts) align reads to a reference genome before counting. Alignment-free methods (Salmon, kallisto) use k-mer matching to a transcriptome, which is much faster and memory-efficient, with comparable accuracy for differential expression.

Which differential expression tool should I use?

Use DESeq2 for most bulk RNA-Seq analyses (robust, well-documented), edgeR for small sample sizes or complex designs, and limma-voom for large datasets or when speed is critical. All three are R/Bioconductor packages and produce similar results for typical experiments.

How many biological replicates do I need for RNA-Seq?

A minimum of 3 biological replicates per condition is strongly recommended. More replicates increase statistical power and reduce false positives. For cell lines, 3-4 is common; for human samples, 10+ is ideal. Technical replicates are not a substitute for biological replicates.

What is the difference between bulk and single-cell RNA-Seq?

Bulk RNA-Seq measures average expression across all cells in a sample. Single-cell RNA-Seq (scRNA-Seq) measures expression in individual cells, revealing cell-type heterogeneity, rare cell populations, and developmental trajectories. scRNA-Seq requires specialized tools like Seurat, Scanpy, and Cell Ranger.

What is TPM and FPKM?

TPM (Transcripts Per Million) and FPKM (Fragments Per Kilobase Million) are normalization methods. TPM is preferred because it's comparable across samples. FPKM is deprecated. For DE analysis, use raw counts with DESeq2/edgeR, not TPM/FPKM.

What is the difference between stranded and unstranded?

Stranded libraries preserve information about which DNA strand the RNA came from. This improves quantification accuracy, especially for overlapping genes and antisense transcripts. Most modern protocols are stranded. Use -s 2 for reverse-stranded, -s 1 for forward-stranded.

How do I handle batch effects?

Include batch in your design formula (~ batch + condition) if it's not confounded with condition. For visualization, use limma::removeBatchEffect or ComBat-seq. Never correct for batch if it's confounded with your biological variable.

What is GSEA?

GSEA (Gene Set Enrichment Analysis) tests whether a predefined gene set is enriched at the top or bottom of a ranked gene list. Unlike ORA, it uses all genes (no cutoff) and is more sensitive. It's implemented in clusterProfiler, fgsea, and the GSEA desktop app.

What is the difference between DESeq2 and edgeR?

Both use negative binomial models. DESeq2 uses median-of-ratios normalization and shrinkage estimation. edgeR uses TMM normalization and offers quasi-likelihood F-tests. They produce highly correlated results. DESeq2 is more commonly used; edgeR is preferred for small sample sizes and complex designs.

# Quick Reference
print("=== Quick Answers ===")
print("""
Q: What is RNA-Seq?
A: Sequence RNA to measure expression.

Q: Best aligner?
A: STAR (accurate), HISAT2 (memory).

Q: Best quantifier?
A: featureCounts (alignment), Salmon (fast).

Q: Best DE tool?
A: DESeq2 (gold standard).

Q: How many replicates?
A: ≥ 3 biological.

Q: TPM vs FPKM?
A: TPM preferred.

Q: Stranded?
A: Yes, use -s 2 for reverse.

Q: Batch effects?
A: Include in design.

Q: Single-cell tools?
A: Seurat (R), Scanpy (Python).

Q: Pathway analysis?
A: GSEA, clusterProfiler.

Q: Free tools?
A: Yes — all major tools are free.
""")

🚀 Future

🎯 Where RNA-Seq Is Going

The future is single-cell, spatial, long-read, AI-augmented, and clinical-grade.

Future Directions

Single-cell: Atlases. Spatial: Location. Long-read: Isoforms. AI: Foundation models. Clinical: Routine. Multi-omics: Integration.

🚀 Single-cell🚀 Spatial🚀 Long-read🚀 AI

# Future of RNA-Seq
print("=== 1. Single-Cell Atlases ===")
print("""
Trend: Cell atlases

Current:
  - Thousands of cells
  - Few tissues

Future:
  - Millions of cells
  - All tissues
  - All species
  - Human Cell Atlas

Examples:
  - Human Cell Atlas
  - Tabula Sapiens
  - Mouse Cell Atlas
  - Fly Cell Atlas

Tools:
  - Seurat
  - Scanpy
  - CellxGene
""")

print("\n=== 2. Spatial Transcriptomics ===")
print("""
Trend: Location + expression

Current:
  - Dissociated cells
  - Lose location

Spatial:
  - Tissue context
  - Cell-cell interactions
  - Tumor microenvironment

Platforms:
  - 10x Visium
  - MERFISH
  - seqFISH
  - Slide-seq
  - Stereo-seq

Tools:
  - Seurat
  - Squidpy
  - Giotto
""")

print("\n=== 3. Long-Read RNA-Seq ===")
print("""
Trend: Full isoforms

Current:
  - Short reads
  - Isoform inference

Long reads:
  - Full-length isoforms
  - Alternative splicing
  - Fusion genes
  - Novel transcripts

Platforms:
  - PacBio Iso-Seq
  - Nanopore

Tools:
  - IsoQuant
  - StringTie2
  - FLAIR
  - Bambu
""")

print("\n=== 4. AI Foundation Models ===")
print("""
Trend: Deep learning

Ideas:
  - Gene expression prediction
  - Cell type annotation
  - Perturbation prediction
  - Drug response

Examples:
  - scBERT
  - scGPT
  - Geneformer
  - scFoundation

Future:
  - Learned representations
  - Transfer learning
  - Zero-shot
  - Multi-modal
""")

print("\n=== 5. Clinical RNA-Seq ===")
print("""
Trend: Routine clinical

Current:
  - Research
  - Specialized

Future:
  - Clinical routine
  - Fusion detection
  - Expression signatures
  - Drug response

Examples:
  - Oncotype DX
  - MammaPrint
  - FoundationOne
  - Caris

Requirements:
  - Accuracy
  - Speed
  - Cost
  - Regulation
""")

print("\n=== 6. Multi-Omics Integration ===")
print("""
Trend: Integrate everything

Current:
  - Single modality
  - Separate analysis

Future:
  - Multi-omics
  - Integration
  - Systems biology

Examples:
  - Genomics + transcriptomics
  - Proteomics + transcriptomics
  - Metabolomics + transcriptomics

Tools:
  - MOFA
  - mixOmics
  - Seurat (multi-modal)
""")

# TIMELINE
print("\n=== Timeline ===")
print("""
Past:
  2006: First RNA-Seq
  2010: TopHat
  2012: STAR
  2013: DESeq2
  2014: Salmon
  2015: 10x Genomics
  2020: Single-cell atlases
  2022: Spatial transcriptomics

Present:
  2025: Multi-omics
  2026: AI foundation models

Future:
  2030: Clinical routine?
  2035: Whole-organism atlases?
  2040: Predictive biology?
""")

print("\n=== Final Thoughts ===")
print("""
RNA-Seq is:
  - Fundamental to transcriptomics
  - Needed for disease research
  - Needed for drug discovery
  - Rapidly improving

The future is:
  - Single-cell
  - Spatial
  - Long-read
  - AI-augmented
  - Clinical-grade
  - Multi-omics

"RNA-Seq reveals the music of the genome.
 Analysis is how we hear it."

Keep exploring. Keep building.
""")