RNA-Seqの発現差解析では、edgeRやDESeq2のようなRNA-Seq専用の統計手法がよく使われます。 これらの手法は、Gene Countsデータの特徴に合わせて特別に設計された統計モデルとなっています。
一方で、 「edgeRを使うべき」 「DESeq2を使うべき」 「RNA-Seqではt検定を使ってはいけない」 といった単純化された断定が浸透しつつあります。 そこで本記事では、具体的なデータを使って、edgeR、DESeq2、そして前処理済みデータに対するt検定の結果を比較します。
今回のデータの特徴：ばらつきの小さいin vitro実験
今回使用しているGSE49110は、biopsyサンプルではなく、 MDA-MB-231細胞を用いた培養細胞ベースのRNA-Seqデータです。
ERRαを標的とするsiRNA処理細胞とcontrol siRNA処理細胞を比較しています。 各条件はduplicateで測定されており、N数は2しかありません。 一方で、同一条件内の繰り返し測定は非常に安定しており、 群内ばらつきが比較的小さいという特徴があります。
現実に、このようなタイプの研究データは数が非常に多いので、 まずは、このようなデータで比較し、結果を検証します。
本記事の目的は、どの手法が優れているかを明らかにしようというものではありません。 それぞれの手法を使ったときに、どのような発現パターンが有意判定されるのか、 あるいは有意判定されないのかを確認し、 手法選択や結果解釈の目安にすることを目的としています。
比較する3つの方法
ここでは、同じsiCとsiE1の比較に対して、次の3つの方法を比較しました。
- edgeRによる発現差解析
- DESeq2による発現差解析
- 前処理・正規化・フィルタリング後のデータに対するt検定
edgeRとDESeq2は、RNA-SeqのCountデータに合わせた統計モデルを使います。 一方、t検定は、前処理・正規化・フィルタリングを行った後の値に対して適用しています。 生のGene Countsや全測定データを対象にt検定を行ったわけではありません。 低Counts領域の不安定な遺伝子を減らし、 一般的な統計手法を適用しやすい状態に整えたうえで比較しています。
前処理・正規化・フィルタリングの具体的な処理は次の通りです。
- Log2変換（ここで測定が0だった場合は欠損値になります）
- 70 percentileでのグローバルノーマライズ
- 正規化後の値に対して、線形数換算で20未満を20に置き換え
- 欠損値を、線形数換算で16で埋める
- siC処理（コントロール）サンプルの平均値に対する比に変換
- Gene Countsが、すべての群において、1サンプル以上で20未満の遺伝子を除去
p値だけで見ると、手法ごとの違いは大きく見える
まず、edgeR、DESeq2、t検定について、 それぞれ p < 0.05 となった遺伝子を比較しました。
個別サンプル同士を比較した散布図だと、低Counts領域では、わずかな測定値の差が大きなfold changeになりやすく、 対角線から広がって、ばらつきが大きく見えます。
一方で、Fig1 の散布図は、siC群とsiE1群の各２回の繰り返しサンプルの平均値を比較しています。 群平均同士の比較では、個別サンプルのランダムなばらつきがある程度平均化されるため、 低Counts領域の点は、個別サンプル同士を比較した場合のようには広がって見えません。
ただし、低Counts領域で群平均の点が対角線付近に集まって見えることは、 その領域の測定値が安定していることを意味するわけではありません。 むしろこの領域では、真の発現差よりもランダムな測定ゆらぎの影響が大きく、 生物学的な発現差として解釈しにくい値が多く含まれていると考えられます。
edgeRやDESeq2は、Gene Counts全体から平均値と分散の関係を推定して有意性を評価します。 しかし、低Counts領域まで中〜高Counts領域と同じ連続的な分散傾向の延長として扱うと、 シグナルよりもノイズの影響が大きい遺伝子まで有意判定の対象に含まれることがあります。（Fig1 左・中）
つまり、edgeRやDESeq2を使う場合でも、 低Counts領域をそのまま解釈対象に含めるわけではありません。 実際の解析では、低発現遺伝子を除外するフィルタリングを併用することが一般的です。 このことからも、低Counts領域は、統計モデルに任せれば自動的に安全に解釈できる領域ではなく、 解析前にどこまでを対象に含めるかを判断すべき領域だと考えられます。
Fig1: edgeR、DESeq2、t検定で p < 0.05 となった遺伝子の分布比較。 黒い点は全遺伝子を示し、色付きの点はそれぞれの方法で有意と判定された遺伝子を示します。 3つの方法で、有意と判定される遺伝子の大まかな分布は似ていますが、低発現領域や判定の境界付近では違いが見られます。
edgeRとDESeq2の差は、主に境界領域に見られる
Fig2 を見ると、3つの手法で重複する遺伝子は多い一方で、 手法に依存して有意判定されたり、されなかったりする遺伝子もあることが分かります。
edgeRとDESeq2の差は、主に有意・非有意の判定が分かれる境界領域に見られます。(Fig2 左) 明らかに差が大きい遺伝子では両者の結果はよく一致しますが、 p値が閾値付近にある遺伝子では、どちらの手法で有意と判定されるかが分かれます。 また、このデータでは、edgeRの方が低Counts領域で有意判定しやすい傾向もみられます。
少サンプルのRNA-Seqデータでは、edgeRがDESeq2より多くの候補遺伝子を検出する傾向を示すことがあると報告されています。 今回の比較でも、edgeRの方がDESeq2より多くの遺伝子を有意と判定しており、 そのような傾向と一致する結果になりました。 ただし、これは常にedgeRの方が多く検出するという意味ではありません。
Fig2: edgeR、DESeq2、t検定で p < 0.05 となった遺伝子の比較。 下段は3つの方法で検出された遺伝子の重なりを示します。下段のベン図で、灰色のオーバーレイされた領域の遺伝子群を、上段の散布図の黒い点で表しています。
edgeRやDESeq2のみで有意判定された遺伝子は、平均値の差が比較的大きいものが多い
edgeRとDESeq2でのみ有意と判定された遺伝子群は、t検定のみで有意と判定された遺伝子群と比べると、 対角線から少し離れた位置にあります。（Fig2 中） これらの遺伝子群は、群内のばらつきが比較的大きく見える遺伝子も含まれています。（Fig2B 中） これは、edgeRやDESeq2が各遺伝子の4サンプルだけのばらつきではなく、 同程度のGene Countsを持つ多数の遺伝子から推定される分散の傾向を利用しているためと考えられます。 このデータセットはもともと群内ばらつきが小さいため、 同程度のGene Countsを持つ遺伝子全体から推定される分散傾向も小さくなりやすいと考えられます。 その結果、個別遺伝子の測定値ではばらつきが大きく見える遺伝子でも、 edgeRやDESeq2のモデル上では有意と判定されたと考えられます。
t検定のみで有意判定された遺伝子は、平均値の差が小さいものが多い
一方で、t検定のみで有意判定された遺伝子は、対角線に非常に近いものが多いです。（Fig2 右） そして、これらの遺伝子は分散が極端に小さいことが分かります。（Fig2B 右）
繰り返し数が2回と極端に少ない場合、測定される値は母集団の確率が高い領域に集まりやすく、 母集団の分散よりもかなり小さな分散として観測される傾向があると考えられます。 その結果、平均値の差が小さくても、各群内のばらつきが非常に小さく観測された場合には有意と判定されやすいことを示しています。 p値だけで判断する限り、この結果は、RNA-Seqデータにはt検定ではなく、edgeRやDESeq2を用いるべきだという主張に整合性を与えるでしょう。
Fig2B: 3つの方法で共通して有意と判定された遺伝子、edgeRとDESeq2でのみ有意と判定された遺伝子、t検定でのみ有意と判定された遺伝子の発現パターン比較。 下段のベン図で、灰色のオーバーレイされた領域の遺伝子群を、上段の線グラフの黒い線で表しています。
前処理・フィルタリングは、低Counts領域のリスク管理である
前処理やフィルタリングを行わずにt検定を適用すると、 低Counts領域の遺伝子が有意と判定されてDEGリストに含まれてしまいます。 この領域では、カウントの離散性、欠損値、低シグナル由来の不安定さの影響を受けやすく、 p値だけで発現差を解釈するには注意が必要です。
これを確かめるため、正規化とlog2変換のみ行った場合に有意判定された遺伝子群と、 Low Signal Cutoff、欠損値の補完、フィルタリングを施したデータで有意判定された遺伝子群を比較しました。
Low Signal Cutoff、欠損値の補完、フィルタリングを行った後にt検定を適用すると、 低Gene Counts領域で有意と判定される遺伝子は減少します。 この結果は、t検定そのものよりも、 低カウント領域の不安定なデータをそのまま検定に入れることが大きなリスクになることを示しています。
edgeRやDESeq2は、RNA-SeqのCountデータに合わせた統計モデルを自動的に適用します。 これは非常に簡便ですが、統計モデル自体の中で、 解析者がどのCounts領域をどの程度まで解釈対象に含めるかを、 生物学的な目的に応じて調整する余地は限られます。 Fig1の左や中の散布図を見ると、２回の繰り返ししかない状況としては大胆に有意判定しているように見えます。 edgeRやDESeq2を使う場合でも、低発現遺伝子に対する何らかのフィルタリングを併用する方が適切だと考えられます。
これに対して、前処理・フィルタリングをデータを見ながら手作業で行う場合は、 解析目的に合わせてLow Signal Cutoffの閾値を柔軟に変更できます。 探索的解析ではやや低めに、厳密なバイオマーカー候補選定では高めにするなど、 解析目的に応じて判断を変えられることには、生物学的解釈上の意味があります。
Fig3: t検定における前処理とフィルタリングの影響。 左は正規化とlog2変換のみを行った後にt検定を適用した結果を示し、中央はさらにLow Signal Cutoff、欠損値の補完、フィルタリングを行った後にt検定を適用した結果を示します。 右は、前処理とフィルタリングによってt検定の有意遺伝子から除外された遺伝子の分布を示します。 これらの除外された遺伝子の多くは低Gene Counts領域に分布しており、前処理とフィルタリングによって、低カウント領域の不安定な遺伝子を有意と判定するリスクを減らせることがわかります。
p値だけでなくfold changeも加えると、実用的なDEGリストに近づく
実際の発現差解析では、p値だけでDEGリストを作成するとは限りません。 fold changeの大きさを条件に加え、生物学的に解釈しやすい変化に絞り込むことが一般的です。 そこでここでは、edgeR、DESeq2、そして前処理・正規化・フィルタリング後のデータに対するt検定について、 p < 0.05 かつ 1.4倍以上の発現差を示す遺伝子を比較しました。
p値だけで判定した場合には、平均値の差が小さい遺伝子も有意と判定されることがあります。 特にt検定では、群内のばらつきが非常に小さい遺伝子では、 fold changeが小さくても有意になりやすい傾向があります。
fold changeの条件を加えると、対角線に近い小さな変動の遺伝子は除外され、 有意遺伝子はより発現差の大きい領域に絞られます。
Fig4: edgeR、DESeq2、t検定で p < 0.05 かつ fold change 1.4倍以上となった遺伝子の分布比較。
fold change条件を加えると、手法間の違いは小さくなる
Fig5では、p < 0.05 に加えて fold change 1.4倍以上という条件を加えた場合に、 3つの方法で得られるDEGリストがどのように変化するかを比較しています。
p値だけで比較した場合には、t検定のみで有意と判定される遺伝子が355個ありました。 これらの遺伝子には、群内のばらつきが非常に小さい一方で、 平均値の変動幅は小さいものが多く含まれていました。 また、edgeRやDESeq2のみで有意判定される遺伝子にも、 有意判定の境界付近に、変動幅の小さい遺伝子が多く含まれていました。
fold change 1.4倍以上という条件を加えると、 これら変動幅の小さな遺伝子群が大量に除かれます。 その結果、edgeR、DESeq2、t検定で得られるDEGリストの重なりは大きくなり、 手法間の違いはp値だけで比較した場合よりも小さくなります。（Fig5）
RNA-Seqの発現変動解析では、p値だけでDEGリストを作るのではなく、 fold changeを組み合わせることが多いです。 したがってこの結果は、現実的なRNA-Seqの発現変動解析においては 手法間の差が小さくなっていることを示しています。
Fig5: p < 0.05 かつ fold change 1.4倍以上という条件で抽出したDEGリストの比較。 fold change条件を加えることで、p値だけで比較した場合に見られた手法間の違いは小さくなり、3つの方法で得られるDEGリストの重なりが大きくなります。
ただし、fold change条件を加えても手法間の違いは完全には消えない
fold change 1.4倍以上という条件を加えることで、 t検定のみで有意と判定されていた小さな変動の遺伝子は大きく減少します。 一方で、edgeRやDESeq2のみで有意と判定された遺伝子群は、まだ多く残っています。（Fig5）
これらの遺伝子は各群内のばらつきが比較的大きいため、t検定では有意と判定されなかったものと考えられます。（Fig5B 左・中） 一方、edgeRやDESeq2ではGene Counts全体から推定される分散傾向も利用するため、 個別遺伝子のプロファイルだけを見るとばらつきが大きく見える遺伝子でも、有意と判定されることがあります。
この結果は、fold change条件を加えても、手法間の違いが完全になくなるわけではないことを示しています。 edgeRやDESeq2で得られたDEGリストであっても、個別遺伝子の発現パターンを確認する必要があります。
Fig5B: p < 0.05 かつ fold change 1.4倍以上で抽出した遺伝子群の発現パターン比較。 fold change条件を加えた後でも、edgeRやDESeq2のみで有意と判定された遺伝子群には、個別サンプルで見ると群内ばらつきが比較的大きいものが残っています。
この結果は、低分散・少サンプルのin vitroデータに強く依存している
ただし、今回使用したのはin vitroの実験データであり、 繰り返しサンプル間のばらつきが非常に小さいデータです。 さらに、各条件の繰り返し数は2回であり、 統計的には極端にサンプル数が少ない比較です。
このようなin vitro実験で、各群の繰り返し数が非常に少ないデータは、 実際のRNA-Seq研究でも少なくありません。 低分散・少サンプルのデータでは、Gene Counts全体から推定される分散傾向が小さくなりやすく、 edgeRやDESeq2で、t検定では有意にならない遺伝子まで有意と判定されることがあります。 これは、このタイプのデータでedgeRやDESeq2を使う際に知っておくべき特徴の一つです。
そこで次回は、biopsyサンプルを用いた中規模のRNA-Seqデータで、 同様にedgeR、DESeq2、前処理済みデータに対するt検定の結果を比較し、 今回見られた傾向がどの程度再現されるのかを検証してみたいと思います。
重要なのは、どの手法を信じるかではなく、データを見て判断すること
今回の比較では、edgeR、DESeq2、t検定の結果は大きく矛盾しているわけではありません。 p値だけで比較すると手法ごとの差は大きく見えますが、 fold change条件を加えると、実用的なDEG候補リストとしての差はかなり小さくなりました。
一方で、手法ごとの差分に含まれる遺伝子を個別に確認すると、 それぞれ異なる特徴が見えてきます。 t検定のみで有意と判定された遺伝子には、 群内ばらつきが非常に小さい一方で、平均値の差は小さいものが多く含まれていました。 このような遺伝子は、fold change条件を加えることで大きく減少します。
逆に、edgeRやDESeq2で有意、t検定では非有意だった遺伝子には、 個別サンプルで見ると群内ばらつきが比較的大きいものが含まれていました。 低分散・少サンプルのin vitroデータでedgeRやDESeq2を使う場合には、 このような遺伝子がDEGリストに含まれやすいことを知ったうえで、 慎重に解釈する必要があります。
Gene Countsを用いたRNA-Seqの発現変動解析では、 edgeRやDESeq2のようにカウントデータ全体の分散構造を考慮できる方法を基本にするのが安全です。 一方で、t検定も、正規化・前処理・フィルタリング・fold change条件・可視化確認を組み合わせれば、 edgeRやDESeq2と大きく矛盾しない、実用的な発現差評価の方法として利用できます。
統計手法は、候補遺伝子を絞り込むための重要な道具です。 しかし、どの遺伝子を解釈対象として残すかを最終的に判断するには、 p値だけでなく、実際のデータを確認する必要があります。
Subio PlatformでDEGリストの違いを可視化して確認する
RNA-Seqデータ解析チュートリアル では、データのインポートから、正規化、フィルタリング、そしてDEG解析までの手順を紹介するだけでなく、 それぞれのステップで何を見て、どのように判断していくかについて、詳細に解説しています。
DEG解析は、Gene CountsをedgeRかDESeq2に取り込んで処理させればいい、と思っている方は多いと思います。 しかし、それでは見落とされる重要なことは多いのです。
また、チュートリアルページには、SSAファイルのダウンロードリンクもあります。 SSAファイルは、Subio Platformにインポートすることで、 実際にそのデータをご自身の手で触ってみることができるファイル共有の仕組みです。 今回の検証も、ぜひご自身で手を動かして確認してみてください。