Example data

The example file 20181113_010_autoQC01.raw used throughout this manuscript contains Fourier-transformed Orbitrap spectra (FTMS) recorded on a Thermo Fisher Scientific Q Exactive HF in positive mode (+). The mass spectrometer was operated in line with a nano UPLC and a nano electrospray source (NSI). MS2 spectra were generated by HCD fragmentation at normalized collision energy (NCE) of 27. All spectra were written to file after applying centroiding (c) and lock mass correction. The analyzed sample consisted of the iRT peptide mix (Biognosys) in a tryptic BSA digest (NEB) and was separated applying a 20 min linear gradient on C18 reversed-phase material at a constant flow rate of 300 nl/min. The file is part of the MassIVE dataset MSV000086542 (Kockmann 2020) and can be obtained through the FTP download link (MD5: a1f5df9627cf9e0d51ec1906776957ab). Individual scans have been assigned a universal spectrum identifier (USI)(Deutsch et al. 2020) by MassIVE.

If the environment variable MONO_PATH does not include a directory containing the RawFileReader .NET assemblies are installed in a directory derived by the rawrr::rawrrAssemblyPath() function.

rawrr::installRawrrExe()

## Overwrite sourceUrl to https://fgcz-ms.uzh.ch/~cpanse/rawrr/dotnet//linux-x64/rawrr

## MD5 d842e4ba5fe901cc4dd62e439e7c6b9b /github/home/.cache/R/rawrr/rawrrassembly/linux-x64/rawrr

## [1] 0

Additional raw data for demonstration and extended testing is available through the Bioconductor data package tartare (Panse and Kockmann 2019).

# fetch via ExperimentHub
library(ExperimentHub)
eh <- ExperimentHub::ExperimentHub()
EH4547 <- normalizePath(eh[["EH4547"]])

(rawfile <- paste0(EH4547, ".raw"))

## [1] "/github/home/.cache/R/ExperimentHub/12f17984537c_4590.raw"

if (!file.exists(rawfile)){
  file.copy(EH4547, rawfile)
}

## [1] TRUE

Use Case I - Analyzing Orbitrap Spectra

H <- rawrr::readFileHeader(rawfile = rawfile)

The Orbitrap detector has been a tremendous success story in MS, since it offers HRAM data on a time scale that is compatible with chromatographic analysis (LC-MS)(Makarov 2000) and is therefore heavily used in bottom-up proteomics. However, analyzing Orbitrap data in R has so far only been possible after raw data conversion to exchange formats like mz(X)ML. Unfortunately, conversion is accompanied by a loss of Orbitrap-specific information. This use case shows how easy it is to work directly with raw-formated Orbitrap data after installing our R package rawrr which applies vendor APIs for data access. We use a complete LC-MS run recorded on a Q Exactive HF Orbitrap by parallel reaction monitoring (PRM)(Gallien et al. 2012) for demonstration purposes (File name: 12f17984537c_4590.raw). The 35 min run resulted in 21881 scans that were written to the file. Already typesetting the above lines uses rawrr functionality, since instrument model, file name, time range of data acquisition, and number of scans is extracted from the binary file header (Note: This manuscript was written in R markdown and combines R code with narration). The respective function is called readFileHeader() and returns a simple R object of type list (see Table 1).

Individual scans or scan collections (sets) can be read by the function readSpectrum() which returns a rawrrSpectrum object or rawrrSpectrumSet. Our package also provides generics for printing and plotting these objects. The following code chunk depicts how a set of scans is read from the raw file (scan numbers were selected based on a database search). The corresponding Figure 1 shows the resulting plot for scan 9594 (USI: mzspec:MSV000086542:20181113_010_autoQC01:scan:9594:LGGNEQVTR/2) assigned to the doubly-charged iRT peptide LGGNEQVTR by MS-GF+ (Score: 144, SpecProb: 1.9e-12, DB E-Value: 4.4e-4, see MassIVE RMSV000000336.1 for details of the search):

i <- c(9594, 11113, 11884, 12788, 12677, 13204, 13868, 14551, 16136, 17193, 17612)
S <- rawrr::readSpectrum(rawfile = rawfile, scan = i)
class(S)

## [1] "rawrrSpectrumSet"

class(S[[1]])

## [1] "rawrrSpectrum"

summary(S[[1]])

## Total Ion Current:    62374688
## Scan Low Mass:    100
## Scan High Mass:   1015
## Scan Start Time (min):    15.42042
## Scan Number:  9594
## Base Peak Intensity:  12894520
## Base Peak Mass:   860.4223
## Scan Mode:    FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]

plot(S[[1]], centroid=TRUE)

Plot of scan number 9594 (USI: mzspec:MSV000086542:20181113_010_autoQC01:scan:9594) showing a centroided tandem mass spectrum of the iRT peptide precursor LGGNEQVTR++ in positive mode. The scan was acquired on an Orbitrap detector including lock mass correction and using a transient of 64 ms (equal to a resolving power of 30’000 at 200 m/z) and an AGC target of 1e5 elementary charges.Peak attributes like m/z, charge (z), and resolution (R) are shown above the peaks.

The plot shows typical Orbitrap peak attributes like resolution (R) and charge (z) above the most intense peaks when centroided data is available and selected. Centroided data also makes it possible to graph spectra using signal-to-noise as response value. This is potentially interesting since Orbitrap detectors follow $S/N \propto charges \cdot \sqrt R$ with S/N being the signal to noise ratio, charges meaning elementary charges in the Orbitrap, and R being the resolving power of the scan (Kelstrup et al. 2014). Signal-to-noise makes judging the signal quantity more intuitive than using arbitrary signal intensity units. Figure 2 shows that all y-ion signals are several ten or even hundred folds above the noise estimate (Note: The spectrum shown in Figure 2 is decorated with y-ions calculated for LGGNEQVTR++ using the protViz package.). By adding some diagnostic information related to the Orbitrap’s advanced gain control (AGC) system, it becomes apparent that the C-trap managed to collect the defined 100,000 charges within 2.8 ms, corresponding to only ~5% of the maximum injection time of 55 ms.

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])

## [1] 175.1190 276.1666 375.2350 503.2936 632.3362 746.3791 803.4006 860.4221

names(yIonSeries) <- paste0("y", seq(1, length(yIonSeries)))
abline(v = yIonSeries, col='#DDDDDD88', lwd=5)
axis(3, yIonSeries, names(yIonSeries))

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 package (Panse and Grossmann 2020). The horizontal blue line indicates a signal to noise ratio of 5. Diagnostic information related to AGC is shown in grey.

In total, the API provides 83 data items for this particular scan covering ion optics settings, mass calibration, and other diagnostic data (use print() to gain an overview). It needs to be mentioned that this list is highly dynamic and changes with instrument model and installed instrument control software version. Our reader method is designed in a flexible way so that also newly introduced items will be passed onto the R environment in the future. In general, access to rawrrSpectrum items are provided by the subsetting operators $ and [[. But for programmatic access we highly recommend using accessor functions like scanNumber() instead of using $scan or [["scan"]], since the accessors will also take care of casting values stored as-returned by the API (mostly base type character). Another argument for favoring the accessor function mechanism is the more restrictive error handling in case the requested information is not available and more descriptive function names without white spaces. Users can create missing accessor functions by applying a function factory named makeAccessor() or request them from the package maintainers. The following code chunk demonstrates the above mentioned aspects:

# value casting
(x <- S[[1]]$scan)

## [1] 9594

class(x)

## [1] "numeric"

(y <- rawrr::scanNumber(S[[1]]))

## [1] 9594

class(y)

## [1] "integer"

# error handling
S[[1]]$`FAIMS Voltage On:`

## NULL

try(rawrr::faimsVoltageOn(S[[1]]))

## Error in rawrr::faimsVoltageOn(S[[1]]) : 
##   FAIMS Voltage On: is not available!

# generating a well-behaved accessor
maxIonTime <- rawrr::makeAccessor(key = "Max. Ion Time (ms):", returnType = "double")
maxIonTime(S[[1]])

## [1] 55

More sophisticated analysis workflows applying rawrr functionalities have also been demonstrated recently. For example, rawrr was used to annotate and compare marker ions found in HCD MS2 spectra for ADP-ribosylated peptides at different collision energies (Gehrig et al. 2020), as well as for the annotation of small molecule spectra after UVPD dissociation (Panse et al. 2020). Such information can be conveniently extracted, since the rawrrSpectrum object provides easy access to normalized and absolute HCD energies.

Use Case II - iRT Regression for System Suitability Monitoring

By applying linear regression, one can convert observed peptide retention times (RTs) into dimensionless scores termed iRT values and vice versa (Escher et al. 2012). This can be used for retention time calibration/prediction. In addition, fitted iRT regression models provide highly valuable information about LC-MS run performance. This example shows how easy it is to perform iRT regression in R by just using the raw measurement data, our package rawrr, and well known base R functions supporting linear modeling. To get a first impression of the data we calculate a total ion chromatogram (TIC) using the readChromatogram() function. Plotting the TIC shows chromatographic peaks between 15 and 28 min that could be of peptidic origin (see Figure 3). Of note, there is also a type = "bpc" option if you prefer a base peak chromatogram (BPC):

message(rawfile)

## /github/home/.cache/R/ExperimentHub/12f17984537c_4590.raw

rawrr::readChromatogram(rawfile = rawfile, type = "tic") |>
  plot()

Total ion chromatogram (TIC) calculated from all MS1-level scans contained in 20181113_010_autoQC01.raw.

rawrr::readChromatogram(rawfile = rawfile, type = "bpc") |>
  plot()

Base peak chromatogram (BPC) calculated from all MS1-level scans contained in 20181113_010_autoQC01.raw.

The initial step of iRT regression is to estimate the empirical RTs of a peptide set with known iRT scores. In the simplest case, this is achieved by computing an extracted ion chromatogram (XIC) for iRT peptide precursors, provided they were spiked into the sample matrix prior to data acquisition. The code chunk below demonstrates how the function readChromatogram() is called on the R command line to return a rawrrChromatogramSet object of the type xic. This object is plotted for visual inspection (see Figure 4 for resulting plot).

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)

## [1] "rawrrChromatogramSet"

plot(C, diagnostic = TRUE)

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.

Be reminded that the intensity traces are not computed within R, for instance, by reading all scans of a raw file and subsequently iterating over a scan subset. Instead, traces are retrieved by a C# method that calls the vendor API. The API takes care of the scan filtering process (checks filter validity and applies the filter). On the R code level there is no need to know a priori which scans match the filter rule or implement vectorized operations (we generate multiple XICs simultaneously here). Only the API-returned output needs to be parsed into rawrrChromatogram objects. By changing the scan filter, one can easily switch between generating precursor traces and fragment ion traces. The following code chunk shows how to create fragment ion chromatograms (y6 to y8) generated from scans that target LGGNEQVTR++ (see Figure 5 for plot output):

rawrr::readChromatogram(rawfile = rawfile,
      mass = yIonSeries[c("y6", "y7", "y8")],
      type = 'xic',
      tol = 10,
      filter = "FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]") |>
  plot(diagnostic = TRUE)

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.

By using the readIndex() function a data.frame that indexes all scans found in a raw file is returned. When subsetting it for our scan type of interest it becomes clear that the example data was recorded using parallel reaction monitoring (PRM) since the parent ion 487.2567 was isolated in regularly-spaced intervals for fragment ion recording:

rawrr::readIndex(rawfile = rawfile) |>
  subset(scanType == "FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]") |>
  head()

##     scan                                                     scanType
## 2      2 FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]
## 24    24 FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]
## 46    46 FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]
## 68    68 FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]
## 90    90 FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]
## 112  112 FTMS + c NSI Full ms2 [email protected] [100.0000-1015.0000]
##       StartTime precursorMass MSOrder charge masterScan dependencyType
## 2   0.006862313      487.2567     Ms2      2         NA             NA
## 24  0.042385827      487.2567     Ms2      2         NA             NA
## 46  0.077909329      487.2567     Ms2      2         NA             NA
## 68  0.113423450      487.2567     Ms2      2         NA             NA
## 90  0.148931680      487.2567     Ms2      2         NA             NA
## 112 0.184455230      487.2567     Ms2      2         NA             NA
##     monoisotopicMz
## 2                0
## 24               0
## 46               0
## 68               0
## 90               0
## 112              0

The delta between consecutive scans is always 22 scans (one PRM cycle) and the depedencyType of the MS2 scans is NA (since it was triggered by an inclusion list).

For regression, we now extract the RTs at the maximum of the fitted intensity traces stored in the rawrrChromatogram object and fit a linear model of the form rt = a + b ⋅ score. In theory, we could do this at the precursor or fragment ion level. For simplicity, we show only the first option.

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]])

The plots show the fitted peak curves in red.

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)

The fitted model can then be inspected using standard procedures. Figure 6, shows a visual inspection by plotting observed RTs as a function of iRT score together with the fitted model regression line. The corresponding R-squared indicates that the RTs behave highly linear. This is expected since the iRT peptides were separated on a 20 min linear gradient from 5% buffer B to 35% buffer B using C18 reversed-phase material (the change rate is therefore constant at 1.5% per minute). The magnitude of the slope parameter (b) is a direct equivalent of this gradient change rate. The intercept (a) is equal to the predicted RT of iRT peptide GAGSSEPVTGLDAK since it was defined to have a zero score on the iRT scale.

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.

Extension

An extended and dynamic version of the above use cases can be found at (https://fgcz-ms.uzh.ch/~cpanse/rawrr/test/functional_test.html). The web page displays spectra and iRT regression models obtained over a set of raw files recorded approximately every 12 hours on different Orbitrap mass spectrometers at the FGCZ (some systems have gone out of service in the meantime). The original purpose of these injections is automated longitudinal system suitability monitoring and quality control. We re-use the resulting raw files to showcase rawrr’s functionality across different Orbitrap instrument models/generations. In order to find the highest-scoring MS2 scan for iRT peptides we now use a simple scoring function implemented in R (it counts the number of matching y-ions), instead of running an external search engine. The web page automatically updates every 30 minutes using the most recent two files per system as input data. Be aware that the code is executed in a full parallel fashion (each core processes one raw file) on a Linux server with network-attached storage.

The R script that renders the html page is also available as supplementary information. The correspoding R markdown file is part of the rawrr package and can be processed locally after downloading a snapshot of the above described input data from MSV000086542 (Kockmann 2020). In summary, this shows how scalable analysis pipelines can be constructed starting from basic building blocks. It demonstrates that rawrr’s data access mechanism works for all types of Orbitrap instrument models.

Benchmark

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)

rawrr:::.benchmark using rawrr::readSpectrum spectra per second.

How to pronounce the package name the right way?

see https://youtu.be/jBc2MniDBYw

What would be the simplest way of extracting the base peak intensity and m/z value for a large batch of files at the MS1 level?

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

##    scan   StartTime       mZ intensity
## 1     1 0.001619751 445.1181   5979308
## 2    23 0.031642766 445.1180   5902284
## 3    45 0.061663615 445.1182   5630335
## 4    67 0.091651065 445.1182   5857274
## 5    89 0.121640750 391.2826   5366674
## 6   111 0.151644970 445.1181   5926550
## 7   133 0.181667770 445.1182   5840996
## 8   155 0.211526280 445.1182   5989230
## 9   177 0.241284530 391.2826  16794288
## 10  199 0.271307600 445.1182   5756987
## 11  221 0.301222200 391.2826   6351876
## 12  243 0.331145000 445.1184   5810201
## 13  265 0.361147270 445.1184   5687840
## 14  287 0.391168050 445.1183   6016346
## 15  309 0.421057700 391.2827   7898428
## 16  331 0.450970250 445.1183   6577130
## 17  353 0.481045180 391.2826   6162342
## 18  375 0.510989030 391.2827   6502305
## 19  397 0.540955750 445.1182   5700832
## 20  419 0.570893130 391.2826   7186618
## 21  441 0.600724580 391.2825   6955724
## 22  463 0.630620300 445.1183   6123972
## 23  485 0.660428770 391.2827  16238321
## 24  507 0.690318400 391.2827   6862950
## 25  529 0.720320350 391.2826   5815317
## 26  551 0.750197620 445.1183   5746604
## 27  573 0.780105920 391.2824   7327204

mZ and intensity vectors have different lengths. What shall I do?

On Windows, the decimal symbol has to be configured as a ‘.’! See also #fgcz/rawDiag/issues/33.

Howto write the rawrr::readFileHeader() output into json file?

#R
.rawrrHeaderJson <- function(inputRawFile,
                             outputJsonFile = tempfile(fileext = ".json")){
  inputRawFile |>
  rawrr::readFileHeader() |>
    rjson::toJSON(indent = 1) |>
    cat(file = outputJsonFile)
}