QC and Sample Size Estimation
2026-11-04
Contents
1 Purpose
This vignette demonstrates how the package can be used to create descriptive statistics or quality control plots for a dataset. First, we read a “proteinGroups.txt” file from the MaxQuant software, and next, we specify a prolfqua configuration. Afterward, the data is transformed and normalized, and plots describing the data are generated. Finally, we estimate how many samples are necessary to properly quantify a two-fold change for \(90%\) of all the proteins with a power of \(0.8\).
2 Loading data
Generate Simulated data. To learn how to create prolfqua LFQData objects see the vignette: creating configurations.
simdata <- prolfqua::sim_lfq_data_peptide_config()
lfqdata <- prolfqua::LFQData$new(simdata$data,simdata$config)
lfqdata$remove_small_intensities()You can convert the data into a data frame in a wide format, where the intensities of each sample occupy their columns.
lfqdata$data_wide()$data[1:3,1:8]## # A tibble: 3 × 8
## protein_Id peptide_Id isotopeLabel A_V1 A_V2 A_V3 A_V4 B_V1
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0EfVhX~0087 ITLb4x1q light 19.4 NA 18.2 18.9 28.9
## 2 0EfVhX~0087 ahQLlQY7 light 26.9 25.2 26.1 25.2 26.1
## 3 0EfVhX~0087 dJkdz7so light NA 16.1 14.7 NA 25.72.1 Visualization of not normalized data
Next we show how the data is distributed before transformation and normalization.
lfqplotter <- lfqdata$get_Plotter()
lfqplotter$intensity_distribution_density()
2.1.1 Visualization of missing data
Often it also makes sense to look into the part of the data that is not quantified in all samples. Therefore several functions for the missingness are available. In the NA_heatmap function it is possible to see if some missing proteins are specific for one particular dilution or if they appear randomly in all samples.
nah <- lfqplotter$na_heatmap()Prints missing heatmap:
nah
Figure 1: Heatmap, black - missing protein intensities, white - present
Also in the next two figures we show how many missing values are found in the respective conditions and how many proteins do not have missing data and how complete the measurements are for each group.
lfqdata$get_Summariser()$plot_missingness_per_group()
Figure 2: # of proteins with 0,1,…N missing values
lfqdata$get_Summariser()$plot_missingness_per_group_cumsum()
Figure 3: Cumulative sum of the # of proteins with 0,1,…N missing values
In the missigness_histogram we can check if there is a dependency between NAs (and the number of NAs for each protein) with respect to the protein intensity.
lfqplotter$missigness_histogram()
Figure 4: Intensity distribution of proteins depending on # of missing values
3 Computing standard deviations, mean and CV.
Other important statistics that can be easily calculated from the object are coefficient of variation, means and standard deviations.
stats <- lfqdata$get_Stats()
prolfqua::table_facade( stats$stats_quantiles()$wide, paste0("quantile of ",stats$stat ))| probs | A | B | Ctrl | All |
|---|---|---|---|---|
| 0.10 | 1.215060 | 1.236542 | 1.129780 | 7.185163 |
| 0.25 | 2.442726 | 1.999453 | 1.607515 | 12.277918 |
| 0.50 | 3.418211 | 3.026062 | 3.349758 | 15.169313 |
| 0.75 | 4.599101 | 4.155821 | 4.885600 | 19.467593 |
| 0.90 | 6.349093 | 5.322302 | 6.008022 | 21.451325 |
stats$density_median()
Figure 5: Distribution of CV’s
In the next figure we show the dependency of the standard deviation with respect to the mean intensity of the proteins.
stdm_raw <- stats$stdv_vs_mean(size = 10000) +
ggplot2::scale_x_log10() +
ggplot2::scale_y_log10()
stdm_raw
Figure 6: Scatter plot of standard deviation vs mean
3.1 Normalize protein intensities
Next, we want to normalize the data by first \(\log_2\) transforming it and then z-scaling. The \(log_2\) stabilizes the variance, while the z-scaling removes systematic differences from the samples.
lt <- lfqdata$get_Transformer()
transformed <- lt$log2()$robscale()$lfq
transformed$is_transformed()## [1] TRUEAgain we want to look into the distribution of the intensities in our samples after normalization.
pl <- transformed$get_Plotter()
pl$intensity_distribution_density()
Figure 7: Normalized intensities
plots matrix scatter plot in the upper right matrix and some statistics in the lower left matrix.
pl$pairs_smooth()
Figure 8: Scatterplot matrix
## NULLplots protein heatmap:
p <- pl$heatmap()
p 
Figure 9: Heatmap, Rows - proteins, Columns - samples
We can also look at the correlation among the samples or look at a PCA plot.
hc <- pl$heatmap_cor()hc
Figure 10: Heatmap based on sample correlations, Rows - samples, Columns - samples
pl$pca()
Figure 11: Principal component analysis for all samples
prints the standard deviations:
stats <- transformed$get_Stats()
prolfqua::table_facade(stats$stats_quantiles()$wide, "Standard deviations")| probs | A | B | Ctrl | All |
|---|---|---|---|---|
| 0.10 | 0.0355196 | 0.0325012 | 0.0294799 | 0.1152932 |
| 0.25 | 0.0587524 | 0.0582281 | 0.0422670 | 0.1677154 |
| 0.50 | 0.0717438 | 0.0853025 | 0.0566645 | 0.2345489 |
| 0.75 | 0.1240478 | 0.1209993 | 0.0720238 | 0.3004859 |
| 0.90 | 0.1687990 | 0.1388753 | 0.1087199 | 0.3253740 |
plots density of the standard deviation
stats$density_median()
Figure 12: Density of the standard deviation
Check for heteroskedasticity. After transformation, the standard deviation should be independent of the mean intensity.
stdm_trans <- stats$stdv_vs_mean(size = 10000) + ggplot2::scale_x_log10() + ggplot2::scale_y_log10()
stdm_trans
Figure 13: Scatter plot of sd vs mean of protein intensity
#gridExtra::grid.arrange(stdm_raw, stdm_trans, nrow=1)3.2 Estimate sample size
For the sample size estimation we will now focus on only one condition and the different bio-reps from the a group In this condition all samples have the same concentration and we want to check variability within this group. Like this we can estimate how many samples are necessary per group to quantify properly a two-fold change (delta of 1) for 90% of all the proteins with a power of 0.8 and significance level of 0.05.
First, we need to filter our data for group_ A only.
transformedA <- transformed$get_copy()
transformedA$set_data(transformedA$data_long() |> dplyr::filter(group_ == "A"))
stats <- transformedA$get_Stats()stats$density(ggstat = "ecdf")
(#fig:plot_sd_ecdf)Empirical cumulative density function of the standard deviation for all the proteins in the dataset.
The table indicates that for this kind of data - if a two-fold change for 90% of all the proteins should be accurately quantified with a power of 0.8 (at 0.05 significance level) one would need to have 10 samples in each group.
sampleSize <- stats$power_t_test_quantiles() |>
dplyr::filter(group_ != "All")
prolfqua::table_facade(sampleSize, "Sample sizes. delta - Effect size, N - samplesize")| group_ | probs | quantiles | sdtrimmed | N_exact | N | delta |
|---|---|---|---|---|---|---|
| A | 0.10 | 0.0355196 | 0.0355196 | 1.526948 | 2 | 0.59 |
| A | 0.25 | 0.0587524 | 0.0587524 | 1.673109 | 2 | 0.59 |
| A | 0.50 | 0.0717438 | 0.0717438 | 1.758418 | 2 | 0.59 |
| A | 0.75 | 0.1240478 | 0.1240478 | 2.172950 | 3 | 0.59 |
| A | 0.90 | 0.1687990 | 0.1687990 | 2.671229 | 3 | 0.59 |
| A | 0.10 | 0.0355196 | 0.0530265 | 1.500024 | 2 | 1.00 |
| A | 0.25 | 0.0587524 | 0.0587524 | 1.521541 | 2 | 1.00 |
| A | 0.50 | 0.0717438 | 0.0717438 | 1.569649 | 2 | 1.00 |
| A | 0.75 | 0.1240478 | 0.1240478 | 1.768204 | 2 | 1.00 |
| A | 0.90 | 0.1687990 | 0.1687990 | 1.961551 | 2 | 1.00 |
| A | 0.10 | 0.0355196 | 0.1060530 | 1.500024 | 2 | 2.00 |
| A | 0.25 | 0.0587524 | 0.1060530 | 1.500024 | 2 | 2.00 |
| A | 0.50 | 0.0717438 | 0.1060530 | 1.500024 | 2 | 2.00 |
| A | 0.75 | 0.1240478 | 0.1240478 | 1.533718 | 2 | 2.00 |
| A | 0.90 | 0.1687990 | 0.1687990 | 1.616344 | 2 | 2.00 |
The table summarises visually how many samples are needed for different expected differences (in log2-scale) and different portions of all proteins.
sampleSize |>
ggplot2::ggplot(ggplot2::aes(x = factor(probs) , y = N)) +
ggplot2::facet_wrap(~delta) +
ggplot2::geom_bar(stat = "identity")
(#fig:plot_sample_size_quantiles)Estimated sample sizes for various FC levels and various quantiles of the standard deviation.
It is also possible to get the statistics for each protein in each group. The information includes group size, number of observations, standard deviation, variance, mean.
stats$stats() |> head()## # A tibble: 6 × 12
## protein_Id peptide_Id isotopeLabel group_ nrReplicates nrMeasured nrNAs
## <chr> <chr> <chr> <chr> <int> <int> <int>
## 1 0EfVhX~0087 ITLb4x1q light A 4 3 1
## 2 0EfVhX~0087 ahQLlQY7 light A 4 4 0
## 3 0EfVhX~0087 dJkdz7so light A 4 2 2
## 4 7cbcrd~5725 D5dQ4nKk light A 4 4 0
## 5 9VUkAq~4703 eIC06D7g light A 4 4 0
## 6 BEJI92~5282 HBkZvdhT light A 4 2 2
## # ℹ 5 more variables: sd <dbl>, var <dbl>, meanAbundance <dbl>,
## # medianAbundance <dbl>, interaction <chr>You can also estimate the sample size needed to perform differential expression analysis for each protein, given the power, delta, and sig.level.
The next figure shows the distribution of needed sample sizes for all proteins with a power of 0.8 at significance level 0.05 for a two fold change (delta = 1).
x <- stats$power_t_test(delta = 1,power = 0.8, sig.level = 0.05)
x <- x |> dplyr::filter(group_ == "A") |> dplyr::arrange(desc(N),protein_Id)
prolfqua::table_facade(x[1:7,], caption = "Sample size for each protein")| protein_Id | peptide_Id | isotopeLabel | group_ | nrReplicates | nrMeasured | nrNAs | sd | var | meanAbundance | medianAbundance | interaction | delta | N_exact | N |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 0EfVhX~0087 | dJkdz7so | light | A | 4 | 2 | 2 | 0.1816806 | 0.0330078 | 3.971469 | 3.971469 | group_A | 1 | 2.023422 | 3 |
| HvIpHG~9079 | opjydeWJ | light | A | 4 | 3 | 1 | 0.2392614 | 0.0572460 | 4.094835 | 3.997946 | group_A | 1 | 2.344670 | 3 |
| JcKVfU~9653 | 9LMQBevj | light | A | 4 | 4 | 0 | 0.1871014 | 0.0350069 | 5.482917 | 5.406550 | group_A | 1 | 2.050429 | 3 |
| 0EfVhX~0087 | ITLb4x1q | light | A | 4 | 3 | 1 | 0.1535287 | 0.0235710 | 4.190351 | 4.178250 | group_A | 1 | 1.892113 | 2 |
| 0EfVhX~0087 | ahQLlQY7 | light | A | 4 | 4 | 0 | 0.0620097 | 0.0038452 | 4.742744 | 4.735435 | group_A | 1 | 1.533666 | 2 |
| 7cbcrd~5725 | D5dQ4nKk | light | A | 4 | 4 | 0 | 0.0508679 | 0.0025875 | 4.956022 | 4.964554 | group_A | 1 | 1.500024 | 2 |
| 9VUkAq~4703 | eIC06D7g | light | A | 4 | 4 | 0 | 0.0945821 | 0.0089458 | 4.159728 | 4.182873 | group_A | 1 | 1.654299 | 2 |
x |> ggplot2::ggplot(ggplot2::aes(x = N)) +
ggplot2::geom_histogram() +
ggplot2::facet_wrap(~delta) +
ggplot2::xlim(0,100)
Figure 14: Distribution of the required sample sizes for two fold change thresholds for all the proteins
The prolfqua package is described in (Wolski et al. 2023).
4 Session Info
sessionInfo()## R version 4.6.1 (2026-06-24)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.5 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] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_4.0.3 prolfqua_0.99.8 dplyr_1.2.1 BiocStyle_2.41.0
##
## loaded via a namespace (and not attached):
## [1] RColorBrewer_1.1-3 jsonlite_2.0.0 shape_1.4.6.1
## [4] magrittr_2.0.5 TH.data_1.1-5 estimability_2.0.0
## [7] jomo_2.7-6 magick_2.9.1 farver_2.1.2
## [10] logistf_1.26.1 nloptr_2.2.1 rmarkdown_2.32
## [13] GlobalOptions_0.1.4 vctrs_0.7.3 Cairo_1.7-0
## [16] minqa_1.2.8 tinytex_0.61 htmltools_0.5.9
## [19] forcats_1.0.1 progress_1.2.3 broom_1.0.13
## [22] mitml_0.4-5 sass_0.4.10 KernSmooth_2.23-27
## [25] bslib_0.12.0 htmlwidgets_1.6.4 plyr_1.8.9
## [28] sandwich_3.1-3 emmeans_2.0.4 plotly_4.12.1
## [31] zoo_1.9-1 cachem_1.1.0 mime_0.13
## [34] lifecycle_1.0.5 iterators_1.0.14 pkgconfig_2.0.3
## [37] Matrix_1.7-6 R6_2.6.1 fastmap_1.2.0
## [40] rbibutils_2.4.1 shiny_1.14.0 clue_0.3-68
## [43] digest_0.6.39 numDeriv_2016.8-1.1 colorspace_2.1-4
## [46] S4Vectors_0.51.10 crosstalk_1.2.2 labeling_0.4.3
## [49] httr_1.4.9 mgcv_1.9-5 compiler_4.6.1
## [52] withr_3.0.3 doParallel_1.0.17 S7_0.2.2
## [55] backports_1.5.1 UpSetR_1.4.1 pan_2.0
## [58] MASS_7.3-66 rjson_0.2.23 tools_4.6.1
## [61] otel_0.2.0 httpuv_1.6.17 nnet_7.3-21
## [64] glue_1.8.1 nlme_3.1-172 promises_1.5.0
## [67] grid_4.6.1 cluster_2.1.8.3 generics_0.1.4
## [70] operator.tools_1.6.3.1 gtable_0.3.6 formula.tools_1.7.1
## [73] tidyr_1.3.2 data.table_1.18.6.1 hms_1.1.4
## [76] utf8_1.2.6 BiocGenerics_0.59.12 ggrepel_0.9.8
## [79] foreach_1.5.2 pillar_1.11.1 limma_3.99.4
## [82] later_1.4.8 circlize_0.4.18 splines_4.6.1
## [85] lattice_0.23-1 survival_3.8-12 tidyselect_1.2.1
## [88] ComplexHeatmap_2.29.0 knitr_1.52 reformulas_0.4.4
## [91] gridExtra_2.3.1 bookdown_0.48 IRanges_2.47.5
## [94] stats4_4.6.1 xfun_0.62 statmod_1.5.2
## [97] matrixStats_1.5.0 stringi_1.8.9 yaml_2.3.12
## [100] boot_1.3-32 evaluate_1.0.5 codetools_0.2-20
## [103] Rfit_0.27.0 tibble_3.3.1 BiocManager_1.30.27
## [106] cli_3.6.6 rpart_4.1.27 xtable_1.8-8
## [109] Rdpack_2.6.6 jquerylib_0.1.4 dichromat_2.0-1
## [112] Rcpp_1.1.2 coda_0.19-4.1 png_0.1-9
## [115] parallel_4.6.1 prettyunits_1.2.0 lme4_2.1-0
## [118] glmnet_5.2 viridisLite_0.4.3 limpa_1.5.2
## [121] mvtnorm_1.4-2 lmerTest_3.2-1 scales_1.4.0
## [124] purrr_1.2.2 crayon_1.5.3 GetoptLong_1.1.1
## [127] rlang_1.3.0 multcomp_1.4-32 mice_3.19.0