Filter features based on feature annotations

Below are some examples of criteria that are used to identify low-confidence. The information is readily available since this was computed by MaxQuant:

• Remove PSMs that are matched to contaminants
• Remove PSMs that are matched to the decoy database
• Keep PSMs that exhibit a high PIF (parental ion fraction), indicative of the purity of a spectrum

We can perform this filtering using the filterFeatures function from QFeatures. filterFeatures automatically accesses the feature annotations and selects the rows that meet the provided condition(s). For instance, Reverse != "+" keeps the rows for which the Reverse variable in the rowData is not "+" (i.e. the PSM is not matched to the decoy database).

scp <- filterFeatures(scp,
                      ~ Reverse != "+" &
                          Potential.contaminant != "+" &
                          !is.na(PIF) & PIF > 0.8)
#> 'Reverse' found in 4 out of 4 assay(s)
#> 'Potential.contaminant' found in 4 out of 4 assay(s)
#> 'PIF' found in 4 out of 4 assay(s)

Filter assays based on detected features

To avoid proceeding with failed runs, another interesting filter is to remove assays with too few features. If a batch contains less than, for example, 150 features we can then suspect something wrong happened in that batch and it should be removed. Using dims, we can query the dimensions (hence the number of features and the number of samples) of all assays contained in the dataset.

dims(scp)
#>      190222S_LCA9_X_FP94BM 190321S_LCA10_X_FP97AG 190321S_LCA10_X_FP97_blank_01
#> [1,]                   283                    318                            60
#> [2,]                    11                     11                            11
#>      190914S_LCB3_X_16plex_Set_21
#> [1,]                          200
#> [2,]                           16

Actually, a QFeatures object can be seen as a three-order array: features × samples × assay. Hence, QFeatures supports three-order subsetting x[rows, columns, assays]. We first select the assays that have sufficient PSMs (the number of rows is greater than 150), and then subset the scp object for the assays that meet the criterion.

keepAssay <- dims(scp)[1, ] > 150
scp <- scp[, , keepAssay]
#> Warning: 'experiments' dropped; see 'drops()'
#> harmonizing input:
#>   removing 11 sampleMap rows not in names(experiments)
#>   removing 11 colData rownames not in sampleMap 'primary'
scp
#> An instance of class QFeatures containing 3 assays:
#>  [1] 190222S_LCA9_X_FP94BM: SingleCellExperiment with 283 rows and 11 columns 
#>  [2] 190321S_LCA10_X_FP97AG: SingleCellExperiment with 318 rows and 11 columns 
#>  [3] 190914S_LCB3_X_16plex_Set_21: SingleCellExperiment with 200 rows and 16 columns

Notice the 190321S_LCA10_X_FP97_blank_01 sample was removed because it did not contain sufficient features, as expected from a blank run. This could also have been the case for failed runs.

Filter features based on SCP metrics

Another type of filtering is specific to SCP. In the SCoPE2 analysis, the authors suggest a filters based on the sample to carrier ratio (SCR), that is the reporter ion intensity of a single-cell sample divided by the reporter ion intensity of the carrier channel (200 cells) from the same batch. It is expected that the carrier intensities are much higher than the single-cell intensities.

The SCR can be computed using the computeSCR function from scp. The function must be told which channels are the samples that must be divided and which channel contains the carrier. This information is provided in the sample annotations and is accessed using the colData, under the SampleType field.

table(colData(scp)[, "SampleType"])
#> 
#>      Blank    Carrier Macrophage   Monocyte  Reference     Unused 
#>          3          3         20          5          3          4

In this dataset, SampleType gives the type of sample that is present in each TMT channel. The SCoPE2 protocole includes 5 types of samples:

• The carrier channels (Carrier) contain 200 cell equivalents and are meant to boost the peptide identification rate.
• The normalization channels (Reference) contain 5 cell equivalents and are used to partially correct for between-run variation.
• The unused channels (Unused) are channels that are left empty due to isotopic cross-contamination by the carrier channel.
• The negative control channels (Blank) contain samples that do not contain any cell but are processed as single-cell samples.
• The single-cell sample channels contain the single-cell samples of interest, that are macrophage (Macrophage) or monocyte (Monocyte).

The computeSCR function expects the following input:

• The QFeatures dataset
• The assay name(s) or index/indices for which the SCR should be computed
• colvar: the variable in the sample annotations (colData) that hold the information used to discriminate sample channels from carrier channels.
• carrierPattern: a string pattern (following regular expression syntax) that identifies the carrier channel in each batch.
• samplePattern: a string pattern (following regular expression syntax) that identifies the samples to divide.

Optionally, you can also provide the following arguments:

• rowDataName: the name of the column in the rowData where to store the computed SCR for each feature.
• sampleFUN: when multiple samples are present in an assay, there are as many SCR as there are samples that need to be summarized to a single value per feature. sampleFUN tells which function to use for summarizing the sample values before computing the SCR; the default is the mean.
• carrierFUN: some designs might include several carriers per run (not the case in this example). Similarly to sampleFUN, carrierFUN tells which function to use for summarizing the carrier values before computing the SCR; the default is the same function as sampleFUN.

The function creates a new field in the rowData of the assays. We compute the average SCR for each PSM and store it in the corresponding rowData, under the MeanSCR column.

scp <- computeSCR(scp,
                  i = 1:3,
                  colvar = "SampleType",
                  carrierPattern = "Carrier",
                  samplePattern = "Macrophage|Monocyte",
                  sampleFUN = "mean",
                  rowDataName = "MeanSCR")

Before applying the filter, we plot the distribution of the average SCR. We collect the rowData from several assays in a single table DataFrame using the rbindRowData function from QFeatures.

rbindRowData(scp, i = 1:3) |>
    data.frame() |>
    ggplot(aes(x = MeanSCR)) +
    geom_histogram() +
    geom_vline(xintercept = c(1/200, 0.1),
               lty = c(2, 1)) +
    scale_x_log10()

The expected ratio between single cells and the carrier is 1/200 (dashed line). We can see that the distribution mode is slightly shifted towards higher ratios with a mode around 0.01. However, there are a few PSMs that stand out of the distribution and have a much higher signal than expected, indicating something wrong happened during the quantification of those PSMs. We therefore filter out PSMs with an average SCR higher than 0.1 (solide line). This is again easily performed using the filterFeatures functions.

scp <- filterFeatures(scp,
                      ~ !is.na(MeanSCR) &
                          MeanSCR < 0.1)
#> 'MeanSCR' found in 3 out of 3 assay(s)

Filter features to control for FDR

Finally, we might also want to control for false discovery rate (FDR). MaxQuant already computes posterior error probabilities (PEP), but filtering on PEPs is too conservative (Käll et al. (2008)) so we provide the pep2qvalue function to convert PEPs to q-values that are directly related to FDR. We here compute the q-values from the PEP (dart_PEP) across all 3 assays. dart_PEP contains the PEP values that have been updated using the DART-ID algorithm (Chen, Franks, and Slavov (2019)). The function will store the results in the rowData, we here asked to name the new column qvalue_PSMs.

scp <- pep2qvalue(scp,
                  i = 1:3,
                  PEP = "dart_PEP",
                  rowDataName = "qvalue_PSMs")

We also allow to compute q-values at peptide or protein level rather than PSM. In this case, you need to supply the groupBy argument. Suppose we want to compute the q-values at protein level, we can fetch the protein information stored under Leading.razor.protein in the rowData. This time, we store the q-values in a new field called qvalue_proteins.

scp <- pep2qvalue(scp,
                  i = 1:3,
                  PEP = "dart_PEP",
                  groupBy = "Leading.razor.protein",
                  rowDataName = "qvalue_proteins")

We can now filter the PSM to control, let’s say, the protein FDR at 1%. This can be performed using filterFeatures because the q-values were stored in the rowData.

scp <- filterFeatures(scp,
                      ~ qvalue_proteins < 0.01)
#> 'qvalue_proteins' found in 3 out of 3 assay(s)

Relative reporter ion intensity

In order to partialy correct for between-run variation, SCoPE2 suggests computing relative reporter ion intensities. This means that intensities measured for single-cells are divided by the reference channel containing 5-cell equivalents. We use the divideByReference function that divides channels of interest by the reference channel. Similarly to computeSCR, we can point to the samples and the reference columns in each assay using the annotation contained in the colData.

We here divide all columns (using the regular expression wildcard .) by the reference channel (Reference).

scp <- divideByReference(scp,
                         i = 1:3,
                         colvar = "SampleType",
                         samplePattern = ".",
                         refPattern = "Reference")

Filter samples of interest

We can subset the cells of interest, that is the negative control samples, the macrophages and the monocytes. This can easily be done by taking advantage of the colData and the subsetting operators. Recall that QFeatures objects support three-order subsetting, x[rows, columns, assays], where columns are the samples of interest.

scp
#> An instance of class QFeatures containing 7 assays:
#>  [1] 190222S_LCA9_X_FP94BM: SingleCellExperiment with 239 rows and 11 columns 
#>  [2] 190321S_LCA10_X_FP97AG: SingleCellExperiment with 277 rows and 11 columns 
#>  [3] 190914S_LCB3_X_16plex_Set_21: SingleCellExperiment with 186 rows and 16 columns 
#>  [4] peptides_190222S_LCA9_X_FP94BM: SingleCellExperiment with 216 rows and 11 columns 
#>  [5] peptides_190321S_LCA10_X_FP97AG: SingleCellExperiment with 256 rows and 11 columns 
#>  [6] peptides_190914S_LCB3_X_16plex_Set_21: SingleCellExperiment with 168 rows and 16 columns 
#>  [7] peptides: SingleCellExperiment with 384 rows and 38 columns
scp <- scp[, scp$SampleType %in% c("Blank", "Macrophage", "Monocyte"), ]

The subsetting removes unwanted samples from all assays. The filtered data set contains the same number of assays with the same number of features, but the number of columns (hence sampled) decreased.

scp
#> An instance of class QFeatures containing 7 assays:
#>  [1] 190222S_LCA9_X_FP94BM: SingleCellExperiment with 239 rows and 8 columns 
#>  [2] 190321S_LCA10_X_FP97AG: SingleCellExperiment with 277 rows and 8 columns 
#>  [3] 190914S_LCB3_X_16plex_Set_21: SingleCellExperiment with 186 rows and 12 columns 
#>  [4] peptides_190222S_LCA9_X_FP94BM: SingleCellExperiment with 216 rows and 8 columns 
#>  [5] peptides_190321S_LCA10_X_FP97AG: SingleCellExperiment with 256 rows and 8 columns 
#>  [6] peptides_190914S_LCB3_X_16plex_Set_21: SingleCellExperiment with 168 rows and 12 columns 
#>  [7] peptides: SingleCellExperiment with 384 rows and 28 columns

Filter based on the median relative intensity

We compute the median relative reporter ion intensity for each cell separately and apply a filter based on this statistic. This procedure recalls that of library size filtering commonly performed in scRNA-Seq data analysis, where the library size is the sum of the counts in each single cell. We compute and store the median intensity in the colData.

medians <- colMedians(assay(scp[["peptides"]]), na.rm = TRUE)
scp$MedianRI <- medians

Looking at the distribution of the median per cell can highlight low-quality cells.

colData(scp) |>
    data.frame() |>
    ggplot() +
    aes(x = MedianRI,
        y = SampleType,
        fill = SampleType) +
    geom_boxplot() +
    scale_x_log10()

The negative control samples should not contain any peptide information and are therefore used to assess the amount of background signal. The graph above confirms that the signal measured in single-cells (macrophages and monocytes) is above the background signal, hence no filtering is needed. Would it not be the case, the same procedure as in the previous section can be used for selecting the cells that have an associated median RI lower that a defined threshold.

Filter based on the median CV

The median CV measures the consistency of quantification for a group of peptides that belong to a protein. We remove cells that exhibit high median CV over the different proteins. We compute the median CV per cell using the computeMedianCV function from the scp package. The function takes the peptides assay and computes the CV for each protein in each cell. To perform this, we must supply the name of the rowData field that contains the protein information through the groupBy argument. We also only want to compute CVs if we have at least 5 peptides per protein. Finally, we also perform a normalization and divide the columns by the median. The computed median CVs are automatically stored in the colData under the name that is supplied, here MedianCV.

scp <- medianCVperCell(scp,
                       i = 1:3,
                       groupBy = "Leading.razor.protein",
                       nobs = 5,
                       norm = "div.median",
                       na.rm = TRUE,
                       colDataName = "MedianCV")

The computed CVs are stored in the colData of the peptides assay and holds the median CV per cell computed using at least 5 observations (peptides). The main interest of computing the median CV per cell is to filter cells with reliable quantification. The negative control samples are not expected to have reliable quantifications and hence can be used to estimate an empirical null distribution of the CV. This distribution helps defining a threshold that filters out single-cells that contain noisy quantification.

getWithColData(scp, "peptides") |>
    colData() |>
    data.frame() |>
    ggplot(aes(x = MedianCV,
               fill = SampleType)) +
    geom_boxplot() +
    geom_vline(xintercept = 0.65)

We can see that the protein quantification for single-cells are much more consistent than for negative control samples. Based on the distribution of the negative controls, we decide to keep the cells that have a median CV lower than 0.65. Note this example is inaccurate because the null distribution is based on only 3 negative controls, but more sets could lead to a better estimation of the CV null distribution.

scp <- scp[, !is.na(scp$MedianCV) & scp$MedianCV < 0.65, ]

We can now remove the negative controls since all QC metrics are now computed.

scp <- scp[, scp$SampleType != "Blank", ]

Normalization

The columns (samples) of the peptide data are first normalized by dividing the relative intensities by the median relative intensities. Then, the rows (peptides) are normalized by dividing the relative intensities by the mean relative intensities. The normalized data is stored in a separate assay. This normalization procedure is suggested in the SCoPE2 analysis and is applied using the sweep method. Beside the dataset and the assay to normalize, the method expects a MARGIN, that is either row-wise (1) or column-wise (2) transformation, the FUN function to apply and STATS, a vector of values to apply. More conventional normalization procedure can be found in ?QFeatures::normalize.

## Divide columns by median
scp <- sweep(scp,
             i = "peptides",
             MARGIN = 2,
             FUN = "/",
             STATS = colMedians(assay(scp[["peptides"]]), na.rm = TRUE),
             name = "peptides_norm_col")
## Divide rows by mean
scp <- sweep(scp,
             i = "peptides_norm_col",
             MARGIN = 1,
             FUN = "/",
             STATS = rowMeans(assay(scp[["peptides_norm_col"]]),  na.rm = TRUE),
             name = "peptides_norm")

Notice each call to sweep created a new assay. Let’s have a look to the current stat of the QFeatures plot:

plot(scp)

Remove peptides with high missing rate

Peptides that contain many missing values are not informative. Therefore, another common procedure is to remove higly missing data. In this example, we remove peptides with more than 99 % missing data. This is done using the filterNA function from QFeatures.

scp <- filterNA(scp,
                i = "peptides_norm",
                pNA = 0.99)

Log-transformation

In this vignette, we perform log2-transformation using the logTransform method from QFeatures. Other log-transformation can be applied by changing the base argument.

scp <- logTransform(scp,
                    base = 2,
                    i = "peptides_norm",
                    name = "peptides_log")

Similarly to sweep, logTransform creates a new assay in scp.

Normalization

Normalization is performed similarly to peptide normalization. We use the same functions, but since the data were log-transformed at the peptide level, we subtract by the statistic (median or mean) instead of dividing.

## Center columns with median
scp <- sweep(scp, i = "proteins",
             MARGIN = 2,
             FUN = "-",
             STATS = colMedians(assay(scp[["proteins"]]),
                                na.rm = TRUE),
             name = "proteins_norm_col")
## Center rows with mean
scp <- sweep(scp, i = "proteins_norm_col",
             MARGIN = 1,
             FUN = "-",
             STATS = rowMeans(assay(scp[["proteins_norm_col"]]),
                              na.rm = TRUE),
             name = "proteins_norm")

Imputation

The protein data contains a lot of missing values.

scp[["proteins_norm"]] |>
    assay() |>
    is.na() |>
    mean()
#> [1] 0.3291962

The average missingness in the proteins assay is around 25 %. Including more samples and hence more batches can increase the missingness up to 70 % as seen for the complete SCoPE2 dataset (Specht et al. (2021)). Whether imputation is beneficial or deleterious for the data will not be discussed in this vignette. But taking those missing value into account is essential to avoid artefacts in downstream analyses. The data imputation is performed using the K nearest neighbors algorithm, with k = 3. This is available from the impute() mehtod. More details about the arguments can be found in ?impute::impute.knn.

scp <- impute(scp,
              i = "proteins_norm",
              name = "proteins_imptd",
              method = "knn",
              k = 3, rowmax = 1, colmax= 1,
              maxp = Inf, rng.seed = 1234)
#> Imputing along margin 1 (features/rows).

Note that after imputation, no value are missing.

scp[["proteins_imptd"]] |>
    assay() |>
    is.na() |>
    mean()
#> [1] 0

Batch correction

A very important step for processing SCP data is to correct for batch effects. Batch effects are caused by technical variation occurring during different MS runs. Since only a small number of single-cells can be acquired at once, batch effects are unavoidable.

The ComBat() function from the sva package can be used to perform batch correction as it is performed in the SCoPE2 analysis. We do not claim that ComBat is the best algorithm for batch correcting SCP data and other batch correcting methods could be used using the same procedure.

We first extract the assay to process.

sce <- getWithColData(scp, "proteins_imptd")
#> Warning: 'experiments' dropped; see 'drops()'

Next, we need to provide a design matrix and the batch annotation to Combat. The design matrix allows to protect variables of interest, in our case SampleType.

batch <- sce$runCol
model <- model.matrix(~ SampleType, data = colData(sce))

We then load and call ComBat and overwrite the data matrix. Recall the data matrix can be accessed using the assay function.

library(sva)
assay(sce) <- ComBat(dat = assay(sce),
                     batch = batch,
                     mod = model)

Finally, we add the batch corrected assay to the QFeatures object and create the feature links.

scp <- addAssay(scp,
                y = sce,
                name = "proteins_batchC")
scp <- addAssayLinkOneToOne(scp,
                            from = "proteins_imptd",
                            to = "proteins_batchC")

For the last time, we plot the overview of the fully processed data set:

plot(scp)

PCA

PCA can be computed using the runPCA method. It returns a SingleCellExperiment object for which the dimension reduction results are stored in the reducedDim slot.

scp[["proteins_batchC"]] <- runPCA(scp[["proteins_batchC"]],
                                   ncomponents = 5,
                                   ntop = Inf,
                                   scale = TRUE,
                                   exprs_values = 1,
                                   name = "PCA")

The computed PCA can be displayed using the plotReducedDim function. The dimred arguments should give the name of the dimension reduction results to plot, here we called it PCA. The samples are colored by type of sample.

plotReducedDim(scp[["proteins_batchC"]],
               dimred = "PCA",
               colour_by = "SampleType",
               point_alpha = 1)

This is a minimalistic example with only a few plotted cells, but the original SCoPE2 dataset contained more than thousand cells.

UMAP

Similarly to PCA, we can compute a UMAP using the runUMAP method. Note however that the UMAP implementation requires a initialization, usually provided by PCA. The previous PCA results are used automatically when supplying dimred = "PCA" (PCA is the name of the dimension reduction result that we supplied in the previous section).

scp[["proteins_batchC"]] <- runUMAP(scp[["proteins_batchC"]],
                                    ncomponents = 2,
                                    ntop = Inf,
                                    scale = TRUE,
                                    exprs_values = 1,
                                    n_neighbors = 3,
                                    dimred = "PCA",
                                    name = "UMAP")

The computed UMAP can be displayed using the plotReducedDim function. The dimred arguments gives the name of the dimension reduction results to plot, here we called it UMAP. The samples are colored by type of sample.

plotReducedDim(scp[["proteins_batchC"]],
               dimred = "UMAP",
               colour_by = "SampleType",
               point_alpha = 1)

The UMAP plot is a very interesting plot for large datasets. A UMAP on this small example dataset is not useful but is shown for illustration.