今回は、遺伝子発現データから差次的発現遺伝子（DEG）を見つけるための強力なツール「edgeR」について解説します。
edgeRは、少ないサンプル数でも信頼性の高い解析が可能で、TMM正規化やGLMモデルなど多彩な統計手法を提供しているのが特徴です。 この記事では、GEOからのデータ取得、前処理、edgeRによる発現変動解析、そして結果の解釈まで、実践的なコードとともにステップバイステップで紹介していきます。
RNA-seqデータを使って統計的に意味のある発現変化を見つけたい方は、ぜひ参考にしてください。
macOS Monterey (12.4), Quad-Core Intel Core i7, Memory 32GB, edgeR 4.4.2
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！
公共データを用いたシングルセル ダイナミクス解析に関する初心者向け技術書を販売中
シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！
edgeRによる発現変動解析
edgeRはRNA-seqデータから差次的発現遺伝子（DEG）を同定するための効率的なツールです。DESeq2と同様にネガティブバイノミアル分布を使用しますが、異なる正規化手法と統計的アプローチを採用しています。edgeRの主な特徴は以下の通りです：
- TMM（Trimmed Mean of M-values）法による効果的なサンプル間の正規化
- タグワイズ分散推定による遺伝子ごとの生物学的変動の正確な評価
- 少ないサンプル数でも信頼性の高い結果を得るための経験ベイズ法
- 一般化線形モデル（GLM）を用いた複雑な実験デザインへの対応
- 正確検定（exact test）やQLF検定など、状況に応じた統計的検定手法
- 計算効率の良いアルゴリズムによる大規模データセットの迅速な処理
この記事では、edgeRを使用して発現変動遺伝子を特定し、様々な手法による結果の違いや解釈方法について説明します。
また前回DEseq2に関する解説も行いました。今回はDEseq2とedgeRの比較も検証してみたいと思います。
https://labo-code.com/bioinformatics/rna-seq-deseq2/
解析準備
GEOからカウントデータを取得
今回の解析に使うカウントデータを取得してください。
# GEOqueryパッケージを読み込む
library(GEOquery)
# GEOからカウントデータを取得
gse <- getGEO("GSE141413", GSEMatrix = FALSE)
# カウントデータを含む補足ファイルをダウンロード
supp_dir <- "data/geo_suppl"
dir.create(supp_dir, recursive = TRUE, showWarnings = FALSE)
count_file <- getGEOSuppFiles("GSE141413", baseDir = supp_dir, makeDirectory = FALSE)
# ダウンロードしたカウントデータを読み込む
counts_path <- file.path(supp_dir, "GSE141413_gene_count.txt.gz")
counts_data <- read.delim(counts_path, header = TRUE)
カウントデータの前処理
取得したカウントデータを解析に適した形式に整理します。このデータセットには発現量データと遺伝子アノテーション情報が含まれています。
# ステップ1: カウントデータの確認
dim(counts_data)
head(counts_data[, 1:8])
# ステップ2: データの構造確認
colnames(counts_data)
str(counts_data)
# ステップ3: 発現データとアノテーション情報の分離
count_cols <- c("Control", "NC_LM", "st2_LM", "IFNG", "st2_IFNG")
annotation_cols <- c("gene_name", "gene_chr", "gene_start", "gene_end", "gene_strand",
"gene_length", "gene_biotype", "gene_description", "tf_family")
# サブセットの作成
counts <- counts_data[, count_cols]
gene_annotation <- counts_data[, annotation_cols]
# ステップ4: 行名の設定
rownames(counts) <- rownames(counts_data)
rownames(gene_annotation) <- rownames(counts_data)
# ステップ5: サンプル情報の作成
sample_info <- data.frame(
sample_id = colnames(counts),
condition = c("Control", "NC_LM", "st2_LM", "IFNG", "st2_IFNG"),
group = c("Control", rep("LM", 2), rep("IFNg", 2)),
treatment = c("None", "NC", "st2", "None", "st2"),
row.names = colnames(counts)
)
# ステップ6: フィルタリングとデータ保存
# 空でない遺伝子のみを保持
counts <- counts[rowSums(counts) > 0, ]
gene_annotation <- gene_annotation[rownames(counts), ]
# 元の条件が複雑すぎてエラーになるため、シンプルな条件に変換
simplified_condition <- ifelse(grepl("Control", sample_info$condition), "Control", "Treatment")
sample_info$simple_condition <- factor(simplified_condition, levels = c("Control", "Treatment"))
これで準備完了です。
edgeRを使ったワークフロー
edgeRは標準的な発現変動解析のワークフローを提供しており、データの正規化から統計的検定まで一連の処理を効率的に実行できます。以下のコードは、edgeRの基本的なワークフローを示しています。
# 必要なライブラリの読み込み
library(edgeR)
# DGEListオブジェクトの作成（シンプル化した条件を使用）
dge <- DGEList(counts = counts, group = sample_info$simple_condition)
# 低発現遺伝子のフィルタリング
keep <- filterByExpr(dge)
dge <- dge[keep, , keep.lib.sizes = FALSE]
# 正規化係数の計算
dge <- calcNormFactors(dge, method = "TMM")
# 分散の推定
dge <- estimateDisp(dge)
# 差次的発現テスト（exactTest）
# Control vs Treatment の比較
et <- exactTest(dge, pair = c("Control", "Treatment"))
# 結果の取得
res_edgeR <- topTags(et, n = Inf)
# 結果の要約
summary(decideTests(et))
# Treatment-Control
# Down 1642
# NotSig 13062
# Up 1335
# 有意なDEGの数
sum(res_edgeR$table$FDR < 0.05) # FDR < 0.05の遺伝子数
# [1] 2977
このコードブロックでは、以下の重要なステップを実行しています：
- DGEListオブジェクトの作成：
DGEList関数はカウントデータと実験グループ情報からedgeR用のデータ構造を作成します。countsパラメータには生のカウント行列を指定します。groupパラメータにはサンプルの条件情報（factor型）を指定します。
- 低発現遺伝子のフィルタリング：
filterByExpr関数は統計的に意味のある発現変動を検出できるレベルで発現している遺伝子を自動的に特定します。- 返される論理ベクトル
keepを使って、保持する遺伝子を選択します。 keep.lib.sizes = FALSEオプションは、フィルタリング後にライブラリサイズを再計算することを指示します。
- 正規化係数の計算：
calcNormFactors関数はTMM（Trimmed Mean of M-values）法を使用してサンプル間の組成バイアスを調整します。method = "TMM"はデフォルト設定で、他にも”RLE”や”upperquartile”などのオプションがあります。- このステップにより、サンプル間のライブラリサイズの違いと組成バイアスが調整されます。
- 分散の推定：
estimateDisp関数は共通分散（common dispersion）、遺伝子ごとの分散（tagwise dispersion）、トレンド分散（trended dispersion）を推定します。- この分散推定はネガティブバイノミアル分布に基づいた正確な統計的検定のために不可欠です。
- 分散推定は経験ベイズ法を使用して、少ないサンプル数でも信頼性の高い推定を行います。
- 差次的発現テスト：
exactTest関数は条件間の発現変動を検定します。pairパラメータには比較する条件の組を指定します（ここでは「Control」と「Treatment」）。- この関数は正確検定（exact test）を実行し、各遺伝子の発現変動の統計的有意性を評価します。
- 結果の取得と要約：
topTags関数は検定結果を整理し、p値で順位付けします。n = Infは全ての遺伝子の結果を取得することを意味します。decideTestsとsummary関数により、統計的に有意な発現変動遺伝子の数をカウントします。- 結果はlog2スケールでの発現変化（logFC）、平均発現量（logCPM）、p値（PValue）、多重検定補正後のp値（FDR）などを含みます。
edgeRのワークフローは、サンプル間の生物学的・技術的変動を適切に考慮し、信頼性の高い発現変動遺伝子の同定を可能にします。
GLM（一般化線形モデル）を用いた複雑なデザイン
edgeRのもう一つの強力な機能は、一般化線形モデル（GLM）を使用して複雑な実験デザインに対応できる点です。GLMアプローチは、バッチ効果や交互作用など、複数の要因を考慮した解析を可能にします。
# デザイン行列の作成
design <- model.matrix(~ simple_condition, data = sample_info)
# 分散の推定
dge <- estimateDisp(dge, design)
# GLMフィッティング
fit <- glmQLFit(dge, design)
# 差次的発現テスト
qlf <- glmQLFTest(fit, coef = "simple_conditionTreatment")
# 結果の取得
res_edgeR_glm <- topTags(qlf, n = Inf)
head(res_edgeR_glm)
# Coefficient: simple_conditionTreatment
# logFC logCPM F PValue FDR
# Bex3 -3.211163 6.009465 199.2349 1.185237e-05 0.02369648
# Med12 6.631471 4.713399 165.3379 1.988956e-05 0.02369648
# ENSMUSG00000097815 5.804860 4.209378 150.5126 2.578899e-05 0.02369648
# Cela1 -3.875430 4.268405 148.4434 2.685287e-05 0.02369648
# Tm4sf19 6.486567 5.554821 147.5358 2.725051e-05 0.02369648
# Ankrd52 5.653018 5.420780 144.2497 2.899722e-05 0.02369648
GLMアプローチを使用した解析では、以下のステップを実行しています：
- デザイン行列の作成：
model.matrix関数は実験デザインを表現する行列を作成します。- 構文
~ simple_conditionはサンプルの条件のみに基づくモデルを指定しています。 - 複雑なデザインの場合、
~ batch + condition + batch:conditionのように複数の要因や交互作用を含めることができます。 - デザイン行列の各列は、モデルの中の要因やそれらの組み合わせを表します。
- 分散の推定：
estimateDisp関数はGLMアプローチ用の分散を推定します。- 指定されたデザイン行列を考慮して分散推定を行います。
- これにより、特定のデザインに対して最適化された分散推定が得られます。
- GLMフィッティング：
glmQLFit関数は準尤度F検定（quasi-likelihood F-test）用のフィッティングを行います。- この手法は小さなサンプルサイズの場合でも過分散を適切に処理します。
- 代替として
glmFit関数もありますが、glmQLFitは統計的検出力と型I誤差のバランスが良いとされています。
- 差次的発現テスト：
glmQLFTest関数は準尤度F検定を実行します。coefパラメータは検定対象の係数（デザイン行列の列）を指定します。ここでは「simple_conditionTreatment」（Treatmentの効果）を指定しています。- 複雑なコントラストをテストする場合は、
contrastパラメータを使用することもできます。
- 結果の取得と表示：
topTags関数で検定結果を取得し、p値で順位付けします。- 結果はlogFC（log2 fold-change）、logCPM（log2 counts-per-million）、F値（F-統計量）、PValue（p値）、FDR（false discovery rate）の情報を含みます。
GLMアプローチはより柔軟で強力な解析を可能にし、exactTestでは扱えない複雑な実験デザインに対応できます。特に、バッチ効果の調整、交互作用の検出、時系列データの解析などに適しています。
edgeRの結果解釈
edgeRによる解析結果は、各遺伝子の発現変動に関する詳細な統計情報を含んでいます。この情報を理解することで、生物学的に重要な発見につながります。
# 結果テーブルの最初の数行を表示
head(res_edgeR$table)
# logFC logCPM PValue FDR
# Cela1 -3.874993 4.268496 1.835498e-16 2.943956e-12
# Rpsa-ps2 -3.752911 4.866253 6.654905e-16 5.336901e-12
# Med12 6.631524 4.713376 3.730481e-15 1.905339e-11
# Bex3 -3.211216 6.009488 5.823562e-15 1.905339e-11
# Tm4sf19 6.486588 5.554808 5.939706e-15 1.905339e-11
# Mgst2 -4.670871 4.522703 1.590844e-14 4.252591e-11
edgeRの結果テーブルに含まれる各カラムは以下の情報を提供します：
- logFC：log2スケールでの発現変化（fold change）。正の値は上方制御（上群＞下群）、負の値は下方制御（上群＜下群）を示します。ここでは「Treatment vs Control」の比較であり、正の値はTreatmentで発現が高いことを意味します。
- logCPM：log2スケールでの平均発現量（counts per million）。この値が高いほど、その遺伝子の発現量が多いことを示します。高発現の遺伝子ほど、一般的に結果の信頼性が高くなります。
- PValue：未調整のp値。各遺伝子の差次的発現の統計的有意性を示します。値が小さいほど、発現変動が統計的に有意であることを示します。
- FDR：偽発見率（false discovery rate）。多重検定補正後のp値で、Benjamini-Hochberg法によって調整されています。通常は、FDR < 0.05または0.1を基準として有意な発現変動遺伝子を選択します。
これらの情報を組み合わせることで、統計的に有意で生物学的にも意味のある発現変動遺伝子を特定することができます。例えば、「FDR < 0.05かつ|logFC| > 1」というフィルタリング条件を設定することで、統計的に有意な2倍以上の発現変化を示す遺伝子を選別できます。
DESeq2とedgeRの比較
DESeq2とedgeRはどちらも広く使われている発現変動解析ツールですが、アルゴリズムや特性に違いがあります。以下のコードでは、両者の結果を比較しています。
# DESeq2とedgeRの結果を比較
# DESeq2結果をデータフレームに変換
res_deseq2_df <- as.data.frame(res)
res_deseq2_df$gene <- rownames(res_deseq2_df)
# edgeR結果をデータフレームに変換
res_edgeR_df <- as.data.frame(res_edgeR)
res_edgeR_df$gene <- rownames(res_edgeR_df)
# 共通の遺伝子を取得
common_genes <- intersect(rownames(res_deseq2_df), rownames(res_edgeR_df))
# DESeq2とedgeRのlog2FCの相関をプロット
plot_data <- data.frame(
DESeq2_log2FC = res_deseq2_df[common_genes, "log2FoldChange"],
edgeR_log2FC = res_edgeR_df[common_genes, "logFC"]
)
ggplot(plot_data, aes(x = DESeq2_log2FC, y = edgeR_log2FC)) +
geom_point(alpha = 0.5) +
geom_smooth(method = "lm", color = "red") +
theme_bw() +
ggtitle("DESeq2 vs edgeR: log2FC の比較") +
xlab("DESeq2 log2FoldChange") +
ylab("edgeR logFC")
# 相関係数の計算
cor(plot_data$DESeq2_log2FC, plot_data$edgeR_log2FC, method = "pearson")
# [1] 0.9918045
このコードブロックでは、DESeq2とedgeRの結果比較を以下のステップで行っています：
- 結果の準備：
- 両方のツールの結果をデータフレームに変換します。
rownamesの情報（通常は遺伝子ID）をカラムとして追加し、後の結合操作を容易にします。
- 共通遺伝子の特定：
intersect関数を使用して、両方のデータセットに共通する遺伝子識別子を特定します。- これにより、直接比較可能な遺伝子セットを作成します。
- 比較データの作成：
- 共通遺伝子について、DESeq2のlog2FoldChangeとedgeRのlogFCの値を含む新しいデータフレームを作成します。
- 両者のlog2スケールでの発現変化量を直接比較することができます。
- 相関の可視化：
ggplot2パッケージを使用して、DESeq2とedgeRのfold changeの散布図を作成します。geom_pointでデータポイントを表示し、alpha = 0.5で透明度を設定して重なりを可視化します。geom_smooth(method = "lm")で線形回帰直線を表示し、全体的な傾向を示します。- 適切なタイトルと軸ラベルを設定して、グラフを分かりやすくします。
- 相関係数の計算：
cor関数でPearson相関係数を計算し、両者の結果の一致度を数値化します。- 1に近い値は、両方のツールが非常に類似した結果を出していることを示します。
一般的な違いは以下の通りです:
計算効率: サンプル数が多い場合、edgeRの方が計算効率が良い傾向がある
正規化方法: DESeq2は中央値比（median-of-ratios）法、edgeRはTMM法を使用
分散推定: 両方ともネガティブバイノミアル分布を使用するが、推定方法が異なる
統計テスト: DESeq2はWald検定、edgeRは正確検定または尤度比検定を使用
フィルタリング: デフォルトのフィルタリング方法が異なる
最後に
いかがだったでしょうか。
edgeRを使うことで、RNA-seqデータから信頼性の高い発現変動遺伝子を効率的に抽出することができます。今回紹介した基本的なワークフローを参考に、ぜひご自身のデータでも差次的発現解析を実践してみてください。
公共データを用いたSingle Cell RNA-seq解析に関する初心者向け技術書を販売中
プログラミング初心者でも始められるわかりやすい解説！
RとSeuratで始めるSingle Cell RNA-seq解析！
公共データを用いたシングルセル ダイナミクス解析に関する初心者向け技術書を販売中
シングルセルデータの高度な解析であるTrajectory解析、RNA Velocity解析、空間トランスクリプトーム解析が環境構築方法から詳しく解説されています！
