一見奇妙な事実
RNA-seq解析を行う際、ふつう我々は「発現量のカウント（read counts）は負の二項分布（Negative Binomial distribution, NB）に従う」と言う。この結論は、DESeq2、edgeR、limma-voom といった主要な差次的発現解析ツールのほぼすべてのドキュメントに登場し、さながら自明の前提のように扱われている。
しかし、なぜか。
私は多くの資料を調べたが、大半の説明は浅すぎる（「データが過分散だから」）か、いきなり数式に飛ぶ（「NBの確率質量関数は……」）か、あるいは「待ち時間」の話にすり替わる（「NBは $r$ 回目の成功までに必要な試行回数を記述する」）——最後の説明はとりわけ混乱を招く。RNA-seqのカウントと「待ち時間」は何の関係もないからだ。
そこで、この話を徹底的に整理したい。最も単純な例から出発し、なぜNBが自然な選択なのかを一歩ずつ導き、かの古典的な待ち時間による定義が実は我々の導出とは無関係であることも示す。
直感的な例から始めよう
ある遺伝子の細胞Aと細胞Bにおける発現量を調べたいとする。
細胞の内部では何が起きているのか。 遺伝子が転写されるには、まずRNAポリメラーゼがプロモーターに結合し、そこから合成が始まる。ポリメラーゼはランダムにプロモーター領域へ到達し離脱するため、転写そのものが確率的な過程である——任意の時間窓の中で合成されるmRNA分子の数はランダムである。
観察時間を固定すれば（たとえば1秒間）、転写によって産生されるmRNAの数はおおむねPoisson分布に従う。理由はこうだ：mRNAの産生は一連の独立な「成功」事象（1回の転写開始）とみなせ、固定された時間窓において、Poisson分布はこの種のカウント過程を記述する自然なモデルだからである。
しかしここで問題が生じる——Poisson分布には強い仮定がある：平均と分散が等しい。
単一細胞であれば、平均と分散の関係はたしかにその通りになりやすい。しかし一群の細胞（たとえば組織サンプル中の数万個の細胞）を比較する場合、状況はまったく異なる。同じ遺伝子であっても、細胞ごとの発現量は大きくばらつく——まったく発現していない細胞もあれば、高発現の細胞もある。この細胞間の異質性（biological variability）は、集団全体の分散を平均よりもはるかに大きくする。
これがいわゆる「過分散」（overdispersion）である。Poisson分布ではこの問題に対処できない。
Gamma-Poisson階層モデル
解決策は階層的な発想から生まれる。
まず認めよう：細胞ごとの真の発現率（$\lambda$）はそもそも異なる。この違いをひとつの分布で記述できる。どの分布を選ぶか。統計学において、Gamma分布は自然な選択である——数学的に扱いやすく、かつ多様な形状の確率変数を記述できるだけの柔軟性を備えているからだ。
こうして我々は二層モデルを構築する：
第一層（細胞間変動）： $\lambda \sim \text{Gamma}(\alpha, \beta)$
ここで $\lambda$ は各細胞の真の発現率（単位時間あたりの平均mRNA産生量）を表す。$\alpha$ と $\beta$ はGamma分布のパラメータである。平均 $\mu = \alpha/\beta$、分散 $= \alpha/\beta^2$。
第二層（細胞内のランダム性）： $\lambda$ を与えられたもとで、カウント $X$ はPoisson分布に従う：$X | \lambda \sim \text{Poisson}(\lambda)$
この層が記述するのは：たとえ二つの細胞がまったく同じ $\lambda$ を持っていたとしても、転写過程そのもののランダム性により、実際に観測されるカウントにはゆらぎが生じる、ということである。
ここで $\lambda$ を積分消去（marginalize out）し、$X$ の周辺分布を求めよう：
$$P(X = k) = \int_0^\infty P(X=k|\lambda) \cdot P(\lambda) \, d\lambda$$この積分の結果は：
$$P(X = k) = \frac{\Gamma(k+\alpha)}{k! \cdot \Gamma(\alpha)} \left(\frac{\beta}{1+\beta}\right)^\alpha \left(\frac{1}{1+\beta}\right)^k$$そしてこれこそが——まさに負の二項分布のパラメトライズ形式である。
通常NBはNB(r, p)と書かれ、ここでの対応関係は：$r = \alpha$、$p = \beta/(1+\beta)$ である。あるいは平均-分散パラメトライゼーションを用いれば：平均 $\mu = r(1-p)/p$、分散 $= \mu + \mu^2/r$。
重要な洞察： 導出された分散は $\mu + \mu^2/r$ であり、$\mu$ ではない。これは分散が平均より大きいことを意味する——モデルは過分散を自然に内包している。そして $r$ が小さいほど、分散の程度は大きくなる。
これこそが、RNA-seq解析において負の二項分布がこれほど自然である理由である：NBは二層の変動——細胞間の異質性（Gamma層）と細胞内のランダム性（Poisson層）——を同時に捕捉する。
Poissonではなぜ駄目なのか
もし無理やりPoissonモデルを使おうとすれば、すべての細胞の $\lambda$ が同一だと仮定する以外に道はない。しかしそれは非現実的である。実データでは、実験条件、組織部位、細胞周期の段階によって、遺伝子発現率は異なる。
Poissonで同一細胞集団のデータに適合させると、分散が深刻に過小評価され、偽陽性率が急上昇する——実際には差がないのに、差があると誤判定されてしまうのである。
負の二項分布は、追加パラメータ $r$（しばしば「dispersion parameter」と呼ばれる）を導入することで、この異質性に適応するだけの柔軟性をモデルに与えている。
「待ち時間」とは何の関係もない
古典的な負の二項分布の教科書はこう教える：NBは「$r$ 回目の成功が起きるまでに必要な試行回数」を記述する、と。たとえばコイン投げで、3回目の表が出るまでに合計何回投げたか。
この定義は、工業品質管理や待ち行列理論などの分野ではたしかに有用である。しかし——それはいま我々が行った Gamma-Poisson 導出とは何の関係もない。
この二つの道筋は、たまたま同じ数学的形式に行き着いただけである。
確率論では、異なる確率過程が同じ分布に収束することが多々ある——これは「分布の同等性」と呼ばれる。NBもその一例だ：「待ち時間」の観点からもNBが得られるし、Gamma-Poisson階層モデルからもNBが得られる。しかし因果の連鎖はまったく異なる。
RNA-seqの文脈で、待ち時間を使ってNBを理解しようとすれば混乱を招くだけである（「我々は何を待っているのか？」）。一方、Gamma-Poissonの枠組みはずっと直感的だ：細胞ごとに発現率が異なり（Gamma）、各細胞内部のカウントはランダムであり（Poisson）、両者が重なり合ってNBが生まれる。
だから次に誰かが「NBの dispersion parameter $r$ は自由度を表す」だの「待ち時間」だのと言うのを聞いたら、警戒してほしい。それは別の視点であり、我々の生物学的応用とは直接のつながりはない。
最後に
統計モデルの意味は、それがどのような物語を語るかによって決まる。
負の二項分布には、数学的に等価な定義・導出が複数存在する。しかしそれらは等価な直感ではない。RNA-seqにとって、Gamma-Poissonの物語が最も実態に即している——我々はたしかに一群の異質な細胞を測定しており、細胞間にはたしかに変動があり、Poissonではたしかにその変動を記述できないからだ。
もうひとつ興味深い点：dispersion parameter $r$ を推定するということは、実質的にこの集団がどれほど異質かを推論していることになる。$r$ が大きければ、分散は平均に近づき、NBはPoissonに縮退する——これは細胞間が比較的均一であることを示す。$r$ が小さければ、分散は平均よりはるかに大きい——異質性が強いことを示す。これはとりわけシングルセルデータにおいて意味を持つ。
だから次にDESeq2やedgeRを走らせるときは、その背後でどのような物語が語られているかを忘れないでほしい。
付録：完全な導出
第一歩：何を計算すべきかを明確にする。
我々には二層モデルがある：
- 第一層：$\lambda$ は Gamma 分布に従う。すなわち $\lambda \sim \text{Gamma}(\alpha, \beta)$
- 第二層：$\lambda$ を与えられたもとで、カウント $X$ は Poisson 分布に従う。すなわち $X | \lambda \sim \text{Poisson}(\lambda)$
求めたいのは $X$ の周辺分布、つまり $\lambda$ を「積分消去」したものである：
$$P(X = k) = \int_0^\infty P(X=k|\lambda) \cdot P(\lambda) \, d\lambda$$第二歩：二つの分布の具体的な形を書く。
Poisson 分布（$\lambda$ を与えられたもとでの $X$ の条件付き分布）：
$$P(X = k | \lambda) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, ...$$Gamma 分布（$\lambda$ の分布）：
$$P(\lambda) = \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta\lambda}, \quad \lambda > 0$$ここで $\Gamma(\alpha)$ は Gamma 関数であり、階乗の一般化である（$\alpha$ が正の整数のとき、$\Gamma(\alpha) = (\alpha-1)!$）。
第三歩：積分式に代入する。
上記の二つの式を周辺分布の公式に代入する：
$$P(X = k) = \int_0^\infty \frac{\lambda^k e^{-\lambda}}{k!} \cdot \frac{\beta^\alpha}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta\lambda} \, d\lambda$$第四歩：被積分関数を整理する。
定数項を積分の外に出し、$\lambda$ を含む項をまとめる：
$$P(X = k) = \frac{\beta^\alpha}{k! \cdot \Gamma(\alpha)} \int_0^\infty \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda$$注意：
- $e^{-\lambda} \cdot e^{-\beta\lambda} = e^{-\lambda(1+\beta)}$
- $\lambda^k \cdot \lambda^{\alpha-1} = \lambda^{k+\alpha-1}$
第五歩：重要な変数変換のテクニック。
積分の中身の形を観察する：$\lambda^{\text{あるべき乗}} \cdot e^{-\lambda \cdot \text{ある定数}}$。これは Gamma 分布の密度関数によく似ている！
Gamma 分布の密度を思い出そう：$\frac{\beta^{\alpha}}{\Gamma(\alpha)} \lambda^{\alpha-1} e^{-\beta\lambda}$。ここで $\alpha$ と $\beta$ はそれぞれ形状パラメータとレートパラメータである。
我々の積分には、前についている正規化定数が欠けている。しかし、それを作り出すことができる！
$\alpha' = k + \alpha$、$\beta' = 1 + \beta$ とおく。すると：
$$\int_0^\infty \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda = \frac{\Gamma(k+\alpha)}{(1+\beta)^{k+\alpha}} \cdot \underbrace{\int_0^\infty \frac{(1+\beta)^{k+\alpha}}{\Gamma(k+\alpha)} \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda}_{= 1}$$積分の中身はまさに $\text{Gamma}(k+\alpha, 1+\beta)$ 分布の密度関数であり、積分結果は 1 になる。
したがって：
$$\int_0^\infty \lambda^{k + \alpha - 1} e^{-\lambda(1+\beta)} \, d\lambda = \frac{\Gamma(k+\alpha)}{(1+\beta)^{k+\alpha}}$$第六歩：元の式に戻す。
$$P(X = k) = \frac{\beta^\alpha}{k! \cdot \Gamma(\alpha)} \cdot \frac{\Gamma(k+\alpha)}{(1+\beta)^{k+\alpha}}$$整理する：
$$P(X = k) = \frac{\Gamma(k+\alpha)}{k! \cdot \Gamma(\alpha)} \cdot \frac{\beta^\alpha}{(1+\beta)^{k+\alpha}}$$$$= \frac{\Gamma(k+\alpha)}{k! \cdot \Gamma(\alpha)} \cdot \left(\frac{\beta}{1+\beta}\right)^\alpha \cdot \left(\frac{1}{1+\beta}\right)^k$$第七歩：負の二項分布と比較する。
標準的な負の二項分布 NB(r, p) の確率質量関数は：
$$P(Y = k) = \frac{\Gamma(k+r)}{k! \cdot \Gamma(r)} (1-p)^r p^k, \quad k = 0, 1, 2, ...$$比較すると：
- $r = \alpha$（形状パラメータ）
- $p = 1/(1+\beta)$（成功確率）
- $1-p = \beta/(1+\beta)$
完全に一致する！
第八歩：平均と分散の確認（任意だが推奨）。
平均：
全期待値の公式を用いて $E[X] = E[E[X|\lambda]] = E[\lambda] = \alpha/\beta$
分散：
全分散の公式を用いて
$$Var(X) = E[Var(X|\lambda)] + Var(E[X|\lambda]) = E[\lambda] + Var(\lambda) = \alpha/\beta + \alpha/\beta^2 = \alpha(1+\beta)/\beta^2$$負の二項分布のパラメータで表現すると：
- 平均 $\mu = r(1-p)/p$
- 分散 $= \mu + \mu^2/r$
我々の導出から：
- $\mu = \alpha/\beta$
- $Var = \alpha/\beta + \alpha/\beta^2 = \mu + \mu^2/\alpha$
したがって $r = \alpha$、つまり Gamma 分布の形状パラメータがそのまま NB の dispersion パラメータである。
導出全体の核心的発想は：
- 階層モデルの二つの分布を書き下す
- 隠れ変数 $\lambda$ について積分する（周辺化）
- Gamma 関数の積分性質を利用して計算を完了する
- 結果を整理すると、それが負の二項分布の形になっていることを確認する
この導出は、負の二項分布がなぜ RNA-seq データを自然に記述できるのかを示している：NB は本質的に Poisson 分布の「加重平均」であり、その重みは Gamma 分布によって与えられる。細胞間の異質性が大きいほど（Gamma 分布の分散が大きいほど）、最終的な分布の overdispersion はより深刻になる。