---
title: "Longitudinal System Suitability Monitoring and Quality Control with MSstatsQC"
author: 
  - "Eralp Dogu <eralp.dogu@gmail.com>"
  - "Sara Taheri <srtaheri66@gmail.com>"
  - "Olga Vitek <o.vitek@neu.edu>"
date: "`r Sys.Date()`"
output: 
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 2
    number_sections: true
vignette: >
  %\VignetteIndexEntry{Longitudinal System Suitability Monitoring and Quality Control with MSstatsQC}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,   
  message = FALSE,  
  fig.width = 7,     
  fig.height = 5     
)
```
# Introduction to MSstatsQC

Liquid chromatography coupled with mass spectrometry (LC-MS) is the gold standard for peptide quantification in complex matrices. However, achieving reproducible results across different laboratories and runs requires rigorous monitoring. This involves two critical components: **System Suitability Tests (SST)** to verify instrumentation performance, and **Quality Control (QC)** to ensure in-process quality assurance.

While the longitudinal nature of SST and QC data is well recognized, traditional statistical methods often fall short in detecting subtle drifts. **MSstatsQC** addresses this gap by translating modern **Statistical Process Control (SPC)** methods—such as simultaneous and time-weighted control charts and change point analysis—into the context of LC-MS experiments.

## Key Features & Innovations

* **MSstatsQC-ML (Machine Learning):** The latest release introduces a cutting-edge machine learning approach optimized for decision-making in complex standard mixtures. By combining supervised tree-based classifiers with experimental design strategies, `MSstatsQC-ML` can simulate unobserved, suboptimal MS runs. This allows for the training of robust models even with limited failure data, enabling users to interpret the root causes of performance loss and design preventive actions.

* **Versatility & Compatibility:** The package supports data processing for missing values (v2.0) and includes a converter function for the standardized **mzQC** format.

* **Proven Application:** The methods are demonstrated using diverse datasets, including SRM-based system suitability data (CPTAC Study 9.1, Site 54), DDA and SRM-based QC datasets (QCloud system), and DIA quantification for iRT peptides (QuiC system).

This vignette summarizes the functionalities of the `MSstatsQC` package, which is available as a stand-alone R package or for integration into automated pipelines.

## Installation

**Note on System Requirements:**
The machine learning module (`MSstatsQC-ML`) relies on the `h2o` package, which requires **Java (JDK 8 or higher)** to be installed on your system. Please ensure Java is installed before proceeding.

To install the stable version from Bioconductor:

```{r install-bioc, eval = FALSE}
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("MSstatsQC")
```

# Input

In order to analyze QC/SST data in `MSstatsQC`, input data must be a .csv file in a "long" format with related columns. This is a common data format that can be generated from spectral processing tools such as Skyline and Panorama AutoQC..

The recommended format includes `Acquired Time`, `Peptide name`, `Annotations` and data for any QC metrics such as `Retention Time`, `Total Peak Area` and `Mass Accuracy` etc. Each input file should include `Acquired Time`, `Peptide name` and `Annotations`. After the `Annotations` column user can parse any metric of interest with a proper column name.

(a) `AcquiredTime`: This column shows the acquired time of the QC/SST sample in the format of MM/DD/YYYY HH:MM:SS AM/PM. 

(b) `Precursor`: This column shows information about Precursor id. Statistical analysis will be conducted for each unique label in this column. 

(c) `Annotations`: Annotations are free-text information given by the analyst about each run. They can be information related to any special cause or any observations related to a particular run. Annotations are carried in the plots provided by `MSstatsQC` interactively.  

(d)-(f) `RetentionTime`, `TotalPeakArea`, `FWHM`, `MassAccuracy`, and `PeakAssymetry`, and other metrics: These columns define a *feature* of a peak for a specific peptide. 

The example dataset is shown below. Each row corresponds to a single time point. Additionally, other inputs such as predefined limits or guide sets are discussed in further steps.

### Example

```{r, eval=TRUE}
# A typical multi peptide and multi metric system suitability dataset
# This dataset was generated during CPTAC Study 9.1 at Site 54
library(MSstatsQC)
data <- MSstatsQC::S9Site54
```

# `MSnbaseToMSstatsQC` functions

`MSnbaseToMSstatsQC` function converts `MSnbase` output to `MSstatQC` format using `QCmetrics` objects. 

### Arguments

* `msfile`: Mass spectrometry raw file that contains quality control measurements. 

### Example

```{r, eval=FALSE}
MSnbaseToMSstatsQC(msfile)
```

# Data processing

Data is processed with `DataProcess()` function to ensure data sanity and efficiently use core and summary `MSstatsQC` funtions. `MSstatsQC` uses a data validation method where slight variations in column names are compansated and converted to the standard `MSstatsQC` format. For example, our data validation function converts column names like `Best.RT`, `best retention time`, `retention time`, `rt` and `best ret` into `BestRetentionTime`. This conversion also deals with case-sensitive typing. 

### Arguments

* `data` : comma-separated (.csv), metric file. It should contain a "Precursor" column and the metrics columns. It should also include "Annotations" for each observation.

### Example

```{r, eval=TRUE}
data <- DataProcess(data)
```

# `MSstatsQC` core functions: control charts

The fuction `XmRChart()` is used to generate individual (X) and moving range (mR), and the function `CUSUMChart()` is used to construct cumulative sum for mean (CUSUMm) and cumulative sum for variability (CUSUMv) control charts for each metric. As a follow up change point estimation procedure `ChangePointEstimator` can be used. 

Metrics (e.g. retention time and peak area) and peptides are chosen within all core functions with 'metric' and 'peptide' arguments. `MSstatsQC` can handle any metrics of interest. User needs to create data columns just after `Annotations` to import metrics into `MSstatsQC` successfully. 

Predefined limits are commonly used in system sutiability monitoring and quality control studies. If the mean and variability of a metric is well known, they can be defined using 'selectMean' and 'selectSD' arguments in core plot functions (e.g. XmRplots function). For example, if mean of retention time is 28.5 minutes, standard deviation is 1 minutes and X chart is used for peptide LVNELTEFAK, we use XmRplot function as follows.

The true values of mean and variability of a metric is typically unknown, and their estimates are obtained from a guide set of high quality runs. Generally, a data gathering and parameter estimation step is required. Within that phase, control limits are obtained to test the hypothesis of statistical control. These thresholds are selected to ensure a specified type I error probability (e.g. 0.0027). Constructing control charts and real time evaluation are considered after achieving this phase. Guide sets are defined with 'L' and 'U' arguments. For example, if retention time of a peptide is monitored and first 20 observations of the dataset are used as a guide set, a plot is constructed as follows. 

# `MSstatsQC` core functions: `XmRChart()`

### Arguments

* `data`: comma-separated (.csv), metric file. It should contain a "Precursor" column and the metrics columns. It should also include "Annotations" for each observation.
* `peptide`: the name of precursor of interest.
* `L`: lower bound of the guide set.
* `U`: upper bound of the guide set.
* `metric`: the name of metric of interest.
* `normalization`: TRUE if data is standardized.
* `ytitle`: the y-axis title of the plot. The x-axis title is by default "Time : name of peptide"
* `type`: the type of the control chart. Two values can be assigned, "mean" or "dispersion". Default is "mean"
* `selectMean`: the mean of a metric. It is used when mean is known. It is NULL when mean is not known.  The default is NULL.
* `selectSD`: the standard deviation of a metric. It is used when standard deviation is known. It is NULL when mean is not known. The default is NULL.

### Example

```{r, eval=TRUE, fig.width=8, fig.height=5}
# An X chart when a guide set (1-20 runs) is used to monitor the mean of retention time
XmRChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "BestRetentionTime", normalization = FALSE, ytitle = "X Chart : retention time", type = "mean", selectMean = NULL, selectSD = NULL)
# An X chart when a guide set (1-20 runs) is used to monitor the mean of total peak area
XmRChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "TotalArea", normalization = FALSE, ytitle = "X Chart : peak area", type = "mean", selectMean = NULL, selectSD = NULL)
```

```{r, eval=TRUE, fig.width=8, fig.height=5}
# An mR chart when a guide set (1-20 runs) is used to monitor the variability of total peak area
XmRChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "TotalArea", normalization = TRUE, ytitle = "mR Chart : peak area", type = "variability", selectMean = NULL, selectSD = NULL)
# An mR chart when a guide set (1-20 runs) is used to monitor the variability of retention time
XmRChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "BestRetentionTime", normalization = TRUE, ytitle = "mR Chart : retention time", type = "variability", selectMean = NULL, selectSD = NULL)
# Mean and standard deviation of LVNELTEFAK is known
XmRChart(data, "LVNELTEFAK", metric = "BestRetentionTime", selectMean = 28.5, selectSD = 0.5)
```

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Mean and standard deviation of LVNELTEFAK is known
XmRChart(data, "LVNELTEFAK", metric = "BestRetentionTime", selectMean = 28.5, selectSD = 0.5)
```
# `MSstatsQC` core functions: `CUSUMChart()`

### Arguments

* `data`: comma-separated (.csv), metric file. It should contain a "Precursor" column and the metrics columns. It should also include "Annotations" for each observation.
* `peptide`: the name of precursor of interest.
* `L`: lower bound of the guide set.
* `U`: upper bound of the guide set.
* `metric`: the name of metric of interest.
* `normalization`: TRUE if data is standardized.
* `ytitle`: the y-axis title of the plot. The x-axis title is by default "Time : name of peptide"
* `type`: the type of the control chart. Two values can be assigned, "mean" or "dispersion". Default is "mean"
* `referenceValue`: the value that is used to tune the control chart for a proper shift size. Recommended setting is 0.5 for standardized data.
* `decisionInterval`: the threshold to detect an out-of-control observation. Recommended setting is 5 for standardized data. 
* `selectMean`: the mean of a metric. It is used when mean is known. It is NULL when mean is not known.  The default is NULL.
* `selectSD`: the standard deviation of a metric. It is used when standard deviation is known. It is NULL when mean is not known. The default is NULL.

### Example

```{r, eval=TRUE, echo =FALSE, fig.width=8, fig.height=5}
# A CUSUMm chart when a guide set (1-20 runs) is used to monitor the mean of retention time
CUSUMChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "BestRetentionTime", normalization = TRUE, ytitle = "CUSUMm Chart : retention time", type = "mean", referenceValue = 0.5, decisionInterval = 5, selectMean = NULL, selectSD = NULL)
# A CUSUMm chart when a guide set (1-20 runs) is used to monitor the mean of total peak area
CUSUMChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "TotalArea", normalization = TRUE, ytitle = "CUSUMm Chart : peak area", type = "mean", referenceValue = 0.5, decisionInterval = 5, selectMean = NULL, selectSD = NULL)
```

```{r, eval=TRUE, echo =FALSE, fig.width=8, fig.height=5}
# A CUSUMv chart when a guide set (1-20 runs) is used to monitor the variability of retention time
CUSUMChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "BestRetentionTime", normalization = TRUE, ytitle = "CUSUMv Chart : retention time", type = "variability", referenceValue = 0.5, decisionInterval = 5, selectMean = NULL, selectSD = NULL)
# A CUSUMv chart when a guide set (1-20 runs) is used to monitor the variability of total peak area
CUSUMChart(data, peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "TotalArea", normalization = TRUE, ytitle = "CUSUMv Chart : peak area", type = "variability", referenceValue = 0.5, decisionInterval = 5, selectMean = NULL, selectSD = NULL)
```
# `MSstatsQC` core functions: `ChangePointEstimator()`

Follow-up change point analysis is helpful to identify the time of a change for each peptide and metric. `ChangePointEstimator()` function is used for the analysis. This function is one of the core functions and uses the same arguments. We recommend using this function after control charts generate an out-of-control observation. For example, retention time of TAAYVNAIEK increases over time as CUSUMm statistics increases steadily after the 20th time point. User can follow-up with `ChangePointEstimator()` function to find the exact time of retention time drift.

### Arguments

* `data`: comma-separated (.csv), metric file. It should contain a "Precursor" column and the metrics columns. It should also include "Annotations" for each observation.
* `peptide`: the name of precursor of interest.
* `L`: lower bound of the guide set.
* `U`: upper bound of the guide set.
* `metric`: the name of metric of interest.
* `normalization`: TRUE if data is standardized.
* `ytitle`: the y-axis title of the plot. The x-axis title is by default "Time : name of peptide"
* `type`: the type of the control chart. Two values can be assigned, "mean" or "dispersion". Default is "mean"
* `selectMean`: the mean of a metric. It is used when mean is known. It is NULL when mean is not known.  The default is NULL.
* `selectSD`: the standard deviation of a metric. It is used when standard deviation is known. It is NULL when mean is not known. The default is NULL.

### Example

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Retention time >> first 20 observations are used as a guide set
XmRChart(data, "TAAYVNAIEK", metric = "BestRetentionTime", type = "mean", L = 1, U = 20)
ChangePointEstimator(data, "TAAYVNAIEK", metric = "BestRetentionTime", type = "mean", L = 1, U = 20)
```

We don't recommend using this function when all the observations are within control limits. In the case of retention time monitoring of LVNELTEFAK, there is no need to further analyse change point.

The time of a variability change can be analyzed with the same fucntion. For example, retention time of YSTDVSVDEVK experiences a drift in the mean of retention time and variability of retention time increases simultaneously. In this case, `ChangePointEstimator()` can be used to identify exact times of both changes.

### Example

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Retention time >> first 20 observations are used as a guide set
XmRChart(data, "YSTDVSVDEVK", metric = "BestRetentionTime", type = "mean", L = 1, U = 20)
ChangePointEstimator(data, "YSTDVSVDEVK", metric = "BestRetentionTime", type = "variability", L = 1, U = 20)
```

# `MSstatsQC` summary functions: river and radar plots

`RiverPlot()` and `RadarPlot()` functions are the summary functions used in `MSstatsQC`. They are used to aggregate results over all analytes for X and mR charts or CUSUMm and CUSUMv charts. `method` argument is used to define the method where the results for multiple peptides are aggregated. For example, if user would like to aggregate information gathered from the X charts of retention time for all analytes, upper panel of `RiverPlot()` show for the increases and decreases in retention time. Next, `RadarPlot()` are used to find out which peptides are affected by the problem. 

If the mean and standard deviation is known, summary functions uses `listMean` and `listSD` arguments. For example, if user monitors retention time and peak assymetry and mean and standard deviations of these metrics are known, arguments will require entering a vector for means and another vector for standard deviations. 

### Example

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Retention time >> first 20 observations are used as a guide set
RiverPlot(data = S9Site54, L = 1, U = 20, method = "XmR")
RiverPlot(data = S9Site54, L = 1, U = 20, method = "CUSUM")
RadarPlot(data = S9Site54, L = 1, U = 20, method = "XmR")
RadarPlot(data = S9Site54, L = 1, U = 20, method = "CUSUM")
```

# `MSstatsQC` summary functions: decision map

`DecisionMap()` functions another summary function used in `MSstatsQC`. It is used to compare aggregated results over all analytes for a certain method such as XmR charts with the user defined criteria. Firstly, user defines the performance criteria and run `DecisionMap()` function to visualize overall performance. This function uses all the arguments of summary plots listed previously. Additionally, the following arguments are used

### Arguments

* `method`: the name of the method prefered. It is either "CUSUM" or "XmR"` interest.
* `peptideThresholdRed`: a threshold that marks percentage of out-of-control peptides. if the percentage is above this threshold, the color is red meaning fail. Default is 0.7.
* `peptideThresholdYellow`: a threshold that marks percentage of out-of-control peptides. if the percentage within this threshold and `peptideThresholdRed`, the color is yellow meaning warning. Default is 0.5.

### Example

```{r, eval=TRUE,fig.width=10, fig.height=5}
# A decision map for Site 54 can be generated using the following script
# Retention time >> first 20 observations are used as a guide set
DecisionMap(data,
    method = "XmR", peptideThresholdRed = 0.25, peptideThresholdYellow = 0.10,
    L = 1, U = 20, type = "mean", title = "Decision map", listMean = NULL, listSD = NULL
)
```

# Use case: longitudinal profiling of DDA with missing values

We analyzed the QCloudDDA dataset. The dataset had many missing values and MSstatsQC 2.0 processed these missing values and generated control charts and summary plots for longitudinal performance assessment.

```{r, eval=TRUE, fig.width=8, fig.height=5}
mydata <- DataProcess(MSstatsQC::QCloudDDA)
# Creating a missing data map
MissingDataMap(mydata)
XmRChart(mydata, "EACFAVEGPK", metric = "missing", type = "mean", L = 1, U = 15)
mydata <- RemoveMissing(mydata)
RiverPlot(mydata[, -9], L = 1, U = 15, method = "XmR")
RadarPlot(mydata[, -9], L = 1, U = 15, method = "XmR")
```

```{r, eval=TRUE, fig.width=8, fig.height=5}
mydata <- DataProcess(MSstatsQC::QCloudDDA)
# Creating a missing data map
MissingDataMap(mydata)
```

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Creating an X chart for missing counts
XmRChart(mydata, "EACFAVEGPK", metric = "missing", type = "mean", L = 1, U = 15)
```

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Removing missing values and analyzing the data
mydata <- RemoveMissing(mydata)
RiverPlot(mydata[, -9], L = 1, U = 15, method = "XmR")
RadarPlot(mydata[, -9], L = 1, U = 15, method = "XmR")
```
# Use case: longitudinal profiling of QC with iRT peptides

We analyzed the QuiCDIA dataset which included longitudinal profiles of 11 iRT peptides. The data comprised two DIA experiments acquired in duplicate with 11 days in between the two measurement sequences. 

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Checking missing values and analyzing the data
MissingDataMap(MSstatsQC::QuiCDIA)
RiverPlot(data = QuiCDIA, L = 1, U = 20, method = "XmR")
RadarPlot(data = QuiCDIA, L = 1, U = 20, method = "XmR")
```

# Use case: longitudinal profiling of QC from an SRM experiment

Following the implementation in Dogu et al. (2017), we analyzed the QCloudSRM dataset. This dataset was previously evaluated by the experts as well-performing for all peptides and metrics, except for the outlying peptide TCVADESHAGCEK.

```{r, eval=TRUE, fig.width=8, fig.height=5}
# Checking missing values and analyzing the data
MissingDataMap(MSstatsQC::QCloudSRM)
RiverPlot(data = QCloudSRM, L = 1, U = 20, method = "CUSUM")
RadarPlot(data = QCloudSRM, L = 1, U = 20, method = "CUSUM")
```

# `MSstatsQC-ML` functions: MSstatsQC-ML.trainR() and MSstatsQC.deployR()

`MSstatsQC-ML.trainR()` and `MSstatsQC.deployR()` functions are the two main functions to use machine learning approach proposed in this package. The two packages needs to be used in an order. The user first needs to train a model with the first function (trainR) and than use this function to predict labels per QC run using the second function (deployR). The first function returns a trained RF model from the guide set automating split, feature engineering, hypertuning and finalizing training based on train/validation/test datasets. It requires a simulation size to generate synthetic runs with random artifacts of deterioration. Deployment function labels each peptide per run with a probability of failure output that is defined to evaluate quality of a run.   

### Example

```{r, eval=TRUE, fig.width=8, fig.height=5}

if (requireNamespace("h2o", quietly = TRUE)) {
  tryCatch({
    h2o::h2o.init(nthreads = -1)
  }, error = function(e) {
    warning("H2O did not start.")
  })
}

utils::data("S9Site54", package = "MSstatsQC", envir = environment())

S9Site54.dataML <- DataProcess(S9Site54[, -c(4, 5, 8)])
colnames(S9Site54.dataML)[1] <- c("idfile")
colnames(S9Site54.dataML)[2] <- c("peptide")
S9Site54.dataML$peptide <- as.factor(S9Site54.dataML$peptide)
S9Site54.dataML$idfile <- as.numeric(S9Site54.dataML$idfile)
S9Site54.dataML <- within(S9Site54.dataML, rm(Annotations, missing))
guide.set <- dplyr::filter(S9Site54.dataML, idfile <= 20)
guide.set <- as.data.frame(guide.set)
Test.set <- S9Site54.dataML

rf_model = MSstatsQC::MSstatsQC.ML.trainR(guide.set, sim.size = 10)

results <- MSstatsQC.ML.deployR(Test.set, guide.set, rf_model)

print(results$DecisionMap)

for (p in results$InterpretationPlots) {
    print(p)
}

```

# Output

Plots created by the core plot functions are generate by `plotly` which is an R package for interactive plot generation. These interactive plots created by `MSstatsQC`, can be saved as an html file using the save widget function. If the user wants to save a static png file, then `export` function can be used. The outputs of other MSstatsQC functions are generated by `ggplot2` package and saving those outputs would require using `ggsave` function.

### Example

```{r, eval=FALSE}
# Saving plots generated by plotly
p <- XmRChart(data,
    peptide = "TAAYVNAIEK", L = 1, U = 20, metric = "BestRetentionTime", normalization = FALSE,
    ytitle = "X Chart : retention time", type = "mean", selectMean = NULL, selectSD = NULL
)
htmlwidgets::saveWidget(p, "Aplot.html")
export(p, file = "Aplot.png")

# Saving plots generated by ggplot2
p <- RiverPlot(data, L = 1, U = 20)
ggsave(filename = "Summary.pdf", plot = p)
# or
ggsave(filename = "Summary.png", plot = p)
```

## Project website
Please use [MSstats.org/MSstatsQC](https://www.MSstats.org/MSstatsQC) and [github repository](https://github.com/eralpdogu/MSstatsQC) and [github page](https://eralpdogu.github.io/MSstatsQC/) for further details about this tool.

## Question and issues
Please use [Google group](https://groups.google.com/forum/#!forum/msstatsqc) if you want to file bug reports or feature requests.

# Citation

Please cite MSstatsQC: 

* Dogu, E., Mohammad-Taheri, S., Abbatiello, S. E., Bereman, M. S., MacLean, B., Schilling, B., & Vitek, O. (2017). MSstatsQC: Longitudinal system suitability monitoring and quality control for targeted proteomic experiments. Molecular & Cellular Proteomics, 16(7), 1335-1347. https://doi.org/10.1074/mcp.M116.064774
* Dogu, E., Taheri, S. M., Olivella, R., Marty, F., Lienert, I., Reiter, L., ... & Vitek, O. (2019). MSstatsQC 2.0: R/Bioconductor package for statistical quality control of mass spectrometry-based proteomics experiments. Journal of proteome research, 18(2), 678-686. https://doi.org/10.1021/acs.jproteome.8b00732 
* Dogu, E., Gupta, S, Olivella, R., Sabido, E. & Vitek, O. (2026). MSstatsQC-ML: A supervised machine learning approach to monitor system suitability and quality control in mass spectrometry-based proteomics. Under Review.

# Session information

```{r si}
sessionInfo()
```



