Inferring differential subcellular localisation in comparative spatial proteomics using BANDLE

The steady-state localisation of proteins provides vital insight into their function. These localisations are context specific with proteins translocating between different sub-cellular niches upon perturbation of the subcellular environment. Differential localisation, that is a change in the steady-state subcellular location of a protein, provides a step towards mechanistic insight of subcellular protein dynamics. Aberrant localisation has been implicated in a number of pathologies, thus differential localisation may help characterise disease states and facilitate rational drug discovery by suggesting novel targets. High-accuracy high-throughput mass spectrometry-based methods now exist to map the steady-state localisation and re-localisation of proteins. Here, we propose a principled Bayesian approach, BANDLE, that uses these data to compute the probability that a protein differentially localises upon cellular perturbation, as well quantifying the uncertainty in these estimates. Furthermore, BANDLE allows information to be shared across spatial proteomics datasets to improve statistical power. Extensive simulation studies demonstrate that BANDLE reduces the number of both type I and type II errors compared to existing approaches. Application of BANDLE to datasets studying EGF stimulation and AP-4 dependent localisation recovers well studied translocations, using only two-thirds of the provided data. Moreover, we potentially implicate TMEM199 with AP-4 dependent localisation. In an application to cytomegalovirus infection, we obtain novel insights into the rewiring of the host proteome. Integration of high-throughput transcriptomic and proteomic data, along with degradation assays, acetylation experiments and a cytomegalovirus intcractome allows us to provide the functional context of these data.


Introduction
The cell is compartmentalised into organelles and sub-cellular niches, allowing many biological processes to occur in synchrony (Gibson, 2009).Proteins are localised to these niches in accordance to their function and thus to shed light on the function of a protein, it is necessary to determine its subcellular location.A number of pathologies have implicated incorrect localisation as a contributing factor, including obesity (Siljee et lF, 2018), cancers (Kau et lF, 2004), neurological disorders (Davies et lF, 2018), as well as multiple others (Laurila and Vihinen, 2009).It is estimated that up to 50% of proteins reside in multiple locations (Christoforou et lF, 2016;Thul et lF, 2017), which complicates the study of their localisations.Community approaches have led to substantial improvements in our understanding of sub-cellular localisation (Thul et lF, 2017;Sullivan et lF, 2018).However, image-based approaches are often low in throughput and high-throughput alternatives are desirable.Furthermore, many biological processes are regulated by re-localisation of proteins, such as transcription factors shuttling from the cytoplasm to the nucleus, which are dicult to map using imaging methods at scale (Plotnikov et lF, 2011).
To simultaneously study the steady-state localisation and re-localisation of proteins, one approach is to couple gentle cell lysis and whole-cell fractionation with high-accuracy mass spectrometry (MS) (Christoforou et lF, 2016;Mulvey et lF, 2017;Geladaki et lF, 2019;Orre et lF, 2019).These approaches have already led to high-resolution subcellular maps of mouse embryonic stem cell (mESC) (Christoforou et lF, 2016), human cell lines (Geladaki et lF, 2019;Orre et lF, 2019), F erevisie (bakers' yeast) (Nightingale et lF, 2019), cyanobacterium (Baers et lF, 2019) and the apicomplexan oxoplsm qondii (Barylyuk et lF, 2020).Dynamic experiments have given us unprecedented insight into HCMV infection (Beltran et lF, 2016), EGF stimulation (Itzhak et lF, 2016), and EGFR inhibition (Orre et lF, 2019).In addition, CRISPR-Cas9 knockouts coupled with spatial proteomics has given insights into AP-4 vesicles (Davies et lF, 2018), as well as AP-5 cargo (Hirst et lF, 2018).In a study by Shin et lF (2020), the golgin long coiled-coil proteins that selectively capture vesicles destined for the Golgi were re-located to the mitochondria by replacing their Golgi targeting domains with a mitochondrial transmembrane domain (Shin et lF, 2020).This allowed the authors to readily observe the vesicle cargo and regulatory proteins that are redirected to the mitochondria, whilst avoiding technical issues that arise because of the redundancy of the golgins and their transient interaction with vesicles.Together, these collections of experiments suggest spatial proteomics can provide unprecedented insight into biological function.
Mass spectrometry-based spatial proteomics currently relies on supervised machine learning methods, such as support vetor mhines, to assign proteins to a sub-cellular niche using marker proteins with known localisations (Gatto et lF, 2010(Gatto et lF, , 2014a)).Advanced computational approaches have also been developed, including novelty detection algorithms (Breckels et lF, 2013;Crook et lF, 2020a) and transfer learning approaches (Breckels et lF, 2016b).These approaches are implemented in the pRoloc software suite (Gatto et lF, 2014b;Breckels et lF, 2016a), which builds on the MSnbase software (Gatto and Lilley, 2012) as part of the Bioconductor project (Gentleman et lF, 2004;Huber et lF, 2015).However, most machine learning methods fail to quantify uncertainty (estimate reliability) in the assignment of a protein to an organelle, which is paramount to obtaining a rich interrogation of the data.Crook et lF (2018) developed a Bayesian model to analyse spatial proteomics data and highlighted that uncertainty quantication can give insights into patterns of multi-localisation.The method is implemented as a tool as part of the Bioconductor project (Crook et lF, 2019a).
In dynamic and comparative experiments; that is, those where we expect re-localisation upon some stimulus to sub-cellular environment, the data analysis is more challenging.The task can no longer be phrased as a supervised learning problem, but the question under consideration is clear: which proteins have dierent sub-cellular niches after cellular perturbation?Procedures to answer this question have been presented by authors (Beltran et lF, 2016;Itzhak et lF, 2016Itzhak et lF, , 2017;;Kennedy et lF, 2020) and reviewed in Crook et lF (2020b).The approach of Itzhak et lF (2016,2017) relies on coupling a multivariate outlier test and a reproducibility score -termed the movement-reproducibility (MR) method.A threshold is then applied to these scores to obtain a list of proteins that re-locate; "moving" proteins.However, these scores can be challenging to interpret, since their ranges dier from one experiment to another and require additional replicates to calibrate the scores.Furthermore, the test ignores the spatial context of each protein, rendering the approach inecient with some applications allowing false discovery rates of up to 23% (Hirst et lF, 2018).Finally, the approach does not quantify uncertainty which is of clear importance when absolute purication of sub-cellular niches is impossible and multi-localising proteins are present.Recently, Kennedy et lF (2020) introduced a computational pipeline for analysing dynamic spatial proteomics experiments by reframing it as a classication task.However, this formulation ignores that some changes in localisation might be shifts in multi-localisation patterns or only partial changes , and it relies on the success of the classication.Furthermore, their approach cannot be applied to replicated experiments and so its applicability is limited.In addition, the authors found that they needed to combine several of the organelle classes together to obtain good results.Finally, the framing of the problem as a classication task only allows a descriptive analysis of the data.These considerations motivate the development of a more sophisticated and reasoned methodology.
We present Bayesian ANalysis of Dierential Localisation Experiments (BANDLE) -an integrative semi-supervised functional mixture model, to obtain the probability of a protein being differentially localised between two conditions.Posterior Bayesian computations are performed using Markov-chain Monte-Carlo and so uncertainty estimates are also available (Gilks et lF, 1995).We associate the term di'erentilly lolised to those proteins which are assigned dierent sub-cellular localisations between two conditions.Then, we refer precisely to this phenomenon as di'erentil lolistion, throughout the text.Hence, our main quantity of interest is the probability that a protein is dierentially localised between two conditions.BANDLE models the quantitative protein proles of each sub-cellular niche in each replicate of each experiment non-parametrically (Crook et lF, 2019b).A rst layer of integration combines replicate information in each experiment to obtain the localisation of proteins within a single experimental condition.To probabilistically integrate the two conditions, we use a joint prior probability distribution on the localisation of proteins in each condition, so that the datasets are modelled together.By examining the dierences in localisations, we can obtain a dierential localisation probability.Two prior distributions are proposed: one using a matrix extension of the Dirichlet Distribution and another, more exible prior, based on Pólya-Gamma augmentation (Polson et lF, 2013;Choi et lF, 2013;Linderman et lF, 2015).
A number of integrative mixture models for systems biology have been proposed including Multiple Dataset Integration (Kirk et lF, 2012), innite tensor factorisation approaches (Banerjee et lF, 2013), Bayesian Consensus Clustering (Lock and Dunson, 2013) and Clusternomics (Gabasova et lF, 2017).The methods have been developed mostly in the context of cancer sub-typing or transcriptional module discovery.Our approach is most similar to Clusternomics, which places a prior on the tensor (outer) product between the mixing proportions; but instead our model denes mixing proportions across datasets -upon which we introduce a prior.Importantly, our approach demands that there is an explicit link between components in each dataset, which can be dicult to assume outside the semi-supervised setting because of a statistical issue known as label-switching, allowing us to identify key model components (Richardson and Green, 1997).
To demonstrate the utility of BANDLE, we rst perform extensive simulations and compare to the MR approach.Our simulation study shows that our approach reduces the number of Type I and Type II errors, and, as a result, can report an increased number of dierentially localised proteins.These simulations also highlight the robustness of our approach to a number of experimental scenarios including batch eects.Our simulation studies also highlight that BANDLE provides interpretative improvements and clearer visualisations, and makes less restrictive statistical assumptions.We then apply our method to a number of datasets with well studied examples of dierential localisation, including EGF stimulation and AP-4 dependent localisation.We recover known biology and provided additional cases of dierential localisation, and suggest that TMEM199 localisation is potentially AP-4 dependent.Finally, we apply BANDLE to a human cytomegalovirus (HCMV) dataset -a case where MR approach is not applicable because the MR approach requires multiple replicates.Integration of high-throughput transcriptomic and proteomic data, along with degradation assays, acetylation experiments and a cytomegalovirus interactome allows us to provide the functional context of these data.In particular, we provide the spatial context of the HCMV interactome data.

The BANDLE workow
To perform statistical analysis of dierential subcellular localisation experiments we developed BAN-DLE.BANDLE is a semi-supervised integrative functional mixture model that allows the computation of a dierential localisation probability.The BANDLE workow, visualised in gure 1, begins with a well dened mass-spectrometry based spatial proteomics experiment.A cellular perturbation of interest is performed alongside control experiments in wild-type cells or another suitable control, depending on the application.To avoid confounding factors, control and treatments experiments should be performed in parallel with identical mass-spectrometry settings (Gatto et lF, 2010(Gatto et lF, , 2014a)).We further recommend checking that the experiment is successful by performing western blots on organelle markers prior to mass-spectrometry analysis, as well as examining the quality of clustering of marker proteins computationally (Gatto et lF, 2010(Gatto et lF, , 2014a(Gatto et lF, , 2019)).Our Bayesian approach BANDLE is applicable to experiments with any number of replicates, as well as several subcellular fractionation approaches or mass-spectrometry methods.BANDLE supports multiple mass-spectrometry-based methods, including label-free and labelled quantitation (e.g.SILAC and isobaric tags), and data-dependent and data-independent acquisition.BANDLE models the proles of the subcellular niches non-parametrically, using Gaussian random elds to capture spatial correlation.Niches are modelled separately across replicates and conditions to allow for variations.Each dataset is hence modelled as a mixture of the dierent subcellular niches.A distribution of 4 localisations is determined for each protein in each condition.Subsequently, this information is shared across the two conditions by using a joint prior on the protein localisation in each condition.A visualisation in the case of two replicates and two conditions is given in Figure 1C.By examining the inferred subcellular localisation across the two conditions and using Bayesian inference, we can compute a dierential localisation probability.To apply the Bayesian model, we rst calibrate the prior based on prior preditive heks (simulation of data from the prior probability distribution) (Gelman et lF, 1995).In all scenarios, we check the prior expectation of the number of dierentially localised proteins and the probability that more than l proteins are dierentially localised.These are reported in the supplement.We then proceed with Bayesian parameter inference using Markov-Chain Monte-Carlo (MCMC), which samples from posterior distribution of the model parameters, (Gilks et lF, 1995) and the checking of convergence.We visualise our results principally using rank plots, where proteins are ranked from those most likely to be dierentially localised or not (Figure 1D).

BANDLE reduces false positives and increases power
To assess the performance of BANDLE and the MR approach, we run a number of simulations allowing us to ascertain the dierence between each method in scenarios where we know the ground truth.Two variations of the MR approach were proposed and we compare to both (see methods section 4.1) We rst start with a real dataset from Drosophila embryos and simulate replicates, as well as 20 protein re-localisations (Tan et lF, 2009).To simulate these datasets a bootstrapping approach is used, coupled with additional noise eects.The rst simulation uses a simple bootstrapping approach, where a niche-specic noise component is included (see supplementary methods).The subsequent simulations start with the basic bootstrapping approach and add additional eects.The second and third simulations add batch eects: random and systematic respectively (see supplementary methods).Whilst the fourth simulation generates misaligned features across datasets by permuting them (fraction swapping) -this models misaligned fractions between replicates (see supplementary methods).The nal simulation includes both batch eects and feature permutations.The simulations are repeated 10 times, where each time we simulate entirely new datasets and relocalisations -this is repeated for each simulation task.We assess the methods on two metricsthe area under the curve (AUC) of the true positive rate and false positive rate for the detection of dierential localised proteins.Furthermore, we determine the number of correctly dierentially localised proteins at xed thresholds (see supplementary methods).
Our proposed method, BANDLE, signicantly outperforms the MR methods with respect to AUC in all scenarios (t-test p < 0.01, Figure 2).Furthermore, it demonstrates that BANDLE is robust to a variety of situations, including batch eects.We used two separate approach to incorporate prior information into BANDLE.The performance of BANDLE based on the Dirichlet prior is already very good and thus it is unsurprising that we do not observe any signicant improvements in AUC by including prior information on correlations captured by the Pólya-Gamma prior.Additional comparisons are made in the supplement where we make similar observations.The improved AUC, which demonstrates improved control of false positives and increased power, translates into increased discovery of dierentially localised proteins.Indeed, BANDLE with the Dirichlet prior discovers around twice as many such re-localising proteins than the MR approach.Allowing prior correlations through the Pólya-Gamma prior demonstrates that additional dierentially localised proteins are discovered.This is an important reality of those performing comparative and dynamic spatial proteomics experiments, since the experiments become more worthwhile with additional biological discoveries.In practice, the authors of the MR approach advocate additional replicates to calibrate which thresholds are used to declare a protein dierentially localised.This assumes that the perturbation of interest does not have a strong eect on the properties of the sub-cellular niches, which restricts applicability.In contrast, BANDLE does not need additional mass-spectrometry experiments to calibrate its probabilistic ranking meaning more discoveries are made at lower cost.
In the following section, we examine the dierences between the approaches in a simulated example.There we focus on the output, interpretation and statistical qualities of each approach, rather than the predictive performance of the methods.Figure 2: Boxplots comparing the performance of the MR approach (2016 and 2017) and our proposed method BANDLE.BANDLE is separated into whether a Dirichlet-based prior was used or if the Polya-Gamma augmentation was applied.Each boxplot corresponds to a dierent simulation scenario.The rst 5 boxplots show BANDLE has signicantly improved AUC in all scenarios.These AUCs are converted (posterior probability 0.999 or M = 2 and R = 0.9) to the correct number of re-localisations and we can see that our method clearly outperforms the MR approach.

BANDLE quanties uncertainty and is straightforward to interpret
In this section, we further explore the application of BANDLE with a Dirichlet prior and the MR approach, focusing on the interpretation and statistical properties of the two methods.Again, we simulate dynamic spatial proteomics data, starting from the Drosophila experiment in the scenario in which the MR method performed best.This is where there are cluster specic noise distributions but no other eects, such as batch eects, were included.Sample PCA plots of the data are presented in gure 3 A.There is a clear pattern of localisations across the data where proteins with known sub-cellular localisations are closer to each other.However, the organelle distributions clearly overlap and in some cases are highly dispersed -a representation of the challenges faced in real data.These data are annotated with 11 sub-cellular niches and 888 proteins are measured across 3 replicates of control and 3 of treatment (totalling 6 experiments).Re-localisations are simulated for 20 proteins.
We rst apply the MR method according to the methods in Itzhak et lF (2016,2017).We provide a brief description of the approach with full details in the methods (see section 4.1).To begin, the dierence proles are computed by subtracting the quantitative values for each treatment from each control.Then the squared Mahalanobis distance is computed to the centre of the data and under a Gaussian assumption the null hypothesis is that these distances follow a Chi-squared distribution, ergo a p-value is obtained.This process is repeated across the 3 replicates and the largest p-value was then cubed and then corrected from multiple hypothesis testing using the Benjamini-Hochberg procedure (Benjamini and Höchberg, 1995).A negative log 10 transform in then performed to obtain the M-score.To produce the R-score, Pearson correlations are computed between each dierence prole for all pairwise combination of dierence proles.The lowest of the three R-scores is reported.The M-score and R-score are plotted against each other (see gure 3 B) and the proteins with high M-score and high R-score are considers "hits".
There are a number of assumptions underlying the MR methodology.Firstly, comparing dierence proles pairwise assumes that the features in both datasets exactly correspond.However, this precludes any stimuli that changes the biochemical properties of the organelles, since changing these properties may result in diering buoyant densities, pelleting of niches at dierent centrifugation speeds or dierential detergent solubility.Thus, whether density-gradient, dierential centrifugation or dierential solubilisation is used for organelle separation this assumption must be carefully assessed.Secondly, the Gaussian assumption ignores the natural clustering structure of the data because of the dierent organelle properties.Indeed, examination of the p-value distributions in a histogram (supplementary gure 4) shows that it clearly deviates from the mixture of distributions expected (p-values are uniformly distributed under the null).The peaking of p-value towards 1 suggests poor distributional assumptions (Holmes and Huber, 2018).Thus perhaps the Chi-squared distribution is a poor t for the statistic of interested.Exploring this further, we t a Chi-squared and Gamma distribution empirically to the statistics using maximum likelihood estimation (MLE) of the parameters.Figure 3 C shows that the Gamma distribution is a better distributional tsuccessfully capturing the tail behaviour of the statistic (log Likelihood ratio: 1355 on 1 degree of freedom).The Chi-squared family is nested in the Gamma family of distributions, so if the theoretical Chi-squared distribution was a good t the distributions would overlap.For a quantitative assessment of model ts we compute the negative log-likelihood of the data given the optimal distributions -the Gamma distribution has a markedly lower negative log-likelihood (log Likelihood Ratio 1355).This provides strong evidence that the underlying Gaussian assumptions are likely violated.Thirdly, it is inappropriate to cube p-values: to combine p-value across experiments one could use Fisher's method (Mosteller and Fisher, 1948;Brown, 1975;Kost and McDermott, 2002) or the Harmonic mean p-value (HMP) (Good, 1958;Wilson, 2019) depending on the context.Indeed, the cube of the p-value is no longer a p-value.To elaborate, if P are a set of p-values, then under the assumption of the null hypothesis P is uniformly distributed; however, the cube is clearly not uniformly distributed.Since we no longer work with p-values, Benjamini-Hochberg correction becomes meaningless in this context.Transforming these values to a "Movement score", conates signicance with eect size which confounds data interpretation.Finally, summarising to a single pair of scores ignores their variability across experimental replicates.
BANDLE rst models each sub-cellular niche non-parametrically (since the underlying functional forms are unknown (Crook et lF, 2019b)).Visualisation of the posterior predictive distributions from these ts for selected sub-cellular niches is given in supplementary gure 3 -we observe a good correspondence between the model and the data.We can see that the dierent sub-cellular niches have contrasting correlation structures and thus niche specic distributions are required.These distributions are specic for each replicate of the experiment and also the two experimental conditions.The information from the replicates, and the control and treatment are combined using an integrative mixture model.Briey, mixing proportions are dened across datasets allowing information to be shared between the control and treatment (see methods for more details).This formulation allows us to compute the probability that a protein is assigned to a dierent sub-cellular niche between the two experiments -the di'erentil lolistion proility.The proteins can then be ranked from most probably dierentially localised to least (gure 3 D).The gure is simple to interpret: the proteins with highest rank are the most likely to have dierentially localised between the experiment, having been condently assigned to dierent sub-cellular niches in the control versus treatment.The proteins with lowest rank are highly unlikely to have moved during the experiment -the localisation are stable.This is important information in itself, especially when combined with other information such as changes in abundance or post translational modication.Figure 3 D (right) shows the 30 proteins with highest rank visualising the uncertainty in the dierential localisation probability (see methods).This ranking allows us to prioritise which proteins to follow up in validation experiments.The ranking can also be mapped onto other experimental data, such as expression or protein-protein interaction data.The probabilistic ranking produced by BANDLE is more closely aligned with the phenomenon of interest.Indeed, if we divide the data into the proteins that were dierentially localised and those that were not.Then from plotting the distribution of the statistics from the respective methods, it is clear that output from BANDLE is most closely associated with re-localisation events (gure 3 E).A Chi-square (orange) and Gamma (blue) t are overlaid (obtained using maximum likelihood estimation).The Gamma distribution clearly captures the tail behaviour.(D) A BANDLE rank plot where proteins are ranked from most to least likely to dierentially localised.The dierentially localisation probability is recorded on the y-axis.(right) A BANDLE rank plot of the top 30 dierentially localised proteins with uncertainty estimates for the dierential localisation probability.Proteins marked in orange were simulated translocations.(E) Violin plots for the dierential localisation probabilities (BANDLE), the M score (MR method) and R score (MR method).The distributiosn are split between dierential localised (movers) and spatially stable proteins.Clearly, the dierential localisation probabilities correlate most closely with the phenomena of interest.

Applications to dierential localisation experiments 2.3.1 Characterising dierential localisation upon EGF stimulation
Having carefully assessed the statistical properties of our approach, BANDLE, and the MR method, we apply these approaches to a number of datasets.First, we consider the hynmi yrgneller wps (DOMs) dataset of Itzhak et lF (2016), exploring the eects of EGF stimulation in HeLa cells.In this experiment, SILAC labelled HeLa cell were cultured and recombinant EGF was added to the culture at a concentration of 20 ng ml −1 (see (Itzhak et lF, 2016)).A total of 2237 complete protein proles were measured across 3 replicates of control and 3 replicates of EGF treated HeLa cells.Principal Component Analysis (PCA) projections of the data can be visualised in the supplement.A quality control assessment was performed using the approach of Gatto et lF (2019).As a result, Nuclear pore complex, peroxisome and Golgi annotations were removed, since the marker proteins of these classes were highly dispersed.
The MR method was applied as described in the methods and the results can be visualised in gure 4 A. 7 proteins are predicted to be dierentially localised using the MR method with the thresholds suggested by Itzhak et lF (2016).These include 3 core proteins of the EGF signalling pathway SHC1, GRB2 and EGFR (Oda et lF, 2005) and other, potentially related, proteins TMEM214, ACOT2, AHNAK, PKN2.Since the MR approach does not provide information about how the functional residency of the proteins change, it is challenging to interpret these results without further analytical approaches.
To quantify uncertainty and gain deeper insight into the perturbation of HeLa cell after EGF stimulation we applied our BANDLE pipeline.Firstly, the rank plots display a characteristic shape suggesting that most proteins are unlikely to be dierentially localised upon EGF stimulation (gure 4 B).Furthermore, we provide uncertainty estimates in the probability that a protein is dierentially localised for selected top proteins (gure 4 C).Furthermore, we visualise the change in localisation for the proteins known to re-localise upon EGF stimulation: SHC1, GRB2 and EGFR (gure 4 E).This is displayed by projecting the posterior localisation probabilities on to the corresponding PCA coordinates.These probabilities are then smoothed using a Nadaraya-Watson kernel estimator (Nadaraya, 1964;Watson, 1964) and visualised as contours.PCA plots of the raw data are found in the supplementary materials.
Given the well-documented interplay between phosphorylation and sub-cellular localisation (Lee et lF, 2012; Christian et lF, 2016; Puertollano et lF, 2018; Balta et lF, 2018), we hypothesised that proteins with the greatest dierential phosphorylation would correlate with proteins that were more likely to be dierentially localised.To this end, we integrated our analysis with a time-resolved phosphoproteomic dataset of EGF stimulation using MS-based quantitation (Köksal et lF, 2018).In their study, cells were harvested at eight dierent time points after EGF stimulation: 0,2,4,8,16,32,64 and 128 minutes.Cells were harvested and protein digested to peptides using trypsin.Peptides corresponding to each time point were labelled with a dierent iTRAQ tag before combining all samples together and quantifying using LC-MS/MS.Immunoprecipitation was used to enrich for phosphorylated tyrosine residues (Possemato et lF, 2017) and the enrichment of phosphosites on serine and threonine residues was performed via immobilized metal anity chromatography (IMAC) (Ficarro et lF, 2002;Moser and White, 2006).
For each phosphopeptide corresponding to a unique protein, we computed the largest log 2 fold 13 change observed across the time course.Given that the changes in localisation occur within 20 minutes, we restricted ourselves to the rst 6 time points (Itzhak et lF, 2016).We then took the top 10 proteins ranked by each of the MR method and BANDLE.These rankings are then correlated with rankings obtained from the changes in phosphorylation.The Spearman rank correlations were recomputed for 5, 000 bootstrap resamples to obtain bootstrap distributions of correlations (see gure 4).We report the mean correlation and the 95% boostrap condence intervals.The correlation between the ranks of the MR method and the phosphoproteomic dataset was ρ S = 0.40 (−0.49, 0.85), whilst the correlation when using the ranking of BANDLE was ρ S = 0.68 (0.02, 0.98).That is, phosphorylated proteins are more likely to be dierential localised and this signal is more clear using the ranking obtained from using BANDLE.Alongside the statistical and interpretable benets of BANDLE, it is clear the approach has the utility to provide insight into localisation dynamics. 14

BANDLE obtains deeper insights into AP-4 dependent localisation
The adaptor protein (AP) complexes are a set of heterotetrameric complexes, which transport transmembrane cargo protein vesicles (Robinson, 2015).The AP1-3 complexes are well characterised: AP-1 mediates the transport of lysosomal hydrolases from the trans-Golgi to the endsomes (Karin et lF, 1997;Hess et lF, 2004); AP-2 has a signicant role in the regulation of endocytosis (Motley et lF, 2003); AP-3 is involved in the sorting of trans-Golgi proteins targeted to the lysosome (Dell'Angelica et lF, 1998).The role of the AP-4 complex has become better understood in recent years (Hirst et lF, 1999(Hirst et lF, , 2013 AP-4 consists of four subunits (β4, ε, µ4 and σ4) forming an obligate complex (Dell'Angelica et lF, 1998).Davies et lF ( 2018) study the functional role of AP-4 using spatial proteomics; in particular, the DOM workow mentioned previously.As part of their study, they use AP-4 CRISPR knockout cells to interrogate the eect on the spatial proteome when AP-4 function has been ablated.
Re-analysis of this subcellular proteomics experiment provides full quantitative measurements for 3926 proteins across two replicates of wild-type cells and two replicates where the β4 subunit has been knocked-out.They also produce two replicates where AP4E1 is knocked-out but this is not considered here for brevity.The data are visualised as PCA plots (see supplement).As in the previous analysis, we run a quality control step removing the actin binding protein and nuclear pore complex annotations (Gatto et lF, 2019).This dataset is particularly challenging to analyse because there are only two replicates for each condition.The value of Bayesian analysis is the ability to provide prior information to regularise, as well as the quantication of uncertainty, which is more critical in data sparse scenarios.
Previous application of the MR methods led to authors to nd that SERINC 1 (Q9NRX5), SERINC 3 (Q13530) were dierentially localised, as well as an altered subcellular distribution for ATG9A (Q7Z3C6) (Davies et lF, 2018).Their results suggest these are cargo proteins of the AP-4 complex that are packaged into vesicles at the trans-Golgi before being transported to the cell periphery.All together their results suggest AP-4 provides spatial regulation of autophagy and that AP-4 neurological pathology is linked to disturbances in membrane tracking in neurons (Mattera et lF, 2017;Davies et lF, 2018).
We apply our method BANDLE to gain further insights into AP-4 dependent localisation.We compute the dierential localisation probability; the associated uncertainty estimates and rank proteins according to this statistic (see gures 5 A and B).Characteristic S shaped plots are observed with most proteins not dierentially localised upon knock-out of AP-4 β4.The results of both SER-INC 1 and 3 are validated, as we compute a dierential localisation probability greater than 0.95 for these proteins.Furthermore, 16 of the top 20 proteins are membrane-bound or membrane-associated proteins (FDR < 0.01 hyper-geometric test).To demonstrate the benet of our probabilistic ranking, we perform two-sided KS rank test against the functional annotations provided in the STRING database (corrected for multiple testing within each functional framework) (Szklarczyk et lF, 2019).
We nd that processes such as ER to Golgi transport and lipid metabolism are more highly ranked than would be expected at random (FDR < 0.01), as well as endosomes and Golgi localisations (FDR < 0.01).Whilst processes associated with translation, ribosome localisation and function appear signicantly lower in the ranking (FDR < 0.01).As expected, this provides a high level overview and evidence for the functional nature of AP-4 in the secretary pathway.
Taking a more precise view on our results, we examine the top 20 dierentially localised proteins in more detail.We compute the Spearman correlation matrix between these proteins and observe strong correlation, suggesting the proteins act in a coordinated way (see gure 5 C).Visualising the data in a heatmap (gure 5 D), after mean and variance normalisation, we observe a highly concordant pattern: most proteins are enriched in fractions 4 and 5.These fractions are obtained from the highest centrifugation speeds and so dierentially pellet light membrane organelles, such as endosomes and lysosomes (Itzhak et lF, 2016;Geladaki et lF, 2019).Again, further evidence for the role of AP-4 dependent localisation dynamics within the secretary pathway.
In gure 5 C, we observed a large cluster of 9 proteins, which included SERINC 1 and 3. Amongst these 9 proteins is SLC38A2, a ubiquitously expressed amino-acid transporter that is widely expressed in the central nervous system and is recruited to the plasma membrane from a pool localised in the trans-Golgi (Hatanaka et lF, 2000;Bevilacqua et lF, 2005;Gonzalez-Gonzalez et lF, 2005;Melone et lF, 2006).Thus, its suggested dierential localisation here provides further evidence for the role of AP-4 as a membrane tracker from the trans-Golgi.Another protein in this cluster is TMEM 199 (Q8N511) a protein of unknown function that is involved in lysosomal degradation (Miles et lF, 2017).Furthermore, it has been implicated in Golgi homoeostasis but the functional nature of this process is unknown (Jansen et lF, 2016).Probing further, we observe that TMEM199 acts in a coordinated fashion with SERINC 1 and 3. Marked re-localisations are observed on PCA plots toward the endo/lysosomal regions (see gure 5 E) and we note that the quantitative proles of SERINC 1, SERINC 3 and TMEM199 act in an analogous way upon AP-4 knockout (see gure 5 F).Our ndings motivate additional studies to elucidate AP-4 dependent localisation and separate these observations from potential clonal artefacts.

Rewiring the proteome under Cytomegalovirus infection 2.4.1 The host spatial-temporal proteome
Human Cytomegalovirus (HCMV) infection is a ubiquitous herpesvirus that burdens the majority of the population (Cannon et lF, 2010).In healthy immune systems, HCMV establishes latent infection following initial viral communication (Reeves et lF, 2005) and reactivation can lead to serious pathology in some imunno-compromised individuals (Boeckh and Nichols, 2004).HCMV has the largest genome of any known human virus, at 236 kbp it encodes for over 170 proteins that modulate almost all aspects of the hosts cellular environment for its benet (Murphy et lF, 2003; Stern-Ginossar et lF, 2012; Jean Beltran and Cristea, 2014).Initial viral infection involves endocytosis of the virion into the cell (Isaacson et lF, 2008), host machinery is then used to transport viral capsids into the nucleus (Ogawa-Goto et lF, 2003).Within the host nucleus viral transcription and genome replication occurs (Milbradt et lF, 2007;Gibson, 2008;Kalejta, 2008).Meanwhile, other viral proteins are targeted to the secretory pathway to inhibit the host immune response and regulate the expression of viral genes (Stamminger et lF, 2002;Feng et lF, 2006;Hwang and Kalejta, 2007;Mitchell et lF, 2009;Cristea et lF, 2010;Li et lF, 2013), rewire signalling pathways (Yurochko, 2008) and modulate metabolism (Yu et lF, 2011).In later phases, the cellular tracking pathways and the secretory organelles are hijacked for the formation of the viral assembly complex (vAC) (Buchkovich et lF, 2010;Moorman et lF, 2010;Alwine, 2012;Das et lF, 2007;Das and Pellett, 2011).Due to the diversity of cellular processes manipulated during HCMV infection, it is often used as a paradigm to analyse virus-host interactions (Weekes et lF, 2014).
There has been a recent urry in applying system-wide proteomic approaches to the HCMV infection model.Weekes et lF (2014) developed quntittive temporl viromis a multiplexed proteomic approach to understand the temporal response of thousands of cellular host and viral proteins.More recently, to discover proteins involved in the innate immune response, a multiplexed proteasomelysosome degradation assay found that more than 100 host proteins are degraded shortly after viral-infection (Nightingale et lF, 2018).Meanwhile, a comprehensive mass spectrometry interactome analysis has identied thousands of host-virus interactions (Nobre et lF, 2019).Furthermore, high-throughput temporal proteomic analysis has revealed the importance of protein acetylation (post-translational modication of lysine amino acids), as an integral component during HCMV infection (Murray et lF, 2018).
Beltran et lF ( 2016) use spatial and temporal proteomics to investigate the response of the human host proteome to HCMV infection.The authors perform subcellular fractionation on uninfected (control) and HCMV infected (treated) cells at 5 dierent time point (24,48,72,96,120) hours post infection (hpi).The authors then used neural networks to classify proteins to sub-cellular niches at each time point in the control and treated cells, allowing a descriptive initial analysis of the data.Proteins with dierential classication at each time point are those that are believed to be dierential localised.However, the challenge of this study is that only a single replicate is produced in each condition.This renders the MR method of (Itzhak et lF, 2016) inapplicable.
Dierential classication is a reasonable approach to probe dierential localisation though it neglects information shared across both experiments and it is not quantitative (i.e.no p-value or posterior probability of change).In the case of single replicates, by sharing information and providing prior information we are able to improve inference and obtain deeper insights.We apply BANDLE to control and HCMV-treated cells at 24 hpi, in the interest of brevity, to explore further the host spatial-temporal proteome.Our analysis reects a potential rewiring of the proteome with many possible proteins dierentially localised on HCMV infection.We highlight an example of dierential localisation with SCARB1 (see gure 6 A), with a localisation in the secretory pathway shifting toward a PM/cytosolic localisation, similar to what has previously been observed (Beltran et lF, 2016).
To obtain global insights into the functional behaviour of the dierentially localised proteins, we performed a Gene Ontology (GO) enrichment analysis.An extensive list of terms is enriched and these can be divided broadly into subcategories such as translation and transcription; transport; viral processes; and immune process (see supplementary materials).These results reect closely the early phase of HCMV infection (Jean Beltran and Cristea, 2014).Pathway enrichment analysis highlights terms related to a viral infection (Viral mRNA Translation, Inuenza Life Cycle, Infectious disease, Innate Immune System, Immune System, MHC class II antigen presentation, Antigen processing-Cross presentation, Host Interactions of HIV factors, HIV Infection) (see supplementary materials).Pathway analysis also reveals known processes that are modulated during HCMV infection, such as membrane tracking (Bozidis et lF, 2010;Niemann et lF, 2014;Zeltzer et lF, 2018), Extracellular matrix organization (Reinhardt et lF, 2006) and rab regulation of tracking (Lu£in et lF, 2018).

Integrating HCMV proteomic datasets to add functional relevance to spatial proteomics data
The spatial information obtained here allows us to perform careful integration with other highresolution proteomic datasets.The degradation screens by Nightingale et lF ( 2018) identied proteins that were actively degraded during HCMV infection but gave no information regarding the spatial location of the targets.To determine the location of host proteins targeted by HCMV for degradation, the BANDLE revised spatial data at 24 hpi was overlapped with proteins that were degraded by the proteasome or lysosome.The subcellular location of the host proteins is displayed for the 24 h timepoint.To determine the spatial granularity of the degradation data we tested whether the proteins assigned to each spatial pattern had a signicantly dierent degradation distribution than the distribution of all proteins in the experiment (t-test).We note that proteins that are dierentially localised are no more likely to be targeted for degradation and those that are not (see supplement).Analysis of changes in protein abundance can be used to generate turnover rates in both HCMVinfected and Mock-infected cells.Comparing these turnover rates allows us to calculate the rescue ratio, which identies proteins that exhibit increased degradation during viral infection compared to baseline.Specically the rescue ratio is obtained by comparing abundance during HCMV infection inhibitor with protein abundance during mock infection ∓ inhibitor.Degradation data from Nightingale et lF ( 2018) are overlaid as a heatmap, showing a − log 10 (p-value) for each inhibitor (see supplementary gures 35 and 37).For proteasomal targeted proteins (MG132 inhibitor), the data highlight a high number of proteins degraded from the mitochondria.The mitochondria act as a signalling platform for apoptosis and innate immunity and it is already well-established that HCMV can subvert these processes to its advantage (Crow et lF, 2016).Furthermore, there is a high degree of protein degradation as one might expect in proteasome fractions (dense cytosol), with an enrichment of proteins recruited from the ER and cytosol (see supplementary gure 35).For lysosomal targeted proteins (leupeptin) there was a high degree of proteins degraded from the mitochondria, cytosol and plasma membrane.There were also several proteins degraded that moved from the cytosol to the dense cytosol (see supplementary material section 16).
Many host proteins are up-or-down regulated upon HCMV infection (Weekes et lF, 2014).We examine more recent abundance data from Murray et lF (2018) at 24 hpi and rst we note that dierentially localised proteins are not more abundant than spatially stable proteins (see supplementary gure 41).However, we see a strong spatial pattern when we overlay the abundance pattern on a heatmap.In supplementary gure 34, we report the mean log2 fold change for proteins stratied according to predicted subcellular localisation.It is important to combine spatial and abundance data, since a dierentially localised protein may not undergo a true translocation event but rather a new pool of proteins is synthesised.Which of these options is true could be further investigated by coupling spatial proteomics workows with time-resolved incorporation of stable isotope labelled amino acids or azido-homoalanine to interrogate the location of newly synthesized proteins.The signicance of these abundance changes is highlighted in supplementary gure 33.For example, there is a signicant decrease in the abundance of the protein recruited to the dense cytosol from the ER (see supplementary gure 42).Some of the larger changes are not signicant because there are too few proteins with the same spatial pattern.We note that FAM3C, a protein involved in platelet degranulation, is upregulated at 24 hpi.Furthermore, FAM3C relocalises from the Golgi to the lysosome, its Golgi localisation is in concordance with the Human Protein Atlas (HPA) (Thul et lF, 2017) and its lysosome relocalisation suggests that it is tracked through the secretory pathway before undergoing degranulation.
Upon integration of the acetylation data of Murray et lF (2018), the spatial patterns are much more nuanced (see gures 6 E and supplementary material section 17).Perhaps surprisingly, we do not observe increased acetylation levels amongst dierentially localised proteins (see supplementary gure 44).The only signicant pattern is for proteins relocalising from the dense cytosol to the cytosol; however, we observe this is driven by a single protein Skp1 (see supplementary gure 45), which shows a 2.5-fold increase in acetylation at 24 hpi for Skp1 and there is an increase in its RNA transcript at 24 hpi (Nightingale et lF, 2018).The Skp1 protein is part of an E3 ubiquitin ligase complex that targets proteins for degradation.E3 ligases are often manipulated by viruses in order to control cellular processes to create a cell states that benet viral replication and survival (Mahon et lF, 2014).It is therefore possible that HCMV is controlling Skp1 activity through acetylation at its C-terminus, leading to its translocation and likely change in function.
The recent publication of the HCMV interactome has provided a wealth of data that gives insights into the function of the 170 canonical and two non-canonical viral protein-coding genes (Nobre et lF, 2019).However, a common diculty with analysing large interactome projects is the ability to reduce the number of false-positive interactions, leading to poor agreement between experimental and computational datasets.This can be controlled through replicates, supervised machine learning and increased statistical stringency; however, background contamination can never be eliminated.If a protein is located in a single location, you would expect true positive interactors to be located in the same subcellular compartment.Therefore, to narrow the list of viral-protein interactors, we overlapped spatial information from Beltran et lF (2016) with the viral interactors from (Nobre et lF, 2019) (Figure 6 C and D, supplementary gures 48 and 49).This provides a far more stringent set of high condence protein-protein interactions.Although, this comes the cost of removing some interactions between proteins that are located in more than one subcellular location and are therefore absent from the spatial dataset from Beltran et lF (2016).
We plot heatmaps to indicate the spatial distribution of the host proteins (gure 6 C and D).The overall distribution is plotted in the heatmap of gure 6 B. Firstly, we are interested in scenarios where the interacting host proteins were more likely to retain their localisation upon HCMV infection (than the computed posterior distribution would have predicted).Thus, for each viral bait, we simulated from a binomial A ∼ Bin(n, p) where p is the posterior probability that a random protein was assigned to the same localisation and n is the number of interactors of that viral bait.We then simulated from this distribution 5, 000 times to obtain a histogram (see supplementary gure 47).Viral baits of interest are those were the observed statistic in the tails of these histograms.
Examples of such cases are shown for viral proteins UL8 and UL70 (see gure 6 D and supplementary gure 48).The majority of UL8 interactors were located in the plasma membrane and cytosol.UL8 is a transmembrane protein that is transiently localised at the cell surface, with a small cytoplasmic pool (Pérez-Carmona et lF, 2018), perfectly mimicking the location of the majority of UL8 interactors.Practically all UL70 interactors were located in the cytosol.Viral UL70 is a primase known to locate to both the nucleus and cytoplasmic compartments during HCMV infection (Shen et lF, 2011).As the nucleus was removed prior to fractionation then one expects only to be able to interrogate cytosolic interactors.An example where the host proteins were spatially diuse was UL148A an elusive viral protein of unknown function, believed to be involved with modulating the innate immune response (Dassa et lF, 2018).UL148A appears to interact with host proteins distributed throughout the cell suggesting it is highly promiscuous (gure 6 C).Perhaps UL148A is a protein with multiple possible functions (Jeery, 2009) making its function hard to pinpoint and such an observation would not be uncommon for viral proteins because of limited genomic size (Cook and Lee, 2013;Copley, 2014).These results illustrate the strength in overlapping spatial proteomics with interactome studies to decrease the number of false positives and focus research on higher condence protein-protein interactions.The entire list of spatially resolved viral protein interactions is shown in the supplementary material.

3 Discussion
We have presented a Bayesian model for comparative and dynamic spatial proteomic experiments.Unlike current approaches, our exible integrative mixture model allows any number of replicate experiments to be included.Furthermore, subcellular proles are modelled separately for each condition and each replicate, allowing cases where the correlation proles dier between experiments.Crucially, our model facilitates the computation of dierential localisation probability, which cannot be performed by other methods in the literature.Furthermore, BANDLE probabilistically assigns proteins to organelles and can model outliers meaning that further supervised machine learning after application of BANDLE is not required.The probabilistic ranking obtained from BANDLE can be used for downstream pathway or GO enrichment analysis, likewise it can be mapped onto other orthogonal high-throughput datasets.
We compared BANDLE to the MR approach of Itzhak et lF (2016,2017).The MR method is not as broadly applicable as BANDLE, and BANDLE does not require additional experiments to interpret the thresholds.In our careful simulation study, we demonstrate reduced Type 1 error and increased power when using our approach.In a further simulation, we demonstrated that BANDLE has more desirable statistical properties than the MR approach, the results are easier to interpret and more information is available.Since we are in a Bayesian framework, our approach also quanties uncertainty allowing us to obtain .
Application of our approach to three dynamic and comparative mass-spectrometry based spatial proteomic experiments demonstrates the broad applicability of our approach.We validated many previously known ndings in the literature, placing condence in these results.When BANDLE was applied to EGF stimulation dataset, we saw increased correlation between our dierential localisation results and a phosphoproteomic timecourse than when compared to the results of the MR approach.
We applied BANDLE to an AP-4 knockout dataset to investigate AP-4 dependant localisation and, as with other studies, we observe SERINC 1 and SERINC 3 are examples AP-4 Cargo.Furthermore, we implicate TMEM199 as potentially overlooked AP-4 cargo, though it remains to rule it out as a potential clonal artefact.We apply BANDLE to a datasets where the MR approach is not applicable -an HCMV infection spatial proteomic dataset.Pathway and GO enrichment results implicate dierentially localised protein in well-studied processes of early viral infection; such as, membrane tracking and immune response.
We then carefully integrated several HCMV proteomic datasets and place a spatial perspective on these data, including proteins targeted for degradation, as well as abundance and acetylation dataset.In addition, we augment a recent HCMV interactome by placing it in its spatial context and note that most host protein interactomes are in the same localisation as their viral bait.This provides an excellent resource for the community and highlights the benet of integrating spatial proteomics and interactomics datasets.
Our analysis here highlights the potential role for post-translational modications (PTMs) and their inuence on localisation.The current datasets are limited because the spatial information is averaged over dierent PTMs.Thus, it is vital to develop methods to obtain spatial PTM information and develop corresponding computational tools to analyse these data.Furthermore, our approach here can only look at a single condition at a time.In the future, more complex spatial proteomics designs will be available that will study multiple perturbations simultaneously.
Overall, dierential localisation experiments seek to add an orthogonal perspective to other assays, such as classical high-throughput dierential abundance testing.Currently, dierential localisation has not been extensively explored in high-throughput.We hope rigorous statistical methods will spur extensive and illuminating applications.An R-package is provided for analysis at https://github.com/ococrook/bandle,building on a suite of packages for analysing spatial proteomics data (Gatto and Lilley, 2012;Gatto et lF, 2014b;Crook et lF, 2019a).

The Movement-Reproducibility method
The movement-reproducibility (MR) method was proposed by Itzhak et lF (2016,2017) and this is our interpretation of their method.We suppose that we are given two spatial proteomics experiments under a single contrast/perturbation/treatment, and denote unperturbed by (d = 1) and (d = 2) for the perturbed condition.Furthermore, assume we measure each condition with r = 1, ..., R biological replicates.Let 1 ] denote the concatenation of replicates for condition 1 and, likewise, for condition 2 we denote where ∆ = [∆ (1) , ..., ∆ (R) ].This assumes that both features and replicates are comparable in some way; that is, a feature in the r th replicate is directly comparable to the same feature in another replicate.Then, for each ∆ r , r = 1, .., R, the squared Mahalanobis distance D M from each protein to the empirical mean is computed using a robust estimate of the covariance matrix -the minimum covariance determination method (Hubert and Debruyne, 2010).Under a Gaussian assumption on ∆ r , D M (p i ) follows a chi-squared distribution with degrees of freedom equal to the dimension of the data G.Then, for each protein and each replicate a p-value is computed, such that there are R such p-values for each protein.These p-values are combined into a score by taking the cube of the largest p-value for each protein, correcting for multiple hypothesis testing using the Benjamini-Hochberg procedure and computing the − log 10 of the resultant value.For Itzhak et lF ( 2016), the p-value is not cubed and simply the largest p-value is taken.The nal score is called the M score.This process means that the computed value can no longer be interpreted as truly derived from a p-value.To maintain this interpretation, one could instead combine p-values using Fisher's method (Mosteller and Fisher, 1948).Furthermore, the authors are, implicitly, concerned with nding ny false positives and as such control over the FWER is desired rather than the FDR.Since FWER ≥ FDR, control of the FDR does not lead to control over the FWER.
A so-called reproducibility (R) score is obtained by rst computing the pearson correlation pairwise between matrices ∆ i , ∆ j , i = j for each protein.A nal R score, for each protein, is obtain by taking the minimum value for each protein.Again this score could have be interpret in a formal testing procedure using a permutation test (Efron, 2012) and furthermore includes an assumption of bivariate normality.Moreover, Pearson's correlation is unresponsive to many non-linear relationships which might be present.
Finally, each protein has an associated pair of scores, referred to as the MR-score.To determine thresholds for these scores the authors take a desired FDR = 0.01.Thus they repeat a control experiment 6 times to determine thresholds M = 2, R = 0.9: a region with no false discoveries.
Repeating the control experiment 6 times is a costly process and likely to be prohibitive for most experiments, particularly for cells that are expensive to culture.Furthermore, since the thresholds are empirically derived, this process needs to be repeated for every new experiment to determine optimal thresholds.

A model for dierential localisation
In the following, we layout our model for BANDLE, along with methods for inference, and approaches for summarising and visualising the output.Firstly, suppose we have two spatial proteomics experiments with unperturbed (d = 1) and perturbed conditions (d = 2).Furthermore, assume we measure each condition with r = 1, ..., R biological replicates.Let 1 ] denote the concatenation of replicates for condition 1 and likewise for condition 2 denotes We introduce the following latent allocation variable z i,d , representing the localisation of protein i in condition d.Thus, if z i,d = k this means that protein i localises to organelle k in dataset d.Given this latent allocation variable, we assume that the data from replicate r = 1, ..., R arises from some component density F (•|θ (r) k ).Hence, letting θ be the set of all component parameters, we can write We assume that biological replicates are independent and so we factorise as follows To couple the two conditions together we assume a joint prior structure for the latent allocation variable in each dataset.To be more precise, we construct a prior for the pair (z i,1 , z i,2 ).We x the possible number of subcellular niches to which a protein may localise to be K.Now, we introduce the matrix Dirichlet distribution, which we denote as MDir(α, K).The concentration parameter α is a K × K matrix, such that for a matrix π, the pdf of the matrix Dirichlet distribution is where B denotes the beta function, α k denotes the k th row of α, and j,k π jk = 1.Thus, we propose the following hierarchical structure π|α ∼ MDir(α, K) (5) where (z i,1 , z i,2 ) ∼ cat(π) means that the prior allocation probabilities are given by The above model is conjugate, and so if where γ j,k = α jk + n j,k .The likelihood models for the data are Gaussian Random Fields, which we elaborate on in the following section.Hence, the conditional posterior of the allocation probabilities are given by

Likelihood Model
The model described in the previous section is presented in a general form, so it could be applied to many dierent modes of data.We describe the model for a single spatial proteomics experiment, since the same model is assumed across all spatial proteomics experiments, that are then subsequently joined together using the approach in the previous section.Though the model is the same across experiments, the parameters are experiment-specic.
We assume that the protein intensity x i at each fraction s j can be described by some regression model with unknown regression function: where µ i is some unknown deterministic function of space and ε ij is a noise variable, which we assume is ε ij ∼ N (0, σ 2 i ).Proteins are grouped together according to their subcellular localisation; such that, all proteins associated to subcellular niche k = 1, ..., K share the same regression model.Hence, we write µ i = µ k and σ i = σ k .Throughout, for clarity, we refer to sub-cellular structures, whether they are organelles, vesicles or large protein complexes, as omponents.The regression functions µ k are unknown and thus we place priors over these functions to represent our prior uncertainty.Protein intensities are spatially correlated and thus we place qussin ndom pield (GRF) priors over these regression functions.We pedantically refer to these as GRF priors rather than qussin roess priors to make the distinction between the 1D spatial process that separates sub-cellular niches and the experimental cellular perturbations, which are potentially temporal in nature.Hence, we write the following which is dened as: Thus, each component is captured by a Gaussian Random Field model and the full complement of proteins as a nite mixture of GRF models.The protein intensity for each experiment might be measured in replicates.For a suciently exible model, we allow dierent regression models across dierent replicates.To be more precise, consider the protein intensity x (r) i for the i th protein measured in replicate r at fraction s (r) j , then we can write the following having assumed that the i th protein is associated to the k th component.The (hyper)parameters for the Gaussian Random Field priors for the r th replicate in experiment d are denoted by θ (r) k,d .We denote by θ the collection of all hyperparameters and the collection of priors for these hyperparameters by G 0 (θ).The loss of conjugacy between the prior on the hyperparameters and likelihood is unavoidable.
The GRF is used to model the uncertainty in the underlying regression functions; however, we have yet to consider the uncertainty that a protein belongs to each of these components.To capture these uncertainties, we can use the model in the previous section, allowing information to be shared across each condition.Following from the previous section, the conditional posterior of the allocation probabilities is where, in the specic case of our likelihood model the probabilities in the terms of the product can be computed using the appropriate GRF.We assume that our GRFs are centred and that the covariance is from the Matern class (Stein, 1999).The Matern covariance is specied as follows where Γ is the gamma function and K v denotes the modied Bessel function of the second kind of order ν > 0. Furthermore, a and ρ are positive parameters of the covariance.a 2 is interpreted as a marginal variance, whilst the non-standard choice of √ 8ν in the denition of the Matern covariance, allows us to interpret ρ as a range parameter and thus ρ is the distance at which the correlation is 0.1 for any ν (Lindgren et lF, 2011).The Matern covariance arises from solutions of the following linear fractional stochastic partial dierential equation (SPDE): where W(u) is spatial Gaussian white noise with unit variance and ∆ is the Laplacian.The parameter ν controls the dierentiability of the resulting sample paths; such that, v is the number of mean-square derivatives.For typical applications, ν is poorly identiable and xed; ν = 1/2 recovers the exponential covariance, whereas taking the limit ν → ∞ one obtains the squared exponential (Gaussian) covariance.We x ν = 2.
A ridge in the marginal likelihood for the marginal variance and range parameters of the Matern covariance makes inference challenging.Indeed, dierent hyperparameters lead to unconditional prior simulations with the same spatial pattern but dierent scales (Rasmussen and Williams, 2006;Fuglstad et lF, 2019).Furthermore, when the intrinsic dimension of the Gaussian random eld is less than four, there is no consistent estimator under in-ll asymptotics for ρ and a.A principled prior, which allows domain expertise to be expressed, is thus desired to enable stable inferences.A number of works considered reference priors for GRFs (Berger et lF, 2001;Paulo et lF, 2005;De Oliveira, 2007;van der Vaart et lF, 2009).Here, we employ a recently introduced collection of weakly-informative priors, which we introduce in the supplementary methods.

Calibration of Dirichlet prior
The following section describes how to calibrate the prior based on expert information and prior predictive checks.Recall the prior on the allocation probabilities is the following p(z i,1 = k, z i,2 = k |π) = π kk . (16) The matrix π has π jk has its (j, k) th entry and π jk is the prior probability that a protein belongs to organelle j in dataset 1 (control) and k in dataset 2 (contrast).The diagonal terms represent the probability that the protein was allocated to the same organelle in each dataset.The non-diagonal terms are the prior probability that the protein was not allocated to the same organelle.Since the number of non-diagonal terms greatly exceeds to the number of diagonal entries, it is important to specify this prior carefully.Recall that the prior is given a matrix Dirichlet distribution with concentration parameter α.Firstly, we are interested in the prior expectation of the number of proteins that are di'erentil lolised ; that is, proteins not allocated to the same organelle in both conditions.Let ρ be the prior probability that a protein is not allocated to the same organelle.Then it follows that p(z i,1 = z i,2 |π) =: ρ = j,k;j =k π jk . (17) By properties of the Dirichlet distribution we have that the the marginal distribution of π jk is given by π jk ∼ B(α jk , α 0 − α jk ), where α 0 = j,k α jk .Thus, the expected value of ρ is computed as follows ( We are further interested in the probability that a certain number of proteins, say q, are dierential localised.Letting N U be the number of unlabelled proteins in the experiment, then the distribution of the prior number of dierential localised proteins is Computing δ is not simple; however, it is straightforward to estimate δ using Monte-Carlo, by simply sampling from Beta distributions: Thus, we recommend calibrating the Dirichlet prior using the above expectation and quantile.It may be important to calibrate several quantiles to ensure sucient mass is placed on desired regions 33 of the probability space.For example, let q 1 < q 2 , then we may desire that δ 1 , below, is not so small to rule out reasonable inferences and that δ 2 < δ 1 is suciently large.These can be computed from the equations below: More precise and informative prior biological knowledge can be specied; for example, should we suspect that some relocalisation events between particular organelles are more likely than others due to the stimuli, these can be encoded into the prior.If we expect more relocalisation events between organelle j and k 1 than organelle j and k 2 , this can be encoded by ensuring Alternatively, if an objective Bayesian analysis is preferred, the Jeery's prior sets α jk = 0.5 for every j, k = 1, .., K.This approach is not generally recommended by the authors, because the diagonal terms of π have a dierent interpretation to the o-diagonal terms.

Dierential localisation probability
The main posterior quantity of interest is the probability that a protein is dierentially localised.This can be approximated from the T Monte-Carlo samples as follows, suppressing notational dependence on all data and parameters for clarity where t denotes the t th of the MCMC algorithm.It is important to note that this quantity is agnostic to the assigned subcellular niche.We notice that the distribution of the number of MCMC observations for which z 1 is not equal to z 2 is given by: Hence, in this case, the Monte-Carlo estimator for χ is simply the maximum likelihood estimator of the probability parameter of the above binomial distribution.As T is given, uncertainty estimates, such as credible intervals, can be obtained from this binomial distribution directly.An alternative, but less computationally ecient approach, to perform uncertainty quantication on the the dierential localisation probability, could use the non-parametric bootstrap on the Monte-Carlo samples.More precisely, rst sample uniformly with replement from {z

Figure 1 :
Figure 1: An overview of the BANDLE workow.(A) A dierential localisation experiment is set-up with a perturbation of interest (B) Mass-spectrometry based spatial proteomics methods are applied to generate abundance proles across the subcellular fraction.(C) BANDLE is applied by rst calibrating the prior.The model is visualised as follows: each dataset is described as a mixture of subcellular niches, modelled as Gaussian random elds.Allocations are obtained for each condition and then integrated between the two conditions, using a joint prior.(D) The major results of BANDLE are represented in a rank plot.(E) Example alluvial plot generated from the results of BANDLE.
Figure 3: (A) Example PCA plots where pointers correspond to proteins.Marker proteins are coloured according to their subcellular niche, whilst proteins with unknown localisation are in grey.Simulated translocations are highlighted in black, where the left corresponds to control and right to the simulated perturbed dataset.(B) An MR-plot showing movement score against reproducibility score.Each pointer corresponds to a protein and orange pointers correspond to simulated translocations and blue otherwise.Teal lines are drawn at suggested thresholds with proteins in the top right corner considered hits.(C) A histogram of the raw statistics underlying the MR method.A Chi-square (orange) and Gamma (blue) t are overlaid (obtained using maximum likelihood estimation).The Gamma distribution clearly captures the tail behaviour.(D) A BANDLE rank plot where proteins are ranked from most to least likely to dierentially localised.The dierentially localisation probability is recorded on the y-axis.(right) A BANDLE rank plot of the top 30 dierentially localised proteins with uncertainty estimates for the dierential localisation probability.Proteins marked in orange were simulated translocations.(E) Violin plots for the dierential localisation probabilities (BANDLE), the M score (MR method) and R score (MR method).The distributiosn are split between dierential localised (movers) and spatially stable proteins.Clearly, the dierential localisation probabilities correlate most closely with the phenomena of interest.

Figure 4 :
Figure 4: (A) An MR-plot where dark green lines are drawn at suggested threshold and hits are highlighted in orange.(B) BANDLE rank plot showing the distribution of dierentially localised proteins.(C) The top dierentially localised proteins from BANDLE plotted with uncertainty estimates.(D) Bootstrap distributions of correlations with a phosphoprotemomic time-course experiment.The BANDLE condence intervals dier signicantly from 0, whilst the MR method do not.(E) PCA plots with (smoothed) localisation probabilities project onto them.Each colour represent an organelle and ellipses represent lines of isoprobability.The inner ellipse corresponds to 0.99 and the proceed line 0.95 with further lines decreasing by 0.05 each time.The annotated proteins demonstrate examples of dierential localisations.EGFR (P005330) clearly relocalises from the PM to endosome, whilst SHC-1 (P29353) and GRB2 (P62993) relocalise from unknown localisation to the lysosome.

Figure 6 :
Figure 6: (A) PCA plots with (smoothed) localisation probabilities projected onto them.Each colour represents a subcellular niche and ellipses represent lines of isoprobability.The inner ellipse corresponds to 0.99 and the proceed line 0.95 with further lines decreasing by 0.05 each time.The relocalisation of SCARB1 is highlighted on the plot (B) The spatial allocation derived from BANDLE where each entry of the heatmap is the number of proteins.Odiagonal entries only include condent dierential localisations with probability > 0.999 (C) UL148A interactome mapped onto the BANDLE determined spatial patterns.(D) UL8 interactome mapped onto the BANDLE determine spatial patterns.(E) A heatmap representation of the mean log2 fold changes in acetylation overlaid on spatial pattern of HCMV infection 24 hpi.
(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is made (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is made (which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity.It is made