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