RNA-seq解析入門⑥：発現変動解析編

📚 この記事について
シリーズの締めくくりです。条件間で発現が変化した遺伝子（DEG）を統計的に検定します。
🔙 前の記事：RNA-seq解析入門⑤：カウント編
🔜 次の記事：機能エンリッチメント・パスウェイ解析
📌 前提：カウントマトリクス（counts.txt）が作れていること

RNA-seq（RNA sequencing）の発現変動解析は、条件間で発現量が変化した遺伝子を統計的に検定する工程です。ここでは3つのツールの使い分け、DESeq2 の実践、結果の読み方、そして Python での選択肢を押さえます。

1. 発現変動解析のゴール
2. 3つのツールと使い分け
3. DESeq2 の使い方
4. 結果を読む：火山プロット
5. Python で発現変動解析を行う場合
6. 注意点
まとめ
関連記事
参考文献

1. 発現変動解析のゴール

カウント数を単純に比べるだけでは、「本当に発現が変化したのか、それとも偶然のばらつきか」を判断できません。発現変動解析（DEA: Differential Expression Analysis）は、統計モデルを使ってその差が統計的に有意かどうかを検定し、DEG（発現変動遺伝子）を同定します。

💡 なぜ負の二項分布を使うのか
RNA-seq のカウントデータは、値が 0 に偏り、かつ同じ条件のサンプル間でもばらつきが大きいという性質があります。平均と分散が等しいポアソン分布ではこのばらつきを表現しきれないため、分散を独立に持てる負の二項分布を使います。DESeq2 と edgeR はどちらもこのモデルを採用しています。

2. 3つのツールと使い分け

ツール	統計モデル	特徴
DESeq2	負の二項分布	サンプル数が少なくても安定。ドキュメントが充実
edgeR	負の二項分布	複雑な実験デザイン・多群比較に強い
limma-voom	正規分布（voom 変換後）	サンプル数が多い大規模データに強い

DESeq2 は各群 2〜3 サンプル程度でも安定した結果が得られるよう設計されています（Love et al., Genome Biology, 2014）。正規化から検定までを一貫したワークフローで実行でき、最初に学ぶツールとして適しています。

edgeR も負の二項分布を使いますが、多群比較や交互作用・共変量を含む複雑なデザインへの対応が柔軟です（Robinson et al., Bioinformatics, 2010）。

limma-voom は、マイクロアレイ向けの limma に voom 変換を組み合わせ、RNA-seq のカウントを扱えるようにした方法です（Law et al., Genome Biology, 2014／Ritchie et al., Nucleic Acids Research, 2015）。サンプル数が多いほど検出力を発揮します。

状況	選択
はじめての発現変動解析・各群3サンプル程度	DESeq2
多群比較・交互作用を含む複雑なデザイン	edgeR
サンプル数が多い（10以上）大規模コホート	limma-voom

3. DESeq2 の使い方

インストール

r
# BiocManager 経由でインストール
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

BiocManager::install("DESeq2")

💡 用語：BiocManager
生命科学向けの R パッケージ（Bioconductor）を管理するツールです。DESeq2 や edgeR は CRAN ではなく Bioconductor で配布されているため、BiocManager を使ってインストールします。

ステップ①：カウントマトリクスの読み込み

r
library(DESeq2)

# HTSeq-count が出力した counts.txt を読み込む
count_data <- read.table(
  "counts.txt",
  header = TRUE,
  row.names = 1,
  sep = "\t"
)

# __no_feature などの特殊行を除外する（必須）
count_data <- count_data[!grepl("^__", rownames(count_data)), ]

ステップ②：サンプル情報の設定

r
# どのサンプルが対照群・処理群かを定義する
col_data <- data.frame(
  condition = factor(c("control", "control", "treated", "treated")),
  row.names = colnames(count_data)
)

ステップ③：DESeqDataSet の作成と解析の実行

r
# DESeqDataSet オブジェクトを作る
dds <- DESeqDataSetFromMatrix(
  countData = count_data,
  colData   = col_data,
  design    = ~ condition    # 比較する条件
)

# 低発現の遺伝子を除外（任意だが推奨）
dds <- dds[rowSums(counts(dds)) >= 10, ]

# 正規化と検定を一括で実行
dds <- DESeq(dds)

引数	意味
`countData`	カウントマトリクス（行＝遺伝子、列＝サンプル）
`colData`	サンプルのメタデータ（条件情報）
`design = ~ condition`	比較する条件。バッチを考慮するなら `~ batch + condition` と書きます
`rowSums >= 10`	全サンプル合計が 10 未満の遺伝子を除外し、ノイズを減らします

ステップ④：結果の取得

r
# treated vs control の比較結果を取得
res <- results(dds, contrast = c("condition", "treated", "control"))

# 調整p値でソート
res <- res[order(res$padj), ]

# CSV に保存
write.csv(as.data.frame(res), "DESeq2_results.csv")

# 有意な DEG を絞り込む（padj < 0.05, |log2FC| > 1）
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(sig_genes)  # DEG の数

4. 結果を読む：火山プロット

横軸に log2 fold change（変化の向きと大きさ）、縦軸に -log10(調整p値)（有意性）を取ります。右上が発現の増えた DEG、左上が減った DEG です。中央付近（変化が小さい、または有意でない）は DEG ではありません。

💡 用語：padj と log2FoldChange
padj（調整p値）は、多重検定を補正した p 値（FDR）です。数万の遺伝子を同時に検定すると偶然の当たりが増えるため、補正が必須です。通常 0.05 未満を有意とします。log2FoldChange は発現変化の大きさで、1 なら 2 倍、-1 なら 2 分の1 を意味します。

可視化や品質確認には次の関数も使えます。

lfcShrink()：log2FoldChange の縮小推定。低カウント遺伝子の過大評価を補正します。
vst() / rlog()：分散安定化変換。PCA やヒートマップの前処理に使います。
plotPCA()：サンプル間の類似性を可視化し、外れサンプルを確認します。
plotMA()：発現量と変化量の関係を示す MA プロットを描きます。

5. Python で発現変動解析を行う場合

DESeq2・edgeR・limma-voom はいずれも R のツールです。Python で解析を進めている場合、PyDESeq2 が現時点で最有力の選択肢です。DESeq2 のワークフローを Python で再実装したもので、scverse エコシステムの一部として整備されています（Muzellec et al., Bioinformatics, 2023）。

bash
# PyDESeq2 のインストール
pip install pydeseq2

💡 R と Python のどちらを選ぶか
bulk RNA-seq の発現変動解析では R が事実上の標準で、DESeq2・edgeR の情報量と論文実績が圧倒的です。この工程だけ R を使うか、Python で通したいなら PyDESeq2 を選ぶ、というのが現実的な判断になります。なお PyDESeq2 と DESeq2 は独立した開発チームによるもので、結果は近いものの完全一致はしません。

6. 注意点

生の p 値ではなく調整p値で判断する：多重検定の補正は必須です。
正規化済みの値を入力しない：DESeq2 には生のカウントを渡します。TPM や FPKM を入力すると統計モデルの前提が崩れます。
生物学的反復が必要：技術的反復ではばらつきを推定できません。各群に最低2、できれば3以上の生物学的反復を用意します。
DEG の数が0でも失敗ではない：本当に差がない場合もあります。PCA でサンプルが条件ごとに分離しているかを先に確認します。

まとめ

発現変動解析は、条件間の差が統計的に有意かを検定し DEG を同定する工程。
少サンプルには DESeq2、複雑なデザインには edgeR、大規模には limma-voom。
DESeq2 は読み込み → サンプル情報 → DESeq() → results() の4ステップ。
padj < 0.05 かつ |log2FoldChange| > 1 が一般的な DEG の基準。
Python で通すなら PyDESeq2。ただし bulk では R が標準。

参考文献

Love, M. I., Huber, W., & Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology, 15(12), 550. doi:10.1186/s13059-014-0550-8
Robinson, M. D., McCarthy, D. J., & Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1), 139–140. doi:10.1093/bioinformatics/btp616
Law, C. W., Chen, Y., Shi, W., & Smyth, G. K. (2014). voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology, 15(2), R29. doi:10.1186/gb-2014-15-2-r29
Ritchie, M. E., Phipson, B., Wu, D., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Research, 43(7), e47. doi:10.1093/nar/gkv007
Muzellec, B., Teleńczuk, M., Cabeli, V., & Andreux, M. (2023). PyDESeq2: a python package for bulk RNA-seq differential expression analysis. Bioinformatics, 39(9), btad547. doi:10.1093/bioinformatics/btad547