## ----style, echo = FALSE, results = 'asis'------------------------------------
BiocStyle::markdown()

## ----setup, include=FALSE-----------------------------------------------------
knitr::opts_chunk$set(dpi=25,fig.width=7)

## ----env, message = FALSE, warning = FALSE, echo = TRUE-----------------------
# options(width = 150)
library("bandle")
library("ggplot2")
library("pRoloc")
library("MSnbase")
library("pheatmap")
library("viridis")
library("dplyr")

setStockcol(NULL)
setStockcol(paste0(getStockcol(),"90")) # make colours transparent

## ----loaddata-----------------------------------------------------------------
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

## ----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")

## 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(10, 60, 250),
                                       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("dl.estimate", "bandle.allocation.pred.pred")],
    fData(msnset_cands[[2]])[, "bandle.allocation.pred.pred"])
colnames(df_cands) <- c("dl.estimate", "0hr_location", "12h_location")

# order by highest differential localisation estimate
df_cands <- df_cands %>% arrange(desc(dl.estimate))
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()

