RNA-Seqの発現差解析では、edgeRやDESeq2のようなGene Countsに対応した統計手法がよく使われます。 これらの手法は、RNA-Seq Countsの性質を考慮して設計されており、 DEG解析の基本となる方法です。
一方で、実際のRNA-Seqデータでは、 in vitro実験かbiopsy由来データか、繰り返し数が多いか少ないか、インプットRNA量が多いか少ないか、paired designかどうかなどのデータタイプによってデータの性質が異なり、その結果として有意判定される遺伝子の発現パターンの特徴が変わります。
本ページでは、edgeR、DESeq2、t検定の原理的な違いを整理したうえで、 代表的な4つのデータタイプにおいて、 それぞれの手法でどのような遺伝子が有意と判定されやすいのかを見ていくことで、これらの手法をどのように使い分ければよいかを考えます。
t検定とedgeR・DESeq2の基本的な違い
t検定とedgeR・DESeq2の大きな違いは、 分散をどのように扱うかにあります。
t検定は、各遺伝子について観測されたサンプル間のばらつきに基づいて、 群間差が大きいかどうかを評価します。 つまり、各遺伝子で観測された群内分散に強く依存します。
一方、edgeRやDESeq2は、 その遺伝子自身のCountsだけでなく、 同程度のCountsを持つ他の遺伝子のばらつきも利用して、 RNA-Seqデータ全体における平均値と分散の関係を推定します。
RNA-Seq Gene Countsには、 発現量が高い領域では相対的なばらつきが小さく、 発現量が低くなるほどばらつきが大きくなる性質があります。 さらに、低発現領域では、 あるサンプルでは測定され、別のサンプルでは測定されないという、 検出の偶然性も強くなります。
edgeRやDESeq2は、このようなCountデータの平均値と分散の関係を考慮して、 発現差を評価します。 そのため、edgeRやDESeq2で有意となった遺伝子を散布図上で見ると、 低Count側に向かって対角線から離れるような、 RNA-Seq Counts特有の分散構造を反映した分布として見えます。ただし、これは低Count領域のすべての差が信頼できるという意味ではありません。 そのため、edgeRやDESeq2を使う一般的なプロトコルでは、低発現遺伝子を解析から除外するのが普通です。
以上の特徴を踏まえて、ここからは代表的な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条件を加えた後でも、 手法間の違いは残ります。 まず、edgeRやDESeq2では有意となる一方でt検定では有意とならない遺伝子は、 観測された群内分散が比較的大きい傾向があります。これは、その遺伝子の実際に観測された分散ではなく、同程度のCountsを持つ他の遺伝子のばらつきも利用して分散が推定されることに由来します。in vitroの低分散データでは、観測されたばらつきよりも小さな分散を割り当てられることがあり、結果としてt検定よりもリベラルな有意判定になることがあります。
逆に、t検定のみで有意となる遺伝子は、 観測された群内分散が極端に小さい傾向があります。 このような遺伝子では、 たまたま繰り返し間の値がそろったことで、 t検定上は小さなp値になりやすくなります。このタイプの遺伝子については、fold changeの閾値を適切に設定して、小さすぎる変動を除外することが有効です。この考え方は、t検定だけでなく、edgeRやDESeq2を使う場合にも共通します。
ケース2：biopsy由来の中規模データ
このケースの詳しい検証は、 biopsy由来の中規模データでedgeR・DESeq2・t検定を比較した事例 で紹介しています。
次に、biopsy由来の中規模RNA-Seqデータを考えます。 このようなデータでは、ケース1と比べてサンプル数が多く、 発現量が十分に高いシグナル領域では、 分散の見積もりが比較的安定します。そのため、発現量が十分に高い領域では、 edgeR、DESeq2、t検定の結果の違いは小さく、 どの手法を用いても主要な発現変動遺伝子は共通して検出できると考えていいでしょう。
一方で、このケースでは発現量が低い領域において、 edgeRやDESeq2のみで有意と判定される遺伝子が多数見られました。 このような遺伝子を詳しく確認すると、 明確な生物学的発現差というより、 サンプル間のダイナミックレンジの違いに由来する 正規化後Gene Countsの下限値のズレと対応していました。
RNA-Seqデータでは、サンプル間で総リード数に数倍の差があることは珍しくありません。 ただし、総リード数の差とダイナミックレンジの差は単純に対応するわけではありません。RNA品質、サンプル組成、特定の高発現遺伝子へのreadの偏り、マッピング率など、複数の要因が重なるためです。いずれにせよダイナミックレンジの差は、正規化後のGene Countsの下限値のズレにつながります。そして、この下限値のズレが一方の群に偏ると、 群間差のように見える原因になります。
edgeRやDESeq2の統計モデルは、 低発現領域のCountsが不安定になりやすいことを考慮して設計されています。 しかし、サンプル間のダイナミックレンジの違いによって、 正規化後Gene Countsの下限値がずれることはモデルに組み込まれていません。
したがって、低Counts領域でedgeRやDESeq2のみ有意となる遺伝子については、 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の平均値と分散の関係、低Countsの不確実性、 モデル全体の分散推定も同時に扱います。
そのため、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手法の差は小さくなる。 |
| biopsy由来の中規模データ | サンプル数はある程度あるが、個体差、組織構成、総リード数、測定下限の違いが混ざる。 | 低Counts領域で、ダイナミックレンジ差に由来する見かけ上の差を有意と判定することがある。 | シグナル領域では解釈しやすいが、低Counts領域や欠損を含む領域には注意が必要。 | ダイナミックレンジが最も狭いサンプルに合わせて、正規化後Gene Countsにlow signal cutoffを適用する。 低Counts領域のDEG抽出方法については、別途検討する。 |
| 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では何が起きるのか
まとめ：データの性質や研究目的に合わせて使い分ける
RNA-SeqのDEG解析では、 Gene Countsを入力として扱えるedgeRやDESeq2を用いることが推奨されています。 しかし、edgeRやDESeq2を使えば、 常に解析結果をそのまま信頼できるわけではありません。 上記の通り、データセットの性質の違いにより、 得られるDEGリストも、その信頼性も大きく変わります。
一方、t検定はRNA-Seqの生のGene Countsにそのまま適用する方法ではありません。 しかし、適切に前処理・正規化・フィルタリングを行い、 シグナル領域に限定して使う場合には、 観測された発現パターンと対応づけて結果を解釈するための有用な比較手法になります。
それぞれの手法がどのような遺伝子を有意と判定しやすいのかを踏まえたうえで、 研究目的に合わせて手法を使い分けるとよいでしょう。 また、必要に応じて複数の手法で得られたDEGリストを組み合わせることも、実務上有効です。
