為什麼 raw counts 不能直接比？
RNA-seq 與 LFQ proteomics 的原始 counts 受 (1) 樣本總定序深度、(2) 基因/蛋白長度、(3) 組成偏差影響。Wagner et al. 2012 證明 RPKM/FPKM 在樣本間不可比（總和不為定值），TPM 才是現代標準。但 TPM 仍不適合 DE 分析——應改用 negative binomial 模型（DESeq2 / edgeR）或 voom 加權 (Law 2014)。
差異基因分析的金標準是 DESeq2 (Love et al. 2014, Genome Biol) + apeglm (Zhu et al. 2018) 的 shrinkage log2FC。蛋白體 LFQ 用 MaxLFQ (Cox et al. 2014, MCP) 配合 Perseus / FragPipe / DIA-NN。
Raw counts from RNA-seq / LFQ proteomics are influenced by (1) total sequencing depth, (2) gene/protein length, (3) composition bias. Wagner 2012 proved RPKM/FPKM are not comparable across samples (sum is not constant); TPM is now the standard. But TPM is still wrong for DE — use negative binomial models (DESeq2 / edgeR) or voom weighting (Law 2014).
The DE gold standard: DESeq2 (Love et al. 2014, Genome Biol) + apeglm (Zhu et al. 2018) shrinkage log2FC. Proteomics LFQ uses MaxLFQ (Cox et al. 2014, MCP) with Perseus / FragPipe / DIA-NN.
一、定量單位辨識卡
| 單位 | 計算 | 可比？ | DE 用 |
|---|---|---|---|
| Raw counts | 原始 read 數 | ❌ | ✓ DESeq2 |
| RPKM / FPKM | counts / (length × total/1e6 / 1000) | ❌ 總和不定 | ❌ |
| TPM | 先除長度 | ✓ 總和固定 | ❌ 無變異模型 |
| DESeq2 norm counts | 中位比例 | ✓ | ✓ |
| vst / rlog | VST | ✓ | PCA |
| MaxLFQ intensity | 配對比例 | ✓ | 蛋白 |
apeglm shrinkage 讓低 count 基因不再吵
調整 base mean 與 dispersion，觀察 raw 與 shrunken log2FC 差異。關鍵：低 count 基因 raw log2FC 變異極大，shrinkage 把它們拉回 0，讓排序更可靠。
Adjust base mean and dispersion. Key: low-count genes have wild raw log2FC; shrinkage pulls them toward zero, making rankings reliable.
MA plot
二、DESeq2 完整流程
library(DESeq2); library(apeglm) # 1. 從 tximport 建 dds（建議 gene-level） library(tximport) txi <- tximport(files=salmon_files, type="salmon", tx2gene=tx2gene) dds <- DESeqDataSetFromTximport(txi, colData=coldata, design= ~ condition) # 2. 預過濾低 count（提速） dds <- dds[rowSums(counts(dds)) >= 10, ] # 3. 跑差異分析 dds <- DESeq(dds) res <- results(dds, contrast=c("condition","treated","control")) # 4. log2FC shrinkage（apeglm 推薦） res_shr <- lfcShrink(dds, coef="condition_treated_vs_control", type="apeglm") # 5. VST 用於 PCA / heatmap vsd <- vst(dds, blind=FALSE) plotPCA(vsd, intgroup="condition") # 6. Volcano plot library(EnhancedVolcano) EnhancedVolcano(res_shr, lab=rownames(res_shr), x="log2FoldChange", y="padj")
from pydeseq2.dds import DeseqDataSet from pydeseq2.ds import DeseqStats dds = DeseqDataSet(counts=counts_df, metadata=meta, design_factors="condition", refit_cooks=True) dds.deseq2() stat_res = DeseqStats(dds, contrast=["condition","treated","control"]) stat_res.summary() # 含 padj (BH-FDR) stat_res.lfc_shrink(coeff="condition_treated_vs_control") # apeglm # PCA (sklearn) from sklearn.decomposition import PCA import numpy as np log_counts = np.log1p(dds.layers["normed_counts"]) pca = PCA(n_components=2).fit_transform(log_counts.T)
三、最常見錯誤
❌ 把 TPM 餵 t-test
TPM 不滿足 NB 假設，違反獨立性。請給 DESeq2 raw counts。
TPM violates NB assumptions and independence. Feed raw counts to DESeq2.
❌ 忽略 shrinkage
未 shrink 的 log2FC 在低 count 區噪聲極大；篩 top hits 會優先選到低 count 基因。
Un-shrunken log2FC is noisy at low counts; top hits become dominated by low-count genes.
❌ n=3 就跑 DE
DESeq2 / edgeR 在 n=3 對 dispersion 估計很差，多為 false positive。2025 文獻建議 biological n ≥ 6，理想 12。
DESeq2 / edgeR estimate dispersion poorly at n=3 → high FP. 2025 reviews recommend biological n ≥ 6, ideally 12.
❌ 蛋白體亂 imputation
MaxLFQ 缺失值常為 MNAR（低豐度未檢出）。先以 valid value 過濾再用 left-shifted Gaussian imputation。
MaxLFQ missing values are usually MNAR (low-abundance, not detected). Filter valid values first, then left-shifted Gaussian imputation.
🎯 章末小測驗
1. 樣本間可比的 RNA-seq 單位？
2. 排序 DE 基因用？
3. RNA-seq DE 最少 n？