In this document we will apply an automated screening for possible candidate signals that show a response to treatment with a compound. For details how to do the experimental part check out our Nature Protocols publication Unger 2021 “Label-free Cell Assays of Compound Uptake and Drug Action using MALDI-TOF Mass Spectrometry”. Briefly, cells were treated with the following concentration of a compound: 0, 0.04, 0.12, 0.37, 1.11, 3.33, 10, 30 uM. After incubation the cells were washed, transferred to a MALDI target plate and matrix was applied. For each concentration 4 spots were applied which means that we also have 4 “measurement replicates” for each concentration.
The most important function in this package is fitCurve(). It does not only fit the dose response curve and compute score values but it also does most of the necessary steps needed to prepare the data for fitting like normalization and alignment. In this tutorial we want to first look at each of this steps individually and then apply the fitCurve() function the helps us combine all of them.
For the sake to keep this package small, the spectra used as an example here were trimmed to 400-900 Da mass-range and stored in a compact form as example data.
library(MALDIcellassay)
library(MALDIquant)
#>
#> This is MALDIquant version 1.22.1
#> Quantitative Analysis of Mass Spectrometry Data
#> See '?MALDIquant' for more information about this package.First we will check the spectra for general quality.
The baseline of the spectra already looks ok but for the sake of this tutorial we will apply a baseline correction. But first we want so save the names of the spectra (which are the concentrations used to treat the cells). Also to get a better overview of the data we compute the mean spectra of each concentration (there are 4 measurement replicates for each concentration) and plot them around m/z 760 as we want to use this signal to normalize and also as lock mass.
conc <- as.numeric(names(Blank2022spec))
spec_prc <- MALDIquant::removeBaseline(Blank2022spec)
avg <- MALDIquant::averageMassSpectra(spec_prc, labels = conc)
MALDIquant::plot(avg[[1]], main = "Overview of mean spectra", xlim = c(755, 765))
for(i in 2:length(avg)) {
MALDIquant::lines(avg[[i]], col = i)
}
legend("topright", legend = paste0(names(avg), "uM"), col = 1:8, lty=1)
We see that the alignment of the spectra is quite good but the intensity
varies for the peak around m/z 760. To get best results its necessary to
normalize the spectra to each other. One way is to add a standard to
each measurement spot (internal standard). This is not always possible
and another option is to use endogenous peaks with high abundance like
the signal around 760 m/z. As an example we will now find out the m/z
deviation to this signal for each spectrum. First we need to detect the
peaks, then we convert the peaks to a data.frame and finally we use
getMzShift() to find the mass shift in each spectrum.
peaks <- MALDIquant::detectPeaks(Blank2022spec, method = "SuperSmoother", SNR = 5)
names(peaks) <- names(Blank2022spec)
mz_shift <- getMzShift(peaks = peaks, targetMz = 760.585, tol = 500)
#> found mz 760.585 in 32 / 32 spectra
#> 15:58 mzshift was -0.0796625 in mean and 0.0828 abs. max.
summary(mz_shift$mzshift)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> -0.08280 -0.08280 -0.08280 -0.07966 -0.08280 -0.05770As we already suspected the shift is quite small. Now lets use it to align the spectra
Now that we have aligned spectra the next step is to normalize them. Also here, a internal standard would be great but as we dont have one we again use the m/z 760 endogenous signal. First we extract the intensities for each spectrum of this signal and then we use it as a normalization factor.
peaks_align <- MALDIquant::detectPeaks(spec_align, method = "SuperSmoother", SNR = 3)
norm <- getNormFactors(peaksdf = peaks2df(peaks_align), targetMz = 760.585, tol = 100)
summary(norm$norm_factor)
#> Min. 1st Qu. Median Mean 3rd Qu. Max.
#> 7949 12842 14492 14099 16145 19570
spec_rdy <- normalizeByFactor(spec_align, norm$norm_factor)As a final check we now plot again the aligned and normalizes mean spectra for each concentration.
avg_rdy <- MALDIquant::averageMassSpectra(spec_rdy, labels = conc)
MALDIquant::plot(avg_rdy[[1]], main = "Overview of mean spectra", xlim = c(755, 765))
for(i in 2:length(avg_rdy)) {
MALDIquant::lines(avg_rdy[[i]], col = i)
}
legend("topright", legend = paste0(names(avg), "uM"), col = 1:8, lty=1)
As expected we find our spectra to be perfectly aligned and normalized to 760.585. Now we can again detect peaks, bin them and compose a intensity matrix of all peaks.
peaks_rdy <- MALDIquant::detectPeaks(avg_rdy, method = "SuperSmoother", SNR = 3)
peaks_rdy <- MALDIquant::binPeaks(peaks_rdy)
intmat <- MALDIquant::intensityMatrix(peaks_rdy, avg_rdy)
dim(intmat)
#> [1] 8 211We get a intensity matrix with 211 peaks in total. The next step is to screen for high variant peaks.
We end up with 20 candidate peaks for which we want to fit a IC50 curve. Lets do the fitting for one exemplary signal.
concLog <- log10(unique(conc))
if(any(concLog == -Inf)) {
concLog[which(concLog == -Inf)] <- (min(concLog[which(!concLog == -Inf)])-1)
}
resp <- nplr::convertToProp(y = intmat[,10])
model <- nplr::nplr(x = concLog, y = resp, useLog = FALSE, npars = 4)
title <- paste0("m/z =", round(as.numeric(colnames(intmat)[12]), 2))
plot(model, main = title)
And with that we are finished! As said during the introduction we did
each step manually. To do what we did in an automated way for each
signal of interest use the following function:
sessionInfo()
#> R version 4.3.2 (2023-10-31 ucrt)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows 10 x64 (build 19045)
#>
#> Matrix products: default
#>
#>
#> locale:
#> [1] LC_COLLATE=C LC_CTYPE=German_Germany.utf8
#> [3] LC_MONETARY=German_Germany.utf8 LC_NUMERIC=C
#> [5] LC_TIME=German_Germany.utf8
#>
#> time zone: Europe/Berlin
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] MALDIquant_1.22.1 MALDIcellassay_0.4.47
#>
#> loaded via a namespace (and not attached):
#> [1] gtable_0.3.4 jsonlite_1.8.7 highr_0.10
#> [4] dplyr_1.1.4 compiler_4.3.2 tidyselect_1.2.0
#> [7] parallel_4.3.2 tidyr_1.3.0 jquerylib_0.1.4
#> [10] scales_1.3.0 yaml_2.3.7 fastmap_1.1.1
#> [13] ggplot2_3.5.0 R6_2.5.1 generics_0.1.3
#> [16] knitr_1.45 forcats_1.0.0 XML_3.99-0.15
#> [19] tibble_3.2.1 nplr_0.1-7 munsell_0.5.0
#> [22] readBrukerFlexData_1.9.1 bslib_0.6.0 pillar_1.9.0
#> [25] rlang_1.1.2 utf8_1.2.4 cachem_1.0.8
#> [28] xfun_0.41 sass_0.4.7 cli_3.6.1
#> [31] withr_2.5.2 magrittr_2.0.3 digest_0.6.33
#> [34] grid_4.3.2 rstudioapi_0.15.0 base64enc_0.1-3
#> [37] lifecycle_1.0.4 vctrs_0.6.4 readMzXmlData_2.8.3
#> [40] evaluate_0.23 glue_1.6.2 MALDIquantForeign_0.14
#> [43] svMisc_1.2.3 colorspace_2.1-0 fansi_1.0.5
#> [46] purrr_1.0.2 rmarkdown_2.25 tools_4.3.2
#> [49] pkgconfig_2.0.3 htmltools_0.5.7