## ----style, echo = FALSE, results = 'asis'------------------------------------ BiocStyle::markdown() ## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(dpi=25,fig.width=7) ## ----eval=FALSE--------------------------------------------------------------- # if (!requireNamespace("BiocManager", quietly=TRUE)) # install.packages("BiocManager") # BiocManager::install("bandle") ## ----ldpkg, message=FALSE----------------------------------------------------- library("bandle") ## ----ldpkg2, message=FALSE---------------------------------------------------- library("pheatmap") library("viridis") ## ----loadpkgdat, message=FALSE, warning=FALSE, fig.width=7, fig.height=6------ 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) ## ----simdata, fig.width=8, fig.height=10-------------------------------------- set.seed(1) tansim <- sim_dynamic(object = tan2009r1, numRep = 6L, numDyn = 100L) ## ----plotsims----------------------------------------------------------------- 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])) ## ----------------------------------------------------------------------------- tansim$lopitrep[[1]] ## ----fitgps, fig.height=10, fig.width=8--------------------------------------- par(mfrow = c(4, 3)) gpParams <- lapply(tansim$lopitrep, function(x) fitGPmaternPC(x, hyppar = matrix(c(10, 60, 250), nrow = 1))) ## ----plotgps, fig.height=10, fig.width=8-------------------------------------- par(mfrow = c(4, 3)) plotGPmatern(tansim$lopitrep[[1]], params = gpParams[[1]]) ## ----sethyppar---------------------------------------------------------------- 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) ## ----runhyppar---------------------------------------------------------------- gpParams <- lapply(tansim$lopitrep, function(x) fitGPmaternPC(x, hyppar = pc_prior)) ## ----setweightprior----------------------------------------------------------- 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) ## ----------------------------------------------------------------------------- predDirPrior$meannotAlloc ## ----------------------------------------------------------------------------- predDirPrior$tailnotAlloc ## ----------------------------------------------------------------------------- hist(predDirPrior$priornotAlloc, col = getStockcol()[1]) ## ----runbandle, message=FALSE, warning=FALSE, error=FALSE, echo = TRUE, results = 'hide'---- 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) ## ----------------------------------------------------------------------------- bandleres ## ----processbandle1----------------------------------------------------------- summaries(bandleres) ## ----processbandle2----------------------------------------------------------- bandleres <- bandleProcess(bandleres) ## ----processbandle3----------------------------------------------------------- summary(summaries(bandleres)) ## ----res---------------------------------------------------------------------- res <- summaries(bandleres) length(res) ## ----seeslots----------------------------------------------------------------- str(res[[1]]) ## ----postest------------------------------------------------------------------ posteriorEstimates(res[[1]]) ## ----getposteriors------------------------------------------------------------ pe1 <- posteriorEstimates(res[[1]]) pe2 <- posteriorEstimates(res[[2]]) ## ----barplotalloc, fig.width=8, fig.height=5---------------------------------- 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") ## ----allocspost--------------------------------------------------------------- pe_alloc1 <- pe1$bandle.probability pe_alloc2 <- pe1$bandle.probability ## ----heatmap_control---------------------------------------------------------- bjoint_control <- bandleJoint(summaries(bandleres)[[1]]) pheatmap(bjoint_control, cluster_cols = FALSE, color = viridis(n = 25)) ## ----heatmap_treatment-------------------------------------------------------- bjoint_treatment <- bandleJoint(summaries(bandleres)[[2]]) pheatmap(bjoint_treatment, cluster_cols = FALSE, color = viridis(n = 25)) ## ----bandpred----------------------------------------------------------------- xx <- bandlePredict(control, treatment, params = bandleres, fcol = "markers") res_control <- xx[[1]] res_treatment <- xx[[2]] ## ----showlength--------------------------------------------------------------- length(res_control) length(res_treatment) ## ----viewdata----------------------------------------------------------------- fvarLabels(res_control[[1]]) ## ----fdata, eval=FALSE-------------------------------------------------------- # fData(res_control[[1]])$bandle.probability # fData(res_control[[1]])$bandle.allocation ## ----setthreshold------------------------------------------------------------- res_control[[1]] <- getPredictions(res_control[[1]], fcol = "bandle.allocation", scol = "bandle.probability", mcol = "markers", t = .99) res_treatment[[1]] <- getPredictions(res_treatment[[1]], fcol = "bandle.allocation", scol = "bandle.probability", mcol = "markers", t = .99) ## ----numtransloc-------------------------------------------------------------- diffloc_probs <- pe1$bandle.differential.localisation head(diffloc_probs, 50) length(which(diffloc_probs[order(diffloc_probs, decreasing = TRUE)] > 0.99)) ## ----extractdiffloc----------------------------------------------------------- plot(diffloc_probs[order(diffloc_probs, decreasing = TRUE)], col = getStockcol()[3], pch = 19, ylab = "Probability", xlab = "Rank", main = "Differential localisation rank plot") ## ----diffloc_boot------------------------------------------------------------- set.seed(1) boot_t <- bootstrapdiffLocprob(params = bandleres, top = 100, Bootsample = 5000, decreasing = TRUE) boxplot(t(boot_t), col = getStockcol()[5], las = 2, ylab = "Probability", ylim = c(0, 1), main = "Differential localisation \nprobability plot (top 100 proteins)") ## ----diffloc_binom------------------------------------------------------------ bin_t <- binomialDiffLocProb(params = bandleres, top = 100, nsample = 5000, decreasing = TRUE) boxplot(t(bin_t), col = getStockcol()[5], las = 2, ylab = "Probability", ylim = c(0, 1), main = "Differential localisation \nprobability plot (top 100 proteins)") ## ----get_pe------------------------------------------------------------------- # 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)) ## ----------------------------------------------------------------------------- EFDR(diffloc_probs, threshold = 0.90) ## ----alluvial, warning=FALSE, message=FALSE, fig.height=8, fig.width=7-------- plotTranslocations(bandleres) ## ----chord, warning=FALSE, message=FALSE, fig.height=7, fig.width=7----------- plotTranslocations(bandleres, type = "chord") ## ----plotafterthresh---------------------------------------------------------- ## 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, ] ## ----plottlres, warning=FALSE, fig.height=8, fig.width=7---------------------- ## 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) ## ----summaryfinal, warning=FALSE, message=FALSE------------------------------- (tbl <- plotTable(res, fcol = "bandle.allocation.pred")) ## ----bandleagain, message=FALSE, warning=FALSE, error=FALSE, echo = TRUE, results = 'hide',fig.height=4, fig.width=4---- 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) ## ----convergence, fig.height=8------------------------------------------------ par(mfrow = c(2, 2)) out <- plotConvergence(bandleres) ## ----sessionInfo-------------------------------------------------------------- sessionInfo()