Bayesian ANalysis of Differential Localisation Experiments (BANDLE) is an integrative semi-supervised functional mixture model, developed by (Crook et al. 2021), to obtain the probability of a protein being differentially localised between two conditions.
In this vignette we walk users through how to install and use the R
(R Development Core Team 2011)
Bioconductor (Gentleman et al. 2004) bandle
package by simulating a well-defined differential localisation
experiment from spatial proteomics data from the pRolocdata
package (Gatto et al. 2014).
The BANDLE method uses posterior Bayesian computations performed using Markov-chain Monte-Carlo (MCMC) and thus uncertainty estimates are available (Gilks, Richardson, and Spiegelhalter 1995). Throughout this vignette we use the term differentially localised to pertain to proteins which are assigned to different sub-cellular localisations between two conditions. One of the main outputs of BANDLE is the probability that a protein is differentially localised between two conditions.
The package can be installed with the BiocManager
package:
if (!requireNamespace("BiocManager", quietly=TRUE))
install.packages("BiocManager")
BiocManager::install("bandle")
and then loaded,
For visualisation we also load the packages,
In this vignette and Crook et al.
(2021), the main data source that we use to study differential
protein sub-cellular localisation are data from high-throughput mass
spectrometry-based experiments. The data from these types of experiments
traditionally yield a matrix of measurements wherein we have, for
example, PSMs, peptides or proteins along the rows, and
samples/channels/fractions along the columns. The bandle
package uses the MSnSet
class as implemented in the
Bioconductor MSnbase
package and thus requires users to import and store their data as a
MSnSet
instance. For more details on how to create a
MSnSet
see the relevant vignettes in pRoloc.
The pRolocdata
experiment data package is a good starting place to look for test data.
This data package contains tens of quantitative proteomics experiments,
stored as MSnSet
s.
To get started with the basics of using bandle
we begin
by generating a simple example dataset which simulates a differential
localisation experiment (please see the second vignette in this package
for a full real-life biological use case). In this example data, the key
elements are replicates, and a perturbation of interest. There is code
within the bandle
package to simulate an example experiment.
In the code chunk below we begin by loading the pRolocdata
package to obtain a spatial proteomics dataset. This will be the basis
of our simulation which will use boostrapping to generate new datasets.
The dataset we have chosen to load is a dataset from 2009
(tan2009r1
). This is data from a early LOPIT experiment
performed on Drosophila embryos (Tan et al.
2009). The aim of this experiment was to apply LOPIT to an
organism with heterogeneous cell types. This experiment used four
isotopes across four distinct fractions and thus yielded four
measurements (features) per protein profile. We visualise the data by
using principal components analysis.
library("pRolocdata")
data("tan2009r1")
## Let's set the stock colours of the classes to plot to be transparent
setStockcol(NULL)
setStockcol(paste0(getStockcol(), "90"))
## Plot the data using plot2D from pRoloc
plot2D(tan2009r1,
main = "An example spatial proteomics datasets",
grid = FALSE)
addLegend(tan2009r1, where = "topleft", cex = 0.7, ncol = 2)
The following code chuck simulates a differential localisation
experiment. It will generate numRep/2
of each a control and
treatment condition. We will also simulate relocalisations for
numDyn
proteins.
## [1] "markers"
The list of the 6 simulated experiments are found in
tansim$lopitrep
. Each one is an MSnSet
instance (the standard data container for proteomics experimental data).
The first 3 are the simulated control experiments (see
tansim$lopitrep[1:3]
), and the following 3 in the list are
the treatment condition simulated experiments (see
tansim$lopitrep[4:6]
). We can plot them using the
plot2D
function from pRoloc
.
plot_title <- c(paste0("Replicate ", seq(3), " condition", " A"),
paste0("Replicate ", seq(3), " condition", " B"))
par(mfrow = c(2, 3))
out <- lapply(seq(tansim$lopitrep), function(z)
plot2D(tansim$lopitrep[[z]], grid = FALSE, main = plot_title[z]))
For understanding, exploring and visualizing individual spatial
proteomics experiments, see the vignettes in pRoloc
and
MSnbase
packages.
## MSnSet (storageMode: lockedEnvironment)
## assayData: 888 features, 4 samples
## element names: exprs
## protocolData: none
## phenoData
## sampleNames: X114 X115 X116 X117
## varLabels: Fractions
## varMetadata: labelDescription
## featureData
## featureNames: P20353 P53501 ... P07909 (888 total)
## fvarLabels: FBgn Protein.ID ... knn.scores (18 total)
## fvarMetadata: labelDescription
## experimentData: use 'experimentData(object)'
## pubMedIds: 19317464
## Annotation:
## - - - Processing information - - -
## Added markers from 'mrk' marker vector. Thu Jul 16 22:53:44 2015
## Performed knn prediction (k=10) Tue Sep 10 05:14:03 2024
## MSnbase version: 1.17.12
bandle
analysisThe main function of the package is bandle
, this uses a
complex model to analyse the data. Markov-Chain Monte-Carlo (MCMC) is
used to sample the posterior distribution of parameters and latent
variables. From which statistics of interest can be computed. Here we
only run a few iterations for brevity but typically one needs to run
thousands of iterations to ensure convergence, as well as multiple
parallel chains.
First, we need to fit non-parametric regression functions to the
markers profiles, upon which we place our analysis. This uses Gaussian
processes. The fitGPmaternPC
function can be used and fits
some default penalised complexity priors (see ?fitGP
),
which works well. However, these can be altered, which is demonstrated
in the next code chunk
par(mfrow = c(4, 3))
gpParams <- lapply(tansim$lopitrep, function(x)
fitGPmaternPC(x, hyppar = matrix(c(10, 60, 250), nrow = 1)))
We apply the fitGPmaternPC
function to each datasets by
calling lapply
over the tansim$lopitrep
list
of datasets. The output of fitGPmaternPC
returns a list of
posterior predictive means and standard deviations. As well as MAP
hyperparamters for the GP.
Note here we the use the default
hyppar = matrix(c(10, 60, 250), nrow = 1)
as a starting
point for fitting the GPs to the marker profile distributions. In the
Crook et al. (2021) manuscript we found
that these values worked well for smaller spatial proteomics datasets.
This was visually assessed by passing these values and visualising the
GP fit using the plotGPmatern
function.
The plotGPmatern
function can be used to plot the
profiles for each class in each replicate condition with the posterior
predictive distributions overlayed with the markers protein
profiles.
For example, to plot the predictive distributions of the first dataset,
The prior needs to form a K*3
matrix. K
corresponds to the number of subcellular classes in the data, and 3
columns for (1) the prior, (2) length-scale amplitude and (3) standard
deviation parameters (see hyppar
in
?fitGPmaternPC
). Increasing these values, increases the
shrinkage. For more details see the manuscript by Crook et al. (2021). We strongly recommend users
start with the recommended hyppar
parameters and change and
assess them as necessary for their dataset by visually evaluating the
fit of the GPs using the plotGPmatern
function.
K <- length(getMarkerClasses(tansim$lopitrep[[1]], fcol = "markers"))
pc_prior <- matrix(NA, ncol = 3, K)
pc_prior[seq.int(1:K), ] <- matrix(rep(c(10, 60, 250),
each = K), ncol = 3)
Now we have generated these complexity priors we can pass them as an
argument to the fitGPmaternPC
function. For example,
By looking at the plot of posterior predictives using the
gpParams
we can see the GP fit looks sensible.
The next step is to set up the matrix Dirichlet prior on the mixing
weights. These weights are defined across datasets so these are slightly
different to mixture weights in usual mixture models. The (i, j)th
entry is the prior probability that a protein localises to organelle
i in the control and j in the treatment. This mean that
off-diagonal terms have a different interpretation to diagonal terms.
Since we expect re-localisation to be rare, off-diagonal terms should be
small. The following functions help set up the priors and how to
interpret them. The parameter q
allow us to check the prior
probability that more than q
differential localisations are
expected.
set.seed(1)
dirPrior = diag(rep(1, K)) + matrix(0.001, nrow = K, ncol = K)
predDirPrior <- prior_pred_dir(object = tansim$lopitrep[[1]],
dirPrior = dirPrior,
q = 15)
The mean number of re-localisations is small:
## [1] 0.3457542
The prior probability that more than q
differential
localisations are expected is small
## [1] 4e-04
The full prior predictive can be visualised as histogram. The prior probability that proteins are allocated to different components between datasets concentrates around 0.
For most use-cases we indeed expect the number of differential
localisations to be small. However, there may be specific cases where
one may expect the either a smaller or larger number of differential
localisations. Users could try testing different values for the
dirPrior
for example, replacing 0.001 with 0.0005 or
smaller, for larger datasets to bring the number of expected
re-localisations inline with the biological expectation, and visa-versa
when we expect the number of proteins to have changed to be higher.
bandle
functionWe are now ready to run the main bandle
function.
Remember to carefully select the datasets and replicates that define the
control and treatment. As a reminder, in this introductory vignette we
have used a small dataset and generated theoretical triplicates of each
theoretical condition. Please see the second vignette in this package
for a more detailed workflow and real biological use-case. In the below
code chunk we run bandle
for only 50 iterations for the
convenience of building the vignette, but typically we’d recommend you
run the number of iterations (numIter
) in the 1000s.
Remember: the first 3 datasets are the first 3 elements of
tansim
and the final 3 elements are the “treatment”
triplicate datasets.
control <- tansim$lopitrep[1:3]
treatment <- tansim$lopitrep[4:6]
bandleres <- bandle(objectCond1 = control,
objectCond2 = treatment,
numIter = 20, # usually 10,000
burnin = 5L, # usually 5,000
thin = 1L, # usually 20
gpParams = gpParams,
pcPrior = pc_prior,
numChains = 1, # usually >=4
dirPrior = dirPrior)
The bandle
function generates an object of class
bandleParams
. The show
method indicates the
number of parallel chains that were run, this should typically be
greater than 4 (here we use 1 just as a demo). Normally, we should also
assess convergence but this omitted for the moment so that we can move
forward with the analysis. Please see the end of the vignette for
convergence plots.
## Object of class "bandleParams"
## Method: bandle
## Number of chains: 1
bandle
outputBefore we can begin to extract protein allocation information and a
list of proteins which are differentially localised between conditions,
we first need to populate the bandleres
object by calling
the bandleProcess
function.
bandleres
objectCurrently, the summary slots of the bandleres
object are
empty. The summaries
function accesses them.
## [[1]]
## An object of class "bandleSummary"
## Slot "posteriorEstimates":
## <S4 Type Object>
## attr(,"elementType")
## [1] "ANY"
## attr(,"elementMetadata")
## `\001NULL\001`
## attr(,"metadata")
## list()
##
## Slot "diagnostics":
## <0 x 0 matrix>
##
## Slot "bandle.joint":
## <0 x 0 matrix>
##
##
## [[2]]
## An object of class "bandleSummary"
## Slot "posteriorEstimates":
## <S4 Type Object>
## attr(,"elementType")
## [1] "ANY"
## attr(,"elementMetadata")
## `\001NULL\001`
## attr(,"metadata")
## list()
##
## Slot "diagnostics":
## <0 x 0 matrix>
##
## Slot "bandle.joint":
## <0 x 0 matrix>
These can be populated as follows
These slots have now been populated
## Length Class Mode
## [1,] 1 bandleSummary S4
## [2,] 1 bandleSummary S4
bandle
resultsWe can save the results by calling summaries
. We see
that it is of length 2. 1 for control and 1 for treatment.
## [1] 2
There are a number of slots,
## Formal class 'bandleSummary' [package "bandle"] with 3 slots
## ..@ posteriorEstimates:Formal class 'DFrame' [package "S4Vectors"] with 6 slots
## .. .. ..@ rownames : chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. ..@ nrows : int 677
## .. .. ..@ elementType : chr "ANY"
## .. .. ..@ elementMetadata: NULL
## .. .. ..@ metadata : list()
## .. .. ..@ listData :List of 7
## .. .. .. ..$ bandle.allocation : Named chr [1:677] "PM" "Golgi" "Ribosome 40S" "PM" ...
## .. .. .. .. ..- attr(*, "names")= chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. .. ..$ bandle.probability : Named num [1:677] 1 0.996 1 0.971 0.989 ...
## .. .. .. .. ..- attr(*, "names")= chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. .. ..$ bandle.outlier : num [1:677] 0 0 0 0 0 0 0 0 0 0 ...
## .. .. .. ..$ bandle.probability.lowerquantile: Named num [1:677] 1 0.991 1 0.947 0.963 ...
## .. .. .. .. ..- attr(*, "names")= chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. .. ..$ bandle.probability.upperquantile: Named num [1:677] 1 0.999 1 0.992 0.998 ...
## .. .. .. .. ..- attr(*, "names")= chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. .. ..$ bandle.mean.shannon : Named num [1:677] 0.00 0.00 1.66e-07 0.00 0.00 ...
## .. .. .. .. ..- attr(*, "names")= chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. .. ..$ bandle.differential.localisation: Named num [1:677] 0 0 0 0 0 0 0 0 0 0 ...
## .. .. .. .. ..- attr(*, "names")= chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## ..@ diagnostics : logi [1, 1] NA
## ..@ bandle.joint : num [1:677, 1:11] 2.25e-97 5.82e-223 1.75e-113 1.77e-126 2.09e-60 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:677] "P20353" "P53501" "Q7KU78" "P04412" ...
## .. .. ..$ : chr [1:11] "Cytoskeleton" "ER" "Golgi" "Lysosome" ...
The main one of interest is the posteriorEstimates
slot,
## DataFrame with 677 rows and 7 columns
## bandle.allocation bandle.probability bandle.outlier
## <character> <numeric> <numeric>
## P20353 PM 0.999998 0
## P53501 Golgi 0.996345 0
## Q7KU78 Ribosome 4... 1.000000 0
## P04412 PM 0.971033 0
## Q7KJ73 Ribosome 6... 0.988720 0
## ... ... ... ...
## Q95TL8 ER 1.000000 0
## P25007 Proteasome 0.999996 0
## P41374 Cytoskelet... 0.999991 0
## Q8SZM1 Peroxisome 0.912287 0
## P07909 Nucleus 0.999999 0
## bandle.probability.lowerquantile bandle.probability.upperquantile
## <numeric> <numeric>
## P20353 0.999996 1.000000
## P53501 0.990644 0.999324
## Q7KU78 1.000000 1.000000
## P04412 0.947008 0.991517
## Q7KJ73 0.963326 0.997742
## ... ... ...
## Q95TL8 1.000000 1.000000
## P25007 0.999992 0.999999
## P41374 0.999980 0.999999
## Q8SZM1 0.768171 0.991368
## P07909 0.999996 1.000000
## bandle.mean.shannon bandle.differential.localisation
## <numeric> <numeric>
## P20353 0.00000e+00 0
## P53501 0.00000e+00 0
## Q7KU78 1.65884e-07 0
## P04412 0.00000e+00 0
## Q7KJ73 0.00000e+00 0
## ... ... ...
## Q95TL8 0.00000e+00 0
## P25007 0.00000e+00 0
## P41374 1.52728e-05 0
## Q8SZM1 0.00000e+00 1
## P07909 0.00000e+00 0
This output object is a data.frame
containing the
protein allocations and associated localisation probabilities (including
the upper and lower quantiles of the allocation probability
distribution), the mean Shannon entropy and the
bandle.differential.localisation
probability.
We create two new objects pe1
and pe2
in
the below code chunk which contain the output of the
posteriorEstimates
slot.
One quantity of interest is the protein allocations, which we can plot as a barplot.
alloc1 <- pe1$bandle.allocation
alloc2 <- pe2$bandle.allocation
par(mfrow = c(1, 2), oma = c(6,2,2,2))
barplot(table(alloc1), col = getStockcol()[2],
las = 2, main = "Protein allocation: control")
barplot(table(alloc2), col = getStockcol()[2],
las = 2, main = "Protein allocation: treatment")
The barplot tells us for this example that bandle
has
allocated the majority of unlabelled proteins to the ER, followed by the
Golgi (irrespective of the posterior probability).
The associated posterior estimates are located in the
bandle.probability
column.
The full allocation probabilities are stored in the
tagm.joint
slot. These can be visualised in a heatmap
We can append the results to our original MSnSet
datasets using the bandlePredict
function.
xx <- bandlePredict(control,
treatment,
params = bandleres,
fcol = "markers")
res_control <- xx[[1]]
res_treatment <- xx[[2]]
The output is a list
of MSnSets
. In this
example, we have 3 for the control and 3 for the treatment.
## [1] 3
## [1] 3
The results are appended to the first
MSnSet
feature data slot for each condition.
## [1] "FBgn" "Protein.ID"
## [3] "Flybase.Symbol" "AccessionNo"
## [5] "EntryName" "AccessionNoAll"
## [7] "EntryNameAll" "No.peptide.IDs"
## [9] "Mascot.score" "No.peptide.quantified"
## [11] "PLSDA" "pd.2013"
## [13] "pd.markers" "markers.orig"
## [15] "markers" "markers.tl"
## [17] "knn" "knn.scores"
## [19] "bandle.allocation" "bandle.probability"
## [21] "bandle.probability.lowerquantile" "bandle.probability.upperquantile"
## [23] "bandle.mean.shannon" "bandle.differential.localisation"
## [25] "bandle.outlier" "bandle.joint"
To access them use the fData
function
It is common practice in supervised machine learning to set a
specific threshold on which to define new assignments/allocations, below
which classifications are left unassigned/unknown. Indeed, we do not
expect the whole subcellular diversity to be represented by the 11
niches defined here, we expect there to be many more, many of which will
be multiply localised within the cell. It is important to allow for the
possibility of proteins to reside in multiple locations (this
information is available in the bandle.joint
slot - see
above for more details on multiple location).
As we are using a Bayesian model the outputs of the classifier are probabilities. This not only allows us to look at the distribution of probabilities over all subcellular classes but also allows us to extract a probability threshold on which we can define new assignments.
The subcellular allocations are located in the
bandle.allocation
column of the fData
slot and
the posteriors are located in the bandle.probability
slot.
We can use the getPredictions
function from the
pRoloc
package to return a set of predicted localisations
according to if they meet a probability threshold.
For example, in the below code chunk we set a 1% FDR for assigning proteins a subcellular nice, below which we leave them unlabelled.
res_control[[1]] <- getPredictions(res_control[[1]],
fcol = "bandle.allocation",
scol = "bandle.probability",
mcol = "markers",
t = .99)
## ans
## Cytoskeleton ER Golgi Lysosome Nucleus
## 13 216 132 8 39
## PM Peroxisome Proteasome Ribosome 40S Ribosome 60S
## 140 8 43 73 50
## mitochondrion unknown
## 82 84
res_treatment[[1]] <- getPredictions(res_treatment[[1]],
fcol = "bandle.allocation",
scol = "bandle.probability",
mcol = "markers",
t = .99)
## ans
## Cytoskeleton ER Golgi Lysosome Nucleus
## 18 191 121 12 42
## PM Peroxisome Proteasome Ribosome 40S Ribosome 60S
## 142 14 48 79 48
## mitochondrion unknown
## 86 87
We may also wish to take into account the probability of the protein
being an outlier and thus use the results in the
bandle.outlier
column of the feature data. We could
calculate the product of the posterior and the outlier (as they are both
probabilities) to obtain a localisation score that takes into account
the outlier model. More details on this are found in the second vignette
of this package.
As previously mentioned the term “differentially localised” is used to pertain to proteins which are assigned to different sub-cellular localisations between two conditions. For the majority of users this is the main output they are keen to extract using the BANDLE method.
Following on from the above example, after extracting posterior
estimates for each condition using the summaries
function
we can also access the differential localisation probability as it is
stored in the bandle.differential.localisation
column of
the data.frames
of pe1
and pe2
,
in the above sections.
The differential localisation probability tells us which proteins are most likely to differentially localise. We can for example, examine how many proteins get a differential probability greater than 0.99 to look for the most confident differentially localised candidates.
## P20353 P53501 Q7KU78 P04412 Q7KJ73 Q9VM65 Q9VCK0
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## B7Z0W3 Q9V415 Q00174 Q9V769 Q27593 Q9V780 P19109
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## Q95NR4 O77047 Q8SXD0 Q9VBV5 Q9VBU5 Q8IA62 Q9Y105
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.06666667
## Q9W2M4 P26308 Q9VP77 Q5U0Z2 Q960E2 Q9VLJ6 Q9VJ39
## 0.00000000 1.00000000 1.00000000 0.00000000 0.00000000 1.00000000 1.00000000
## Q9VTX8 Q9VTZ5 B7Z0C1 Q9VRJ4 M9PCB7 P46150 A1ZBH5
## 0.00000000 0.33333333 0.00000000 0.00000000 0.00000000 0.06666667 0.00000000
## Q9W3M7 A8DZ29 Q9VN14 Q9VZL3 M9PC99 Q86BP3 Q9W3N9
## 0.00000000 0.73333333 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000
## Q7K5M6 P16620 P48375 Q9VMB9 Q9VI55 Q9VU58 Q9V498
## 0.00000000 0.00000000 0.00000000 0.00000000 0.00000000 0.40000000 0.00000000
## Q9V496
## 1.00000000
## [1] 90
We find there are 90 proteins above this threshold.
This can also be seen on a rank plot
plot(diffloc_probs[order(diffloc_probs, decreasing = TRUE)],
col = getStockcol()[3], pch = 19, ylab = "Probability",
xlab = "Rank", main = "Differential localisation rank plot")
In-line with our expectations, the rank plot indicates that most proteins are not differentially localised.
bootstrapdiffLocprob
functionWe can examine the top n
proteins (here we use an
example of top = 100
) and produce bootstrap estimates of
the uncertainty (note here the uncertainty is likely to be
underestimated as we did not produce many MCMC samples). These can be
visualised as ranked boxplots.
binomDiffLoc
functionInstead of applying the bootstrapdiffLocprob
we could
use the binomDiffLoc
function to obtain credible intervals
from the binomial distribution.
There are many ways we could obtain probability estimates from either of the above methods. We could, for example, take the mean of each protein estimate, or compute the cumulative error (there is not really a false discovery rate in Bayesian statistics) or we could threshold on the interval to reduce the number of differential localisations if you feel the model has been overconfident.
# more robust estimate of probabilities
dprobs <- rowMeans(bin_t)
# compute cumulative error, there is not really a false discovery rate in
# Bayesian statistics but you can look at the cumulative error rate
ce <- cumsum(1 - dprobs)
# you could threshold on the interval and this will reduce the number of
# differential localisations
qt <- apply(bin_t, 1, function(x) quantile(x, .025))
Instead of estimating the false discovery rate we can estimate the expected false discovery rate from the posterior probabilities at a particular threshold. This mean that for fixed threshold, we compute the expected proportion of false discoveries. Here is an example below. We can see that setting a probability threshold of 0.9 leads to an expected false discovery rate of less than 0.5%
## [1] 0.004166667
We can visualise the changes in localisation between conditions on an
alluvial plot using the plotTranslocations
function
Or alternatively, on a chord (circos) diagram
These visualisations are showing the change in class label between
the two conditions (as assigned by bandle
i.e. the result
stored in bandle.allocation
). The results are taken
directly from bandleres
and thus no thresholding on the
class label and posterior to allow for proteins to be left “unassigned”
or unknown, is conducted. Furthermore, there is no thresholding on the
bandle.differential.localisation
probability.
It would be better to re-plot these figures with some thresholds on
the above quantities, to get a better representation of what is moving.
The easiest way to do this is to pass the MSnSets
output
after performing bandlePredict
and
getPredictions
.
For example, first let’s identify which proteins get a high differential localisation probability,
## identify which proteins get a high differential localisation probability
ind <- which(fData(res_control[[1]])$bandle.differential.localisation > 0.99)
## create two new MSnSets with only these proteins
res_dl_control <- res_control[[1]][ind, ]
res_dl_treatment <- res_treatment[[1]][ind, ]
Now we can plot only these 90 proteins which are deemed to
move/differentially localise. We also specify where the prediction
results are located
e.g. fcol = "bandle.allocation.pred"
.
## specify colour palette
mycols <- c(getStockcol()[seq(getMarkerClasses(res_control[[1]]))], "grey")
names(mycols) <- c(getMarkerClasses(res_control[[1]]), "unknown")
## Create a list of the datasets for plotTranslocations
res <- list(res_dl_control, res_dl_treatment)
plotTranslocations(res, fcol = "bandle.allocation.pred", col = mycols)
## 90 features in common
## ------------------------------------------------
## If length(fcol) == 1 it is assumed that the
## same fcol is to be used for both datasets
## setting fcol = c(bandle.allocation.pred,bandle.allocation.pred)
## ----------------------------------------------
We can also use the function plotTable
to display a
summary table of the number of proteins that have changed in location
between conditions.
## Condition1 Condition2 value
## 1 ER Golgi 1
## 2 ER Nucleus 1
## 3 ER PM 3
## 4 ER Peroxisome 7
## 5 ER Proteasome 4
## 6 ER Ribosome 40S 4
## 7 ER Ribosome 60S 1
## 9 ER unknown 6
## 10 ER Cytoskeleton 2
## 11 ER Lysosome 1
## 12 Golgi ER 1
## 13 Golgi Nucleus 2
## 18 Golgi Ribosome 60S 2
## 19 Golgi mitochondrion 2
## 25 Nucleus PM 1
## 26 Nucleus Peroxisome 1
## 35 PM Golgi 1
## 38 PM Proteasome 2
## 39 PM Ribosome 40S 4
## 42 PM unknown 2
## 44 PM Lysosome 1
## 46 Peroxisome Golgi 2
## 59 Proteasome PM 2
## 61 Proteasome Ribosome 40S 2
## 62 Proteasome Ribosome 60S 1
## 63 Proteasome mitochondrion 3
## 64 Proteasome unknown 1
## 68 Ribosome 40S Golgi 1
## 69 Ribosome 40S Nucleus 1
## 70 Ribosome 40S PM 1
## 74 Ribosome 40S mitochondrion 1
## 76 Ribosome 40S Cytoskeleton 1
## 77 Ribosome 40S Lysosome 2
## 83 Ribosome 60S Proteasome 2
## 84 Ribosome 60S Ribosome 40S 2
## 85 Ribosome 60S mitochondrion 1
## 94 mitochondrion Proteasome 1
## 97 mitochondrion unknown 1
## 98 mitochondrion Cytoskeleton 1
## 100 unknown ER 2
## 101 unknown Golgi 2
## 102 unknown Nucleus 1
## 103 unknown PM 1
## 105 unknown Proteasome 3
## 106 unknown Ribosome 40S 1
## 109 unknown Cytoskeleton 1
In this section, we demonstrate how to visually assess convergence of the MCMC algorithm. In the chunk below, we use 4 chains so that we can assess the convergence of the method.
control <- tansim$lopitrep[1:3]
treatment <- tansim$lopitrep[4:6]
bandleres <- bandle(objectCond1 = control,
objectCond2 = treatment,
numIter = 50, # usually 10,000
burnin = 10L, # usually 5,000
thin = 1L, # usually 20
gpParams = gpParams,
pcPrior = pc_prior,
numChains = 4,
dirPrior = dirPrior)
We then use the plotConvergence
function which plots the
ranks of the total number of outliers within in each of the chains. If
convergence has been reached these plots should be uniform. As we can
see these plots are not uniform and are skewed towards the extremes, so
they have not converged. Clearly one of the chains has higher values on
average than the other chains producing the skew. Running the algorithm
for more iterations (typically 10,000) should produce convergence. This
is not done here for brevity of the vignette.
bandle
parametersThe bandle
function has a significant number of
parameters to allow flexible and bespoke analysis. Here, we describe
these parameters in more detail to allow user to make decisions on the
level of flexibility they wish to exploit.
objectCond1
. This is a list of MSnSets
containing the first condition.
objectCond2
. This is a list of MSnSets
containing the second condition.
fcol
indicates the feature column in the
MSnSets
that indicated the markers. Proteins that are not
markers should be labels unknown
. The default is
markers
.
hyperLearn
is the algorithm used to learn the
hyperparameters of the Gaussian processes. For speed the default is an
optimization algorithm called “LBFGS”, however is users want to perform
uncertainty quantification on these parameters we can use Markov-chain
Monte Carlo (MCMC) methods. This is implemented using the
Metropolis-Hastings algorithm. Though this latter methodology provides
more information, it is much more costly. The analysis is expected to
take several days rather than hours.
numIter
is the number of MCMC iterations for the
algorithm. We typically suggest around 10,000 iterations is plenty for
convergence. Though some cases may take longer. If resources are
constrained, we suggest 4,000 iterations as acceptable. A minimum number
of iterations is around 1,000 though at this level we expect the
posterior estimates to suffer considerably. If possible more parallel
chains should be run in this case by changing numChains
to,
say, 9. The more chains and iterations the more computationally
expensive the algorithm. The time taken for the algorithm scales roughly
linearly in the number of iterations
burnin
is the number of samples that should be
discarded from the beginning of the chain due to the bias induced by the
starting point of the algorithm. We suggest sensible burnin
values to be roughly 10 − 50% of the
number of iterations
thin
reduces auto-correlation in the MCMC samples.
The default is 5, which means every 5th
sample is taken. If memory requirements are an issue, we suggest to
increase the thinning amount. Though above 20, you will see a decrease in
performance.
u
and v
represent the prior
hyperparameters of the proportion of outliers. This is modelled using a
Beta(u,v)
with u = 2
and v = 10
a
default. This suggest that roughly $\frac{u}{u
= V} = 16%$ of proteins are believed to be outliers and that it
is quite unlikely that more than 50 of
proteins are outliers. Users can examine the quantiles of the
Beta(u,v)
distribution if they wish to place a more bespoke
prior. For example, increasing u
will increase the number
of a prior believed outliers.
lambda
is a ridge parameter used for numerical
stability and is set to 0.01. If you
experience the algorithm fails due to numerical issue then you can set
this value larger. If you require values above 1 it is likely that there are other issues
with the analysis. We suggest checking the method is appropriate for
your problem and opening issue detailing the problems.
gpParams
results from fitting Gaussian proccess
(Gaussian random fields). We refer the users to those functions. The
default is NULL
which will fit GPs internally but we
recommend setting these outside the bandle function because it leads to
more stable results.
hyperIter
if the hyperparameters of the GP are
learnt using MH algorithm then this is the frequency at which these are
updated relative to the bandle algorithm. By default this is unused, but
if hyperLearn
is set to MH
then this proceed
at every 20 iterations.
hyperMean
is the mean of the log normal prior used
on the hyperparameters. Though by default this is not used unless
PC
is set to false
hyperSd
is the standard deviation of the log normal
prior used on the hyperparameters. The default is c(1,1,1)
for the 3 hyperparameters, increasing these values increases the
uncertainty in the prior values of the hyperparameters.
seed
is the random-number seed.
pg
indicates whether or not to use the Polya-Gamma
(PG) prior. The default is false and a Dirichlet prior is used instead.
If set to true the pg
is used. In which case a default PG
prior is used. This prior attempts to match the default Dirichlet prior
that is used when PG prior is set to false. The PG prior is more
computationally expensive but can provide prior information on
correlations
pgPrior
is by default NULL. We suggest using the
pg_prior
function to help set this parameter and the
documentation therein. This function uses an empirical approach to
compute a sensible default.
tau
is a parameter used by the Polya-Gamma prior and
we refer to BANDLE manuscript for details. By default it is only used if
pg
prior is true, when the default becomes
0.2
. At this value the pg
prior is similar to
the Dirichlet prior but with information on correlations.
dirPrior
is the Dirichlet matrix prior on the
correlations. This should be provided as a K by K matrix, where K is the
number of subcellular niches. The diagonal component should represent
the prior belief that organelles do not re-localise (same compartment),
where as the off-diagonal terms represent the prior terms of
re-localisation. The prior_pred_dir
can be used to provide
a prior predictive check based on the provided prior. It is recommended
that the off-diagonal terms are at least two orders of magnitude smaller
than the diagonal terms. An example is given in the vignette.
maternCov
is this true the covariance function is
the matern covariance, otherwise a Gaussian covariance is used.
PC
indicates whether a penalised complexity (PC) is
used. The default is true and otherwise log normal priors are
used.
pcPrior
is a numeric of length 3 indicating the
parameters of the PC prior. The prior is placed on the parameters of
length-scale, amplitude, and variance in that order. The default values
are 0.5, 3, 100, and increasing the
value increases the shrinkage towards straight-lines with zero
variance.
nu
which defaults to 2 is the smoothness of the
matern covariance. By increasing nu
you encourage smoother
solutions. nu
should be an integer, though for values of
nu
above 3, we have observed numerical
instability.
propSd
is the standard deviation of the random-walk
update used in the MH algorithm. We do not recommend changing this
unless you are familiar with Bayesian analysis. The default is
c(0.3,0.1,0.05)
for the 3 hyperparameters. Changing these
will alter the efficiency of the underlying samplers.
numChains
is the number of parrallel chains and
defaults to 4. We recommend using as much processing resources as you
have and frequently have used 9 in practise.
BPPARAM
is the BiocParallel back-end which defaults
to BiocParallel::bpparam()
. We refer you to the BiocParallel
package for details on setting this dependent on your computing
system.
All software and respective versions used to produce this document are listed below.
## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so; LAPACK version 3.12.0
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## time zone: Etc/UTC
## tzcode source: system (glibc)
##
## attached base packages:
## [1] stats4 stats graphics grDevices utils datasets methods
## [8] base
##
## other attached packages:
## [1] pRolocdata_1.43.0 viridis_0.6.5 viridisLite_0.4.2
## [4] pheatmap_1.0.12 bandle_1.9.0 pRoloc_1.45.1
## [7] BiocParallel_1.39.0 MLInterfaces_1.85.0 cluster_2.1.6
## [10] annotate_1.83.0 XML_3.99-0.17 AnnotationDbi_1.67.0
## [13] IRanges_2.39.2 MSnbase_2.31.1 ProtGenerics_1.37.1
## [16] mzR_2.39.0 Rcpp_1.0.13 Biobase_2.65.1
## [19] S4Vectors_0.43.2 BiocGenerics_0.51.1 BiocStyle_2.33.1
##
## loaded via a namespace (and not attached):
## [1] splines_4.4.1 filelock_1.0.3
## [3] tibble_3.2.1 hardhat_1.4.0
## [5] preprocessCore_1.67.0 pROC_1.18.5
## [7] rpart_4.1.23 lifecycle_1.0.4
## [9] httr2_1.0.3 doParallel_1.0.17
## [11] globals_0.16.3 lattice_0.22-6
## [13] MASS_7.3-61 MultiAssayExperiment_1.31.5
## [15] dendextend_1.17.1 magrittr_2.0.3
## [17] limma_3.61.9 plotly_4.10.4
## [19] sass_0.4.9 rmarkdown_2.28
## [21] jquerylib_0.1.4 yaml_2.3.10
## [23] MsCoreUtils_1.17.1 DBI_1.2.3
## [25] buildtools_1.0.0 RColorBrewer_1.1-3
## [27] lubridate_1.9.3 abind_1.4-5
## [29] zlibbioc_1.51.1 GenomicRanges_1.57.1
## [31] purrr_1.0.2 mixtools_2.0.0
## [33] AnnotationFilter_1.29.0 nnet_7.3-19
## [35] rappdirs_0.3.3 ipred_0.9-15
## [37] circlize_0.4.16 lava_1.8.0
## [39] GenomeInfoDbData_1.2.12 ggrepel_0.9.6
## [41] listenv_0.9.1 gdata_3.0.0
## [43] maketools_1.3.0 parallelly_1.38.0
## [45] ncdf4_1.23 codetools_0.2-20
## [47] DelayedArray_0.31.11 xml2_1.3.6
## [49] shape_1.4.6.1 tidyselect_1.2.1
## [51] farver_2.1.2 UCSC.utils_1.1.0
## [53] matrixStats_1.4.1 BiocFileCache_2.13.0
## [55] jsonlite_1.8.8 caret_6.0-94
## [57] e1071_1.7-14 ggalluvial_0.12.5
## [59] survival_3.7-0 iterators_1.0.14
## [61] foreach_1.5.2 segmented_2.1-2
## [63] tools_4.4.1 progress_1.2.3
## [65] lbfgs_1.2.1.2 glue_1.7.0
## [67] prodlim_2024.06.25 gridExtra_2.3
## [69] SparseArray_1.5.32 xfun_0.47
## [71] MatrixGenerics_1.17.0 GenomeInfoDb_1.41.1
## [73] dplyr_1.1.4 withr_3.0.1
## [75] BiocManager_1.30.25 fastmap_1.2.0
## [77] fansi_1.0.6 digest_0.6.37
## [79] timechange_0.3.0 R6_2.5.1
## [81] colorspace_2.1-1 gtools_3.9.5
## [83] lpSolve_5.6.20 biomaRt_2.61.3
## [85] RSQLite_2.3.7 utf8_1.2.4
## [87] tidyr_1.3.1 generics_0.1.3
## [89] hexbin_1.28.4 data.table_1.16.0
## [91] recipes_1.1.0 FNN_1.1.4
## [93] class_7.3-22 prettyunits_1.2.0
## [95] PSMatch_1.9.0 httr_1.4.7
## [97] htmlwidgets_1.6.4 S4Arrays_1.5.7
## [99] ModelMetrics_1.2.2.2 pkgconfig_2.0.3
## [101] gtable_0.3.5 timeDate_4032.109
## [103] blob_1.2.4 impute_1.79.0
## [105] XVector_0.45.0 sys_3.4.2
## [107] htmltools_0.5.8.1 MALDIquant_1.22.3
## [109] clue_0.3-65 scales_1.3.0
## [111] png_0.1-8 gower_1.0.1
## [113] knitr_1.48 reshape2_1.4.4
## [115] coda_0.19-4.1 nlme_3.1-166
## [117] curl_5.2.2 GlobalOptions_0.1.2
## [119] proxy_0.4-27 cachem_1.1.0
## [121] stringr_1.5.1 parallel_4.4.1
## [123] mzID_1.43.0 vsn_3.73.0
## [125] pillar_1.9.0 grid_4.4.1
## [127] vctrs_0.6.5 pcaMethods_1.97.0
## [129] randomForest_4.7-1.1 dbplyr_2.5.0
## [131] xtable_1.8-4 evaluate_0.24.0
## [133] mvtnorm_1.3-1 cli_3.6.3
## [135] compiler_4.4.1 rlang_1.1.4
## [137] crayon_1.5.3 future.apply_1.11.2
## [139] LaplacesDemon_16.1.6 mclust_6.1.1
## [141] QFeatures_1.15.2 affy_1.83.0
## [143] plyr_1.8.9 stringi_1.8.4
## [145] munsell_0.5.1 Biostrings_2.73.1
## [147] lazyeval_0.2.2 Matrix_1.7-0
## [149] hms_1.1.3 bit64_4.0.5
## [151] future_1.34.0 ggplot2_3.5.1
## [153] KEGGREST_1.45.1 statmod_1.5.0
## [155] highr_0.11 SummarizedExperiment_1.35.1
## [157] kernlab_0.9-33 igraph_2.0.3
## [159] memoise_2.0.1 affyio_1.75.0
## [161] bslib_0.8.0 sampling_2.10
## [163] bit_4.0.5