DESeq2の使い方｜nf-core/rnaseqのカウント行列から発現変動遺伝子を検出する

カウント行列が手に入ったら、次は「どの遺伝子が群間で変動したか」を統計的に判定する段階です。DESeq2 は、そのための標準的なRパッケージです。

この記事では、nf-core/rnaseqでFASTQからカウント行列を作る記事の続きとして、そのカウント行列をDESeq2に渡し、発現変動遺伝子のリストを得るまでを通しで説明します。nf-coreを使っていなくても構いません。 featureCountsなど他の方法で作った生カウント行列でも手順は同じで、変わるのは「読み込み方」の1節だけです。

掲載しているコードと数値・図は、すべて手元のR 4.5.1／DESeq2 1.50.2で実際に実行した出力です。題材にはBioconductorの airway パッケージを使うので、FASTQもマッピングも不要で、この記事のコードをそのまま上から順に実行できます。

DESeq2とは：なぜ生カウントを渡すのか

TPM・FPKMを渡してはいけない

最初につまずきやすいのがここです。DESeq2に渡すのは、正規化されていない生のカウント値です。TPMやFPKM、RPKMを渡してはいけません。

理由は、DESeq2が内部で正規化と統計モデリングを同時に行っているからです。DESeq2は各サンプルのライブラリサイズをsize factorとして推定し、その上で「カウントが整数で、平均が大きいほど分散も大きい」という性質を負の二項分布でモデル化します。すでに長さや深さで割られた値を渡すと、この前提が壊れます。 TPMは小数で、しかもサンプルごとに合計が100万に固定されているため、「読まれた回数」という情報がすでに失われています。

Salmonの推定カウントが小数を含むのは問題ありません。tximport が長さ情報と一緒に渡し、DESeq2側がオフセットとして扱うためです。「割り算で正規化済みの値」がだめであって、「小数」がだめなのではありません。

負の二項分布とdispersion

RNA-seqのカウントは、同じ条件のサンプル間でもばらつきます。このばらつきの大きさをdispersion（分散パラメータ）と呼びます。

問題は、生物学的反復が3〜4個しかない実験では、遺伝子ごとにdispersionを正確に推定できないことです。DESeq2は、全遺伝子から「発現量とdispersionの関係」の傾向線を学習し、各遺伝子の推定値をその線に引き寄せます（shrinkage）。反復が少なくても極端な推定値に振り回されにくいのは、この仕組みのためです。

環境を準備する

RとRStudioが入っている前提で進めます。DESeq2はCRANではなくBioconductorのパッケージなので、install.packages("DESeq2") ではインストールできません。

install.packages("BiocManager")

BiocManager::install(c("DESeq2", "apeglm", "airway", "pheatmap"))

apeglm は後で使うfold changeの補正に、airway はこの記事の題材データに、pheatmap はヒートマップの描画に使います。読み込んでエラーが出なければ完了です。

library(DESeq2)
library(airway)

入力をそろえる

DESeq2に必要なのは次の2つだけです。

• カウント行列 — 行が遺伝子、列がサンプルの整数行列

• colData — 行がサンプル、列が群分けやバッチなどの属性を表す表

この2つのサンプルの並び順が一致していることが絶対条件です。ここがずれたまま実行すると、エラーは出ないのに結果が完全に無意味になります。

nf-core/rnaseq の出力から読み込む

nf-core/rnaseqの star_salmon 経路を使った場合、results/star_salmon/ に統合済みの行列が出ています。最も簡単なのは、パイプラインが出力する SummarizedExperiment をそのまま読むことです。

se <- readRDS("results/star_salmon/salmon.merged.gene_counts_length_scaled.rds")

dds <- DESeqDataSet(se, design = ~ condition)

TSVから読む場合は、1列目が遺伝子ID、2列目が遺伝子名になっている点に注意します。カウント本体は3列目以降です。

tab <- read.delim("results/star_salmon/salmon.merged.gene_counts_length_scaled.tsv",
                  row.names = 1, check.names = FALSE)

cts <- as.matrix(tab[, -1])       # gene_name 列を落とす
mode(cts) <- "integer"            # DESeq2は整数を要求する

Salmonの quant.sf から自分で集計する場合は tximport を使いますが、nf-coreを使っているなら上の統合済みファイルで足ります。

featureCounts のカウント行列から読み込む

featureCountsの出力は、先頭6列がアノテーション情報（Geneid、Chr、Start、End、Strand、Length）で、7列目以降がサンプルごとのカウントです。

fc <- read.delim("counts.txt", comment.char = "#", row.names = 1)

cts <- as.matrix(fc[, 6:ncol(fc)])   # 先頭のアノテーション列を落とす

colnames(cts) <- sub("\\.bam$", "", basename(colnames(cts)))

featureCountsは列名にBAMのフルパスをそのまま入れるので、最後の行で短いサンプル名に整えています。

colData を用意する

colData はサンプル情報の表です。行名がカウント行列の列名と一致している必要があります。

coldata <- data.frame(
  condition = factor(c("control", "control", "control",
                       "treated", "treated", "treated")),
  row.names = colnames(cts)
)

coldata$condition <- relevel(coldata$condition, ref = "control")

relevel() は比較の基準をどちらに置くかの指定です。これを忘れると、Rはアルファベット順で最初の水準を基準にします。control と treated なら偶然うまくいきますが、KO と WT だと KO が基準になり、fold changeの向きが期待と逆になります。

この記事で使うデータ

以降は airway パッケージのデータを使います。4人のドナー由来の気道平滑筋細胞に、ステロイド（デキサメタゾン）を処理した群と未処理群があり、同じ細胞株から処理・未処理がペアで取られた8サンプルという構成です。

data(airway)
se <- airway
se$dex <- relevel(se$dex, ref = "untrt")

as.data.frame(colData(se)[, c("cell", "dex")])

              cell   dex
SRR1039508  N61311 untrt
SRR1039509  N61311   trt
SRR1039512 N052611 untrt
SRR1039513 N052611   trt
SRR1039516 N080611 untrt
SRR1039517 N080611   trt
SRR1039520 N061011 untrt
SRR1039521 N061011   trt

cell が細胞株（＝ドナー）、dex が処理の有無です。この「同じドナーからペアで取られている」構造は、後のデザイン式で効いてきます。

DESeq2を実行する

DESeqDataSet を作る

カウント行列とcolDataからDESeq2のオブジェクトを作ります。行列から作る場合は DESeqDataSetFromMatrix()、SummarizedExperiment から作る場合は DESeqDataSet() です。

# 行列から作る場合
dds <- DESeqDataSetFromMatrix(countData = cts,
                              colData   = coldata,
                              design    = ~ condition)

# airway（SummarizedExperiment）から作る場合
dds <- DESeqDataSet(se, design = ~ cell + dex)

design は「発現量の違いを何で説明するか」を書く数式です。ここでは ~ cell + dex として、細胞株の違い（ドナー間の個体差）を考慮した上で処理の効果を見ています。デザイン式は後の節で詳しく扱います。

低発現遺伝子をフィルタする

ほとんど読まれていない遺伝子は、検定してもノイズにしかなりません。先に落としておきます。

nrow(dds)
# [1] 63677

keep <- rowSums(counts(dds) >= 10) >= 4
dds  <- dds[keep, ]

nrow(dds)
# [1] 16139

63,677遺伝子のうち47,538が除かれ、16,139遺伝子まで絞られました。 条件は「カウントが10以上のサンプルが4つ以上ある遺伝子を残す」です。この 4 は、比較したい群のうち小さい方のサンプル数に合わせます。 片群にだけ発現している遺伝子を落とさないための設定です。

DESeq() を実行する

正規化、dispersion推定、検定までを1行で実行します。

dds <- DESeq(dds)

estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

この6行が、DESeq2が実際にやっていることです。順に、ライブラリサイズの推定、遺伝子ごとのdispersion推定、発現量とdispersionの関係の学習、その線への引き寄せ、そしてモデルの当てはめと検定です。

推定されたsize factorは確認できます。

round(sizeFactors(dds), 3)

SRR1039508 SRR1039509 SRR1039512 SRR1039513 SRR1039516 SRR1039517 SRR1039520 SRR1039521 
     1.022      0.896      1.174      0.668      1.175      1.399      0.917      0.945 

1.0前後に収まっていれば、サンプル間の読み取り深度に極端な偏りはありません。0.5を下回る、あるいは2を超えるサンプルがある場合は、そのサンプルのシーケンス量やQCを疑ってください。

結果を読む

results() で結果を取り出す

比較したい対比（contrast）を指定して結果を取り出します。

res <- results(dds, contrast = c("dex", "trt", "untrt"), alpha = 0.05)

summary(res)

out of 16139 with nonzero total read count
adjusted p-value < 0.05
LFC > 0 (up)       : 2190, 14%
LFC < 0 (down)     : 1891, 12%
outliers           : 0, 0%
low counts         : 313, 1.9%

contrast = c("dex", "trt", "untrt") は「dex 列について、trt を untrt と比べる」という意味です。後ろに書いた方が基準になります。alpha = 0.05 は、後述する独立フィルタリングの最適化に使われる閾値なので、実際に使うFDR水準と合わせておきます。

結果、padj < 0.05 の遺伝子が4,081個（上昇2,190・低下1,891）検出されました。

log2FoldChange と padj の意味

results() が返す主要な列は次のとおりです。

絞り込みに使うのは pvalue ではなく padj です。16,139遺伝子を同時に検定しているので、p < 0.05 だけで選ぶと偶然の当たりが数百個単位で混ざります。

padj順に並べると、上位はこうなりました。

library(org.Hs.eg.db)

res$symbol <- mapIds(org.Hs.eg.db, keys = rownames(res),
                     column = "SYMBOL", keytype = "ENSEMBL", multiVals = "first")

head(res[order(res$padj), c("symbol", "baseMean", "log2FoldChange", "padj")], 5)

                 symbol   baseMean log2FoldChange          padj
ENSG00000152583 SPARCL1   997.9612       4.570944 5.641707e-132
ENSG00000165995  CACNB2   495.5675       3.287083 5.641707e-132
ENSG00000120129   DUSP1  3411.9661       2.943859 1.273761e-124
ENSG00000101347  SAMHD1 12712.9456       3.763021 2.968002e-124
ENSG00000189221    MAOA  2343.5733       3.349721 3.609424e-119

DUSP1 はステロイド応答で古くから知られる遺伝子で、既知の生物学と一致しています。解析が正しく回ったかどうかは、統計値ではなく「既知の応答遺伝子が上位に来ているか」で確認するのが実務的です。

lfcShrink でfold changeを補正する

results() が返す log2FoldChange には落とし穴があります。発現量が低い遺伝子ほど、fold changeが極端な値になりやすいのです。カウントが5と50なら10倍ですが、この差はほぼ偶然でも起こります。

lfcShrink() は、この不確かさを考慮してfold changeを0方向に引き寄せます。

resLFC <- lfcShrink(dds, coef = "dex_trt_vs_untrt", type = "apeglm")

coef に指定できる名前は resultsNames(dds) で確認できます。

resultsNames(dds)
# [1] "Intercept"               "cell_N061011_vs_N052611"
# [3] "cell_N080611_vs_N052611" "cell_N61311_vs_N052611" 
# [5] "dex_trt_vs_untrt"

効果はMAプロットで一目で分かります。横軸が平均発現量、縦軸がfold changeです。

▲ 左が補正なし、右が lfcShrink(type = "apeglm") 後。赤が padj < 0.05 の遺伝子。左図の左端（低発現側）で上下に飛び散っている点が、右図では中央に引き戻されているのが分かります。

補正後の値は、遺伝子を順位づけしたり、fold changeの閾値で絞ったりするときに使います。

volcanoプロットで全体を眺める

横軸にfold change、縦軸に -log10(padj) を取ると、全体像が一目で分かります。

library(ggplot2)

vd <- as.data.frame(resLFC)
vd <- vd[!is.na(vd$padj), ]
vd$status <- ifelse(vd$padj < 0.05 & vd$log2FoldChange >  1, "up",
             ifelse(vd$padj < 0.05 & vd$log2FoldChange < -1, "down", "ns"))

ggplot(vd, aes(log2FoldChange, -log10(padj), color = status)) +
  geom_point(size = 1.1, alpha = 0.6) +
  geom_vline(xintercept = c(-1, 1), linetype = "dashed") +
  geom_hline(yintercept = -log10(0.05), linetype = "dashed") +
  theme_bw()

▲ lfcShrink 後の値で描画。padj < 0.05 かつ |log2FC| > 1 は、上昇411・低下338でした。

ここで数字がずれることに気づいてください。補正前の res で同じ条件（padj < 0.05 かつ |log2FC| > 1）を数えると951個ですが、補正後は749個に減ります。低発現で極端なfold changeを示していた遺伝子が、shrinkageによって閾値の内側に戻ったためです。どちらが正しいという話ではなく、どちらの値で数えたのかを記録しておくことが重要です。

サンプルの品質を確認する

発現変動遺伝子のリストを見る前に、サンプルそのものが信頼できるかを確認します。この順序を逆にすると、外れサンプルが作った偽のシグナルを一生懸命解釈することになります。

分散を安定化してから可視化する

カウントデータをそのままPCAにかけると、発現量の大きい遺伝子だけで結果が決まってしまいます。可視化の前に vst() で分散を安定化します。

vsd <- vst(dds, blind = FALSE)

blind = FALSE は「デザイン式を知った上で変換する」という指定です。QCのために群情報を一切使いたくない場合は TRUE にしますが、通常のQCでは FALSE で構いません。

PCAで群構造と外れサンプルを見る

plotPCA(vsd, intgroup = c("dex", "cell"))

▲ 色が処理の有無、形が細胞株。PC1（48%）が処理の効果、PC2（24%）が細胞株の違いに対応しているのが読み取れます。

この図は、解析を進めてよいという判断材料になります。見るべき点は3つです。

• 群がPC1またはPC2で分離しているか — 分離していれば、処理の効果がデータ全体に現れている

• 群から大きく外れたサンプルがないか — 1点だけ遠く離れていれば、そのサンプルのQCを再確認する

• 群以外の構造が見えていないか — ここではPC2に細胞株の差が出ており、デザイン式に cell を入れるべき根拠になっている

サンプル間距離で確認する

PCAは上位2軸しか見ていないので、全体の距離関係も確認します。

library(pheatmap)
library(RColorBrewer)

sampleDists <- dist(t(assay(vsd)))
m <- as.matrix(sampleDists)
rownames(m) <- paste(vsd$cell, vsd$dex, sep = " / ")
colnames(m) <- NULL

pheatmap(m,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colorRampPalette(rev(brewer.pal(9, "Blues")))(255))

▲ 濃いほど距離が近い。未処理4サンプルと処理4サンプルが、それぞれ別のクラスタにきれいに分かれています。

処理の有無ではなく、バッチや測定日でクラスタが分かれていたら、それはデザイン式に入れるべき交絡因子です。

dispersionの当てはまりを見る

モデルが妥当かどうかは、dispersionプロットで確認できます。

plotDispEsts(dds)

▲ 黒が遺伝子ごとの生の推定値、赤が当てはめた傾向線、青が縮小後の最終値。黒い点が赤い線の周りに散らばり、青い点が線に引き寄せられているのが正常な形です。

赤い線から系統的に外れていたり、青い点が線に乗りすぎていたりする場合は、デザイン式かサンプル構成に問題がある可能性があります。

デザイン式を組む

交絡因子を式に入れる

デザイン式は、発現量の違いを何で説明するかの宣言です。airwayのデータでは同じドナーから処理・未処理をペアで取っているため、cell を入れています。

design(dds) <- ~ cell + dex

順序には意味があります。一番右に、効果を知りたい変数を置きます。 results() は既定で一番右の変数の効果を返すためです。

cell を入れた効果は数字に出ます。同じデータ・同じフィルタで式だけを変えると、こうなりました。

cell を落とすと検出数が4,081から2,773へ約3割減ります。ドナー間の個体差が説明されないまま誤差に押し込まれ、処理の効果が埋もれるためです。ペア構造や既知のバッチがあるなら、必ず式に入れてください。

結果を書き出して次へ進む

絞り込んで保存します。

resOrdered <- resLFC[order(resLFC$padj), ]

sig <- subset(resOrdered, padj < 0.05 & abs(log2FoldChange) > 1)
nrow(sig)

write.csv(as.data.frame(resOrdered), "deseq2_all_genes.csv")
write.csv(as.data.frame(sig),        "deseq2_significant.csv")

全遺伝子の結果も一緒に保存してください。 有意な遺伝子だけを残すと、後から閾値を変えたくなったときや、特定の遺伝子が有意でなかったことを示したいときに、解析をやり直すことになります。

再現のために残しておく情報は次の4点です。

1. カウント行列の由来（nf-coreの版とパラメータ、またはfeatureCountsのコマンド）

2. 低発現フィルタの条件

3. デザイン式と、relevel() で指定した基準水準

4. sessionInfo() の出力（RとDESeq2のバージョン）

ここまでで遺伝子リストが得られました。次は、そのリストが生物学的に何を意味するかを調べる段階です。共発現ネットワークから疾患に関わるモジュールを探す方法は、別記事にまとめています。

関連記事WGCNA解析 part 1：RNA-seqデータを使って、疾患に絡む遺伝子ネットワークを可視化する →

ケース別：よく出るエラーと対処

よくある質問

Q. DESeq2とedgeR、limma-voomはどう使い分ける？ どれも生カウントを入力とする標準的な手法で、同じデータならおおむね似た結果になります。 DESeq2とedgeRは負の二項分布を仮定し、反復が少ない実験に強いのが特徴です。limma-voomはカウントを重み付き線形モデルに乗せる方式で、サンプル数が多いとき（数十以上）に高速で、複雑なデザイン式を扱いやすい利点があります。迷ったら、反復3〜6程度の一般的な実験ならDESeq2で問題ありません。

Q. TPMをDESeq2に渡してはいけない理由をもう一度 DESeq2は「カウントが整数で、平均が大きいほど分散も大きい」という性質を前提にモデルを組み、内部でsize factorによる正規化まで行います。TPMはすでにサンプルごとの合計が固定されており、「何回読まれたか」という情報が失われています。 分散の構造も壊れるため、統計的な前提が成立しません。ヒートマップや発現量の比較にTPMを使うのは適切ですが、検定の入力にはなりません。

Q. nf-coreの gene_counts.tsv と gene_counts_length_scaled.tsv のどちらを使う？ gene_counts_length_scaled.tsv です。転写産物長の構成がサンプル間で違う場合の影響を補正した上で、カウントのスケールを保っています。gene_counts.tsv でも動きますが、アイソフォームの使い分けが群間で変わる遺伝子で偏りが残ります。SummarizedExperiment のRDSがあるなら、それを DESeqDataSet() に渡すのが最も確実です。

Q. 低発現フィルタの閾値はどう決める？ この記事では「カウント10以上のサンプルが4つ以上」としました。4 は比較したい群のうち小さい方のサンプル数に合わせます。 片群にだけ強く発現する遺伝子を落とさないためです。10 は目安で、シーケンス深度が浅いなら下げます。厳密な正解はありませんが、決めた条件を記録しておくことが重要です。

Q. padj が NA になる遺伝子があるのはなぜ？ DESeq2の独立フィルタリングによって、検出力が低く検定しても意味がないと判断された遺伝子が除外されるためです。除外することで検定数が減り、残りの遺伝子のFDR補正が緩くなります。エラーではないので、0や1で埋めないでください。極端な外れ値を持つ遺伝子（Cook距離が大きい）も NA になることがあります。

Q. pvalue で絞ってはいけない？ いけません。数万遺伝子を同時に検定しているため、p < 0.05 だけでは偶然有意になったものが大量に混ざります。16,139遺伝子なら、本当は何も変動していなくても800個前後がp < 0.05 になります。 必ず padj（FDR補正後）で絞ってください。

Q. fold changeの閾値は必要？ 必須ではありません。padj だけで絞るのが統計的には素直で、|log2FC| > 1 のような閾値は「生物学的に意味のある大きさ」という別の基準です。閾値を使うなら lfcShrink() 後の値で判断してください。補正前の値は低発現遺伝子で過大になっています。

Q. lfcShrink の type は何を選べばいい？ apeglm が現在の推奨です。coef で係数名を指定する必要があり、resultsNames(dds) で確認できます。contrast を使った複雑な比較で apeglm が使えない場合は ashr を選びます。古い normal は後方互換のために残されているだけなので、新しく書くコードでは選ばないでください。

Q. 3群以上を比較したい results() の contrast で比較したいペアを指定します。contrast = c("condition", "B", "A") のように書けば、A を基準に B を比較できます。「どれか1つでも違う遺伝子」を探したい場合は、DESeq(dds, test = "LRT", reduced = ~ 1) で尤度比検定を使います。この場合の log2FoldChange は特定の2群間の値なので、解釈に注意してください。

Q. バッチ効果はどう扱う？ デザイン式に入れます（~ batch + condition）。カウントを事前に補正してからDESeq2に渡してはいけません。 DESeq2はモデルの中でバッチを考慮するので、二重に補正することになります。可視化のためだけにバッチを除いた行列が欲しい場合は、limma::removeBatchEffect() をvst変換後のデータに適用します。あくまで図のためであって、検定には使いません。

Q. model matrix is not full rank が出た バッチと群が完全に交絡しているか、デザイン式に不要な列が入っています。前者は実験デザインの問題で、解析では救えません。 例えばcontrolを全部1日目、treatedを全部2日目に測定した場合、日付の効果と処理の効果は原理的に区別できません。後者なら、片方の群にしか水準が存在しない列をデザイン式から外します。

Q. 生物学的反復が各群2サンプルしかない 実行はできますが、dispersionの推定が不安定になり、検出力が大きく落ちます。 DESeq2は全遺伝子から情報を借りて補正するので破綻はしませんが、結果は探索的なものとして扱ってください。3反復以上が最低ラインで、4以上あると安定します。テクニカルレプリケート（同じサンプルを2回シーケンス）は生物学的反復の代わりになりません。

Q. 結果に遺伝子名（シンボル）を付けたい org.Hs.eg.db（ヒトの場合）で変換します。BiocManager::install("org.Hs.eg.db") の後、mapIds(org.Hs.eg.db, keys = rownames(res), column = "SYMBOL", keytype = "ENSEMBL") です。ENSEMBL IDにバージョン番号（.15 など）が付いている場合は、sub("\\..*$", "", rownames(res)) で先に落としてください。 付いたままだと変換に失敗して全部 NA になります。

Q. StringTie・Ballgownの結果からDESeq2に渡せる？ そのままでは渡せません。Ballgownが扱うのはFPKMなどの連続値なので、DESeq2の入力にはなりません。StringTieを使うなら、-e オプション付きで実行した結果を prepDE.py でカウント行列に変換する経路があります。ただし遺伝子レベルの発現変動が目的なら、SalmonやfeatureCountsでカウントを取る方が素直です。

Q. 図の見た目を論文用に整えたい plotPCA(..., returnData = TRUE) でデータフレームを取り出せば、あとはggplot2で自由に描けます。volcanoは EnhancedVolcano パッケージが手軽です。ヒートマップは pheatmap か ComplexHeatmap を使います。どの図もvst変換後のデータから描き、生カウントからは描かないようにしてください。

関連記事

カウント行列を作るところからやり直したい場合は、こちらが上流にあたります。

関連記事nf-core/rnaseqの使い方｜NextflowでRNA-seq解析を自動化する →

アラインメントの中身を手で確かめたい場合はこちらです。

関連記事超初心者向け！！RNA-seq解析シリーズ④HISAT2でマッピングする →

転写産物・アイソフォーム単位の発現差を見たい場合は、Ballgownを使う経路になります。

関連記事超初心者向け！！RNA-seq解析シリーズ⑥ Ballgownで発現差解析 →

得られた遺伝子リストから共発現ネットワークへ進む場合はこちらです。

関連記事WGCNA解析 part 1：RNA-seqデータを使って、疾患に絡む遺伝子ネットワークを可視化する →

公式情報

外部リンクDESeq2 公式vignette（Analyzing RNA-seq data with DESeq2） ↗https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

外部リンクLove MI, Huber W, Anders S. (2014) Genome Biology — DESeq2原著論文 ↗https://doi.org/10.1186/s13059-014-0550-8

外部リンクBioconductor RNA-seq workflow（rnaseqGene） ↗https://bioconductor.org/packages/release/workflows/vignettes/rnaseqGene/inst/doc/rnaseqGene.html