空間クラスタリング手法による多重突然変異のある疾患遺伝子のファインマッピング

方法

I人の個体より構成されるディプロイドデータにおいて表現型 $y_i, i=1,...,I$ で、ハプロタイプ $h_{ij},j=1,2$ である。 $h_{ij} \in \{1,...,H\}$ で、 $H$ はunique haplotypeの総数。

連続的表現型について、次のモデルを使用する。

$y_i=\sum^2_{j=1} \gamma_{c_{h_{ij}}} + \epsilon_i$

バイナリの表現型について、次のモデルを使用する。

$Probit [Pr(y_i=1)]=\sum^2_{j=1} \gamma_{c_{h_{ij}}} ...(1-1)$

ここで $c_{h_{ij}} \in \{1,...,C\}$ で、Cはクラスターの総数。(1-1)式はadditive modelである。dominant modelとしたいときには、

$Probit[Pr(y_i=1)]=max\{\gamma_{c_{h_{i1}}},\gamma_{c_{h_{i2}}}\}$

とすればよいだろう。

クラスターcのリスク $\gamma_c \sim N(\alpha,\delta_{\gamma}^2)$ とモデル化する。プライアーをあてはめてよいが、今回著者らは全ての場合で $\alpha=0, \delta_{\gamma}^2=1$ としている。これはつまりプロビットの式が

$p_i=Pr(y_i=1) \sim Unif(0,1)$

となるが、ケースコントロールスタディではこれが適切なモデルであると思われる。

これの証明

正規分布の累積分布関数を $F$ とする。(1-1)式は

$F^{-1}(p_i)=\gamma_c$

両辺にFをかけて

$p_i=F(\gamma_c) ...(2-1)$

累積分布関数の定義より、 $\gamma_c \sim N(0,1)$ より

$P(\gamma_c \leq y)=F(y) ...(2-2)$

さらに $p_i \sim Unif(0,1)$ は、次のように書ける。

$P(p_i \leq x)=x ...(2-3)$

すると、(2-3)式の左辺に(2-1)、(2-2)を順々に代入して

$P(p_i \leq x)=P(F(\gamma_c) \leq x)=P(\gamma_c \leq F^{-1}(x))=F(F^{-1}(x))=x$

したがって(2-3)式が証明された。

ひとつの表現型毎にひとつのハプロタイプが対応する場合、(1-1)式を次のように書き換えることができる。

$Probit [Pr(y_i=1)]=\gamma_{c_{h_{ij}}} ...(1-2)$

ボロノイ分割構造 Voronoi Tessellation Structure

確率的にアサインした「中心」ハプロタイプをもとにハプロタイプリスクのクラスタリングを行う。
- ここで $\mathcal{T}=(t_1,...,t_C)$ はC個の中心をあらわす。
- 中心ハプロタイプは、クラスターが由来する先祖ハプロタイプと考えられる。
- 中心ハプロタイプは、unobserved valuesであってよい。ただしマーカーがきわめて多い場合は観察データとしたほうがよいだろう
決定論的に、それぞれのハプロタイプを、最も近しい中心にアサインする。その「距離」の定義は後述。同じ距離に二つの中心が存在する場合は、リストの最初のほうを選択する。
ハプロタイプ領域 $R_c$ を次のように定義する。

$R_c=\{h \in \mathcal{H} : \|| h-t_c \|| \lt \|| h-t_{c'} \|| \forall c' \ne c\}$

が、後述する「距離」である。
- $\mathcal{H}$ はデータセットを構成するユニークハプロタイプのセットである。
- $R_c$ は、他のどの中心よりもクラスターcの中心に近い観察ハプロタイプの集合である。

これがボロノイ分割と呼ばれる手法である。中心についてはMetropolis-Hastingsアルゴリズムで求めた。

類似性の測度と遺伝子マッピング

類似性の測度は、クラスターc内の想定上の機能的突然変異、 $x_c$ により定める。 $x_c$ は簡便のため観察されたマーカー部位に限る。

ハプロタイプと中心との類似性は上の二つのハプロタイプのthe shared length identical by state (IBS)とする。
- $m_r$ は $x_c$ の右側でハプロタイプ $h$ と中心 $t_c$ とが最初に相同でなくなる部位である。
- $m'_r$ は $m_r$ のすぐ左のマーカーである。
次のように定義する。

$R_{ht_c}(x_c)=\frac{m_r+m'_r}{2}$

左については同様に $L_{ht_c}(x_c)$ とすると、類似性は次のように表す。

$w_{ht_c}=R_{ht_c}(x_c)-L_{ht_c}(x_c)$

ベイズ的マルコフ連鎖モンテカルロ(MCMC)推定法

Gibbs samplingを用いる。 (1-2)式は $y_i^* \sim N(\gamma_{c_h_i},1)$ を用いて次のようになる。

$y_i^*=\gamma_{c_h_i}+\epsilon_i...(5)$

ここで $\epsilon \sim N(0,1)$ である。

Gibbs samplerのそれぞれのサイクルでを発生させる。
- $y_i=1$ なら $y_i^* \gt 0$ とし、一方 $y_i=0$ なら $y_i^* \lt 0$ とする。

これで標準的MCMCが動く。

クラスター数が任意の数に固定されないようランダム化するため、revesible jump methodsを用いた。
- 通常使われる方法ではなくDenison and Holmes(2001)による。すなわちクラスター数が変化するごとに、(5)式の $\gamma_c$ の値をintegrate outする（積分して0のようにするの意か？）。
$\mathbf{y}^*=(y_1^*,...,y_I^*)$ 、 $\mathbf{\gamma}=(\gamma_1,...,\gamma_C)$ と定義する。
もう一度言うけど、 $\mathcal{H}$ はデータセットのユニークハプロタイプの空間で、 $\mathcal{T}$ は現時点でのハプロタイプ中心のセットである。
次の式をたてる。

$f(\mathbf{y}^*|\mathcal{T},\mathcal{H})=\int^{\infty}_{-\infty} ... \int^{\infty}_{-\infty} f(\mathbf{y}^*|\mathbf{\gamma},\mathcal{T},\mathcal{H}) f(\mathbf{\gamma})d_{\gamma_1}...d_{\gamma_C}...(6)$

新しいクラスターを追加したいときは、通常新たに $\gamma_c$ をたてる必要があるが、Denison and Holmesの方法によれば、 $\gamma$ 値のintegrate outした値にもとづく尤度構成にもとづくことができる。
Gibbs samplerのそれぞれのサイクルの最後において、式(5)にもとづき、それぞれの $\gamma_c$ のfull条件付き確率分布を使用して、新たな $\gamma$ のセットを発生させる。
これらのパラメータがupdateされたら、上記のとおり $y_i^*$ もupdateする。モデル式に対応した全てのパラメータの結合確率分布は

$f(\mathbf{y}^*|\mathbf{y},\mathbf{\gamma},\mathbf{c}_h)f(\mathbf{\gamma}|C)f(\mathbf{c}_h|C,\mathcal{T},\mathcal{H},\mathbf{x})f(\mathcal{T})f(\mathbf{x})$ .

詳しい話

(Appendix Aより)

Reversible Jump MCMC

prior assumption: 全てのハプロタイプ中心のlikelyは同等で、クラスターCの数は $[C_{min},C_{max}]$ の範囲に制限される。
中心の数についてはuninformative prior structureを用い、中心の提案と削減は同等の確率とし、Metropolis-Hastingsのupdateにおいてcenter proposalと事前確率分布をうまくcancel outできた。
(6)式を用いて、新しい中心を提案し、中心値の新しいセットを $\mathcal{T}^{new}$ とすることで新しいクラスターの提案(Birth Move)が可能である。現状のセットは $\mathcal{T}^{old}$ となる。
提案を受け入れるかどうかを、次のMetropolis-Hastings比によって決定する。

$R=min \left( 1,\frac{f(\mathbf{y}^*|\mathcal{T}^{new},\mathcal{H})}{f(\mathbf{y}^*|\mathcal{T}^{old},\mathcal{H})} \right)...(A1)$ .

連続変数を用いず離散的な中心のセットを扱っているので、(A1)式にヤコビアンは必要ないことに注目してほしい。
クラスターの削除はランダムにひとつ選択し、削除するというものである。
さらに、あるひとつの現在存在するクラスターを選択し、中心ハプロタイプを変化させるという状態を導入する。

クラスター配置に影響を与えないupdate

式(5)のとおりに $y_i^*$ と $\gamma_c$ をupdateする。

クラスター配置に影響するupdate

$z_i=y_i^*$ とし、 $I \times C$ 行列 $\mathbf{X}$ を定義する。 $c_h_i=j$ なら $\mathbf{X}_{ij}=1$ であり、それ以外なら0。式(5)はつぎのようになる。

$z=X\mathbf{\gamma}+\epsilon$

prior $\gamma \sim \mathbf{N(\mu,\Sigma)}$ とすると $f(\gamma|z) \sim \mathbf{N(\mu^*,\Sigma^*)}$ で、このとき

$\mathbf{\mu^*}=(\mathbf{X'X+\Sigma^{-1})^{-1}(\Sigma^{-1}\mu+X'}z)$

$\mathbf{\Sigma^*=(X'X+\Sigma^{-1})^{-1}}$

パラメータの提案は次の周辺尤度に沿うMetropolis-Hasting比で評価される。

$f(z|\mathbf{\mu,\Sigma},\mathcal{T,H}) \propto \frac{|\Sigma^*|^{\frac{1}{2}}exp(-b^*)}{|\Sigma|^{\frac{1}{2}}(2\pi)^{\frac{I}{2}}}$

ここで

$b^*=\frac{(\mathbf{z'z+\mu'\Sigma^{-1}\mu-(\mu^*)'(\Sigma^*)^{-1}\mu^*})}{2}$

われわれのpriorは $\gamma \sim \mathbf{N(0,I)}$ であるよ。

次の5つのupdateのうちどれかを、均一乱数にしたがって選択し行う。いずれも式(6)に基づく提案の受け入れor拒絶をおこなうためのupdateである。

Update $x_c$: 非標準的方法により機能的変異 $x_c$ の場所をサンプリングする。確率0.75で、現在の値のランダムな摂動に基づき $x_c$ の新しい値を提案する、すなわち $x_c^{new}=x_c^{old}+\Delta$ で $\Delta$ はサンプラーのburn inフェイズで変化する。