この記事について
バイオインフォマティクス Advent Calendar 2019 の25日目の記事です。
『ガウス過程と機械学習』（持橋・大羽 2019）を読み、ガウス過程がバイオインフォマティクスへどう応用されているか気になったので、調べた内容をまとめました。
ガウス過程は確率的生成モデルの一種です。ガウス過程の定義、導出、式変形、推定アルゴリズムなどについては『ガウス過程と機械学習』（持橋・大羽 2019）を読んでいただければと思います。ガウス過程一般についてはネット上に良質な資料があるので（例：持橋先生の講義資料）、この記事ではガウス過程の説明は控えたいと思います🙇♀️
転写制御解析での応用
遺伝子発現量の時系列データが与えられた時に、既知制御因子（転写因子など）の「活性の時系列」を推定する、という問題設定に対し、 Neil D. Lawrence や Antti Honkela がガウス過程を積極的に取り入れてきました。
Lawrence+, Adv. Neural Inf. Process. Syst., 2007 では、TFの活性（濃度）の時系列をガウス過程としてモデル化し、観測モデルが線型および非線形の微分方程式の時の標的遺伝子の発現量の観測モデルを定式化しました。これにより、遺伝子発現量の時系列データが与えられた時に、既知標的遺伝子の感受性を推定を計算効率よくできたと主張しています。続く Gao+, Bioinformatics, 2008 では、ミカエリスメンテン式や転写抑制の場合での定式化を提案しました。
また、Honkela+, PNAS, 2010 では、遺伝子発現データのみから転写因子のターゲット遺伝子を推定する（Prioritizationが目的）を提案しています。さらに、Titsias+, BMC Syst. Biol., 2012では、問題設定を複数の転写因子があった場合に拡張させ、転写因子のターゲット遺伝子を推定する手法を提案しています。
[Gao+, Bioinformatics, 2008] の例
例えば、Gao+, Bioinformatics, 2008 では、mRNA量の時系列データと転写制御関係が与えられた時に上流制御因子の活性の時系列を推定する問題にガウス過程を応用しています。
時刻 $t$ での遺伝子$j$のmRNA量 $x_j(t)$ が上流転写因子によって活性化されているとします。この時、
- 遺伝子$j$のベースラインの発現量 $B_j$
- 遺伝子$j$に対する上流転写因子の影響の大きさ $S_j$
- 遺伝子$j$のmRNAの分解率 $D_j$
- 上流の転写因子の時刻 $t$ での活性 $f(t)$
とし、転写活性化を線型微分方程式として表すと、$x_j(t)$ の観測モデルは以下のようになります。
$$x_j(t) = \frac{B_j}{D_j} + S_j \int_{0}^{t} e^{-D_j(t-u)} f(u) du$$
ここで、 $f(t)$ をガウス過程でモデリングします。すなわち、
$$f(t) \sim N(0, K)$$
$B_j$, $S_j$, $D_j$ を推定し、
パラメータを推定した後に、$f(t)$ の予測分布
さらに、サンプリングしていない時刻も含む全時刻での遺伝子$j$のmRNA量 $x_j(t)$ を得られます。
著者らは、ガウス過程によるモデル化のメリットについて、（時間の関数の離散化をしないことで）パラメータ数を抑えられる、非線形化（ミカエリスメンテン式、抑制モデル）にも開かれているといった主張をしています。
ChIP-seq データ解析への応用
Neil D. Lawrence や Antti Honkela は Pol II ChIP-seqの解析に対してもガウス過程を導入しています。
wa Maina+, PLoS Comput. Biol., 2014 では、Pol II ChIP-seq時系列データのgene bodyに対するシグナルに対し、時間方向と空間（ゲノム上のbin）報告での共分散を同時に考慮した convolved Gaussian process を適用しました。これにより、転写開始・伸長活性を反映するPol II のシグナルの動態にから転写速度や活性を推定し、遺伝子を転写の動態に応じて分類することができます。時間方向だけでなく、空間方向も同時に考慮するのが面白いです。
Honkela+, PNAS, 2015では、RNAの転写活性とmRNA量の遅れ（RNA production delay）の時間変化をガウス過程でモデル化するアプローチが提案されています。Pol II ChIP-seqとRNA-seqの時系列データからそれぞれ転写活性とmRNA発現量の時間変化を定量して入力データとし、RNA production delayを推定することで、
シングルセルRNA-seq解析への応用
擬時間推定への応用
PLOS Comput Biol. 2016、Reid+, Bioinformatics, 2016、Ahmed+, Bioinformatics, 2018では、GPLVMを用いた細胞の擬時間推定手法が提案されています。
また、Boukouvalas+, Genome Biology, 2018では、BGP (branching Gaussian process)という方法を提案しています。これは、時系列の混合モデルを仮定し、あらかじめ細胞の擬時間と大まかな分岐パターンが与えられた時に、遺伝子ごとの分岐パターンを推定します。
遺伝子発現の細胞間分散の解析への応用
Buettner+, Nature Biotechnology, 2015では、GPLVMをscRNA-seqデータ解析に応用したscLVMという手法を提案しています。scLVMでは、scRNA-seqデータにおいて各遺伝子の細胞間での変動が細胞ごとに異なる潜在因子の値の違いによって説明できると仮定し、全遺伝子についてガウス過程回帰を行って潜在因子ごとの寄与度を調べたり、特定の潜在因子（細胞周期など）の影響を除去します。
まず、GPLVMを用いて細胞ごとの潜在変数の値（およびそれらの細胞間共分散行列）を推定します。あらかじめ設定したマーカー遺伝子群に対して、低ランクの共分散行列をフィッティングします。具体的には、細胞数を $N$ 、ある潜在因子（細胞周期など）を $h$ 、$h$に関連するマーカー遺伝子群を $g_h$ 、 $\mathbf{Y_h}$ を各細胞での発現量を表す $N \times |g_h|$の行列、$\mathbf{C}$を細胞ごとの$Q$個の既知の共変量を表す$N \times Q$の行列、$\mathbf{X}$を細胞ごとの$K$個の未知の共変量を表す$N \times K$の行列とし、 $\mathbf{Y_h}$ を以下の線型生成モデルで表します。
$$\mathbf{Y_h} = \mu + \mathbf{C}\mathbf{U} + \mathbf{X}\mathbf{W} + \mathbf{\psi}$$
ここで、既知および未知の共変量の重み $\mathbf{U}$ 、 $\mathbf{W}$ の事前分布を独立正規分布と仮定します。
$$p(\mathbf{C}) = \prod_{q=1}^{Q} N(\mathbf{u_q}|0, \sigma_u^2 \mathbf{I})$$
$$p(\mathbf{W}) = \prod_{k=1}^{Q} N(\mathbf{w_k}|0, \mathbf{I})$$
ここで、 $\mathbf{C}$ と $\mathbf{W}$ について $\mathbf{Y_h}$ の尤度を周辺化することで、以下のガウス過程の表現が得られます。
$$p(\mathbf{Y_h}|\sigma^2_u, \nu^2, \mathbf{C},\mathbf{X}) = \prod_{g=1}^{|g_h|} N(\mathbf{y_g}|\mu_g\mathbf{1}, \sigma^2_u\mathbf{C}\mathbf{C}^{\mathrm{T}} + \mathbf{X}\mathbf{X}^{\mathrm{T}} + \nu \mathbf{I} )$$
データ $\mathbf{Y_h}$ および $\mathbf{C}$ は与えられているので、 $\sigma^2_u$ 、 $\nu^2$ 、 $\mathbf{X}$ を推定します。
ここで推定された、$\Sigma_h = \mathbf{X}\mathbf{X}^{\mathrm{T}}$ を用いて、scLVMではさらに、各遺伝子の細胞間変動に対する潜在因子ごとの寄与度を調べたり、特定の潜在因子（細胞周期など）の影響を発現量から除去するといったことを行います。
まとめると、潜在因子と関連する特定の遺伝子群（細胞周期関連遺伝子など）についてGPLVMで細胞ごとの潜在変数（$\mathbf{X}$）を得て、その潜在変数の共分散行列（$\mathbf{\Sigma}_h$）に基づいて全遺伝子についてガウス過程回帰を行って潜在因子ごとの寄与度を調べたり、特定の潜在因子（細胞周期など）の影響を除去します。
なお、同じ著者らによって Buettner+, Genome Biology, 2017 で提案されている f-scLVM では、（定式化は似ているものの）細胞の潜在変数を推定するのにGPLVMを用いないアプローチになっています。ちなみに、scLVMおよびf-scLVMの論文の最終著者の Oliver Stegle は後述する Fusi+, PLoS Comput. Biol., 2012 で Neil D. Lawrence と共著者です。
空間トランスクリプトームデータでの応用
最近、空間トランスクリプトームというジャンルの技術が出てきました。これは遺伝子発現量データに二次元や三次元の空間座標が紐づいているものです。実験技術的には、イメージングベースのもの（MERFISH、SeqFISH+など）やシーケンシングベース（Visiumなど）があります。
空間トランスクリプトームデータへの応用では、マルチカーネル学習のような枠組みで、各カーネルの寄与を解釈する方法が提案されています。
例えば、Svensson+, Nature Methods, 2018では、SpatialDEという手法が提案されています。
SpatialDEでは、各遺伝子について、空間中の各座標における発現量をガウス過程でモデル化しています。
この際、共分散行列は、座標のペア同士が近いほど大きくなるカーネルから成る $\sigma_s^2 \cdot\mathbf{\Sigma}$ と 座標間で独立な共分散行列である $\sigma_s^2 \cdot\delta \cdot \mathbf{I}$の和になっています。
$$P(y|\mu,\sigma_s^2,\delta,\mathbf{\Sigma})=N(y|\mu \cdot \mathbf{1}, \sigma_s^2 \cdot (\mathbf{\Sigma} + \delta \cdot \mathbf{I}))$$
これにより、SpatialDEでは「距離が近いと発現量が近くなる」パターンを持ちつつ変動している遺伝子を探索することができます。
また、SpatialDEと同じグループによる Arnol+, Cell Reports, 2018 ではSVCAという方法が提案されています。SVCAでは、ガウス過程の共分散行列を、内在的要因（細胞型など）・細胞間相互作用・環境要因（距離が近いと発現量が近くなる）のカーネルの重み付き和としています。これらの重みを推定することで、各遺伝子がどのような要因によって変動しているかを推定します。
$$P(Y) = N(0, K_{intrinsic} + K_{cell-cell interation} + K_{environment} + \delta^2_{\epsilon} I_n)$$
$$ K_intrinsic = \sigma_I^2 XX^{T} $$
$$ K_{cell-cell interation} = \sigma_{C-C}^2 ZXX^{T}Z^{T} \mathrm{where} Z_{i,j} = exp(\frac{-d_{i,j}^2}{2l^2}) $$
$$ K_{environment} = \sigma_{E} exp{\frac{-d_{i,j}^2}{2l^2}} $$
遺伝統計学への応用
Fusi+, PLoS Comput. Biol., 2012 では、eQTL（expression QTL）解析にガウス過程回帰を応用したPANAMA modelを提案しています。各サンプルにおける遺伝子発現量が、$K$個のSNP（$N \times K$ の行列 $S$ で表現）とQ個の潜在因子（$N \times Q$ の行列 $X$で表現）で決まるというモデルです（下記）。
$$Y=\mathbf{\mu} + \mathbf{S}\mathbf{V} + \mathbf{X}\mathbf{W} + \mathbf{\epsilon}$$
著者らは、交絡因子の影響を除去する既存手法ではtrans eQTLを交絡因子と混同する危険がある一方、提案手法では交絡因子とtrans eQTLを同時に推定することで、コホート間での再現性があがったと主張しています。
論文では明示的にはガウス過程とは書いてませんが、マテメソを追うと（１）線型回帰モデルをガウス分布として設定し、（２）重みパラメータの事前分布をガウス分布と仮定し、（３）重みパラメータを周辺化して多変量ガウス分布を得ており、典型的なガウス過程回帰の導出になっています。
DNase-seq解析への応用
DNase I-seq は、DNase IをゲノムDNAに供した後にDNAシーケンシングを行うことで、ゲノムの各領域のクロマチンアクセシビリティを計測する実験技術です。
ゲノム上である転写因子のモチーフ出現位置に関してリード（やその端点）の分布を積み上げると、モチーフ内の位置に応じたリード分布を示すことが知られています。このような現象はDNase I フットプリントと呼ばれ、ここからこのようなモチーフ出現位置に特異的なリード分布を示す場所には転写因子結合がしているのではないか、とDNase-seqからの転写因子結合部位の推定に使われています。
Sherwood+, Nat. Biotechnol., 2014では、DNase-seqデータによる転写因子結合の推定にガウス過程を用いるPIQという手法提案しています。PIQデハ、ゲノム上のモチーフ出現位置からの相対的な位置 $i$ におけるリードの生成モデルを以下のように定式化します。
まず、観測モデルとして、位置 $i$ におけるリード数 $x_i$ がポワソン分布に従うと仮定します:
$$x_i \sim Poisson(\exp(\mu_i))$$
この時、${\mu_i}$ はガウス過程に従うと仮定します:
$$\mu_i \sim N(\mu_0, \Sigma)$$
$$\Sigma_{(i,j)} = \sigma_0 \cdot k_{|i-j|}$$
さらに、モチーフ出現位置 $m$ に転写因子が結合・非結合の２状態を表す潜在変数 $I_m$ （結合していれば１、していなければ0）を導入し、 $I_m$ の値によって $\mu_i$ が異なる混合モデルになっています。そして、$I_m$ の事後確率を推定することで、 $m$ に転写因子が結合しているかを決定します。
PIQでは、転写因子結合予測だけでなく、時系列DNase-seqデータによる結果をいわゆるパイオニア因子を明らかにするという応用にも利用しています。
配列設計への応用
Jokinen+, Bioinformatics, 2018 では、少数の実験データを多数の分子シミュレーションの結果で補五つ変異によるタンパク質の安定性変化を予測するパイプライン mGPfusion を提案しいます。その中でガウス過程を用いて予測モデルを構築しています。
まとめ
ガウス過程のバイオインフォマティクスへの応用について、概要だけ紹介しました。
今の所、時間、ゲノム座標、空間、擬時間、トランスクリプトームの類似度を表すカーネルを用いることで、ガウス過程回帰で変動解析を行ったり、GPLVMで未知の（交絡）因子を推定するといったアプローチが多い気がしました。
今後、シングルセルマルチオミクスデータや、時間・空間的に高密度にサンプリングされた環境DNAデータの解析といった、他の分野にも使われるかもしれないと思いました。また、遺伝子型の効果がエピジェネティックな効果でよく見えない時 に、DNAメチル化のサンプル間類似度をカーネルで表して遺伝子型の効果を分離できないかなと思いました。
ただ、必ずしも「ガウス過程じゃないといけないのか」という点に答えられていない論文もあり、どんな問題設定ならガウス過程が合っているのかはもっと深く考えないといけないのではないかと思いました（私自身もこの点はまだ掴み切れていないところがあるので、もっと勉強したいと思います）。