Aurora blog

バイオインフォ・バイオテクノロジーにまつわる情報

シングルセル解析⑧： 細胞の運命を確率的に推定する (Plantir, CellRank)

Bioinfo Single-cell Omics

問題設計

CellRankは、Fabian TheisとDana Pe'erのグループから最近報告された、シングルセルデータから細胞が将来的に辿る状態遷移の可能性を確率的に求める手法。「細胞Aは細胞Bよりも細胞Cへと辿り着く可能性が高い」といった推定ができる。

www.nature.com

Plantir

CellRankのアイデアはDana Pe'erのグループから2019年に発表されたPlantirから発展したものだと思われる。Plantirは、細胞の状態遷移にマルコフ性がある (現在の状態のみが次の状態に影響し、過去の状態は影響しない) と仮定し、上記の課題に挑んでいる。手法の概要は以下の通り*1。

  1. Diffusion mapでscRNA-seqデータを次元圧縮
  2. Top nのDiffusion components (DCs) を使いkNN graph G_{E}を求める
  3. 各細胞のPseudo-time \tau (解析者が選んだ起点からのkNN graphにおける最小パス数) を求める
  4. Pseudo-timeをもとに G_{E}を有向グラフ G_{D}に変換する
    • \tau_{i} - \tau_{j} \leq \sigma_{i}の場合は細胞 iから細胞 jへのエッジと考える*2
  5. 解析者が選んだ終点からのEdgeを G_{D}から除外 *3
  6. G_{D}から細胞間の状態遷移確率行列 Pを求める
  7. 状態遷移確率行列 Pを使ったランダムウォーク

結果として下図dのような結果が得られる。例えば、図dのTrajectoryの(2)に位置する細胞は細胞状態A/B/Cのいずれにも同程度分化できるポテンシャルがあるが、(6)の細胞はB/Cにのみ分化する、といったような解釈ができる。

引用元：Setty et al. “Palantir characterizes cell fate continuities in human hematopoiesis” bioRxiv 385328; doi: https://doi.org/10.1101/385328. CC-BY 4.0

www.nature.com

CellRank

CellRankはマルコフ連鎖での細胞状態遷移のモデリングにRNA velocity*4を利用する。RNA velocityを使うことによって、(i) 解析者がTrajectoryの始点・終点を決める必要がなく、(ii) Pseudo-timeに依存したPlantirのように「始点から終点へ」という単一方向の状態遷移を仮定する必要がなくなる。したがって、より複雑な細胞状態遷移の過程を調べるのに適した手法だと言える。RNA velocityに関しては以前の記事を参考のこと。手法の概要は以下の通り。7-8のあたりの展開は正直理解できていない。

  1. kNN graphを作る
    • Plantirと違いDiffusion map以外の次元圧縮法を適用可能 (Default: PCA)
  2. Cell-cell adjacency matrix A \in R^{N \times N}
    • 細胞 iと jがいずれかのkNNに含まれていない場合: A_{ij} = 0
    • それ以外: A_{ij} = \rm{dist}(x_{i}, x_{j})
      • dist: デフォルトはユークリッド距離
      • x_{i}: first top n PCs of cell i (Default: top 30 PCs)
  3. Connectivity kernel P^{s} \in R^{N \times N}
    • P^{s} = row-normalized A
  4. Velocity kernel P^{v} \in R^{N \times N}
    • 下図b参照
    • P^{v}_{ij} = \frac{ \rm{exp} (\sigma c_{ij}) }{ \sum_{k \in K_{i}} \rm{exp} (\sigma c_{ik}) }
      • c_{ij} = \rm{corr} (v_{i}, s_{ij})
        • v_{i}: 細胞 iのVelocity
        • s_{ij} = x_{i} - x_{j}: 細胞 iと jの遺伝子発現の差
      • K_{i}: 細胞 iのkNN
      • \sigma = 1 / \rm{median}_{k \in K_{i}}(c_{ik})
        • VelocityがkNNのどの方向にも向いていない場合に P^{v}_{ij}を小さくする用途のパラメータ
  5. Merged kernel P \in R^{N \times N}
    • P = (1 - \lambda) P^{v} + \lambda P^{s}
      • 論文では \lambda=0.2が安定した結果を出すと記載されている
      • カーネルの設定はユーザーに委ねられている*5。ベストプラクティスが特に無いので \lambda=0.2を使うしかなさそう。
  6. Macrostate identification
    • 単一細胞レベルのデータであるカーネル Pはノイズが大きいので、細胞集団(Macrostateと呼ぶ)レベルの情報へと変換する
    • GPCCA (Generalized Perron Cluster Cluster Analysis)
      1. カーネル PをSchur分解
        • P = QRQ^{T}
        • Q \in R^{N \times N}: ユニタリ行列 ( QQ^{*} = I) 各列がMacrostateに対応する基底ベクトル
        • R \in R^{N \times N}: Pの固有値 \lambda_{1}...\lambda_{l}が対角成分にある上三角行列
      2. 以下のいずれかの基準で上位のMacrostateを選ぶ
        • Eigengap heuristics (固有値に閾値を設ける方法) *6
        • minChi algorithm
        • Crispness
      3. 選ばれたMacrostateのみを含む基底ベクトル行列を \hat{Q} \in R^{N \times n_{s}}と呼ぶ
  7. Macrostate membership \chi \in R^{N \times n_{s}}
    • \chi = \hat{Q}A: \chi_{is}は細胞 iのState sとの類似度を表す
      • A: Rotation matrix *7
    • Aは以下の f_{n_{s}}(A)を最小化する最適化問題を解き求める
      • f_{n_{s}}(A) = n_{s} - \rm{trace}( \hat{D}^{-1} \chi^{T} D \chi)
      • \hat{Q}^{T} D Q = I
      • 制約条件
        • \chi_{is} \leq 0
        • \sum_{s} \chi_{is} = 1
      • これにより \chiの各列 (各Macrostateに該当する細胞のパターン) は直交に近づく
  8. Macrostate-level transition matrix P^{c} \in R^{n_{s} \times n_{s}}
    • P^{c} = ( \chi^{T} D \chi )^{-1} (\chi^{T} D P \chi)
  9. Terminal stateを推定する
    • 分化の終末にいるMacrostateは他の細胞への状態遷移が少ないと考えられる。したがってCellRankでは遷移行列 P^{c}の情報をもとに終末にあるMacrostateを以下のように定める。
    • State i is terminal if P^{c}_{ii} > \epsilon
      • \epsilon: Threshold (Default: 0.96)
  10. Initial stateを推定する
    1. Stationary distribution \pi \in R^{N}を求める
      • \pi: 定常状態にある細胞状態組成
      • \pi^{T}P = \pi^{T}
      • \sum_{i} \pi^{T}_{i} = 1
    2. \piをMacrostate-levelに変換
      • \pi_{p} = \chi^{T} \pi
    3. Initial state
      • 定常状態を「一定以上の時間が経過した後の細胞状態組成」と仮定すると、Initial stateの細胞は定常状態では少なくなると考えられる。したがってCellRankでは以下の基準で選択する。Initial stateの数 nはユーザーが決める必要がある。
      • State i is initial if \pi_{i} is in top n smallest values
  11. Terminal cells R_{t}を定義
    • Terminal state t ごとにMembership \chi_{t}の値が高い細胞をTerminal cells R_{t}とする
  12. 各細胞の運命を推定する
    • Transient cellとRecurrent cellという2つの状態を考える。Transient cellからはRecurrent cellへ遷移できるのに対し、反対の遷移は生じないとすると、Transition matrix Pは以下のように表せる
      • このとき A=(I - Q)^{-1}Sの A_{ij}がTransient cell iがRecurrent cellのうち jへ最初に到着する確率を表すことが証明されているらしい*8
    • CellRankではTerminal cell R_{t}をRecurrent cellと見なし、その他の細胞(Transient cell)への遷移をゼロに変換。上述の Aを求め、各細胞がどのTerminal stateへ辿り着くか (Absorption probabilityと呼んでいる)を推定する
 
P = \begin{bmatrix}
\hat{P}&0\\
S&Q
\end{bmatrix}

引用元：Lange et al. “ CellRank for directed single-cell fate mapping” bioRxiv 2020.10.19.345983; doi: https://doi.org/10.1101/2020.10.19.345983 CC-BY 4.0

*1:簡略化しているので詳細は論文を要参照

*2:σは

*3:ランダムウォークの際に終点で収束するように

*4:過去記事: https://auroratummy.hatenablog.com/entry/2021/12/07/011122

*5:https://cellrank.readthedocs.io/en/stable/kernels_and_estimators.html#Initialize-a-kernel

*6:https://cellrank.readthedocs.io/en/stable/kernels_and_estimators.html#Compute-a-matrix-decomposition

*7:前述のAdjacency matrixと違うので注意

*8:Tolverらの教科書が引用されていた。http://web.math.ku.dk/noter/filer/stoknoter.pdf