## ----style, echo = FALSE, results = 'asis'------------------------------------ BiocStyle::markdown() knitr::opts_chunk$set(fig.wide = TRUE, fig.retina = 3, error=FALSE) ## ----HexSticker, echo=FALSE, out.width="50%", eval=TRUE----------------------- knitr::include_graphics("rawrr_logo.png") ## ----installRuntime, echo=TRUE------------------------------------------------ rawrr::installRawrrExe() ## ----tartareEH4547, warning=FALSE, message=FALSE, eval=TRUE------------------- # fetch via ExperimentHub library(ExperimentHub) eh <- ExperimentHub::ExperimentHub() EH4547 <- normalizePath(eh[["EH4547"]]) (rawfile <- paste0(EH4547, ".raw")) if (!file.exists(rawfile)){ file.copy(EH4547, rawfile) } ## ----check, eval=TRUE, echo=FALSE--------------------------------------------- stopifnot(file.exists(rawfile)) ## ----readFileHeader, eval=TRUE------------------------------------------------ H <- rawrr::readFileHeader(rawfile = rawfile) ## ----------------------------------------------------------------------------- try({ i <- c(9594, 11113, 11884, 12788, 12677, 13204, 13868, 14551, 16136, 17193, 17612) S <- rawrr::readSpectrum(rawfile = rawfile, scan = i) class(S) class(S[[1]]) summary(S[[1]]) plot(S[[1]], centroid=TRUE) }) ## ----checkScan, echo=FALSE, message=FALSE, eval=FALSE------------------------- # # AGC-related # S[[1]]$`AGC:` # S[[1]]$`AGC Target:` # # Injection time-related # S[[1]]$`Ion Injection Time (ms)` # S[[1]]$`Max. Ion Time (ms):` # # Resolving power-related # S[[1]]$`FT Resolution:` # # HCD-related # S[[1]]$`HCD Energy:` # S[[1]]$`HCD Energy eV:` ## ----plotSN, eval = TRUE, fig.cap = "Spectrum plot using signal-to-noise option. The vertical grey lines indicate the *in-silico* computed y-ions of the peptide precusor LGGNEQVTR++ as calculated by the [protViz]( https://CRAN.R-project.org/package=protViz) package [@protViz]. The horizontal blue line indicates a signal to noise ratio of 5. Diagnostic information related to AGC is shown in grey."---- plot(S[[1]], centroid=TRUE, SN = TRUE, diagnostic = TRUE) # S/N threshold indicator abline(h = 5, lty = 2, col = "blue") # decorate plot with y-ion series of target peptide LGGNEQVTR++ (yIonSeries <- protViz::fragmentIon("LGGNEQVTR")[[1]]$y[1:8]) names(yIonSeries) <- paste0("y", seq(1, length(yIonSeries))) abline(v = yIonSeries, col='#DDDDDD88', lwd=5) axis(3, yIonSeries, names(yIonSeries)) ## ----------------------------------------------------------------------------- # value casting (x <- S[[1]]$scan) class(x) (y <- rawrr::scanNumber(S[[1]])) class(y) # error handling S[[1]]$`FAIMS Voltage On:` try(rawrr::faimsVoltageOn(S[[1]])) # generating a well-behaved accessor maxIonTime <- rawrr::makeAccessor(key = "Max. Ion Time (ms):", returnType = "double") maxIonTime(S[[1]]) ## ----TIC, fig.cap="Total ion chromatogram (TIC) calculated from all MS1-level scans contained in 20181113_010_autoQC01.raw."---- message(rawfile) rawrr::readChromatogram(rawfile = rawfile, type = "tic") |> plot() ## ----BPC, fig.cap="Base peak chromatogram (BPC) calculated from all MS1-level scans contained in 20181113_010_autoQC01.raw."---- rawrr::readChromatogram(rawfile = rawfile, type = "bpc") |> plot() ## ----plotrawrrchromatogram, fig.cap="Extracted ion chromatograms (XIC) for iRT peptide precursors. Each XIC was calculated using a tolerance of 10 ppm around the target mass and using only MS1-level scans.", error=TRUE---- try({ iRTmz <- c(487.2571, 547.2984, 622.8539, 636.8695, 644.8230, 669.8384, 683.8282, 683.8541, 699.3388, 726.8361, 776.9301) names(iRTmz) <- c("LGGNEQVTR", "YILAGVENSK", "GTFIIDPGGVIR", "GTFIIDPAAVIR", "GAGSSEPVTGLDAK", "TPVISGGPYEYR", "VEATFGVDESNAK", "TPVITGAPYEYR", "DGLDAASYYAPVR", "ADVTPADFSEWSK", "LFLQFGAQGSPFLK") C <- rawrr::readChromatogram(rawfile, mass = iRTmz, tol = 10, type = "xic", filter = "ms") class(C) plot(C, diagnostic = TRUE) }) ## ----fragmentIonTraces, fig.cap="Extracted ion chromatograms (XICs) for LGGNEQVTR++ fragment ions y6, y7, and y8. Each target ion (see legend) was extracted with a tolerance of 10 ppm from all scans matching the provided filter.", fig.retina=1---- rawrr::readChromatogram(rawfile = rawfile, mass = yIonSeries[c("y6", "y7", "y8")], type = 'xic', tol = 10, filter = "FTMS + c NSI Full ms2 487.2567@hcd27.00 [100.0000-1015.0000]") |> plot(diagnostic = TRUE) ## ----readIndex---------------------------------------------------------------- rawrr::readIndex(rawfile = rawfile) |> subset(scanType == "FTMS + c NSI Full ms2 487.2567@hcd27.00 [100.0000-1015.0000]") |> head() ## ----fittedAPEX, fig.cap="The plots show the fitted peak curves in red."------ par(mfrow = c(3, 4), mar=c(4,4,4,1)) rtFittedAPEX <- C |> rawrr:::pickPeak.rawrrChromatogram() |> rawrr:::fitPeak.rawrrChromatogram() |> lapply(function(x){ plot(x$times, x$intensities, type='p', ylim=range(c(x$intensities,x$yp)), main=x$mass); lines(x$xx, x$yp, col='red'); x}) |> sapply(function(x){x$xx[which.max(x$yp)[1]]}) ## a simple alternative to derive rtFittedAPEX could be rt <- sapply(C, function(x) x$times[which.max(x$intensities)[1]]) ## ----iRTscoreFit, error=TRUE-------------------------------------------------- try({ iRTscore <- c(-24.92, 19.79, 70.52, 87.23, 0, 28.71, 12.39, 33.38, 42.26, 54.62, 100) fit <- lm(rtFittedAPEX ~ iRTscore) }) ## ----iRTscoreFitPlot, fig.small=TRUE, echo=FALSE, fig.cap="iRT regression. Plot shows observed peptide RTs as a function of iRT scores and fitted regression line of corresponding linear model obtained by ordinary least squares (OLS) regression."---- # iRTscoreFitPlot plot(rtFittedAPEX ~ iRTscore, ylab = 'Retention time [min]', xlab = "iRT score", pch=16, frame.plot = FALSE) abline(fit, col = 'grey') abline(v = 0, col = "grey", lty = 2) legend("topleft", legend = paste("Regression line: ", "rt =", format(coef(fit)[1], digits = 4), " + ", format(coef(fit)[2], digits = 2), "score", "\nR2: ", format(summary(fit)$r.squared, digits = 4)), bty = "n", cex = 0.75) text(iRTscore, rt, iRTmz, pos=1,cex=0.5) ## ----benchmark, message = FALSE, fig.cap="`rawrr:::.benchmark using rawrr::readSpectrum` spectra per second.", fig.small=TRUE---- set.seed(1) seq(1, 4) |> lapply(FUN=function(x){rawrr::sampleFilePath() |> rawrr:::.benchmark()}) |> Reduce(f = rbind) -> S boxplot(count / runTimeInSec ~ count, data = S, log ='y', sub = paste0("Overall runtime took ", round(sum(S$runTimeInSec), 3), " seconds."), xlab = 'number of random generated scan ids') legend("topleft", c(Sys.info()['nodename'], Sys.info()['sysname']), cex = 1) ## ----faq1--------------------------------------------------------------------- # use sample.raw file f <- rawrr::sampleFilePath() (rawrr::readIndex(f) |> subset(MSOrder == "Ms"))$scan |> sapply(FUN = rawrr::readSpectrum, rawfile=f) |> parallel::mclapply(FUN = function(x){ idx <- x$intensity == max(x$intensity); data.frame(scan=x$scan, StartTime=x$StartTime, mZ=x$mZ[idx], intensity=x$intensity[idx] )}) |> base::Reduce(f=rbind) ## ----json, eval=FALSE--------------------------------------------------------- # #R # .rawrrHeaderJson <- function(inputRawFile, # outputJsonFile = tempfile(fileext = ".json")){ # inputRawFile |> # rawrr::readFileHeader() |> # rjson::toJSON(indent = 1) |> # cat(file = outputJsonFile) # } ## ----sessioninfo, echo=FALSE-------------------------------------------------- sessionInfo()