RNA-Seqの発現差解析では、edgeRやDESeq2のようなGene Countsに対応した統計手法が推奨されています。 これらの手法は、RNA-Seq Countsの性質を考慮して設計されており、 「少数のサンプルでも信頼性の高い解析ができる」といわれています。
しかし、マイクロアレイの黎明期からこのオミクス業界では、「みんなが使う手法」を盲目的に使うという悪しき伝統があります。 「少数のサンプルでも信頼性の高い解析ができる」ようにするための工夫には高い価値があるのですが、万能薬のような魔法はやはり存在しないのです。
現実のRNA-Seqのデータは、 in vitro実験かbiopsy由来データか、繰り返し数が多いか少ないか、インプットRNA量が多いか少ないか、paired designかどうかなどによって、データの性質が大きく異なります。当然ながら、データの性質が異なれば、ある手法を使うことでうまくいくこともあれば、うまくいかないこともあるでしょう。
本ページでは、edgeR、DESeq2、t検定の原理的な違いを整理したうえで、 代表的な4つの実験データタイプにおいて、 それぞれの手法でどのような遺伝子が有意と判定されやすいのかを見ていくことで、これらの手法をどのように使い分ければよいかを考えます。
t検定とedgeR・DESeq2の基本的な違い
t検定とedgeR・DESeq2の大きな違いは、 分散をどのように扱うかにあります。
ただし、edgeRやDESeq2について一般的に説明されている内容や、 そこから受けるイメージと、実際の計算や結果との間には、 3つのギャップがあります。
t検定は、各遺伝子について実際に観測されたサンプル間のばらつきに基づいて、 群間差が大きいかどうかを評価します。 つまり、各遺伝子で観測された群内分散を、 その遺伝子の検定に直接使用します。
一方、edgeRやDESeq2では、 各遺伝子で観測された分散がそのまま検定に使われるわけではありません。 同程度のGene Countsを持つ多数の遺伝子のばらつきを利用して、 RNA-Seqデータ全体における平均値と分散の関係を推定し、 個々の遺伝子の分散をその傾向へ近づけます。 この処理は、一般にdispersion shrinkageと呼ばれます。
一般的なRNA-Seq発現差解析の説明では、 edgeRやDESeq2の特徴として、 Gene Countsを負の二項分布でモデル化することが強調されます。 確かに、低Count領域では、 あるサンプルでは測定され、別のサンプルでは測定されないという検出の偶然性や、 平均値に対して大きなばらつきが生じやすく、 Countデータに適したモデルを使うことには意味があります。
しかし、標準的な解析では、 このような低発現遺伝子の多くがフィルタリングによって除かれます。 そのため、フィルタリング後のデータでは、 負の二項分布を採用したこと自体よりも、 edgeRやDESeq2が行うdispersion shrinkageの方が、 結果に大きな影響を与えると考えられます。
一般的な説明では負の二項分布が強調されますが、 実際の計算を大きく左右するdispersion shrinkageの影響は、 それほど明確には説明されていません。 これが、一般的な説明と現実の第一のギャップです。
dispersion shrinkageでは、 個別の遺伝子で観測された分散が、 同程度のGene Countsを持つ他の遺伝子より大きい場合、 検定に使用される分散が小さい方向へ補正されることがあります。 その結果、実際に観測された分散を直接使うt検定と比べて、 edgeRやDESeq2の方が検定を通過しやすくなる遺伝子があります。
これは、edgeRやDESeq2はt検定よりも保守的であるという 一般的なイメージとは異なります。 観測された分散が小さい遺伝子では、 分散が大きい方向へ補正されて検定を通過しにくくなる一方、 観測された分散が大きい遺伝子では、 分散が小さい方向へ補正され、 t検定よりも検定を通過しやすくなる場合があります。
したがって、edgeRやDESeq2が、 個々の遺伝子について常にt検定より保守的に判定するわけではありません。 これが第二のギャップです。
さらに、dispersion shrinkageは、 反復数が極端に少なく、 個々の遺伝子の分散を安定して見積もれない場合だけに行われる処理ではありません。 ある程度のサンプル数があり、 個々の遺伝子で観測された分散の信頼性が高くなったと考えられる場合でも、 この処理は行われます。
サンプル数が増えれば、 個々の遺伝子で観測された分散の重みは大きくなります。 しかし、同程度のGene Countsを持つ他の遺伝子から推定された分散傾向も、 引き続き検定結果に影響します。
そのため、この処理を単に 「遺伝子間の情報共有による分散推定の安定化」 と説明するだけでは、 実際に何が行われ、 結果がどのように変わる可能性があるのかを十分に表していません。
個々の遺伝子で観測された大きな分散が、 その遺伝子の本当の生物学的な不均一性を反映している場合でも、 その分散が小さい方向へ補正される可能性があります。 このような場合まで含めて、 この処理を単純に「分散推定の安定化」と呼ぶことが適切なのかは、 慎重に考える必要があります。 これが第三のギャップです。
ただし、ここまで述べたdispersion shrinkageの特徴とは別に、 edgeRやDESeq2がRNA-Seq発現差解析にもたらした重要な改善もあります。 それは、TPMやFPKMではなく、 Gene Countsを使って発現差を評価する流れを定着させたことです。 これは、edgeRとDESeq2の非常に大きな功績だと思います。
TPMやFPKMは、遺伝子長で正規化することにより、 Gene Countsが持つ平均値と分散の関係を変えてしまいます。 これは、RNA-Seqに限らず、 実験測定データの解析で蓄積されてきた理論とノウハウをなげうつことに等しく、 信頼性の低いDEG解析の原因となっていました。 TPM、FPKM、Gene Countsの違いについては、 こちらの記事で詳しく説明しています。
以上の特徴を踏まえて、ここからは代表的な4つのデータタイプについて、 edgeR、DESeq2、t検定の結果がどのように異なり、 どのように使い分けられるのかを見ていきます。
ケース1：低分散・少数サンプルのin vitroデータ
このケースの詳しい検証は、 低分散・少数サンプルのin vitroデータでedgeR・DESeq2・t検定を比較した事例 で紹介しています。
in vitroの実験では、個体差に由来するばらつきがないため、 群内分散が非常に小さくなりやすいです。 さらに、繰り返し数が2など極端に少ない場合には、 たまたま2つの値が近いだけで、 その遺伝子の分散は小さく見積もられがちになることが予想されます。
もちろんこのような実験デザインは、 統計学的モデルの想定とはかけ離れています。 しかし、現実のRNA-Seqデータではよく見られる条件でもあります。 そのため、このような場合にそれぞれの手法でどのような遺伝子が有意と判定されやすいのかを理解しておく必要があります。
実際に計算してみた結果、測定される群内分散が非常に小さいことで、いずれの手法を用いても、 変動幅が極めて小さい遺伝子まで有意と判定されました。
そのため、このようなケースでは、p値だけでDEGを抽出するのではなく、 fold change条件を組み合わせることが推奨されます。 これにより、 わずかな変動幅で有意になっていた遺伝子を除けます。そして、fold change条件を組み合わせると、3手法で得られるDEGリストの違いはかなり小さくなります。
ただし、fold change条件を加えた後でも、 手法間の違いは残ります。
t検定のみで有意となる遺伝子は、 観測された群内分散が極端に小さい傾向 があります。 このような遺伝子では、 たまたま繰り返し間の値がそろったことで、 t検定上は小さなp値になりやすくなります。ただ、この種のエラーは、fold change条件と組み合わせることで簡単に対処できます。そしてこれは、edgeRやDESeq2を使う場合でも同様に効きます。
次に、edgeRやDESeq2では有意となる遺伝子は、 観測された群内分散が比較的大きい傾向があります。これは、その遺伝子の実際に観測された分散ではなく、同程度のCountsを持つ他の遺伝子のばらつきも利用して分散が推定されることに由来します。in vitroの低分散データでは、観測されるばらつきが小さい遺伝子が多いことで、結果として、t検定よりも甘めに有意と判定される遺伝子が出てきます。
ケース2：biopsy由来の中規模データ
このケースの詳しい検証は、 biopsy由来の中規模データでedgeR・DESeq2・t検定を比較した事例 で紹介しています。
次に、biopsy由来の中規模RNA-Seqデータを考えます。 このようなデータでは、ケース1と比べてサンプル数が多く、 発現量が十分に高いシグナル領域では、 分散の見積もりが比較的安定します。そのため、発現量が十分に高い領域では、 edgeR、DESeq2、t検定の結果の違いは小さく、 どの手法を用いても主要な発現変動遺伝子は共通して検出できると考えていいでしょう。
一方で、このケースでは発現量が低い領域において、 edgeRやDESeq2のみで有意と判定される遺伝子が多数見られました。 このような遺伝子を詳しく確認すると、 明確な生物学的発現差というより、 サンプル間のダイナミックレンジの違いに由来する 正規化後Gene Countsの下限値のズレと対応している可能性が高い場合があることが分かりました。
RNA-Seqデータでは、サンプル間で総リード数に数倍の差があることは珍しくありません。 総リード数の差とダイナミックレンジの差は単純に比例するわけではありませんが、ダイナミックレンジの差が正規化後のGene Countsの下限値のズレを引き起こすことは、ある程度の規模のRNA-Seq実験データでは常に見られる問題です。そして、この下限値のズレが一方の群に偏ると、 群間差のように見える原因になります。
edgeRやDESeq2の統計モデルは、 低発現領域のCountsが不安定になりやすいことを考慮して設計されています。 しかし、サンプル間のダイナミックレンジの違いによって、 正規化後Gene Countsの下限値がずれることはモデルに組み込まれていません。 したがって、低Count領域でedgeRやDESeq2のみ有意となる遺伝子については、 p値だけで判断するのは危険です。
ただし、edgeRやDESeq2の標準的な解析プロトコルでは、 検定前に低発現遺伝子のフィルタリングなどの処理が入るため、この問題が表面化するのは、ダイナミックレンジの差が極端に大きい場合に限られます。 とはいえ、このようなRNA-Seqデータの性質は常に意識しておく必要があります。さもなければ、ダイナミックレンジの差に由来するDEG候補をそのまま信用してしまうリスクがあります。
edgeRやDESeq2と同様に、t検定を用いる場合も適切な前処理とフィルタリングは必須です。低Count領域の測定値を適切に処理していれば 、3つの手法間の差は小さくなりました。 つまり、どの手法を選ぶかだけでなく、 検定前に適切な前処理とフィルタリングを行うことが、 信頼性の高いDEG解析にとって重要だということです。
ただし、この場合は低Count領域でON/OFFの動きをしている遺伝子を見逃してしまいがちです。これらの遺伝子については、測定値の有無を条件に抽出することで、低Count領域の不安定性やダイナミックレンジの差を受けてしまうp値に頼らずに、安全に拾うことができます。
一方で、ケース1と異なり、 ある程度サンプル数がある中規模データでは、 fold change条件を必ずしも一律に加える必要はありません。 発現量が十分に高い領域で、 低変動ながら一貫して有意となる遺伝子は、 生物学的な小さな変化を反映している可能性があります。ただし、低変動DEGの解釈は、大きな変動を示す遺伝子よりは慎重に行うことが推奨されます。
ケース3：paired designのデータ
このケースの詳しい検証は、 paired designでedgeR・DESeq2・paired t-testを比較した事例 で紹介しています。
3つ目のケースは、paired designのRNA-Seqデータです。 たとえば、同じ患者からcontrolとcaseのサンプルが得られている場合、paired t-testでは、個体差をキャンセルしたうえで、 各ペア内でcontrolからcaseへどの方向にどれだけ変化したかを評価します。 そのため、対応のない2群比較よりも、paired t-testの方が一般的には検出力が高くなります。
一方、edgeRやDESeq2でも、 patientなどのペア情報をdesignに入れることで、 個体差を考慮することができます。 ただし、edgeRやDESeq2は、 単にペア内の変化方向を見るだけではありません。 Gene Countsの平均値と分散の関係、低Countの不確実性、 モデル全体の分散推定も同時に扱います。
そのため、edgeRやDESeq2でpaired designを考慮したからといって、 必ずしも有意遺伝子数が増えるとは限りません。 場合によっては、対応のない2群比較よりも有意遺伝子数が減ることもあります。
paired t-testとedgeR・DESeq2のどちらが正しいかを論じるより、結果にこのような違いがあることを知っておくことが、実務上重要です。
ケース4：low-input RNA-Seqデータ
このケースの詳しい検証は、 low-input RNA-SeqでedgeR・DESeq2・t検定を比較した事例 で紹介しています。
4つ目のケースは、low-input RNA-Seqです。 この場合は、どの統計手法を使うか以前に、 データそのものがDEG解析に耐える再現性を持っているかを確認する必要があります。
100 cellsのような極端に低いインプット量では、 繰り返しサンプル間の再現性がほとんど見られず、 発現差として大きく見える要因の多くが、 偶然増幅されたか、あるいは偶然検出されなかったかに由来していました。
このようなデータでは、 edgeRやDESeq2の統計モデルが想定する通常のbulk RNA-Seqの分散構造とは大きく異なります。 そのため、edgeRやDESeq2で有意判定された遺伝子であっても、 その多くが生物学的な発現差ではなく、 増幅や検出の有無に由来する見かけ上の差を反映していても不思議ではありません。
これはedgeRやDESeq2が誤っているというより、 入力データの構造が、 通常のbulk RNA-Seqを前提とした分散モデルから大きく外れているために生じる当然の結果です。
1000 cellsになると、 散布図上では少しだけRNA-Seqらしい分散構造が見え始めます。 つまり、高発現領域ではばらつきが小さく、 低発現領域に向かうほどばらつきが広がるという、 edgeRやDESeq2が想定する平均値と分散の関係に近い形が見えるようになります。
しかし、分散構造が少し見え始めることと、 edgeRやDESeq2によるDEG解析の結果が信用できるようになることは別です。 1000 cellsでも、通常のbulk RNA-Seqと比べればデータはまだ不安定であり、 増幅や検出の偶然性は残ります。
このデータでは、100k cellsの結果を参照データとして使うことができました。 1000 cellsの結果を100k cellsの結果と比較すると、 一貫した変化として確認できたのは、 Gene Countsが5000以上の非常に高発現な遺伝子に限られていました。
つまり、1000 cellsで見かけ上edgeRやDESeq2が想定するような分散構造に近づいて見えても、DEG解析の結果をそのまま信頼できるとは限りません。低インプットのRNA-Seqデータから得られたDEGリストは、非常に高発現の遺伝子に限定して見るなど、極めて慎重に扱う必要があります。
データタイプ別の実践ルール
| データタイプ | 主な問題 | edgeR / DESeq2での注意点 | t検定での注意点 | 実践ルール |
|---|---|---|---|---|
| 低分散・少数サンプルのin vitroデータ | 群内分散が極端に小さく見える。nが少ないため、遺伝子ごとの分散推定が小さくなりがち。 | 同程度のCountsを持つ遺伝子群から分散を推定するため、観測分散が大きい遺伝子でも有意になることがある。このような場合、t検定より有意判定が甘くなる。 | 観測された群内分散が極端に小さい遺伝子が有意になりやすい。 | p値だけでなくfold change条件を併用する。fold changeを加えると3手法の差は小さくなる。 edgeR / DESeq2の場合、観測されたばらつきの大きい遺伝子をDEGリストから除外することも検討する。 |
| biopsy由来の中規模データ | サンプル数はある程度あるため、分散推定は安定してくる。ただし、個体差、組織構成、総リード数、ダイナミックレンジなどの違いが、DEGの信頼性を不安定にする要因となりうる。 | ノイズとの境界領域で、ダイナミックレンジ差に由来する見かけ上の差を有意と判定することがある。 観測分散が大きい遺伝子でも、有意になることがある。 |
ノイズとの境界領域で、ダイナミックレンジ差に由来する見かけ上の差を有意と判定することがある。 観測された群内分散が極端に小さい遺伝子が有意になりやすい。 |
正規化後Gene Countsの分布を確認し、最もダイナミックレンジが狭いサンプル群に合わせてlow signal cutoffを設定する。シグナル領域では通常のDEG解析を行い、低Count領域ではp値に頼らず、測定値の有無に基づいてON/OFF型の遺伝子を別枠で抽出する。 |
| paired design | 個体差を考慮しないと、患者間のベースライン差が群間差の評価に混ざる。 | patientをdesignに入れて個体差を考慮できるが、Countモデル全体の不確実性も評価するため、必ずしも有意遺伝子数が増えるとは限らない。 | ペア内の差を直接評価するため、患者内で一貫した小さな変化を拾いやすい。ただし、これらが完全に生物学的に意味のないノイズとは言い切れない。 | ペア構造は必ず考慮する。paired t-testのみの遺伝子も、患者内で一貫した小変化としてクラスタリングや機能解析で確認する。 |
| low-input RNA-Seq | 増幅や検出の偶然性が大きく、通常のbulk RNA-Seqの分散構造から外れやすい。 | 有意判定された遺伝子が、増幅・検出の有無に由来する見かけ上の差を反映していることがある。 | 観測値の不安定さをそのまま反映しやすく、p値やfold changeが大きく見えても信頼できるとは限らない。 | p値だけでDEGリストを信用しない。見た目よりもさらに保守的に判断する。 |
詳しい検証記事
- 事例1：低分散・少数サンプルのin vitroデータ
- 事例2：biopsy由来の中規模データ
- 事例3：paired designで何が変わるのか
- 事例4：low-input RNA-Seqでは何が起きるのか
まとめ：データの性質や研究目的に合わせて使い分ける
t検定は、edgeRやDESeq2によって置き換えられるべき 古く未熟な方法というわけではありません。 より正確には、 t検定は各遺伝子で実際に観測された分散を重視する方法であり、 edgeRやDESeq2は遺伝子間で分散情報を共有することによって、 分散推定の安定性を重視する方法です。 どちらの方法が適しているかは、 データの性質と解析目的によって変わります。
サンプル数が少なく、 低Count遺伝子が多く、 個々の遺伝子の分散推定誤差が結果を大きく左右する場合には、 dispersion shrinkageは合理的です。 ただし、ケース1とケース2で見たように、 p値だけで比較した場合には手法間の違いが見られても、 fold change条件を組み合わせると、 最終的に抽出されるDEGリストは互いに近いものになります。
一方、十分なサンプル数があり、 データの品質が高く、 低Count遺伝子が適切に除外され、 効果量や遺伝子固有のばらつきも重要な情報として扱いたい場合には、 分散推定を安定させるために遺伝子間で情報を共有するメリットは小さくなり、 個々の遺伝子で観測された分散を全体傾向へ近づけるデメリットが相対的に大きくなります。 そのような条件では、各遺伝子で観測された分散を直接使う方法に十分な合理性があります。
したがって、本来問うべきなのは、
「RNA-Seqではt検定を使うべきではない？」ではなく、
どの程度、遺伝子間で分散情報を共有することが妥当な状況なのか
ということです。
dispersion shrinkageは、
「分散推定の不確実性とバイアスの交換」と考えることができます。
言い換えると、
分散推定の不安定さを減らす代わりに、
個々の遺伝子で観測された分散を、
同程度のGene Countsを持つ多数の遺伝子から推定された全体傾向へ近づける処理です。
この交換が有利に働く状況もあれば、
遺伝子固有のばらつきを重要な情報として扱いたい場合には、
不利に働く状況もあります。
実務的には、 どちらか一方を常に選ぶのではなく、 サンプル数、低Count遺伝子の割合、データ品質、 観測された分散や効果量、 そして解析目的を確認しながら使い分ける必要があります。 t検定もedgeRやDESeq2も、 解析者が目的に応じて選択するツールに過ぎません。
edgeRやDESeq2と、t検定では、甘い有意判定になりやすい発現パターンが異なります。 t検定は、観測分散が小さい遺伝子が有意判定されやすいのですが、 fold change条件を組み合わせることで対処可能です。 一方、edgeRやDESeq2では、観測分散が大きい遺伝子でも有意と判定されることがあります。 変動の大きい遺伝子が重要ではないと考えられる場面では、 CVが大きい遺伝子をDEGリストから除外するなどの対応を検討してもいいでしょう。 また、複数の手法で得られたDEGリストを比較し、 共通部分や手法ごとの特徴を可視化して確認することは、 単一のp値リストに頼るよりも、実務上有用な判断材料になります。
また、一段視点をあげて、検定だけでなく、その前後関係も考慮して判断することが 実務上重要になります。 さらに、どの手法を使う場合でも、 解析結果がどのような特徴を持つのかを可視化しながら理解することが、 DEG解析の基本になります。
関連トピック
ChatGPTを使ったedgeR・DESeq2解析（RによるRNA-Seqチュートリアル）
