deseq2-differential-expression

DESeq2 による差次的発現解析（R/Bioconductor）

概要

DESeq2 は、バルク RNA-seq のカウントデータを用いた差次的遺伝子発現解析のための Bioconductor R パッケージです。遺伝子ごとに負の二項一般化線形モデルを当てはめ、遺伝子間の経験ベイズ縮小を用いて分散パラメーターを推定し、Wald 検定（二群の場合）または尤度比検定（複雑なデザインの場合）によって差次的発現を検定します。DESeq2 は RNA-seq の差次的発現解析における R のゴールドスタンダードであり、Salmon（tximeta／tximport）、featureCounts、HTSeq からシームレスにデータをインポートできる Bioconductor とのネイティブな連携機能を備えています。

使用する場面

• バルク RNA-seq のカウントデータから、2 つの実験条件（処理群と対照群、疾患群と健常群など）間で差次的に発現する遺伝子を特定する場合
• バッチ効果や共変量を考慮した多因子デザイン（例：~ batch + condition）を解析する場合
• 交互作用項（例：時間 × 処理）を含む複雑な仮説や、尤度比検定（LRT）を用いた縮小モデルを検定する場合
• tximeta または tximport を介して Salmon の疑似アラインメント出力をインポートし、転写産物レベルの不確実性を伝播する場合
• apeglm を用いて LFC を縮小し、遺伝子のランキングリスト、ボルケーノプロット、下流のパスウェイ解析に利用する場合
• 時系列実験や、モデル比較が必要な 2 水準を超えるデザインを実施する場合
• SummarizedExperiment、clusterProfiler、EnhancedVolcano との連携が必要な R／Bioconductor エコシステムで作業する場合
• Python ベースのパイプラインで同じ統計モデルを使用する場合は、代わりに pydeseq2-differential-expression を使用してください
• TMM 正規化、準尤度 F 検定、または TREAT 検定による負の二項差次的発現解析には edgeR を使用してください
• 差次的発現解析の後、DESeq2 の結果を論文品質のプロットにするには omics-plotting SKILL を使用してください
• 差次的発現解析の後、パスウェイレベルで結果を解釈するには gseapy-gene-enrichment SKILL を使用してください

前提条件

• R パッケージ：DESeq2（Bioconductor）、tximeta または tximport（Salmon のインポート）、apeglm（LFC の縮小）、pheatmap、ggplot2、EnhancedVolcano
• データ要件：生の（正規化されていない）整数カウント行列（遺伝子を行 × サンプルを列）と、列名が一致するサンプルメタデータのデータフレーム。Salmon を使用する場合は、サンプルごとの quant.sf ファイルと転写産物から遺伝子への対応表が必要です
• 環境：R ≥ 4.2、Bioconductor ≥ 3.16

# Install Bioconductor packages
if (!require("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install(c("DESeq2", "tximeta", "tximport", "apeglm",
                        "EnhancedVolcano"))
install.packages(c("pheatmap", "ggplot2", "dplyr"))

クイックスタート

カウント行列から 2 群の比較を 20 行未満で実行できます。

library(DESeq2)

# Load count matrix (genes x samples) and metadata
counts <- as.matrix(read.csv("counts.csv", row.names = 1))
coldata <- read.csv("metadata.csv", row.names = 1)
coldata$condition <- factor(coldata$condition)

# Build DESeqDataSet, run full pipeline
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData   = coldata,
                               design    = ~ condition)
dds <- dds[rowSums(counts(dds)) >= 10, ]   # pre-filter
dds <- DESeq(dds)

# Extract results (treated vs. control)
res <- results(dds, contrast = c("condition", "treated", "control"),
               alpha = 0.05)
summary(res)
# LFC shrinkage for visualization
res_shrunk <- lfcShrink(dds, contrast = c("condition", "treated", "control"),
                         type = "apeglm", coef = 2)

ワークフロー

ステップ 1：カウント行列とサンプルメタデータを準備する

遺伝子 × サンプルのカウント行列と colData データフレームから DESeqDataSet を作成します。カウント行列の列名は colData の行名と一致している必要があります。

library(DESeq2)
library(dplyr)

# Load raw count matrix (genes as rows, samples as columns)
counts <- as.matrix(read.csv("featureCounts_matrix.csv", row.names = 1))
# Or from featureCounts tab-delimited output — skip first comment line
# fc <- read.table("featurecounts_output.txt", header = TRUE, skip = 1, row.names = 1)
# counts <- as.matrix(fc[, 7:ncol(fc)])  # columns 7+ are sample counts

# Load sample metadata
coldata <- data.frame(
    condition = factor(c("control", "control", "control",
                          "treated", "treated", "treated")),
    batch     = factor(c("A", "A", "B", "A", "B", "B")),
    row.names = colnames(counts)
)

# Build DESeqDataSet
dds <- DESeqDataSetFromMatrix(countData = counts,
                               colData   = coldata,
                               design    = ~ condition)

cat("Samples:", ncol(dds), "\n")
cat("Genes:", nrow(dds), "\n")
cat("Condition levels:", levels(coldata$condition), "\n")

ステップ 2：tximeta を介して Salmon からインポートする

リードの定量に Salmon を使用した場合は、tximeta を使って転写産物レベルの推定値をインポートします。これにより、転写産物の長さと GC バイアスを考慮した適切なオフセット補正が行われます。

library(tximeta)

# Build a coldata with a "files" column pointing to quant.sf files
quant_dirs <- file.path("salmon_output", coldata$sample_id, "quant.sf")
coldata$files <- quant_dirs
coldata$names <- coldata$sample_id

# Import with tximeta — automatically fetches transcript metadata from Ensembl
se <- tximeta(coldata)           # SummarizedExperiment with transcript-level data

# Summarize to gene level (requires Ensembl or custom txdb)
gse <- summarizeToGene(se)       # gene-level SummarizedExperiment

# Build DESeqDataSet from the SummarizedExperiment
dds <- DESeqDataSet(gse, design = ~ condition)

cat("Gene-level SE dimensions:", dim(gse), "\n")
cat("Assay names:", assayNames(gse), "\n")
# Assay "counts" = estimated counts; "abundance" = TPM; "length" = eff. lengths

ステップ 3：事前フィルタリングと品質管理

統計的検出力を高め、多重検定の負担を軽減するため、カウント数が非常に少ない遺伝子を除外します。分散安定化カウントに対する PCA でサンプル品質を確認します。

library(ggplot2)

# Pre-filter: keep genes with at least 10 reads total
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep, ]
cat("Genes after pre-filtering:", nrow(dds), "\n")

# Variance-stabilizing transformation for QC visualization
vsd <- vst(dds, blind = TRUE)    # blind = TRUE for exploratory QC

# PCA plot
pca_data <- plotPCA(vsd, intgroup = c("condition", "batch"), returnData = TRUE)
percent_var <- round(100 * attr(pca_data, "percentVar"))

ggplot(pca_data, aes(PC1, PC2, color = condition, shape = batch)) +
    geom_point(size = 3) +
    xlab(paste0("PC1: ", percent_var[1], "% variance")) +
    ylab(paste0("PC2: ", percent_var[2], "% variance")) +
    ggtitle("PCA — Variance-Stabilized Counts") +
    theme_bw()

ggsave("pca_plot.pdf", width = 6, height = 5)
cat("Saved pca_plot.pdf\n")

ステップ 4：DESeq() を実行する — 正規化、分散推定、モデル当てはめ

DESeq() は 3 つのステップを順に実行します。（1）中央値比によるサイズファクターの推定、（2）遺伝子ごとの分散と縮小分散の推定、（3）負の二項 GLM の当てはめと Wald 統計量の計算。

# Run full DESeq2 pipeline
dds <- DESeq(dds)

# Inspect size factors (library normalization)
cat("Size factors:\n")
print(sizeFactors(dds))
# Expected range: 0.3–3.0; outliers indicate library quality issues

# Dispersion plot — should show decreasing dispersion trend with mean
plotDispEsts(dds, main = "Dispersion Estimates")

# Inspect fitted model coefficients
resultsNames(dds)
# [1] "Intercept" "condition_treated_vs_control"

ステップ 5：結果を抽出して FDR 補正を適用する

results() を使って Wald 検定の統計量を抽出します。明確さのため、コントラストを明示的に指定してください。独立フィルタリングにより、検出可能な遺伝子の数を最大化できます。

# Extract results for the contrast of interest
res <- results(dds,
               contrast          = c("condition", "treated", "control"),
               alpha             = 0.05,   # FDR threshold for independent filtering
               lfcThreshold      = 0,      # Test H0: |LFC| = 0 (Wald test)
               independentFilter = TRUE)   # Maximize power via adaptive filtering

# Summary: how many genes are up/down/NA
summary(res)
# out of N with nonzero total read count:
# adjusted p-value < 0.05: ...
# LFC > 0 (up): ...
# LFC < 0 (down): ...
# outliers [Cook's distance]: ...
# low counts [below independent filtering threshold]: ...

# Convert to data frame and sort by adjusted p-value
res_df <- as.data.frame(res)
res_df <- res_df[order(res_df$padj, na.last = TRUE), ]

cat("Significant genes (padj < 0.05):", sum(res_df$padj < 0.05, na.rm = TRUE), "\n")
cat("Upregulated:", sum(res_df$padj < 0.05 & res_df$log2FoldChange > 0, na.rm = TRUE), "\n")
cat("Downregulated:", sum(res_df$padj < 0.05 & res_df$log2FoldChange < 0, na.rm = TRUE), "\n")

write.csv(res_df, "deseq2_all_results.csv")

ステップ 6：apeglm による LFC の縮小

適応型 t 事前分布（apeglm）を用いて、log2 フォールドチェンジを縮小します。MA プロット、ボルケーノプロット、遺伝子ランキングには縮小後の LFC を使用してください。有意性の判定には使用しないでください（padj は縮小前の Wald 検定から得られます）。

library(apeglm)

# apeglm shrinkage — specify coefficient name from resultsNames(dds)
res_shrunk <- lfcShrink(dds,
                         coef = "condition_treated_vs_control",
                         type = "apeglm",
                         res  = res)    # pass original results for speed

# Compare shrunk vs. unshrunk LFC
cat("Max |LFC| unshrunk:", max(abs(res$log2FoldChange), na.rm = TRUE), "\n")
cat("Max |LFC| shrunk:  ", max(abs(res_shrunk$log2FoldChange), na.rm = TRUE), "\n")

# Export shrunk results
res_shrunk_df <- as.data.frame(res_shrunk)
res_shrunk_df <- res_shrunk_df[order(res_shrunk_df$padj, na.last = TRUE), ]
write.csv(res_shrunk_df, "deseq2_shrunk_results.csv")

# Filter significant with effect size threshold
sig_genes <- res_shrunk_df[
    !is.na(res_shrunk_df$padj) &
    res_shrunk_df$padj < 0.05 &
    abs(res_shrunk_df$log2FoldChange) > 1, ]
cat("Significant (padj<0.05, |LFC|>1):", nrow(sig_genes), "\n")
write.csv(sig_genes, "deseq2_significant.csv")

ステップ 7：可視化 — MA プロット、ボルケーノプロット、ヒートマップ

計算と可視化の分離（必須）： DESeq2 による計算とプロットを同じコードブロックにまとめないでください。ステップ 6（またはステップ 5）で結果を CSV に保存し、このステップではその CSV を読み込んでプロットのみを作成します。こうすることで、負荷の高い計算を再実行せずに、プロットを再生成したりスタイルを変更したりできます。 プロット作成： ボルケーノプロットでは、エクスポートした DEG CSV に対して skills/data-visualization/omics-plotting/SKILL.md を読み、「Volcano」のレシピに従ってください（→ figures/volcano_plot.pdf）。以下のヒートマップは、引き続きネイティブ R（pheatmap）で作成します。

# Volcano: use the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) "Volcano" recipe on the exported DEG CSV -> figures/volcano_plot.pdf

library(pheatmap)
# ---Heatmap of top 50 DE genes ---
top50 <- head(rownames(sig_genes[order(sig_genes$padj), ]), 50)
mat   <- assay(vsd)[top50, ]            # variance-stabilized expression
mat   <- mat - rowMeans(mat)            # center by gene mean

annotation_col <- data.frame(
    Condition = coldata$condition,
    Batch     = coldata$batch,
    row.names = colnames(mat)
)
# Use the code in the omics-plotting SKILL (`skills/data-visualization/omics-plotting/SKILL.md`) for publication-quality heatmap plots or use below code to generate a heatmap
pheatmap(mat,
         annotation_col  = annotation_col,
         cluster_rows    = TRUE,
         cluster_cols    = TRUE,
         show_rownames   = TRUE,
         fontsize_row    = 6,
         filename        = "figures/heatmap_top50.pdf",
         width           = 8,
         height          = 10)

ステップ 8：多因子デザインと交互作用項

バッチ補正、ペアサンプル、または 2 つの変数間の交互作用効果を扱うデザイン向けです。

# --- Batch-corrected two-group comparison ---
dds_batch <- DESeqDataSetFromMatrix(countData = counts,
                                     colData   = coldata,
                                     design    = ~ batch + condition)
dds_batch <- dds_batch[rowSums(counts(dds_batch)) >= 10, ]
dds_batch <- DESeq(dds_batch)

res_batch <- results(dds_batch,
                      contrast = c("condition", "treated", "control"),
                      alpha    = 0.05)
cat("Significant (batch-corrected):", sum(res_batch$padj < 0.05, na.rm = TRUE), "\n")

# --- Interaction design: test whether treatment effect differs between genotypes ---
# H0: the effect of treatment is the same in WT and KO
coldata$genotype <- factor(coldata$genotype)  # "WT" or "KO"

dds_int <- DESeqDataSetFromMatrix(
    countData = counts, colData = coldata,
    design = ~ genotype + condition + genotype:condition
)
dds_int <- DESeq(dds_int)
resultsNames(dds_int)
# Interaction coefficient: "genotypeKO.conditiontreated"

# Extract interaction term — genes where treatment effect differs by genotype
res_int <- results(dds_int, name = "genotypeKO.conditiontreated", alpha = 0.05)
cat("Genes with significant interaction:", sum(res_int$padj < 0.05, na.rm = TRUE), "\n")

# --- Likelihood Ratio Test for complex designs (e.g., time series) ---
# LRT compares full model vs. reduced model
dds_lrt <- DESeq(dds_int, test = "LRT", reduced = ~ genotype + condition)
res_lrt  <- results(dds_lrt)
cat("Genes with genotype:condition interaction (LRT):",
    sum(res_lrt$padj < 0.05, na.rm = TRUE), "\n")

主要なパラメーター

パラメーター	関数	デフォルト	範囲／選択肢	効果
design	DESeqDataSetFromMatrix()	（必須）	~ var または ~ cov + var	GLM の式を指定します。対象の変数より前にバッチ／共変量を配置します
contrast	results()	NULL	c(var, level1, level2) または係数名	比較を定義します。明確にするため、常に明示的に指定してください
alpha	results()	0.1	0.01–0.10	独立フィルタリングに使用する FDR 閾値です。padj の値は変わりません
lfcThreshold	results()	0	0–2	H0：
independentFilter	results()	TRUE	TRUE／FALSE	棄却可能な仮説の数を最大化するための適応的な平均カウントフィルタリングです
type	lfcShrink()	"apeglm"	"apeglm"、"ashr"、"normal"	縮小事前分布です。apeglm は最速で、最も適切に校正されています
coef	lfcShrink()	（必須）	resultsNames() の係数名	縮小するモデル係数です。resultsNames() の出力と一致する必要があります
test	DESeq()	"Wald"	"Wald"、"LRT"	Wald：2 群間の比較。LRT：複雑なデザインのモデル比較
reduced	DESeq()	NULL	検定対象の項を含まない式	LRT に必須です。帰無モデルを指定します
blind	vst()／rlog()	TRUE	TRUE／FALSE	QC／探索には TRUE、分散情報を使用する後続解析には FALSE を指定します

主要な概念

負の二項モデルと分散

DESeq2 は、各遺伝子のカウントを、平均 = µ、分散パラメーター = α の負の二項分布としてモデル化します。平均 µ は、対数線形モデルを介して実験デザインに依存します。重要な革新は分散の縮小です。経験ベイズ法を用いて、遺伝子ごとの分散推定値を全遺伝子にわたる適合トレンドに向けて縮小し、遺伝子間の情報を活用して、サンプル数が少ない場合の推定値を安定化します。

サイズファクターによる正規化

DESeq2 は比の中央値法を使用します。各サンプルについて、各遺伝子のカウントを全サンプルにわたるその遺伝子の幾何平均で割った比を計算し、すべての遺伝子にわたる比の中央値をサイズファクターとします。この方法は差次的発現遺伝子の影響を受けにくく、差次的発現遺伝子の割合に関する仮定も必要ありません。

Wald 検定と尤度比検定

• Wald 検定（デフォルト）：単一の係数 β = 0 を検定します。高速で、2 群間の比較や多因子デザインにおける特定のコントラストに適しています
• LRT：尤度比を用いて完全モデルと reduced モデルを比較します。多水準因子（例：5 つの時点）のいずれかの水準に効果があるかを検定する場合や、交互作用項をまとめて検定する場合に適しています

Cook 距離による外れ値検出

DESeq2 は、発現の外れ値の可能性を検出するため、サンプルごと、遺伝子ごとに Cook 距離を計算します。いずれかのサンプルで Cook 距離が閾値を超える遺伝子は、結果で自動的に NA に置き換えられます。refit = TRUE（デフォルト）の場合、外れ値サンプルを除外し、その遺伝子についてモデルを再適合します。

よく使われるレシピ

レシピ：1つのフィット済みモデルから複数のコントラストを取得する

複数の処理群を効率よく比較できます。DESeq()を一度だけ実行し、複数のコントラストを抽出します。

# Fit once
dds_multi <- DESeqDataSetFromMatrix(counts, coldata, design = ~ condition)
dds_multi <- dds_multi[rowSums(counts(dds_multi)) >= 10, ]
dds_multi <- DESeq(dds_multi)

# Extract contrasts — all vs. "control"
contrasts_list <- list(
    A_vs_ctrl = c("condition", "treatment_A", "control"),
    B_vs_ctrl = c("condition", "treatment_B", "control"),
    C_vs_ctrl = c("condition", "treatment_C", "control")
)

results_list <- lapply(names(contrasts_list), function(nm) {
    res <- results(dds_multi, contrast = contrasts_list[[nm]], alpha = 0.05)
    res_shrunk <- lfcShrink(dds_multi, contrast = contrasts_list[[nm]],
                             type = "apeglm", res = res)
    write.csv(as.data.frame(res_shrunk), paste0("results_", nm, ".csv"))
    cat(nm, "— significant:", sum(res$padj < 0.05, na.rm = TRUE), "\n")
    as.data.frame(res_shrunk)
})
names(results_list) <- names(contrasts_list)

レシピ：ペアサンプルのデザイン

被験者をブロッキング因子として含めることで、被験者内のばらつきを取り除きます。

# Paired design: each subject has pre- and post-treatment samples
coldata$subject   <- factor(c("S1","S1","S2","S2","S3","S3"))
coldata$timepoint <- factor(c("pre","post","pre","post","pre","post"))

dds_paired <- DESeqDataSetFromMatrix(counts, coldata,
                                      design = ~ subject + timepoint)
dds_paired <- dds_paired[rowSums(counts(dds_paired)) >= 10, ]
dds_paired <- DESeq(dds_paired)

res_paired <- results(dds_paired,
                       contrast = c("timepoint", "post", "pre"),
                       alpha    = 0.05)
cat("Significant (paired):", sum(res_paired$padj < 0.05, na.rm = TRUE), "\n")

レシピ：縮小LFC × -log10(pvalue)を使ってGSEA用の遺伝子を順位付けする

GSEAまたはgseapyのprerank関数に入力する、事前順位付け済みの遺伝子リストを作成します。

# Build ranking metric: signed -log10(pvalue) with LFC direction
res_rank <- as.data.frame(res_shrunk)
res_rank <- res_rank[!is.na(res_rank$pvalue) & !is.na(res_rank$log2FoldChange), ]
res_rank$rank_metric <- sign(res_rank$log2FoldChange) *
                         (-log10(res_rank$pvalue + 1e-300))
res_rank <- res_rank[order(res_rank$rank_metric, decreasing = TRUE), ]

# Export as two-column TSv for fgsea or gseapy prerank
write.table(data.frame(gene  = rownames(res_rank),
                        score = res_rank$rank_metric),
            "gsea_ranked_list.rnk",
            sep = "\t", quote = FALSE, row.names = FALSE, col.names = FALSE)
cat("Ranked gene list saved:", nrow(res_rank), "genes\n")

レシピ：フィット済みDESeq2オブジェクトを保存して再読み込みする

後続のセッションで、時間のかかるDESeq()の実行を避けます。

# Save fitted DDS object
saveRDS(dds, "dds_fitted.rds")
cat("Saved fitted DESeqDataSet to dds_fitted.rds\n")

# Reload in a new session
dds_loaded <- readRDS("dds_fitted.rds")
# Continue with results(), lfcShrink(), etc.
res_reloaded <- results(dds_loaded, contrast = c("condition", "treated", "control"))

期待される出力

ファイル	説明
deseq2_all_results.csv	全遺伝子の完全な結果テーブル：baseMean、log2FoldChange、lfcSE、stat、pvalue、padj
deseq2_shrunk_results.csv	apeglmでLFCを縮小した結果。順位付けと可視化に使用します
deseq2_significant.csv	フィルタリング済みの結果：padj < 0.05かつ
pca_plot.pdf	分散安定化したカウントのPCA。バッチ構造や外れ値の確認に使用します
dispersion_plot.pdf	遺伝子ごとの分散とフィットされた分散。トレンド周辺に密集した分布が見られるはずです
figures/volcano_plot.pdf	有意性とLFCのしきい値を表示したボルケーノプロット
figures/heatmap_top50.pdf	上位50個のDE遺伝子を階層クラスタリングしたヒートマップ
gsea_ranked_list.rnk	パスウェイ濃縮解析用の事前順位付け済み遺伝子リスト（fgsea/gseapy）
dds_fitted.rds	チェックポイント／再開用にシリアライズしたDESeqDataSet

トラブルシューティング

問題	原因	解決策
Error: design contains one or more variables with all samples having the same value	デザイン因子のレベルが1つしかありません	table(coldata$condition)を確認し、すべての因子レベルが存在することを確認します
Model matrix not full rank	デザインが交絡しています（例：バッチと条件が完全に相関している）	table(coldata$batch, coldata$condition)を実行し、交絡している変数を除外するか、サンプルを追加します
ほとんど、またはすべてのpadj = NA	独立フィルタリング（平均カウントが低い）またはCookの外れ値によって遺伝子がフラグ付けされています	summary(res)でフィルタリングされた遺伝子数を確認し、alphaを緩和するか、データ品質を確認します
生物学的な効果が強いにもかかわらず、有意な遺伝子がありません	実験の検出力不足、またはグループ内のばらつきが大きい	グループあたりn ≥ 3であることを確認し、PCAで外れ値を確認します。p値のヒストグラムを調べます（0付近にピークがあるはずです）
Error in lfcShrink: 係数が見つかりません	coef名がresultsNames(dds)と一致していません	resultsNames(dds)を実行し、係数の文字列を正確にコピーしてcoef =に指定します
apeglmが収束警告で失敗します	LFC推定値が極端です（データが疎、または完全分離）	このようなケースにより頑健なtype = "ashr"に切り替えます
サイズファクターが非常に大きい（> 5）	ライブラリサイズが極端に異なるか、正規化済みカウントを誤って使用しています	colSums(counts(dds))を確認し、入力がアライナー／カウンターの生の整数カウントであることを確認します
分散プロットで、トレンドから大きく離れた外れ値の集まりが見られます	ライブラリのノイズまたは失敗、あるいはメタデータのグループ分けが不適切です	PCAを確認し、外れ値のサンプルを除外します。design式が生物学的グループと一致していることを確認します

参考資料

• DESeq2 Bioconductorページ — パッケージのホーム、ビネット、リファレンスマニュアル
• DESeq2ビネット — すべてのデザインタイプを網羅した包括的な実例
• Love et al. (2014) Genome Biology — DESeq2の原著論文：RNA-seqデータのフォールドチェンジと分散のモデレート推定
• Zhu et al. (2019) Bioinformatics — DESeq2向けのapeglmによるLFC縮小
• tximeta Bioconductorページ — 転写産物レベルのメタデータを含むSalmonのインポート
• EnhancedVolcano Bioconductorページ — DE結果用の出版品質のボルケーノプロット