この記事では、カウントマトリクスをもとに「どの遺伝子が条件間で発現量が変化しているか」を統計的に検定する「発現変動解析」について解説します。DESeq2・edgeR・limma-voomの3ツールの特徴と使い分け、DESeq2を使った実践コードを紹介します。また、DESeq2はR環境で実行されるため、現在使用ユーザーが多い、Pythonでの発現変動解析を行う場合の選択肢も紹介します。
発現変動解析とは？なぜ必要なのか
カウントで得られた遺伝子ごとのリード数を比較するだけでは、「この遺伝子は本当に発現が変化したのか、それとも偶然のばらつきなのか」を判断できません。発現変動解析（DEA: Differential Expression Analysis）では、統計モデルを使ってその差が統計的に有意かどうかを検定します。
FASTQファイル→トリミング→マッピング→カウント→発現変動解析（DESeq2）
主なツールの紹介と比較
発現変動解析で使われる主なツールはすべてR言語で動作します。
| ツール名 | 統計モデル | 特徴 |
|---|---|---|
| DESeq2 | 負の二項分布 | サンプル数が少ない場合に強い・初心者向け |
| edgeR | 負の二項分布 | 柔軟な実験デザインに対応・高機能 |
| limma-voom | 正規分布（voom変換後） | マイクロアレイ由来・大規模サンプルに強い |
💡 負の二項分布とは？
RNA-seqのカウントデータは「0が多い・ばらつきが大きい」という特徴があります。このような性質を適切に扱える統計モデルが負の二項分布です。DESeq2とedgeRはどちらもこのモデルを採用しています。
DESeq2
Bioconductorで提供されているR パッケージで、RNA-seq発現変動解析の世界標準ツールのひとつです。サンプル数が少ない（各群2〜3サンプル程度）場合でも安定した結果が得られるよう設計されています。
メリット
- 少サンプルでも安定した検定ができる
- ドキュメント・チュートリアルが充実している
- 正規化・検定を一貫したワークフローで実行できる
デメリット
- サンプル数が多い場合はedgeRより保守的になることがある
- 複雑な実験デザインにはやや不向き
このツールを選ぶ理由：RNA-seq解析を初めて学ぶ方や、生物学的サンプルが少ない実験（各群3サンプル程度）に最適です。
edgeR
DESeq2と同様に負の二項分布を使うツールですが、複数群の比較や複雑な実験デザイン（交互作用・共変量など）への対応が柔軟です。DESeq2と並んで最も引用されているツールのひとつです。
メリット
- 複雑な実験デザインにも対応できる
- 多群比較が得意
- 論文での使用実績が非常に多い
デメリット
- DESeq2よりコードがやや複雑
- 初心者には最初とっつきにくい
このツールを選ぶ理由：2群比較だけでなく多群比較や交互作用など、複雑な実験デザインを扱う場合に向いています。
limma-voom
もともとマイクロアレイ解析用に開発されたlimmaパッケージに、RNA-seqのカウントデータをvoom変換で対応させた方法です。サンプル数が多い大規模データで特に力を発揮します。
メリット
- サンプル数が多い場合に高い検出力を発揮
- マイクロアレイとRNA-seqを統一的に扱える
- limmaの豊富な機能をそのまま使える
デメリット
- サンプル数が少ないとDESeq2・edgeRに劣る場合がある
- voom変換の概念が初心者にはわかりにくい
このツールを選ぶ理由：サンプル数が10以上ある大規模コホートや、マイクロアレイとRNA-seqを混在して解析したい場合に適しています。
ツール選びのまとめ
| 状況 | おすすめの選択 |
|---|---|
| はじめて発現変動解析をする・各群3サンプル程度 | DESeq2 |
| 複雑な実験デザイン・多群比較 | edgeR |
| サンプル数が多い（10以上）・大規模コホート | limma-voom |
DESeq2の基本的な使い方
インストール
# BiocManagerでインストール
# BiocManagerでインストール if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager") BiocManager::install("DESeq2")
💡 BiocManagerとは？
バイオインフォマティクス向けRパッケージを管理するツールです。DESeq2やedgeRなど、CRAN（通常のRパッケージ配布場所）には登録されていない生命科学系のパッケージはBiocManagerを使ってインストールします。
ステップ① カウントマトリクスの読み込み
HTSeq-countで作成した counts.txt を読み込みます。
library(DESeq2) # カウントデータの読み込み count_data <- read.table( "counts.txt", header = TRUE, row.names = 1, sep = "\t" ) # __no_feature などの特殊行を除外 count_data <- count_data[!grepl("^__", rownames(count_data)), ]
ステップ② サンプル情報の設定
どのサンプルが「処理群」で、どのサンプルが「対照群」かを定義します。
# サンプル情報（条件）を定義
# サンプル情報（条件）を定義 col_data <- data.frame( condition = factor(c("control", "control", "treated", "treated")), row.names = colnames(count_data) )
ステップ③ DESeqDataSetの作成と解析実行
# DESeqDataSetオブジェクトの作成
# 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未満の遺伝子を解析前に除外します。ノイズの多い低発現遺伝子を減らすための前処理です。 |
ステップ④ 結果の取得と保存
# 結果の取得（treated vs control の比較）
# 結果の取得（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 の数を確認
💡 padj・log2FoldChange とは？
padj（調整済みp値）は多重検定補正を行ったp値です。通常 0.05 未満を有意とします。log2FoldChange は発現変化の大きさを対数スケールで表したもので、値が1なら「2倍に増加」、-1なら「2倍に減少」を意味します。
発展的なオプション（参考）
- lfcShrink()log2FoldChangeのシュリンケージ（縮小推定）を行う。低カウント遺伝子の過大推定を補正でき、 MA プロットや可視化に有効。
- vst() / rlog()カウントデータを分散安定化変換する。PCAやヒートマップなど可視化の前処理に使用する。
- plotMA()発現量と変化量の関係を示すMAプロットを描画する。DEGの全体像を把握するのに便利。
- plotPCA()サンプル間の類似性を主成分分析（PCA）で可視化する。外れサンプルの確認に使う。
- alpha = 0.1results() の有意水準を変更する（デフォルトは0.1）。論文によって0.05に変更することが多い。
Pythonで発現変動解析をする場合
DESeq2・edgeR・limma-voomはいずれもR言語のツールです。Pythonで発現変動解析を行いたい場合は、以下のツールが選択肢になります。
🐍 pydeseq2（おすすめ）
DESeq2のアルゴリズムをPythonで再実装したパッケージです。scikit-learnなどのPythonエコシステムと統合しやすく、Pythonに慣れた方には最もとっつきやすい選択肢です。pip install pydeseq2 でインストールできます。
🐍 diffxpy
Pythonネイティブの発現変動解析パッケージです。単細胞RNA-seq（scRNA-seq）にも対応しており、複雑な統計モデルを柔軟に組めます。
📌 現実的なアドバイス：
バイオインフォマティクスの発現変動解析はRが業界標準であり、DESeq2・edgeRの情報量・サポート・論文実績は圧倒的です。Python環境で解析を行っている方は、この工程だけRを使うかpydeseq2を選択するのが現実的です。
まとめ
- 発現変動解析は条件間の発現変化を統計的に検定し、DEGを同定する作業
- 初心者・少サンプルにはDESeq2、複雑なデザインにはedgeR、大規模サンプルにはlimma-voomが向いている
- DESeq2は「データ読み込み → サンプル情報設定 → DESeq() 実行 → 結果取得」の4ステップで完結する
- padj < 0.05 かつ |log2FoldChange| > 1 が一般的なDEG判定の基準
- Pythonで行う場合はpydeseq2が現時点での最有力候補
次のステップ：DEGが得られたら、GOエンリッチメント解析やKEGGパスウェイ解析などの機能解析へと進みましょう。
コメント