## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()
## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(dpi=25,fig.width=7)
## ----loaddata, message=FALSE--------------------------------------------------
library("pRolocdata")
data("thpLOPIT_unstimulated_rep1_mulvey2021")
data("thpLOPIT_unstimulated_rep3_mulvey2021")
data("thpLOPIT_lps_rep1_mulvey2021")
data("thpLOPIT_lps_rep3_mulvey2021")
## ----summarydata--------------------------------------------------------------
thpLOPIT_unstimulated_rep1_mulvey2021
thpLOPIT_lps_rep1_mulvey2021
## ----ldpkg, message=FALSE-----------------------------------------------------
library("bandle")
library("pheatmap")
library("viridis")
library("dplyr")
library("ggplot2")
## ----datadim------------------------------------------------------------------
dim(thpLOPIT_unstimulated_rep1_mulvey2021)
dim(thpLOPIT_unstimulated_rep3_mulvey2021)
dim(thpLOPIT_lps_rep1_mulvey2021)
dim(thpLOPIT_lps_rep3_mulvey2021)
## ----cmnprots-----------------------------------------------------------------
thplopit <- commonFeatureNames(c(thpLOPIT_unstimulated_rep1_mulvey2021, ## unstimulated rep
thpLOPIT_unstimulated_rep3_mulvey2021, ## unstimulated rep
thpLOPIT_lps_rep1_mulvey2021, ## 12h-LPS rep
thpLOPIT_lps_rep3_mulvey2021)) ## 12h-LPS rep
## ----listmsnsets--------------------------------------------------------------
thplopit
## ----exampledata, fig.height=10, fig.width=10---------------------------------
## create a character vector of title names for the plots
plot_id <- c("Unstimulated 1st rep", "Unstimulated 2nd rep",
"12h-LPS 1st rep", "12h-LPS 2nd rep")
## Let's set the stock colours of the classes to plot to be transparent
setStockcol(NULL)
setStockcol(paste0(getStockcol(), "90"))
## plot the data
par(mfrow = c(2,2))
for (i in seq(thplopit))
plot2D(thplopit[[i]], main = plot_id[i])
addLegend(thplopit[[4]], where = "topleft", cex = .75)
## ----fitgps-------------------------------------------------------------------
gpParams <- lapply(thplopit, function(x) fitGPmaternPC(x))
## ----lengthmrk----------------------------------------------------------------
(mrkCl <- getMarkerClasses(thplopit[[1]], fcol = "markers"))
## ----sethyppar----------------------------------------------------------------
K <- length(mrkCl)
pc_prior <- matrix(NA, ncol = 3, K)
pc_prior[seq.int(1:K), ] <- matrix(rep(c(1, 60, 100),
each = K), ncol = 3)
head(pc_prior)
## ----runhyppar----------------------------------------------------------------
gpParams <- lapply(thplopit,
function(x) fitGPmaternPC(x, hyppar = pc_prior))
## ----plotgps, fig.height=10, fig.width=8--------------------------------------
par(mfrow = c(4, 3))
plotGPmatern(thplopit[[1]], gpParams[[1]])
## ----setweightprior-----------------------------------------------------------
set.seed(1)
dirPrior = diag(rep(1, K)) + matrix(0.001, nrow = K, ncol = K)
predDirPrior <- prior_pred_dir(object = thplopit[[1]],
dirPrior = dirPrior,
q = 15)
## -----------------------------------------------------------------------------
predDirPrior$meannotAlloc
## -----------------------------------------------------------------------------
predDirPrior$tailnotAlloc
## ----fig.height=4, fig.width=6------------------------------------------------
hist(predDirPrior$priornotAlloc, col = getStockcol()[1])
## ----try, fig.height=4, fig.width=6-------------------------------------------
set.seed(1)
dirPrior = diag(rep(1, K)) + matrix(0.0005, nrow = K, ncol = K)
predDirPrior <- prior_pred_dir(object = thplopit[[1]],
dirPrior = dirPrior,
q = 15)
predDirPrior$meannotAlloc
predDirPrior$tailnotAlloc
hist(predDirPrior$priornotAlloc, col = getStockcol()[1])
## ----runbandle, message=FALSE-------------------------------------------------
control <- list(thplopit[[1]], thplopit[[2]])
treatment <- list(thplopit[[3]], thplopit[[4]])
bandleres <- bandle(objectCond1 = control,
objectCond2 = treatment,
numIter = 50, # usually 10,000
burnin = 5L, # usually 5,000
thin = 1L, # usually 20
gpParams = gpParams,
pcPrior = pc_prior,
numChains = 1, # usually >=4
dirPrior = dirPrior,
seed = 1)
## -----------------------------------------------------------------------------
bandleres
## ----processbandle------------------------------------------------------------
bandleres <- bandleProcess(bandleres)
## -----------------------------------------------------------------------------
summary(summaries(bandleres))
## -----------------------------------------------------------------------------
length(summaries(bandleres))
## ----getposteriors------------------------------------------------------------
pe1 <- posteriorEstimates(summaries(bandleres)[[1]])
pe2 <- posteriorEstimates(summaries(bandleres)[[2]])
head(pe1)
## ----barplotalloc, fig.width=9, fig.height=6----------------------------------
par(mfrow = c(1, 2), oma = c(6, 2, 2, 2))
barplot(table(pe1$bandle.allocation), col = getStockcol()[2],
las = 2, main = "Control: Protein allocation",
ylab = "Number of proteins")
barplot(table(pe2$bandle.allocation), col = getStockcol()[2],
las = 2, main = "Treatment: Protein allocation")
## ----meanbandleprob, fig.width=9, fig.height=6--------------------------------
par(mfrow = c(1, 2), oma = c(6, 2, 2, 2))
boxplot(pe1$bandle.probability ~ pe1$ bandle.allocation,
col = getStockcol()[2], xlab = "",
ylab = "BANDLE probability (mean)",
las = 2, main = "Control: Probability distribution\n by allocation class")
boxplot(pe2$bandle.probability ~ pe1$ bandle.allocation,
col = getStockcol()[2], xlab = "", ylab = "",
las = 2, main = "Treatment: Probability distribution\n by allocation class")
## ----predictlocation----------------------------------------------------------
## Add the bandle results to a MSnSet
xx <- bandlePredict(control,
treatment,
params = bandleres,
fcol = "markers")
res_0h <- xx[[1]]
res_12h <- xx[[2]]
## ----showappended, eval=FALSE-------------------------------------------------
# fvarLabels(res_0h[[1]])
# fvarLabels(res_0h[[2]])
#
# fvarLabels(res_12h[[1]])
# fvarLabels(res_12h[[2]])
## ----thresholddata------------------------------------------------------------
## threshold results using 1% FDR
res_0h[[1]] <- getPredictions(res_0h[[1]],
fcol = "bandle.allocation",
scol = "bandle.probability",
mcol = "markers",
t = .99)
res_12h[[1]] <- getPredictions(res_12h[[1]],
fcol = "bandle.allocation",
scol = "bandle.probability",
mcol = "markers",
t = .99)
## ----threshold2---------------------------------------------------------------
## add outlier probability
fData(res_0h[[1]])$bandle.outlier.t <- 1 - fData(res_0h[[1]])$bandle.outlier
fData(res_12h[[1]])$bandle.outlier.t <- 1 - fData(res_12h[[1]])$bandle.outlier
## threshold again, now on the outlier probability
res_0h[[1]] <- getPredictions(res_0h[[1]],
fcol = "bandle.allocation.pred",
scol = "bandle.outlier.t",
mcol = "markers",
t = .99)
res_12h[[1]] <- getPredictions(res_12h[[1]],
fcol = "bandle.allocation.pred",
scol = "bandle.outlier.t",
mcol = "markers",
t = .99)
## ----appendtosecond-----------------------------------------------------------
## Add results to second replicate at 0h
res_alloc_0hr <- fData(res_0h[[1]])$bandle.allocation.pred.pred
fData(res_0h[[2]])$bandle.allocation.pred.pred <- res_alloc_0hr
## Add results to second replicate at 12h
res_alloc_12hr <- fData(res_12h[[1]])$bandle.allocation.pred.pred
fData(res_12h[[2]])$bandle.allocation.pred.pred <- res_alloc_12hr
## ----plotmyres, fig.height=14, fig.width=5------------------------------------
par(mfrow = c(5, 2))
plot2D(res_0h[[1]], main = "Unstimulated - replicate 1 \n subcellular markers",
fcol = "markers")
plot2D(res_0h[[1]], main = "Unstimulated - replicate 1 \nprotein allocations (1% FDR)",
fcol = "bandle.allocation.pred.pred")
plot2D(res_0h[[2]], main = "Unstimulated - replicate 2 \nsubcellular markers",
fcol = "markers")
plot2D(res_0h[[2]], main = "Unstimulated - replicate 2 \nprotein allocations (1% FDR)",
fcol = "bandle.allocation.pred.pred")
plot2D(res_0h[[1]], main = "12h LPS - replicate 1 \nsubcellular markers",
fcol = "markers")
plot2D(res_0h[[1]], main = "12h LPS - replicate 1 \nprotein allocations (1% FDR)",
fcol = "bandle.allocation.pred.pred")
plot2D(res_0h[[2]], main = "12h LPS - replicate 2 \nsubcellular markers",
fcol = "markers")
plot2D(res_0h[[2]], main = "12h LPS - replicate 2 \nprotein allocations (1% FDR)",
fcol = "bandle.allocation.pred.pred")
plot(NULL, xaxt='n',yaxt='n',bty='n',ylab='',xlab='', xlim=0:1, ylim=0:1)
addLegend(res_0h[[1]], where = "topleft", cex = .8)
## ----numtransloc--------------------------------------------------------------
diffloc_probs <- pe1$bandle.differential.localisation
## ----cutoffres----------------------------------------------------------------
length(which(diffloc_probs[order(diffloc_probs, decreasing = TRUE)] > 0.95))
## ----extractdifflocp, fig.height=4, fig.width=6-------------------------------
plot(diffloc_probs[order(diffloc_probs, decreasing = TRUE)],
col = getStockcol()[2], pch = 19, ylab = "Probability",
xlab = "Rank", main = "Differential localisation rank plot")
## ----diffloc_binom------------------------------------------------------------
set.seed(1)
bin_t <- binomialDiffLocProb(params = bandleres, top = 500,
nsample = 500, decreasing = TRUE)
## ----get_pe-------------------------------------------------------------------
qt <- apply(bin_t, 1, function(x) quantile(x, .025))
## ----candidates---------------------------------------------------------------
candidates <- names(which(qt > 0.95))
head(candidates)
## ----chkallsame---------------------------------------------------------------
all(featureNames(res_0h[[1]]) == featureNames(res_0h[[2]]))
all(featureNames(res_0h[[1]]) == featureNames(res_12h[[1]]))
all(featureNames(res_12h[[1]]) == featureNames(res_12h[[2]]))
## ----addtomsn-----------------------------------------------------------------
dl.estimate <- qt[candidates]
fn <- featureNames(control[[1]])
cmn <- fn %in% names(dl.estimate)
## Add results to the 0h time-point (control)
for (i in seq(res_0h)) {
## create column called "dl.estimate" in the data
mcol <- "dl.estimate"
fData(res_0h[[i]])[, mcol] <- NA
fData(res_0h[[i]])[cmn, mcol] <- dl.estimate
## create column called "dl.candidate" in the data
mcol <- "dl.candidate"
fData(res_0h[[i]])[, mcol] <- "unknown"
fData(res_0h[[i]])[cmn, mcol] <- "DL candidate"
}
## Add results to the 12h time-point (treatment)
for (i in seq(res_12h)) {
## create column called "dl.estimate" in the data
mcol <- "dl.estimate"
fData(res_12h[[i]])[, mcol] <- NA
fData(res_12h[[i]])[cmn, mcol] <- dl.estimate
## create column called "dl.candidate" in the data
mcol <- "dl.candidate"
fData(res_12h[[i]])[, mcol] <- "unknown"
fData(res_12h[[i]])[cmn, mcol] <- "DL candidate"
}
## ----alluvial, warning=FALSE, message=FALSE-----------------------------------
msnset_cands <- list(res_0h[[1]][candidates, ],
res_12h[[1]][candidates, ])
## ----dataframeres-------------------------------------------------------------
# construct data.frame
df_cands <- data.frame(
fData(msnset_cands[[1]])[, c("bandle.differential.localisation",
"dl.estimate",
"bandle.allocation.pred.pred")],
fData(msnset_cands[[2]])[, "bandle.allocation.pred.pred"])
colnames(df_cands) <- c("bandle.differential.localisation", "dl.estimate",
"0hr_location", "12h_location")
# order by highest differential localisation estimate
df_cands <- df_cands %>% arrange(desc(bandle.differential.localisation))
head(df_cands)
## ----plotres, fig.height=8, fig.width=8---------------------------------------
## set colours for organelles and unknown
cols <- c(getStockcol()[seq(mrkCl)], "grey")
names(cols) <- c(mrkCl, "unknown")
## plot
alluvial <- plotTranslocations(msnset_cands,
fcol = "bandle.allocation.pred.pred",
col = cols)
alluvial + ggtitle("Differential localisations following 12h-LPS stimulation")
## ----tbllocs------------------------------------------------------------------
(tbl <- plotTable(msnset_cands, fcol = "bandle.allocation.pred.pred"))
## ----plotlysos, fig.height=8, fig.width=8-------------------------------------
secretory_prots <- unlist(lapply(msnset_cands, function(z)
c(which(fData(z)$bandle.allocation.pred.pred == "Golgi Apparatus"),
which(fData(z)$bandle.allocation.pred.pred == "Endoplasmic Reticulum"),
which(fData(z)$bandle.allocation.pred.pred == "Plasma Membrane"),
which(fData(z)$bandle.allocation.pred.pred == "Lysosome"))))
secretory_prots <- unique(secretory_prots)
msnset_secret <- list(msnset_cands[[1]][secretory_prots, ],
msnset_cands[[2]][secretory_prots, ])
secretory_alluvial <- plotTranslocations(msnset_secret,
fcol = "bandle.allocation.pred.pred",
col = cols)
secretory_alluvial + ggtitle("Movements of secretory proteins")
## ----plotdfprof---------------------------------------------------------------
head(df_cands)
## ----protprof, fig.height=7, fig.width=8--------------------------------------
par(mfrow = c(2, 1))
## plot lysosomal profiles
lyso <- which(fData(res_0h[[1]])$bandle.allocation.pred.pred == "Lysosome")
plotDist(res_0h[[1]][lyso], pcol = cols["Lysosome"], alpha = 0.04)
matlines(exprs(res_0h[[1]])["B2RUZ4", ], col = cols["Lysosome"], lwd = 3)
title("Protein SMIM1 (B2RUZ4) at 0hr (control)")
## plot plasma membrane profiles
pm <- which(fData(res_12h[[1]])$bandle.allocation.pred.pred == "Plasma Membrane")
plotDist(res_12h[[1]][pm], pcol = cols["Plasma Membrane"], alpha = 0.04)
matlines(exprs(res_12h[[1]])["B2RUZ4", ], col = cols["Plasma Membrane"], lwd = 3)
title("Protein SMIM1 (B2RUZ4) at 12hr-LPS (treatment)")
## ----plotpcacands, fig.height=6, fig.width=9----------------------------------
par(mfrow = c(1, 2))
plot2D(res_0h[[1]], fcol = "bandle.allocation.pred.pred",
main = "Unstimulated - replicate 1 \n predicted location")
highlightOnPlot(res_0h[[1]], foi = "B2RUZ4")
plot2D(res_12h[[1]], fcol = "bandle.allocation.pred.pred",
main = "12h-LPS - replicate 1 \n predicted location")
highlightOnPlot(res_12h[[1]], foi = "B2RUZ4")
## ----sessionInfo--------------------------------------------------------------
sessionInfo()