# Rnaseq Analysis

> Run bulk RNA-seq differential expression from a gene x sample count matrix using DESeq2 (R), pydeseq2, or edgeR — normalization (median-of-ratios/TPM), Wald/LRT testing, BH-adjusted p-values, volcano/MA plots, and GSEA/ORA. Use when doing RNA-seq DE, comparing treatment vs control expression, building a FASTQ-to-DESeq2 pipeline, or asked about TPM/RPKM/FPKM, count matrices, log2FoldChange, padj, or STAR/Salmon/featureCounts/tximport.

- Skill: `pavel-kravchenko/rnaseq-analysis` (Agent Skill)
- Install (CLI): `npx skillmds@latest add pavel-kravchenko/rnaseq-analysis`
- Raw SKILL.md: https://api.skillmd.com/api/skills/pavel-kravchenko/rnaseq-analysis/raw
- Safety review: pending
- Works with: Claude Code, Claude.ai, OpenAI Codex
- Category: DevOps & Infra
- Author: pavel-kravchenko (https://skillmd.com/u/pavel-kravchenko)
- Updated: 2026-09-17
- Page: https://skillmd.com/skills/pavel-kravchenko/rnaseq-analysis

---


# RNA-seq Differential Expression Analysis

## When to Use
- Testing which genes differ between conditions (treatment vs. control, tumor vs. normal) from RNA-seq count data
- Deciding between alignment-based (STAR + featureCounts) and alignment-free (Salmon/kallisto + tximport) quantification
- Choosing/justifying a normalization unit (raw counts vs. TPM/RPKM/FPKM) for a given use case
- Building volcano plots, MA plots, or PCA/clustering QC from a count matrix
- Following up DE results with gene set enrichment (GSEA/ORA)

## Version Compatibility
DESeq2 ≥1.42 (Bioconductor 3.18+), apeglm ≥1.24, edgeR ≥4.0, R ≥4.3, pydeseq2 ≥0.4, Python ≥3.10, pandas ≥2.0, numpy ≥1.26, scipy ≥1.11, scikit-learn ≥1.3, statsmodels ≥0.14.

## Prerequisites
- R: `BiocManager::install(c("DESeq2","apeglm"))`. Python: `pip install pydeseq2 pandas numpy scipy scikit-learn statsmodels matplotlib seaborn`.
- CLI tools (if starting from FASTQ): STAR or HISAT2, Subread (featureCounts), Salmon, and `tximport` (R) for alignment-free counts.
- Concepts: count matrix layout (genes x samples), negative-binomial count model, Benjamini-Hochberg FDR.

## Workflow

```text
FASTQ -> QC (FastQC/fastp) -> Alignment (STAR/HISAT2) -> featureCounts/HTSeq
                                   OR
                          -> Pseudoalignment (Salmon/kallisto) -> tximport
-> Count matrix (genes x samples) -> Normalization -> DESeq2/edgeR -> Volcano/MA -> GSEA/ORA
```

| Unit | Formula | Use case |
|------|---------|----------|
| RPKM/FPKM | `(C / L) / N x 1e9` | Single/paired-end, within-sample only |
| TPM | `(C/L) / sum(Cj/Lj) x 1e6` | Cross-sample comparison (sums to 1M) |
| DESeq2 size factor | `median(count / gene_geomean)` per sample | Robust normalization for DE testing |

C = read count, L = gene length (bp), N = total mapped reads. **Never feed TPM/RPKM/FPKM into DESeq2/edgeR — use raw counts; report TPM only for visualization.**

**Goal:** get a per-sample scale factor that corrects for library size and composition without being skewed by a few very highly expressed genes.
**Approach:** DESeq2's median-of-ratios — divide each gene's count by its across-sample geometric mean, then take the per-sample median of those ratios.

```python
import numpy as np
import pandas as pd

def deseq2_size_factors(count_matrix: pd.DataFrame) -> pd.Series:
    """Compute DESeq2-style median-of-ratios size factors.

    count_matrix: genes (rows) x samples (columns), raw integer counts.
    Only genes with nonzero counts in every sample contribute to the
    geometric mean (required since log(0) is undefined).
    """
    nonzero_mask = (count_matrix > 0).all(axis=1)
    filtered = count_matrix.loc[nonzero_mask]
    geo_means = np.exp(np.log(filtered).mean(axis=1))
    return pd.Series(
        {s: np.median(filtered[s] / geo_means) for s in filtered.columns}
    )

def counts_to_tpm(counts: np.ndarray, lengths: np.ndarray) -> np.ndarray:
    """Convert raw counts to TPM (length-normalize, then scale to 1e6)."""
    rate = counts / lengths
    return rate / rate.sum() * 1e6

size_factors = deseq2_size_factors(counts_df)
normalized_counts = counts_df.div(size_factors, axis=1)
```

**Goal:** call differentially expressed genes with correct FDR control.
**Approach:** for production use DESeq2 (R) or pydeseq2 (Python) — both fit a negative-binomial GLM per gene and borrow strength across genes for dispersion. A quick non-parametric fallback (Mann-Whitney + BH) is fine for prototyping only.

```r
library(DESeq2)
dds <- DESeqDataSetFromMatrix(countData = count_matrix,
                               colData = sample_info, design = ~condition)
dds <- DESeq(dds)                                    # size factors + dispersion + Wald test
res <- results(dds, contrast = c("condition", "Treatment", "Control"))
res_shrunk <- lfcShrink(dds, coef = "condition_Treatment_vs_Control", type = "apeglm")
sig <- subset(res_shrunk, padj < 0.05 & abs(log2FoldChange) > 1)
summary(res)
```

```python
from pydeseq2.dds import DeseqDataSet
from pydeseq2.ds import DeseqStats

dds = DeseqDataSet(counts=count_matrix, metadata=sample_info, design_factors="condition")
dds.deseq2()
stat_res = DeseqStats(dds, contrast=["condition", "Treatment", "Control"])
stat_res.summary()
results_df = stat_res.results_df  # baseMean, log2FoldChange, pvalue, padj
```

```python
from scipy import stats
from statsmodels.stats.multitest import multipletests

def simple_de(count_df: pd.DataFrame, ctrl_samples, treat_samples) -> pd.DataFrame:
    """Prototype-only DE test: BH-corrected Mann-Whitney on DESeq2-normalized counts.
    Underpowered vs. DESeq2/edgeR (no dispersion borrowing) — do not use for publication.
    """
    sf = deseq2_size_factors(count_df)
    norm = count_df.div(sf, axis=1)
    rows = []
    for gene in count_df.index:
        c, t = norm.loc[gene, ctrl_samples], norm.loc[gene, treat_samples]
        lfc = np.log2((t.mean() + 1) / (c.mean() + 1))
        _, p = stats.mannwhitneyu(c, t, alternative="two-sided")
        rows.append({"gene": gene, "log2FC": lfc, "pvalue": p,
                     "baseMean": (c.mean() + t.mean()) / 2})
    df = pd.DataFrame(rows).set_index("gene")
    df["padj"] = multipletests(df["pvalue"], method="fdr_bh")[1]
    return df.sort_values("pvalue")
```

**Goal:** get from raw reads to a count matrix.
**Approach:** align + count (splice-aware, gene-level counts) or pseudoalign + import (fast, transcript-level then aggregated).

STAR + featureCounts (alignment-based):
```bash
STAR --runMode genomeGenerate --genomeDir star_index/ --genomeFastaFiles genome.fa --sjdbGTFfile genes.gtf
STAR --genomeDir star_index/ --readFilesIn R1.fastq R2.fastq --outSAMtype BAM SortedByCoordinate --quantMode GeneCounts
featureCounts -a genes.gtf -o counts.txt -T 4 -p --countReadPairs *.bam
```

Salmon (alignment-free; aggregate transcript counts to genes with `tximport` in R):
```bash
salmon index -t transcriptome.fa -i salmon_index
salmon quant -i salmon_index -l A -1 R1.fastq -2 R2.fastq -o sample_quant --validateMappings
```

## Pitfalls
- **TPM/RPKM/FPKM for DE testing** — ratio-of-ratios artifacts and no variance model; always pass raw counts to DESeq2/edgeR.
- **Library composition bias** — one highly-expressed gene can suppress the apparent expression of everything else; DESeq2 median-of-ratios and edgeR TMM are robust, plain CPM/TPM are not.
- **Skipping dispersion shrinkage** — DESeq2/edgeR borrow information across genes to stabilize low-count variance estimates; a per-gene t-test or Mann-Whitney misses this and is underpowered, especially at n<3/group.
- **Too few replicates** — use >=3 biological replicates per condition (5+ preferred); replicates improve power far more than deeper sequencing.
- **No LFC shrinkage** — raw log2FoldChange is noisy for low-count genes; use `lfcShrink(..., type="apeglm")` before ranking/plotting.
- **Not filtering low counts** — genes with near-zero counts across samples inflate the multiple-testing burden and destabilize dispersion estimates; filter before testing.
- **Confusing biological with technical replicates** — technical replicates (same sample, resequenced) understate true variance and inflate significance.

## See Also
- `bio-differential-expression-deseq2-basics` — DESeq2 API details and design formulas
- `bio-differential-expression-edger-basics` — edgeR/TMM alternative workflow
- `bio-rna-quantification-tximport-workflow` — Salmon/kallisto to gene-level counts
- `bio-pathway-analysis-gsea` — downstream gene set enrichment on DE results