## ----style, echo = FALSE, results = 'asis'------------------------------------ BiocStyle::markdown() ## ----setup, include=FALSE----------------------------------------------------- knitr::opts_chunk$set(dpi=40,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") library("gridExtra") ## ----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 replicate 1", "Unstimulated replicate 2", "12h-LPS replicate 1", "12h-LPS replicate 2") ## 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) ## ----lengthmrk---------------------------------------------------------------- (mrkCl <- getMarkerClasses(thplopit[[1]], fcol = "markers")) ## ----------------------------------------------------------------------------- K <- length(mrkCl) ## ----sethyppar---------------------------------------------------------------- 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=10------------------------------------- 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=3, fig.width=4------------------------------------------------ hist(predDirPrior$priornotAlloc, col = getStockcol()[1]) ## ----try, fig.height=3, fig.width=4------------------------------------------- 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]]) params <- bandle(objectCond1 = control, objectCond2 = treatment, numIter = 10, # usually 10,000 burnin = 5L, # usually 5,000 thin = 1L, # usually 20 gpParams = gpParams, pcPrior = pc_prior, numChains = 4, # usually >=4 dirPrior = dirPrior, seed = 1) ## ----------------------------------------------------------------------------- params ## ----gelman------------------------------------------------------------------- calculateGelman(params) ## ----tracedens, fig.height=8, fig.width=8------------------------------------- plotOutliers(params) ## ----pressure, echo=FALSE, fig.cap="", out.width = "100%"--------------------- knitr::include_graphics("figs/traceplot.png") ## ----rmchains----------------------------------------------------------------- params_converged <- params[-c(1, 4)] ## ----viewnew------------------------------------------------------------------ params_converged ## ----processbandle------------------------------------------------------------ params_converged <- bandleProcess(params_converged) ## ----------------------------------------------------------------------------- bandle_out <- summaries(params_converged) ## ----------------------------------------------------------------------------- length(bandle_out) class(bandle_out[[1]]) ## ----results = "hide"--------------------------------------------------------- bandle_out[[1]]@posteriorEstimates bandle_out[[1]]@diagnostics bandle_out[[1]]@bandle.joint ## ----predictlocation---------------------------------------------------------- ## Add the bandle results to a MSnSet xx <- bandlePredict(control, treatment, params = params_converged, 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) ## ----------------------------------------------------------------------------- ## Remove the markers from the MSnSet res0hr_unknowns <- unknownMSnSet(res_0h[[1]], fcol = "markers") res12h_unknowns <- unknownMSnSet(res_12h[[1]], fcol = "markers") ## ----------------------------------------------------------------------------- res1 <- fData(res0hr_unknowns)$bandle.allocation.pred.pred res2 <- fData(res12h_unknowns)$bandle.allocation.pred.pred res1_tbl <- table(res1) res2_tbl <- table(res2) ## ----fig.height=4, fig.width=8------------------------------------------------ par(mfrow = c(1, 2)) barplot(res1_tbl, las = 2, main = "Predicted location: 0hr", ylab = "Number of proteins") barplot(res2_tbl, las = 2, main = "Predicted location: 12hr", ylab = "Number of proteins") ## ----allocspost, fig.height=4, fig.width=8------------------------------------ pe1 <- fData(res0hr_unknowns)$bandle.probability pe2 <- fData(res12h_unknowns)$bandle.probability par(mfrow = c(1, 2)) boxplot(pe1 ~ res1, las = 2, main = "Posterior: control", ylab = "Probability") boxplot(pe2 ~ res2, las = 2, main = "Posterior: treatment", ylab = "Probability") ## ----------------------------------------------------------------------------- res0hr_mixed <- unknownMSnSet(res0hr_unknowns, fcol = "bandle.allocation.pred.pred") res12hr_mixed <- unknownMSnSet(res12h_unknowns, fcol = "bandle.allocation.pred.pred") ## ----------------------------------------------------------------------------- nrow(res0hr_mixed) nrow(res12hr_mixed) ## ----------------------------------------------------------------------------- fn1 <- featureNames(res0hr_mixed) fn2 <- featureNames(res12hr_mixed) ## ----plot_mixed, fig.height=10, fig.width=10---------------------------------- g <- vector("list", 9) for (i in 1:9) g[[i]] <- mcmc_plot_probs(params_converged, fn1[i], cond = 1) do.call(grid.arrange, g) ## ----plot_mixed2, fig.height=10, fig.width=10--------------------------------- g <- vector("list", 9) for (i in 1:9) g[[i]] <- mcmc_plot_probs(params_converged, fn1[i], cond = 2) do.call(grid.arrange, g) ## ----------------------------------------------------------------------------- head(fData(res0hr_mixed)$bandle.joint) ## ----heatmap_control---------------------------------------------------------- bjoint_0hr_mixed <- fData(res0hr_mixed)$bandle.joint pheatmap(bjoint_0hr_mixed, cluster_cols = FALSE, color = viridis(n = 25), show_rownames = FALSE, main = "Joint posteriors for unlabelled proteins at 0hr") ## ----heatmap_treatment-------------------------------------------------------- bjoint_12hr_mixed <- fData(res12hr_mixed)$bandle.joint pheatmap(bjoint_12hr_mixed, cluster_cols = FALSE, color = viridis(n = 25), show_rownames = FALSE, main = "Joint posteriors for unlabelled proteins at 12hr") ## ----------------------------------------------------------------------------- dl <- diffLocalisationProb(params_converged) head(dl) ## ----cutoffres---------------------------------------------------------------- length(which(dl[order(dl, decreasing = TRUE)] > 0.95)) ## ----extractdifflocp, fig.height=4, fig.width=6------------------------------- plot(dl[order(dl, decreasing = TRUE)], col = getStockcol()[2], pch = 19, ylab = "Probability", xlab = "Rank", main = "Differential localisation rank plot") ## ----------------------------------------------------------------------------- candidates <- names(dl) ## ----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", "bandle.allocation.pred.pred")], fData(msnset_cands[[2]])[, "bandle.allocation.pred.pred"]) colnames(df_cands) <- c("differential.localisation", "0hr_location", "12h_location") # order by highest differential localisation estimate df_cands <- df_cands %>% arrange(desc(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()