# Signal drift and batch effect correction for mass spectrometry

Andris Jankevics1* and Ralf Johannes Maria Weber1**

1Phenome Centre Birmingham, University of Birmingham, UK

*a.jankevics@bham.ac.uk
**r.j.weber@bham.ac.uk

#### 2022-04-26

#### Package

pmp 1.8.0

# 1 Introduction

This vignette demonstrates how to apply Quality Control-Robust Spline Correction (QC-RSC) (Kirwan et al., 2013) algorithm for signal drift and batch effect correction within/across a multi-batch direct infusion mass spectrometry (DIMS) and liquid chromatography mass spectrometry (LCMS) datasets.

Please read “Signal drift and batch effect correction and mass spectral quality assessment” vignette to learn how to assess your dataset and details on algorithm itself.

# 2 Installation

You should have R version 4.0.0 or above and Rstudio installed to be able to run this notebook.

Execute following commands from the R terminal.

    install.packages("gridExtra")

    if (!requireNamespace("BiocManager", quietly = TRUE))
        install.packages("BiocManager")
    BiocManager::install("pmp")

Load the required libraries into the R environment

    library(S4Vectors)
    library(SummarizedExperiment)
    library(pmp)
    library(ggplot2)
    library(reshape2)
    library(gridExtra)

# 3 Dataset

In this tutorial we will be using a direct infusion mass spectrometry (DIMS) dataset consisting of 172 samples measured across 8 batches and is included in `pmp` package as `SummarizedExperiemnt` class object `MTBLS79`. More detailed description of the dataset is available from Kirwan et al. (2014), MTBLS79 and R man page.

    help ("MTBLS79")

    data("MTBLS79")

    class <- MTBLS79$Class
    batch <- MTBLS79$Batch
    sample_order <- c(1:ncol(MTBLS79))

    # Input data structure
    MTBLS79
    #> class: SummarizedExperiment
    #> dim: 2488 172
    #> metadata(0):
    #> assays(1): ''
    #> rownames(2488): 70.03364 70.03375 ... 569.36369 572.36537
    #> rowData names(0):
    #> colnames(172): batch01_QC01 batch01_QC02 ... Batch08_C09 Batch08_QC39
    #> colData names(4): Batch Sample_Rep Class Class2

    class[1:10]
    #>  [1] "QC" "QC" "QC" "C"  "S"  "C"  "QC" "S"  "C"  "S"
    batch[1:10]
    #>  [1] "1" "1" "1" "1" "1" "1" "1" "1" "1" "1"
    sample_order[1:10]
    #>  [1]  1  2  3  4  5  6  7  8  9 10

# 4 Missing values

Current implementation of `QCRSC` algorithm does support missing values in the input data object, but we would recommend to filter out features which were net reproducibly measured across quality control (QC) sample. In this example we will use 80% detection threshold.

    data <- filter_peaks_by_fraction(df=MTBLS79, classes=class, method="QC",
        qc_label="QC", min_frac=0.8)

# 5 Applying signal drift and batch effect correction

Function `QCRSC` should be used to apply signal batch correction.

Argument `df` should be `SummarizedExperiment` object or matrix-like R data structure with all `numeric()` values.

Argument `order` should be `numeric()` vector containing sample injection order during analytical measurement and should be the same length as number of sample in the input object.

Argument `batch` should be `numeric()` or `character()` vector containing values of sample batch identifier. If all samples were measured in 1 batch, then all values in the `batch` vector should be identical.

Values for `classes` should be character vector containing sample class labels. Class label for quality control sample has to be `QC`.

    corrected_data <- QCRSC(df=data, order=sample_order, batch=batch,
        classes=class, spar=0, minQC=4)
    #> The number of NA and <= 0 values in peaksData before QC-RSC: 15330

# 6 Visual comparison of the results

Function ‘sbc_plot’ provides visual comparison of the data before and after correction. For example we can check output for features ‘1’, ‘5’, and ‘30’ in peak matrix.

    plots <- sbc_plot (df=MTBLS79, corrected_df=corrected_data, classes=class,
        batch=batch, output=NULL, indexes=c(1, 5, 30))
    plots
    #> [[1]]

    #>
    #> [[2]]

    #>
    #> [[3]]

The scores plots of principal components analysis (PCA) before and after correction can be used to asses effects of data correction.

In this example, probabilistic quotient normalisation (`PQN`) method is used to normalise data, k-nearest neighbours (`KNN`) for missing value imputation and `glog` for data scaling. All functions are availiable as a part of pmp package.

See Di Guida et al. (2016) for a more detailed review on common pre-processing steps and methods.

    manual_color = c("#386cb0", "#ef3b2c", "#7fc97f", "#fdb462", "#984ea3",
        "#a6cee3", "#778899", "#fb9a99", "#ffff33")

    pca_data <- pqn_normalisation(MTBLS79, classes=class, qc_label="QC")
    pca_data <- mv_imputation(pca_data, method="KNN", k=5, rowmax=0.5,
        colmax=0.5, maxp=NULL, check_df=FALSE)
    pca_data <- glog_transformation(pca_data, classes=class, qc_label="QC")

    pca_corrected_data <- pmp::pqn_normalisation(corrected_data, classes=class,
        qc_label="QC")
    pca_corrected_data <- pmp::mv_imputation(pca_corrected_data, method="KNN", k=5,
        rowmax=0.5, colmax=0.5, maxp=NULL, check_df=FALSE)
    #> Warning in knnimp(x, k, maxmiss = rowmax, maxp = maxp): 8 rows with more than 50 % entries missing;
    #>  mean imputation used for these rows
    pca_corrected_data <- pmp::glog_transformation(pca_corrected_data,
        classes=class, qc_label="QC")

    pca_data <- prcomp(t(assay(pca_data)), center=TRUE, scale=FALSE)
    pca_corrected_data <- prcomp(t(assay(pca_corrected_data)),
        center=TRUE, scale=FALSE)

    # Calculate percentage of explained variance of the first two PC's
    exp_var_pca <- round(((pca_data$sdev^2)/sum(pca_data$sdev^2)*100)[1:2],2)
    exp_var_pca_corrected <- round(((pca_corrected_data$sdev^2) /
        sum(pca_corrected_data$sdev^2)*100)[1:2],2)

    plots <- list()

    plotdata <- data.frame(PC1=pca_data$x[, 1], PC2=pca_data$x[, 2],
        batch=as.factor(batch), class=class)

    plots[[1]] <- ggplot(data=plotdata, aes(x=PC1, y=PC2, col=batch)) +
        geom_point(size=2) + theme(panel.background=element_blank()) +
        scale_color_manual(values=manual_color) +
        ggtitle("PCA scores, before correction") +
        xlab(paste0("PC1 (", exp_var_pca[1] ," %)")) +
        ylab(paste0("PC2 (", exp_var_pca[2] ," %)"))

    plots[[2]] <- ggplot(data=plotdata, aes(x=PC1, y=PC2, col=class)) +
        geom_point(size=2) + theme(panel.background=element_blank()) +
        scale_color_manual(values=manual_color) +
        ggtitle("PCA scores, before correction") +
        xlab(paste0("PC1 (", exp_var_pca[1] ," %)")) +
        ylab(paste0("PC2 (", exp_var_pca[2] ," %)"))

    plotdata_corr <- data.frame(PC1=pca_corrected_data$x[, 1],
        PC2=pca_corrected_data$x[, 2], batch=as.factor(batch), class=class)

    plots[[3]] <- ggplot(data=plotdata_corr, aes(x=PC1, y=PC2, col=batch)) +
        geom_point(size=2) +
        theme(panel.background=element_blank()) +
        scale_color_manual(values=manual_color) +
        ggtitle("PCA scores, after correction") +
        xlab(paste0("PC1 (", exp_var_pca_corrected[1] ," %)")) +
        ylab(paste0("PC2 (", exp_var_pca_corrected[2] ," %)"))

    plots[[4]] <- ggplot(data=plotdata_corr, aes(x=PC1, y=PC2, col=class)) +
        geom_point(size=2) +
        theme(panel.background=element_blank()) +
        scale_color_manual(values=manual_color) +
        ggtitle("PCA scores, after correction") +
        xlab(paste0("PC1 (", exp_var_pca_corrected[1] ," %)")) +
        ylab(paste0("PC2 (", exp_var_pca_corrected[2] ," %)"))

    grid.arrange(ncol=2, plots[[1]], plots[[2]], plots[[3]], plots[[4]])

# 7 Session information

    sessionInfo()
    #> R version 4.2.0 RC (2022-04-19 r82224 ucrt)
    #> Platform: x86_64-w64-mingw32/x64 (64-bit)
    #> Running under: Windows Server x64 (build 20348)
    #>
    #> Matrix products: default
    #>
    #> locale:
    #> [1] LC_COLLATE=C
    #> [2] LC_CTYPE=English_United States.utf8
    #> [3] LC_MONETARY=English_United States.utf8
    #> [4] LC_NUMERIC=C
    #> [5] LC_TIME=English_United States.utf8
    #>
    #> attached base packages:
    #> [1] stats4    stats     graphics  grDevices utils     datasets  methods
    #> [8] base
    #>
    #> other attached packages:
    #>  [1] gridExtra_2.3               reshape2_1.4.4
    #>  [3] ggplot2_3.3.5               SummarizedExperiment_1.26.0
    #>  [5] Biobase_2.56.0              GenomicRanges_1.48.0
    #>  [7] GenomeInfoDb_1.32.0         IRanges_2.30.0
    #>  [9] S4Vectors_0.34.0            BiocGenerics_0.42.0
    #> [11] MatrixGenerics_1.8.0        matrixStats_0.62.0
    #> [13] pmp_1.8.0                   BiocStyle_2.24.0
    #>
    #> loaded via a namespace (and not attached):
    #>  [1] sass_0.4.1             splines_4.2.0          jsonlite_1.8.0
    #>  [4] foreach_1.5.2          bslib_0.3.1            assertthat_0.2.1
    #>  [7] BiocManager_1.30.17    highr_0.9              doRNG_1.8.2
    #> [10] GenomeInfoDbData_1.2.8 yaml_2.3.5             impute_1.70.0
    #> [13] pillar_1.7.0           lattice_0.20-45        glue_1.6.2
    #> [16] digest_0.6.29          XVector_0.36.0         randomForest_4.7-1
    #> [19] colorspace_2.0-3       htmltools_0.5.2        Matrix_1.4-1
    #> [22] plyr_1.8.7             pkgconfig_2.0.3        magick_2.7.3
    #> [25] bookdown_0.26          zlibbioc_1.42.0        purrr_0.3.4
    #> [28] scales_1.2.0           itertools_0.1-3        tibble_3.1.6
    #> [31] mgcv_1.8-40            generics_0.1.2         farver_2.1.0
    #> [34] ellipsis_0.3.2         withr_2.5.0            cli_3.3.0
    #> [37] magrittr_2.0.3         crayon_1.5.1           evaluate_0.15
    #> [40] fansi_1.0.3            nlme_3.1-157           tools_4.2.0
    #> [43] lifecycle_1.0.1        stringr_1.4.0          munsell_0.5.0
    #> [46] rngtools_1.5.2         DelayedArray_0.22.0    missForest_1.5
    #> [49] pcaMethods_1.88.0      compiler_4.2.0         jquerylib_0.1.4
    #> [52] rlang_1.0.2            grid_4.2.0             RCurl_1.98-1.6
    #> [55] iterators_1.0.14       bitops_1.0-7           labeling_0.4.2
    #> [58] rmarkdown_2.14         gtable_0.3.0           codetools_0.2-18
    #> [61] DBI_1.1.2              R6_2.5.1               knitr_1.38
    #> [64] dplyr_1.0.8            fastmap_1.1.0          utf8_1.2.2
    #> [67] stringi_1.7.6          parallel_4.2.0         Rcpp_1.0.8.3
    #> [70] vctrs_0.4.1            tidyselect_1.1.2       xfun_0.30

# References

Di Guida, R., Engel, J., Allwood, J.W., et al. (2016) Non-targeted uhplc-ms metabolomic data processing methods: A comparative investigation of normalisation, missing value imputation, transformation and scaling. Metabolomics, 12 (5): 93. doi:10.1007/s11306-016-1030-9.

Kirwan, J.A., Weber, R.J., Broadhurst, D.I., et al. (2014) Direct infusion mass spectrometry metabolomics dataset: A benchmark for data processing and quality control. Scientific data, 1: 140012. doi:10.1038/sdata.2014.12.

Kirwan, J., Broadhurst, D., Davidson, R., et al. (2013) Characterising and correcting batch variation in an automated direct infusion mass spectrometry (dims) metabolomics workflow. Analytical and bioanalytical chemistry, 405 (15): 5147–5157. doi:10.1007/s00216-013-6856-7.