ガウス過程のバイオインフォマティクスへの応用

## この記事について

[バイオインフォマティクス Advent Calendar 2019](https://qiita.com/advent-calendar/2019/bioinfo) の25日目の記事です。

[『ガウス過程と機械学習』（持橋・大羽 2019）](https://www.kspub.co.jp/book/detail/1529267.html)を読み、ガウス過程がバイオインフォマティクスへどう応用されているか気になったので、調べた内容をまとめました。

ガウス過程は確率的生成モデルの一種です。ガウス過程の定義、導出、式変形、推定アルゴリズムなどについては[『ガウス過程と機械学習』（持橋・大羽 2019）](https://www.kspub.co.jp/book/detail/1529267.html)を読んでいただければと思います。ガウス過程一般についてはネット上に良質な資料があるので（例：[持橋先生の講義資料](https://www.ism.ac.jp/~daichi/lectures/H26-GaussianProcess/gp-lecture2-daichi.pdf)）、この記事ではガウス過程の説明は控えたいと思います🙇‍♀️






<!-- 今回話したい点は、ガウス過程回帰とGPLVM（Gaussian process (GP) latent variable models）です。 -->

<!-- 以下ではどんなモデル化をしているかについてのみ説明します。 -->


<!-- ### ガウス過程回帰

ガウス過程回帰の場合、データを入力 $x_i={xid}$と出力の$y_i$の組、カーネル行列（グラム行列）を$K$とすると、ハイパーパラメータ $\theta$ です。


$y=N(0|K)$
$K = {k(x_i,x_j)}$

この場合、カーネルのハイパーパラメータ $\theta$を推定します。

### GPLVM
また、GPLVMの場合は、データ$y_i$が与えられた時に、潜在変数$x_i={xid}$があり、カーネル行列（グラム行列）を$K$を通じて


$y=N(0|K)$
$K = {k(x_i,x_j)}$

となります。この場合、潜在変数 $x$ とカーネルのハイパーパラメータ $\theta$の両方を推定します。
 -->


## 転写制御解析での応用

遺伝子発現量の時系列データが与えられた時に、既知制御因子（転写因子など）の「活性の時系列」を推定する、という問題設定に対し、 [Neil D. Lawrence](https://scholar.google.com/citations?user=r3SJcvoAAAAJ) や [Antti Honkela](https://scholar.google.com/citations?user=XsyLs6AAAAAJ) がガウス過程を積極的に取り入れてきました。

[Lawrence+, Adv. Neural Inf. Process. Syst., 2007](https://papers.nips.cc/paper/3119-modelling-transcriptional-regulation-using-gaussian-processes) では、TFの活性（濃度）の時系列をガウス過程としてモデル化し、観測モデルが線型および非線形の微分方程式の時の標的遺伝子の発現量の観測モデルを定式化しました。これにより、遺伝子発現量の時系列データが与えられた時に、既知標的遺伝子の感受性を推定を計算効率よくできたと主張しています。続く [Gao+, Bioinformatics, 2008](https://doi.org/10.1093/bioinformatics/btn278) では、ミカエリスメンテン式や転写抑制の場合での定式化を提案しました。

また、[Honkela+, PNAS, 2010](https://www.ncbi.nlm.nih.gov/pmc/articles/PMC2867914/) では、遺伝子発現データのみから転写因子のターゲット遺伝子を推定する（Prioritizationが目的）を提案しています。さらに、[Titsias+, BMC Syst. Biol., 2012](https://bmcsystbiol.biomedcentral.com/articles/10.1186/1752-0509-6-53)では、問題設定を複数の転写因子があった場合に拡張させ、転写因子のターゲット遺伝子を推定する手法を提案しています。


### [Gao+, Bioinformatics, 2008] の例

例えば、[Gao+, Bioinformatics, 2008](https://doi.org/10.1093/bioinformatics/btn278) では、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)$$

<!-- 入力は
- サンプリング時刻 $t$
- 時刻 $t$ での遺伝子$j$のmRNA量 $x_j(t)$ -->

$B_j$, $S_j$, $D_j$ を推定し、
パラメータを推定した後に、$f(t)$ の予測分布

さらに、サンプリングしていない時刻も含む全時刻での遺伝子$j$のmRNA量 $x_j(t)$　を得られます。

著者らは、ガウス過程によるモデル化のメリットについて、（時間の関数の離散化をしないことで）パラメータ数を抑えられる、非線形化（ミカエリスメンテン式、抑制モデル）にも開かれているといった主張をしています。


<!--  Lawrence, N. D., Sanguinetti, G. & Rattray, M. Modelling transcriptional regulation using Gaussian Processes. Adv. Neural Inf. Process. Syst. 19 785–792 (2007). -->
<!-- 
要は:  TFの濃度と既知ターゲット遺伝子のTFに対する感度を発現データから推定
質的差分:  計算効率がよいのでスケールできる、アドホックな離散化が必要ない、内挿がいらない
技術の肝: Linear / non-linear response model を定式化。TFの濃度をGPとし、感度とmRNA分解速度とベースライン発現レベルをハイパーパラメータとした
有効性の検証: p53のデータに適用、先行研究と推定結果が近いことを確認。
議論: p53のデータに適用、先行研究と推定結果が近いことを確認。
次に読むべき論文:  Barenco+2006, Rogers+2006 -->

<!-- Gao, P., Honkela, A., Rattray, M. & Lawrence, N. D. Gaussian process modelling of latent chemical species: Applications to inferring transcription factor activities. Bioinformatics 24, i70–i75 (2008).
Gaussian process modelling of latent chemical species: applications to inferring transcription factor activities
[Gao+, Bioinformatics, 2008]

遺伝子発現量の時系列データから上流転写因子の活性を予測する。
制御ネットワークが既知であることを想定している。 -->

<!-- 要は:  ターゲット遺伝子の発現データからTFの活性を推定する問題を、GPで縛ったLatent chemical species をベイズ周辺化する方法を提案
質的差分: なるべく仮定を少なくした。パラメータ数爆発しない。非線形化（ミカエリスメンテン式、抑制）
技術の肝: GPによるモデル化（時間の関数の離散化をしないことでパラメータの数を抑えられる)、非線形化にも開かれている
有効性の検証: 推定結果がウェスタンブロットや先行研究と目視で似ていること
議論: ターゲット遺伝子は与えられているが、そのTFだけに制御されているかはわからない
次に読むべき論文: Lawrence,N.D. et al. (2007), Barenco+2006, Alon2006 -->




<!-- 
 Honkela, A. et al. Model-based method for transcription factor target identification with limited data. Proc. Natl. Acad. Sci. 107, 7793–7798 (2010).
要は:  遺伝子発現データのみから転写因子のターゲット遺伝子を推定する（Prioritizationが目的）
質的差分: 発現データのみを使用
技術の肝: 転写因子タンパク質の発現プロファイルはGPで推定して積分消去する 
有効性の検証: ChIP-chipやTF KOの発現データでターゲット遺伝子を定義し、提案手法のランクが高い遺伝子と重なり。TFBSからの距離の分布。
議論: 恒常発現の遺伝子は低発現を除いたら除けた。ターゲット遺伝子が多様な発現プロファイルを示す時、単にTFの発現との相関をとるよりも精度が高かった 
次に読むべき論文: 8, 11,14, 15, 18 -->


<!-- TFタンパク質の時系列を推定
TFから各遺伝子への影響のパラメータを推定
入力は時系列発現データ、働いているTF(複数)、TFの既知ターゲット。出力はTFのターゲット遺伝子のリスト。 -->

<!-- ### Titsias, M. K., Honkela, A., Lawrence, N. D. & Rattray, M. Identifying targets of multiple co-regulating transcription factors from expression time-series by Bayesian model comparison. BMC Syst. Biol. 6, 53 (2012).

要は:  入力は時系列発現データ、働いているTF(複数)、TFの既知ターゲット。出力はTFのターゲット遺伝子のリスト。
質的差分: 複数TFのターゲット遺伝子を効率的に探せる（教師ありの問題設定をGPでやる方法はあったが、単一のtfのみだった。複数tfでシグモイド関数の入ったODEを解くものがあった。それらを組み合わせた）
技術の肝: ODEとGPと非線形なシグモイド関数による転写制御モデル。既知ターゲットの発現データから、TFタンパク質の時系列を推定。次にそれを用いて、TFによる各遺伝子への影響のパラメータを推定する。 
有効性の検証: 人工データ。実データでChIPとの重なりが大きいこと
議論: 
次に読むべき論文: 21, 18, 15-17, 5 -->


## ChIP-seq データ解析への応用
Neil D. Lawrence や Antti Honkela は Pol II ChIP-seqの解析に対してもガウス過程を導入しています。

<!-- 
 wa Maina, C. et al. Inference of RNA Polymerase II Transcription Dynamics from Chromatin Immunoprecipitation Time Course Data. PLoS Comput. Biol. 10, e1003598 (2014). -->

[wa Maina+, PLoS Comput. Biol., 2014](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1003598) では、Pol II ChIP-seq時系列データのgene bodyに対するシグナルに対し、時間方向と空間（ゲノム上のbin）報告での共分散を同時に考慮した convolved Gaussian process を適用しました。これにより、転写開始・伸長活性を反映するPol II のシグナルの動態にから転写速度や活性を推定し、遺伝子を転写の動態に応じて分類することができます。時間方向だけでなく、空間方向も同時に考慮するのが面白いです。



<!-- 要は:  Pol II ChIP-Seqの時系列データからPol IIの転写動態を推定した
質的差分: 転写速度の変化を推定できる（一定性を仮定しない）
技術の肝: convolved GPで時間と空間の関係をモデル化。遅れも推定。
有効性の検証: 人工データで転写速度の遅れの推定が正解に近い。実データはホルモン刺激後の時系列Pol II ChIP-Seqから、transcriptionall engaged geneを抽出、転写速度と転写活性を推定した。転写活性による遺伝子クラスタリング。
議論: モデルへのフィッティングでPoll II が動かない遺伝子もわかる
次に読むべき論文: -->



<!-- ### Honkela, A. et al. Genome-wide modeling of transcription kinetics reveals patterns of RNA production delays. Proc. Natl. Acad. Sci. 112, 201420404 (2015). -->
[Honkela+, PNAS, 2015](https://www.pnas.org/content/112/42/13115)では、RNAの転写活性とmRNA量の遅れ（RNA production delay）の時間変化をガウス過程でモデル化するアプローチが提案されています。Pol II ChIP-seqとRNA-seqの時系列データからそれぞれ転写活性とmRNA発現量の時間変化を定量して入力データとし、RNA production delayを推定することで、

<!-- 
要は:  transcription kinetics が重要であるかという問いをゲノムワイドに検討した
質的差分: RNA processingに由来する遅れの存在をゲノムワイドに評価した
技術の肝: 遅れ（RNA production delay）をパラメータとしていれた。ベイズでパラメータ推定。GPで時間発展を記述。
有効性の検証: シミュレーションデータでパラメータの推定。
議論: RNA production delayとintron length, intron retentionなどの関係を考察。
次に読むべき論文: 1,5,6, 9,10,11, 4,7,8 -->


## シングルセルRNA-seq解析への応用

### 擬時間推定への応用
[PLOS Comput Biol. 2016](https://doi.org/10.1371%2Fjournal.pcbi.1005212)、[Reid+, Bioinformatics, 2016](https://doi.org/10.1093%2Fbioinformatics%2Fbtw372)、[Ahmed+, Bioinformatics, 2018](https://academic.oup.com/bioinformatics/article/35/1/47/5047752)では、GPLVMを用いた細胞の擬時間推定手法が提案されています。


また、[Boukouvalas+, Genome Biology, 2018](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-018-1440-2)では、BGP (branching Gaussian process)という方法を提案しています。これは、時系列の混合モデルを仮定し、あらかじめ細胞の擬時間と大まかな分岐パターンが与えられた時に、遺伝子ごとの分岐パターンを推定します。


### 遺伝子発現の細胞間分散の解析への応用

[Buettner+, Nature Biotechnology, 2015](https://www.nature.com/articles/nbt.3102)では、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](https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1334-8) で提案されている f-scLVM では、（定式化は似ているものの）細胞の潜在変数を推定するのにGPLVMを用いないアプローチになっています。ちなみに、scLVMおよびf-scLVMの論文の最終著者の [Oliver Stegle](https://scholar.google.com/citations?user=ClSXZ4IAAAAJ) は後述する [Fusi+, PLoS Comput. Biol., 2012](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002330) で Neil D. Lawrence と共著者です。



<!-- 次に、$\Sigma_h = \mathbf{X}\mathbf{X}^{\mathrm{T}}$ が求まったら、今度は $\Sigma_h$ を用いて遺伝子発現量の細胞間変動を分解します。ここで、 $\mathbf{\psi}_{e}$ をredisual noise、 $\mathbf{\psi}_{t}$ を technical noise とすると、発現量は以下のようにモデル化されます。

$$\mathbf{y_g} = \mu \mathbf{I} + \sum_{h=1}^{H} \mathbf{u}_h + \mathbf{\psi}_{e} + \mathbf{\psi}_{t}$$

ここで $p(\mathbf{u}_h)=N(\mathbf{u}_h|\mathbf{0}, \sigma^2_{g,h}\mathbf{\Sigma}_h)$ 、 $\mathbf{\psi}_{e} \sim N(\mathbf{0}, \delta^2_g \mathbf{I})$ 、 $\mathbf{\psi}_{t} \sim N(\mathbf{0}, \nu^2_g \mathbf{I})$ とし、$\mathbf{u}_h$ 、 $\mathbf{\psi}_{e}$ 、 $\mathbf{\psi}_{t}$ について周辺化すると、以下の式が得られます。


$$p(\mathbf{y_g} | \mu_g, \{\sigma^2_{g,h}, \Sigma_{h}\}_{h=1}^{H}, \nu_g^{2}, \delta_{g}^2] )  = N(\mathbf{y}_g| \mu_g \mathbf{1}, \sum_{h=1}^{H} \sigma^2_{g,h}\Sigma_h + \nu^2_{g} \mathbf{I} + \delta^2_{g} \mathbf{I})$$

上の式で、$\{\sigma^2_{g,h}\}$ 、 $\nu_{g}^2$ 、 $\delta_{g}^2$ を推定することで、遺伝子 $g$ の細胞間変動に対してどの潜在因子の影響が大きいかといった解釈ができます。 -->



<!-- 


scLVM estimates the cell-to-cell similarity matrix due to hidden factors such as cell cycle by fitting a low-rank covariance matrix to a set of predefined marker genes. The fitting process allows known covariates (if available) to be accounted for. In some instances, it may also be desirable to account for multiplicative interactions between the identified factors (see Section 3.1.2).

予め細胞周期に関連する遺伝子を入れておく

GPLVM だけど

To estimate the components of variance, scLVM employs a linear mixed model that is fitted to the expression levels of each gene, decomposing sources of variation.

線型混合モデルによる、変動の原因の分解。

$$Y_h = $$ -->


<!-- ### Computational analysis of cell-to-cell heterogeneity in single-cell RNA-sequencing data reveals hidden subpopulations of cells.  -->




<!-- BGP: identifying gene-specific branching dynamics from single-cell data with a branching Gaussian process -->





<!-- 


A robust nonlinear low-dimensional manifold for single cell RNA-seq data

https://www.biorxiv.org/content/10.1101/443044v1

GPLVM -->



<!-- #### GrandPrix: scaling up the Bayesian GPLVM for single-cell data 
https://academic.oup.com/bioinformatics/article/35/1/47/5047752
[Ahmed, Bioinformatics, 2018]
GPLVM
細胞サンプリング時刻情報つきのscRNA-seqに対する擬時間推定で
高速化した。
擬時間上の時点間でのカーネル -->




<!-- 
GPLVMとは関係ない #### f-scLVM: scalable and versatile factor analysis for single-cell RNA-seq
[Buettner+, Genome Biology, 2017]
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-017-1334-8 -->



<!-- データの変動を、潜在変数・共変量によって説明するという発想自体は、GPLVM以外でもあります。例えば、[Surrogate variable analysis](https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.0030161)とか。 -->

<!-- scLVM is closely related to previous approaches that correct for hidden confounding factors in gene expression data 21,48 and the inference employed to fit hidden factors builds on the PANAMA model30.

Related approaches have been proposed to account for relatedness between individuals in the context of expression Quantitative Trait Loci (eQTL) studies30, where an individual-to-individual covariance is inferred to explain the heterogeneity in gene expression levels between individuals rather than cells. -->

<!-- 21  Leek, J.T. & Storey, J.D. Capturing heterogeneity in gene expression studies by surrogate variable analysis. PLoS Genet. 3, e161 (2007). https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.0030161 -->

<!-- 48 Gagnon-Bartsch, J.A. & Speed, T.P. Using control genes to correct for unwanted variation in microarray data. Biostatistics 13, 539–552 (2012).
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3577104/ -->

<!-- Lippert, C. et al. FaST linear mixed models for genome-wide association studies. Nat. Methods 8, 833–835 (2011). -->




<!-- ### Bayesian Gaussian Process Latent Variable Models for pseudotime inference in single-cell RNA-seq data

Kieran Campbell, Christopher Yau
https://www.biorxiv.org/content/10.1101/026872v1
[Campbell, bioRxiv, 2015] -->


## 空間トランスクリプトームデータでの応用

最近、空間トランスクリプトームというジャンルの技術が出てきました。これは遺伝子発現量データに二次元や三次元の空間座標が紐づいているものです。実験技術的には、イメージングベースのもの（MERFISH、SeqFISH+など）やシーケンシングベース（Visiumなど）があります。

空間トランスクリプトームデータへの応用では、マルチカーネル学習のような枠組みで、各カーネルの寄与を解釈する方法が提案されています。

例えば、[Svensson+, Nature Methods, 2018](https://dx.doi.org/10.1038/nmeth.4636)では、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](https://dx.doi.org/10.1016/j.celrep.2019.08.077) では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}} $$



<!-- SpatialDE: identification of spatially variable genes -->


<!-- ### Modeling Cell-Cell Interactions from Spatial Molecular Data with Spatial Variance Component Analysis

SVCA () -->





<!-- SpatialDE: identification of spatially variable genes
[Svensson+, Nature Methods, 2018]

https://dx.doi.org/10.1038/nmeth.4636 -->


## 遺伝統計学への応用

[Fusi+, PLoS Comput. Biol., 2012](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002330) では、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を同時に推定することで、コホート間での再現性があがったと主張しています。

論文では明示的にはガウス過程とは書いてませんが、マテメソを追うと（１）線型回帰モデルをガウス分布として設定し、（２）重みパラメータの事前分布をガウス分布と仮定し、（３）重みパラメータを周辺化して多変量ガウス分布を得ており、典型的なガウス過程回帰の導出になっています。


<!-- 
30. Fusi, N., Stegle, O. & Lawrence, N.D. Joint modelling of confounding factors and prominent genetic regulators provides increased accuracy in genetical genomics studies. PLoS Comput. Biol. 8, e1002330 (2012).
https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002330

PANAMA model: 発現量を分解する線型モデル -->


## DNase-seq解析への応用

DNase I-seq は、DNase IをゲノムDNAに供した後にDNAシーケンシングを行うことで、ゲノムの各領域のクロマチンアクセシビリティを計測する実験技術です。

ゲノム上である転写因子のモチーフ出現位置に関してリード（やその端点）の分布を積み上げると、モチーフ内の位置に応じたリード分布を示すことが知られています。このような現象はDNase I フットプリントと呼ばれ、ここからこのようなモチーフ出現位置に特異的なリード分布を示す場所には転写因子結合がしているのではないか、とDNase-seqからの転写因子結合部位の推定に使われています。

[Sherwood+, Nat. Biotechnol., 2014](https://doi.org/10.1038/nbt.2798)では、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データによる結果をいわゆるパイオニア因子を明らかにするという応用にも利用しています。

<!-- 

モチーフ出現位置 $m$ に転写因子が結合しているかどうかを表す潜在変数 $I_m$ （結合していれば１、していなければ0）を設定します。また、 $y_m$ を$m$の中心点の座標、 $W$ をモチーフ出現位置周辺で考慮する幅

$I_m=1$ かつ $|y_m-j| \leq W$ の場合は、

$$\hat{\mu}_i^{+} = \mu_i^{+} + \beta^{+}_{i-j,l}  : |y_m -j| \leq W \mathrm{ and } I_m =1$$

$y_m=0$ の場合は、

$$\hat{\mu}_{i}^{+} = \mu_i^{+}$$

また、 $I_j$ の事前分布は、モチーフ出現位置 $j$ のモチーフスコア $s_j$ とリードカウント $c_j$ で定義されます。

$$\log( P (I_j=1) ) = f(s_j) + g(c_j)$$

これによって、 -->




<!-- PIQ 
Discovery of directional and nondirectional pioneer transcription factors by modeling DNase profile magnitude and shape
[Sherwood+, Nat. Biotechnol., 2014]  -->

<!-- 10.1038/nbt.2798 -->



## 配列設計への応用
<!-- ### mGPfusion: predicting protein stability changes with Gaussian process kernel learning and data fusion -->
[Jokinen+, Bioinformatics, 2018](https://academic.oup.com/bioinformatics/article/34/13/i274/5045756) では、少数の実験データを多数の分子シミュレーションの結果で補五つ変異によるタンパク質の安定性変化を予測するパイプライン mGPfusion を提案しいます。その中でガウス過程を用いて予測モデルを構築しています。

<!-- 最近の強化学習やGANでの材料設計に比べてどうなのかはわからない。 -->

<!-- # TODO -->
<!-- 
J. Hensman, P. Papastamoulis, P. Glaus, A. Honkela, and M. Rattray. Fast and accurate approximate inference of transcript expression from RNA-seq data. Bioinformatics (2015). -->


<!-- P. Glaus, A. Honkela*, and M. Rattray*. Identifying differentially expressed transcripts from RNA-seq data with biological variation. Bioinformatics 28(13):1721-1728 (2012). -->

## まとめ

ガウス過程のバイオインフォマティクスへの応用について、概要だけ紹介しました。

今の所、時間、ゲノム座標、空間、擬時間、トランスクリプトームの類似度を表すカーネルを用いることで、ガウス過程回帰で変動解析を行ったり、GPLVMで未知の（交絡）因子を推定するといったアプローチが多い気がしました。

今後、シングルセルマルチオミクスデータや、時間・空間的に高密度にサンプリングされた環境DNAデータの解析といった、他の分野にも使われるかもしれないと思いました。また、[遺伝子型の効果がエピジェネティックな効果でよく見えない時](https://nosumi.exblog.jp/20278470/) に、DNAメチル化のサンプル間類似度をカーネルで表して遺伝子型の効果を分離できないかなと思いました。

ただ、必ずしも「ガウス過程じゃないといけないのか」という点に答えられていない論文もあり、どんな問題設定ならガウス過程が合っているのかはもっと深く考えないといけないのではないかと思いました（私自身もこの点はまだ掴み切れていないところがあるので、もっと勉強したいと思います）。