Comparing Two Groups with prolfqua

Witold E. Wolski

Contents

1 Purpose

This vignette demonstrates how two conditions, e.g., treatment vs. control, can be compared and the differences statistically tested. We use simulated protein abundances of an LFQ experiment with two treatment groups (A, B) and a control group (Ctrl). After first examining the data using QC plots and then normalizing the data, we compare each treatment group with the control. The output of the comparison is the difference in the mean intensities for quantified proteins (log2 fold-change) in each group, along with statistical parameters such as degrees of freedom, standard errors, p-value and the FDR.

2 Simulated protein abundances

library(prolfqua)

sim_lfq_data_protein_config() simulates protein abundances in long format, one row per protein and sample, and returns them together with a matching AnalysisConfiguration. To read your own data, see the vignette Creating configurations.

xx <- prolfqua::sim_lfq_data_protein_config(Nprot = 100)
xx

## $data
## # A tibble: 1,200 × 8
##    sample sampleName group_ isotopeLabel protein_Id  abundance qValue nrPeptides
##    <chr>  <chr>      <chr>  <chr>        <chr>           <dbl>  <dbl>      <dbl>
##  1 A_V1   A_V1       A      light        0EfVhX~3967      20.4      0          2
##  2 A_V1   A_V1       A      light        0YSKpy~2865      18.2      0          2
##  3 A_V1   A_V1       A      light        0m5WN4~6025      16.4      0          7
##  4 A_V1   A_V1       A      light        3QLHfm~8938      23.0      0          1
##  5 A_V1   A_V1       A      light        3QYop0~7543      23.1      0          7
##  6 A_V1   A_V1       A      light        76k03k~7094      25.5      0          4
##  7 A_V1   A_V1       A      light        7QuTub~1867      18.4      0          4
##  8 A_V1   A_V1       A      light        7cbcrd~7351      27.1      0          2
##  9 A_V1   A_V1       A      light        7soopj~5352      19.5      0          2
## 10 A_V1   A_V1       A      light        7zeekV~7127      20.6      0          8
## # ℹ 1,190 more rows
## 
## $config
## <AnalysisConfiguration>
##   Public:
##     annotation_vars: function () 
##     bin_resp: 
##     clone: function (deep = FALSE) 
##     factor_depth: 1
##     factor_keys: function () 
##     factor_keys_depth: function () 
##     factors: list
##     file_name: sample
##     get_response: function () 
##     hierarchy: list
##     hierarchy_depth: 1
##     hierarchy_keys: function () 
##     hierarchy_keys_depth: function (names = TRUE) 
##     id_vars: function () 
##     ident_q_value: qValue
##     ident_score: 
##     is_response_transformed: FALSE
##     isotope_label: isotopeLabel
##     min_peptides_protein: 2
##     norm_value: NULL
##     nr_children: nrPeptides
##     opt_mz: 
##     opt_rt: 
##     opt_se: 
##     pattern_contaminants: NULL
##     pattern_decoys: NULL
##     pop_response: function () 
##     sample_name: sampleName
##     sep: ~
##     set_response: function (col_name) 
##     value_vars: function () 
##     work_intensity: abundance

Create the LFQData class instance and remove zero or very small intensities, which most quantification software reports for missing values.

lfqdata <- prolfqua::LFQData$new(xx$data, xx$config)
lfqdata$remove_small_intensities()
lfqdata$factors()

## # A tibble: 12 × 3
##    sample  sampleName group_
##    <chr>   <chr>      <chr> 
##  1 A_V1    A_V1       A     
##  2 A_V2    A_V2       A     
##  3 A_V3    A_V3       A     
##  4 A_V4    A_V4       A     
##  5 B_V1    B_V1       B     
##  6 B_V2    B_V2       B     
##  7 B_V3    B_V3       B     
##  8 B_V4    B_V4       B     
##  9 Ctrl_V1 Ctrl_V1    Ctrl  
## 10 Ctrl_V2 Ctrl_V2    Ctrl  
## 11 Ctrl_V3 Ctrl_V3    Ctrl  
## 12 Ctrl_V4 Ctrl_V4    Ctrl

You can always convert the data into wide format.

lfqdata$data_wide()$data[1:3,1:7]

## # A tibble: 3 × 7
##   protein_Id  isotopeLabel  A_V1  A_V2  A_V3  A_V4  B_V1
##   <chr>       <chr>        <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 0EfVhX~3967 light         20.4  23.8  22.3  21.5  19.1
## 2 0YSKpy~2865 light         18.2  17.7  16.7  17.7  18.2
## 3 0m5WN4~6025 light         16.4  NA    16.1  18.2  24.3

2.1 Visualization of not normalized data

After this first setting up of the analysis we show now how to normalize the proteins and the effect of normalization. Furthermore we use some functions to visualize the missing values in our data.

lfqplotter <- lfqdata$get_Plotter()
density_nn <- lfqplotter$intensity_distribution_density()

2.1.1 Visualization of missing data

lfqplotter$na_heatmap()

Figure 1: Heatmap of missing protein intensities, black - missing, white - present

lfqdata$get_Summariser()$plot_missingness_per_group()

Figure 2: # of proteins with 0,1,…N missing values

lfqplotter$missigness_histogram()

Figure 3: Intensity distribution of proteins depending on # of missing values

2.1.2 Computing standard deviations, mean and CV.

Other important statistics such as coefficient of variation, means and standard deviations can be easily calculated using the get_Stats function and visualized with a violin plot.

stats <- lfqdata$get_Stats()
stats$violin()

Figure 4: Violin plots of CVs in the different groups and among all groups

prolfqua::table_facade( stats$stats_quantiles()$wide, paste0("quantile of ",stats$stat ))

Table 1: quantile of CV

probs	A	B	Ctrl	All
0.10	2.014611	1.993591	1.762209	3.451035
0.25	2.711702	2.956419	2.882970	4.370724
0.50	3.655067	3.700140	4.261204	5.296948
0.75	5.338899	5.462486	5.904953	7.411380
0.90	7.541564	7.118907	7.471854	9.912513

stats$density_median()

Figure 5: Distribution of CV’s for top 50% and bottom 50% proteins by intensity

2.2 Normalize protein intensities and show diagnostic plots

We normalize the data by \(\log_2\) transforming and then \(z-scaling\).

lt <- lfqdata$get_Transformer()
transformed <- lt$log2()$robscale()$lfq
transformed$is_transformed()

## [1] TRUE

pl <- transformed$get_Plotter()
density_norm <- pl$intensity_distribution_density()

gridExtra::grid.arrange(density_nn, density_norm)

Figure 6: Distribution of intensities before and after normalization

pl$pairs_smooth()

Figure 7: Scatterplot matrix

## NULL

p <- pl$heatmap_cor()

p

Figure 8: Heatmap, Rows - proteins, Columns - samples

3 Fitting a linear model

For fitting linear models to the transformed intensities for all our proteins we have to first specify the model function and define the contrasts that we want to calculate.

transformed$factors()

## # A tibble: 12 × 3
##    sample  sampleName group_
##    <chr>   <chr>      <chr> 
##  1 A_V1    A_V1       A     
##  2 A_V2    A_V2       A     
##  3 A_V3    A_V3       A     
##  4 A_V4    A_V4       A     
##  5 B_V1    B_V1       B     
##  6 B_V2    B_V2       B     
##  7 B_V3    B_V3       B     
##  8 B_V4    B_V4       B     
##  9 Ctrl_V1 Ctrl_V1    Ctrl  
## 10 Ctrl_V2 Ctrl_V2    Ctrl  
## 11 Ctrl_V3 Ctrl_V3    Ctrl  
## 12 Ctrl_V4 Ctrl_V4    Ctrl

formula_Condition <-  strategy_lm("transformedIntensity ~ group_")

# specify model definition
modelName  <- "Model"
contr_spec <- c("AvsC" = "group_A - group_Ctrl",
"BvsC" = "group_B - group_Ctrl")

Here we have to build the model for each protein.

mod <- prolfqua::build_model(
  transformed$data_long(),
  formula_Condition,
  subject_id = transformed$hierarchy_keys() )

In this plot we can see what factors in our model are mostly responsible for the adjusted p-values calculated from an analysis of variance.

mod$anova_histogram("FDR")

## $plot

Figure 9: Distribution of adjusted p-values (FDR)

## 
## $name
## [1] "Anova_p.values_Model.pdf"

One can also look at which proteins show different abundances in any of the groups by looking at the FDR values of the analysis of variance.

aovtable <- mod$get_anova()
head(aovtable)

## # A tibble: 6 × 10
##   protein_Id  factor    Df  Sum.Sq  Mean.Sq F.value  p.value isSingular nr_coef
##   <chr>       <chr>  <int>   <dbl>    <dbl>   <dbl>    <dbl> <lgl>        <int>
## 1 0EfVhX~3967 group_     2 0.0626  0.0313     5.18  0.0416   FALSE            3
## 2 0YSKpy~2865 group_     2 0.00753 0.00376    0.605 0.569    FALSE            3
## 3 0m5WN4~6025 group_     2 0.324   0.162     29.4   0.000205 FALSE            3
## 4 3QLHfm~8938 group_     2 0.0637  0.0319    13.2   0.00213  FALSE            3
## 5 3QYop0~7543 group_     2 0.00121 0.000603   0.178 0.840    FALSE            3
## 6 76k03k~7094 group_     2 0.0264  0.0132     3.06  0.0968   FALSE            3
## # ℹ 1 more variable: FDR <dbl>

dim(aovtable)

## [1] 98 10

xx <- aovtable |> dplyr::filter(FDR < 0.2)
signif <- transformed$get_copy()
signif$set_data(signif$data_long() |> dplyr::filter(protein_Id %in% xx$protein_Id))
hmSig <- signif$get_Plotter()$heatmap()

hmSig

Figure 10: Heatmap for proteins with FDR < 0.2 in the analysis of variance

4 Compute contrasts

Next we do calculate the statistics for our defined contrasts for all the proteins. For this we can use the Contrasts function.

contr <- prolfqua::Contrasts$new(mod, contr_spec)
v1 <- contr$get_Plotter()$volcano()

Alternatively, we can moderate the variance using the empirical Bayes method implemented in ContrastsModerated.

contr <- prolfqua::ContrastsModerated$new(contr)
contrdf <- contr$get_contrasts()

Moderation shrinks the variance of each protein towards the variance of all proteins, which stabilises the p-values of proteins whose variance is estimated from few observations. The next figure compares the volcano plots without and with moderation.

plotter <- contr$get_Plotter()
v2 <- plotter$volcano()
gridExtra::grid.arrange(v1$FDR,v2$FDR, ncol = 1)

Figure 11: Volcano plot, Left panel - no moderation, Right panel - with moderation

plotter$ma_plotly()

Figure 12: MA plot showing the dependency of mean abundance with respect to the difference

myProteinIDS <- head(unique(contrdf$protein_Id), 3)
dplyr::filter(contrdf, protein_Id %in% myProteinIDS)

## # A tibble: 6 × 16
##   modelName estimate_type protein_Id  contrast    diff avgAbd
##   <chr>     <chr>         <chr>       <chr>      <dbl>  <dbl>
## 1 WaldTest  observed      0EfVhX~3967 AvsC      0.0758   4.43
## 2 WaldTest  observed      0YSKpy~2865 AvsC     -0.0176   4.16
## 3 WaldTest  observed      0m5WN4~6025 AvsC     -0.260    4.21
## 4 WaldTest  observed      0EfVhX~3967 BvsC     -0.115    4.33
## 5 WaldTest  observed      0YSKpy~2865 BvsC     -0.0625   4.13
## 6 WaldTest  observed      0m5WN4~6025 BvsC      0.174    4.42
## # ℹ 10 more variables: std.error.unmoderated <dbl>, df.unmoderated <int>,
## #   std.error <dbl>, statistic <dbl>, df <dbl>, p.value <dbl>, conf.low <dbl>,
## #   conf.high <dbl>, sigma <dbl>, FDR <dbl>

5 Contrasts with missing value imputation

Use this method if there are proteins with no observations in one of the groups. The lm_impute facade refits the linear model for these proteins after imputing the missing observations with the limit of detection (LOD), and borrows the variance from the proteins that could be fitted.

mC <- prolfqua::ContrastsLMImputeFacade$new(transformed, "~ group_", contr_spec)$contrast
colnames(mC$get_contrasts())

##  [1] "modelName"             "estimate_type"         "protein_Id"           
##  [4] "contrast"              "diff"                  "avgAbd"               
##  [7] "std.error.unmoderated" "df.unmoderated"        "std.error"            
## [10] "statistic"             "df"                    "p.value"              
## [13] "conf.low"              "conf.high"             "sigma"                
## [16] "FDR"

Finally we are merging the results and give priority to the results where we do not have missing values in one group.

merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)$merged
plotter <- merged$get_Plotter()
tmp <- plotter$volcano()
tmp$FDR

Figure 13: Volcano plots for the two contrasts with missing value imputation from the lm_impute model

Look at proteins which could not be fitted using the linear model, if any.

merged <- prolfqua::merge_contrasts_results(prefer = contr,add = mC)

moreProt <- transformed$get_copy()
moreProt$set_data(moreProt$data_long() |> dplyr::filter(protein_Id %in% merged$more$contrast_result$protein_Id))
moreProt$get_Plotter()$raster()

# here we do not get anything because there is nothing imputed!

6 Gene Set Enrichment Analysis (GSEA)

After computing contrasts, proteins can be ranked by log2 fold change or t-statistic and subjected to gene set enrichment analysis. The typical workflow is:

1. Extract UniProt IDs from the contrast results using prolfqua::get_uniprot_id_from_fasta_header()
2. Map UniProt IDs to Entrez Gene IDs using AnnotationDbi::mapIds() with a species-specific annotation package
3. Create a named ranked list (statistic values named by gene IDs)
4. Run GSEA using clusterProfiler::gseGO() or clusterProfiler::gseKEGG()

For details see:

• clusterProfiler book – comprehensive guide to enrichment analysis in R
• clusterProfiler Bioconductor page
• AnnotationDbi – for mapping between gene identifier types

The prolfqua package is described in (Wolski et al. 2022).

7 Session Info

sessionInfo()

## R version 4.6.1 Patched (2026-06-24 r90190)
## Platform: x86_64-apple-darwin20
## Running under: macOS Ventura 13.7.8
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.6-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
## 
## locale:
## [1] C/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: America/New_York
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] prolfqua_0.99.2  dplyr_1.2.1      BiocStyle_2.41.0
## 
## loaded via a namespace (and not attached):
##   [1] Rdpack_2.6.6           gridExtra_2.3.1        sandwich_3.1-3        
##   [4] rlang_1.3.0            magrittr_2.0.5         clue_0.3-68           
##   [7] GetoptLong_1.1.1       multcomp_1.4-32        otel_0.2.0            
##  [10] matrixStats_1.5.0      compiler_4.6.1         mgcv_1.9-4            
##  [13] png_0.1-9              vctrs_0.7.3            pkgconfig_2.0.3       
##  [16] shape_1.4.6.1          crayon_1.5.3           fastmap_1.2.0         
##  [19] magick_2.9.1           backports_1.5.1        labeling_0.4.3        
##  [22] utf8_1.2.6             promises_1.5.0         rmarkdown_2.32        
##  [25] nloptr_2.2.1           UpSetR_1.4.1           tinytex_0.61          
##  [28] purrr_1.2.2            xfun_0.62              glmnet_5.1            
##  [31] jomo_2.7-6             logistf_1.26.1         cachem_1.1.0          
##  [34] jsonlite_2.0.0         progress_1.2.3         later_1.4.8           
##  [37] pan_2.0                prettyunits_1.2.0      broom_1.0.13          
##  [40] parallel_4.6.1         cluster_2.1.8.3        R6_2.6.1              
##  [43] stringi_1.8.9          bslib_0.12.0           RColorBrewer_1.1-3    
##  [46] limma_3.99.2           boot_1.3-32            rpart_4.1.27          
##  [49] jquerylib_0.1.4        estimability_2.0.0     Rcpp_1.1.2            
##  [52] bookdown_0.48          iterators_1.0.14       knitr_1.52            
##  [55] zoo_1.9-1              IRanges_2.47.5         httpuv_1.6.17         
##  [58] Matrix_1.7-6           splines_4.6.1          nnet_7.3-21           
##  [61] tidyselect_1.2.1       dichromat_2.0-1        yaml_2.3.12           
##  [64] doParallel_1.0.17      codetools_0.2-20       lattice_0.23-1        
##  [67] tibble_3.3.1           plyr_1.8.9             shiny_1.14.0          
##  [70] withr_3.0.3            S7_0.2.2               coda_0.19-4.1         
##  [73] evaluate_1.0.5         survival_3.8-12        circlize_0.4.18       
##  [76] pillar_1.11.1          BiocManager_1.30.27    mice_3.19.0           
##  [79] KernSmooth_2.23-27     foreach_1.5.2          stats4_4.6.1          
##  [82] reformulas_0.4.4       plotly_4.12.1          generics_0.1.4        
##  [85] hms_1.1.4              S4Vectors_0.51.10      ggplot2_4.0.3         
##  [88] scales_1.4.0           minqa_1.2.8            xtable_1.8-8          
##  [91] glue_1.8.1             lazyeval_0.2.3         emmeans_2.0.4         
##  [94] tools_4.6.1            data.table_1.18.6.1    lme4_2.0-6            
##  [97] forcats_1.0.1          mvtnorm_1.4-2          Cairo_1.7-0           
## [100] grid_4.6.1             tidyr_1.3.2            crosstalk_1.2.2       
## [103] rbibutils_2.4.1        colorspace_2.1-4       nlme_3.1-171          
## [106] formula.tools_1.7.1    cli_3.6.6              viridisLite_0.4.3     
## [109] ComplexHeatmap_2.29.0  gtable_0.3.6           sass_0.4.10           
## [112] digest_0.6.39          operator.tools_1.6.3.1 BiocGenerics_0.59.12  
## [115] ggrepel_0.9.8          TH.data_1.1-5          htmlwidgets_1.6.4     
## [118] rjson_0.2.23           farver_2.1.2           htmltools_0.5.9       
## [121] lifecycle_1.0.5        httr_1.4.9             mime_0.13             
## [124] GlobalOptions_0.1.4    mitml_0.4-5            statmod_1.5.2         
## [127] MASS_7.3-66

References

Wolski, Witold E., Paolo Nanni, Jonas Grossmann, Maria d’Errico, Ralph Schlapbach, and Christian Panse. 2022. “Prolfqua: A Comprehensive R-package for Proteomics Differential Expression Analysis.” bioRxiv, ahead of print. https://doi.org/10.1101/2022.06.07.494524.