scEpathの論文のメソッドのざっくりした解説，というかメモ．エネルギー最高のクラスタから単純にエネルギーが下がっていくって方針を取っているので，複数の経路を再現するとかはできない．悲しみ．

scEnergyの導入

遺伝子遺伝子間相互作用gene-gene interactionのネットワークを作る．つまり，遺伝子が$n$あるとして，そのグラフを考える．$i$番目の遺伝子と$k$番目の遺伝子が繋がっていれば$a_{ik}=1$，そうでなければ$a_{ik}=0$とする（隣接行列adjacency matrix）．詳しくは，$m$個の細胞での$i$番目の遺伝子と$k$番目の遺伝子の発現をそれぞれ，$x_i=(x_{i1},\ldots,x_{im}),x_k=(x_{k1},\ldots,x_{km})$とする．この時，

aik={1(∣cor(xi,xk)∣>τ)0(∣cor(xi,xk)∣≤τ).(1)\begin{aligned} a_{ik}= \begin{cases} 1 & (\mid\text{cor}(x_i,x_k)\mid>\tau) \\ 0 & (\mid\text{cor}(x_i,x_k)\mid\leq\tau) \end{cases}. \end{aligned}\tag{1}

カノニカル分布（正準分布）にいついて．系がエネルギー$E_i\ (i=0,1,\ldots)$を取るとき，系のエネルギーが$E_j$である確率$p_j$は

pj=exp⁡(−βEj)∑jexp⁡(−βEj)\begin{aligned} p_j=\frac{\exp(-\beta E_j)}{\sum_{j}\exp(-\beta E_j)} \end{aligned}

で与えられる．特に，分母の$\sum_j\exp(-\beta E_j)$を分配関数$Z$と呼ぶ．

scEpathでは$\beta=1$としている 1．つまり，遺伝子発現が$y$である細胞が状態$j$にある確率は，細胞が取るステージの総数，もっと言えば，分化していくときに異種とみなす細胞の数を$m$とすれば，

pj(y)=exp⁡(Ej(y))∑j=1mexp⁡(−Ej(y))(2)\begin{aligned} p_j(y)=\frac{\exp(E_j(y))}{\sum_{j=1}^{m}\exp(-E_j(y))} \end{aligned}\tag{2}

となる．ステージは適当なunsupervised clusteringでのclusterなどを使う．この場合，$m$はクラスターの数．

$E_j(y)$（遺伝子発現が$y$である$j$番目の細胞のエネルギー）については，状態$j$の細胞での$i$番目の遺伝子の発現レベルを$x_{ij}$として，標準化発現レベル$y_{ij}$を

yij=xij−xjminxjmax−xjmin\begin{aligned} y_{ij}=\frac{x_{ij}-x_j^\text{min}}{x_j^\text{max}-x_j^\text{min}} \end{aligned}

で定義する．ただし，$x_j^\text{max},x_j^\text{min}$は$j$番目の細胞の中にある全遺伝子の発現レベルのうちの最大値，最小値である．

$E_j(y)$（遺伝子発現が$y$である$j$番目の細胞のエネルギー）については，$j$番目の細胞での$i$番目の遺伝子の標準化発現レベルを$y_{ij}\ (0\leq y_{ij}\leq1)$を使って，

Ej(y)=∑inEij(y)=−∑inyijln⁡yij∑k∈N(i)ykj.\begin{aligned} E_j(y)=\sum_i^nE_{ij}(y)=-\sum_i^ny_{ij}\ln\frac{y_{ij}}{\sum_{k\in N(i)}y_{kj}}. \end{aligned}

$N(i)$は，$(1)$で考えたネットワークで，$i$番目の遺伝子の隣のある遺伝子の集合．さらに，これを標準化した

Ej^(y)=(Ej(y)/Eˉ(y))21+(Ej(y)/Eˉ(y))2\begin{aligned} \hat{E_j}(y)=\frac{(E_j(y)/\bar{E}(y))^2}{1+(E_j(y)/\bar{E}(y))^2} \end{aligned}

を使う．ただし，$\bar{E}(y)$は，全細胞（仮に$m’$個とする）についての$E_j(y)$の平均である：

Eˉ(y)=1m′∑j=1m′Ej(y).\begin{aligned} \bar{E}(y)=\frac{1}{m'}\sum_{j=1}^{m'}E_j(y). \end{aligned}

メタ細胞

各クラスタ（主成分分析を行い，２成分からなっている）2でクラスタのエネルギーの$\theta_1=80$%を占める細胞たちをメタ細胞metacellとして定義する．クラスタの数を$N$個とすれば，それぞれのメタ細胞のエネルギー$E_k^\text{M}$は，トリム平均trimeanを取る（外れ値を弾くことができる）．すなわち，$k$番目のメタ細胞のエネルギーの四部位数を下から$Q_1,Q_2,Q_3$として，

EkM=12(Q2+Q1+Q32).\begin{aligned} E_k^\text{M}=\frac{1}{2}\left(Q_2+\frac{Q_1+Q_3}{2}\right). \end{aligned}

$(2)$と同様に3すれば，ある細胞が$k$番目のメタ細胞にある確率 $p_k^\text{M}$は，

pkM=exp⁡(−EkM)∑j=1Nexp⁡(−EjM)\begin{aligned} p_k^\text{M}=\frac{\exp(-E_k^\text{M})}{\sum_{j=1}^N\exp(-E_j^\text{M})} \end{aligned}

で与えられる．ここで，$k$番目のメタ細胞に属する細胞のエネルギーは大体，$E_k^\text{M}$だと仮定している4．もちろん，メタ細胞，つまり，ある細胞が取りうる状態は$N$個あるので分配関数は$N$個の和になっている．また，ある細胞が$k$番目のメタ細胞に残る確率は$p_k^\text{M}$，それ以外のメタ細胞へ遷移する確率は$1-p_k^\text{M}$である．

Markov過程

時間的に定常なMarkov過程において，時刻が$1$経過して，$i$から$j$になる確率は，

pij=Pr⁡(Xn=j∣Xn−1=i)\begin{aligned} p_{ij}=\Pr(X_n=j\mid X_{n-1}=i) \end{aligned}

で与えられる．ここで，あるベクトル$\boldsymbol{\pi}$を考える．$\boldsymbol{\pi}$が

πj=∑iπipij,∑jπj=1\begin{aligned} \pi_j=\sum_i\pi_ip_{ij},\quad\sum_j\pi_j=1 \end{aligned}

を満たす時，$\boldsymbol{\pi}$を定常分布stationary distributionと呼ぶ．

主成分分析で得られた２次元空間で，$k$番目のメタ細胞から$l$番目のメタ細胞に遷移する確率は，これらの距離に反比例すると仮定する．メタ細胞間の隣接行列を，距離によって次のように定義する：

Gkl=exp⁡(−∥zk−zl∥24ε2).\begin{aligned} G_{kl}=\exp\left(-\frac{\|z_k-z_l\|^2}{4\varepsilon^2}\right). \end{aligned}

ただし，$z_k$は２次元空間におけるメタ細胞の中心，$\varepsilon$はメタ細胞間の距離の最大値である．$G_{kl}$を規格化すれば遷移行列が得られる：

G~klasym=Gkl∑j=1NGkj.\begin{aligned} \tilde{G}_{kl}^\text{asym}=\frac{G_{kl}}{\sum_{j=1}^NG_{kj}}. \end{aligned}

遷移行列$\tilde{G}_{kl}^\text{asym}$は非対称であるが，次に示す定常状態$\boldsymbol{\pi}_k\ (1\leq k\leq N)$を持つ：

πk=∑i=1NGkj∑k=1N∑j=1NGkj.\begin{aligned} \boldsymbol{\pi}_k=\frac{\sum_{i=1}^NG_{kj}}{\sum_{k=1}^N\sum_{j=1}^NG_{kj}}. \end{aligned}

また，$G_{kl}$が対称であることに注意すれば，容易に分かるように，

πkG~klasym=πlG~lkasym\begin{aligned} \pi_k\tilde{G}_{kl}^\text{asym}=\pi_l\tilde{G}_{lk}^\text{asym} \end{aligned}

が成立する．よって，

G~klsym=πkπlG~klasym\begin{aligned} \tilde{G}_{kl}^\text{sym}=\sqrt{\frac{\pi_k}{\pi_l}}\tilde{G}_{kl}^\text{asym} \end{aligned}

を定義すれば，これは対称な遷移行列となる．

遷移確率

メタ細胞の項の最後にも述べたことに注意すると，$k$番目のメタ細胞から$l$番目のメタ細胞に遷移する確率行列$T$は，

Tkl={(1−pkM)G~klsymk≠lpkMk=l (no transition)\begin{aligned} T_{kl}= \begin{cases} (1-p_k^\text{M})\tilde{G}_{kl}^\text{sym} & k\neq l\\ p_k^\text{M} & k=l\ (\text{no transition}) \end{cases} \end{aligned}

となる．

遷移確率行列$T$は２つのメタ細胞間がお互いに遷移する確率($T_{kl},T_{lk}$)を含む．つまり，この状態ではメタ細胞の集合は双方向性のあるグラフになる．しかし，実際に細胞が時間発展していく場合はエネルギーが減少する経路のみを取るので，グラフは一方向性なはずである．よって，Wilcoxonの順位和検定によって，２つのメタ細胞を比較する．この場合の帰無仮説$\text{H}_0$は２つのメタ細胞間のエネルギー（のトリム平均）に差がないことである．scEpathでは，p値が$\alpha=0.01$以下の場合は帰無仮説が棄却される．これらのことを考慮すれば，向き付けを考慮したメタ細胞のグラフの遷移確率行列$W$は

Wkl={(1−pkM)G~klsymk≠l,p<α,EkM>ElM;k≠l,p≥αpkMk=l (no transition)\begin{aligned} W_{kl}= \begin{cases} (1-p_k^\text{M})\tilde{G}_{kl}^\text{sym} & k\neq l, p<\alpha,E_k^\text{M}>E_l^\text{M}; \\ & k\neq l,p\geq\alpha\\ p_k^\text{M} & k=l\ (\text{no transition}) \end{cases} \end{aligned}

となる．$p$は$k$番目と$l$番目のメタ細胞を検定した時の$p$値．帰無仮説が採用された場合は，$k$番目と$l$番目のメタ細胞のエッジは双方向性となる．これによって，メタ細胞をノードとするグラフはエネルギーが増加しない遷移のみを許す有向グラフとなる．その遷移確率行列は$W$である．

scEpathを求める

細胞の実際の経路を求める場合は，最もエネルギーが高いメタ細胞から，遷移確率が最も高いような経路を辿ることになる．これは最もエネルギーが高いメタ細胞を根とし，各エッジの重みを$1-W$として構成した最小全域木minimum directed spanning tree (MDST)を見つけるのに等価である．

メタ細胞の中心をあらかじめ求めておき，その中心から$\theta_2=0.75$分位数より内側にある細胞をコア細胞とする．そして，先程求めた経路をコア細胞に限定してprincipal curveでフィッティングし，それを$[0,1]$にスケール変換する．そして，psudotime reconstruction score (PRS)を，他のデータを使って次のように定める：

PRS=c−c′c+c′.\begin{aligned} \text{PRS}=\frac{c-c'}{c+c'}. \end{aligned}

ただし，$c$は一致した細胞の数，$c^\prime$は一致しなかった細胞の数である．

pseudotime依存的なTF

pseudotime依存的な遺伝子を見つけるため，pseudotimeを10個に分割する．さらに，各部分での遺伝子発現のトリム平均を取る．各遺伝子の発現の平均を３次スプライン曲線で滑らかにする．細胞の順番をランダムに並び替えたものと比べて，標準偏差が$0.5$以上，Bonferroni調節p値が$0.01$以下のものをpseudotime依存的な遺伝子とする．

既存の転写因子(TF)の中から，pseudotime依存的な遺伝子を探し，細胞系列で発言が有意に変化している（系列に含まれるクラスター間で比べた際，Bonferroni調節p値が$0.01$以下，少なくとも閾値倍以上の変化がある）ものを，pseudotime依存的なTFとする．

さらに，TFのネットワークを構成し，Spearmanの順位相関係数を使ってp値が$10^{-5}$以下のものは相関がないとし，そうでなければ，相関があるとして結ぶ．この操作を$90$%の細胞で1000回繰り返し，どの操作でも結ばれたTFは，真に関係があると結論づける．

  1. 逆温度$\beta$を必ずしも$1$にする必要はない？ ↩

  2. これはunsupervised clusteringとかで得られた結果から行う？ ↩

  3. 正準分布の考え方．つまり，ある細胞が$k$番目のメタ細胞に属する，と言う状態になる確率を議論している． ↩

  4. この妥当性を検証するのは面白いかも ↩