1 The problem

Heatmaps, PCA and clustering need a complete data matrix. Label-free data are not complete: a protein is often not observed in some samples, and often not in any sample of one group. While such proteins are frequently below the limit of detection (LOD), complete absence can also reflect true biological on/off expression, which makes them candidates of high biological interest. PCA drops every protein with a missing value, and with it exactly these proteins.

library(prolfqua)
library(ggplot2)
sim <- sim_lfq_data_protein_config(Nprot = 100, weight_missing = 0.5, seed = 1234)
lfq <- LFQData$new(sim$data, sim$config)

wide <- lfq$data_wide(as.matrix = TRUE)$data
empty_group <- lfq$data_long() |>
  dplyr::group_by(protein_Id, group_) |>
  dplyr::summarise(empty = all(is.na(abundance)), .groups = "drop") |>
  dplyr::filter(empty)
c(
  proteins = nrow(wide),
  complete = sum(stats::complete.cases(wide)),
  with_empty_group = dplyr::n_distinct(empty_group$protein_Id)
)

##         proteins         complete with_empty_group 
##              100               23                5

The proteins with an empty group:

lfq$data_long() |>
  dplyr::filter(protein_Id %in% empty_group$protein_Id) |>
  dplyr::select(protein_Id, sampleName, abundance) |>
  tidyr::pivot_wider(names_from = sampleName, values_from = abundance) |>
  knitr::kable(digits = 1)

protein_Id	A_V1	A_V2	A_V3	A_V4	B_V1	B_V2	B_V3	B_V4	Ctrl_V1	Ctrl_V2	Ctrl_V3	Ctrl_V4
DoWup2~8058	17.9	NA	18.9	NA	NA	NA	NA	NA	NA	NA	NA	17
LMQBev~8107	NA	16.0	NA	NA	NA	15.8	16.0	15.8	NA	NA	NA	NA
ZHBkZv~0394	NA	NA	NA	NA	NA	NA	NA	13.3	NA	NA	12.6	NA
cnCfON~9769	NA	15.9	14.3	14.4	NA	NA	NA	NA	NA	15.8	17.5	NA
umjHSH~2090	NA	NA	NA	NA	NA	NA	15.8	16.5	NA	NA	17.4	NA

2 Fold changes are a separate question

To test for differences you do not need to impute anything yourself: build_contrast_analysis(method = "lm_impute") also estimates fold changes for proteins with an empty group, refitting them at the LOD and borrowing their variance from the other proteins. See Comparing two groups and Comparing modelling methods on one dataset. Here we reuse the models of that fit to fill the data matrix.

3 Model-based imputation

impute_from_model() fills every missing cell with the prediction of its protein model and keeps the observed values. The summary tells which route each protein took:

complete: no missing values, nothing to fill
fitted: the model fitted on the observed values predicts the missing cells from their group
lod_refit: a protein whose fit failed, typically because a group is empty, is refitted with its missing values set to the LOD and observed values below the LOD raised to it; the empty group gets a value at the LOD

contrasts <- c(AvsC = "group_A - group_Ctrl", BvsC = "group_B - group_Ctrl")
fa <- build_contrast_analysis(lfq, "~ group_", contrasts, method = "lm_impute")
filled <- impute_from_model(fa$model, lfq)
table(filled$summary$route)

## 
##  complete    fitted lod_refit 
##        23        70         7

MissingHelpers estimates the LOD from the data. Groups in which a protein is observed in only one sample, or in the fewest samples, are the closest the data comes to the detection limit; the LOD is the median of their mean intensities. Pass lod to build_contrast_analysis() to use a LOD of your own.

MissingHelpers$new(lfq$data_long(), lfq$get_config())$get_lod()

## [1] 15.97769

The proteins with an empty group after imputation; compare with the table above:

filled$lfqdata$data_long() |>
  dplyr::filter(protein_Id %in% empty_group$protein_Id) |>
  dplyr::select(protein_Id, sampleName, abundance) |>
  tidyr::pivot_wider(names_from = sampleName, values_from = abundance) |>
  knitr::kable(digits = 1)

protein_Id	A_V1	A_V2	A_V3	A_V4	B_V1	B_V2	B_V3	B_V4	Ctrl_V1	Ctrl_V2	Ctrl_V3	Ctrl_V4
DoWup2~8058	17.9	17.6	18.9	17.6	16.0	16.0	16.0	16.0	16.4	16.4	16.4	17.0
LMQBev~8107	16.0	16.0	16.0	16.0	16.0	15.8	16.0	15.8	16.0	16.0	16.0	16.0
ZHBkZv~0394	16.0	16.0	16.0	16.0	16.0	16.0	16.0	13.3	16.0	16.0	12.6	16.0
cnCfON~9769	16.0	15.9	14.3	14.4	16.0	16.0	16.0	16.0	16.5	15.8	17.5	16.5
umjHSH~2090	16.0	16.0	16.0	16.0	16.2	16.2	15.8	16.5	16.9	16.9	17.4	16.9

4 Heatmap before and after

We cluster the proteins once, on the imputed data, and draw both heatmaps with that row order, so each row shows the same protein before and after.

before <- lfq$data_wide(as.matrix = TRUE)
after <- filled$lfqdata$data_wide(as.matrix = TRUE)
row_order <- rownames(after$data)[hclust(dist(t(scale(t(after$data)))))$order]

ht_before <- plot_heatmap(
  before$data, before$annotation, lfq$factor_keys(), lfq$sample_name(),
  na_fraction = 1, row_order = row_order
)
ht_before

Figure 1: Before imputation: missing cells are grey

ht_after <- plot_heatmap(
  after$data, after$annotation, lfq$factor_keys(), lfq$sample_name(),
  na_fraction = 1, row_order = row_order
)
ht_after

Figure 2: After imputation: every cell holds a value

5 PCA before and after

Before imputation PCA can only use the proteins without missing values; after imputation it uses all proteins.

n_complete <- sum(filled$summary$route == "complete")
gridExtra::grid.arrange(
  lfq$get_Plotter()$pca() +
    ggtitle(paste0("Before imputation: ", n_complete, " complete proteins")),
  filled$lfqdata$get_Plotter()$pca() +
    ggtitle(paste0("After imputation: ", nrow(filled$summary), " proteins")),
  ncol = 1
)

Figure 3: PCA before (top) and after (bottom) imputation

6 Caveats

After imputation the groups separate more cleanly, in the heatmap and in the PCA. Part of that is circular: each filled value is the model’s estimate for its group, so imputation adds the group structure of the model to the data. The filled values also carry no noise. Use the filled matrix to look at the data; for tests use the contrasts of lm_impute, which account for the imputation.

7 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] 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-5            
##  [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.2            
##  [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.4           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      scales_1.4.0          
##  [88] minqa_1.2.8            xtable_1.8-8           glue_1.8.1            
##  [91] emmeans_2.0.4          tools_4.6.1            data.table_1.18.6.1   
##  [94] lme4_2.1-0             forcats_1.0.1          mvtnorm_1.4-2         
##  [97] Cairo_1.7-0            grid_4.6.1             tidyr_1.3.2           
## [100] crosstalk_1.2.2        rbibutils_2.4.1        colorspace_2.1-4      
## [103] nlme_3.1-172           formula.tools_1.7.1    cli_3.6.6             
## [106] viridisLite_0.4.3      ComplexHeatmap_2.29.0  gtable_0.3.6          
## [109] sass_0.4.10            digest_0.6.39          operator.tools_1.6.3.1
## [112] BiocGenerics_0.59.12   ggrepel_0.9.8          TH.data_1.1-5         
## [115] htmlwidgets_1.6.4      rjson_0.2.23           farver_2.1.2          
## [118] htmltools_0.5.9        lifecycle_1.0.5        httr_1.4.9            
## [121] mime_0.13              GlobalOptions_0.1.4    mitml_0.4-5           
## [124] statmod_1.5.2          MASS_7.3-66

Filling the gaps: model-based imputation for heatmaps and PCA

2026-11-04

Contents