Contents

1 Overview

ProBatchFeatures is an extension of proBatch based on QFeatures (1), which extends its functionality by storing each processing stage (for example, raw data → log-transformation → normalization → batch effect correction → …) as a new SummarizedExperiment assay, and logging every step. It maintains full compatibility with QFeatures functions, and facilitates quick benchmarking of data preprocessing pipelines.

This vignette demonstrates how to:

  1. Construct a ProBatchFeatures object from matrix-formatted quantification data.

  2. Apply and record standard preprocessing steps (filtering features with too many missing values, log-transformation, normalization, and batch-effect correction) in a chained workflow.

  3. Perform diagnostic checks, including PCA, t-SNE, and UMAP.

  4. Retrieve processed assay data and the corresponding preprocessing pipeline for benchmarking or reuse.

2 Setup

2.1 Installation

To install the latest version of proBatch package, use BiocManager package:

if (!requireNamespace("BiocManager", quietly = TRUE)) {
    install.packages("BiocManager")
}
BiocManager::install("proBatch")

2.2 Loading required packages

This vignette uses dplyr, tibble, and ggplot2 (alongside proBatch) for data manipulation and visualization.

library(proBatch)

library(dplyr)
library(tibble)
library(ggplot2)
library(QFeatures)

2.3 Loading and preparing the example dataset

As an example, we will use the E. coli dataset described in the FedProt paper (2):

data("example_ecoli_data", package = "proBatch")

The dataset consists of data from five centers that were measured and analyzed independently. For simplicity, data from all centers have been merged into a single dataset, and a subset containing 400 protein groups and their corresponding peptides is used in this example.

# Extracting data
all_metadata <- example_ecoli_data$all_metadata
all_precursors <- example_ecoli_data$all_precursors
all_protein_groups <- example_ecoli_data$all_protein_groups
all_precursor_pg_match <- example_ecoli_data$all_precursor_pg_match

# Keeping only essential
rm(example_ecoli_data)

3 Constructing a ProBatchFeatures object

The dataset contains two hierarchical levels. While additional levels can be added in a manner similar to QFeatures, this example focuses on the peptide and protein group levels. The all_precursor_pg_match dataframe stores the mapping between these levels. We begin by using the ProBatchFeatures() constructor to create an object representing the peptide-level data:

# Building from a (features x samples) matrix
pbf <- ProBatchFeatures(
    data_matrix = all_precursors, # matrix of peptide intensities (features x samples)
    sample_annotation = all_metadata, # data.frame of sample metadata
    sample_id_col = "Run", # column in metadata that matches sample IDs
    level = "peptide" # label this assay level as "peptide"
)

pbf # show basic info about the ProBatchFeatures object
#> An instance of class ProBatchFeatures (type: bulk) with 1 set:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 1463 rows and 118 columns 
#>   Processing chain: unprocessed data (raw)

Assay names follow the convention “::”, for example, “peptide::raw”, or “protein::median_on_log”.

Next, we add the protein-level data as a second assay in the same ProBatchFeatures instance. Internally, this is achieved by creating a SummarizedExperiment for the protein-level data and adding it to the existing object using QFeatures functionality. The same sample annotation should be used for the new assay to ensure correct column alignment:

# Adding proteins as a new level and mapping them to precursors
#    all_precursor_pg_match has columns: "Precursor.Id", "Protein.Ids"
pbf <- pb_add_level(
    object = pbf,
    from = "peptide::raw",
    new_matrix = all_protein_groups,
    to_level = "protein", # will name "protein::raw" by default
    mapping_df = all_precursor_pg_match,
    from_id = "Precursor.Id",
    to_id = "Protein.Ids",
    map_strategy = "as_is"
)

pbf
#> An instance of class ProBatchFeatures (type: bulk) with 2 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 1463 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 400 rows and 118 columns 
#>   Processing chain:
#>   [1] add_level(protein)_byVar
#>   Steps logged: 1 (see get_operation_log())

If a precursor (peptide) maps to multiple protein groups, map_strategy specifies how to resolve multimappers. Available options are "first", "longest", and "as_is". The first two strategies select a single ProteinID per peptide: "first" takes the first match in the mapping data frame, while "longest" prefers the match with the most semicolon-separated identifiers, then the longest string, then the first occurrence. "as_is" requires a one-to-one mapping between identifiers.

# Check
validObject(pbf) # should be TRUE
#> [1] TRUE
assayLink(pbf, "protein::raw") # verify link summary
#> AssayLink for assay <protein::raw>
#> [from:peptide::raw|fcol:ProteinID|hits:1463]

# Remove unused variables to clean up the workspace
rm(all_metadata, all_precursor_pg_match, all_precursors, all_protein_groups)

4 Processing pipeline with diagnostics

With the data now loaded into a ProBatchFeatures object with peptide- and protein-level assays, we can demonstrate an example of typical preprocessing workflow. This pipeline includes filtering low-quality features, applying a log2-transformation, median normalization, and batch effect correction.

4.1 Step 1 - filtering features with too many missing values

Features detected in only a small subset of samples are unreliable and may obscure systematic effects in subsequent analyses. To examine the proportion of missing values across datasets, pb_nNA() function is used:

Function pb_nNA() summarizes missingness per assay:

# pb_nNA returns a DataFrame with counts of NAs per feature/sample
na_counts <- pb_nNA(pbf)

# na_counts contains info about total number of NAs and % in the data:
na_counts[["nNA"]] %>% as_tibble()
#> # A tibble: 2 × 3
#>   assay          nNA   pNA
#>   <chr>        <int> <dbl>
#> 1 peptide::raw 78870 0.457
#> 2 protein::raw  9163 0.194

Here nNA is a number of missing values in the assay and pNA is its proportion. It is possible to inspect the per-feature and per-sample breakdown for each assay:

# check NAs per sample:
head(na_counts[["peptide::raw"]]$nNAcols) %>% as_tibble()
#> # A tibble: 6 × 4
#>   assay        name       nNA   pNA
#>   <chr>        <chr>    <int> <dbl>
#> 1 peptide::raw lab_A_S1   496 0.339
#> 2 peptide::raw lab_A_S2   483 0.330
#> 3 peptide::raw lab_A_S3   524 0.358
#> 4 peptide::raw lab_A_S4   525 0.359
#> 5 peptide::raw lab_A_S5   541 0.370
#> 6 peptide::raw lab_A_S6   555 0.379

For each assay the following tables are available:

print(names(na_counts[["peptide::raw"]]))
#> [1] "nNA"     "nNArows" "nNAcols"

Visualisations help explore batch-specific patterns of missingness. Firstly, we define a color scheme which will be used in all diagnostic plots:

# color scheme for lab and condition labels
color_scheme <- sample_annotation_to_colors(
    pbf,
    sample_id_col = "Run",
    factor_columns = c("Lab", "Condition"),
    numeric_columns = NULL
)

Displaying features that contain at least one missing as a heatmap:

plot_NA_heatmap(
    pbf, # by default the last created assay
    show_rownames = FALSE, show_row_dend = FALSE,
    color_by = "Lab"
)

To display all features, set drop_complete = FALSE. When the number of features or samples exceeds 5000, a random subset of 5000 rows or columns is used for visualization by default. To display the entire dataset, set use_subset = FALSE.

To display multiple assays in one plot, specify them using pbf_name parameter:

plot_NA_heatmap(
    pbf,
    pbf_name = c("peptide::raw", "protein::raw"),
    show_rownames = FALSE, show_row_dend = FALSE,
    color_by = "Lab"
)

The overlap between peptide and protein identifications across samples can be visualized as follows:

plot_NA_frequency(
    pbf,
    pbf_name = c("peptide::raw", "protein::raw"),
    show_percent = FALSE
)

In the resulting plot, particularly for peptides, batch-specific patterns are evident, reflecting the 23-24 samples contributed by each laboratory.

Here, we apply a simple filter that retains only features quantified in at least 40% of the samples. This threshold can be adjusted if necessary.

pbf <- pb_filterNA(
    pbf, # without specifying, the filter will be applied to all assays in the object
    inplace = TRUE, # if false, filtered matrix will be saved as a new assay.
    pNA = 0.6 # the maximal proportion of missing values per feature
)
pbf
#> An instance of class ProBatchFeatures (type: bulk) with 2 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 331 rows and 118 columns 
#>   Processing chain:
#>   [1] add_level(protein)_byVar filterNA                 filterNA                
#>   Steps logged: 3 (see get_operation_log())

After filtering, the ProBatchFeatures object stores fewer features, which sharpens downstream diagnostics.

plot_NA_heatmap(
    pbf,
    pbf_name = c("peptide::raw", "protein::raw"),
    show_rownames = FALSE, show_row_dend = FALSE,
    color_by = "Lab", # currently only one level is supported
    border_color = NA # pheatmap parameter to remove cell borders
)

However, as the plots indicate, the remaining patterns of missingness are still associated with the lab.
To remove the residual systematic effects from the data, we further perform normalization and batch-effect correction.

4.2 Step 2 - log2-transformation

Proteomics data are commonly log-transformed to stabilize variance for statistical analysis and data visualization. Here, we perform a log2-transformation of both peptide- and protein-level data:

pbf <- log_transform_dm(
    pbf,
    log_base = 2, offset = 1,
    pbf_name = "protein::raw"
)

pbf <- log_transform_dm(
    pbf,
    log_base = 2, offset = 1,
    pbf_name = "peptide::raw"
)

pbf
#> An instance of class ProBatchFeatures (type: bulk) with 2 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 331 rows and 118 columns 
#>   Processing chain:
#>   [1] add_level(protein)_byVar filterNA                 filterNA                ; [4] log2                     log2                    
#>    - peptide: log2_on_raw
#>    - protein: log2_on_raw
#>   Steps logged: 5 (see get_operation_log())

After log-transformation, pbf records the virtual targets peptide::log2_on_raw and protein::log2_on_raw. Because log-transformation is fast to recompute, these results are replayed when requested instead of being stored as assays by default. Set store_fast_steps = TRUE or supply an explicit final_name to materialize them.

The intensity distributions of each sample after log-transformation can be visualized as boxplots to identify potential batch-specific shifts:

plot_boxplot(
    pbf,
    pbf_name = c(
        "peptide::log2_on_raw", "protein::log2_on_raw"
    ), # plot multiple on one
    sample_id_col = "Run",
    order_col = "Run",
    batch_col = "Lab",
    color_by_batch = TRUE,
    color_scheme = color_scheme,
    base_size = 7,
    plot_ncol = 1, # how many columns use for plotting multy-assay plot
    plot_title = c( # title for each assay
        "Peptide level - log2 scale",
        "Protein level - log2 scale"
    )
)

Additionally, the intensity distribution of peptides and proteins with and without NAs can be plotted:

plot_NA_density(
    pbf,
    pbf_name = c("peptide::log2_on_raw", "protein::log2_on_raw")
)

The comparison confirms that log2 stabilises the variance for both peptides and proteins, yet the cross-lab shifts in median intensity persist and will need to be addressed downstream.

Next, we will perform a Principal Component Analysis (PCA) to get a high-level overview of the sample clustering. We expect to see a strong batch effect, where samples cluster by their center of origin (“Lab”) rather than their biological condition (“Condition”).

pca1 <- plot_PCA(
    pbf,
    pbf_name = "protein::log2_on_raw",
    sample_id_col = "Run",
    color_scheme = color_scheme,
    color_by = "Lab",
    shape_by = "Condition",
    fill_the_missing = NULL, # remove incomplete feature rows for this plot
    plot_title = "NA rows removed, protein, log2",
    base_size = 10, point_size = 3, point_alpha = 0.5
)

pca2 <- plot_PCA(
    pbf,
    pbf_name = "protein::log2_on_raw",
    sample_id_col = "Run",
    color_scheme = color_scheme,
    color_by = "Condition",
    shape_by = "Lab",
    fill_the_missing = 0, # default value is -1
    plot_title = "NA replaced with 0, protein, log2",
    base_size = 10, point_size = 3, point_alpha = 0.5
)

library(gridExtra)
grid.arrange(pca1, pca2, ncol = 2, nrow = 1)

Similarly, the PCA plots clearly demonstrate that the primary source of variance in the raw data is the lab.

4.3 Step 3 - median normalization

To make the samples more comparable by removing the effect of sample load on intensities, we apply median normalization:

pbf <- pb_transform(
    object = pbf,
    from = "peptide::log2_on_raw",
    steps = "medianNorm"
)

pbf <- pb_transform(
    object = pbf,
    from = "protein::log2_on_raw",
    steps = "medianNorm"
)

pbf
#> An instance of class ProBatchFeatures (type: bulk) with 4 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 331 rows and 118 columns 
#>  [3] peptide::medianNorm_on_log2_on_raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [4] protein::medianNorm_on_log2_on_raw: SummarizedExperiment with 331 rows and 118 columns 
#>   Processing chain:
#>   [1] add_level(protein)_byVar filterNA                 filterNA                ; [4] log2                     log2                     medianNorm              ; [7] medianNorm              
#>    - peptide: log2_on_raw, medianNorm_on_log2_on_raw
#>    - protein: log2_on_raw, medianNorm_on_log2_on_raw
#>   Steps logged: 7 (see get_operation_log())

pb_transform() appends peptide::medianNorm_on_log2_on_raw and protein::medianNorm_on_log2_on_raw. We can now plot and compare these assays before and after normalization:

# boxplot
plot_boxplot(
    pbf,
    sample_id_col = "Run",
    order_col = "Run",
    batch_col = "Lab",
    color_by_batch = TRUE,
    color_scheme = color_scheme,
    base_size = 7,
    pbf_name = c("protein::log2_on_raw", "protein::medianNorm_on_log2_on_raw"),
    plot_ncol = 1,
    plot_title = c(
        "Before Median Normalization, protein level",
        "After Median Normalization, protein level"
    )
)

# plot PCA plot
plot_PCA(
    pbf,
    pbf_name = c("protein::log2_on_raw", "protein::medianNorm_on_log2_on_raw"),
    sample_id_col = "Run",
    color_scheme = color_scheme,
    color_by = "Lab",
    shape_by = "Condition",
    fill_the_missing = NULL, # remove incomplete feature rows for these plots
    plot_title = c(
        "Before Median Normalization, protein level",
        "After Median Normalization, protein level"
    ),
    base_size = 10, point_size = 3, point_alpha = 0.5
)

Although median normalization made intensity distributions more similar, the residual lab-specific batch effects remain evident on the PCA plots and need to be corrected.

4.4 Step 4 - batch-effect correction

The user can choose from two Core-registered batch-effect correction steps: combat, backed by ComBat (3) from sva, and limmaRBE, backed by removeBatchEffect() from limma (4). Both steps can be applied and logged with pb_transform().

Correction methods use the common missing-value policies "error", "keep", "drop_features", and "fill". The default, "error", stops when values are missing. "keep" passes missing values to a method that supports them, "drop_features" removes incomplete feature rows, and "fill" requires an explicit numeric fill_value. Here we choose "drop_features" so both correction backends receive complete feature rows.

# sample annotations used by both correction methods
sa <- colData(pbf) %>% as.data.frame()

Apply limma’s removeBatchEffect() to log-transformed and normalized data:

# limma::removeBatchEffect method
pbf <- pb_transform(
    pbf,
    from = "protein::log2_on_raw",
    steps = "limmaRBE",
    params_list = list(list(
        sample_annotation = sa,
        batch_col = "Lab",
        covariates_cols = c("Condition"),
        fill_the_missing = "drop_features"
    ))
)

pbf <- pb_transform(
    pbf,
    from = "protein::medianNorm_on_log2_on_raw",
    steps = "limmaRBE",
    params_list = list(list(
        sample_annotation = sa,
        batch_col = "Lab",
        covariates_cols = c("Condition"),
        fill_the_missing = "drop_features"
    ))
)

Because batch-effect correction is not a fast step, two new protein-level assays are materialized:

  • protein::limmaRBE_on_log2_on_raw
  • protein::limmaRBE_on_medianNorm_on_log2_on_raw

Both contain only complete feature rows because this example requested "drop_features".

print(pbf)
#> An instance of class ProBatchFeatures (type: bulk) with 6 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 331 rows and 118 columns 
#>  [3] peptide::medianNorm_on_log2_on_raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [4] protein::medianNorm_on_log2_on_raw: SummarizedExperiment with 331 rows and 118 columns 
#>  [5] protein::limmaRBE_on_log2_on_raw: SummarizedExperiment with 244 rows and 118 columns 
#>  [6] protein::limmaRBE_on_medianNorm_on_log2_on_raw: SummarizedExperiment with 244 rows and 118 columns 
#>   Processing chain:
#>   [1] add_level(protein)_byVar filterNA                 filterNA                ; [4] log2                     log2                     medianNorm              ; [7] medianNorm               limmaRBE                 limmaRBE                
#>    - peptide: log2_on_raw, medianNorm_on_log2_on_raw
#>    - protein: log2_on_raw, limmaRBE_on_log2_on_raw, medianNorm_on_log2_on_raw, limmaRBE_on_medianNorm_on_log2_on_raw
#>   Steps logged: 9 (see get_operation_log())

Alternatively, batch effect correction can be performed using ComBat:

# sva::ComBat wrapped via pb_transform()
pbf <- pb_transform(
    pbf,
    from = "protein::log2_on_raw",
    steps = "combat",
    params_list = list(list(
        sample_annotation = sa,
        batch_col = "Lab",
        sample_id_col = "Run",
        par.prior = TRUE,
        fill_the_missing = "drop_features"
    ))
)

pbf <- pb_transform(
    pbf,
    from = "protein::medianNorm_on_log2_on_raw",
    steps = "combat",
    params_list = list(list(
        sample_annotation = sa,
        batch_col = "Lab",
        sample_id_col = "Run",
        par.prior = TRUE,
        fill_the_missing = "drop_features"
    ))
)

The ComBat method from the sva package cannot be applied to data containing missing values. We therefore requested fill_the_missing = "drop_features" explicitly. To replace missing values with a constant instead, use fill_the_missing = "fill" together with a finite numeric fill_value; omitting a policy retains the safer "error" default.

print(pbf)
#> An instance of class ProBatchFeatures (type: bulk) with 8 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 331 rows and 118 columns 
#>  [3] peptide::medianNorm_on_log2_on_raw: SummarizedExperiment with 867 rows and 118 columns 
#>  ...
#>  [6] protein::limmaRBE_on_medianNorm_on_log2_on_raw: SummarizedExperiment with 244 rows and 118 columns 
#>  [7] protein::combat_on_log2_on_raw: SummarizedExperiment with 244 rows and 118 columns 
#>  [8] protein::combat_on_medianNorm_on_log2_on_raw: SummarizedExperiment with 244 rows and 118 columns 
#>   Processing chain:
#>    [1] add_level(protein)_byVar filterNA                 filterNA                ;  [4] log2                     log2                     medianNorm              ;  [7] medianNorm               limmaRBE                 limmaRBE                ; [10] combat                   combat                  
#>    - peptide: log2_on_raw, medianNorm_on_log2_on_raw
#>    - protein: log2_on_raw, combat_on_log2_on_raw, limmaRBE_on_log2_on_raw, medianNorm_on_log2_on_raw, combat_on_medianNorm_on_log2_on_raw, limmaRBE_on_medianNorm_on_log2_on_raw
#>   Steps logged: 11 (see get_operation_log())

Similarly, two new protein-level ComBat assays are materialized:

  • protein::combat_on_log2_on_raw
  • protein::combat_on_medianNorm_on_log2_on_raw

The corrected assays contain fewer protein groups than the filtered raw data whenever incomplete features were dropped.

All steps can be applied to peptide-level data in the same way:

pbf <- pb_transform(
    pbf,
    from = "peptide::medianNorm_on_log2_on_raw",
    steps = "limmaRBE",
    params_list = list(
        list(
            sample_annotation = sa,
            batch_col = "Lab",
            covariates_cols = c("Condition"),
            fill_the_missing = "drop_features"
        )
    )
)

pbf <- pb_transform(
    pbf,
    from = "peptide::medianNorm_on_log2_on_raw",
    steps = "combat",
    params_list = list(
        list(
            sample_annotation = sa,
            batch_col = "Lab",
            sample_id_col = "Run",
            par.prior = TRUE,
            fill_the_missing = "drop_features"
        )
    )
)

pbf
#> An instance of class ProBatchFeatures (type: bulk) with 10 sets:
#> 
#>  [1] peptide::raw: SummarizedExperiment with 867 rows and 118 columns 
#>  [2] protein::raw: SummarizedExperiment with 331 rows and 118 columns 
#>  [3] peptide::medianNorm_on_log2_on_raw: SummarizedExperiment with 867 rows and 118 columns 
#>  ...
#>  [8] protein::combat_on_medianNorm_on_log2_on_raw: SummarizedExperiment with 244 rows and 118 columns 
#>  [9] peptide::limmaRBE_on_medianNorm_on_log2_on_raw: SummarizedExperiment with 152 rows and 118 columns 
#>  [10] peptide::combat_on_medianNorm_on_log2_on_raw: SummarizedExperiment with 152 rows and 118 columns 
#>   Processing chain:
#>    [1] add_level(protein)_byVar filterNA                 filterNA                ;  [4] log2                     log2                     medianNorm              ;  [7] medianNorm               limmaRBE                 limmaRBE                ; [10] combat                   combat                   limmaRBE                ; [13] combat                  
#>    - peptide: log2_on_raw, medianNorm_on_log2_on_raw, combat_on_medianNorm_on_log2_on_raw, limmaRBE_on_medianNorm_on_log2_on_raw
#>    - protein: log2_on_raw, combat_on_log2_on_raw, limmaRBE_on_log2_on_raw, medianNorm_on_log2_on_raw, combat_on_medianNorm_on_log2_on_raw, limmaRBE_on_medianNorm_on_log2_on_raw
#>   Steps logged: 13 (see get_operation_log())

Two new peptide-level assays are materialized:

  • peptide::limmaRBE_on_medianNorm_on_log2_on_raw
  • peptide::combat_on_medianNorm_on_log2_on_raw

As at the protein level, each contains only feature rows that were complete across samples under the requested policy.

4.5 Step 5 - assess processed assays

The four new assays (protein::limmaRBE_on_… and protein::combat_on_…) are now stored alongside the previous processing steps. We can explore these assays using PCA to verify that the batch effect has been successfully removed and compare them to the data before batch effect correction.

4.5.1 Principal component analysis

The progression from log\(_2\)-transformed to batch-corrected assays demonstrates the expected reduction in lab-specific variance and a clearer separation by biological condition in the final panel.

plot_PCA(
    pbf,
    pbf_name = c(
        "protein::log2_on_raw",
        "protein::medianNorm_on_log2_on_raw",
        "protein::limmaRBE_on_medianNorm_on_log2_on_raw",
        "protein::combat_on_medianNorm_on_log2_on_raw"
    ),
    sample_id_col = "Run",
    color_scheme = color_scheme,
    color_by = "Lab",
    shape_by = "Condition",
    fill_the_missing = NULL,
    base_size = 8, point_size = 3, point_alpha = 0.5
)

Similarly, for peptide-level data:

plot_PCA(
    pbf,
    pbf_name = c(
        "peptide::log2_on_raw",
        "peptide::medianNorm_on_log2_on_raw",
        "peptide::limmaRBE_on_medianNorm_on_log2_on_raw",
        "peptide::combat_on_medianNorm_on_log2_on_raw"
    ),
    sample_id_col = "Run",
    color_scheme = color_scheme,
    color_by = "Lab",
    shape_by = "Condition",
    fill_the_missing = NULL,
    base_size = 8, point_size = 3, point_alpha = 0.5
)

4.5.2 t-SNE and UMAP

Core also provides plot_TSNE() and plot_UMAP() for matrices and ProBatchFeatures objects. Like PCA, these functions treat matrix columns as samples and align colData() by sample_id_col. Their plotting-family fill_the_missing control is deliberately different from the correction policy above: the default numeric value -1 fills missing values, while NULL (or FALSE) removes incomplete feature rows. Both choices emit a warning when missing values are present.

A single selected assay returns a static ggplot. Supply random_state explicitly when reproducible UMAP coordinates are required; omitting it leaves randomness under the caller’s control:

umap_protein <- plot_UMAP(
    pbf,
    pbf_name = "protein::combat_on_medianNorm_on_log2_on_raw",
    sample_id_col = "Run",
    color_by = "Lab",
    shape_by = "Condition",
    color_scheme = color_scheme,
    fill_the_missing = NULL,
    n_neighbors = 10,
    random_state = 2026,
    plot_title = "UMAP after protein-level ComBat correction"
)
umap_protein

For multiple assays, pbf_name = NULL selects assays in object order, whereas an explicit vector preserves the requested order. Shared arguments are supplied normally, and assay_args is a named list keyed by selected assay for settings that differ between assays. The next example requests a specific order, assigns per-assay t-SNE seeds and backend settings, and asks for both the arranged grob and its named plot list:

embedding_assays <- c(
    "protein::limmaRBE_on_medianNorm_on_log2_on_raw",
    "protein::combat_on_medianNorm_on_log2_on_raw"
)
n_embedding_samples <- ncol(pb_assay_matrix(pbf, embedding_assays[[1L]]))
tsne_perplexity <- min(
    10,
    max(1, floor((n_embedding_samples - 1) / 3))
)
tsne_args <- setNames(
    list(
        list(
            perplexity = tsne_perplexity,
            max_iter = 300,
            random_seed = 2026
        ),
        list(
            perplexity = tsne_perplexity,
            max_iter = 300,
            random_seed = 2027
        )
    ),
    embedding_assays
)

tsne_protein <- plot_TSNE(
    pbf,
    pbf_name = embedding_assays,
    sample_id_col = "Run",
    color_by = "Lab",
    shape_by = "Condition",
    color_scheme = color_scheme,
    fill_the_missing = NULL,
    plot_title = setNames(
        c("limma corrected", "ComBat corrected"),
        embedding_assays
    ),
    assay_args = tsne_args,
    return_gridExtra = TRUE,
    plot_ncol = 2
)
names(tsne_protein$plots) # same order as embedding_assays
#> [1] "protein::limmaRBE_on_medianNorm_on_log2_on_raw"
#> [2] "protein::combat_on_medianNorm_on_log2_on_raw"
grid::grid.draw(tsne_protein$grob)

random_seed controls t-SNE reproducibility. If a requested perplexity is too large for an assay’s sample count, Core reduces it to the supported maximum and warns.

Static rendering is the default and does not require Plotly. To make the embeddings interactive, set use_plotlyrender = TRUE. One assay returns one plotly object, while multiple assays return a named ordered list. Set return_subplots = TRUE to combine a multi-assay result; subplot_ncol and share_axes control its layout. The example below shows both return forms:

embedding_assays <- c(
    "protein::limmaRBE_on_medianNorm_on_log2_on_raw",
    "protein::combat_on_medianNorm_on_log2_on_raw"
)
umap_args <- setNames(
    list(
        list(n_neighbors = 10, random_state = 2026),
        list(n_neighbors = 15, random_state = 2027)
    ),
    embedding_assays
)

umap_plotly_list <- plot_UMAP(
    pbf,
    pbf_name = embedding_assays,
    sample_id_col = "Run",
    color_by = "Lab",
    shape_by = "Condition",
    color_scheme = color_scheme,
    fill_the_missing = NULL,
    assay_args = umap_args,
    use_plotlyrender = TRUE
)
names(umap_plotly_list) # same order as embedding_assays

umap_subplot <- plot_UMAP(
    pbf,
    pbf_name = embedding_assays,
    sample_id_col = "Run",
    color_by = "Lab",
    shape_by = "Condition",
    color_scheme = color_scheme,
    fill_the_missing = NULL,
    assay_args = umap_args,
    use_plotlyrender = TRUE,
    return_subplots = TRUE,
    subplot_ncol = 2,
    share_axes = c(x = TRUE, y = TRUE)
)
umap_subplot

4.5.3 Hierarchical clustering

Hierarchical clustering provides another way to explore sample similarities. Here we display clustering results for protein::combat_on_medianNorm_on_log2_on_raw and protein::limmaRBE_on_medianNorm_on_log2_on_raw assays and colour samples by lab and condition:

plot_hierarchical_clustering(
    pbf,
    pbf_name = "protein::combat_on_medianNorm_on_log2_on_raw",
    sample_id_col = "Run",
    label_font = 0.6,
    color_list = color_scheme
)

plot_hierarchical_clustering(
    pbf,
    pbf_name = "protein::limmaRBE_on_medianNorm_on_log2_on_raw",
    sample_id_col = "Run",
    label_font = 0.6,
    color_list = color_scheme,
    fill_the_missing = NULL
)

4.5.4 Principal Variance Component Analysis (PVCA)

PVCA quantifies the contribution of known factors to the observed variance (5) which can be compared across raw, normalised, and ComBat-corrected assays:

plot_PVCA(
    pbf,
    pbf_name = c(
        "protein::medianNorm_on_log2_on_raw",
        "protein::combat_on_medianNorm_on_log2_on_raw"
    ),
    sample_id_col = "Run",
    technical_factors = c("Lab"),
    biological_factors = c("Condition"),
    fill_the_missing = NULL,
    base_size = 7,
    plot_ncol = 2, # the number of plots in a row
    variance_threshold = 0 # the percentile value of weight each of the
    # covariates needs to explain
    # (the rest will be lumped together)
)

5 Inspecting operation log

Each call to pb_transform() or pb_filterNA() creates a log entry that records the source assay, the applied operation, and the resulting assay name. Reviewing this log ensures that complex processing pipelines remain transparent and reproducible.

get_operation_log(pbf) %>%
    as_tibble()
#> # A tibble: 13 × 7
#>    step                 fun   from  to    params       timestamp           pkg  
#>    <chr>                <chr> <chr> <chr> <list>       <dttm>              <chr>
#>  1 add_level(protein)_… addA… pept… prot… <named list> 2026-08-02 23:42:11 proB…
#>  2 filterNA             filt… pept… pept… <named list> 2026-08-02 23:42:16 proB…
#>  3 filterNA             filt… prot… prot… <named list> 2026-08-02 23:42:16 proB…
#>  4 log2                 log2  prot… prot… <named list> 2026-08-02 23:42:17 proB…
#>  5 log2                 log2  pept… pept… <named list> 2026-08-02 23:42:17 proB…
#>  6 medianNorm           medi… pept… pept… <list [0]>   2026-08-02 23:42:26 proB…
#>  7 medianNorm           medi… prot… prot… <list [0]>   2026-08-02 23:42:26 proB…
#>  8 limmaRBE             limm… prot… prot… <named list> 2026-08-02 23:42:32 proB…
#>  9 limmaRBE             limm… prot… prot… <named list> 2026-08-02 23:42:33 proB…
#> 10 combat               comb… prot… prot… <named list> 2026-08-02 23:42:33 proB…
#> 11 combat               comb… prot… prot… <named list> 2026-08-02 23:42:33 proB…
#> 12 limmaRBE             limm… pept… pept… <named list> 2026-08-02 23:42:33 proB…
#> 13 combat               comb… pept… pept… <named list> 2026-08-02 23:42:33 proB…

If you need a compact name for the current assay chain (for example, to label benchmark results), use pb_pipeline_name(pbf, "protein::combat_on_medianNorm_on_log2_on_raw").

pb_pipeline_name(pbf, "protein::combat_on_medianNorm_on_log2_on_raw")
#> [1] "combat_on_medianNorm_on_log2_on_filterNA_on_add_level(protein)_byVar_on_filterNA_on_raw"

5.1 Extracting data matrices from pbf object

It is possible to extract a specific data matrix from a ProBatchFeatures object using the assay name:

extracted_matrix <- pb_assay_matrix(
    pbf,
    assay = "peptide::raw"
)
# show
extracted_matrix[1:5, 1:5]
#>                 lab_A_S1 lab_A_S2 lab_A_S3 lab_A_S4 lab_A_S5
#> ELGNWKDFIEVMLR3       NA       NA       NA       NA       NA
#> DFIEVMLR2             NA       NA       NA       NA       NA
#> AAADLISR2       14647900 13463600  8518130 12172600 11612300
#> AGIGINAGR2      11785300 11025600  9713860  9612410 13795300
#> AEQYLLENETTK2   30835300 29444100 19451100 20090000 22758800
rm(extracted_matrix)

For consistency with the main ProBatch vignette, assay can also be extracted in the long format. In this case, the sample_id_col needs to be specified:

extracted_long <- pb_as_long(
    pbf,
    sample_id_col = "Run",
    pbf_name = "peptide::raw"
)

# show
extracted_long[1:3, ]
#>     feature_label      Run Intensity   Lab Condition
#> 1 ELGNWKDFIEVMLR3 lab_A_S1        NA lab_A       Pyr
#> 2 ELGNWKDFIEVMLR3 lab_A_S2        NA lab_A       Pyr
#> 3 ELGNWKDFIEVMLR3 lab_A_S3        NA lab_A       Glu
rm(extracted_long)

6 Session info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 LTS
#> 
#> Matrix products: default
#> BLAS:   /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so 
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
#>  [3] LC_TIME=en_GB              LC_COLLATE=C              
#>  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
#>  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
#>  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
#> 
#> time zone: America/New_York
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats4    stats     graphics  grDevices utils     datasets  methods  
#> [8] base     
#> 
#> other attached packages:
#>  [1] QFeatures_1.23.1            MultiAssayExperiment_1.39.0
#>  [3] SummarizedExperiment_1.43.0 Biobase_2.73.2             
#>  [5] GenomicRanges_1.65.1        Seqinfo_1.3.0              
#>  [7] IRanges_2.47.2              S4Vectors_0.51.6           
#>  [9] BiocGenerics_0.59.10        generics_0.1.4             
#> [11] MatrixGenerics_1.25.0       matrixStats_1.5.0          
#> [13] gridExtra_2.3.1             knitr_1.51                 
#> [15] proBatch_2.1.3              ggplot2_4.0.3              
#> [17] tibble_3.3.1                dplyr_1.2.1                
#> [19] BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>   [1] splines_4.6.1           ggplotify_0.1.3         preprocessCore_1.75.0  
#>   [4] XML_3.99-0.23           rpart_4.1.27            lifecycle_1.0.5        
#>   [7] Rdpack_2.6.6            fastcluster_1.3.0       edgeR_4.11.4           
#>  [10] doParallel_1.0.17       lattice_0.22-9          MASS_7.3-66            
#>  [13] backports_1.5.1         magrittr_2.0.5          limma_3.69.2           
#>  [16] Hmisc_5.2-6             sass_0.4.10             rmarkdown_2.31         
#>  [19] jquerylib_0.1.4         yaml_2.3.12             wesanderson_0.3.7      
#>  [22] otel_0.2.0              askpass_1.2.1           reticulate_1.46.0      
#>  [25] MsCoreUtils_1.25.4      DBI_1.3.0               minqa_1.2.8            
#>  [28] RColorBrewer_1.1-3      lubridate_1.9.5         abind_1.4-8            
#>  [31] Rtsne_0.17              purrr_1.2.2             AnnotationFilter_1.37.0
#>  [34] yulab.utils_0.2.4       nnet_7.3-20             rappdirs_0.3.4         
#>  [37] sva_3.61.0              genefilter_1.95.0       pheatmap_1.0.13        
#>  [40] umap_0.2.10.0           RSpectra_0.16-2         annotate_1.91.0        
#>  [43] codetools_0.2-20        DelayedArray_0.39.3     tidyselect_1.2.1       
#>  [46] farver_2.1.2            lme4_2.0-6              viridis_0.6.5          
#>  [49] dynamicTreeCut_1.63-1   base64enc_0.1-6         jsonlite_2.0.0         
#>  [52] Formula_1.2-5           survival_3.8-9          iterators_1.0.14       
#>  [55] foreach_1.5.2           tools_4.6.1             Rcpp_1.1.2             
#>  [58] glue_1.8.1              SparseArray_1.13.2      BiocBaseUtils_1.15.1   
#>  [61] xfun_0.60               mgcv_1.9-4              ggfortify_0.4.19       
#>  [64] HDF5Array_1.41.0        withr_3.0.3             BiocManager_1.30.27    
#>  [67] fastmap_1.2.0           pvca_1.53.0             rhdf5filters_1.25.3    
#>  [70] boot_1.3-32             openssl_2.4.2           digest_0.6.39          
#>  [73] gridGraphics_0.5-1      timechange_0.4.0        R6_2.6.1               
#>  [76] colorspace_2.1-3        dichromat_2.0-1         RSQLite_3.53.3         
#>  [79] h5mread_1.5.0           utf8_1.2.6              tidyr_1.3.2            
#>  [82] data.table_1.18.4       httr_1.4.8              htmlwidgets_1.6.4      
#>  [85] S4Arrays_1.13.0         pkgconfig_2.0.3         gtable_0.3.6           
#>  [88] blob_1.3.0              S7_0.2.2                impute_1.87.0          
#>  [91] XVector_0.53.0          htmltools_0.5.9         bookdown_0.47          
#>  [94] ProtGenerics_1.45.0     clue_0.3-68             scales_1.4.0           
#>  [97] png_0.1-9               reformulas_0.4.4        corrplot_0.95          
#> [100] rstudioapi_0.19.0       reshape2_1.4.5          checkmate_2.3.4        
#> [103] nlme_3.1-170            nloptr_2.2.1            cachem_1.1.0           
#> [106] rhdf5_2.57.3            stringr_1.6.0           parallel_4.6.1         
#> [109] foreign_0.8-91          AnnotationDbi_1.75.2    vsn_3.81.0             
#> [112] pillar_1.11.1           grid_4.6.1              vctrs_0.7.3            
#> [115] xtable_1.8-8            cluster_2.1.8.3         htmlTable_2.5.0        
#> [118] evaluate_1.0.5          tinytex_0.60            magick_2.9.1           
#> [121] cli_3.6.6               locfit_1.5-9.12         compiler_4.6.1         
#> [124] rlang_1.3.0             crayon_1.5.3            labeling_0.4.3         
#> [127] fs_2.1.0                affy_1.91.0             plyr_1.8.9             
#> [130] stringi_1.8.7           viridisLite_0.4.3       WGCNA_1.74             
#> [133] BiocParallel_1.47.0     Biostrings_2.81.6       lazyeval_0.2.3         
#> [136] Matrix_1.7-6            bit64_4.8.2             Rhdf5lib_2.1.0         
#> [139] KEGGREST_1.53.6         statmod_1.5.2           rbibutils_2.4.1        
#> [142] igraph_2.3.3            memoise_2.0.1           affyio_1.83.0          
#> [145] bslib_0.11.0            bit_4.6.0

7 Citation

To cite this package, please use:

citation("proBatch")
#> To cite proBatch in publications use:
#> 
#>   Cuklina J, Lee CH, Williams EG, Sajic T, Collins BC, Rodriguez
#>   Martinez M, Sharma VS, Wendt F, Goetze S, Keele GR, Wollscheid B,
#>   Aebersold R, Pedrioli PGA. Diagnostics and correction of batch
#>   effects in large-scale proteomic studies: a tutorial. Molecular
#>   Systems Biology 17(8), e10240 (2021).
#>   https://doi.org/10.15252/msb.202110240
#> 
#> A BibTeX entry for LaTeX users is
#> 
#>   @Article{,
#>     title = {Diagnostics and correction of batch effects in large-scale proteomic studies: a tutorial},
#>     author = {Jelena Cuklina and Chloe H. Lee and Evan G. Williams and Tatjana Sajic and Ben C. Collins and Maria Rodriguez-Martinez and Varun S. Sharma and Fabian Wendt and Sandra Goetze and Gregory R. Keele and Bernd Wollscheid and Ruedi Aebersold and Patrick G. A. Pedrioli},
#>     journal = {Molecular Systems Biology},
#>     year = {2021},
#>     volume = {17},
#>     number = {8},
#>     pages = {e10240},
#>     doi = {10.15252/msb.202110240},
#>     url = {https://doi.org/10.15252/msb.202110240},
#>   }

References

1. Gatto, L. & Vanderaa, C. QFeatures: Quantitative features for mass spectrometry data. (2025). at <https://github.com/RforMassSpectrometry/QFeatures>

2. Burankova, Y. et al. Privacy-preserving multicenter differential protein abundance analysis with fedprot. Nature Computational Science 1–14 (2025).

3. Johnson, W. E., Li, C. & Rabinovic, A. Adjusting batch effects in microarray expression data using empirical bayes methods. Biostatistics 8, 118–127 (2006).

4. Ritchie, M. E. et al. Limma powers differential expression analyses for rna-sequencing and microarray studies. Nucleic Acids Research 43, e47–e47 (2015).

5. Bushel, P. pvca: Principal variance component analysis (pvca). (2025). doi:10.18129/B9.bioc.pvca