はじめに
この記事はkyoto.bioinfoアドベントカレンダー2022の記事です。
1細胞シーケンシング（1細胞RNA-seqや1細胞ATAC-seqなど）を行うことで、1細胞レベルでの遺伝子発現量データやクロマチンアクセシビリティデータを得ることができます。
このようなデータに基づいて、ある遺伝子の発現量をノックダウンなどにより変化させたとき、他の遺伝子の発現がどのように変化するかを予測する(in silico perturbation)手法がいくつか開発されています。
in silico perturbationを行うことで、特定の生命現象に重要な遺伝子を効率的に見つけ、検証実験につなげられる可能性があります（あくまで可能性ではありますが）。
この記事では、これらの手法がどのような仕組みで遺伝子発現変化を予測するのかを簡単に説明します。（簡潔にするため、記述が不正確な可能性があります。詳細な点については、各ツールの論文やドキュメントを読んでいただければと思います。）
今回紹介する手法は以下の二つです。
- 1細胞RNA-seqデータのみ用いる手法 (Dynamo)
- 1細胞RNA-seqデータと1細胞ATAC-seqデータを用いる手法 (SCENIC+)
これらの他に、遺伝子摂動（遺伝子発現のノックダウンや活性化）を行なった1細胞RNA-seqデータを用いる手法（CPA (Lotfollahi et al., 2021), PerturbNet (Yu et al., 2022), GEARS (Roohani et al., 2022) など）も開発されてきています。ただ、これらの手法は一般的な摂動なしの1細胞RNA-seqデータには適用できないため、今回は触れません。
Dynamo
この手法は Qiu et al., 2022 で提案されています。Dynamoは代謝ラベリングを行なった1細胞RNA-seqデータの解析を主眼に開発されていますが、代謝ラベリングなしの1細胞RNA-seqデータにも適用することができます。
Dynamoは遺伝子発現空間のベクトル場をデータから推定し、推定されたベクトル場を用いて in silico perturbation など様々な解析を行う手法です。
in silico perturbation を行う場合の流れは以下のようになります。
RNA velocityの推定
1細胞RNA-seqを行うことで、各細胞 $c$ の、遺伝子発現量ベクトル $\mathbf{x}_c $ が得られます。$\mathbf{x}_c $ は遺伝子数次元のベクトルです。 (La Manno et al., 2018)で提案された手法を用いることで、この $\mathbf{x}_c $ に加え、発現変化速度 $\frac{d\mathbf{x}_c}{dt}$（RNA velocity）もデータから推定します。
RNA velocity推定の仕組みは直感的には以下のようになります。上述の $\mathbf{x}_c $ で表される「遺伝子発現量」とは成熟mRNAの数を（おおよそ）意味します。一方で、ゲノムのイントロン領域にマップされたリードを用いることで、イントロンが除去される前のmRNA前駆体(pre-mRNA)の数を推定することができます。mRNA前駆体の数は微小時間後の成熟mRNAの数の情報を含んでいます。そのため、各細胞における成熟mRNAの数とmRNA前駆体の数のデータを用いることで、成熟mRNAの数の変化速度 $\frac{d\mathbf{x}_c}{dt}$ を計算することができます。
ベクトル場の推定
各細胞 $c$ について、遺伝子発現量 $\mathbf{x}_c $ とその変化速度 $\frac{d\mathbf{x}_c}{dt}$ が得られました。このデータから、遺伝子発現量ベクトル $\mathbf{x}$ を受け取ってその変化速度 $\frac{d\mathbf{x}}{dt}$ を返す関数 $\mathbf{f}$ （ベクトル場）を学習します。
\frac{d\mathbf{x}}{dt} = \mathbf{f}(\mathbf{x})
いま、$\mathbf{x}$ は遺伝子数次元のベクトルであり、発現変動が大きい遺伝子だけに絞っても数百~数千次元の高次元ベクトルとなります。また、$\mathbf{x}_c $、$\frac{d\mathbf{x}_c}{dt}$ には計測ノイズや推定誤差が含まれています。そのため、関数 $\mathbf{f}$ をうまく推定するために、$\mathbf{f}$ を再生核ヒルベルト空間に属す関数であると仮定し、以下のように表します。
\mathbf{f(x)} = \sum_{j=1}^{m}\Gamma(\mathbf{x}, \tilde{\mathbf{x}}_j)\mathbf{c}_j
$\tilde{\mathbf{x}}_j$ はコントロール点となる細胞の発現量ベクトルです。コントロール点となる細胞は、全細胞の中から一部（デフォルトでは全細胞の5%）をサンプリングして用います。これは全細胞をコントロール点とすると計算に時間がかかるからです。$\mathbf{c}_j$ はコントロール点ごとに与えられた遺伝子数次元の重みベクトルです。
$\Gamma$ は $w$ を幅パラメータとするガウスカーネルです。
\Gamma(\mathbf{x}, \tilde{\mathbf{x}}) = \exp(-w(\mathbf{x} - \tilde{\mathbf{x}})^2)
$\mathbf{f}$ をこのように定式化することで、関数 $\mathbf{f}$ を推定する問題は、係数 $\mathbf{c}_j$ を推定する問題になります。係数 $\mathbf{c}_j$ を推定する方法として、sparse VFCアルゴリズム (Ma et al., 2013)を用いています。この方法では、なめらかな関数 $\mathbf{f}$ を推定するために、EMアルゴリズムを用いてoutlierとなる細胞を検出して排除しつつ、$\mathbf{f}$ のノルムが小さくなるような $\mathbf{c}_j$ を推定します。
in silico perturbation
各遺伝子の発現量を微小量 ($dx_1, dx_2, \dots, dx_n$) 変化させたときの、遺伝子 $i$ の変化速度 $f_i$ の変化は以下のように表されます。
df_i = \frac{\partial f_i}{\partial x_1} dx_1 + \frac{\partial f_i}{\partial x_2} dx_2 + \cdots + \frac{\partial f_i}{\partial x_n} dx_n
これをベクトル化すると、
\begin{bmatrix}
df_1 \\
df_2 \\
\cdots \\
df_n
\end{bmatrix}
=
\begin{bmatrix}
\frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\
\frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n} \\
\cdots & \cdots & \cdots & \cdots \\
\frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \cdots & \frac{\partial f_n}{\partial x_n} \\
\end{bmatrix}
\begin{bmatrix}
dx_1 \\
dx_2 \\
\cdots \\
dx_n
\end{bmatrix}
上式右辺の行列は、ベクトル場 $\mathbf{f}$ のヤコビ行列 $\mathbf{J}$ です。無限小の変化を有限の摂動に置き換えると、上式は以下のようになります。
\Delta \mathbf{f} = \mathbf{J} \Delta \mathbf{x}
ここで、ヤコビ行列 $\mathbf{J}$ は、ガウスカーネルと係数 $\mathbf{c}_j$ を用いて以下のように計算することができます。
\mathbf{J} = \frac{\partial \mathbf{f}(\mathbf{x})}{\partial \mathbf{x}} = -2w \sum_{j=1}^{m} \Gamma(\mathbf{x}, \tilde{\mathbf{x}}_j) \mathbf{c}_j (\mathbf{x} - \tilde{\mathbf{x}}_j)^{\mathrm{T}}
これで、遺伝子発現の摂動 $\Delta \mathbf{x}$ による全遺伝子の発現変化速度の変化 $\Delta \mathbf{f}$ を計算することができます。$\Delta \mathbf{f}$ を積分することで、遺伝子発現の摂動 $\Delta \mathbf{x}$ による一定時間後の全遺伝子の発現変化を予測します。
ここまで、ベクトル場推定も in silico perturbation も遺伝子発現空間で直接行う（$\mathbf{x}$の次元が遺伝子数となる）方法を紹介しました。一方で、Dynamoでは、まず発現量ベクトルをPCAで次元削減し、低次元のPCA空間でベクトル場推定及び in silico perturbation を行い、その結果を元の遺伝子発現空間に戻す、という方法も提案されています。その場合、$\mathbf{x}$の次元はPCA空間の次元数になります。
適用例
元論文(Qiu et al., 2022)のfigure 7に、ヒト造血幹細胞の1細胞RNA-seqデータ（代謝ラベリングあり）を用いて in silico perturbation を行なった結果が載っています。in silico perturbation のやり方は Dynamoのチュートリアル に載っています。
(Sinha et al., 2022)では、トナカイの角の1細胞RNA-seqデータ（代謝ラベリングなし）を用いて in silico perturbation を行なっています (figure 6)。
SCENIC+
この手法は González-Blas et al., 2022 で提案されています（2022年12月22日現在では、bioRxivのプレプリントです）。
SCENIC+は1細胞RNA-seq（遺伝子発現）データだけでなく、1細胞ATAC-seq（クロマチンアクセシビリティ）データも用いる手法です。
1細胞で遺伝子発現量とクロマチンアクセシビリティを同時計測（10x multiome や SHARE-seq など）したデータと、遺伝子発現量とクロマチンアクセシビリティを別々に計測したデータの両方に適用することができます。ただし、後者の場合は1細胞RNA-seqデータの各細胞と1細胞ATAC-seqデータの各細胞に、共通の細胞種ラベルがついている（つまり、1細胞RNA-seqのどの細胞と1細胞ATAC-seqのどの細胞が同じ細胞種かがわかっている）必要があります。
SCENIC+は1細胞RNA-seqデータと1細胞ATAC-seqデータに基づいて遺伝子制御ネットワークを推定し、このネットワークに基づいて in silico perturbation を含む様々な解析を行う手法です。
SCENIC+の前身としてSCENIC (Aibar et al., 2017)という手法があります。SCENICは1細胞RNA-seqデータのみから遺伝子制御ネットワークを推定します。SCENIC+はSCENICを拡張し、1細胞ATAC-seqのデータも使うことで、エンハンサーを考慮して遺伝子制御ネットワークを推定する手法になります。
SCENIC+は遺伝子制御ネットワークに基づいて in silico perturbation を行います。そのため、遺伝子制御ネットワークの推定方法とネットワークに基づいた in silico perturbation の方法を順に説明します。
遺伝子制御ネットワークの推定
ここでは転写因子による遺伝子制御のみを考えます。転写因子はゲノム上のエンハンサー及びプロモーター領域に結合することで遺伝子の発現を制御します。そのため、以下の情報が得られれば、遺伝子制御ネットワークを推定できると考えられます。
- ゲノム上のエンハンサー領域及びそのエンハンサーが制御しているターゲット遺伝子
- 各エンハンサーに結合する転写因子
- エンハンサーに結合した転写因子の活性（実際にターゲット遺伝子を制御しているか）
これらの情報を、1細胞RNA-seqデータと1細胞ATAC-seqデータから推定します。
エンハンサー候補領域の推定
まず、1細胞ATAC-seqデータからエンハンサー候補領域を推定します。ここでは、細胞種間でアクセシビリティが変化している領域をエンハンサーとして推定します。
ただし、1細胞ATAC-seqデータには、データが非常にスパースであるという問題があります。これは、mRNAと違い、1細胞中にはゲノムDNAは2コピーしか存在しないことに由来します。そのため、PycisTopicという手法を用いてエンハンサー候補を推定します。
PycisTopicはLatent Dirichlet Allocation (LDA)によるトピックモデリングを行う手法です。これは、cisTopic (González-Blas et al., 2019)という既存のR実装の手法をPythonで実装し直して高速化し、かつQCなどの機能を付け足した手法です。
PycisTopicは、各細胞のトピック分布と各トピックのゲノム領域分布をそれぞれ推定します。ここでは、co-accessible（和訳不明、同時にopen/closeするという意味）なゲノム領域がトピックとして表されています。これらの確率分布を用いて、元の1細胞ATAC-seqデータの欠損値補完 (imputation) を行います。
補完されたデータを用いて、細胞種間でアクセシビリティが変動している領域 (Differentially Accessible Regions; DARs) を検出します。このDARsがエンハンサー候補となります。また、各トピックに含まれているゲノム領域もエンハンサー候補とします。
なお、トピックの数はデータから自動決定されます。トピック数は通常、細胞数よりずっと少ない数（多くて数十個）になります。1細胞レベルでそのまま解析するのではなく、少数のトピックに一度まとめてから解析することで、前述のスパース性の問題を克服する方法だと言えます。
エンハンサー候補領域への転写因子の結合の推定
PycisTopicにより、エンハンサー候補領域を推定できました。次は、PycisTargetという手法を用いて、このエンハンサー候補領域に結合しうる転写因子を推定します。PycisTargetは、データベースから取得した転写因子結合モチーフのデータを用いて、エンハンサー候補領域における結合モチーフのエンリッチメント解析を行います。
PycisTargetは、HOMER (Heinz et al., 2010) や cisTarget (Verfaillie et al., 2015) などのモチーフエンリッチメント解析手法を組み合わせた方法です。
これにより、転写因子がどのエンハンサー候補領域に結合しうるかを推定します。
エンハンサーを介した遺伝子制御ネットワークの推定
ここまでで、エンハンサー領域と、その領域に結合しうる転写因子が推定できました。
これらの情報と、PycisTopicで補完されたクロマチンアクセシビリティデータ、（今まで使っていなかった）遺伝子発現データを用いて、エンハンサーを介した遺伝子制御ネットワーク (enhancer-driven gene regulatory network; eGRN) を推定します。
まず、各エンハンサー領域が、どの遺伝子を制御しているかを推定します。各エンハンサー領域のアクセシビリティ（PycisTopicで補完された値）から各遺伝子の発現量を予測する予測器（勾配ブースティングによる回帰モデル）を作ります。この予測器を用いて、ゲノム上で遺伝子 $g$ から一定距離内に存在し、かつ遺伝子 $g$ の発現量を予測するのに重要なエンハンサー領域を推定します。これらのエンハンサー領域が、その遺伝子 $g$ を制御しているとみなします。
ここまでで、各転写因子がどのエンハンサーに結合し、各エンハンサーがどの遺伝子を制御しているかが推定できました。エンハンサーを介して、各転写因子のターゲット遺伝子候補が得られました。
最後に、各転写因子が実際にターゲット遺伝子候補を制御しているのか（活性があるか）を、遺伝子発現量データを用いて推定します。各転写因子の発現量から各遺伝子の発現量を予測する予測器を先ほどと同様にして作ります。この予測器を用いて、遺伝子 $g$ の発現量を予測するのに重要な転写因子を推定します。
この情報を用いて、先ほどの各転写因子のターゲット遺伝子候補のうち、実際に制御していそうなもののみを残し、他を消します。
これでようやく遺伝子制御ネットワークの推定ができました。
in silico perturbation
推定された遺伝子制御ネットワークを用いて、in silico perturbation を行います。
ネットワークから、遺伝子 $g$ を制御している転写因子群 $TF_1, TF_2, \dots, TF_m$ を得ます。この転写因子群 $TF_1, TF_2, \dots, TF_m$ の発現量から、（同一細胞における）遺伝子 $g$ の発現量を予測する予測器 $h_g$ を、ランダムフォレスト回帰を用いて（各遺伝子 $g$ について）作成します。
ある転写因子の発現量を変化させたときの全遺伝子の発現量変化は以下のように計算します。まず、変化させた転写因子発現量を入力として、予測器 $h_g$ で全遺伝子（転写因子含む）の発現量を予測します。これにより変化した転写因子発現量を再び入力として、予測器 $h_g$ で全遺伝子の発現量を予測します。以下これを繰り返します。
この in silico perturbation 方法自体は、（1細胞RNA-seqデータに加え）遺伝子制御ネットワークさえあればできます。そのため、SCENIC+以外の手法で推定した遺伝子制御ネットワークを用いても同様の解析を行うことができます（うまくいくかは別として）。
適用例
元論文(González-Blas et al., 2022)のfigure 4に、ヒトのメラノーマ細胞株に対する1細胞RNA-seq, 1細胞ATAC-seqデータ（別々に計測）を用いて in silico perturbation を行なった結果が載っています。in silico perturbation のやりかたはSCENIC+のチュートリアルに載っています。
まとめ
in silico perturbation を行う上では、遺伝子間の制御関係（ある遺伝子の発現を変化させた際にどの遺伝子がどのような影響を受けるか）のモデリングが重要になります。
Dynamoは発現変化速度（RNA velocity）データを用いて、各遺伝子の発現量と発現変化速度の間の関係（ベクトル場）をモデル化しています。一方、SCENIC+は遺伝子制御ネットワークを用いて、転写因子の発現量と各遺伝子の発現量の間の関係をモデル化しています。
筆者の知る限りではこれらの手法のきちんとした比較はまだなされていません。しかし、両者を同一データに適用することは可能なため、今後比較研究が現れてくるのではないかと思います。
また、これらの手法を組み合わせた方法（発現変化速度データと遺伝子制御ネットワークを両方活用する方法）も今後開発されそうです。