EdgeRと一般化線形モデルまわり

ちょっと仕事も含めて、EdgeRの一般化線形モデルをじっくり見る必要があったので、そのあたりのメモ。

EdgeRは以前見たときは、1因子はExact Test をしていて、多因子の時が一般化線形モデルという話だと思ったけれど、最近は1因子でも一般化線形モデルという事らしい。DESeqもそうなっているので、まぁそういうことなのだろう。

一般化線形モデルはリンク関数を使って、応答変数を変換し、線形予測子でモデルを作るようにしたもので、応答変数が正規分布をしてなくてもいいというところで応用の幅が広い（ただし指数分布族である必要はある）。

で、いきなりEdgeRでのモデルを見てみる。マニュアル18ページに出てる。ちなみに2015年10月8日にReviseされたバージョン。

$\begin{align} \log\mu_{gi} = {\bf x}^T_i \beta_g + \log N_i　\end{align}$

ここで $\mu_{gi}$ は $g$ 番目の遺伝子、 $i$ 番目のサンプルの観測値の期待値、 ${\bf x}^T_i$ はデザイン行列、 $\beta_g$ は係数。そして $\log N_i$ はサンプル $i$ 番目のライブラリサイズ、つまり総リード数。ここで、ふと思ったわけです、この $\log N_i$ ってなんでそこに居るの、と。普通教科書に出てくる一般化線形モデルは、以下

$\displaystyle \log\mu_{gi} = {\bf x}^T_i \beta_g$

みたいな感じ。で、とりあえず分かりやすくするために左辺をもとの値に戻してみるかと、式変形してみると

$\begin{align} \log\mu_{gi}&={\bf x}^T_i \beta_g + \log N_i \\\ \exp(\log\mu_{gi})&=\exp({\bf x}^T_i \beta_g + \log N_i) \\\ \mu_{gi}&=\exp({\bf x}^T_i \beta_g)N_i \\\ \end{align}$

で、マニュアルの14ページあたりに、以下の記載があるわけです。

$\begin{align} E(y_{gi})=\mu_{gi}=N_i \pi_{gi} \end{align}$

ここで、 $y_{gi}$ が実際の観測されたリード数、その期待値が、総リード数に割合 $\pi_{gi}$ をかけたものであらわされていて（ $\Sigma^G_{g=1}\pi_{gi}=1$ ）なんだ、そういうことだったんですか、Robinsonさんと著者の人を思っていました。まぁ自明だから説明省くよ、って事なのかもしれませんが、なるほど、総リード数の割合のところを推定するようなモデルなのだなーとなるほどなるほどと思っていました。
バイオインフォのアルゴリズムは、当たり前といえば当たり前なのでしょうが、教科書どおりとはならないので、有る意味解読的に読まないといけませんが、おや？と思ったところが分かると、パズルが解けたような面白さですね。と、教科書には載ってないと書いて、よくよくいろいろな本を見てみると、一般化線形モデルの教科書的な久保先生のみどり本には、47ページ脚注に説明があり、138ページで同様な形のGLMが紹介されていました。あはは。。。GLM奥深し。。大学院で習った記憶はあるものの、あまり専門ではないのですが、これを機に勉強しなおそうかなとか思ったり。

まぁSlow and Steady でいきませう。