edgeRで使われる負の二項分布（Negative binomial distribution）とは何か

 

はじめに

RNA-Seq（RNAシーケンシング）実験では、まず複数のRNAサンプルからRNAを抽出し、シーケンシングを行います。次に、得られた大量の短い配列データ（リード）をゲノムにマッピングし、各遺伝子に由来するリードの数をカウントします。このカウントデータが、遺伝子の発現量の指標となります。サンプル $i$ の遺伝子 $g$ から得られたリード数を $y_{gi}$ とします。

なぜ負の二項モデルが必要なのか？

RNA-Seqのカウントデータは、ばらつきを持っています。このばらつきには、主に2つの要因が考えられます。

1. 技術的なばらつき (Technical variation): 同じRNAサンプルを繰り返しシーケンスした場合に生じる測定誤差です。これは、リードカウントがある程度ポアソン分布に従うとされています。ポアソン分布では、分散と平均が等しくなります。
2. 生物学的なばらつき (Biological variation): 生物学的な反復サンプル間（例えば、異なる個体や異なる培養細胞など）での遺伝子発現の真の変動です。同じ処理条件下であっても、個体差などにより遺伝子の発現レベルはサンプルごとに異なります。

もし技術的なばらつきだけしか存在しないのであれば、ポアソン分布で十分にデータを説明できます。しかし、実際には生物学的なばらつきが無視できないほど大きく、ポアソン分布で仮定されるよりも大きなばらつき（過分散）を示すことが一般的です。

そこで、この過分散をうまく捉えるために負の二項分布 (Negative Binomial distribution) が用いられます。

そもそも違う条件のものを比べるんだから差が出るのは当たり前ではないのか？

例えば、

• 新しい薬を病気のグループに投与する　vs　投与しないグループ
• 特定の栄養素を与えた植物　vs　与えなかった植物

こういった実験では、「薬を投与すれば症状が良くなるはず」「栄養素を与えれば成長が良くなるはず」といったように、条件の違いによって何らかの「差」が出ることを期待しています。

では、なぜわざわざ「負の二項モデル」のような少し複雑な話が出てくるのでしょうか？

それは、「見えた差」が本当に「条件の違いによる本物の差」なのか、それとも「たまたま起こった偶然の差」なのかを区別したいからです。

実験で何か「差」が見つかったとしても、その差が、

• 本当に「条件の違い（例：薬の効果）」によって生じた意味のある差なのか？
• それとも、「測定の誤差が大きめに出たサンプル」が集まった結果として見えているだけの、偶然の差なのか？**

この2つをしっかり区別する必要があるのです。

「負の二項モデル」のような統計的な方法は、まずこの「偶然起こりうるばらつきの大きさ」をきちんと評価します。その上で、観測された「差」が、その偶然のばらつきの範囲を大きく超えるほど確かなものなのかどうかを判断するのに役立ちます。

負の二項モデルの基本

遺伝子のカウント数が、実験のたびにどれくらいブレるか、その「ばらつき具合」を予測するために、以下のな数式を使います。
$$var(y_{gi})={\mu }_{gi}+{\varphi }_g{\mu }_{gi}^2$$

• $y_{gi}$ (ワイ・ジーアイ)：ある遺伝子 $g$ の、あるサンプル $i$ におけるカウント数
• ${\mu }_{gi}$ (ミュー・ジーアイ)：遺伝子 $g$ の平均値
• $var(y_{gi})$ (ヴァリアンス・ワイ・ジーアイ)：これが、遺伝子 $g$ の分散値。この数字が大きいほど、カウント数のブレが大きい。
• ${\varphi }_g$ (ファイ・ジー): 生物学的なばらつきを表す項です。${\varphi }_g$ はディスパージョン（分散係数） と呼ばれるパラメータで、遺伝子 $g$ に固有の生物学的なばらつきの大きさを反映します。

この式 $var(y_{gi})={\mu }_{gi}+{\varphi }_g{\mu }_{gi}^2$ は、大きく分けて2つの足し算でできています。これは、分散値が2種類のブレから成り立っていることを示しています。最初の部分「${\mu }_{gi}$」：これは「機械のせいで起こる小さなブレ」に相当します。次の部分「${\varphi }_g{\mu }_{gi}^2$」：これが「生き物の個性による大きなブレ」に相当します。

この式の両辺を ${\mu }_{gi}^2$ で割ると、変動係数の二乗 (CV${}^2$) の関係が見えてきます。

$$C{V}^2(y_{gi})=\frac{1}{{\mu }_{gi}}+{\varphi }_g$$

ここで、

• $1/{\mu }_{gi}$ は、技術的なばらつきに起因するCV${}^2$です。
• ${\varphi }_g$ は、生物学的なばらつきに起因するCV${}^2$であり、生物学的変動係数 (Biological Coefficient of Variation, BCV) の二乗 (${\text{BCV}}_g^2$) にあたります。つまり、${\text{BCV}}_g=\sqrt{{\varphi }_g}$ です。

BCVは、生物学的反復サンプル間での遺伝子の真の存在量（発現レベル）が、平均に対してどれくらいばらついているかを示す指標です。シーケンスの深さ（リード総数）をどれだけ増やしても、この生物学的なばらつき自体は減りません。特に発現量の高い遺伝子では、技術的なばらつきの相対的な影響が小さくなるため、BCVが不確実性の主要な要因となります。したがって、遺伝子の発現変動を正確に評価するためには、このBCVをいかに信頼性高く推定するかが非常に重要になります。

edgeRにおけるBCVの推定

edgeRでは、このBCV、すなわちディスパージョンパラメータ ${\varphi }_g$ を推定するために、いくつかの段階的なアプローチを取ることができます。

1. 共通ディスパージョン (Common Dispersion): 全ての遺伝子でBCVが等しい（つまり、${\varphi }_g$ が全ての遺伝子で共通の値 $\varphi$ を持つ）と仮定する最も単純なモデルです。
2. トレンドディスパージョン (Trended Dispersion): 遺伝子の平均発現量（期待カウント数）に応じてBCVが系統的に変動する傾向（例えば、低発現の遺伝子ほどBCVが大きいなど）を考慮し、そのトレンドに沿ったディスパージョンを推定します。
3. タグワイズディスパージョン (Tagwise Dispersion; 遺伝子ごとのディスパージョン): 各遺伝子に固有のディスパージョン ${\varphi }_g$ を推定します。これが最も柔軟なアプローチであり、edgeRでは強く推奨されています。なぜなら、遺伝子によっては反復間のばらつきが非常に小さいものもあれば、非常に大きいものもあるため、遺伝子ごとに異なるBCVを許容することで、より現実に即したモデル化が可能になるからです。これにより、発現変動遺伝子のランキングにおいて、反復間で一貫して変動している遺伝子をより高く評価できます。

ただし、特にサンプル数が少ない場合、各遺伝子のデータだけで ${\varphi }_g$ を安定して推定するのは困難です。そこで、edgeRでは経験的ベイズ (Empirical Bayes) の戦略を用います。これは、個々の遺伝子の情報だけでなく、全遺伝子から得られる情報を「借りてくる」ことで、遺伝子ごとのディスパージョン推定値をより安定させる手法です。具体的には、個々の遺伝子から計算されたタグワイズディスパージョンの値を、共通ディスパージョンやトレンドディスパージョンで示される全体の傾向に向かって「縮約 (squeeze)」させます。この縮約の度合いは、データの信頼性などに基づいて統計的に決定されます。

準負の二項モデル (Quasi Negative Binomial Model)

さらにedgeRでは、負の二項モデルを拡張した準負の二項 (Quasi-Negative Binomial, QL) モデルも利用できます。これは、遺伝子特有のさらなるばらつきを考慮に入れるためのものです。このモデルでは、カウント $y_{gi}$ の分散は以下のように表されます。

$$var(y_{gi})={\sigma }_g^2({\mu }_{gi}+\varphi {\mu }_{gi}^2)$$

ここで、

• $\varphi$: NBディスパージョン。これは、全ての遺伝子に共通する、あるいは発現量のトレンドに従う、生物システム全体としてのばらつきを表します（タグワイズディスパージョンではなく、共通またはトレンドディスパージョンが用いられます）。
• ${\sigma }_g^2$: QLディスパージョン。これは遺伝子 $g$ に特有のばらつきで、全体の傾向（$\varphi$ で捉えられるばらつき）からさらに上振れしたり下振れしたりする部分をモデル化します。

観測されるカウントのばらつきの増加は、$\varphi$ または ${\sigma }_g^2$ の推定値の増加としてモデル化されます。

遺伝子特有のQLディスパージョン ${\sigma }_g^2$ の推定も、特に反復サンプル数が少ない場合には困難が伴います。そのため、ここでも経験的ベイズのアプローチが用いられます。まず、生のQLディスパージョン推定値に対して、平均発現量に依存するトレンドをフィットさせます。そして、生の推定値をこのトレンドに向かって縮約することで、より安定したEB (Empirical Bayes) 推定値を得て、これらを仮説検定に用います。これにより、推定の不確実性を減らし、検定力を向上させることができます。

まとめ

edgeRで使われる負の二項モデルは、RNA-Seqデータにおける技術的なばらつきと生物学的なばらつきの両方を考慮し、特に生物学的変動係数 (BCV) を適切にモデル化・推定することに重点を置いています。遺伝子ごとのディスパージョン（タグワイズディスパージョン）を経験的ベイズ法によって安定的に推定し、さらに準負の二項モデルを導入することで、より現実に近い複雑なばらつき構造を捉え、信頼性の高い発現変動解析を実現しています。