RNA-Seqの発現差解析では、DESeq2やedgeRのように、 負の二項分布を前提とした統計モデルが広く使われています。 一方で、t検定については、 「RNA-Seqには使うべきではない」 「カウントデータには不適切」 と説明されることが少なくありません。
この説明は、入門的な注意喚起としては有用です。 特に、生のcount値をそのままt検定にかけるような誤用を防ぐという意味では、 重要な役割を持っています。
ただし、 「負の二項分布モデルの方がt検定よりもRNA-Seqのデータ構造を的確に捉えている」 という説明には、前提があります。 それは、生のcount値、あるいは低発現遺伝子や分布の歪みが十分に処理されていないデータを対象にする場合です。
生のcount値では、値は0以上の整数であり、 低発現領域には0や1が多く含まれます。 また、平均値と分散の関係や過分散も無視できません。 このようなデータに対して、正規分布を前提とするt検定をそのまま適用することは危険です。 その意味で、負の二項分布モデルは、 RNA-Seqの生カウントデータの性質をより適切に扱うためのモデルだと言えます。
しかし、適切な正規化、前処理、低発現遺伝子のフィルタリングを行い、 データ分布やサンプル間の関係を確認した後のデータでは、 この説明をそのまま当てはめることはできません。 その時点では、解析対象となるデータは、 すでに生のcount値そのものではなく、 統計解析に適用しやすい形に整えられたデータになっています。
そのようなデータに対しては、 「負の二項分布モデルの方が本質的に正しい」 「t検定はRNA-Seqには不適切」 という単純な優劣の議論は成り立ちにくくなります。 どちらも、データを解釈するための一つの数学モデルであり、 それぞれ異なる前提と限界を持つ方法として扱うべきです。
重要なのは、どの統計手法を使うかを先に決めることではなく、 データの状態を確認し、 前処理やフィルタリングの影響を理解し、 その手法が何を仮定し、何を評価しているのかを踏まえて結果を解釈することです。
なぜt検定は避けられるようになったのか
RNA-Seq解析でt検定が避けられるようになった背景には、 単純な数学的優劣だけでなく、実務上の理由があります。
RNA-Seqが普及し始めた時期には、 低発現遺伝子、0や1の多いカウントデータ、 ライブラリサイズの違い、過分散などを十分に考慮しないまま、 生カウントや不適切に処理されたデータに対して 一般的な統計検定を適用する例が少なくありませんでした。
このようなデータに対して単純なt検定を適用すると、 低カウント領域の不安定な値や、 サンプル間のライブラリサイズ差、 平均値に依存する分散の違いなどによって、 結果が大きく歪む可能性があります。
さらに、RPKMやFPKMの普及も、この問題を分かりにくくした要因の一つです。 RPKMやFPKMは、ライブラリサイズや遺伝子長で補正された値であるため、 一見すると「正規化済みの発現量」として、そのまま統計解析に使えるように見えます。
しかし、RPKMやFPKMに変換しても、 低カウント値そのものの信頼性が高くなるわけではありません。 raw countが1と2の違いであっても、変換後には2倍の差として見えることがあります。 このような差は、測定上の不確かさやサンプリングノイズを強く含んでいる可能性があります。
つまり、RPKMやFPKMは、低カウント領域のばらつきを解決するのではなく、 その不安定さを「正規化済みの連続値」のように見せてしまうことがあります。 その結果、RPKMやFPKMに対してt検定などの一般的な統計検定を適用すると、 低発現遺伝子由来の不安定な差を、過剰に意味のある差として扱ってしまう危険があります。
その後、TPMが登場し、FPKMやRPKMに代わって使われる場面が増えました。 TPMには、FPKMよりもサンプル間で発現量を比較しやすいという利点があります。 しかし、TPMはここで議論している問題を本質的に解決したわけではありません。
TPMに変換しても、低カウント由来の不安定性、 発現量に依存する分散、少サンプル条件でのばらつきといった問題がなくなるわけではありません。 そのため、TPMをそのまま一般的な統計検定へ入力すればよい、という話にはなりません。 また、TPMを負の二項分布モデルの入力として用いることは、 モデルの前提から外れており、不適切です。
そのような誤用を避けるために、 DESeq2やedgeRのようなRNA-Seq専用の統計モデルが広く推奨されるようになりました。 これらのツールは、低発現遺伝子の扱い、ライブラリサイズの補正、 発現量に依存する分散、少サンプルでの分散推定などを、 一連の枠組みの中で扱うことができます。
この流れを後押ししたのが、Gene Countsをそのまま入力として扱えるという実務上の効率性です。 前処理を施していないGene Countsをそのまま単純なt検定にかけると、 大きな誤りを生みやすいのに対して、 RNA-Seq専用モデルはGene Countsを入力として受け取り、 一連の解析手順の中で処理できるため、 大量のデータを標準的な手順で処理するパイプラインに組み込みやすいという利点がありました。
この性質は、RNA-Seq解析の標準化に大きく貢献しました。 しかし同時に、 「Gene Countsを入力すれば標準モデルが処理してくれる」という安心感によって、 データの状態を十分に確認しないまま結果だけを受け取る解析を増やす要因にもなりました。
この意味で、 「RNA-Seqではt検定をそのまま使ってはいけない」 という説明は、初心者に対する安全側の注意喚起としては合理的です。 しかし、それがいつの間にか 「t検定はRNA-Seq解析では本質的に劣った手法である」 という単純な理解に変わってしまうと、話は不正確になります。
edgeR、DESeq2、t検定で有意判定されやすい遺伝子の特徴を理解し、使い分ける
ここまで見てきたように、生のcount値やRPKM/FPKM、TPMをそのまま統計検定に使うことには注意が必要です。
しかし、低発現遺伝子を除外し、分布やサンプル間の関係を確認したうえで、
正規化や対数変換を行ったデータでは、議論の前提が変わります。
この段階のデータは、すでに生のcount値そのものではなく、
統計解析に適用しやすい形に整えられています。
そのようなデータに対しては、 「t検定か、負の二項分布か」という名前だけで、 どちらが本質的に正しいかを決めることはできません。 どの統計モデルも、本来の生物学的状態をそのまま正確に表しているわけではなく、 特定の前提のもとで差を評価するための近似だからです。
こうなると、実務上の焦点は、 理論上の優劣ではなく、 その手法を使うと、どのような遺伝子が拾われやすく、 どのような遺伝子が拾われにくいのかを理解することです。 （関連ケーススタディ：edgeR・DESeq2・t検定の使い分け｜原理とデータタイプ別の実践ルール）そのうえで、 生物学的背景、臨床モデル、実験系の性質、そして実験生物学者としての経験や勘をもとに、 目的の遺伝子がどのように見えるはずかを考えることが重要になります。
たとえば、疾患モデルで一貫して小さく変化する遺伝子を重視したいのか、 処理によって大きくON/OFFする遺伝子を探したいのか、 一部のサンプルで強く反応するサブタイプ特異的な遺伝子を見つけたいのかによって、 適した手法や閾値の考え方は変わります。
発現量が高く、サンプル間で一貫して小さな差を示す遺伝子、 Fold Changeは大きいが低発現で測定値が不安定な遺伝子、 一部のサンプルだけで大きく外れる遺伝子などは、 統計手法や前処理の違いによって扱われ方が変わることがあります。
したがって、統計手法は「正解を出す装置」ではなく、 研究目的に照らして候補遺伝子を選び出すための道具です。 実際の解析では、その手法がどのような遺伝子リストを生みやすいのか、 そしてそのリストが研究目的や生物学的仮説に合っているのかを確認することが欠かせません。
一方で、生物学者が本当に探している遺伝子が、常に統計的に有意な発現差として現れるとは限りません。とくに、原因遺伝子や、原因遺伝子でなくとも上流のパスウェイに関連する遺伝子は、健常群や疾患群で高発現と低発現が入り混じる可能性が考えられます。（参照： t検定が理論的には問題を含みながらも、実用上はうまく働くことがある理由 ）このような遺伝子を抽出したい場合は、有意差の強い遺伝子だけでなく、P値が高めで通常のDEGリストからは外れやすい遺伝子群も探索対象としてクラスタリングで確認する、という使い方が有効なことがあります。
DESeq2やedgeRでも、偽陽性が問題になることはある
DESeq2やedgeRは、RNA-SeqのGene Countsを用いた発現変動解析で広く使われている標準的な手法です。 しかし、これらの手法を使えば、常に偽陽性が適切に抑えられるわけではありません。
実際に、ヒト集団由来の大規模RNA-Seqデータを対象にした研究では、 DESeq2やedgeRで、想定よりも多くの偽陽性が検出されることが報告されています。 この研究では、条件ラベルを入れ替えたpermutation解析を行い、 本来は群間差がないはずのデータでも、DESeq2やedgeRが多数の発現変動遺伝子を検出することを示しています。 （参照：Exaggerated false positives by popular differential expression methods when analyzing human population samples）
この問題の原因として、主に議論されているのは、 負の二項分布モデルへの適合不良、外れ値の影響、FDR制御のずれなどです。 その後の議論でも、外れ値をwinsorizationで抑えると偽陽性が減ることや、 シミュレーションで正規化の影響を無視すると人工的な偽陽性が生じることが指摘されています。
これらの議論では、 正規化後のヒストグラム左側に現れる低Counts領域の下限値のずれ、 つまりサンプルごとの実質的な検出限界の違いは、主要な原因としては扱われていません。 しかし、GEOに登録されている複数のデータセットを用いた私たちの解析では、 この問題が解析結果に影響しうることを確認しています。 （ケーススタディ：No.426（中規模データで一般的に観察されるレベルの事例）、No.403（解析結果を大きく歪めるレベルの事例））
私たちが言いたいのは、t検定を否定することでも、DESeq2やedgeRを否定することでもありません。 どの手法を使う場合でも、その手法が前提としているデータの見方と、 実際にどのような発現パターンを持つ遺伝子が抽出されやすいのかを理解したうえで、 データの状態と研究目的に合わせて適切に使い分けることです。
現代統計モデルの便利さが、バッチ効果を見えにくくすることもある
DESeq2やedgeRのような現代的な統計モデルは非常に有用です。 しかし、Gene Countsを入力すれば、 標準的なパイプラインが自動的に発現差解析の結果を返してくれるという便利さは、 別の問題も生みます。それは、データの状態を十分に確認しないまま、標準手法として使われてしまうことです。
その結果、バッチ効果やサンプル構成の偏りを含んだままの解析結果が、 統計的に整ったDEGリストとして扱われてしまう危険が、 むしろ高まったという側面があります。 バッチ効果と判断する前に、PCAや発現分布を確認する重要性 については、別のケーススタディでも詳しく解説しています。
手作業でデータを見ることは、古いやり方ではない
手作業でデータを確認しながら解析することは、 単に古い、非効率な方法ではありません。 むしろ、バッチ効果や外れ値、サンプル構成の偏りに気づきながら解析を進めるうえで、必要なプロセスです。
データを段階的に確認しながら進める解析では、 正規化前後の分布、低発現遺伝子の影響、PCAでのサンプル分離、 クラスタリングでのまとまり、各遺伝子の発現パターンを確認できます。
たとえば、PCAで条件ではなく実験日ごとにサンプルが分かれて見えた場合でも、単純に補正アルゴリズムを適用すればいいわけではありません。 確認すべきことは多岐にわたり、その過程を自動処理で代替することはできません。
統計モデルは万能ではないので、生物学的に重要そうな、詳しく確認すべき候補遺伝子を自動的に排除することもあります。
手作業でデータを見ることは、統計モデルを否定するためではありません。統計モデルによる自動化と、研究者自身による視覚的な確認を組み合わせることが、現実的で信頼性の高い解析戦略になります。
