MSstatsQC: Longitudinal System Suitability Monitoring and Quality Control for Targeted Proteomic Experiments

Overview and Notation

Denote i = 1,2, …, I the index of a metric (such as retention time or peak area), j = 1,2, …, J the index of a peptide, and t = 1,2, …, T the time point that an SST or QC run is acquired. Denote X_ijt the value of the metric for this peptide and time point. Because most discussion below focuses on one metric and one peptide at a time, we will omit the indices of metrics and peptides for simplicity, and denote the values in a time point as X_t.

Because X_t is subject to stochastic variation, it is a random variable. Denote μ and σ the expected value and the standard deviation of X_t. It is also useful to define the quantity Z_t = , called Z score, which standardizes X_t. The Z score is independent of the measurement units of the metric. It fluctuates around the expected value of 0, and has the standard deviation of 1. Note that μ is a standard deviation and not a standard error, and therefore the Z score is different from a statistic used, e.g. in a t test. The true values of μ and σ are typically unknown, and their estimates μ̂ and σ̂ are obtained from a guide set of high quality runs.

Example Data Sets for Longitudinal QC Monitoring: Simulated QC

Three simulated data sets were designed to evaluate longitudinal QC monitoring algorithms for SRM proteomic workflows in situations where the true nature of the problem is known. The simulations were based on the time course QC profiles of proteins in blood presented in tutorials for SProCop (20) and Panorama AutoQC (21). Ten measurements (runs 5–15) were used to estimate peak area mean (μ̂) and standard deviation (σ̂) of peptide VYVEELKPTPEGDLEILLQK. This peptide was stable in these runs, and produced good quality data (20). Next, peak areas measurements at 50 time points were simulated from Normal distribution X_t ∼ N(μ̂, σ̂²). Finally, three types of disturbances were created. The Simulated Data Set 1 introduced three outlying observations with a larger mean level of peak area by simulating X_t∼N(μ̂ + 3σ̂, σ̂²) at t = 30, 40, 50. The Simulated Data Set 2 introduced a small systematic increase (i.e. step shift) in peak area measurements by simulating X_t∼N(μ̂ + σ̂, σ̂²) at 26 ≤ t ≤ 50. The Simulated Data Set 3 imitated a drift in peak areas, and introduced an incremental increase by simulating X_t ∼ N(μ̂ + βσ̂(t − 25), σ̂²) at 26 ≤ t ≤ 50, with β = 0.1. We used these data sets to illustrate the ability of SPC to detect isolated outliers and subtle measurement drifts.

Example Data Sets for Longitudinal SST Monitoring: CPTAC Study 9.1

CPTAC study 9.1 (sites 54, 56A, 65, and 86) (8) were used to evaluate longitudinal SST monitoring algorithms for SRM proteomic workflows. The study designed a reference sample with 15 peptides to characterize instrument performance in terms of peak area, retention time and FWHM metrics. The study acquired guide sets from this sample, as well as longitudinal measurements collected during two months in 11 labs. The number of time points in a lab varied between 36 and 89. The data sets were processed with Skyline 3.6.1.10690 to generate a table of quantitative metrics for input to MSstatsQC. We used these data sets to evaluate the ability of SPC to discover abnormalities such as retention time and peak area drifts (Site 54), retention time fluctuation (Site 86), and chromatographic shape problems (Site 65), and to illustrate effects of guide set selection on the performance of detecting a retention time drift (Site 56A). Guide sets are manually determined from in-control observations for each site.

Previously Used SPC Methods in Proteomics: Control Charts for Monitoring Mean

The fundamental SPC tool for monitoring is a control chart, first proposed by Shewhart (13, 14), and later transferred to clinical chemistry by Levey and Jennings (15). A control chart summarizes QC measurements during time (e.g. peak area in Fig. 1A and 1B). Decisions from the control charts are made by comparing the observations to control limits. In some cases, the control limits can be pre-defined by industry standards for the purposes of SST, e.g. “mass accuracy must be less than 2 ppm” or “retention time measurements must be within 2 min of the target retention time.” In many other cases, however, the control limits are empirically derived, as we discuss below.

To make conclusions on the original (i.e. not standardized) scale of the metric, the X chart (Fig. 1A) monitors the mean of X_t and sets control limits in the form (μ̂ ± Lσ̂). This expression assumes that the observations are independent and Normally distributed. A typical value of L is 3, which corresponds to ±3 standard deviations from the mean of a metric, and results in the false alarm probability α = 0.0027. This approach is used in SProCop (20). Equivalently, the Z chart (Fig. 1B) monitors the standardized measurements Z_t = , and sets the control limits to ±3. If an observation is within the control limits, it is considered in-control. If an observation exceeds the limits, it is out-of-control. The goal of the next steps is to find the earliest time point when the variation occurred, and investigate and eliminate its root cause.

SPC Method Advocated In This Manuscript: Simultaneous Monitoring of Mean and Variability With XmR and ZmR Charts

Although X charts can detect shifts the mean value of a metric, they can fail detecting shifts in its variation. Yet shifts in variation can be a major issue in mass spectrometric experiments, e.g. in cases of undesirable between-run fluctuation in retention time of an analyte. Monitoring variation of a metric along with its mean allows us to detect these problems. To this end, the X chart is supplemented with another control chart that displays moving ranges mR_t = |X_t+1 − X_t|, t ≥ 2 (Fig. 1C) (16). Values of mR are sensitive to changes between consecutive measurements, and therefore to changes in the variability of a metric. Control limits for the mR chart are set to (0, Lm̅R̅). The constant L is often set to L = 3.268, matching the false alarm probability α = 0.0027 of the X and Z charts for moving ranges of length 2. The value of m̅R̅ is estimated as the average of all the moving ranges in the guide set.

A pair (X chart, mR chart) of a metric is called XmR chart, and a pair (Z chart, mR chart) is called ZmR chart. In XmR and ZmR charts the estimate of variation σ̂ is typically replaced to m̅R̅, and the constant L in the control limits of the X chart is set to L = 2.66, to match the false alarm probability α = 0.0027 (16). The XmR and ZmR charts detect observations that are out-of-control in terms of both the mean of the metric, and its variation.

SPC Methods Advocated in This Manuscript: Detection of Subtle Shifts With CUSUM_mand CUSUM_v

Although XmR and ZmR charts detect large undesirable deviations, they miss subtler problems, such as small isolated deviations of values of a metric, or slow sustained drifts. As an attempt to remedy this situation, Westgard rules (17) are used. These rules are used as plug-ins to X control chart. For example, Rule 10x raises a red flag when 10 consecutive retention time measurements fall on one side of their respective μ̂ for a peptide. Unfortunately, these rules do not meet the needs of mass spectrometry-based proteomics. They are somewhat arbitrary, and do not scale beyond a small number of metrics and analytes. The differences in failure times and types of the problems between the analytes complicate the use of these rules even further. We advocate another approach that uses a family of cumulative sum (CUSUM) control charts (33) detects such subtle shifts effectively by monitoring time-weighted measurements.

The CUSUM_m chart (16, 34) (Fig. 2C) monitors shifts in the metric's mean. It has two components, CUSUM_m+ and CUSUM_m− that separately consider positive and negative changes. For example, a slow but steady increase in retention time is reflected by an increased CUSUM_m+. For at time point (t ≥ 1), CUSUM_m+,t and CUSUM_m−,t are defined as cumulative sums of sequential measurements, i.e. CUSUM_m+,t = max (0,X_t − k + CUSUM_{m+,(t − 1)}), and CUSUM_m−,t = max(0,−X_t − k + CUSUM_{m−,t − 1}), with starting values CUSUM_m+,0 = CUSUM_m−,0 = 0. Here k is a parameter (also called slack value), which affects the sensitivity of the approach. The control limit h is defined such that P(0 < CUSUM_m+ < h) = P(0 < CUSUM_m− < h) = 1 − . If CUSUM_m+,t or CUSUM_m−,t exceed the control limit h, the process is considered out-of-control.

Here we advocate a modified version of CUSUM that simplifies the interpretation. Specifically, we replace X_t with its standardized version Z_t, and set k = 0.5 and h ≅ 5. In this definition, CUSUM_m+,t = max (0,Z_t − 0.5 + CUSUM_{m+,(t − 1)}), and CUSUM_m−,t = max (0,−Z_t − 0.5 + CUSUM_{m−,t − 1}). It can be shown that this approach detects 1σ deviations from the metric mean, with the false alarm rate of α = 0.0027 that matches the approaches above.

Similarly, CUSUM_ν charts (Fig. 2D) detect slow changes in variation. For example, if retention time becomes increasingly unstable, the value of CUSUM_ν+ increases. CUSUM_ν+ is constructed from the transformation V_t = (16, 34). Hawkins (1981) (35) showed that is very close to Normal distribution with mean 0.822 and standard deviation 0.349, and therefore V_t follows approximately standard Normal distribution. Similarly to CUSUM_m, define CUSUM_ν+,t = max (0, V_t − k + CUSUM_ν+,(t−1)), and CUSUM_ν−,t = max(0,−V_t − k + CUSUM_ν−,t−1) We suggest the same parameters k = 0.5 and h ≅ 5 to match α = 0.0027. If CUSUM_ν+,t or CUSUM_ν−,t exceeds the control limit h, the variation of the process is out-of-control.

SPC Method Advocated in this Manuscript: Change Point Analysis

The goal of longitudinal profiling is to detect the earliest time point, called change point, at which a problem occurs. However, in some cases, especially in presence of slow sustained drifts, the problem can only be detected with a delay. Change point analysis, performed after an out-of-control observation occurs, helps pinpoint the change point more exactly (e.g. supplemental Fig. S2). Many statistical approaches for change point analysis exist (36).

Here we advocate detecting change point based on maximum likelihood estimation of (1) a step shift change model for the mean (37), and (2) a step shift change model for the variability (38). This approach has been proven effective in conjunction with control charts (37⇓–39). Assume that a step change in the mean of a metric occurs after the change point τ (τ = 1, …, T). Then, where N(.) denotes Normal distribution, μ₁ = μ + δσ is the post-change mean and δ is the magnitude of a change. The estimate of the change point τ̂ is obtained by maximizing the likelihood L(τ,μ₁|X) = . This is equivalent to maximizing C_t = (T − t + 1) (X̄_T,t − μ̂)² where X̄_T,t =Σ_t′=t^T X_t′/(T − t′ + 1) is the reverse cumulative mean. The change point estimate is therefore τ̂ = argmax_1≤t<T (C_t).

Similarly, the change point for variation is estimated by maximizing the likelihood L(τ,σ₁|X) = , or, equivalently, maximizing D_t = SS_T,t/2σ̂² − ((T − t + 1)/2)log_e(SS_T,t/((T − t + 1)σ̂²)) − (T − t + 1)/2, where SS_T,t = Σ_t′=t^T (X_t′ − μ̂)² is the reverse cumulative sum of squares of T − t most recent observations. The change point estimate is τ̂ = argmax_1≤t<T (D_t) The change point estimate for variability helps user identify the time of a change in the variability of a metric.

SPC Method Advocated in this Manuscript: Summary Plots and Decision Maps

A distinctive feature of mass spectrometric experiments is the ability to simultaneously quantify multiple metrics of multiple analytes, for the purposes of biological investigations, but also for the purposes of SST or QC. However, extending the control charts to highly multivariate settings is non-trivial. The metrics and analytes produce multiple charts, which exponentially increase the comparisons to the control limits, and inflate the number of false alarms.

Given the complexity of the situation, we advocate a visual representation of individual longitudinal profiles of metrics and analytes, as well as simple ad hoc summaries (pass, warning, fail) of the overall performance, as described below.

A decision map (e.g. Fig. 4A and 4B) presents the most general overview of the overall performance per method (e.g. X chart and mR chart). It takes as input user-defined pass/warning/fail criteria for a control chart, and aggregates the conclusions from each metric over all the analytes and summaries obtained for the method. For example, the user-defined criteria can generate a warning signal if the proportion p_t of analytes with out-of-control values in at least one metric (e.g. retention time) is 0.10 ≤ p_t < 0.25, and generate a fail signal if p_t ≥ 0.25. With this approach, each metric generates its own pass/warning/fail signals at each time point for each control chart. The signals are visually summarized, e.g. with a heatmap. Although this approach does not guarantee optimal decisions, it is a step toward making more objective, reproducible, and documentable decisions in multivariate situations.

A more detailed Look into the results of the control charts is required to efficiently identify the root causes of the problems. A river plot and a radar plot provide such more detailed insight. The river plot (e.g. Fig. 4C) details the performance of each metric in terms of the mean and variability of a metric, aggregated over the analytes. For example, in the case of retention time or total peak area in Fig. 4C, the river plots distinguish the proportions of out-of-control peptides with positive and negative deviations in mean and variation (colored lines), present trends over time (x axis), and indicate the estimated change point time for each deviation type and each analyte (red dots). Sloping curves and large numbers of red dots indicate worsening of the performance.

A radar plot (e.g. Fig. 4D) investigates the extent of the problems at the peptide level. For each metric, each summary and each peptide, it displays the number of time points that exceeded the control limits throughout the duration of the experiment. Different types of deviations (i.e. increases and decreases in the mean or variability) are color-coded. Large colored areas indicate worsening of the performance.

Implementation in Open-source Software

We implemented the methods discussed in this manuscript in an R package MSstatsQC. MSstatsQC is available free and open-source under Artistic 2.0 license (same license as, e.g. the project Bioconductor (40)). The package takes as input SST or QC measurements in a tabular format produced by data processing tools such as Skyline, or any other. The package can be accessed through its online graphical user interface, where input files of at most 5MB are uploaded on a remote server hosted by RStudio. Alternatively, the package can be installed on a local server, and used via its web interface or via a command line. The size of the input files can be increased to 100MB in that case. This is useful for organizations that wish to perform all the analyses locally, or integrate SPC into broader automated pipelines. The implementations, examples, experimental data sets, link to the source code on Github are available at msstats.org/msstatsqc.