## ----style, echo = FALSE, results = 'asis'------------------------------------ BiocStyle::markdown() ## ----environment, cache=FALSE, echo=FALSE------------------------------------- suppressPackageStartupMessages(library("ggplot2")) suppressPackageStartupMessages(library("MSnbase")) suppressPackageStartupMessages(library("zoo")) suppressPackageStartupMessages(require("Rdisop")) suppressPackageStartupMessages(require("pRolocdata")) suppressPackageStartupMessages(require("pRoloc")) suppressPackageStartupMessages(require("msdata")) library("grid") suppressPackageStartupMessages(library("BiocParallel")) ## ----include_forword, echo=FALSE, results="asis"------------------------------ cat(readLines("./Foreword.md"), sep = "\n") ## ----include_bugs, echo=FALSE, results="asis"--------------------------------- cat(readLines("./Bugs.md"), sep = "\n") ## ---- eval=FALSE-------------------------------------------------------------- # library("doParallel") # registerDoParallel(3) ## using 3 slave nodes # register(DoparParam(), default = TRUE) # # ## rest of script comes below ## ----readdata, echo=TRUE, cache=FALSE, tidy=FALSE----------------------------- file <- dir(system.file(package = "MSnbase", dir = "extdata"), full.names = TRUE, pattern = "mzXML$") rawdata <- readMSData(file, msLevel = 2, verbose = FALSE) ## ----writeMSData, cache = FALSE----------------------------------------------- writeMSData(rawdata, file = paste0(tempfile(), ".mzML"), copy = TRUE) ## ----MSnExp, cache=FALSE, echo=TRUE------------------------------------------- library("MSnbase") itraqdata head(fData(itraqdata)) ## ----experimentsize, echo=FALSE, cache=FALSE---------------------------------- sz <- sum(sapply(assayData(itraqdata), object.size)) + object.size(itraqdata) sz <- as.numeric(sz) sz <- round(sz/(1024^2), 2) ## ----Spectrum, cache=FALSE, echo=TRUE----------------------------------------- sp <- itraqdata[["X1"]] sp ## ----accessors, cache=FALSE, echo=TRUE---------------------------------------- peaksCount(sp) head(peaksCount(itraqdata)) rtime(sp) head(rtime(itraqdata)) ## ----ReporterIons------------------------------------------------------------- iTRAQ4 TMT16 ## ----Chromatogram------------------------------------------------------------- f <- c(system.file("microtofq/MM14.mzML", package = "msdata")) mtof <- readMSData(f, mode = "onDisk") mtof_tic <- chromatogram(mtof) mtof_tic ## ----Chromatogram-continue---------------------------------------------------- mtof_tic[1, 1] head(intensity(mtof_tic[1, 1])) head(rtime(mtof_tic[1, 1])) mz(mtof_tic[1, 1]) ## ----Chromatogram-bpc--------------------------------------------------------- mtof_bpc <- chromatogram(mtof, aggregationFun = "max") ## ----msmap, eval=FALSE-------------------------------------------------------- # ## downloads the data # library("rpx") # px1 <- PXDataset("PXD000001") # mzf <- pxget(px1, 7) # # ## reads the data # ms <- openMSfile(mzf) # hd <- header(ms) # # ## a set of spectra of interest: MS1 spectra eluted # ## between 30 and 35 minutes retention time # ms1 <- which(hd$msLevel == 1) # rtsel <- hd$retentionTime[ms1] / 60 > 30 & # hd$retentionTime[ms1] / 60 < 35 # # ## the map # M <- MSmap(ms, ms1[rtsel], 521, 523, .005, hd, zeroIsNA = TRUE) ## ----msmaplaod, echo=FALSE---------------------------------------------------- mrda <- dir(system.file(package = "MSnbase", dir = "extdata"), full.names = TRUE, pattern = "M.rda$") mrda2 <- dir(system.file(package = "MSnbase", dir = "extdata"), full.names = TRUE, pattern = "M2.rda$") load(mrda) load(mrda2) ## ----msmaphow----------------------------------------------------------------- M ## ----mapheat, fig.cap = "Heat map of a chunk of the MS data."----------------- plot(M, aspect = 1, allTicks = FALSE) ## ----map3d, fig.cap = "3 dimensional represention of MS map data."------------ plot3D(M) ## ----msmap2, eval=FALSE------------------------------------------------------- # i <- ms1[which(rtsel)][1] # j <- ms1[which(rtsel)][2] # M2 <- MSmap(ms, i:j, 100, 1000, 1, hd) ## ----m2----------------------------------------------------------------------- M2 ## ----map3d2, fig.cap = "3 dimensional represention of MS map data. MS1 and MS2 spectra are coloured in blue and magenta respectively."---- plot3D(M2) ## ----spectrumPlot, fig.keep='high', fig.cap = "Raw MS2 spectrum with details about reporter ions."---- plot(sp, reporters = iTRAQ4, full = TRUE) ## ----bsaSelect, eval=TRUE, echo=FALSE----------------------------------------- ## bsasel <- fData(itraqdata)$ProteinAccession == "BSA" bsasel <- 1:3 ## ----subset, echo=TRUE-------------------------------------------------------- sel <- fData(itraqdata)$ProteinAccession == "BSA" bsa <- itraqdata[sel] bsa as.character(fData(bsa)$ProteinAccession) ## ----msnexpPlot, fig.keep='last', fig.cap = "Experiment-wide raw MS2 spectra. The y-axes of the individual spectra are automatically rescaled to the same range. See section \\@ref(sec:norm) to rescale peaks identically."---- plot(bsa, reporters = iTRAQ4, full = FALSE) + theme_gray(8) ## ----chromPlot, fig.keep='high', fig.cap = "Base peak chromatogram."---------- plot(mtof_bpc) ## ----iddf0-------------------------------------------------------------------- library("msdata") f <- "TMT_Erwinia_1uLSike_Top10HCD_isol2_45stepped_60min_01-20141210.mzid" idf <- msdata::ident(full.names = TRUE, pattern = f) iddf <- readMzIdData(idf) str(iddf) ## ----dups1, echo=FALSE-------------------------------------------------------- iddf[grep("scan=5291", iddf[, "spectrumID"]), c("spectrumID", "sequence", "DatabaseAccess")] ## ----dups2, echo=FALSE-------------------------------------------------------- iddf[grep("scan=4936", iddf[, "spectrumID"]), c("spectrumID", "sequence", "modName", "modLocation")] ## ----ideda-------------------------------------------------------------------- table(iddf$isDecoy) table(iddf$chargeState) ## ----idvis-------------------------------------------------------------------- library("ggplot2") ggplot(data = iddf, aes(x = MS.GF.RawScore, colour = isDecoy)) + geom_density() + facet_wrap(~chargeState) ## ----------------------------------------------------------------------------- iddf <- filterIdentificationDataFrame(iddf) ## ----rediddf------------------------------------------------------------------ iddf2 <- reduce(iddf, key = "spectrumID") ## ----dups3, echo=FALSE-------------------------------------------------------- iddf2[grep("scan=5291", iddf2[, "spectrumID"]), c("spectrumID", "sequence", "DatabaseAccess")] ## ----dups4, echo=FALSE-------------------------------------------------------- iddf2[grep("scan=4936", iddf2[, "spectrumID"]), c("spectrumID", "sequence", "modName", "modLocation")] ## ----msnexpIdentification, echo=TRUE, cache=FALSE, tidy=FALSE----------------- ## find path to a mzXML file quantFile <- dir(system.file(package = "MSnbase", dir = "extdata"), full.name = TRUE, pattern = "mzXML$") ## find path to a mzIdentML file identFile <- dir(system.file(package = "MSnbase", dir = "extdata"), full.name = TRUE, pattern = "dummyiTRAQ.mzid") ## create basic MSnExp msexp <- readMSData(quantFile, verbose = FALSE) head(fData(msexp), n = 2) ## ----msnexpIdentification2, echo=TRUE, cache=FALSE, tidy=FALSE---------------- msexp <- addIdentificationData(msexp, id = identFile) head(fData(msexp), n = 2) ## ----msnexpIdentification3, echo=TRUE, cache=FALSE, tidy=FALSE---------------- idSummary(msexp) ## ----fragplot0---------------------------------------------------------------- itraqdata2 <- pickPeaks(itraqdata, verbose=FALSE) i <- 14 s <- as.character(fData(itraqdata2)[i, "PeptideSequence"]) ## ----fragplot, fig.cap = "Annotated MS2 spectrum."---------------------------- plot(itraqdata2[[i]], s, main = s) ## ----msnexpIdentification4, echo=TRUE, cache=FALSE, tidy=FALSE---------------- fData(msexp)$sequence msexp <- removeNoId(msexp) fData(msexp)$sequence idSummary(msexp) ## ----calculateFragments, echo=TRUE, cache=FALSE, tidy=FALSE------------------- calculateFragments("ACEK", type = c("a", "b", "c", "x", "y", "z")) ## ----msnexpcalculateFragments, echo=TRUE, cache=FALSE, tidy=FALSE------------- pepseq <- fData(msexp)$sequence[1] calculateFragments(pepseq, msexp[[1]], type=c("b", "y")) ## ----removePeaks, echo=TRUE, cache=FALSE-------------------------------------- experiment <- removePeaks(itraqdata, t = 400, verbose = FALSE) ionCount(itraqdata[["X55"]]) ionCount(experiment[["X55"]]) ## ----spectrum-clean-plot, echo=FALSE, fig.keep='high', fig.cap = "Same spectrum before (left) and after setting peaks <= 400 to 0."---- p1 <- plot(itraqdata[["X55"]], full = TRUE, plot = FALSE) + theme_gray(5) p2 <- plot(experiment[["X55"]], full = TRUE, plot = FALSE) + theme_gray(5) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) print(p1,vp=vplayout(1,1)) print(p2,vp=vplayout(1,2)) ## ----clean, echo=TRUE, cache=FALSE-------------------------------------------- peaksCount(itraqdata[["X55"]]) peaksCount(experiment[["X55"]]) experiment <- clean(experiment, verbose = FALSE) peaksCount(experiment[["X55"]]) ## ----preprosp, cache=FALSE, echo=FALSE---------------------------------------- int <- c(0,1,1,3,1,1,0,0,0,1,3,7,3,1,0) mz <- c(113.9,114.0,114.05,114.1,114.15,114.2,114.25, 114.3,114.35,114.4,114.42,114.48,114.5,114.55,114.6) ppsp <- new("Spectrum2",intensity=int,mz=mz,centroided=FALSE) p1 <- plot(ppsp, full = TRUE, plot = FALSE) + theme_gray(5) + geom_point(size=3,alpha=I(1/3)) + geom_hline(yintercept=3,linetype=2) + ggtitle("Original spectrum") p2 <- plot(removePeaks(ppsp,t=3), full=TRUE, plot = FALSE) + theme_gray(5) + geom_point(size=3,alpha=I(1/3)) + geom_hline(yintercept=3,linetype=2) + ggtitle("Peaks < 3 removed") p3 <- plot(clean(removePeaks(ppsp,t=3)), full = TRUE, plot = FALSE) + theme_gray(5) + geom_point(size=3,alpha=I(1/3)) + geom_hline(yintercept=3,linetype=2) + ggtitle("Peaks < 3 removed and cleaned") ## ----preprocPlot, echo=FALSE, fig.cap = "This figure illustrated the effect or the `removePeaks` and `clean` methods. The left-most spectrum displays two peaks, of max height 3 and 7 respectively. The middle spectrum shows the result of calling `removePeaks` with argument `t=3`, which sets all data points of the first peak, whose maximum height is smaller or equal to `t` to 0. The second peak is unaffected. Calling `clean` after `removePeaks` effectively deletes successive 0 intensities from the spectrum, as shown on the right plot."---- grid.newpage() pushViewport(viewport(layout = grid.layout(3, 1))) print(p1, vp=vplayout(1, 1)) print(p2, vp=vplayout(2, 1)) print(p3, vp=vplayout(3, 1)) ## ----trimMz, echo=TRUE, cache=FALSE------------------------------------------- range(mz(itraqdata[["X55"]])) experiment <- filterMz(experiment, mzlim = c(112,120)) range(mz(experiment[["X55"]])) experiment ## ----reporters, echo=TRUE, cache=FALSE---------------------------------------- mz(iTRAQ4) width(iTRAQ4) ## ----simplesp, cache=FALSE, echo=FALSE, fig.keep='none'----------------------- int <- c(0, 1, 1, 3, 1, 1, 0) mz <- c(113.9, 114.0, 114.05, 114.1, 114.15, 114.2, 114.25) ssp <- new("Spectrum2", intensity = int, mz = mz, centroided = FALSE) p <- plot(ssp, full = TRUE, plot = FALSE) p <- p + theme_gray(5) ## ----quantitationPlot, echo=FALSE, fig.cap = "The different quantitation methods. See text for details."---- grid.newpage() pushViewport(viewport(layout = grid.layout(2, 2))) vplayout <- function(x, y) viewport(layout.pos.row = x, layout.pos.col = y) print(p + ggtitle("Quantitation using 'sum'") + geom_point(size = 3, alpha = I(1/3), colour = "red"), vp = vplayout(1, 1)) print(p + ggtitle("Quantitation using 'max'") + geom_point(aes(x = 114.1, y = 3), alpha = I(1/18), colour = "red", size = 3), vp = vplayout(1, 2)) print(p + ggtitle("Trapezoidation and strict=FALSE") + geom_polygon(alpha = I(1/5), fill = "red"), vp = vplayout(2, 1)) print(p + ggtitle("Trapezoidation and strict=TRUE") + geom_polygon(aes(x = c(NA, 114.05, 114.05, 114.1, 114.15, 114.15, NA), y = c(NA, 0, 1, 3, 1, 0, NA)), fill = "red", alpha = I(1/5)), vp = vplayout(2,2)) ## ----quantify, echo=TRUE, cache=FALSE, tidy=FALSE----------------------------- qnt <- quantify(experiment, method = "trap", reporters = iTRAQ4, strict = FALSE, verbose = FALSE) qnt head(exprs(qnt)) ## ----filterNA, echo=TRUE------------------------------------------------------ table(is.na(qnt)) qnt <- filterNA(qnt, pNA = 0) sum(is.na(qnt)) ## ----removeNa, echo=TRUE, eval=FALSE------------------------------------------ # whichRow <- which(is.na((qnt))) %% nrow(qnt) # qnt <- qnt[-whichRow, ] ## ----readSRMData, echo = TRUE, warning = FALSE-------------------------------- fl <- proteomics(full.names = TRUE, pattern = "MRM") srm <- readSRMData(fl) srm ## ----srm-precursorMz---------------------------------------------------------- head(precursorMz(srm)) head(productMz(srm)) ## ----pheplus1, echo=TRUE, cache=FALSE----------------------------------------- library(Rdisop) ## Phenylalanine immonium ion Fim <- getMolecule("C8H10N") getMass(Fim) isotopes <- getIsotope(Fim) F1 <- isotopes[2, 2] F1 ## ----purityCorrect, echo=TRUE, cache=FALSE, tidy = FALSE---------------------- impurities <- matrix(c(0.929, 0.059, 0.002, 0.000, 0.020, 0.923, 0.056, 0.001, 0.000, 0.030, 0.924, 0.045, 0.000, 0.001, 0.040, 0.923), nrow = 4) qnt.crct <- purityCorrect(qnt, impurities) head(exprs(qnt)) head(exprs(qnt.crct)) ## ----impute------------------------------------------------------------------- ## an example MSnSet containing missing values data(naset) table(is.na(naset)) ## number of NAs per protein table(fData(naset)$nNA) x <- impute(naset, "min") processingData(x) table(is.na(x)) ## ----miximp, echo = FALSE, fig.cap = "Mixed imputation method. Black cells represent presence of quantitation values and light grey corresponds to missing data. The two groups of interest are depicted in green and blue along the heatmap columns. Two classes of proteins are annotated on the left: yellow are proteins with randomly occurring missing values (if any) while proteins in brown are candidates for non-random missing value imputation."---- x <- impute(naset, "zero") exprs(x)[exprs(x) != 0] <- 1 suppressPackageStartupMessages(library("gplots")) heatmap.2(exprs(x), col = c("lightgray", "black"), scale = "none", dendrogram = "none", trace = "none", keysize = 0.5, key = FALSE, RowSideColors = ifelse(fData(x)$randna, "orange", "brown"), ColSideColors = rep(c("steelblue", "darkolivegreen"), each = 8)) ## ----miximp2------------------------------------------------------------------ x <- impute(naset, method = "mixed", randna = fData(naset)$randna, mar = "knn", mnar = "min") x ## ----normalise, echo=TRUE, cache=FALSE---------------------------------------- qnt.max <- normalise(qnt, "max") qnt.sum <- normalise(qnt, "sum") qnt.quant <- normalise(qnt, "quantiles") qnt.qrob <- normalise(qnt, "quantiles.robust") qnt.vsn <- normalise(qnt, "vsn") ## ----normPlot, echo=FALSE, fig.cap = "Comparison of the normalisation *MSnSet* methods. Note that vsn also glog-transforms the intensities."---- .plot <- function(x, ttl = NULL) { boxplot(exprs(x), main = ifelse(is.null(ttl), processingData(x)@processing[2], ttl), cex.main = .8, cex.lab = .5, cex.axis = .5, cex = .8) grid() } oldmar <- par()$mar par(mfrow=c(3,2),mar=c(2.9,2.9,2.9,1)) .plot(qnt, ttl = "Non-normalised data") .plot(qnt.max, ttl = "Maximum") .plot(qnt.sum, ttl = "Sum") .plot(qnt.quant, ttl = "Quantile") .plot(qnt.qrob, ttl = "Robust quantile") .plot(qnt.vsn, ttl = "vsn") ## ----cvdata, echo=FALSE, cache=FALSE------------------------------------------ sd1 <- apply(log2(exprs(qnt))+10,1,sd) mn1 <- apply(log2(exprs(qnt))+10,1,mean) cv1 <- sd1/mn1 sd2 <- apply(exprs(qnt.vsn)+10,1,sd) mn2 <- apply(exprs(qnt.vsn)+10,1,mean) cv2 <- sd2/mn2 dfr <- rbind(data.frame(rank=order(mn1),cv=cv1,norm="raw"), data.frame(rank=order(mn2),cv=cv2,norm="vsn")) library("zoo") ## rmed1 <- rollapply(cv1,7,function(x) median(x,na.rm=TRUE)) ## rmed2 <- rollapply(cv2,7,function(x) median(x,na.rm=TRUE)) ## ## Calling directly rollapply.zoo to make it zoo_1.6-4 compatible. ## The above requires zoo >= 1.7-0, which is as of 15 March 2011 ## not yet available on CRAN (only on r-forge). rmed1 <- zoo:::rollapply.zoo(cv1,7,function(x) median(x,na.rm=TRUE)) rmed2 <- zoo:::rollapply.zoo(cv2,7,function(x) median(x,na.rm=TRUE)) dfr2 <- rbind(data.frame(x=seq(1,70,by=(70/length(rmed1))), y=rmed1,norm="raw"), data.frame(x=seq(1,70,by=(70/length(rmed1))), y=rmed2,norm="vsn")) p <- qplot(rank,cv,data=dfr,col=norm) + geom_line(data=dfr2,aes(x=x,y=y,colour=norm)) + theme_gray(7) ## ----cvPlot, echo=FALSE, fig.cap = "CV versus signal intensity comparison for log2 and vsn transformed data. Lines indicate running CV medians."---- print(p) ## ----prepareMsnsetNormPlot, cache=FALSE, echo=FALSE, keep.fig='none'---------- p <- plot(normalise(experiment[bsasel], "max"), reporters = iTRAQ4, full = FALSE, plot = FALSE) p <- p + theme_gray(7) ## ----msnexpNormPlot, echo=FALSE, fig.cap = "Experiment-wide normalised MS2 spectra. The y-axes of the individual spectra is now rescaled between 0 and 1 (highest peak), as opposed to figure \\@ref(fig:msnexpPlot)."---- print(p) ## ----makeGroups1,echo=FALSE,cache=FALSE--------------------------------------- gb <- fData(qnt)$ProteinAccession ## ----makeGroups2,echo=TRUE,cache=FALSE---------------------------------------- gb <- fData(qnt)$ProteinAccession table(gb) length(unique(gb)) ## ----combineFeatures, echo=TRUE, cache=FALSE---------------------------------- qnt2 <- combineFeatures(qnt, groupBy = gb, method = "median") qnt2 ## ----count-------------------------------------------------------------------- sc <- quantify(msexp, method = "count") ## lets modify out data for demonstration purposes fData(sc)$DatabaseAccess[1] <- fData(sc)$DatabaseAccess[2] fData(sc)$DatabaseAccess sc <- combineFeatures(sc, groupBy = fData(sc)$DatabaseAccess, method = "sum") exprs(sc) ## ----labelfree---------------------------------------------------------------- fData(msexp)[, c("DatabaseAccess", "nprot")] ## ----SIn---------------------------------------------------------------------- siquant <- quantify(msexp, method = "SIn") processingData(siquant) exprs(siquant) ## ----compms2plot, fig.cap = "Comparing two MS2 spectra."---------------------- centroided <- pickPeaks(itraqdata, verbose = FALSE) (k <- which(fData(centroided)[, "PeptideSequence"] == "TAGIQIVADDLTVTNPK")) mzk <- precursorMz(centroided)[k] zk <- precursorCharge(centroided)[k] mzk * zk plot(centroided[[k[1]]], centroided[[k[2]]]) ## ----compareSpectra, tidy=FALSE----------------------------------------------- compareSpectra(centroided[[2]], centroided[[3]], fun = "common") compareSpectra(centroided[[2]], centroided[[3]], fun = "cor") compareSpectra(centroided[[2]], centroided[[3]], fun = "dotproduct") ## ----msnexpcompareSpectra----------------------------------------------------- compmat <- compareSpectra(centroided, fun="cor") compmat[1:10, 1:5] ## ----dendo-------------------------------------------------------------------- plot(hclust(as.dist(compmat))) ## ----incompdiss, echo=TRUE, cache=FALSE, tidy=FALSE--------------------------- iTRAQ5 incompdiss <- quantify(itraqdata, method = "trap", reporters = iTRAQ5, strict = FALSE, verbose = FALSE) head(exprs(incompdiss)) ## ----incompdissPlot, echo=FALSE, warning=FALSE, fig.cap ="Boxplot and scatterplot comparing intensities of the 4 reporter ions (or their sum, on the right) and the incomplete dissociation specific peak."---- dfr <- cbind(expand.grid("feature" = featureNames(incompdiss), "reporters" = sub("iTRAQ5\\.", "", sampleNames(incompdiss)), stringsAsFactors=FALSE), "intensity" = as.vector(exprs(incompdiss))) repsum <- rowSums(exprs(incompdiss)[, 1:4]) dfr2 <- data.frame(iTRAQ1to4 = repsum, iTRAQ5 = exprs(incompdiss)[,5]) p1 <- ggplot(data = dfr, aes(x = reporters,y = log10(intensity))) + geom_boxplot() + theme_gray(6) p2 <- ggplot(data = dfr2, aes(x = log10(iTRAQ1to4), y = log10(iTRAQ5))) + geom_point(alpha = I(1/2)) + geom_abline(intercept = 0, slope = 1, linetype = "dotted") + stat_smooth(method = lm, se = FALSE) + xlab(label = expression(log[10]~sum~114~to~117)) + ylab(label = expression(log[10]~145)) + theme_gray(6) grid.newpage() pushViewport(viewport(layout = grid.layout(1, 2))) print(p1, vp = vplayout(1, 1)) print(p2, vp = vplayout(1, 2)) ## ----makeexp12, echo=TRUE, cache=FALSE, tidy = FALSE-------------------------- exp1 <- quantify(itraqdata, reporters = iTRAQ4, verbose = FALSE) sampleNames(exp1) centroided(rawdata) <- FALSE exp2 <- quantify(rawdata, reporters = iTRAQ4, verbose = FALSE) sampleNames(exp2) ## ----updateFnames, echo=TRUE, cache=FALSE------------------------------------- head(featureNames(exp1)) exp1 <- updateFeatureNames(exp1) head(featureNames(exp1)) head(featureNames(exp2)) exp2 <- updateFeatureNames(exp2) head(featureNames(exp2)) ## ----comb1, echo=TRUE, cache=FALSE-------------------------------------------- exp12 <- combine(exp1, exp2) dim(exp1) dim(exp2) dim(exp12) ## ----make2exps2, echo=TRUE, cache=FALSE, tidy=FALSE--------------------------- set.seed(1) i <- sample(length(itraqdata), 35) j <- sample(length(itraqdata), 35) exp1 <- quantify(itraqdata[i], reporters = iTRAQ4, verbose = FALSE) exp2 <- quantify(itraqdata[j], reporters = iTRAQ4, verbose = FALSE) exp1 <- droplevels(exp1) exp2 <- droplevels(exp2) table(featureNames(exp1) %in% featureNames(exp2)) exp1 <- combineFeatures(exp1, groupBy = fData(exp1)$ProteinAccession) exp2 <- combineFeatures(exp2, groupBy = fData(exp2)$ProteinAccession) head(featureNames(exp1)) head(featureNames(exp2)) ## ----renameSamples, echo=TRUE, cache=FALSE------------------------------------ sampleNames(exp1) exp1 <- updateSampleNames(exp1) sampleNames(exp1) sampleNames(exp1) <- c("Ctrl1", "Cond1", "Ctrl2", "Cond2") sampleNames(exp2) <- c("Ctrl3", "Cond3", "Ctrl4", "Cond4") ## ----fdatanames, echo=TRUE, cache=FALSE--------------------------------------- stopifnot(all(fvarLabels(exp1) == fvarLabels(exp2))) fData(exp1)["BSA", 1:4] fData(exp2)["BSA", 1:4] ## ----renameFvars, echo=TRUE, cache=FALSE-------------------------------------- exp1 <- updateFvarLabels(exp1) exp2 <- updateFvarLabels(exp2) head(fvarLabels(exp1)) head(fvarLabels(exp2)) ## ----combine, echo=TRUE, cache=FALSE------------------------------------------ exp12 <- combine(exp1, exp2) dim(exp12) pData(exp12) exprs(exp12)[25:28, ] exp12 ## ----split-------------------------------------------------------------------- data(dunkley2006) head(pData(dunkley2006)) split(dunkley2006, dunkley2006$replicate) ## or, defining the appropriate annotation variable name dun <- split(dunkley2006, "replicate") ## ----unsplit------------------------------------------------------------------ dun2 <- unsplit(dun, pData(dunkley2006)$replicate) compareMSnSets(dunkley2006, dun2) ## ----avg---------------------------------------------------------------------- library("pRolocdata") data(tan2009r1) data(tan2009r2) data(tan2009r3) msnl <- MSnSetList(list(tan2009r1, tan2009r2, tan2009r3)) avgtan <- averageMSnSet(msnl) head(exprs(avgtan)) head(fData(avgtan)$disp) head(fData(avgtan)$nNA) ## ----plot2Davg, fig.cap = "PCA plot of the averaged *MSnSet*. The point sizes are proportional to the dispersion of the protein quantitation across the averaged data."---- disp <- rowMax(fData(avgtan)$disp) disp[disp == 0] <- max(disp) range(disp) library("pRoloc") plot2D(avgtan, cex = 3 * disp) ## ----sessioninfo, echo=FALSE-------------------------------------------------- sessionInfo()