Contents

1 Overview

normScore provides tools to evaluate and rank normalization methods for omics data using a composite multi-metric scoring framework.

The package is designed to help compare different normalization strategies based on several complementary criteria, including variability, correlation structure, MA-trend behavior, relative log expression consistency, and sample-level intensity consistency.

In this vignette:

  1. We explain the main basics of normScore metric.
  2. We show a minimal workflow using simulated data.
  3. We provide an example of integration of normScore with SummarizedExperiment and NormalyzerDE R Packages.

2 Methods

2.1 General rationale

A good normalization method should ideally:

  • reduce unwanted technical variability
  • preserve biological structure
  • avoid introducing systematic distortions
  • maintain sample comparability
  • stabilize distributions across samples

Because these goals are not fully captured by a single number, normScore uses a multi-item scoring approach.

Lower final scores indicate better normalization performance relative to the other methods under comparison.

2.2 Inputs required by normScore

The main function of the package expects three key inputs:

  1. a named list of normalized datasets
  2. a sample annotation table with group membership
  3. the raw data matrix

Each normalized dataset should contain the same features and samples, arranged in the same order.

The raw data are used to derive the log-transformed baseline normalization and to compute a correction factor related to global sample intensity variation.

2.3 The six scoring items

The composite score is based on six item-wise metrics.

2.3.1 Item 1. Pooled coefficient of variation (PCV)

This item summarizes within-group variability across samples.

For each group, the coefficient of variation is computed feature-wise, and the mean value is summarized across groups. Lower values indicate lower relative dispersion within biological groups.

This item captures whether a normalization method improves sample consistency without reference to differential structure.

2.3.2 Item 2. Within-group sample correlation

This item evaluates pairwise sample correlations within each group.

For each group, all unique pairwise sample-sample correlations are extracted, and a summary statistic based on the median and interquartile range is used. The score is transformed so that lower values indicate better performance.

This item reflects whether samples from the same group become more similar after normalization.

2.3.3 Item 3. MA-trend deviation

This item assesses whether the relationship between fold change and average expression shows systematic distortion.

For two selected groups, an MA-like summary is constructed using:

  • logFC: the difference between group means
  • AveExpr: the average of group means

A linear model is fitted, and the area between the fitted line and the expected horizontal line logFC = 0 is computed. An additional shape correction factor penalizes undesirable trends in the spread of logFC across the expression range.

Lower values indicate a flatter and more stable MA trend.

2.3.4 Item 4. Mean-SD trend deviation

This item evaluates whether sample standard deviations depend systematically on sample mean intensity.

For each sample, the mean and standard deviation are computed. Samples are ordered by mean intensity, and a linear model is fitted using standard deviation as a function of sample order.

The slope of this fitted trend is then converted into an area-based metric measuring deviation from horizontality.

Lower values indicate a weaker dependency between sample mean and sample variability.

2.3.5 Item 5. RLE consistency

This item is inspired by the relative log expression (RLE) concept.

For each feature, the median across samples is used as a reference. Relative expression values are obtained by dividing each observation by the feature-wise median, and the median relative expression is then computed sample-wise.

The final score is the mean absolute percentage error of these sample-wise medians relative to 1.

Lower values indicate that sample-wise medians are closer to the expected reference value.

2.3.6 Item 6. Total intensity consistency

This item compares the distributional center and spread of sample intensity profiles.

For each sample, the first quartile, median, and third quartile are computed. The score is based on the sum of the mean absolute percentage errors of these three summaries relative to their global sample medians.

Lower values indicate more consistent sample distributions.

2.4 Scaling and aggregation

Because the six items are on different scales, each one is transformed using min-max scaling before aggregation.

After scaling, the six items are summed into a total score.

This makes the final ranking easier to interpret and ensures that no single item dominates only because of its numerical range.

2.5 Correction applied to the log-transformed baseline

The log-transformed raw data are included in the normalization list as a baseline reference.

A correction factor derived from the coefficient of variation of total raw sample intensities is applied specifically to the "Log" method. This is used to account for the extent of global intensity imbalance already present in the raw data.

2.6 Additional adjustments and weighting

During the development of the scoring framework, two additional adjustments were introduced to improve the balance between items and avoid dominance by specific components.

First, the contribution of the within-group correlation metric (Item 2) was reduced. This item was observed to have relatively low discriminative power across normalization methods, often yielding very similar values. To prevent it from contributing disproportionately to the final score despite its limited ability to differentiate methods, its scaled value was down-weighted by a factor of 0.1.

Second, the correction factor applied to the log-transformed baseline normalization (Item 0) was moderated. This factor is derived from the coefficient of variation of total sample intensities in the raw data and is intended to account for global intensity imbalance.

However, without adjustment, this correction could have an excessive influence on the final score, particularly for datasets with very low or very high values of the correction factor. To mitigate this effect, the magnitude of the correction was reduced by scaling it down (e.g., dividing or adjusting its impact), ensuring that it contributes to the ranking without dominating it.

These adjustments were introduced to maintain a balanced contribution across items and to improve the stability and interpretability of the final ranking.

3 Basic usage

3.1 Simulating example data

To keep the package lightweight and reproducible, normScore includes a data simulation function that generates proteomics-like two-group datasets.

# BiocManager::install("normScore")
library(normScore)

simData <- simulateData(
    nProteins = 500,
    nPerGroup = 4,
    sampleShiftSd = 0.12,
    sampleShiftCap = 1,
    sampleSdStrength = 1.2,
    sampleSdRho = 0.3,
    rhoWithin = 0.6,
    rhoBetween = 0.3,
    loadingSd = 0.45,
    sigmaHi = 0.03,
    sigmaLo = 0.4,
    gammaSigma = 3.8,
    kMnar = 1.6,
    addMissing = FALSE,
    seed = 123
)

The simulated object contains three main elements:

names(simData)
#> [1] "logData"  "rawData"  "metadata"
  • logData: simulated data on the log2 scale
  • rawData: simulated data on the raw scale
  • metadata: sample annotation with group membership

We can inspect the sample metadata:

head(simData$metadata)
#>   Samples Groups
#> 1    G1_1     G1
#> 2    G1_2     G1
#> 3    G1_3     G1
#> 4    G1_4     G1
#> 5    G2_1     G2
#> 6    G2_2     G2

3.2 Building a list of normalized datasets

The main function of the package, normScore::normScore(), expects a named list of normalized datasets. Each element of the list should be a matrix or data frame with the same rows (features) and columns (samples).

For illustration purposes, we create a few simple variants of the simulated log-scale data. In a real analysis, these would correspond to different normalization methods.

# Median normalization
MedianMat <- NormalyzerDE::medianNormalization(simData$rawData)
colnames(MedianMat) <- colnames(simData$rawData)
rownames(MedianMat) <- rownames(simData$rawData)

# Quantile normalization
QuantileMat <- limma::normalizeQuantiles(
    simData$logData
)
colnames(QuantileMat) <- colnames(simData$logData)
rownames(QuantileMat) <- rownames(simData$logData)


# CyclicLoess normalization
CyclicNorm <- limma::normalizeCyclicLoess(
    simData$logData, method="fast",  adaptive.span=FALSE
)


# All together
normalizedDataList <- list(
    Log = simData$logData,
    Median = MedianMat,
    Quantile = QuantileMat, 
    CyclicLoess = CyclicNorm
)

names(normalizedDataList)
#> [1] "Log"         "Median"      "Quantile"    "CyclicLoess"

3.3 Running normScore()

We now compute the normalization ranking.

result <- normScore(
    normalizedDataList = normalizedDataList,
    groupData = simData$metadata,
    rawData = simData$rawData,
    returnDetails = TRUE,
    nBoot = 150
)
#> Reference group set to G2.
#> Alternative group set to G1.

Important considerations! Before running normScore, verify that:

  • all assays contain the same rows and columns;
  • feature and sample names are unique;
  • the order of groupData$Samples matches the assay columns;
  • missing values are represented consistently across assays;
  • every sample has a valid group assignment.

3.4 Final ranking

The finalRanking element contains the final composite score for each normalization method.

result$finalRanking
#>    Quantile CyclicLoess      Median         Log 
#>    0.000000    1.221309    3.389992    6.202657

Lower values indicate better normalization performance according to the scoring framework implemented in the package.

3.5 Detailed scores

The detailRanking element contains the scaled item-wise scores used to build the final ranking.

result$detailRanking
#>                 Item1 Item2      Item3     Item4     Item5     Item6    Total
#> Quantile    0.0000000   0.0 0.00000000 0.0000000 0.0000000 0.0000000 0.000000
#> CyclicLoess 0.2043477   0.1 0.09279472 0.2479291 0.3949613 0.1812763 1.221309
#> Median      0.9599428   0.0 0.95657880 0.6217586 0.1067698 0.7449419 3.389992
#> Log         1.0000000   0.0 1.00000000 1.0000000 1.0000000 1.0000000 5.000000
#>             TotalxItem0
#> Quantile       0.000000
#> CyclicLoess    1.221309
#> Median         3.389992
#> Log            6.202657

The six item scores correspond to:

  1. pooled coefficient of variation
  2. within-group sample correlation
  3. MA-plot trend deviation
  4. mean-SD trend deviation
  5. RLE consistency
  6. total intensity consistency

These are scaled and combined into a final score.

3.6 Bootstrap summary

If returnDetails = TRUE, the function also computes bootstrap-based summary scores and confidence intervals.

result$bootstrapScore
#>   normalization meanNormScore     ll95      ul95
#> 1      Quantile      0.000000 0.000000  0.000000
#> 2   CyclicLoess      3.524893 1.456764  6.607371
#> 3        Median      9.802567 4.350631 18.035879
#> 4           Log     19.443317 8.683720 35.546264

This output contains:

  • the normalization method name
  • the mean bootstrap normScore
  • the lower 95% confidence bound
  • the upper 95% confidence bound

3.7 Bootstrap plot

The bootstrap summary can also be visualized directly:

plotBootstrapNormScore(result)

4 Diagnostic plots

Individual diagnostic plots from state-of-art normalization criteria can be visualized with the following function:

allPlots <- plotNormScoreDiagnostics(
    normalizedDataList = normalizedDataList,
    groupData = simData$metadata,
    rawData = simData$rawData
)
#> Reference group set to G2.
#> Alternative group set to G1.

names(allPlots)
#> [1] "I0_SystematicBias" "I1_PCV"            "I2_WithinGroupCor"
#> [4] "I3_MAplot"         "I4_MeanSD"         "I5_RLE"           
#> [7] "I6_TotalIntensity"

4.1 Item 0 - Systematic Bias Penalty

In the following barplot, each bar represents the mean of total intensity or protein quantity for each sample from the raw data, that is, without any transformations or normalizations applied. The magnitude of the systematic bias in the current data can be assessed using this figure.

If the mean of total intensity across samples shows large differences (for example, a three-fold difference between the highest and the lowest values), then the variability in the data is significant. In this case, a more robust normalization method than just the logarithm transformation is likely to be needed.

allPlots[["I0_SystematicBiasPenalty"]]
#> NULL

4.2 Item 1. Pooled Coefficient of Variance (PCV)

The coefficient of variation (CV) indicates the dispersion of all data around the average of each variable, i.e., each protein. Therefore, it is interpreted as a percentage. High CV values indicate widely scattered data. In this app, the CV is estimated for each group and referred to as the pooled coefficient of variation (PCV). A lower PCV is expected for normalized data against only log-transformed data, which means that normalization reduces the variability within each group.

Once PCV is calculated, the results are displayed in two formats: If the number of groups is fewer than 5, the figure will show the point estimate of the mean and its 95% confidence interval. However, with 5 or more groups, a box plot will be displayed to depict the distribution of PCV across the different groups.

allPlots[["I1_PCV"]]

4.3 Item 2. Within-group sample correlation

We expect a strong correlation of protein groups intensities between samples within the same group, as they reflect the same condition. In this section, coefficients of correlation are estimated by pairs of samples within each study group, using a non-parametric approach: Spearman’s correlation. The distribution of these correlations is then plotted for each normalization.

The correlations are also estimated by normalization type. The best normalization method will be the one that increases intragroup correlations the most, indicating a reduction in intragroup variation.

allPlots[["I2_WithinGroupCorrelation"]]
#> NULL

4.4 Item 3. MA-trend deviation

As most of the proteins will not be differentially expressed between groups, a logFC close to 0 is expected, so most of the points should fall along this line on Y-axis. Moreover, it is more difficult to properly estimate proteins at lower concentrations, and as a consequence, higher dispersion is expected when the average of protein quantities in both groups are small (X-axis). In contrast, proteins with higher concentrations in both groups will be accurately estimated, so less dispersion is expected at the ends of the X-axis.

allPlots[["I3_MAplot"]]

4.5 Item 4. Mean-SD trend deviation

In this type of plot, samples are sorted by their average quantity, and then the standard deviation is represented versus the previously established average order.

How should this graph be interpreted?

There should be no correlation between the standard deviation and the mean of protein quantities through the different samples. For this reason, the best normalization would be the one in which the points are randomly distributed on the plot.

allPlots[["I4_MeanSD"]]

4.6 Item 5. RLE consistency

For each protein, the median of intensity across all samples is estimated, and then the intensity of each protein is divided by its median. This allows us to obtain the relative expression of a protein from a specific sample relative to all samples. Finally, this data is log-transformed, and its distribution (by sample) is displayed using a boxplot.

How should this graph be interpreted?

Most of the proteins should show a log relative expression value close to 0, as the majority of the protein quantification values should match their median across all samples (ratio = 1, and log(1) = 0). For this reason, the boxplots across different samples should be centered around 0, and the best normalization would be the one that achieves this outcome.

allPlots[["I5_RLE"]]

4.7 Item 6. Total intensity consistency

Using all protein quantities in each sample, their distribution is represented through box plots once all normalizations have been applied to the data.

How should this graph be interpreted?

Effective removal of systematic error will result in homogeneous distributions of protein quantities across all samples.

allPlots[["I6_TotalIntensity"]]

5 Interoperability with other packages

SummarizedExperiment is a standard Bioconductor container for storing high-dimensional molecular data together with feature and sample metadata. A single object can contain several assays, making it convenient for keeping the original data and multiple normalized versions aligned. In addition, NormalizerDE is a gold standard option to assess normalization methods.

Due to the relevance of both packages in relation to normScore goals, we show in this last section how to use them together.

# BiocManager::install("normScore")
library(normScore)

# BiocManager::install("SummarizedExperiment")
library(SummarizedExperiment)

# BiocManager::install("S4Vectors")
library(S4Vectors)

# BiocManager::install("NormalyzerDE")
library(NormalyzerDE)

5.1 Example data

For illustration, we generate a small log-transformed abundance matrix using simulateData function from normScore data package.

Rows represent proteins and columns represent samples in logData and rawData items from the output simData object.

simData <- simulateData(
    nProteins = 500,
    nPerGroup = 4,
    sampleShiftSd = 0.12,
    sampleShiftCap = 1,
    sampleSdStrength = 1.2,
    sampleSdRho = 0.3,
    rhoWithin = 0.6,
    rhoBetween = 0.3,
    loadingSd = 0.45,
    sigmaHi = 0.03,
    sigmaLo = 0.4,
    gammaSigma = 3.8,
    kMnar = 1.6,
    addMissing = FALSE,
    seed = 456
)

rawData <- simData$rawData

5.2 Build a SummarizedExperiment object

The original and normalized matrices are stored as separate assays. Sample-level information is stored in colData, whereas optional feature-level annotations can be stored in rowData.

sampleData <- S4Vectors::DataFrame(
    sample = simData$metadata$Samples,
    group = simData$metadata$Groups
)

se <- SummarizedExperiment::SummarizedExperiment(
    assays = list(rawData = rawData),
    colData = sampleData,
    metadata = list(
        sample = "sample",
        group = "group")
)

5.3 Normalization using NormalyzerDE

Using normalyzerDE different normalization methods are simultaneously applied to the input data.

normalyzerObject <- NormalyzerDE::getVerifiedNormalyzerObject(
    jobName = "normScore_example",
    summarizedExp = se,
    noLogTransform = FALSE
)
#> Input data checked. All fields are valid.
#> Sample check: More than one sample group found
#> Sample replication check: All samples have replicates
#> No RT column found, skipping RT processing

normalyzerResults <- NormalyzerDE::normMethods(
    normalyzerObject,
    normalizeRetentionTime = FALSE
)
#> No RT column specified (column named 'RT') or option not specified Skipping RT normalization.

5.4 Prepare normScore inputs

normScore evaluates a named list of matrices. We therefore extract the normalized assays from the SummarizedExperiment object.

normalizedDataList <- methods::slot(
    normalyzerResults,
    "normalizations"
)

Setting rownames as some normalized matrices from normalyzerDE lost protein rownames:

normalizedDataList <- lapply(
    normalizedDataList,
    function(x) {
        rownames(x) <- rownames(rawData)
        x
    }
)

The group information is constructed from colData. The sample order must match the matrix columns.

groupData <- data.frame(
    Samples = sampleData$sample,
    Groups = sampleData$group
)

stopifnot(
    identical(groupData$Samples, colnames(rawData)),
    all(vapply(
        normalizedDataList,
        function(x) identical(colnames(x), groupData$Samples),
        logical(1)
    ))
)

groupData
#>   Samples Groups
#> 1    G1_1     G1
#> 2    G1_2     G1
#> 3    G1_3     G1
#> 4    G1_4     G1
#> 5    G2_1     G2
#> 6    G2_2     G2
#> 7    G2_3     G2
#> 8    G2_4     G2

5.5 Run normScore

The extracted objects can now be supplied to normScore.

normScoreResults <- normScore(
    normalizedDataList = normalizedDataList,
    groupData = groupData,
    rawData = rawData
)
#> Reference group set to G2.
#> Alternative group set to G1.

Depending on the final argument names of the exported normScore::normScore() function, this call may need to be adapted to the current package API.

The returned object can be inspected directly:

names(normScoreResults)
#> [1] "finalRanking"   "detailRanking"  "bootstrapScore"

normScoreResults$finalRanking
#>   Quantile        RLR   CycLoess     median        VSN        Log         GI 
#> 0.03903797 0.31206666 0.48460383 2.62167567 2.97196236 3.26402195 4.96302977 
#>       mean 
#> 4.96484976

5.6 Store results in the SummarizedExperiment metadata

Analysis-level information that is not naturally associated with individual features or samples can be stored in metadata(se).

metadata(se)$normScore <- normScoreResults

The results can then be recovered without separating them from the data:

storedResults <- metadata(se)$normScore
names(storedResults)
#> [1] "finalRanking"   "detailRanking"  "bootstrapScore"

A compact summary table may also be stored separately:

metadata(se)$normalizationRanking <- normScoreResults$score

5.7 Add a selected normalized assay

Once a normalization method has been selected, its matrix can remain as an assay of the same object. For example, if median normalization is selected:

assay(se, "Selected") <- normalizedDataList$Quantile

assayNames(se)
#> [1] "rawData"  "Selected"

The selected assay can then be used in downstream Bioconductor workflows:

selectedData <- assay(se, "Selected")

dim(selectedData)
#> [1] 500   8

5.8 Working with an existing SummarizedExperiment

For an existing object, the minimum workflow is:

assayNames(se)
colnames(se)
colData(se)

rawData <- assay(se, "Unnormalized")

normalizationMethods <- c("Mean", "Median", "Quantile")

normalizedDataList <- lapply(
    normalizationMethods,
    function(method) assay(se, method)
)
names(normalizedDataList) <- normalizationMethods

groupData <- data.frame(
    Samples = colnames(se),
    Groups = as.character(colData(se)$group)
)

normScoreResults <- normScore(
    normalizedDataList = normalizedDataList,
    groupData = groupData,
    rawData = rawData
)

metadata(se)$normScore <- normScoreResults

6 Session information

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] NormalyzerDE_1.31.0         SummarizedExperiment_1.43.0
#>  [3] Biobase_2.73.2              GenomicRanges_1.65.4       
#>  [5] Seqinfo_1.3.2               IRanges_2.47.5             
#>  [7] S4Vectors_0.51.9            BiocGenerics_0.59.12       
#>  [9] generics_0.1.4              MatrixGenerics_1.25.0      
#> [11] matrixStats_1.5.0           normScore_0.99.1           
#> [13] BiocStyle_2.41.0           
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6          xfun_0.60             bslib_0.12.0         
#>  [4] ggplot2_4.0.3         rstatix_1.1.0         lattice_0.23-1       
#>  [7] vctrs_0.7.3           tools_4.6.1           tibble_3.3.1         
#> [10] vsn_3.81.0            pkgconfig_2.0.3       Matrix_1.7-6         
#> [13] RColorBrewer_1.1-3    S7_0.2.2              lifecycle_1.0.5      
#> [16] compiler_4.6.1        farver_2.1.2          statmod_1.5.2        
#> [19] tinytex_0.60          carData_3.0-6         htmltools_0.5.9      
#> [22] sass_0.4.10           yaml_2.3.12           Formula_1.2-6        
#> [25] preprocessCore_1.75.1 car_3.1-5             tidyr_1.3.2          
#> [28] ggpubr_1.0.0          pillar_1.11.1         jquerylib_0.1.4      
#> [31] MASS_7.3-66           affy_1.91.0           DelayedArray_0.39.6  
#> [34] cachem_1.1.0          limma_3.99.0          magick_2.9.1         
#> [37] boot_1.3-32           abind_1.4-8           tidyselect_1.2.1     
#> [40] digest_0.6.39         purrr_1.2.2           dplyr_1.2.1          
#> [43] bookdown_0.48         labeling_0.4.3        cowplot_1.2.0        
#> [46] fastmap_1.2.0         grid_4.6.1            cli_3.6.6            
#> [49] SparseArray_1.13.2    magrittr_2.0.5        S4Arrays_1.13.0      
#> [52] dichromat_2.0-1       broom_1.0.13          withr_3.0.3          
#> [55] backports_1.5.1       scales_1.4.0          rmarkdown_2.32       
#> [58] XVector_0.53.0        affyio_1.83.0         otel_0.2.0           
#> [61] ggsignif_0.6.4        evaluate_1.0.5        knitr_1.52           
#> [64] rlang_1.3.0           Rcpp_1.1.2            glue_1.8.1           
#> [67] BiocManager_1.30.27   jsonlite_2.0.0        R6_2.6.1