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解析では本質的に劣った手法である」 という単純な理解に変わってしまうと、話は不正確になります。
前処理後のデータでは、統計手法の優劣よりも出力傾向が重要になる
ここまで見てきたように、生のcount値やRPKM/FPKM、TPMをそのまま統計検定に使うことには注意が必要です。
しかし、低発現遺伝子を除外し、分布やサンプル間の関係を確認したうえで、
正規化や対数変換を行ったデータでは、議論の前提が変わります。
この段階のデータは、すでに生のcount値そのものではなく、
統計解析に適用しやすい形に整えられています。
そのようなデータに対しては、 「t検定か、負の二項分布か」という名前だけで、 どちらが本質的に正しいかを決めることはできません。 どの統計モデルも、本来の生物学的状態をそのまま正確に表しているわけではなく、 特定の前提のもとで差を評価するための近似だからです。
こうなると、実務上の焦点は、 理論上の優劣ではなく、 その手法を使うと、どのような遺伝子が拾われやすく、 どのような遺伝子が拾われにくいのかを理解することです。 そのうえで、 生物学的背景、臨床モデル、実験系の性質、そして実験生物学者としての経験や勘をもとに、 目的の遺伝子がどのように見えるはずかを考えることが重要になります。
たとえば、疾患モデルで一貫して小さく変化する遺伝子を重視したいのか、 処理によって大きくON/OFFする遺伝子を探したいのか、 一部のサンプルで強く反応するサブタイプ特異的な遺伝子を見つけたいのかによって、 適した手法や閾値の考え方は変わります。
発現量が高く、サンプル間で一貫して小さな差を示す遺伝子、 Fold Changeは大きいが低発現で測定値が不安定な遺伝子、 一部のサンプルだけで大きく外れる遺伝子などは、 統計手法や前処理の違いによって扱われ方が変わることがあります。
したがって、統計手法は「正解を出す装置」ではなく、 研究目的に照らして候補遺伝子を選び出すための道具です。 実際の解析では、その手法がどのような遺伝子リストを生みやすいのか、 そしてそのリストが研究目的や生物学的仮説に合っているのかを確認することが欠かせません。
統計モデルの誤り率は、生物学的な差そのものではない
p値やFDRは、統計モデルの前提のもとで計算された指標です。 それらは、モデル内で定義された誤り率や有意性を表すものであり、 生物学的な差そのものを直接示すものではありません。
たとえば、統計的には有意であっても、 Fold Changeが非常に小さい差が、 研究目的に照らして意味のある変化とは限りません。 逆に、p値がそれほど小さくなくても、 サンプル数、ばらつき、測定限界、外れ値、バッチ効果を考慮すると、 慎重に見るべき遺伝子が含まれることもあります。
このように、統計手法の理論的な前提と、実際のデータ解析での使われ方は、 必ずしも単純には一致しません。 t検定が理論的には問題を含みながらも、実用上はうまく働くことがある理由 については、別の記事でも詳しく解説しています。
また、bulk RNA-Seqのデータには、 protein-coding genesだけでなく、non-coding genesも含まれます。 通常のbulk RNA-Seqでは、リードカウントの多くはprotein-coding genesに割り当てられる一方で、 低カウント領域には、測定限界に近い値、低発現遺伝子、non-coding genes、 アノテーションやライブラリ調製の影響を受けやすい測定値が混在します。
つまり、低カウント領域は単に「発現量が低い遺伝子」の集まりではありません。 そこには、測定上の不確かさ、遺伝子種別の違い、ライブラリ調製の特性、 アノテーションの不確実性など、性質の異なる要因が重なっています。
このような複雑なデータを一つの統計モデルに入れたからといって、 生物学的な解釈まで自動的に保証されるわけではありません。 統計モデルは、データを解釈するための道具であって、 生物学的判断そのものではありません。
そのため、このような領域で重要になるのは、 統計モデルに判断を任せることではなく、 解析者自身が、研究目的に照らして「注目すべき遺伝子はどのような発現パターンとして見えるはずか」を考えることです。 たとえば、Fold差は小さいけれど分散も小さく、群間で一貫した差を示す遺伝子を重視するのか、 分散は大きいものの、特定のサブタイプや反応群で二値的に変化する遺伝子を候補として残すのかによって、 適した手法や閾値の考え方は変わります。
現代統計モデルの便利さが、バッチ効果を見えにくくすることもある
DESeq2やedgeRのような現代的な統計モデルは非常に有用です。 しかし、gene countを入力すれば、 標準的なパイプラインが自動的に発現差解析の結果を返してくれるという便利さは、 別の問題も生みます。
現代的な統計モデルそのものが問題なのではありません。 問題は、gene countを入力すればパイプラインが自動的に結果を返してくれるという便利さによって、 データの状態を十分に確認しないまま、標準手法として使われてしまうことです。 その結果、バッチ効果やサンプル構成の偏りを含んだままの解析結果が、 統計的に整ったDEGリストとして扱われてしまう危険が、 むしろ高まったという側面があります。 バッチ効果と判断する前に、PCAや発現分布を確認する重要性 については、別のケーススタディでも詳しく解説しています。
負の二項分布モデルは、 ライブラリサイズやカウントデータのばらつきを扱うためのモデルです。 しかし、実験日、試薬ロット、測定施設、処理順序、 サンプル調製条件などに由来するバッチ効果を、 自動的に見つけて取り除いてくれるわけではありません。
バッチ効果は、単なるランダムノイズではなく、 サンプル群や実験条件と結びついた系統的な歪みであることが多くあります。 もし、コントロール群と処理群が別の日に調製されていたり、 別のロットで測定されていたりすれば、 統計モデルはその差を生物学的な差として扱ってしまう可能性があります。
バッチ効果をモデル内で扱うには、 それを実験デザインやモデル式に明示的に組み込む必要があります。 さらに、条件とバッチが完全に重なっている場合には、 統計モデルだけで生物学的な差とバッチ効果を分離することはできません。
つまり、現代的な統計モデルを使ったからといって、 実験デザインやサンプル構成の問題が解決されるわけではありません。 入力データに構造的な偏りがあれば、 高度な統計モデルは、その偏りを含んだまま、 整った形式の結果を返してしまいます。
手作業でデータを見ることは、古いやり方ではない
手作業でデータを確認しながら解析することは、 単に古い、非効率な方法ではありません。 むしろ、PCA、クラスタリング、ヒストグラム、散布図、発現分布、 サンプル間の関係を確認する過程で、 バッチ効果や外れ値、サンプル構成の偏りに気づくことがあります。
パイプラインにgene countを入力して、 DEGリストだけを受け取る解析では、 その途中でデータの状態を確認する機会が少なくなります。 その結果、統計的には整った形式の結果が得られても、 その結果がどのようなデータ構造に支えられているのかを見落としやすくなります。
一方で、データを段階的に確認しながら進める解析では、 正規化前後の分布、低発現遺伝子の影響、PCAでのサンプル分離、 クラスタリングでのまとまり、各遺伝子の発現パターンを確認できます。 この過程は、統計モデルでは自動的に代替できません。
たとえば、PCAで条件ではなく実験日ごとにサンプルが分かれている場合、 それは発現差解析の前に確認すべき重要な情報です。 ヒストグラムでサンプル間の分布が大きく異なっている場合も、 単純に統計検定へ進む前に、その原因を考える必要があります。
このような確認を行わずに、 いきなりDEGリストだけを見ると、 バッチ効果、外れ値、正規化の影響、低発現遺伝子の混入などに気づかないまま、 統計的に有意な遺伝子を生物学的に意味のある遺伝子として扱ってしまう危険があります。
手作業でデータを見ることは、 統計モデルを否定するためではありません。 むしろ、統計モデルの結果を正しく解釈するために必要な作業です。
重要なのは、どの手法を信じるかではなく、データを見て判断すること
この記事は、t検定を推奨するための記事ではありません。 また、DESeq2やedgeRのような現代的なRNA-Seq統計モデルを否定するための記事でもありません。
問題にしたいのは、 「t検定はダメ」 「現代的な統計モデルを使えばよい」 という一方的な説明だけでは、 RNA-Seqデータ解析の判断として不十分だという点です。
生のcount値や十分に前処理されていないデータに対しては、 負の二項分布モデルの方がRNA-Seqのデータ構造を適切に扱いやすい、 という説明は妥当です。 しかし、適切に正規化・前処理・フィルタリングされ、 分布やサンプル間の関係を確認した後のデータでは、 単純に「負の二項分布の方が正しい」「t検定は不適切」とは言えません。
現代統計モデルであっても、一般的なt検定であっても、 どの統計手法であっても、結果はデータの状態、前処理、 フィルタリング、モデルの前提、しきい値、サンプル構成、 バッチ効果、外れ値などに依存します。 したがって、危険なのは特定の統計手法そのものではなく、 データの状態や手法の前提を確認しないまま、 やみくもに適用することです。
統計モデルは、答えを自動的に与える装置ではありません。 それは、データを解釈するための道具です。 重要なのは、どの手法を使ったかではなく、その結果がデータの分布、サンプル間の関係、 発現パターン、生物学的な文脈と整合しているかを確認することです。
RNA-Seq解析では、統計手法を選ぶ前に、まずデータを見る必要があります。
そして、統計結果を得た後にも、再びデータに戻って確認する必要があります。
データを見て、前提を確認し、結果を問い直す。
それが、RNA-Seq発現差解析において、
統計手法の名前以上に重要な作業です。