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

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.

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.

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)

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

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):

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):

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