Palantirの論文の解説．元論文はほとんど式の解説とかなくて，ほぼsupplementaryに数式による詳細が載っている．内容としては，コンピュータにヤバヤバな計算をさせる．サンプルの細胞の経路を複数調べたりできるけど，細胞の数を減らしたりする操作（scEpathならメタ細胞を作ったりした）が無いから，計算量はお察し．

データ幾何学（多様体）

データが２次元曲面に分布しているとする．これは３次元空間$\boldsymbol{R}^3$（多様体）に埋め込まれた部分多様体$L^2$（もう少し詳しく言えば超曲面）と考えられる．データを解析する場合は，３次元空間でやるより，これを２次元に直した方が圧倒的に効率が良い．よって，次元を落とすわけだが，この際にデータ間の距離を不変にしたい．この場合，距離は部分多様体での測地線距離とする．

diffusion mapは，データを低次元のユークリッド空間に埋め込む方法．Palantirではdiffusion mapを使う．

Nearest Neighbor Graph

nearest neighbor graphは似た細胞同士を結んだグラフ．当然のことながら，この長い経路は細胞の成長におけるトラジェクトリーに一致する．以下では，$k$-nearest neighbor graph，つまり$k$個の近い細胞を結んだグラフ$G_X\in\boldsymbol{R}^{N\times N}$を考える．

$N$個の細胞と$M$個の遺伝子のデータとして$X\in\boldsymbol{R}^{N\times M}$が与えられたとする．$\boldsymbol{R}^M$での距離をaffinityに変換するため，細胞$i$と細胞$j$について，次のadaptive Gaussian kernelを定義する：

K(xi,xj)=12π(σi+σj)exp⁡(−12(xi−xj)t(xi+xj)σi−σj).(1)K(\boldsymbol{x} _ i,\boldsymbol{x} _ j) = \frac{1}{\sqrt{2\pi(\sigma_i+\sigma_j)}} \exp\left(-\frac{1}{2}\frac{(\boldsymbol{x} _ i-\boldsymbol{x} _ j)^t(\boldsymbol{x} _ i+\boldsymbol{x} _ j)}{\sigma_i-\sigma_j}\right). \tag{1}

ただし，$\sigma_i$は，細胞$i$について，$l$番目に近い細胞までの距離である($l<k$)．

$\boldsymbol{x} _ i$ は細胞$i$での遺伝子発現を表す（行）ベクトルである．これによって，affinity行列 $K\in\boldsymbol{R}^{N\times N}$ が得られる．さらに，$K$のラプラシアン行列$T\in\boldsymbol{R}^{N\times N}$を計算する（$X$上での$k(x,y)$の$y$についての和で$k(x,y)$を割る）．$T_{ij}$は規格化されており，細胞$i$から細胞$j$へ１回で遷移できる確率を表わしている（Markov過程での遷移確率行列）．

$T$の固有ベクトルをdiffusion componentとする．さらに，対応する固有値を$\lambda_1,\lambda_2,\ldots$とする．diffusion componentの順序は固有値が降順になるようにする．

元々のグラフ$G_X$から，ノイズを減らしたグラフ$G_E\in\boldsymbol{R}^{N\times N}$を考える．pseudotimeは，このグラフ$G_E$での最短経路の長さによって計算される．また，pseudotimeを決める際，waypointを定めておく．pseudotimeはwaypointによる重み付けスコアに基づいて順番が決まる．すなわち，細胞に最も近いwaypointが，細胞の位置を決める際に最も高い重み付けを与える．多様体全体にwaypointが存在するのが望ましい（詳しくは後述）．

多様体$E\in R^{N\times L}$を考える($L<M$)．つまり，diffusion componentを$L$番目まで取ってきて，$N$個の細胞について$L$種類のdiffusion componentに沿った遺伝子発現をデータにしている．細胞$i$と細胞$j$のdiffusion distanceを次のように定義する：

DDt(ei,ej)2=∑l=1Lλl2t(ei(l)−ej(l))2.\begin{aligned} \text{DD} _ t(e_i,e_j)^2=\sum_{l=1}^L\lambda_l{}^{2t}(e_i{}^{(l)}-e_j{}^{(l)})^2. \end{aligned}

$t$はグラフに沿った（ランダムウォーキングの）ステップの数，$e_i{}^{(l)}$はdiffusion component$l$に沿った細胞$i$の埋め込みである（$l$番目の固有ベクトルを基底とした場合の成分，つまり遺伝子発現）．さらに，multi-scale distanceを次のように定義する：

MS(ei,ej)2=∑t=1∞∑l=1Lλl2t(ei(l)−ej(l))2.\begin{aligned} \text{MS}(e_i,e_j)^2=\sum_{t=1}^\infty\sum_{l=1}^L\lambda_l{}^{2t}(e_i{}^{(l)}-e_j{}^{(l)})^2. \end{aligned}

定義から，$1>\lambda_1>\lambda_2>\cdots>\lambda_L>0$なので，multi-scale distanceは次のように書くことができる：

MS(ei,ej)2=∑l=1L(λl1−λl)2(ei(l)−ej(l))2.\begin{aligned} \text{MS}(e_i,e_j)^2=\sum_{l=1}^L\left(\frac{\lambda_l}{1-\lambda_l}\right)^2(e_i{}^{(l)}-e_j{}^{(l)})^2. \end{aligned}

waypointの決定

Max-min samplingは既に存在するwaypointとの最短距離の最大値を新しいwaypointとする．これによって，多様体全体にwaypointを置くことができる．

$\boldsymbol{E}^{(l)}$を$l$番目のdiffusion componentとする．まず，$\text{WS}^{(l)}$を$N$個の細胞から適当に１つ選んだ細胞とする（初期値）．$j\in\text{WS}^{(l)}$に対し，diffusion componentに沿った，waypointまでの距離は次のように計算される：

wdij=(ei(l)−ej(l))2.\begin{aligned} \text{wd} _ {ij}=(e_i{}^{(l)}-e_j{}^{(l)})^2. \end{aligned}

細胞$i$に対し，waypoint distanceの最小値は次のように計算される：

mdi=min⁡(wdij).\begin{aligned} \text{md} _ i=\min(\text{wd} _ {ij}). \end{aligned} 先述の通り，既存のwaypointからの最短距離が最大になる点がwaypointの集合に加えられる：

WS(l)=⋃(WS(l),argmax⁡({mdi})).\begin{aligned} \text{WS}^{(l)}=\bigcup(\text{WS}^{(l)},\text{arg}\max(\{\text{md} _ i\})). \end{aligned}

満足な数のwaypointが得られるまでこれを繰り返し，これを全成分に関して行えば良い．得られた$L$個の集合の和集合を考えれば，最終的に$nW$個のwaypointの集合$\text{WS}$が得られる．

pseudotimeの計算

Palantirでは，スタートの細胞は多様体の境界にある，つまりあるdiffusion componentに関する端点であるとする．ここで，境界にある細胞を

C=⋃l=1L(argmin⁡E(l),argmax⁡E(l))\begin{aligned} C=\bigcup_{l=1}^L(\text{arg}\min\boldsymbol{E}^{(l)},\text{arg}\max\boldsymbol{E}^{(l)}) \end{aligned}

と表す．ユーザーが設定したスタートの細胞$s$に最も近いような$C$の要素をpseudotimeの開始点$s’$とする：

s′=argmin⁡i∈CMS(es,ei).\begin{aligned} s'=\text{arg}\min_{i\in C}\text{MS}(e_s,e_i). \end{aligned}

pseudotimeの初期値$\tau_i{}^{(0)}$は細胞$s’$から細胞$i$への最短経路の距離とする．これを全細胞に対して考える：${\tau_i{}^{(0)}}$．

waypoint $w$を通って，より離れた（pseudotimeが大きい）細胞$i$に行く場合，その経路の長さは，$w$までのpseudotimeと$w,i$間の距離の和になる．もし，$i$より離れたwaypointを通って$i$に行く経路という場合は，$V_{wi}$は$\tau_w$から$D_{wi}$を引けばスタートから$i$までの適切な距離が求まる．以上のことを踏まえると，$D_{wi}$を細胞$i$からのwaypoint $w$への最短経路の距離として，細胞$i$のwaypoint $w$に対するperspective（上で言った「長さ」）は次のように計算される：

Vwi={τw(0)+Dwi (τi(0)>τw(0))τw(0)−Dwi (otherwise).\begin{aligned} V_{wi}= \begin{cases} \tau_w{}^{(0)}+D_{wi}\ &(\tau_i{}^{(0)}>\tau_w{}^{(0)})\\ \tau_w{}^{(0)}-D_{wi}\ &(\text{otherwise})\\ \end{cases}. \end{aligned}

あるwaypointを経由して細胞$i$に到達する経路のperspective（長さ）を計算した．よって，細胞$i$に至る経路の長さを求めるには全てのwaypointについて平均しなければならない．$i$に近いようなwaypointを通る経路は尤もらしいので，平均を取る時に重みを大きくする．逆に，$i$から離れたwaypointを通るような経路は，平均を取る時に重みを小さくする．具体的には，waypointのperspectiveの重み付け平均は，waypointと細胞の距離に反比例するような，指数関数による重み付けを使う．そのために，重み付け行列$W\in R^{nW\times N}$を次のように定める：

Wwi=exp⁡(−Dwi2/σ)∑k=1Nexp⁡(−Dwk2/σ).\begin{aligned} W_{wi}=\frac{\exp(-{D_{wi}{}^2}/\sigma)}{\sum_{k=1}^{N}\exp(-{D_{wk}{}^2}/\sigma)}. \end{aligned}

$\sigma$は距離行列$D$の標準偏差である．重み付け平均は

τi(1)=∑w∈WSVwiWwi\begin{aligned} \tau_i{}^{(1)}=\sum_{w\in\text{WS}}V_{wi}W_{wi} \end{aligned}

と計算される．これによって，順次$\boldsymbol{\tau^{(0)}},\boldsymbol{\tau^{(1)}},\ldots$が得られ，最終的にpseudotime $\boldsymbol{\tau}$に収束する．

最終状態

waypointを結んだようなグラフ$G’_E\in G_E$ を考える．$G’_E$ は無向グラフであるが，pseudotime $\boldsymbol{\tau}$を使って有向グラフ$G_D\in\boldsymbol{R}^{N\times N}$を定義する：

GDij={GE′ijτi<τj;τi>τj, τi−τj<σi0τi>τj, τi−τj>σi.\begin{aligned} G_D{} _ {ij}= \begin{cases} G'_E{} _ {ij}\quad &\tau_i<\tau_j;\\ & \tau_i>\tau_j,\,\tau_i-\tau_j<\sigma_i\\ 0 \quad & \tau_i>\tau_j,\,\tau_i-\tau_j>\sigma_i \end{cases}. \end{aligned}

次に，$(1)$で使ったkernelを使って，affinity行列$Z\in\boldsymbol{R}^{nW\times nW}$を作る．さらに，これを

Pij=Zij∑kZik\begin{aligned} P_{ij}=\frac{Z_{ij}}{\sum_kZ_{ik}} \end{aligned}

によってMarkov過程の遷移確率行列$P$にする．$P_{ij}$は細胞が１回の遷移で状態$i$から状態$j$に到達する確率である．

最終状態は，それ以上分化しない状態なので，最終状態の集合$\text{TS}$が与えられれば，Markov過程$P$から，

Aij=0 (i∈TS,j=1,…,nW)\begin{aligned} A_{ij}=0\ (i\in\text{TS},j=1,\ldots,nW) \end{aligned}

によって吸収的Markov過程$A$を定義できる．

Markov過程の定常分布$\boldsymbol{\pi}$があるとする： $\boldsymbol{\pi}=P*\boldsymbol{\pi}$．さらに，平均絶対偏差は次のように計算される：

sc=Median(∣πi−Median(π)∣).\begin{aligned} \text{sc}=\text{Median}(\mid\pi_i-\text{Median}(\boldsymbol{\pi})\mid). \end{aligned}

これを使って，この分布における外れ値を次のように計算する：

TScands={i∣πi>Gppf(0.9999,Median(π),sc)}.\begin{aligned} \text{TS}^\text{cands}=\{i\mid\pi_i>\text{Gppf}(0.9999,\text{Median}(\boldsymbol{\pi}),\text{sc})\}. \end{aligned}

ただし，$\text{Gppf}(p,\mu,\sigma^2)$は平均$\mu$，分散$\sigma^2$の正規分布関数での$p$パーセント点である．この$\text{TS}^\text{cands}$のうち，diffusion componentの端点であるものを最終状態とする：

TS=⋂(TScands,C).\begin{aligned} \text{TS}=\bigcap(\text{TS}^\text{cands},C). \end{aligned}

分化ポテンシャル

Markov過程でのランダムウォークにより，細胞が辿る経路が分かる．それぞれの細胞について，吸収的最終状態$b$に到達する確率を表す分岐確率ベクトル$\boldsymbol{B} _ i$を求める．吸収的Markov過程Aは次のように表すことができる：

A=(QR0E)\begin{aligned} A= \begin{pmatrix} Q & R\\ 0 & E \end{pmatrix} \end{aligned}

$Q$は中間状態での遷移確率を表す$(nW-b)\times(nW-b)$行列，$R$は中間状態から最終状態への遷移確率を表す$(nW-b)\times b$行列，$E$は単位行列である．fundamental行列$F$を次のように定義する：

F∗(E−Q)=E,F=(E−Q)−1.\begin{aligned} F*(E-Q)=E,\quad F=(E-Q)^{-1}. \end{aligned}

$Q_{ij}$は中間状態$i$から中間状態$j$に遷移する確率を表すので，$F_{ij}$は中間状態$j$から（いつか）中間状態$i$に到達する確率を表す．次に，分化確率を

B=F∗R\begin{aligned} B=F*R \end{aligned}

によって計算することができる．$R _ {ij} $ は中間状態$i$から最終状態$j$（$nW-b+j$番目の状態）を表すので，$B_{ij}=\sum_kF_{ki}R_{kj}$は中間状態$i$から（いつか）最終状態$j$に到達する確率を表わしている．

$\boldsymbol{B} _ {i}$ は $\sum _ j B _ {ij}=1$ となる多項分布である．また，最終状態から最終状態への分岐確率は，明らかに

Bij={1i=j0i≠j\begin{aligned} B_{ij}= \begin{cases} 1\quad& i=j\\ 0\quad& i\neq j \end{cases} \end{aligned}

となる．これはwaypointの分岐確率なので，$W$がwaypointによる全細胞への投影の重み付けであったことを思い出せば，全細胞に関する遷移確率

Bij=∑w∈WSBwjWwi\begin{aligned} B_{ij}=\sum_{w\in\text{WS}}B_{wj}W_{wi} \end{aligned}

が得られる．

分岐確率ベクトル$\boldsymbol{B} _ i$のエントロピー

Si=−∑jBijlog⁡Bij\begin{aligned} S_i=-\sum_jB_{ij}\log B_{ij} \end{aligned}

を各状態での分化ポテンシャルとする．