### R code from vignette source 'vignettes/pRoloc/inst/doc/pRoloc-tutorial.Rnw' ################################################### ### code chunk number 1: env ################################################### library("knitr") opts_chunk$set(fig.align = 'center', fig.show = 'hold', par = TRUE, prompt = TRUE, eval = TRUE, stop_on_error = 1L, comment = NA) options(replace.assign = TRUE, width = 55) suppressPackageStartupMessages(library("MSnbase")) suppressWarnings(suppressPackageStartupMessages(library("pRoloc"))) suppressPackageStartupMessages(library("pRolocdata")) suppressPackageStartupMessages(library("class")) set.seed(1) ################################################### ### code chunk number 2: libraries ################################################### library("MSnbase") library("pRoloc") library("pRolocdata") ################################################### ### code chunk number 3: setcols ################################################### setStockcol(paste0(getStockcol(), 70)) ################################################### ### code chunk number 4: readCsvData0 ################################################### ## The original data for replicate 1, available ## from the pRolocdata package f0 <- dir(system.file("extdata", package = "pRolocdata"), full.names = TRUE, pattern = "pr800866n_si_004-rep1.csv") csv <- read.csv(f0) ################################################### ### code chunk number 5: showOrgCsv ################################################### head(csv, n=3) ################################################### ### code chunk number 6: readCsvData1 ################################################### ## The quantitation data, from the original data f1 <- dir(system.file("extdata", package = "pRolocdata"), full.names = TRUE, pattern = "exprsFile.csv") exprsCsv <- read.csv(f1) ## Feature meta-data, from the original data f2 <- dir(system.file("extdata", package = "pRolocdata"), full.names = TRUE, pattern = "fdataFile.csv") fdataCsv <- read.csv(f2) ## Sample meta-data, a new file f3 <- dir(system.file("extdata", package = "pRolocdata"), full.names = TRUE, pattern = "pdataFile.csv") pdataCsv <- read.csv(f3) ################################################### ### code chunk number 7: showExprsFile ################################################### head(exprsCsv, n=3) ################################################### ### code chunk number 8: showFdFile ################################################### head(fdataCsv, n=3) ################################################### ### code chunk number 9: showPdFile ################################################### pdataCsv ################################################### ### code chunk number 10: makeMSnSet ################################################### tan2009r1 <- readMSnSet(exprsFile = f1, featureDataFile = f2, phenoDataFile = f3, sep = ",") tan2009r1 ################################################### ### code chunk number 11: showSubset ################################################### smallTan <- tan2009r1[1:5, 1:2] dim(smallTan) exprs(smallTan) ################################################### ### code chunk number 12: rmtan ################################################### ## remove manullay constructred data rm(tan2009r1) data(tan2009r1) ################################################### ### code chunk number 13: loadTan1 ################################################### data(tan2009r1) ################################################### ### code chunk number 14: lookAtTan ################################################### getMarkers(tan2009r1, fcol = "markers.orig") getMarkers(tan2009r1, fcol = "PLSDA") ################################################### ### code chunk number 15: markers ################################################### pRolocmarkers() head(pRolocmarkers("dmel")) table(pRolocmarkers("dmel")) ################################################### ### code chunk number 16: realtiveQuants ################################################### summary(rowSums(exprs(tan2009r1))) ################################################### ### code chunk number 17: norm ################################################### ## create a small illustrative test data data(itraqdata) itraqdata <- quantify(itraqdata, method = "trap", reporters = iTRAQ4, verbose = FALSE, parallel = FALSE) ## the quantification data head(exprs(itraqdata), n = 3) summary(rowSums(exprs(itraqdata))) ## normalising to the sum of feature intensitites itraqnorm <- normalise(itraqdata, method = "sum") processingData(itraqnorm) head(exprs(itraqnorm), n = 3) summary(rowSums(exprs(itraqnorm))) ################################################### ### code chunk number 18: featurenames ################################################### head(featureNames(tan2009r1)) ################################################### ### code chunk number 19: showplotdist (eval = FALSE) ################################################### ## ## indices of the mito markers ## j <- which(fData(tan2009r1)$markers.orig == "mitochondrion") ## ## indices of all proteins assigned to the mito ## i <- which(fData(tan2009r1)$PLSDA == "mitochondrion") ## plotDist(tan2009r1[i, ], ## markers = featureNames(tan2009r1)[j]) ################################################### ### code chunk number 20: plotdist1 ################################################### par(mfrow = c(2,2), mar = c(4, 4, 2, 0.1)) cls <- getStockcol()[1:4] plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "mitochondrion"), ], markers = featureNames(tan2009r1) [which(fData(tan2009r1)$markers.orig == "mitochondrion")], mcol = cls[1]) title(main = "mitochondrion") plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "ER/Golgi"), ], markers = featureNames(tan2009r1) [which(fData(tan2009r1)$markers.orig == "ER")], mcol = cls[2]) title(main = "ER") plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "ER/Golgi"), ], markers = featureNames(tan2009r1) [which(fData(tan2009r1)$markers.orig == "Golgi")], mcol = cls[3]) title(main = "Golgi") plotDist(tan2009r1[which(fData(tan2009r1)$PLSDA == "PM"), ], markers = featureNames(tan2009r1) [which(fData(tan2009r1)$markers.orig == "PM")], mcol = cls[4]) title(main = "PM") ################################################### ### code chunk number 21: plot2d ################################################### plot2D(tan2009r1, fcol = "PLSDA", fpch = "markers.orig") addLegend(tan2009r1, fcol = "PLSDA", where = "bottomright", bty = "n", cex = .7) ################################################### ### code chunk number 22: foi0 ################################################### foi1 <- FeaturesOfInterest(description = "Feats of interest 1", fnames = featureNames(tan2009r1[1:10])) description(foi1) foi(foi1) ################################################### ### code chunk number 23: pRoloc-tutorial.Rnw:656-660 ################################################### foi2 <- FeaturesOfInterest(description = "Feats of interest 2", fnames = featureNames(tan2009r1[880:888])) foic <- FoICollection(list(foi1, foi2)) foic ################################################### ### code chunk number 24: plotfoi ################################################### plot2D(tan2009r1, fcol = "PLSDA") addLegend(tan2009r1, fcol = "PLSDA", where = "bottomright", bty = "n", cex = .7) highlightOnPlot(tan2009r1, foi1, col = "black", lwd = 2) highlightOnPlot(tan2009r1, foi2, col = "purple", lwd = 2) legend("topright", c("FoI 1", "FoI 2"), bty = "n", col = c("black", "purple"), pch = 1) ################################################### ### code chunk number 25: guiinstall (eval = FALSE) ################################################### ## library("BiocInstaller") ## biocLite("pRolocGUI") ## library("pRolocGUI") ## pRolocVis(tan2009r1) ################################################### ### code chunk number 26: plot2dnull ################################################### plot2D(tan2009r1, fcol = NULL) ################################################### ### code chunk number 27: plotKmeans ################################################### kcl <- MLearn( ~ ., tan2009r1, kmeansI, centers=5) plot(kcl, exprs(tan2009r1)) ################################################### ### code chunk number 28: plotHclust ################################################### hcl <- MLearn( ~ ., tan2009r1, hclustI(distFun = dist, cutParm = list(k = 5))) plot(hcl, labels = FALSE) ################################################### ### code chunk number 29: plotPam ################################################### pcl <- MLearn( ~ ., tan2009r1, pamI(dist), k = 5) plot(pcl, data = exprs(tan2009r1)) ################################################### ### code chunk number 30: svmParamOptim (eval = FALSE) ################################################### ## params <- svmOptimisation(tan2009r1, fcol = "markers.orig", ## times = 10, xval = 5, verbose = FALSE) ################################################### ### code chunk number 31: loadParams ################################################### fn <- dir(system.file("extdata", package = "pRoloc"), full.names = TRUE, pattern = "params.rda") load(fn) rm(fn) ################################################### ### code chunk number 32: showParams ################################################### params ################################################### ### code chunk number 33: params ################################################### plot(params) levelPlot(params) ################################################### ### code chunk number 34: f1count ################################################### f1Count(params) getParams(params) ################################################### ### code chunk number 35: svmRes0 (eval = FALSE) ################################################### ## ## manual setting of parameters ## svmres <- svmClassification(tan2009r1, fcol = "markers.orig", ## sigma = 1, cost = 1) ################################################### ### code chunk number 36: svmRes ################################################### ## using default best parameters svmres <- svmClassification(tan2009r1, fcol = "markers.orig", assessRes = params) processingData(svmres) tail(fvarLabels(svmres), 4) ################################################### ### code chunk number 37: getPredictions ################################################### p1 <- getPredictions(svmres, fcol = "svm") minprob <- median(fData(svmres)$svm.scores) p2 <- getPredictions(svmres, fcol = "svm", t = minprob) table(p1, p2) ################################################### ### code chunk number 38: predscoresPlot ################################################### boxplot(svm.scores ~ svm, data = fData(svmres), ylab = "SVM scores") abline(h = minprob, lwd = 2, lty = 2) ################################################### ### code chunk number 39: medscores1 ################################################### ## including marker scores (ts1 <- tapply(fData(svmres)$svm.scores, fData(svmres)$svm, median)) ################################################### ### code chunk number 40: medscores2 ################################################### ## ignoring markers scores (i.e. scores == 1) (ts2 <- tapply(fData(svmres)$svm.scores, fData(svmres)$svm, function(x) { ## assuming scores of 1 are markers scr <- median(x[x != 1]) ## in case no proteins were assigned to an organelle, ## scr would be NA. Setting these cases to 1. ifelse(is.na(scr), 1, scr) })) ################################################### ### code chunk number 41: medscorepreds ################################################### getPredictions(svmres, fcol = "svm", t = ts2) ################################################### ### code chunk number 42: svmresfig ################################################### ptsze <- exp(fData(svmres)$svm.scores) - 1 plot2D(svmres, fcol = "svm", fpch = "markers.orig", cex = ptsze) addLegend(svmres, fcol = "svm", where = "bottomright", bty = "n", cex = .5) ################################################### ### code chunk number 43: runPhenoDisco (eval = FALSE) ################################################### ## pdres <- phenoDisco(tan2009r1, GS = 10, times = 100, fcol = "PLSDA") ################################################### ### code chunk number 44: loadpdres ################################################### fn <- dir(system.file("extdata", package = "pRoloc"), full.names = TRUE, pattern = "pdres.rda") load(fn) rm(fn) ################################################### ### code chunk number 45: phenoDiscoFvar ################################################### processingData(pdres) tail(fvarLabels(pdres), 3) ################################################### ### code chunk number 46: pdresfig ################################################### plot2D(pdres, fcol = "pd") addLegend(pdres, fcol = "pd", ncol = 2, where = "bottomright", bty = "n", cex = .5) ################################################### ### code chunk number 47: pdmarkers ################################################### getMarkers(tan2009r1, fcol = "pd.markers") ################################################### ### code chunk number 48: weights ################################################### w <- table(fData(tan2009r1)[, "pd.markers"]) (w <- 1/w[names(w) != "unknown"]) ################################################### ### code chunk number 49: pdsvmParams (eval = FALSE) ################################################### ## params2 <- svmOptimisation(tan2009r1, fcol = "pd.markers", ## times = 10, xval = 5, ## class.weights = w, ## verbose = FALSE) ################################################### ### code chunk number 50: loadParams2 ################################################### fn <- dir(system.file("extdata", package = "pRoloc"), full.names = TRUE, pattern = "params2.rda") load(fn) rm(fn) ################################################### ### code chunk number 51: pdsvm ################################################### tan2009r1 <- svmClassification(tan2009r1, params2, class.weights = w, fcol = "pd.markers") ################################################### ### code chunk number 52: pdres2fig ################################################### ptsze <- exp(fData(tan2009r1)$svm.scores) - 1 plot2D(tan2009r1, fcol = "svm", cex = ptsze) addLegend(tan2009r1, fcol = "svm", where = "bottomright", ncol = 2, bty = "n", cex = .5)