Semi-Supervised Non-Parametric Bayesian Modelling of Spatial Proteomics

TL;DR

This paper introduces a semi-supervised non-parametric Bayesian model using Gaussian process mixtures to analyze spatial proteomics data, capturing complex correlation structures within sub-cellular niches efficiently.

Contribution

It develops a novel Gaussian process mixture model with an efficient Hamiltonian-within-Gibbs sampler and covariance matrix decomposition for spatial proteomics analysis.

Findings

Effective modeling of protein localization with Gaussian processes.

Reduced computational complexity via tensor decomposition and specialized algorithms.

Open-source R package implementation available.

Abstract

Understanding sub-cellular protein localisation is an essential component to analyse context specific protein function. Recent advances in quantitative mass-spectrometry (MS) have led to high resolution mapping of thousands of proteins to sub-cellular locations within the cell. Novel modelling considerations to capture the complex nature of these data are thus necessary. We approach analysis of spatial proteomics data in a non-parametric Bayesian framework, using mixtures of Gaussian process regression models. The Gaussian process regression model accounts for correlation structure within a sub-cellular niche, with each mixture component capturing the distinct correlation structure observed within each niche. Proteins with a priori labelled locations motivate using semi-supervised learning to inform the Gaussian process hyperparameters. We moreover provide an efficient Hamiltonian-within-Gibbs sampler for our model. As in other recent work, we reduce the computational burden associated with inversion of covariance matrices by exploiting the structure in the covariance matrix. A tensor decomposition of our covariance matrices allows extended Trench and Durbin algorithms to be applied it order to reduce the computational complexity of inversion and hence accelerate computation. A stand-alone R-package implementing these methods using high-performance C++ libraries is available at: https://github.com/ococrook/toeplitz

Figures

Equations

$$x_i(t_j)={\mu }_i(t_j)+{ϵ}_{ij},$$

$${\mu }_k\sim GP(m_k(t),C_k(t,t^{\prime })).$$

$$p(\mathbf{x}|\mathbf{\pi },\mathbf{\theta })=\sum _{k=1}^K{\pi }_{k}F(\mathbf{x}|{\mathbf{\theta }}_k),$$

$$P(z_1,\dots ,z_n|\mathbf{\pi })=\prod _{k=1}^K{\pi }_k^{n_k},$$

$$P(z_1,\dots ,z_n|\alpha )=\frac{\Gamma (\alpha )}{\Gamma (n+\alpha )}\prod _{i=1}^K\frac{\Gamma (n_i+\alpha /K)}{\Gamma (\alpha /K)}.$$

$$P(z_i=k|z_{-i},\alpha )=\frac{n_{-i,k}+\alpha /K}{N-1+\alpha },$$

$$P(z_i=k|z_{-i})\propto \frac{n_{-i,k}+\alpha /K}{N-1+\alpha }F({\mathbf{x}}_i|{\mathbf{\theta }}_k).$$

$$\mathbf{\pi }|z_1,...,z_n,\alpha \sim Dir(\alpha /K+n_1,...,\alpha /K+n_K).$$

$$P(z_i=k|\mathbf{\pi })\propto {\pi }_{k}F({\mathbf{x}}_i|{\mathbf{\theta }}_k).$$

$$C_k(t_i,t_j)=a_k^2\exp \left(-\frac{\Vert t_i-t_j{\Vert }_2^2}{l_k}\right).$$

$$X_k|{\mathbf{\mu }}_k,{\sigma }_k\sim \mathcal{N}({\mathbf{\mu }}_k,{\sigma }_k^2{I}_D).$$

$$\begin{array}{rl}{X}_k(t_1),...,X_k(t_D)& |{\mu }_k,{\sigma }_k\sim \\ & \mathcal{N}({\mu }_k(t_1),...,{\mu }_k(t_D),...,{\mu }_k(t_1),...,{\mu }_k(t_D),{\sigma }_k^2{I}_{n_{k}D}),\end{array}$$

$${\mu }_k(t_1),...,{\mu }_k(t_D),...,{\mu }_k(t_1),...,{\mu }_k(t_D)|a_k,l_k\sim \mathcal{N}(0,C_k),$$

$$C_k=\left[\begin{array}{cccc}{A}_k& A_k& \dots & A_k\\ A_k& A_k& \dots & A_k\\ ⋮& ⋮& \ddots & ⋮\\ A_k& A_k& \dots & A_k\end{array}\right].$$

$$\begin{array}{r}{X}_k(t_1),...,X_k(t_D)|{\mathbf{\theta }}_k\sim \mathcal{N}(0,C_k+{\sigma }_k^2{I}_{n_{k}D}),\end{array}$$

$$P(z_i=k|z_{-i})\propto \frac{n_{-i,k}+\alpha /K}{K-1+\alpha }\int p({\mathbf{x}}_i|{\mu }_k)p({\mu }_k|{\mathbf{\theta }}_{\mathbf{k}},X_{-i,k}(\tau ))\mathrm{d}{\mu }_k,$$

$$C=\left[\begin{array}{cccc}A+{\sigma }^2{I}_D& A& \dots & A\\ A& A+{\sigma }^2{I}_D& \dots & A\\ ⋮& ⋮& \ddots & ⋮\\ A& A& \dots & A+{\sigma }^2{I}_D\end{array}\right].$$

$$C={\sigma }^2{I}_{nD}+B,$$

$$B=J_n\otimes A,$$

$$Q=(I_D+{\sigma }^{-2}nA),$$

$$\begin{array}{rl}{C}^{-1}& ={\sigma }^{-2}{I}_{nD}-{\sigma }^{-4}{J}_n^T\otimes (ZA)\\ & ={\sigma }^{-2}{I}_{nD}-\frac{1}{n{\sigma }^2}{J}_n^T\otimes (I-Z)\end{array}$$

$$\det (C)=({\sigma }^2{)}^{nD}\det (I_D+{\sigma }^{-2}nA).$$

$$p({\mathbf{x}}_i|\mathbf{\pi },\mathbf{\theta })=\sum _{k=1}^K{\pi }_{k}F({\mathbf{x}}_i|{\mathbf{\theta }}_k{)}^{{\varphi }_i}G(x_i|\Phi {)}^{1-{\varphi }_i},$$

$$\begin{array}{r}{X}_k(t_1),...,X_k(t_D)|{\mathbf{\theta }}_{\mathbf{k}}\sim N(0,C_k+{\sigma }_k^2{I}_{n_{k}D}).\end{array}$$

$$\begin{array}{rl}\log p(X_k|\tau ,& {\mathbf{\theta }}_{\mathbf{k}})\\ & =-\frac{1}{2}{X}_k(\tau ){\left(C_k+{\sigma }_k^2{I}_{n_{k}D}\right)}^{-1}{X}_k(\tau {)}^T-\frac{1}{2}\log |C_k+{\sigma }_k^2{I}_{n_{k}D}|-\frac{n_{k}D}{2}\log 2\pi .\end{array}$$

$$\Lambda =\frac{p({\theta }_{t+1}|X_k(\tau ))}{p({\theta }_t|X_k(\tau ))}=\frac{p(X_k(\tau )|{\theta }_{t+1})p_0({\theta }_{t+1})}{p(X_k(\tau )|{\theta }_t)p_0({\theta }_t)}.$$

$${\left\{H_{ik}^{(t)}=-\sum _{k=1}^K{p}_{ik}^{(t)}\log p_{ik}^{(t)}\right\}}_{t=1}^T,$$

$$H_{ik}\approx \frac{1}{T}\sum _{t=1}^T{H}_{ik}^{(t)}.$$

$$X_k|{\mathbf{\mu }}_k,{\sigma }_k\sim \mathcal{N}({\mathbf{\mu }}_k,{\sigma }_k^2{I}_D),$$

$${\mathbf{\mu }}_k|a_k,l_k\sim GP(0,C_k).$$

Full text

1 Introduction

For a protein to make appropriate interactions with binding partners and substrates, it must localise to the correct sub-cellular compartment (Gibson,, 2009). Furthermore, there is mounting evidence implicating aberrant protein localisation in disease, including cancer and obesity (Olkkonen and Ikonen,, 2006, Laurila and Vihinen,, 2009, Luheshi et al.,, 2008, De Matteis and Luini,, 2011, Cody et al.,, 2013, Kau et al.,, 2004, Rodriguez et al.,, 2004, Latorre et al.,, 2005, Shin et al.,, 2013, Siljee et al.,, 2018). Mapping the sub-cellular location of proteins using high-resolution spatial proteomic approaches are thus of high utility in the characterisation of therapeutic targets and in determining pathobiological mechanisms (Cook and Cristea,, 2019). To interrogate the sub-cellular locations of thousands of proteins per experiment, recent advances in high-throughput spatial proteomics (Christoforou et al.,, 2016, Mulvey et al.,, 2017, Geladaki et al.,, 2019), followed by rigorous data analysis (Gatto et al.,, 2010) can be applied. The methodology relies on the observation that organelles, macro-molecular complexes and, more generally, sub-cellular niches are characterised by density-gradient profiles without the necessity to create homogeneous preparations of major sub-cellular components (De Duve and Beaufay,, 1981).

Mass spectrometry(MS)-based spatial proteomics experiments begin with gentle lysis of the cell in such a way that maintains the integrity of their organelles. A diverse set of methods are available to separate cellular content, including equilibrium density separation (Dunkley et al.,, 2004, 2006, Christoforou et al.,, 2016) and differential centrifugation (Itzhak et al.,, 2016, Geladaki et al.,, 2019, Orre et al.,, 2019), amongst others (Parsons et al.,, 2014, Heard et al.,, 2015). In the LOPIT (Dunkley et al.,, 2004, 2006, Sadowski et al.,, 2006) and hyperLOPIT (Christoforou et al.,, 2016, Mulvey et al.,, 2017) approaches, cell lysis is proceeded by the separation of sub-cellular components along a continuous density gradient based on their buoyant density. Discrete fractions along this gradient are then collected, and protein distributions revealing organelle specific correlation profiles within the fractions are achieved using high accuracy MS.

Sophisticated data analysis methods for spatial proteomics have been developed (Breckels et al.,, 2013, Gatto et al., 2014a, , Breckels et al., 2016a, , Crook et al.,, 2018, Gatto et al.,, 2019), along with detailed work flows (Breckels et al., 2016b, ), which have benefited from significant contributions to the R programming language (R Core Team,, 2017) and the Bioconductor project (Gentleman et al.,, 2004, Huber et al.,, 2015) through MS and proteomics packages (Gatto and Lilley,, 2012, Gatto et al., 2014b, ). Applications of such methods have led to organelle-specific localisation information of proteins in many systems (Dunkley et al.,, 2006, Tan et al.,, 2009, Hall et al.,, 2009, Breckels et al.,, 2013), including mouse pluripotent stem cells (Christoforou et al.,, 2016) and cancer cell lines (Thul et al.,, 2017). MS based spatial proteomics has gained in popularity in recent years with several recent applications across many different cell lines (Christoforou et al.,, 2016, Beltran et al.,, 2016, Jadot et al.,, 2017, Itzhak et al.,, 2017, Mendes et al.,, 2017, Hirst et al.,, 2018, Davies et al.,, 2018, Orre et al.,, 2019, Nightingale et al.,, 2019). This motivates the development of a unified statistical framework for spatial proteomics. Furthermore, with the absence of a mechanistic model for the data, quantifying uncertainty in systems biology applications is of paramount importance and, as yet, such a model has not be applied to existing datasets (Kirk et al.,, 2015).

The current goal of computational methods is to assign proteins with unknown localisation to known sub-cellular niches. It is important to note, however, that not all proteins can be robustly assigned to single locations, since many proteins function in multiple cellular compartments, they may reside in uncharacterised organelles or they may translocate between multiple locations all leading to uncertainty in assignments. Recently, Crook et al., (2018) demonstrated the importance of uncertainty quantification in spatial proteomics analysis. This study developed a generative model of these data and computed posterior distributions of protein localisation probabilities demonstrating the variety of reasons for uncertain protein localisation. This study, however, failed to model known features of the biochemical fractionation process. Sub-cellular niches are not purified in a single fraction along the density-gradient but have distinct quantitative profiles reflective of their biophysical properties. Since the exact quantitative profiles depends heavily on experimental design, the functions describing these profiles are unknown. This motivates a non-parametric Bayesian approach to analysing these data in order to infer the unknown functions and quantify the uncertainty in these functions.

We assume each quantitative protein profile can be described by some unknown function, with the uncertainty in this function captured using a Gaussian process (GP) prior. Each sub-cellular niche is described by distinct density-gradient profiles, which display a non-linear structure with no particular parametric assumption being suitable. The contrasting density-gradient profiles are captured as components in a mixture of Gaussian process regression models. Gaussian process regression models have been applied extensively and we refer to Rasmussen, (2004) and Rasmussen and Williams, (2006) for the general theory. In molecular biology and functional genomics the focus of many applications has been on expression time-series data, where sophisticated models have been developed (Kirk and Stumpf,, 2009, Cooke et al.,, 2011, Kalaitzis and Lawrence,, 2011, Kirk et al.,, 2012, Hensman et al.,, 2013, Strauß et al.,, 2019). We remark that many of these applications consider unsupervised clustering problems. In contrast, here we have (partially) labelled data (proteins with location known prior to our experiments) and so we may consider semi-supervised approaches. We explore inference of GP hyperparameters in two ways: firstly, an empirical Bayes approach in which the hyperparameters are optimised by maximising a marginal likelihood; secondly, by placing priors over these GP hyperparameters and performing fully Bayesian inference using labelled and unlabelled data.

A number of computational aspects need to be considered if inference is to be applied to spatial proteomics data. The first is that correlation in the GP hyperparameters can lead to slow exploration of the posterior, thus we use Hamiltonian evolutions to propose global moves through our probability space (Duane et al.,, 1987) avoiding random walk nature evident in traditional symmetric random walk proposals (Metropolis et al.,, 1953, Beskos et al.,, 2013). Hamiltonian Monte-Carlo (HMC) has been explored previously for hyperparameter inference in GP regression (Williams and Rasmussen,, 1996), and here we show that HMC can be up to an order of magnitude more efficient than a Metropolis-Hastings approach. Furthermore, a particular costly computation in our model is the computation of the marginal likelihood (and its gradient) associated with each mixture component, which involves the inversion of a large covariance matrix - even storage of such matrix can be challenging. We demonstrate that a tensor decomposition of the covariance matrix allows application of fast matrix algorithms for covariance inversion and low memory storage (Zhang et al.,, 2005).

2 Methods

2.1 Model specification

In our experiment, we make discrete observations along a continuous density-gradient, ${\bf x}_{i}=\left[x_{i}(t_{1}),...,x_{i}(t_{D})\right]$, where $x_{i}(t_{j})$ indicates the measurement of protein $i$ at fraction $t_{j}$ along the gradient. We assume that protein intensity $x_{i}$ varies smoothly with the distance along the density-gradient. We further assume observations are equally spaced. Thus, the regression model for each protein is

[TABLE]

where $\mu_{i}$ is an unknown deterministic function of space and $\epsilon_{ij}$ a noise variable. We assume that $\epsilon_{ij}\sim_{iid}\mathcal{N}(0,\sigma_{i}^{2})$, for simplicity and remark that more elaborate noise models could be chosen but at additional computational cost and greater model complexity. Proteins are grouped together according to their sub-cellular localisation, with all proteins associated with sub-cellular niche $k=1,...,K$ sharing the same regression model; that is, $\mu_{i}=\mu_{k}$ and $\sigma_{i}=\sigma_{k}$. For clarity, we refer to sub-cellular structures, whether that be organelles, vesicles or large multi-protein complexes, as components. Thus proteins associated with component $k$ can be modelled as i.i.d draws from a multivariate Gaussian random variable with mean vector $\bm{\mu}_{k}=\left[\mu_{k}(t_{1}),...,\mu_{k}(t_{D})\right]$ and covariance matrix $\sigma^{2}_{k}I_{D}$. To perform inference for unknown $\bm{\mu}_{k}$, as is typical for spatial correlated data (Gelfand et al.,, 2005, Steel and Fuentes,, 2010), we specify a Gaussian Process prior for each $\mu_{k}$

[TABLE]

Each component is thus captured by a Gaussian process regression model and the full complement of proteins as a finite mixture of Gaussian process regression models.

2.2 Finite mixture models

This section provides a brief review of finite mixture models (see, for example (Lavine and West,, 1992, Fraley and Raftery,, 2007) for more details) . Finite mixture models are of the form,

[TABLE]

where $K$ is the number of mixture components, $\pi_{k}$ are the mixture proportions, and $F({\bf x}|\bm{\theta}_{k})$ are the component densities. We assume each component density to have the same parametric form, but with component specific parameters, $\bm{\theta}_{k}$. We denote the prior for these unknown component parameters by $G_{0}(\bm{\theta})$. We suppose that we have a collection of $n$ data points, $X=\{{\bf x}_{1},\ldots,{\bf x}_{n}\}$ that we seek to model using Equation (3). We associate with each of these data points a component indicator variable, $z_{i}\in\{1,\ldots,K\}$, which indicates which component generated observation ${\bf x}_{i}$. Given the mixing proportions, the joint prior distribution of these indicators is multinomial with parameter vector $\bm{\pi}=[\pi_{1},\ldots,\pi_{K}]$,

[TABLE]

where $n_{k}$ is the number of data points $x_{i}$ for which $z_{i}=k$. If we assign the mixture proportions a symmetric Dirichlet prior with concentration parameter $\alpha/K$, then we may marginalise the $\pi_{k}$ in order to yield the following joint distribution for the indicators (Murphy,, 2012),

[TABLE]

For Gibbs sampling, we require the conditional priors for a single indicator, $z_{i}$, given all of the others, $z_{-i}$. These are given by (Murphy,, 2012),

[TABLE]

where $n_{-i,k}$ is the number of observations, excluding ${\bf x}_{i}$, that are associated with component $k$. If we are given the parameters, $\bm{\theta}_{k}$, associated with each of the components then we may combine the above conditional priors with the likelihoods, $F({\bf x}_{i}|{\bf\bm{\theta}}_{k})$, in order to obtain the conditional posterior:

[TABLE]

An alternative to integrating out the mixture proportions is to sample them directly from the posterior, which leads to increased posterior variance (Gelfand and Smith,, 1990, Casella and Robert,, 1996) but can be computational advantageous. Conjugacy of the Dirichlet prior and multinomial likelihood means that the posterior distribution of the mixing proportions is also Dirichlet,

[TABLE]

In this situation the conditional posterior becomes

[TABLE]

2.3 Gaussian Process priors

A Gaussian Process (GP) is a continuous stochastic process such that any finite collection of these random variables is jointly Gaussian. A Gaussian Process prior is uniquely specified by a mean function $m$ and covariance function $C$, which determine the mean vectors and covariance matrices of the associated multivariate Gaussian distributions. To elaborate, assuming a GP prior for $\mu_{k}$ means that for indices $t_{1},...,t_{D}$, the joint prior of $\bm{\mu}_{k}=\left[\mu_{k}(t_{1}),...,\mu_{k}(t_{D})\right]^{T}$, is multivariate Gaussian with mean vector $\bm{m}_{k}=\left[m_{k}(t_{1}),...,m_{k}(t_{D})\right]$ and covariance matrix $C_{k}(i,j)=C_{k}(t_{i},t_{j})$. Given no prior belief about symmetry or periodicity in our deterministic function, we assume our GP is centred with squared exponential covariance function

[TABLE]

2.4 Marginalising the unknown function

Having adopted a GP prior with component specific parameters $a_{k}$ and $l_{k}$ for each unknown function $\mu_{k}$, we let observations associated with component $k$ be denoted by $X_{k}=\{x_{1},...,x_{n_{k}}\}$. Our model tells us that

[TABLE]

Then, we can write this as

[TABLE]

where $\mu_{k}(t_{1}),...,\mu_{k}(t_{D})$ is repeated $n_{k}$ times. Our GP prior tell us

[TABLE]

where $C_{k}$ is an $n_{k}D\times n_{k}D$ matrix. This matrix is organised into $n_{k}\times n_{k}$ square blocks each of size $D$. The $(i,j)^{th}$ block of $C_{k}$ being $A_{k}$, where $A_{k}$ is the covariance function for the $k^{th}$ component evaluated at $\tau=\{t_{1},...,t_{D}\}$.

[TABLE]

Letting $\bm{\theta}_{k}=\left\{a_{k},l_{k},\sigma^{2}_{k}\right\}$, we can then marginalise $f_{k}$ to obtain,

[TABLE]

thus avoiding inference of $\mu_{k}$. Let $X_{k}(\tau)$ denote the vector of length $n_{k}\times D$ equal to $\left(x_{1}(t_{1}),...,x_{1}(t_{D}),\dots,x_{n_{k}}(t_{1}),...,x_{n_{k}}(t_{D})\right)$. Then we may rewrite equation 7 by marginalising $\mu_{k}$ to obtain:

[TABLE]

where $X_{-i,k}(\tau)$ is equal to $X_{k}(\tau)$ with observation $i$ removed.

2.5 Tensor decomposition of the covariance matrix for fast inference

Our covariance matrix has a particularly simple structure allowing us to exploit extended Trench and Durbin algorithms for fast matrix computations (Zhang et al.,, 2005). Recall we are interested in the inversion of matrices of the following form

[TABLE]

Note that $A$ is a positive symmetric matrix of size $D\times D$ and furthermore it is Toeplitz (constant diagonal and perisymmetric). Let $J_{n}$ denote an $n\times n$ matrix of ones. It is clear that we can write $C$ in the following form:

[TABLE]

where

[TABLE]

and $\otimes$ denotes the Kronecker (tensor) product. Let us write

[TABLE]

which is a $D\times D$ Toeplitz matrix, for which the inverse and determinant can be inverted in $O(D^{2})$ operations (Durbin,, 1960, Trench,, 1964). If we denote the inverse of $Q$ by $Z$ it follows (see supplementary for full derivation) that:

[TABLE]

and

[TABLE]

Thus, the inversion of $C$ requires only the inversion of a $D\times D$ matrix, which can be performed in $O(D^{2})$ computations, this should be compared with a naïve inversion of $C$ requiring $O((nD)^{3})$ computations, which represents significant savings. The determinant can also be obtained in $O(D^{2})$ operations (Zhang et al.,, 2005). The method of Zhang et al., (2005) can be seen as a special case of our situation where $n=1$. Step by step algorithms for computing this inverse and determinant can also be found in the supplementary materials. We note that Strauß et al., (2019) also exploit the block matrix structure of the covariance matrix efficiently, using a more general approach to compute block matrix inversions and determinants that works also in the case of hierarchical GP models, for which Hensman et al., (2013) had found an alternative way of performing efficient likelihood computations.

2.6 Sampling the underlying function

Whilst it is often mathematically convenient to marginalise the unknown function $\mu_{k}$ from a computational perspective it is not always advantageous to do so. To be precise, marginalising $\mu_{k}$ induces dependencies among the observations; that is, we cannot exploit the conditional independence structure given the underlying function $\mu_{k}$. After marginalising, Gibbs moves must be made sequentially for each protein in turn and this can slow down computation.

The alternative approach is to sample the underlying function and exploit conditional independence. Once a sample is obtained from the GP posterior on $\mu_{k}$, conditional independence allows us to compute the likelihood for all proteins at once, exploiting vectorisation. If there are a particularly large number of observation in each component it is also possible to parallelize computation over the components $k=1,...,K$.

2.7 Modelling outliers

Crook et al., (2018) demonstrated that many proteins are not captured well by any known sub-cellular component. This could be because of yet undiscovered biological novelty, technical variation or a manifestation of some proteins residing in multiple localisations. Modelling outliers in mixture models can be challenging (Hennig,, 2004, Cooke et al.,, 2011, Coretto and Hennig,, 2016). Here, we take the approach of Crook et al., (2018). Briefly, we introduce a further binary latent variable $\phi$ so that for each protein ${\bf x}_{i}$ we have a $\phi_{i}$ indicating whether ${\bf x}_{i}$ is modelled by one of the known components or an outlier component. The augmented model becomes the following

[TABLE]

where $G$ is density of the outlier component. In our case, we specify $G$ as the density of a multivariate T distribution with degrees of freedom $\kappa=4$, mean ${\bf M}$ and scale matrix $V$. ${\bf M}$ is taken as the empirical global mean of the data and the scale matrix $V$ as half the empirical covariance of the data. These choices are motivated by considering a Gaussian component with the same mean and covariance but with heavier tails to better capture dispersed proteins. We remark that other choices of G and parameters may be suitable and can be tailored to the application at hand. In typical Bayesian fashion, we specify a prior for $\phi$ as $p_{0}(\phi_{i}=0)=\epsilon$, where $\epsilon\sim\mathcal{B}(u,v)$. All hyperparameter choices are stated in the appendix.

2.8 Gaussian process hyperparameter inference

2.8.1 Supervised approach: optimising the hyperparameters

Inference of the hyperparameters $\bm{\theta_{k}}$ can be dealt with in several ways. The first is to learn them using only the labelled data (i.e. data that pertains to proteins with well documented sub-cellular locations). Using the labelled data for each component constitutes maximise the marginal likelihood of the hyperparameters with respect to the data. These hyperparameters are then fixed throughout the inference of the unlabelled data. The marginal likelihood can be obtained quickly by recalling that

[TABLE]

Thus the log marginal likelihood is given by

[TABLE]

For convenience of notation set $\hat{C}_{k}=C_{k}+\sigma^{2}_{k}I_{n_{k}D}$. To maximise the marginal likelihood given equation 25, we find the partial derivatives with respect to the parameters (Rasmussen,, 2004). Hence, we can use a gradient based optimisation procedure. Positivity constraints on $a^{2}_{k},l_{k},\sigma_{k}^{2}$ are dealt with by re-parametrisation and so, dropping the dependence on $k$ for notational convenience, and abusing notation, we set $l=\exp(\theta_{1})$, $a^{2}=\exp(2\theta_{2})$ and $\sigma^{2}=\exp(2\theta_{3})$.

Application of the quasi-Newton L-BFGS algorithm (Liu and Nocedal,, 1989) for numerical optimisation of the marginal likelihood with respect to the hyperparameters is now straightforward. The L-BFGS can only find a local optimum and so we initialise over a grid of values. We terminate the algorithm when successive iterations of the gradient are less than $10^{-8}$. We make extensive use of high performance R packages to interface with C++ (Eddelbuettel and Francois,, 2011, Eddelbuettel and Sanderson,, 2014).

2.8.2 Semi-supervised hyperparameter inference

The advantage of adopting a Bayesian approach to hyperparameter inference is that we can quantify uncertainty in these hyperparameters. Uncertainty quantification in GP hyperparameter inference is important, since different hyperparameters can have a strong effect on the GP posterior (Rasmussen,, 2004). Furthermore, we consider a semi-supervised approach to hyperparameter inference. By a semi-supervised approach we mean that a posterior distribution for the hyperparameters can be inferred using both the labelled and unlabelled data, rather than just the labelled data.

Consider at some iteration of our MCMC algorithm the data associated to the $k^{th}$ component $X_{k}$. We can partition this data into the unlabelled (U) and labelled data (L); in particular, $X_{k}=\left[X^{(L)}_{k},X^{(U)}_{k}\right]$. To clarify, the indicators $z_{i}$ are known for $X^{(L)}_{k}$ prior to any inference, whilst allocations $z_{i}$ for $X^{(U)}_{k}$ are sampled at each iteration of our MCMC algorithm. If we believe our labelled data $X^{(L)}_{k}$ are true representatives of the distribution of that component, it is computationally advantageous just to consider the labelled data when performing hyperparameter inference. However, there could be a sampling bias in the labelled data and so the labelled data alone is insufficient to explain the variability in the data. A semi-supervised approach allows the posterior distribution of the hyperparamters to reflect the uncertainty in the component allocations $z_{i}$ and therefore improve our abilities to predict allocations and quantify uncertainty in allocations.

2.8.3 Semi-Supervised approach: hyperparameter inference using MH

In a Bayesian framework, we treat the hyperparameters as random variables and place hyperpriors overs them. Positivity constraints motivate working with the $\log$ of the hyperparameters and using, for example, standard normal priors (Neal,, 1997). Unfortunately loss of conjugacy between the prior on the hyperparameters $G_{0}(\bm{\theta})$ and the likelihood $F({\bf x}|\bm{\theta})$ is unavoidable, and hence we use a Metropolis-Hastings step or Hamiltonian Monte-Carlo step for inference. The Metropolis-Hastings sampler can be summarised as follows:

Metropolis-Hastings algorithm with random walk proposals: Suppose $\theta_{t}$ is the most recently sampled value. Sample a value $\xi\sim N(0,1)$, setting $\theta_{t+1}=\theta_{t}+\xi$ and compute the Metropolis ratio

[TABLE]

This ratio can be computed in $\log$ form using equation 25. Then sample a uniform random number $u\sim U[0,1]$ if $\log(\Lambda)\geq\log(u)$ set $\theta_{t+1}=\theta_{t}+\xi$, otherwise $\theta_{t+1}=\theta_{t}$.

2.8.4 Semi-Supervised approach: hyperparamter inference using HMC

To avoid the random walk nature of the MH sampler, we also consider a Hamiltonian Monte-Carlo approach, which exploits the geometry of the space to provide more efficient proposals (Duane et al.,, 1987, Horowitz,, 1991, Neal et al.,, 2011, Girolami and Calderhead,, 2011). In short, Hamiltonian Monte-Carlo allows us to construct Hamiltonian evolutions $H(\bm{x},\bm{p})$ such that the resulting dynamics efficiently explore a target distribution $p(\bm{x})$. We augment our probability distribution with an auxiliary momentum component $\bm{p}$. An MCMC algorithm can then be constructed to sample from the required distribution, where proposals are made using Hamiltonian evolutions. Full details in the case of the hyperparameters of a Gaussian process are discussed in the supplement. In previous sections, we saw we can exploit a tensor decomposition to accelerate computation of the likelihood and similar formulae are available to accelerate computation of the gradient for use in L-BFGS and Hamiltonian Monte Carlo. These formulae can be found in the supplement.

2.8.5 An overview of the MCMC algorithm for posterior Bayesian computation

In our model $G_{0}(\bm{\theta})$ and $F({\bf x}|\bm{\theta})$ are non-conjugate, which means the integral in equation 16 cannot be obtained analytically. A Gibbs sampling scheme with either an additional Metropolos-Hastings or Hamiltonian Monte Carlo update is used. Each iteration of the MCMC algorithm includes a sampled value for the component indicators, outlier components and current values of the hyperparameters. We also keep track of associated posterior probabilities and marginal likelihoods as appropriate. Furthermore, we can sample the hyperparameters every $T$ iterations of the MCMC algorithm to accelerate computations.

2.9 Summarising uncertainty in posterior localisation probabilities

Summarising uncertainty quantified by Bayesian analysis in an interpretable way can be challenging. As always, we can summarise uncertainty using credible intervals or regions (Gelman et al.,, 1995). One particularly challenging quantity of interest to summarise is the uncertainty in posterior allocations. Whilst, each individual allocation of a protein to a sub-cellular niche can be summarised by a credible interval it is not clear what is the best way to summarise the posterior over all possible localisations for each individual protein. As in previous work (Crook et al.,, 2018), we propose to summarise this uncertainty in an information-theoretic approach by computing the Shannon entropy of the localisation probabilities (Shannon,, 1948) at each iteration of the MCMC algorithm

[TABLE]

where $p_{ik}^{(t)}$ is the probability that protein $i$ belong to component $k$ at iteration $t$. We can then summarise this by a Monte-Carlo average:

[TABLE]

We note that, larger values of a Shannon entropy correspond to greater uncertainty in allocations.

3 Results

3.1 Case Study I: Drosophila melanogaster embryos

3.1.1 Application

The first case study is the Drosophila melanogaster (common fruit fly) embryos (Tan et al.,, 2009), in which we compare the supervised and semi-supervised approaches for updating the model hyperparameters. In particular, we explore the effect on the component specific noise term $\sigma^{2}$, by adopting different inference approaches. For each sub-cellular niche, we learn the hyperparameters by either maximising their marginal likelihood or sampling from their posterior using MCMC. The posterior distribution for the hyperparameters can either be found solely using the labelled data for each component or by making use of labelled and unlabelled data.

Figure 1 demonstrates several phenomena. Reassuringly, the estimates of the noise parameters $\sigma_{k}^{2}$ for $k=1,...,K$ obtained by using the L-BFGS algorithm to maximise the marginal likelihood coincide with the posterior distributions of the noise parameters, inferred using only the labelled data for each component. However, when we perform inference in a semi-supervised way, by using both the labelled and unlabelled data to make inferences, we make several important observations.

Firstly, in many cases, the posterior using both the labelled and unlabelled data is shifted right towards [math]. Recalling that we are working with the log of the hyperparameters, this indicates that the noise parameters is smaller when solely using the labelled data. This is likely a manifestation of experimental bias, since it is reasonable to believe that proteins with known prior locations are those which have less variable localisations and are therefore easier to experimentally validate. A semi-supervised approach is able to overcome these issues, by adapting to proteins in a dense region of space. In some cases the shift is pronounced, with posteriors of the parameters using labelled and unlabelled data found in the tails of the posterior only using the labelled distribution. Furthermore, we notice shrinkage in the posterior distribution of the noise parameter in the semi-supervised setting. The reduction in variance reduces our uncertainty about the underlying true value of $\sigma_{k}^{2}$ for $k=1,...,K$. This variance reduction is observed in most cases even when these is little difference in the mean of the posteriors.

The primary goal of spatial proteomics is to predict the localisation of unknown proteins from data. Our modelling approach allows the allocation probability of each protein to each component to be used to predict the localisation of unknown proteins. Proteins may reside in multiple locations and some sub-cellular niches are challenging to separate because of confounding biochemical properties, leading to uncertainty in a proteins localisation. Thus adopting a Bayesian approach and quantifying this uncertainty is of great importance. Our methods allow point-estimates as well as interval estimates to be obtained for the posterior localisation probabilities. Figure 2 demonstrates the results of applying our method. Each protein in this PCA plot is scaled according to mean of the Monte-Carlo samples from the posterior localisation probability.

Further visualisation of the model and data are possible. We plot two representative example of gradient-density profiles for two components the endoplasmic reticulum (ER) and the nucleus, in figure 3. We plot both the labelled proteins, in colour, which were assigned to each component before our analysis. In grey, for both components, we plot the unlabelled proteins which have been allocated to these components probabilistically. We observe that they have the same gradient-density shape as the labelled proteins - in line with our beliefs about the underlying biology: that proteins from the same components should co-fractionate and therefore have similar density gradient profiles. In addition, we overlay the posterior predictive distribution for these components and observe they represent the data well.

3.1.2 Sensitivity analysis for hyper-prior specification

We use the Drosophila melanogaster dataset to test for sensitivity of the hyper-prior specification. To test for sensitivity, we see if predictive performance is affected by changes in the choice of hyper-prior. The following cross-validation schema assesses whether predictive performance is affected by choice of hyper-prior. We split the labelled data for each experiment into class-stratified training $(80\%)$ and test $(20\%)$ partitions, with the separation formed at random. The true classes of the test profiles are withheld from the classifier, whilst MCMC is performed. This $80/20$ data stratification is performed $100$ times in order produce a distribution of scores. We compare the ability of the methods to probabilistically infer the true classes using the quadratic loss, also referred to as the Brier score (Gneiting and Raftery,, 2007). Thus a distribution of quadratic losses is obtained for each method, with the preferred method minimising the quadratic loss. Each method is run for $10,000$ MCMC iterations with $1000$ iterations for burn-in. We vary the mean of the standard normal hyper-prior for each hyperparameter in turn for a grid of values $\tilde{m}=(0,-1,-2,-3,-4)$, keeping the hyper-prior for the other variable held the same as a standard normal distribution. The results are displayed in figure 4.

We observe only minor sensitivity to the choice of hyper-prior, with no significant difference in performance noted (KS test, threshold = $0.01$). Sensitivity analysis for hyperparameters of GPs is vital, since these hyperparameters have a strong effect on the posterior of the GP (Rasmussen,, 2004). The observed lack of sensitivity in our case is advantageous, since prior information can be included without fear of over fitting. However, practitioners should always take care when specifying priors, especially for variance/covariance parameters as many authors have noted sensitivity of Bayesian models to these parameters (Gelman et al.,, 1995, Lunn et al.,, 2000, Gelman et al.,, 2006, Wang and Dunson,, 2011, Schuurman et al.,, 2016)

3.2 Case Study II: mouse pluripotent embryonic stems cells

3.2.1 Application

Our main case study is the mouse pluripotent E14TG2a stem cell dataset of Christoforou et al., (2016). This dataset contains 5032 quantitative protein profiles, and resolves 14 sub-cellular niches. We first plot the density-gradient profiles of the marker proteins for each sub-cellular niche in figure 5. We fit a Gaussian process prior regression model for each sub-cellular niche with the hyperparameters found by maximising the marginal likelihood.

A table of unconstrained log hyperparameter values found by maximising the marginal likelihood is found in the supplement. Alternatively, placing standard normal priors on each of the log hyperparameters and using a Metropolis-Hastings update we can infer the distributions over these hyperparameters. We perform $20,000$ iterations for each sub cellular niche and discard $15,000$ iterations for burn-in and proceed to thin the remaining samples by $20$. We summarise the Monte-Carlo sample by the expected value as well as the $95\%$ equi-tailed credible interval, which can also be found in the supplement.

We go further to predict proteins with unknown localisation to annotated components using our proposed mixture of GP regression models. As before, we adopt a semi-supervised approach to hyperparameter inference. Again we place standard normal hyper-priors on the log of the hyperparameters. We run our MCMC algorithm for $20,000$ iterations with half taken as burnin and thin by $5$, as well as using HMC to update the hyperparameters. The PCA plot in figure 6 visualises our results. Each pointer represent a single protein and is scaled either to the probability of membership to the coloured component (left) or scaled with the Shannon entropy (right). In these plots we observe regions of high-probability and confidence to each organelle, as well as obtaining a global view of uncertainty. In this example, we observe regions of uncertainty, as measured by the Shannon entropy, concentrating where components overlap. We also observe uncertainty in regions where there is no dominant component. This Bayesian analysis provides a wealth of information on the global patterns of protein localisation in mouse pluripotent embryonic stem cells.

3.3 Assessing predictive performance

We compare the predictive performance of the methods proposed here, as well as against the fully Bayesian TAGM model of Crook et al., (2018), where sub-cellular niches are described by multivariate Gaussian distributions rather than GPs. The following cross-validation schema is used to compare the classifiers. We split the labelled data for each experiment into class-stratified training $(80\%)$ and test $(20\%)$ partitions, with the separation formed at random. The true classes of the test profiles are withheld from the classifier, whilst MCMC is performed. This $80/20$ data stratification is performed $100$ times in order produce a distribution of scores. We compare the ability of the methods to probabilistically infer the true classes using the quadratic loss, also referred to as the Brier score (Gneiting and Raftery,, 2007). Thus a distribution of quadratic losses is obtained for each method, with the preferred method minimising the quadratic loss. Each method is run for $10,000$ MCMC iterations with $1000$ iterations for burn-in. For fair comparison we held priors the same across all datasets. Prior specifications are stated in the supplement.

We compare across $5$ different spatial proteomics datasets across three different organisms. The datasets we compare our methods on are Drosophila melanogaster embryos from Tan et al., (2009), the mouse pluripotent embroyonic stem cell dataset of Christoforou et al., (2016), the HeLa cell line dataset of Itzhak et al., (2016), the mouse primary neuron dataset of Itzhak et al., (2017) and finally a CRISPR-CAS9 knock-out coupled to spatial proteomics analysis dataset (AP5Z1-KO1) of Hirst et al., (2018). The results are found in figure 7.

We see that our in four out five datasets there is an improvement of the GP models over the TAGM model (Kolmogorov-Smirnov (KS) two-sample test $p<0.0001$ ), because the GP model is provided with more explicit correlation structure of the data. The empirical Bayes slightly method outperforms the fully Bayesian approach in three of the data sets ((KS) two-sample test $p<0.01$ ). These are the mouse pluripotent embryonic stem cell dataset, the HeLa data set of Itzhak et al., (2016) and the HeLA AP5Z1 knock-out dataset of Hirst et al., (2018). We observe, the TAGM model outperforms the GP methods in the Itzhak et al., (2016) dataset. The authors of this study used differential centrifugation to separate cellular content and curated a “large protein complex” class. This class could contain multiple sub-cellular structures such as ribosomes, as well as cytosolic and nuclear proteins. In any case, our modelling assumptions are violated in both models and this is issue is exacerbated by parametrising the covariance structure. One solution to this would be to model this mixture of large protein complexes as its own class. However, as this class contains a quite diverse set of sub-cellular compartments, it is difficult to predict behaviour. This class could be itself a mixture of GPs, however the number of components of the class would be unknown and this would have to be carefully modelled, perhaps using reversible jump methods (Richardson and Green,, 1997) or Dirichlet process approaches (Escobar and West,, 1995).

4 Discussion

This article presents semi-supervised non-parametric Bayesian methods to model spatial proteomics data. Sub-cellular niches display unique signatures along a density-gradient and we exploit this information to construct GP regression models for each niche. The full complement of sub-cellular proteins in then described as mixture of GP regression models, with outliers captured by an additional component in our mixture. This provides cell biologist with a fully Bayesian method to analyse spatial proteomics data in the non-parametric framework that more closely reflects the biochemical process used to generate the data. This greatly increases model interpretation and allows us to make more biological sound inferences from our model.

We compared the proposed semi-supervised models to the state-of-the-art model on $5$ different spatial proteomics datasets. Modelling the correlation structure along the density-gradient leads to competitive predictive performance over state-of-the-art models. Empirical Bayes procedures perform either equally well or better than the fully Bayesian approach, at the loss of uncertainty quantification in the hyperparameters. Though this performance improvement should not be over interpreted, since cross-validation assessment is only performed on the labelled data and will not reflect any biased sampling mechanisms that could be at play.

To accelerate computation in our model, we note that the structure of our covariance matrix admits a tensor decomposition, which can be exploited so that fast algorithms for matrix inversion of toeplitz matrices can be employed. These decomposition can then be used to derive formulae for fast computation of the likelihood and gradient of a GP. A stand-alone R-package implementing these methods using high-performance C++ libraries is available at https://github.com/ococrook/toeplitz. These algorithms and associate formulae are useful to those outside the spatial proteomics community to anyone using GPs with equally spaced observations, even in the unsupervised case.

We demonstrated that in the presence of labelled data there are two approach to hyperparameter inference. This first is to use empirical-Bayes to optimise the hyperparameters; the other a fully-Bayesian approach, taking into account the uncertainty in these hyperparameters. We propose to use HMC to update these hyperparameters, since highly correlated hyperparameters can induce high autocorrelation and exacerbate issues with random-walk MH updates. We demonstrate that, in the situation presented here, HMC updates can be up to an order of magnitude more efficient than MH updates. We further explored the sensitivity of our model to hyper-prior specification, which gives practitioners good default choices.

In two case-studies, we highlighted the value of taking a semi-supervised approach to hyperparameter inference, allowing us to explore the uncertainty in our hyperparameters. In a fully Bayesian approach the uncertainty in the hyperparamters is reflected in the uncertainty of the localisation of proteins to components. Quantifying uncertainty provide cell biologists with a wealth of information to make quantifiable inference about protein sub-cellular localisation.

We plan to disseminate our method via the Bioconductor project (Gentleman et al.,, 2004, Huber et al.,, 2015) and include our code in pRoloc package (Gatto et al., 2014b, ). The pRoloc package includes methods for visualisation, processing data and disseminating code in a unified framework. All spatial proteomics data used here is freely available within the Bioconductor package pRolocdata (Gatto et al.,, 2018).

One potential source of uncertainty in protein localisation is that they can be residents of multiple sub-cellular compartments. We believe that by proposing a model which more closely reflects the underlying biochemical rationale for the experiment we can facilitate models which can infer proteins with multiple locations with greater confidence. This is the subject of further work.

5 Supplementary

5.1 GP Prior, posteriors and predictive distributions

Denoting our unknown regression function $\mu_{k}$, we let observations associated with component $k$ be denoted by $X_{k}=\{x_{1},...,x_{n_{k}}\}$. Our model tells us that

[TABLE]

where

[TABLE]

The posterior distribution follows from normal theory:

[TABLE]

where

[TABLE]

The mean and covariance functions for the posterior predictive distribution are given by:

[TABLE]

5.2 Derivation of tensor-Toeplitz decomposition inverse

Let $J_{n}$ denote a $n\times n$ matrix of ones. Recall that we can write $C$ in the following form:

[TABLE]

where

[TABLE]

and we have denoted $\otimes$ as the Kronecker (tensor) product. Let $e_{n}$ denote a column vector of ones of length $n$. It is easy to see that $J_{n}=e_{n}e_{n}^{T}$. Trivially, we can write $A=I_{D}A$ and this leads to the following factorisation

[TABLE]

where the second equality follows from the mixed-product property of the Kronecker product. Observing that $e_{n}\otimes I_{D}$ is a matrix of size $nD\times D$, and $e_{n}^{T}\otimes A$ is matrix of size $D\times nD$. We thus arrive at the following factorisation:

[TABLE]

which is in the following form

[TABLE]

Matrices of this form have a simple formula for their inverse (Woodbury Identity):

[TABLE]

In our case $R$ is trivially its own inverse and the inverse of $M$ requires only a single computation. Thus the only challenge is to invert $(R^{-1}+VM^{-1}U)$. However, consider the following computations

[TABLE]

Recall that $A$ is a $D\times D$ Toeplitz matrix and so it is easy to see that $R^{-1}+VM^{-1}U$ is also Toeplitz and efficient algorithms exist for inverting them. Denote this inverse by Z and so

[TABLE]

where the last follows from the following computations, denoting $Z^{-1}=Q$

[TABLE]

Thus the inversion of $C$ requires only the inversion of a $D\times D$ matrix, which can be performed in $O(D^{2})$ computations, this should be compared with a naïve inversion of $C$ requiring $O((nD)^{3})$ computations, which represents significant savings. We also need the determinant of $C$ and the calculation is straightforward using an elementary determinant lemma.

[TABLE]

As before the term in the determinant is Toeplitz and efficient algorithm exists for calculating this determinant.

5.3 Matrix algorithms

We state here the require algorithm to invert the covariance matrix $C=\sigma^{2}I_{nD}+J_{n}\otimes A$ for a Toeplitz matrix $A$. The algorithms are a minor modification of the algorithms found in Zhang et al., (2005) to handle the Tensor product.

5.4 Derivative of the marginal likelihood

The derivatives of the marginal likelihood given in equation 25 are given by (Rasmussen,, 2004)

[TABLE]

The partial derivatives of the covariance functions can obtained in a straightforward manner and once evaluated at observations can be structured into blocks just as in equation 17. Letting $\hat{A}_{k}$ be the diagonal blocks of the covariance matrix in equation 17. The corresponding diagonal blocks of derivative are given in equation 47. Blocks not on the diagonal are similar and do not include the derivative with respect to $\theta_{3}$.

[TABLE]

5.5 Tensor decompositions for derivatives of the marginal likelihood

In this appendix we derive formulae for the derivative of the marginal likelihood exploiting the block structure of our matrices. We first make some preliminary manipulations. We set the following notation $\partial_{\theta_{j}}=\frac{\partial}{\partial_{\theta_{j}}}$. First we note that

[TABLE]

We recall the following

[TABLE]

and hence the following is true

[TABLE]

We then note that $\partial_{\theta_{j}}J_{n}=0$ and so the following algebraic manipulations hold

[TABLE]

We recall that

[TABLE]

and so

[TABLE]

It is obvious that

[TABLE]

and so

[TABLE]

Whence it follows that

[TABLE]

Recall that

[TABLE]

where $S_{rs}=\frac{(t_{r}-t_{s})^{2}}{l}$. We now derive formulae for the derivatives of the marginal likelihood and we denote $A\odot B$ has the Hadamard (element-wise) product of matrices $A$ and $B$.

Proposition 1**.**

The derivative of the marginal likelihood in 46 with respect to $\theta_{1}$ is given by

[TABLE]

where

[TABLE]

Proof.

We observe the following manipulations, which follow from our preliminary manipulations

[TABLE]

where the third line follows from the second because

[TABLE]

For the trace term, recall that the trace of a product of two matrices is the sum of the Hadamard product of those two matrices. That is

[TABLE]

Applying the mixed product property, we see that the following manipulations hold

[TABLE]

Hence,

[TABLE]

Thus the derivative of the log marginal likelihood is

[TABLE]

Then we can substitute $ZA=(I-Z)\frac{\sigma^{2}}{n}$ to obtain the required result. ∎

Proposition 2**.**

The derivative of the marginal likelihood in 46 with respect to $\theta_{2}$ is given by

[TABLE]

where

[TABLE]

Proof.

As in the previous proposition we observe:

[TABLE]

For the trace term, as for $\theta_{1}$ we proceed as follows

[TABLE]

Hence,

[TABLE]

Thus the derivative of the log marginal likelihood is

[TABLE]

Then we can substitute $ZA=(I-Z)\frac{\sigma^{2}}{n}$ to obtain the required result. ∎

Proposition 3**.**

The derivative of the marginal likelihood in 46 with respect to $\theta_{1}$ is given by

[TABLE]

where

[TABLE]

Proof.

We note that $\partial_{\theta_{3}}\hat{C}_{k}=2\sigma^{2}I_{nD}$ is a scalar multiple of the identity matrix and thus commutes. Hence, we need only compute $\hat{C}_{k}^{-1}\hat{C}_{k}^{-1}$ and the trace term. Note the following algebraic manipulations:

[TABLE]

The compute the trace we note that the following follows directly from the tensor decomposition of $\hat{C}_{k}^{-1}$:

[TABLE]

Substituting the formulae shows the desired result it now clear. ∎

In practice, we never need to compute or even store the full $nD\times nD$ inverse matrix $C^{-1}$, since we can only need to keep track of summaries of the data matrix rather than the full data matrix itself. This is demonstrated in the following proposition.

Proposition 4**.**

Let

[TABLE]

be a $D\times n$ matrix. Let $Y_{i}=\sum_{j}X_{i,j}$ be the sum of the $i^{th}$ row of X and written concisely $Y=Xe_{n}$, where $e_{n}$ is a $n\times 1$ vector of ones. We write $J_{n}$ to be the $n\times n$ matrix of ones. Let $R$ be any $D\times D$ matrix. Then the following holds

[TABLE]

where $vec(X)$ denotes the vectorisation of $X$; that is, the $Dn\times 1$ vector formed by stacking columns of $X$.

Proof.

Firstly, observe the following standard algebraic manipulations

[TABLE]

Thus, using the above, it follows that

[TABLE]

as required. ∎

5.6 Hamiltonian Monte-Carlo for GP hyperparameters

The Hamiltonian can be decomposed into potential and kinetic energies $H(\bm{x},\bm{p})=U(\bm{x})+K(\bm{p})$. The canonical distribution is then given by:

[TABLE]

The distribution of momentum component is chosen as a Gaussian distribution with diagonal covariance matrix $M=diag(m_{1},...,m_{r})$ and thus the distribution and kinetic energies are given by

[TABLE]

It is easy to see from the canonical distribution that $U(\bm{x})=-\log(p(\bm{x}))$ is the required choice for the potential. In practice, we need to simulate from Hamiltonian dynamics. Hamilton’s equations are given by a coupled system:

[TABLE]

Such a system is called symplectic and thus a numerical schema which is a symplectic integrator is required to simulate the required dynamics (Neal et al.,, 2011). The leapfrog algorithm is the standard choice (MacKay,, 2003). This algorithm does not exactly conserve energy and so a Metropolis accept/reject step is required is remove the induced bias (Beskos et al.,, 2013). An MCMC algorithm can then be constructed to sample from the required distribution, where proposals are made using Hamiltonian evolutions. Recall, we are required to simulate the Hamiltonian evolutions. To simulate an evolution over time $T$, take $L$ steps of size $\delta$ such that $L\delta\geq T$. One step of the leapfrog algorithm of size $\delta$ for Hamilton’s dynamics starting at time $t$ is given by the following

[TABLE]

We can now summarise the HMC algorithm to sample $n$ samples from a target distribution $p(\bm{x})$.

Set $t=0$ 2. 2.

Sample a position value from the prior $\bm{x}_{0}\sim p_{0}$ 3. 3.

Do until $t=n$

1. (a)

Set $t=t+1$ 2. (b)

Sample an initial momentum variable $\bm{p}_{0}\sim p(\bm{p})$ 3. (c)

Set $\bm{x}_{0}=\bm{x}_{t-1}$ 4. (d)

Run algorithm 82 for $L$ step of size $\delta$ and obtain proposal states $\bm{x}_{*}$ and $\bm{p}_{*}$ 5. (e)

Compute the Metropolis ratio

[TABLE] 6. (f)

Sample $u\sim U[0,1]$ if $\Lambda>u$ set $\bm{x}_{t}=\bm{x}_{*}$, else $\bm{x}_{t}=\bm{x}_{t-1}$

We can now specify the details for sampling the hyperparameters of a Gaussian Process with standard normal hyperpriors. Using a squared exponential covariance function and re-parametrising, as before, we first specify our target distribution $p(\bm{x})=p(\bm{\theta}|X(\tau))\propto p(X(\tau)|\bm{\theta})p_{0}(\bm{\theta})$. Now considering

[TABLE]

the first term can be computed by marginalising and is recognised as the marginal likelihood given in equation 25. Recalling that we have a standard normal prior the negative log prior and its gradient is given by is given by

[TABLE]

where $\bm{\theta}=(\theta_{1},\theta_{2},\theta_{3}).$ Hence, we can write down the gradient of the potential energy using the above and equation 46. We further reintroduce the dependence on $k$,

[TABLE]

We recall that $\nabla\hat{C}_{k}$ can be computed from equation 47. Thus we have everything we need to simulate Hamiltonian dynamics to explore our target distribution. In practice, we make a few standard adaptations to the above algorithm as detailed in (Neal et al.,, 2011). We sample $\delta$ from a uniform distribution on $\mathcal{U}\left[a,b\right]$, as well as using a partial momentum refreshment with parameter $\alpha$. More specifically, given $\bm{p}$ from the previous iteration of the HMC algorithm and a sample $n\sim p_{0}(\bm{p})$ set $\bm{p}^{\prime}$ as

[TABLE]

5.7 Assessing Convergence

To assess convergence of our MCMC algorithms we visualise trace plots of our chains. In addition, we run two chains in parallel and we look at the potential scale reduction factors (the $\hat{R}$ statistic) and their upper $95\%$ confidence limits (Gelman and Rubin,, 1992, Brooks and Gelman,, 1998) using the coda R package (Plummer et al.,, 2006). We note that values of the $\hat{R}$ statistics far from $1$ indicate non-convergence. For the Drosophila example we run our a Hamiltonian-within-Gibbs sampler for $20,000$ iterations, performing a Hamiltonian Monte Carlo move to update the hyperparameters every $50$ iterations. We monitor the hyperparameters of the GPs carefully, since correlations can lead to slow exploration. As a representative example, we plot, in figure 8, the two parallel chains for the length-scale of GP associated to the ER. We see that mixing is rapidly achieved and that the upper $95\%$ confidence limit of $\hat{R}\approx 1$, indicating convergence. We contrast this with using a Metropolis-within-Gibbs sampler, in which we perform and Metropolis-Hastings move every $10$ iterations. We see from figure 9 the random walk nature of the parameters. In this case the upper $95\%$ confidence limit of $\hat{R}\approx 1.05$, thus our chain has most likely converged but exploration of the probability space is clearly slow. We also asses convergence using parallel chains in the set-up of section 3.2. We monitor $\epsilon$ the mixing weight of the outlier component, as an example, which can be seen in figure 10. The upper $95\%$ confidence limit of $\hat{R}\approx 1.01$, indicating convergence.

5.8 Tables of hyperparameters

Tables of hyperparameters and hyperparameter distributions for the mouse pluripotent stem cell data.

5.9 Prior Specifications

The priors for the comparison between classifiers in section 3.3 are as follows. The normal-inverse-Wishart prior for the multivariate Gaussian distributions was the following: the mean was set as the empirical mean of whole data, the shrinkage was set to $0.01$, the degrees of freedom was set to be the number of variables plus 2, the scale matrix was set to the identity matrix. For the GP we placed standard normal prior on each log hyperparameter. For all methods the Beta prior for the outlier component prior weight was set to be $\mathcal{B}(2,10)$ and the mixing proportions for each component was given symmetric Dirichlet prior with $\alpha=1$.