はじめに

Trajectory inferenceはscRNA-seqデータから細胞状態の遷移パターンを推定する解析手法である。一般的にはこれらの手法は、細胞の遺伝子発現データの分布から、どの細胞種（細胞サブセット）が連続的な関係にあるかを推定することができるが、「どの細胞がどの細胞へ向かっているのか」といった状態遷移の方向性については推定できず、解析者の生物学的な知識で解釈する *1 必要があった。

引用元：Wikipedia: Trajectory inference. https://en.wikipedia.org/wiki/Trajectory_inference CC-BY 4.0

その中で2018年にシングルセルデータから細胞状態遷移の方向性を推定する手法RNA velocityがLa Mannoらによって提案された。後述の通り、scRNA-seqデータ以外のデータを必要とせず、転写物に含まれるイントロンに着目するだけで推定できる画期的なアイデアであったが、しばしば誤った結果が得られることもあり使いづらい印象だった。RNA velocityの現状と課題についてのレビューが最近Theis labから公開されていたので改めてまとめてみる。

La Manno et al. 2018. Nature

RNA velocityのオリジナル論文。Gieole La Manno、Sten Linnarsson、Peter Kharchenkoのグループによる報告。RNA velocityでは「イントロンが取り除かれていない転写物は転写してからの時間が間もない」→「イントロンを含む転写産物（=新しい転写産物）が多い遺伝子ほど転写が活発に行われている」→「その遺伝子の発現量は次の時点で上昇する」という仮定を置く（下図B）ことで、イントロンを含有する転写産物の量から「近い未来の遺伝子発現の変化（細胞状態遷移）の方向性」を推定する。イントロンの情報を利用する本手法は多様なscRNA-seqデータに利用できる。論文では、Smart-seq2などの転写産物の全長を解析する手法だけでなく、Chromiumなどの3'端（または5'端）のみを解析する手法においてもイントロンが検出できることが示されている（下図A）。

RNA velocityではイントロンを含む転写産物 (Unspliced mRNA) と含まない転写産物 (Spliced mRNA) について遺伝子ごとに以下のようなモデルを置く（下図b）。

$u(t)$ : 時間tにおけるUnspliced mRNAの量
$s(t)$ : 時間tにおけるSpliced mRNAの量
$\alpha(t)$ : 時間tにおける転写スピード
$\beta(t)$ : 時間tにおけるスプライシングスピード
$\gamma(t)$ : 時間tにおけるmRNA分解スピード
$\frac{du}{dt} = \alpha(t) - \beta(t)u(t)$ : Unspliced mRNA量の変化
$\frac{ds}{dt} = \beta(t)u(t) - \gamma(t)s(t)$ : Spliced mRNA量の変化

引用元：La Manno et al. “RNA velocity in single cells” bioRxiv. 2018. https://doi.org/10.1101/206052 CC-BY 4.0

RNA velocityでは以下の仮定を置くことで、さらにシンプルなモデルを考える。

転写スピード $\alpha$ とmRNA分解スピード $\gamma$ は時間によらず一定である
スプライシングスピード $\beta$ は遺伝子・時間によらず一定である
スプライシングスピード $\beta$ を $\beta = 1$ とする*2
$\frac{du}{dt} = \alpha - u(t)$ : Unspliced mRNA量の変化
$\frac{ds}{dt} = u(t) - \gamma s(t)$ : Spliced mRNA量の変化

ある遺伝子の発現量が安定した状態 (定常状態: $\frac{du}{dt} = 0$ , $\frac{ds}{dt} = 0$ ) の場合は以下の関係が成り立つ。従って定常状態にあるデータを用いればSpliced/Unspliced mRNAの量から転写スピード $\alpha$ と分解スピード $\gamma$ を推定することができる。

$\alpha = u(t)$
$u(t) = \gamma s(t)$

RNA velocityでは定常状態にあるSpliced/Unspliced mRNAのデータ $s$ , $u$ を使い最小二乗法でパラメータ $\gamma$ を求める。遺伝子によってはエクソン・イントロンの情報が完全ではないのでOffset項 $o$ を設けてその影響を考慮する。RNA velocityでは $u$ と $s$ は以下のように求められる。安定したデータを出すために近傍k個の細胞を考慮する。

$u \sim \gamma * s + o$
- $U, S$ : Unspliced/Spliced read count of the gene in k-nearest neighbor cells
- $N_{u}, N_{s}$ : Total counts of unspliced/spliced reads in k-nearest neighbor cells

ここで定常状態にあるデータをどのように見つけるか、が問題になるが、 $s$ , $u$ が最大・最小に近い細胞において (上図Dの右上・左下) に位置する細胞が定常状態 ( $\frac{du}{dt} = 0$ , $\frac{ds}{dt} = 0$ )) に近いと仮定する。RNA velocityでは $s$ , $u$ が上位 $x$ %に位置する細胞を用いて上記の回帰を行う。この手法 (Extreme quantile regression) は全データを使った回帰と比較して安定した結果が出ることを示している *3。求めたパラメータ $\gamma$ と $\alpha$ を使うことで定常状態にない細胞における各遺伝子の転写スピード $\frac{du}{dt}, \frac{ds}{dt}$ を計算できる。このように求めた転写スピード $\frac{ds}{dt}$ を使って各細胞がどの細胞状態へと近づいているかを推定する。

www.nature.com

Bergen et al. 2020. Nat Biotechnol

Fabian Theisグループから発表されたRNA velocityの派生手法。Extreme quantile regressionにおける定常状態がデータに含まれている、という仮定に従わない*4場合にも適用できる手法を考案。ここでは定常状態のデータを使った線形回帰ではなく、 $u(t), s(t)$ のカーブにフィッティングする方法をとる。

$\frac{du}{dt} = \alpha(t) - \beta(t)u(t)$
$\frac{ds}{dt} = \beta(t)u(t) - \gamma(t)s(t)$

オリジナル論文と同じ上式を変形し、本論文は以下のモデルを考える。オリジナルの手法は転写スピード $\alpha$ が時間によらず一定であると仮定していたが、このモデルでは $\alpha$ に4つのState (Induction ( $k=1$ ), Repression ( $k=0$ ), Steady state at maximum ( $k=ss_1$ ), Steady state at minimum ( $k=ss_0$ ) があり、各細胞が潜在的にいずれかのState $k$ ・時間 $t$ に該当すると考える。*5

$u(t) = u_{0}e^{-\beta\tau} + \frac{\alpha^{(k)}}{\beta} (1 - e^{-\beta\tau})$
$s(t) = s_{0}e^{-\gamma\tau} + \frac{\alpha^{(k)}}{\gamma} (1 - e^{-\gamma\tau}) + \frac{\alpha^{(k)} - \beta u_{0}}{\gamma - \beta} (e^{-\gamma\tau} - e^{-\beta\tau})$
$\tau = t - t_{0}^{(k)}$

この手法では各細胞のパラメータ ( $k$ , $t$ ) と遺伝子発現のパラメータ ( $\alpha^{(k)}$ , $\beta$ , $\gamma$ ) を推定する問題を解く。現実のデータ $u^{obs}$ と $s^{obs}$ と近い値を返すモデルのパラメータを探す。ここでは実測値との残差 $e$ を $x=(u, s)$ の符号付きユークリッド距離として、これができるだけ小さくなるような調整を行う。

$e_i = sign(s_i^{obs} - \hat{s}_{t_i}) \cdot |x_{i}^{obs} - \hat{x}_{t_i}(\theta)|$ : 細胞iの残差

ここでは残差が正規分布に従うと考え、以下の尤度 $L(\theta)$ が最大 (=negative log-likelihood $l(\theta)$ が最小) になるパラメータを探索する。

$\theta = ( \alpha^{(k)}, \beta, \gamma )$
$L(\theta) = \frac{1}{\sqrt{2\pi}\sigma} exp(\frac{1}{2n} \sum_{i}^{n} \frac{|x_{i}^{obs} - \hat{x}_{t_i}(\theta)|^{2}}{\sigma^{2}})$
$l(\theta) = \frac{1}{2} log(2\pi\sigma^{2}) + \frac{1}{2 \sigma^{2}} \frac{1}{n} \sum_{i}^{n} |x_{i}^{obs} - \hat{x}_{t_i}(\theta)|^{2}$

各細胞の状態パラメータ ( $k$ , $t$ ) と遺伝子発現のパラメータ ( $\alpha^{(k)}$ , $\beta$ , $\gamma$ ) を解析的に推定することは困難なので、本手法では以下のようなEMアルゴリズムを使う。

E step: $\theta$ を固定して $t$ を最適化する
M step: を固定して (パラメータkの最適化と同義)とを最適化する
- Nelder-Mead法で最適化をするらしい *6

引用元：Bergen et al. “Generalizing RNA velocity to transient cell states through dynamical modeling” bioRxiv. 2019. https://doi.org/10.1101/820936 CC-BY 4.0

www.nature.com

Bergen et al. 2021. Mol Sys Biol

同グループから最近発表されたレビュー論文。Fig. 2で紹介されているRNA velocityがうまく機能しない場合の事例が興味深い。以下のようなケース・遺伝子はうまく上記のモデルが適用できないらしい。

細胞種ごとに遺伝子発現のKineticsが異なる場合
- 一般性が高そうな印象。全ての細胞を使ってKineticsを推定するのではなく、ある程度状態遷移の流れを見たい細胞集団に絞って解析をする方が良いと思われる。
Transcriptional boost (ある状態時に発現が急上昇する) の場合
- 時間によらず一定と仮定されているKinetics *7 が一定でない場合はうまくいかない
安定した状態の細胞しかいない場合
- PBMCのような分化済みの細胞が中心の場合は発現変動のKineticsを推定しようとしてもランダムな結果が出やすい

引用元：Bergen et al. “RNA velocity—current challenges and future perspectives” Mol Syst Biol. 2021. doi: 10.15252/msb.202110282 CC-BY 4.0

pubmed.ncbi.nlm.nih.gov

使い方

ツールvelocyto.pyを使って、本解析のインプットであるUnspliced mRNAとSpliced mRNAのリード数をカウントする。ChromiumのデータであればCellrangerのマッピング結果(bam file)をそのまま転用することができる。リピート領域のデータはUCSC genome browserでダウンロード可能。GTFファイルはCellrangerのリファレンスゲノムに付随しているものを使うことを推奨。velocyto run10xというオプションも用意されているが、"This is an older version of cellranger"というエラーが出てきたので、開発者が推奨*8するようにvelocyto runを使用した。

velocyto run \
   --mask repeat_mask.gtf \  # リファレンスゲノムのリピート領域の情報
   --samtools-threads 4 \      # samtools sortで使うスレッド数
   --samtools-memory 16000 \    # samtools sort 1スレッドあたり使うメモリ数 (MB)
   path/to/cellranger_output \ # 入力データ Cellrangerの出力 `outs`などを含むフォルダを指定する
   path/to/gtf_file                    # gtfファイル

計算が終わるとSpliced/Unspliced mRNAのカウントデータが格納されたloom形式のファイルが得られる。これを解析するためのツールもvelocyto.pyに用意されているが、ツールの使いやすさと機能の充実度 *9 を考えるとscVeloを使う方が便利な印象。ツールの使い方については本家のTutorialが充実しているのでそちらを参照。

*1:細胞Aは未分化の細胞であるから状態遷移の起点にあるべき、というような解釈

*2:スプライシングスピードあたりの転写スピード、mRNA分解スピードを求める形になる

*3:Supplementary Note 2, Fig. 2

*4:対象の細胞種・細胞分化で発現が上昇・下降し続ける遺伝子など

*5:遺伝子発現の変動は転写調節による影響が強いはずで、オリジナル論文の仮定には明らかに無理がある、という彼らの主張だと思われる。

*6:https://www.kimoton.com/entry/20200624/1592944325が分かりやすそう

*7:Bergenらのモデルでは4状態を仮定

*8:https://github.com/velocyto-team/velocyto.py/issues/121

*9:上述のモデルをすべて使うことが可能