Biological system identification using deep generative models for high-resolution omics data

Abstract

生命活動を支配するシステムの全貌解明は、生命科学の究極的な目標の一つと言える。現在、一細胞、空間オミクスの発展により、特定の生命活動コンテキストにおける現象の網羅的な写像を得ることができるようになってきた。これらからその背後にある生体システムを同定する上で、深層生成モデルが注目を集めている。本総説では、まず、生成モデルが生体システムの同定になぜ有用なのか、並びに、その限界を深層生成モデルがどのように拡張するかを、細胞の遺伝子発現を生成する低次元の因子、細胞状態の推定問題を通して概説する。次に、細胞状態ごとの摂動に対する応答の予測、確率的ダイナミクスの推定、空間オミクスに対する微小環境の状態学習を通して、次元削減に留まらない、深層生成モデルの応用可能性を示す。最後に、さらなるオミクス技術、深層学習技術の発展が、高解像度のオミクスデータに対する深層生成モデルに何をもたらすのか、その展望を述べる。

序論

我々、ヒトをはじめとする生物の実体となる生体組織は、その構成要素たる細胞一つをとっても無数の多様な生体分子の組み合わせによって構築されている。それらが組み上げられた生体組織ともなると、無限の自由度を持つようにも思われるが、現実の生命は、ヒトから生まれる個体はヒトの形態を再構築し、疾患においても決まった時間経過で容体が変容するなど、類似のパターンを地球上のあちらこちらで再現している。このような秩序をもたらすものの一つは、遺伝子間の転写制御ネットワーク等、生体組織の構成要素を結びつける生体システムである。それゆえ、そのシステムの同定は生命の理解と制御、医療の発展において重要な意味を持っているが、そのシステムは、分子、細胞、組織、個体、集団といった様々な時空間スケール、多数の構成要素で駆動しており、それら全てを詳らかにすることはまだ現実的とは言えない段階である。

近年、この生体システムを解き明かす上で、強力な観測技術となっているのが、一細胞・空間オミクス技術である［1, 2, 3］。その中でも最も普及が進んでいる、転写産物を定量する一細胞・空間トランスクリプトーム技術は、細胞ごと、あるいは組織内の空間座標ごとに、数千から数万にも及ぶ遺伝子の発現を同時に定量することが可能となる。転写産物という、生命システムの一側面、かつ一時刻のスナップショットではあるものの、多数の遺伝子にわたり網羅的に定量が可能となったことで生命システムをデータから同定する解析技術の開発が盛んに行われている［4, 5］。

このような情報解析技術において、重要なアプローチの一つが、深層生成モデル［6］となる。深層生成モデルは、データの生成過程のモデル化とその生成過程に対する推論のそれぞれを深層学習により行うことで、複雑な生成過程やシステムの特性を柔軟にモデル化し、効率的に推論する。一細胞や空間オミクスといった複雑な生命データへの応用において、深層生成モデルの高い柔軟性は、単純な生成過程の表現だけでなく、数理生物学モデルとの融合による複雑なシステムの解析を可能にする。

本総説では、まず生成モデルがシステムの同定において有用であることを簡単な数理的記述を交えて議論する。次に、細胞の低次元表現を獲得するVariational Auto Encoder（VAE）［7］を中心に、既存の生成モデル的アプローチと比較し、深層生成モデルの有用性を明らかにする。さらに、VAEで得られる細胞状態を活用した摂動予測［8］から、そのダイナミクス［9］や微小環境状態の表現学習［10］についても触れ、深層生成モデルの応用の広がりを示す。最後に、これらのモデルの有用性と今後期待される発展について議論し、深層生成モデルの可能性を総括する。

隠れ状態推定の有用性

生成モデルは、その名前にあるようにデータの生成過程をモデル化したものである。生成モデルと対比して語られる識別モデルでは、データに対して、細胞種のラベル等、知りたいものをモデルが直接予測するが、生成モデル的アプローチでは知りたいものからデータが生成されるというモデルを立てることで、ベイズ推定の枠組みにより、データが与えられた元で知りたいものの確率分布を推定する。具体的には、データXが生成される過程を興味関心のある確率変数Z、パラメータθに依存する形P_θ（X|Z）で記述し、事前分布P（z）を設定する。そうすることで、XとZの同時分布P_θ（X, Z）＝P_θ（X|Z）P（Z）が定式化されるが、これを活用することでZの事後分布P_θ（Z|X）や周辺化尤度P_θ（X）を、解析的に、あるいは近似的に求める。前者は、Zの確率的な推定結果そのものであり、後者は、θについて最大化することで、データの実現確率を最大化する最尤推定値$\widehat{\theta }$、つまりデータと一貫した生成モデルを得るために用いることができる。この枠組みは、生体システムを確率変数Zやパラメータθによってモデル化した上で、そのシステムのもとに一細胞・空間オミクスデータが生成されるというモデルを立てることにより、データの背後にある生体システムそのものの推定をデータ駆動的に行うことを可能にする。

例えば、Piersonらは、ある細胞の遺伝子発現プロファイルxが少数の因子zの線形変換したものを平均とする確率分布から生成されるという生成モデルを仮定することで、遺伝子発現プロファイルを形作る潜在的な因子zの推定を行った［11］。具体的には、因子zは標準正規分布に従い、遺伝子発現 x′がzの線型変換Az＋μを平均、diag（（σ₁, ..., σ_D））を共分散行列とする多次元正規分布から得られると仮定し、そこにベルヌーイ分布に基づき一定の確率πで観測が0になるDrop outが入り、観測データxが生成されるという生成モデルを立てる。

  

$$\begin{array}{c}P(z)=N(z|0,I)\\ P_{A,\mu ,\sigma }(x\text{'}|z)=N\left(x\text{'}|\mathit{Az}+\mu ,\mathit{diag}(\sigma )\right)\\ P_{\pi }(x|x\text{'})=\mathit{\pi \delta }(x)+(1-\pi )\delta (x-x\text{'})\end{array}$$

この場合は、μ, A, σ, πといったモデルのパラメータが与えられた下で確率変数 x′, zの事後分布P（x′, z|x）を解析的に求めることができる。その一方で、パラメータの周辺化対数尤度を解析的に導出することは出来ないため、パラメータの最適解を得ることが問題となる。ここで、x′, zの事後分布で、x, x′, zの対数同時分布を周辺化すれば、パラメータの最尤推定を可能とするEMアルゴリズム（期待値最大化法）［12］のQ関数を求めることができる。

EMアルゴリズムでは、このQ関数の計算とパラメータについてのQ関数の最大化を行うことで、対数尤度の極大値をとるパラメータを得る。これにより、データの生成確率が高い生成モデルの元でのx′, zの事後分布を求めることができる。今回の例では、細胞ごとの発現データxから、遺伝子発現プロファイルに対応する低次元の因子z、因子zと各遺伝子の関係性を表現するAをデータと一貫した形で求めることができる。このように、生成モデルにおいては、データの生成過程に関わると仮定したシステムそのものをデータ駆動的に推定することができる。

その一方で、これまでの生成モデルの最適化アプローチでは、対数同時分布log P（X, Z）が、Zについてシンプルな形であることを要求するものが多い。例えば、上述のEMアルゴリズムでは、潜在変数zの事後分布による対数同時分布の期待値計算を行う。これを効率的に計算するには、解析的な導出が望ましいが、対数同時分布のzに対する依存が複雑な場合には、サンプリングに基づく数値計算が必要となる。また、事後分布からのサンプリングを可能とするマルコフ連鎖モンテカルロ法［13］は、そのような制約は必要としないものの、パラメータや確率変数が増大するに従い適切な確率変数のサンプリングが困難となる。

そのため、たとえば、今回のような遺伝子発現 xに対応する潜在因子zの推定においても、遺伝子発現への変換を線型写像と仮定する等、モデルの単純化が必要となり、たとえばニューラルネットワーク等の非線形関数によって潜在因子zが遺伝子発現 xに対応づけられるといった複雑なモデルとなると、既存の方法では最適化、並びに事後分布の計算が困難になってしまう場合が多く存在した。このような課題を解決し、深層学習モデルを含む、柔軟な生成モデル設計を可能とする理論的枠組みが、深層生成モデル、特にVAEにおいて導入された変分推論を応用した枠組みとなる。

VAEによる細胞状態の学習

深層生成モデル自体は、Generative Adversarial Network等を含むより広範な概念である［6］が、本総説では一細胞・空間オミクスへの応用で積極的に用いられている変分推論を活用したアプローチについて紹介する。特に、ここではこのアプローチの原点とも言えるVAEによる細胞状態の学習アルゴリズム［14］の紹介を通して、その柔軟性と応用可能性を示す。本総説では、生成モデルの観点からVAEを解説するが、同じくJSBi Bioinfomatics Reviewの大野による総説では、より実装に即した形で解説がなされているため、そちらも参照されたい［15］。

VAEは、自己符号化器を冠するその名の通り、高次元データに対する自己回帰的な低次元表現の獲得を生成モデルの償却変分推論の枠組みにより行う（図1）。まず、遺伝子gの発現 x_gが、細胞の低次元因子zのNeural Network（NN）による変換 λ_θ（z）_gを平均とするPoisson分布から観測され、zは標準正規分布に従うと仮定する。

  

$$\begin{array}{c}P(z)=N(z|0,I)\\ P_{\theta }(x|z)={\mathrm{\Pi }}_{\mathrm{g}}\text{Poisson}\left(x_g|{\lambda }_{\theta }{(z)}_g\right)\end{array}$$

図１：非観測の低次元因子（細胞状態） zによる非線形変換をもとに高次元の遺伝子発現の確率分布が決まる生成モデル。解析的に事後分布P(z|x)を求めることができないが、償却変分推論では、xを入力とする深層学習モデルにより実装された変分事後分布で近似する。

この生成モデルは、前章において最適化が困難であると論じたzとxの間の依存関係が複雑な生成モデルとなっている。しかし、以下に述べる変分推論では、対数尤度の下限値値の最大化を通して、最尤推定の近似値とその近似値の下での事後分布を得ることができる。まず、任意の確率分布q(z)を導入すると、対数尤度の下限値がJensenの不等式により以下のように導ける。

  

$$\begin{array}{c}\log P_{\theta }(x)=\log \int P_{\theta }(x|z)P(z)\mathit{dz}=\log \int P_{\theta }(x|z)P(z)\frac{q(z)}{q(z)}\mathit{dz}\\ \ge \int q(z)\log \frac{P_{\theta }(x|z)P(z)}{q(z)}\mathit{dz}={𝔼}_{q(z)}[\log P_{\theta }(x|z)]-D_{\mathit{KL}}\left(q(z)|P(z)\right)\end{array}$$（1）

この下限は、Evidence Lower Bound（ELBO）という変分推論における目的関数となる。第一項は、元の観測 xをzの非線形変換 λ(z)からどれほど観測されそうか、つまり、λ_θ(z)の再構成性能を定量するものとなっており、第二項はq(z)と標準正規分布のKL divergenceとなっておりq(z)が可能な限り単純な分布となるような正則化項として働く。さらにELBOは、以下のように変形するともう一つの観点から理解することもできる。

  

$${𝔼}_{q(z)}[\log P_{\theta }(x|z)]-D_{\mathit{KL}}\left(q(z)|P(z)\right)=\log P_{\theta }(x)-D_{\mathit{KL}}\left(q(z)|P_{\theta }(z|x)\right)$$

このように、ELBOは、周辺化された対数尤度と生成モデルの事後分布とq(z)の負のKL divergenceの和としても捉えることができる。そのため、θについてELBOを最大化すれば、周辺化対数尤度が向上したより現実のデータ分布に適合した生成モデルとなる一方で、q(z)については事後分布と一致するときにKL Divergenceが0となるため、q(z)＝P_θ(z|x)がELBOを最大化するq(z)となる。これらのことから、ELBOの最大化は、生成モデルのデータへの適合化と事後分布P(z|x)に一致するq(z)の獲得という、当初の目的を達成するものになっていることがわかる。

VAEでは、xを入力とし、変分パラメータφによりパラメトライズされるNNの出力${\mu }_{\varphi }(z),{\sigma }_{\varphi }^2(z)$を平均、分散とする正規分布$N\left(z|{\mu }_{\varphi }(z),\mathit{diag}\left({\sigma }_{\varphi }^2(z)\right)\right)$によりqを実装する。これにより、ニューラルネットワークに十分な表現力があり、データも十分にある場合には、正規分布の中で最もP(z|x)に類似したq_φ(z|x)を得ることができる。

実際の最適化においては、式（1）の第二項の標準正規分布とq(z|x)のKL divergenceは解析的に計算でき、深層学習の自動微分を活用した最適化が可能となる一方で、第一項の再構成項はqについての期待値計算が含まれており解析的な計算ができずθ，φについての最適化を行うことができない。そこで、第一項を以下のようにモンテカルロ近似により近似計算を行う。

  

$${𝔼}_{q_{\varphi }(z|x)}[\log P_{\theta }(x|z)]\mathrm{\approx }\frac{1}{L}\sum _{l=1}^L\log P_{\theta }(x|\overline{z_l})$$

ここで、Lはモンテカルロ近似の際のサンプリング回数であり¹、$\overline{z_l}\sim q_{\varphi }(z|x)$であるが、${ϵ}_l\sim \mathit{𝒩}(0,I)$として$\overline{z_l}={\mu }_{\varphi }(x)+\text{diag}\left({\sigma }_{\varphi }(x)\right){ϵ}_l$と計算することで、φによる自動微分を実行可能とする。このようにパラメータ依存の確率分布からのサンプリングをパラメータ依存の決定的な変換とパラメータ非依存的な確率的サンプリングに分解することで自動微分を可能とする技術をリパラメトライゼーショントリック［7］といい、深層生成モデルを様々な確率分布へと応用可能とする技術である。最終的にVAEでは以下のような目的関数の確率的勾配降下法を用いた最小化により$\widehat{\theta }$，$\widehat{\varphi }$を得る。

  

$${\mathit{𝐿}}_{\theta ,\varphi }(x)=-\frac{1}{L}\sum _{l=1}^L\log P_{\theta }(x|\overline{z_l})+D_{\mathit{KL}}\left(q_{\varphi }(z|x)|P(z)\right)$$

この最適化結果を用いることにより、$\widehat{z}={\widehat{\mu }}_{\varphi }(x)$により、遺伝子発現の再構築を可能な低次元表現zを得るほか、観測発現量xの観測ノイズによるスパース性を補完したような再構成結果$\widehat{x}={\lambda }_{\widehat{\theta }}(\widehat{z})$を得ることが可能となる。このようにVAEでは、低次元表現zからニューラルネットワークによる変換により遺伝子発現プロファイル${\lambda }_{\widehat{\theta }}(\widehat{z})$が生成されるという、複雑な生成過程を仮定したにも関わらず、xに応じて柔軟に変化する変分事後分布をNNで実装したことで、モデルの最適化と事後分布の推定を同時に実現している。

実は、変分事後分布はELBO最大化の観点だけから言えば、q(z)を観測データ点xごとに独立の平均パラメータと分散パラメータを設定し、最適化すれば同様にP_θ(z|x)に最も近い正規分布を得ることができる（同時にその下でのθの最適化も行うことができる）。しかし、この方法では、サンプルサイズに応じて最適化対象が増大してしまい、新規のデータに対するzを推定する際には最適化を再度実行しないといけない上、深層学習の最適化で一般的に有効なミニバッチを使用した学習では各ステップで大多数のパラメータが更新されないという問題がある。

その一方で、q(z)をNNによりx依存的にモデル化することにより、xに応じてP(z|x)を柔軟に近似することを可能にしつつ、NNを使い回す（償却）ことで新規のデータに対するzの前向きな推論やミニバッチ学習を使用した最適化を実現可能にしている。このような枠組みを償却変分推論［16］と言い、この後紹介するような一細胞、空間オミクスへの応用を含め、VAEを超えたより広範な確率モデルで有効なアプローチとなっている。以降の章では、償却変分推論を活用した一細胞・空間オミクス解析手法を複数紹介し、その応用性の高さを示す。

細胞状態を介した摂動応答の予測

VAEの重要な性質の一つに潜在次元の各次元が別々の情報を表現するdisentanglementがある［17］。これは、（1）の第二項の標準正規分布と変分事後分布のKL divergenceにより、標準正規分布が各次元独立の正規分布であることから、q(z|x)についても各次元が独立となりやすくなるような正則化を受けるためである。

このような性質を、摂動に対する遺伝子発現の応答予測に応用したものが、scGen［8］である。特定の刺激が細胞にもたらされると、細胞内の酵素や転写因子の存在により多段階のリン酸化パスウェイを介した発現制御や、遺伝子制御ネットワークによる再帰的な制御により、非線形な遺伝子発現の変動が引き起こされる。そのため、刺激に対する遺伝子発現の応答は、細胞集団によって異なるものとなり、応答が観測されていない細胞集団に対する遺伝子発現の変動を予測することは困難なタスクとなる。

Lotfollahiらは、このような課題をVAEにより学習される潜在空間を使って解決することを試みた。具体的には、細胞iの遺伝子発現をx_i、摂動の有無をc_i、細胞クラスタをk_iとするような一細胞トランスクリプトームデータを考える。ただし、c＝0,1はそれぞれ摂動の無し、有り、に対応し、最適化時にPerturbationを行った集団を使用しない細胞種k_w/ooptの細胞については、全ての細胞がc＝0となる。まず、X＝｛x_i|i＝1, ..., N｝を使用して、VAEの訓練を行い細胞ごとの潜在状態を推定する。ここから、c＝0の細胞のみを用いて条件c＝0,1のそれぞれについて最大の細胞種の個数$N_{\mathit{kmax}}^c$に合わせて全ての細胞のアップサンプリングを行った細胞集団に対して条件cにおける平均的な細胞状態ベクトルを算出し、その差分により、摂動の有無による細胞状態の変化を求めた。この差分を使用し、perturbationが存在しない細胞クラスタに属する細胞のperturbationをした場合の細胞状態、さらにはVAEのデコーダーを使用することで、perturbation時の遺伝子発現をシミュレートする（図2）。

図２：特定の細胞集団（青）のみに摂動後の発現プロファイルが得られている場合、遺伝子発現空間での摂動による差分により、他の集団における摂動を予測することは難しい。scGenでは、細胞状態の空間における差分により、摂動後の細胞状態と発現プロファイルがより良い推定となることを示した。

Lotfollahiらは、このようなアプローチが、遺伝子発現空間における差分を使用したperturbation予測と比べて、摂動時の発現プロファイルを学習に使用していない細胞集団におけるperturbation時の発現プロファイルの予測をより高精度に行えることを示した²。scGenの成功は、遺伝子発現空間では非線形の変動となるような発現変化が潜在空間内では並行移動として表せることを示唆している。これは、VAEの潜在空間がdisentanglementされていることにより、何らかの要因によって引き起こされる遺伝子発現の変動が潜在空間の特定の次元によって表現されるためであると考えられる。さらに、VAEが備えるEncoderとDecoderは、遺伝子発現空間と細胞状態空間の双方向的な変換を可能とし、scGenは、細胞状態空間内での変化を遺伝子発現の空間に落としこむことに成功した。これらのVAE、並びに、その潜在空間の性質は、以降の章で紹介するように摂動による発現変動に留まらないさまざまな生命現象へと応用可能なものとなっている。

細胞状態ダイナミクスの不確実性

一細胞トランスクリプトームによって、発生過程から、がんをはじめとする疾患等に渡る様々な生物学的な系で多様な細胞の状態が網羅的な発現プロファイルとともに明らかになりつつある［19, 20］。それでは、このような多様な状態はどのようにして生じるのだろうか。ここに、一細胞トランスクリプトームをはじめとするオミクス解析技術の限界がある。トランスクリプトームという転写産物の総体を得るためには、対象となる細胞を破壊する必要がある。そのため、基本的には細胞が破壊されるその時点での発現プロファイルのみが得られるため、それがどのような変化し、現在の多様性を生み出しているかを明らかにすることは容易ではない。

この課題に情報解析技術によってアプローチしたものが、RNA velocityである［21］。著者のLa Mannoらは、一細胞トランスリプトームとして観測された転写産物を、Intronの有無に基づいて未成熟RNA uと成熟RNA sとに分離して定量した。この二つは、反応速度論に基づくスプライシングの方程式により、成熟RNAの時間変化と結びつけることができる。

  

$$\frac{\mathit{ds}}{\mathit{dt}}=\mathit{\beta u}-\mathit{\gamma s}$$

この方程式は、成熟RNAの時間変化が、未成熟RNAのスプライシングに起因する増大と成熟RNAの分解に起因する減少により決定されるという数理モデルになっており、β, γはそれぞれスプライシングと分解の速度を示すパラメータとなる。RNA velocityは、この数理モデルを活用することで、一細胞ごとの未成熟、成熟RNAの定量から発現変化の速度を予測する。さらに、予測された発現変化パターンと細胞状態間の発現パターンの差分の類似性に基づき、それぞれの細胞状態がどの細胞状態に遷移するか、つまりは細胞分化の流れを推定する。

その後も多数の改良手法［22, 23, 24］が提案されていたものの、いずれの手法においても発現変化の推定は、遺伝子ごと独立に決定論的な推定を行うものとなっており、細胞状態そのものがどのようにして多様化していくのか、その原点を辿ることは依然として容易ではなかった。そこで、筆者らは、VAEにより推定される細胞状態そのものの確率的なダイナミクスを推定することで、その不確実性から細胞状態が多様化する起点を同定、探索するための方法論VICDYFを開発した［9］（図3）。

図３：細胞状態のダイナミクスを確率分布で捉える。細胞状態の分岐を不確実性の高まりとして検出する。VICDYFでは、成熟RNAを反映した細胞状態をVAEにより学習した上で、細胞状態空間における変化を変分事後分布からサンプリング、遺伝子発現の時間変化に変換し、スプライシング数理モデルにより未成熟RNAの平均発現量を推定する。

具体的には、まず、VAEの潜在空間において、細胞状態zとともに、微小時間δtの間の細胞状態変化dを仮定する。

  

$$\begin{array}{c}P(z)=\mathit{𝒩}(z|0,I)\\ P(d)=\mathit{𝒩}(d|0,0.01I)\end{array}$$

ここで、dの分散パラメータを0.01としているのは、微小時間における変化がzの多様性と比べて十分に小さくなるようにするためである。このz, dを観測データs, uと結びつけることにより、細胞状態zとそのダイナミクスdの変分事後分布を推定する。

まず、遺伝子g成熟RNA s_gは、細胞状態zをNNにより変換したλ_θ(z)_gを平均パラメータとするPoisson分布に従うと仮定する。

  

$$P(s|z)={\mathrm{\Pi }}_{\mathrm{g}}\text{Poisson}\left(s_g|{\lambda }_{\theta }{(z)}_g\right)$$

この定式化は、細胞状態が現在の成熟RNA sの発現プロファイルを反映したものとなることを仮定したものとなる。

次に、dをどのように観測と結びつけるかであるが、ここで、デコーダーλ_θ(z)とスプライシングの方程式を利用する。

まず、dが与えられた元では微小時間前後の成熟RNAの発現プロファイルがλ(z－d)，λ(z＋d)により計算できるため、成熟RNAの時間変化ds/dtを微小時間前後の有限差分により近似することができる。

  

$$\frac{\mathit{ds}}{\mathit{dt}}\mathrm{\approx }\mathit{\delta s}=\frac{{\lambda }_{\theta }(z+d)-{\lambda }_{\theta }(z-d)}{2\mathit{\delta t}}$$

次に、スプライシングの数理モデルをもとにして、uの平均パラメータλ_uを以下のようにモデル化する。

  

$$\begin{array}{c}\mathit{\delta s}={\mathit{\beta \lambda }}_u-{\mathit{\gamma \lambda }}_{\theta }(z)\\ \Rightarrow {\lambda }_u=\frac{\mathit{\delta s}+{\mathit{\gamma \lambda }}_{\theta }(z)}{\beta }\end{array}$$

遺伝子gの未成熟RNA u_gはこのλ_u(z,d)_gを平均とするポアソン分布に従うと仮定する。

  

$$P(u|z,d)={\mathrm{\Pi }}_{\mathrm{g}}\text{Poisson}\left({\lambda }_u{(z,d)}_g\right)$$

zの確率的な変化dはスプライシング数理モデルのもとに、観測データの生成に関わる形で定式化された。VICDYFでは、VAEと同様に償却変分推論を用いて、ELBOの最大化によりデコーダーパラメータの最適化とデータ、そして数理モデルに一貫したz, dの変分事後分布 $q_{\varphi }(z|x)=N\left(z|{\mu }_{\varphi }^z(x),\mathit{diag}\left({\sigma }_{\varphi }^z(x)\right)\right),q_{\varphi }(d|z)=N\left(z|{\mu }_{\varphi }^d(z),\mathit{diag}\left({\sigma }_{\varphi }^d(z)\right)\right)$の推定を行う［9］。ここで、q_φ(d|z)は、細胞状態zにおけるダイナミクスを確率分布として推定したものとなっており、平均的な変化${\mu }_{\varphi }^d(z)$だけではなく、その不確実性${\sigma }_{\varphi }^d(z)$を評価することが可能となる。

筆者らは、このVICDYF法により、造血幹細胞の分化の経路に応用しダイナミクスの不確実性が大きい細胞状態を抽出することで、B細胞と形質樹状細胞の分岐が起こる細胞状態の特定とその細胞状態ではLFA-1の発現の多寡が形質樹状細胞への分化傾向と相関することを見出した［9］。このように、VICDYFでは深層生成モデルにより一細胞トランスクリプトームデータから推定される細胞状態空間に時間変化の概念を取り入れることで、細胞分化の不確実性を推定する方法論を構築した。

空間トランスクリプトームに対する微小環境表現推定

ここまで、一細胞トランスクリプトームに対する深層生成モデルの応用により、遺伝子発現空間とNNによって行き来できる細胞状態空間を中心として、クラスタリング等に有用な低次元表現の獲得にとどまらず、特定刺激に対する細胞の応答や、確率的な時間変化など、様々な生命現象をデータと一貫した形で推論できることを示してきた。その一方で、最近の空間オミクスの発展により、細胞レベルの遺伝子発現は、空間的な座標も伴って観測できるようになってきている［3, 25, 26］。それでは、このような空間構造を伴う遺伝子発現を考慮した時に、データの背後にどのような潜在状態を仮定するべきであろうか。

その上で、有用になるのが、微小環境、またはニッチの概念である。例えば、腫瘍が、発癌、転移等、その病態を進行させる際には、その進行はがん細胞のみによって引き起こされるものではなく、複数の細胞種が織りなす細胞間の相互作用のもとに駆動されている［27］。そのため、これらの現象を理解するためには、細胞単体の多様性をそれぞれ別々に捉えるのではなく、複数種類の細胞の組み合わせによって構成されるローカルな共同体全体、微小環境の違いや変化を捉える必要がある。このような微小環境の状態を推定するためにも、深層生成モデルが用いられつつある［10, 28］。

Birkらは、一細胞レベルの空間オミクスデータから空間近傍グラフを構築し、グラフ生成を活用することで、空間一貫性と遺伝子発現に応じた多様性を両立した微小環境状態を推定する手法、NicheCompassの開発を行った［10］。NicheCompassでは、細胞ごとの発現プロファイルx_i(i＝1, ..., N)とその空間近傍の平均をとったx′_iの他に、その空間座標をもとに作られる空間近傍グラフを入力にする。

  

$$A_{\mathit{ij}}=\{\begin{array}{cc}1& \left(j\in \mathrm{ℰ}(i,k)\right)\\ 0& \text{otherwise}\end{array}$$

ただし、ℰ(i, k)は、細胞iの空間的に近くからk番目までの細胞群である。このようなデータに対して、VAEと同様に標準正規分布に従う潜在状態zを仮定する。

  

$$P(z_i)=\mathit{𝒩}(z_i|0,I)$$

この潜在状態からx_i, x′_iが生成されることを仮定する。

  

$$\begin{array}{c}P(x_i|{\mathrm{z}}_{\mathrm{i}})={\mathrm{\Pi }}_{\mathrm{g}}\text{NegativeBinomial}\left(x_{\mathit{ig}}|\mu {(z_i)}_g,{\alpha }_g\right)\hfill \\ P({x\text{'}}_i|{\mathrm{z}}_{\mathrm{i}})={\mathrm{\Pi }}_{\mathrm{g}}\text{NegativeBinomial}\left({x\text{'}}_{\mathit{ig}}|\mu \text{'}{(z_i)}_g,{\alpha \text{'}}_g\right)\hfill \end{array}$$

これにより、その細胞、さらにはその近傍の発現プロファイルに一貫したzとなるような過程となっている。ただし、一般的な深層生成モデルと異なる点は、μ(z_i), μ′(z_i)にはGene Programを活用した線形モデルを使用しており、z_iの各次元がGene Programのactivationに対応するなど、解釈性の向上が図られている。

  

$$\mathrm{\mu }({\mathrm{z}}_{\mathrm{i}})=\text{softmax}(P\circ W{z}_i)$$

ただし、∘はアダマール積（要素ごとのせき）を示し、Pは各遺伝子の遺伝子とプログラムの所属関係を表す行列、Wは、Gene Programの所属遺伝子への影響を表す重みとなっている（μ′(z_i)も同様）。Pはデータベースから構築される一方、Wは最適化可能なパラメータとなっている。さらに、Nichecompassのユニークな点は、この潜在状態zから隣接行列が生成されるとすることである。具体的には、以下のようなベルヌーイ分布により各Edgeの有無が決まることを仮定している。

  

$$P(A_{\mathit{ij}}|z_i,z_j)=\text{Bernoulli}\left(\text{cosine\_similarity}({\mathrm{z}}_{\mathrm{i}},z_j)\right)$$

これにより、A_ijが1となるz_i, z_jは類似したものになる、つまりは空間的に近傍な細胞の潜在状態が一貫したものになるという性質が与えられる。NicheCompassでは、Encoderにおいてもこの空間隣接行列Aを活用したGraph Attention Layer［29］を使用することで、空間の近接性を反映した潜在状態の変分推論を行なっている。

Brikらは、NicheCompassの有用性を空間的一貫性、近傍関係の保存性、遺伝子発現に対する忠実性等、様々な尺度で評価し、既存手法に対する優位性を示したほか、マウス胚発生、乳がん組織、マウス全脳など多様なデータセットに対して適用し、ニッチの同定やGene Programの空間パターンとそれに紐づく細胞間のクロストークを明らかにした［10］。このように、深層生成モデルによる潜在状態推定は、細胞の状態のみならず、微小環境状態をはじめ生体組織の様々な階層に対して応用が可能であり、償却変分事後分布のスケーラビリティと柔軟性により、大規模データに対しても解釈可能性の高い生成モデルの仮定と最適化による知見抽出を可能とする。

深層生成モデルの有用性と今後の展望

本総説では、深層生成モデルというアプローチが一細胞・空間オミクスデータの背後に存在する生体システムを同定する上で、いかに有用であるか、具体的な応用例を定式化とともに示してきた。その有用性は、一般には、推定対象を観測データと結びつける際の自由度の高さ、そして、償却変分推論による潜在状態の効率的、かつ柔軟な推定によるところが大きい。

その一方で、一細胞・空間オミクスの分野において、深層生成モデルが他の深層学習アプローチと比べて大きな成功を収める背景としては、観測のノイズが大きく推定対象たる生体システムの状態と現実の観測データの間の不確実性が大きいと想定されること、細胞状態やそのダイナミクスに代表されるように推定の対象が訓練データとして与えられない教師なしの枠組みに該当する場合が多いためと考えられる。また、一細胞・空間オミクスという組織内の多様性を捉える上で強力な観測を活用するにあたって、その多様性、つまりは細胞状態、微小環境等の潜在状態、によって条件づけられた生体システムの推定は、その多様性によって生体システムがどう変容するか、というこれまでの生物学では解析が困難であった新たな次元を切り開いている。

一細胞・空間オミクスの深層生成モデルには、まだまだ数多くの発展の余地が残されていると考えている。一つには、配列モデルの取り込みである。近年、タンパク質配列をはじめ、DNA、RNAといった核酸配列においても大規模言語モデルを応用した基盤モデルの創出が相次いでおり［30, 31, 32］、配列分子の機能推定等での有用性が証明されつつある。これらの配列情報と、細胞状態等、一細胞空間オミクスで推定される潜在状態のコンテキストとの関連性は、転写制御をはじめとした様々なシステムにおいて重要な意味を持つと考えられる。

また、近年では大きな規模なデータセットで訓練された一細胞・空間オミクスの基盤モデルも提供され、下流のタスクで優れた精度を示している［33, 34］。これらの汎用的な細胞状態をもとにした、推論を行うことができれば、データセットを跨いだ遺伝子間制御の比較等、新たな展開を迎えることが期待される。一細胞・空間オミクスは、その観測技術は未だ発展を続けている上、急速な普及に伴いデータセットの大規模化も進んでおり、個体間での細胞群、微小環境の違いも議論できるようになりつつある［35, 36］。今後は、このような観測技術に加え、分子配列解析等の周辺分野の進歩により、分子から個体、集団に渡る生体システムの有り様とその多様性が探求できるようになる未来が期待される。

謝辞

本論文の執筆にあたり、原稿や執筆に対する助言をいただいた、国立がん研究センターの廣瀬遥香博士と西村和也博士に謝意を述べたい。

脚注

1　VAEの場合、L＝1が設定されることが多い。

2　ただし、Liらによる複数のデータ、複数の手法を対象とした性能検証［18］においては、scGenを含む深層学習ベースの手法は、Perturbationの有無による平均値の違いを用いた場合に比べて低い性能を示しており、深層生成モデルによる非線形な応答の再現が有効になる場面は限定的とする議論もある。

References

• [1]  Aldridge S, Teichmann SA. Single cell transcriptomics comes of age. Nat Commun. 2020; 11: 4307. doi:10.1038/s41467-020-18158-5
• [2]  Fischer DS, Villanueva MA, Winter PS, Shalek AK. Adapting systems biology to address the complexity of human disease in the single-cell era. Nat Rev Genet. 2025; 1-18. doi:10.1038/s41576-025-00821-6
• [3]  Moses L, Pachter L. Museum of spatial transcriptomics. Nat Methods. 2022; 19: 534-546. doi:10.1038/s41592-022-01409-2
• [4]  Heumos L, Schaar AC, Lance C, Litinetskaya A, Drost F, Zappia L, et al. Best practices for single-cell analysis across modalities. Nat Rev Genet. 2023; 24: 550-572. doi:10.1038/s41576-023-00586-w
• [5]  Kleino I, Frolovaitė P, Suomi T, Elo LL. Computational solutions for spatial transcriptomics. Comput Struct Biotechnol J. 2022; 20: 4870-4884. doi:10.1016/j.csbj.2022.08.043
• [6]  Bond-Taylor S, Leach A, Long Y, Willcocks CG. Deep generative modelling: A comparative review of VAEs, GANs, normalizing flows, energy-based and autoregressive models. arXiv [cs.LG]. 2021. Available: http://arxiv.org/abs/2103.04922
• [7]  Kingma DP, Welling M. Auto-Encoding Variational Bayes. arXiv [stat.ML]. 2013. Available: http://arxiv.org/abs/1312.6114
• [8]  Lotfollahi M, Wolf FA, Theis FJ. scGen predicts single-cell perturbation responses. Nat Methods. 2019; 16: 715-721. doi:10.1038/s41592-019-0494-8
• [9]  Nagaharu K, Kojima Y, Hirose H, Minoura K, Hinohara K, Minami H, et al. A bifurcation concept for B-lymphoid/plasmacytoid dendritic cells with largely fluctuating transcriptome dynamics. Cell Rep. 2022; 40: 111260. doi:10.1016/j.celrep.2022.111260
• [10]  Birk S, Bonafonte-Pardàs I, Feriz AM, Boxall A, Agirre E, Memi F, et al. Large-scale characterization of cell niches in spatial atlases using bio-inspired graph learning. bioRxiv. 2024. p. 2024. 02. 21. 581428. doi:10.1101/2024.02.21.581428
• [11]  Pierson E, Yau C. ZIFA: Dimensionality reduction for zero-inflated single-cell gene expression analysis. Genome Biol. 2015; 16: 241. doi:10.1186/s13059-015-0805-z
• [12]  Dempster AP, Laird NM, Rubin DB. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc. 1977; 39: 1-38. Available: https://www.jstor.org/stable/2984875
• [13]  Robert C, Casella G. Monte Carlo Statistical Methods. 2nd ed. New York, NY: Springer; 2005. Available: https://search.worldcat.org/title/55488084
• [14]  Grønbech CH, Vording MF, Timshel PN, Sønderby CK, Pers TH, Winther O. scVAE: variational auto-encoders for single-cell gene expression data. Bioinformatics. 2020; 36: 4415-4422. doi:10.1093/bioinformatics/btaa293
• [15]  大野健太. 深層生成モデルによるテーブルデータ生成と仮想人体生成モデルへの応用. JSBi Bioinformatics Review. 2024; 5: 16-27. doi:10.11234/jsbibr.2024.2
• [16]  Margossian CC, Blei DM. Amortized variational inference: When and why? arXiv [stat.ML]. 2023. Available: http://arxiv.org/abs/2307.11018
• [17]  Higgins I, Matthey L, Pal A, Burgess C, Glorot X, Botvinick M, et al. beta-VAE: Learning Basic Visual Concepts with a Constrained Variational Framework. International Conference on Learning Representations. 2017. Available: https://openreview.net/pdf?id=Sy2fzU9gl
• [18]  Li C, Gao H, She Y, Bian H, Chen Q, Liu K, et al. Benchmarking AI Models for In Silico Gene Perturbation of Cells. Bioinformatics. bioRxiv; 2024. Available: https://www.biorxiv.org/content/10.1101/2024.12.20.629581v2.full.pdf
• [19]  Farrell JA, Wang Y, Riesenfeld SJ, Shekhar K, Regev A, Schier AF. Single-cell reconstruction of developmental trajectories during zebrafish embryogenesis. Science. 2018; 360. doi:10.1126/science.aar3131
• [20]  Ji AL, Rubin AJ, Thrane K, Jiang S, Reynolds DL, Meyers RM, et al. Multimodal analysis of composition and spatial architecture in human squamous cell carcinoma. Cell. 2020; 182: 497-514.e22. doi:10.1016/j.cell.2020.05.039
• [21]  La Manno G, Soldatov R, Zeisel A, Braun E, Hochgerner H, Petukhov V, et al. RNA velocity of single cells. Nature. 2018; 560: 494-498. doi:10.1038/s41586-018-0414-6
• [22]  Bergen V, Lange M, Peidli S, Wolf FA, Theis FJ. Generalizing RNA velocity to transient cell states through dynamical modeling. Nat Biotechnol. 2020; 38: 1408-1414. doi:10.1038/s41587-020-0591-3
• [23]  Gao M, Qiao C, Huang Y. UniTVelo: temporally unified RNA velocity reinforces single-cell trajectory inference. Nat Commun. 2022; 13: 6586. doi:10.1038/s41467-022-34188-7
• [24]  Cui H, Maan H, Vladoiu MC, Zhang J, Taylor MD, Wang B. DeepVelo: deep learning extends RNA velocity to multi-lineage systems with cell-specific kinetics. Genome Biol. 2024; 25: 27. doi:10.1186/s13059-023-03148-9
• [25]  Eng C-HL, Lawson M, Zhu Q, Dries R, Koulena N, Takei Y, et al. Transcriptome-scale super-resolved imaging in tissues by RNA seqFISH. Nature. 2019; 568: 235-239. doi:10.1038/s41586-019-1049-y
• [26]  Cho C-S, Xi J, Si Y, Park S-R, Hsu J-E, Kim M, et al. Microscopic examination of spatial transcriptome using Seq-Scope. Cell. 2021; 184: 3559-3572.e22. doi:10.1016/j.cell.2021.05.010
• [27]  Seferbekova Z, Lomakin A, Yates LR, Gerstung M. Spatial biology of cancer evolution. Nat Rev Genet. 2023; 24: 295-313. doi:10.1038/s41576-022-00553-x
• [28]  Haviv D, Remšík J, Gatie M, Snopkowski C, Takizawa M, Pereira N, et al. The covariance environment defines cellular niches for spatial inference. Nat Biotechnol. 2025; 43: 269-280. doi:10.1038/s41587-024-02193-4
• [29]  Veličković P, Cucurull G, Casanova A, Romero A, Liò P, Bengio Y. Graph Attention Networks. arXiv [stat.ML]. 2017. Available: http://arxiv.org/abs/1710.10903
• [30]  山口秀輝, 齋藤裕. タンパク質の言語モデル. JSBi Bioinformatics Review. 2023; 4: 52-67. doi:10.11234/jsbibr.2023.1
• [31]  Dalla-Torre H, Gonzalez L, Mendoza-Revilla J, Lopez Carranza N, Grzywaczewski AH, Oteri F, et al. Nucleotide Transformer: building and evaluating robust foundation models for human genomics. Nat Methods. 2025; 22: 287-297. doi:10.1038/s41592-024-02523-z
• [32]  Zhang Y, Lang M, Jiang J, Gao Z, Xu F, Litfin T, et al. Multiple sequence alignment-based RNA language model and its application to structural inference. Nucleic Acids Res. 2024; 52: e3. doi:10.1093/nar/gkad1031
• [33]  Cui H, Wang C, Maan H, Pang K, Luo F, Duan N, et al. scGPT: toward building a foundation model for single-cell multi-omics using generative AI. Nat Methods. 2024; 21: 1470-1480. doi:10.1038/s41592-024-02201-0
• [34]  Hao M, Gong J, Zeng X, Liu C, Guo Y, Cheng X, et al. Large-scale foundation model on single-cell transcriptomics. Nat Methods. 2024; 21: 1481-1491. doi:10.1038/s41592-024-02305-7
• [35]  Nathan A, Asgari S, Ishigaki K, Valencia C, Amariuta T, Luo Y, et al. Single-cell eQTL models reveal dynamic T cell state dependence of disease loci. Nature. 2022; 606: 120-128. doi:10.1038/s41586-022-04713-1
• [36]  Yazar S, Alquicira-Hernandez J, Wing K, Senabouth A, Gordon MG, Andersen S, et al. Single-cell eQTL mapping identifies cell type-specific genetic control of autoimmune disease. Science. 2022; 376: eabf3041. doi:10.1126/science.abf3041

著者略歴

小嶋 泰弘 2020年、東京大学大学院新領域創成科学研究科メディカル情報生命専攻で木立尚孝准教授の元、博士を取得。その後、名古屋大学、東京医科歯科大学で、西川博嘉教授、島村徹平教授のもとで研究・教育活動に従事した。2023年より、国立がん研究センターにて計算生命科学ユニットを主宰し、AI・数理モデリングを活用した一細胞・空間オミクス解析技術の開発を行っている。 https://www.ncc.go.jp/jp/ri/division/computational_life_science/index.html