RNA-Seqの発現差解析では、edgeRやDESeq2のようなRNA-Seq専用の統計手法がよく使われます。 これらの手法は、Gene Countsデータの特徴に合わせて特別に設計された統計モデルとなっています。
一方で、 「edgeRを使うべき」 「DESeq2を使うべき」 「RNA-Seqでは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が、すべての群において、１サンプル以上で20未満の遺伝子を除去
p値だけで見ると、手法ごとの違いは大きく見える
まず、edgeR、DESeq2、t検定について、 それぞれ p < 0.05 となった遺伝子を比較しました。
Fig1 の散布図では、siC群とsiE1群の繰り返しサンプルの平均値を比較しています。 そのため、個別サンプル同士を比較した場合とは見え方が異なります。
個別サンプル同士を比較した散布図だと、低Counts領域では、わずかな測定値の差が 大きなFold Changeになりやすく、対角線から広がって、ばらつきが大きく見えます。
一方で、Fig1 の散布図は、siC群とsiE1群の繰り返しサンプルの平均値を比較しています。 群平均同士の比較では、個別サンプルのランダムなばらつきがある程度平均化されるため、 低Counts領域の点は、個別サンプル同士を比較した場合のようには広がって見えません。 今回のデータは各群2サンプルの平均値なので、ばらつきが強く収束するわけではありませんが、 サンプル数が増えるほどランダムなばらつきは平均化され、 低Counts領域の点はより対角線付近に集まりやすくなると考えられます。
edgeRやDESeq2は、Countデータの平均値と分散の関係を考慮した統計モデルによって有意差を判定します。 一方、t検定では、群内のばらつきと群間差を評価します。 そのため、同じデータを使っていても、有意と判定される遺伝子は異なります。Fig2 を見ると、３つの手法で重複する遺伝子が多かったですが、手法に依存して有意判定されたり、されなかったりする遺伝子もあることが分かりました。
Fig1: edgeR、DESeq2、t検定で p < 0.05 となった遺伝子の分布比較。 黒い点は、それぞれの方法で有意と判定された遺伝子を示します。
edgeRとDESeq2の差は、主に境界領域に見られる
edgeRとDESeq2の差は、主に有意・非有意の判定が分かれる境界領域に見られます。 明らかに差が大きい遺伝子では両者の結果はよく一致しますが、 p値が閾値付近にある遺伝子では、 どちらの手法で有意と判定されるかが分かれます。（Fig2 左）
少サンプルのRNA-Seqデータでは、 edgeRがDESeq2より多くの候補遺伝子を検出する傾向を示すことがあると報告されています。 今回の比較でも、edgeRの方がDESeq2より多くの遺伝子を有意と判定しており、 そのような傾向と一致する結果になりました。ただし、その差は小さいです。 また、これは常にedgeRの方が多く検出するという意味ではありません。
Fig2: edgeR、DESeq2、t検定で p < 0.05 となった遺伝子の差分比較。 edgeRとDESeq2の差は主に判定の境界領域に見られ、 t検定のみで有意と判定された遺伝子には、対角線から大きく離れたものも含まれていました。
edgeRとDESeq2の両方で有意、t検定では非有意と判定された遺伝子
今回の結果で、edgeRとDESeq2の両方で有意と判定され、 t検定では有意と判定されなかった遺伝子を確認しました。（Fig2 中央）
edgeRやDESeq2では、個々の遺伝子で観測されたばらつきだけを使って有意差を判定するのではなく、 同程度のCounts値を持つ遺伝子全体のばらつきも利用して分散を推定します。 そのため、Countsが高い領域の遺伝子は、 その遺伝子の実際の観測値にある程度ばらつきがあっても、 有意と判定されやすくなります。
これに対し、t検定では、各遺伝子について実際に観測された群内ばらつきと群間差を直接評価します。 そのため、単純に繰り返しサンプル間のばらつきが大きければ、有意と判定されにくくなります。
この違いを頭に入れて、Fig2B をご覧ください。Fig2B を見ると、3つすべてで有意判定された遺伝子群（左）とくらべて、edgeRとDESeq2だけで有意判定された遺伝子群は、ばらつきが大きいという特徴がみられました（右）。今回のデータはin vitro実験で全体としてばらつきが非常に小さいため、 RNA-Seq専用モデルでは有意判定を通過しやすかったと考えられます。
Fig2B: 3手法すべてで有意判定された遺伝子群と、edgeRまたはDESeq2側に偏って有意判定された遺伝子群の比較。 3手法で共通して検出された遺伝子群に比べて、edgeRやDESeq2側に偏る遺伝子群では、 群内のばらつきが大きいものが含まれていました。
統計学的に考えると、 繰り返しサンプル数が少ないときは、偶然ばらつきが小さく見える可能性が高くなります。 言い換えると、少ない繰り返しの中でも大きなばらつきが観測されるなら、 それは外れ値をたまたま観測した可能性より、その遺伝子のばらつきが非常に大きいという特徴を表している可能性の方が高いということになります。
このような個別の遺伝子の観測値のばらつきより、同程度のCounts値を持つ遺伝子群全体のばらつきに依存するというedgeRやDESeq2のモデル設計が、生物学的解釈において必ずしも望ましい方向に働くとは限らないかもしれません。
t検定では有意、edgeR・DESeq2では非有意だった遺伝子
次に、t検定では有意と判定され、 edgeRとDESeq2の両方では有意と判定されなかった遺伝子を確認しました。（Fig3）
非常におもしろいのは、P値の閾値付近（散布図でいうと、対角線に近いところ）の遺伝子群と、大きく発現変動している（散布図でいうと対角線から離れている）遺伝子群という、二つのグループがあるということです。おそらく対角線に近いグループの方は、上述のとおり測定値の実際の分散によるものと考えられます。ここでは、対角線から離れた方（Fig3 左で黒く示した）遺伝子群について検証します。
Line Graphで個別に確認すると、 それぞれのグループ内で再現性が高く、 かつグループ間では大きく離れている遺伝子でした。 t検定のみで有意と判定された遺伝子の一部は、 低品質なばらつきの大きい遺伝子ではなく、 生物学的に確認すべき候補として残す価値がある遺伝子を含んでいることが分かりました。このような発現パターンであれば、 t検定で非常に低いP値をもって有意と判定されるのは自然で、実際そうなっています。
これらの遺伝子がなぜedgeRおよびDESeq2では有意判定から外れたのかを考えます。 DESeq2には、Cook's distanceに基づいて極端なcount outlierを検出し、 その遺伝子のP値をNAにしたり、 条件によってoutlier countを置換して再フィットしたりする仕組みがあります。 edgeRでも、外れ値や過剰なばらつきの影響を抑えるためのrobustな推定や quasi-likelihoodに基づく判定が使われます。
そのため、個別に見ると群内の繰り返し測定が安定し、 群間差が大きく見える遺伝子であっても、 RNA-Seq専用モデル上では、外れ値的な発現パターンや分散構造として扱われ、 有意判定から外れることがあります。 今回のデータでも、このようなモデル上の扱いが影響した可能性があります。
DEG解析にedgeRやDESeq2を使ううえで、 この特徴は特に注意して理解しておくべきポイントです。
Fig3: t検定では有意と判定され、edgeRおよびDESeq2では有意判定から外れた遺伝子の例。 散布図上で対角線から大きく離れ、Line Graphでも群内の再現性が高く、 群間差が大きい遺伝子が含まれています。
前処理・フィルタリングは、低Counts領域のリスク管理である
前処理やフィルタリングを行わずにt検定を適用すると、 低Counts領域の遺伝子が有意と判定されてDEGリストに紛れ込んでしまいます。一方で、正規化後にLow Signal Cutoff未満の値を閾値まで引き上げる処理を行い、 さらに常に低Countsの測定値ばかりで解析に資すると考えにくい遺伝子をフィルタリングすれば、 このような低Counts領域の候補を大きく減らすことができると考えられます。これを確かめるため、正規化と対数変換のみ行った場合に有意判定された遺伝子群と、適切な前処理とフィルタリングを施したデータで有意差判定された遺伝子群を比較しました。(Fig4)
処理前後のDEGリストの差分を見ると、除かれた遺伝子は低Counts領域に集中しており、 前処理とフィルタリングが、低Counts領域の信頼性の低さに対するリスク管理として機能していることが分かります。
edgeRやDESeq2は、RNA-SeqのCountデータに合わせた統計モデルを自動的に適用します。 これは非常に簡便ですが、統計モデル自体の中で、 解析者がどのCounts領域をどの程度まで解釈対象に含めるかを、 生物学的な目的に応じて調整する余地は限られます。
これに対して、前処理・フィルタリングをデータを見ながら手作業で行う場合は、 解析目的に合わせてLow Signal Cutoffの閾値を柔軟に変更できます。 探索的解析ではやや低めに、厳密なバイオマーカー候補選定では高めにするなど、 解析目的に応じて判断を変えられることには、生物学的解釈上の意味があります。
Fig4: 前処理・フィルタリングなしのt検定と、前処理・フィルタリング後のt検定の比較。 処理によってDEGリストから除かれた遺伝子は低Counts領域に集中しており、 低Counts領域の不安定な候補を減らすリスク管理として機能していることが分かります。
p値だけでなくFold Changeも加えると、実用的なDEGリストに近づく
実際の発現差解析では、p値だけでDEGリストを作成するとは限りません。 Fold Changeの大きさを条件に加え、 生物学的に解釈しやすい変化に絞り込むことが一般的です。そこでここでは、edgeR、DESeq2、そして前処理・正規化・フィルタリング後のデータに対するt検定について、 p < 0.05 かつ 1.4倍以上の発現差を示す遺伝子を比較しました。（Fig5）
Fig5: edgeR、DESeq2、t検定で、p < 0.05 かつ1.4倍以上の発現差を示した遺伝子の比較。 p値だけでなくFold Change条件を加えることで、 実用的なDEG候補に近い形で手法間の違いを確認しています。
その結果、p値だけで比較した場合に比べて、 3つの手法で共通して検出される遺伝子の割合が多くなり、 手法間の違いはきわめて小さくなりました（Fig6 ベン図）。これは、Fig2で手法による違いが対角線付近の遺伝子に集中していたことから、これらの遺伝子が除外された結果として考えれば、きわめて整合的です。
t検定のみで有意と判定された遺伝子には、 散布図上で対角線から大きく離れた外側の遺伝子があることは、前述のとおりです。これらの遺伝子群は、Fold Changeを条件に加えても影響を受けませんので、そのまま残っていました。 これらの遺伝子は、Fold Changeが十分に大きく、群間差が明瞭に見える遺伝子群です。（Fig6 右）
一方で、edgeRやDESeq2で有意、t検定では非有意だった遺伝子にも特徴が見られました。 Line Graphで確認すると、これらの遺伝子には、 群内の繰り返しサンプル間のばらつきが比較的大きいものが含まれていました。（Fig6 左）
この結果は、 「RNA-Seqではt検定を使ってはいけない」という単純な理解では不十分であることを示しています。 一般的には、「edgeRやDESeq2の方が保守的な結果を出すから、t検定よりこれらの手法を使うべきだと言われている」というイメージが先行していると思います。しかし、今回の検証結果は真逆で、edgeRやDESeq2の方がリベラルな結果だったのです。
Fig6: p < 0.05 かつ1.4倍以上の発現差を示す遺伝子に絞ったうえで、 edgeR/DESeq2側で検出されやすい遺伝子と、t検定側で検出されやすい遺伝子を比較しました。 左はedgeRまたはDESeq2で有意と判定され、t検定では有意と判定されなかった遺伝子、 右はt検定で有意と判定され、edgeRおよびDESeq2では有意と判定されなかった遺伝子を示します。
ただし、今回使用したのはin vitroの実験データであり、 繰り返しサンプル間のばらつきが非常に小さいデータです。 さらに、各条件の繰り返し数は2回であり、 統計的には極端にサンプル数が少ない比較です。 この条件が今回の結果に大きく影響していると考えられるため、 biopsyサンプルや臨床検体のように、 ばらつきが大きくN数も多いデータでは、 異なる結果になる可能性があります。
したがって、次回はbiopsyサンプルを用いた中規模のRNA-Seqデータで、 同様にedgeR、DESeq2、前処理済みデータに対するt検定の結果を比較し、 今回見られた傾向がどの程度再現されるのかを検証してみたいと思います。
一方で、このようなin vitro実験で、各群の繰り返し数が非常に少ないデータは、 実際のRNA-Seq研究でも少なくありません。 今回の結果だけを見ると、このような低ばらつき・少サンプルのin vitroデータでは、 edgeRやDESeq2よりも、前処理済みデータに対するt検定の方が、 生物学的に納得しやすい候補を拾えているように見えます。
つまり、そのようなケースでは、 少なくともedgeRやDESeq2だけに頼るのではなく、 t検定も行ってみる価値が高いと思われます。 そのうえで、それぞれの違いを比べながら候補遺伝子を判断するとよいでしょう。
重要なのは、どの手法を信じるかではなく、データを見て判断すること
今回の比較では、edgeR、DESeq2、t検定の結果は大きく矛盾しているわけではありません。 p値だけで比較すると手法ごとの差は大きく見えますが、 Fold Change条件を加えると、実用的なDEG候補リストとしての差はかなり小さくなりました。
一方で、手法ごとの差分に含まれる遺伝子を個別に確認すると、 それぞれ異なる特徴が見えてきます。 edgeRやDESeq2で有意、t検定では非有意だった遺伝子には、群内ばらつきが比較的大きいものが含まれていました。 逆に、t検定では有意、edgeRやDESeq2では非有意だった遺伝子には、 群内の再現性が高く、群間差も大きいものが含まれていました。
この結果は、 「edgeRが正しい」 「DESeq2が正しい」 「t検定が正しい」 という単純な話ではありません。 統計モデルの違いによって、候補遺伝子の拾われ方は変わります。RNA-Seqの発現差解析では、 手法名だけで結果を判断するのではなく、 p値、Fold Change、Countsの大きさ、群内のばらつき、 フィルタリング条件、そして個々の遺伝子の発現パターンを合わせて確認することが重要です。 統計手法は、候補遺伝子を絞り込むための重要な道具です。 しかし、どの遺伝子を解釈対象として残すかを最終的に判断するには、 p値だけでなく、実際のデータを確認する必要があります。
Subio PlatformでDEGリストの違いを可視化して確認する
RNA-Seqデータ解析チュートリアル では、データのインポートから、正規化、フィルタリング、そしてDEG解析までの手順を紹介するだけでなく、それぞれのステップで何を見て、どのように判断していくかについて、詳細に解説しています。
DEG解析は、Gene CountsをedgeRかDESeq2に取り込んで処理させればいい、と思っている方は多いと思いますが、それでは見落とされる重要なことは多いのです。この機会に、よろしければチュートリアルをお試しください。
また、チュートリアルページには、SSAファイルのダウンロードリンクもあります。SSAファイルは、Subio Platformにインポートすることで、実際にそのデータをご自身の手で触ってみることができるファイル共有の仕組みです。今回の検証も、ぜひご自身で手を動かして確認してみてください。