Overview

This script builds Figure 3: missingness, completeness, and repeatability checks at the feature and sample level, prior to normalization and batch correction (covered separately in Figure 4).

Two panels are adapted based on what this dataset actually supports:

library(dplyr)
library(ggplot2)
library(patchwork)
proteomics_raw <- readxl::read_excel("Copy_of_Proteomics_data_original_1.xlsx",
                                      sheet = "proteins_imputed_80 percent_V3") %>%
  janitor::clean_names() %>%
  mutate(batch = stringr::str_extract(file_name, "^B\\d+"))

metabolomics_raw <- readxl::read_excel("NMR_Metabolomics_Data__main_.xlsx", sheet = "Sheet1") %>%
  janitor::clean_names()

cat("Proteomics QC/pool-labelled samples found:",
    nrow(proteomics_raw %>% filter(grepl("QC|pool", linked_id, ignore.case = TRUE) |
                                      grepl("QC|pool", file_name, ignore.case = TRUE))), "\n")
## Proteomics QC/pool-labelled samples found: 0
cat("Metabolomics QC/pool-labelled samples found:",
    nrow(metabolomics_raw %>% filter(grepl("QC|pool", nmr_id, ignore.case = TRUE))), "\n")
## Metabolomics QC/pool-labelled samples found: 0

Panel A: Metabolomics missingness

zero_summary <- metabolomics_raw %>%
  select(-nmr_id, -group) %>%
  summarise(across(everything(), ~ sum(. == 0, na.rm = TRUE))) %>%
  tidyr::pivot_longer(everything(), names_to = "metabolite", values_to = "n_zero") %>%
  mutate(pct_zero = 100 * n_zero / nrow(metabolomics_raw),
         status = ifelse(pct_zero > 20, "Excluded", "Retained")) %>%
  arrange(desc(pct_zero))
zero_summary$metabolite <- factor(zero_summary$metabolite, levels = zero_summary$metabolite)

pA <- ggplot(zero_summary, aes(x = metabolite, y = pct_zero, fill = status)) +
  geom_col(width = 0.7) +
  geom_hline(yintercept = 20, linetype = "dashed", color = "grey30", linewidth = 0.5) +
  scale_fill_manual(values = c(Excluded = "#C44E52", Retained = "#4C72B0")) +
  labs(x = NULL, y = "% below detection", title = "A: Metabolomics missingness", fill = NULL) +
  theme_minimal(base_size = 13) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5, size = 10),
        axis.text.y = element_text(size = 11),
        plot.title = element_text(size = 15, face = "bold"), legend.position = "top",
        legend.text = element_text(size = 11))
pA

Panel B: Proteomics sample coverage by batch

batch_counts <- proteomics_raw %>% count(batch)
pB <- ggplot(batch_counts, aes(x = batch, y = n, fill = batch)) +
  geom_col(width = 0.6) +
  geom_text(aes(label = n), vjust = -0.4, size = 4) +
  scale_fill_manual(values = c(B1 = "#4C72B0", B2 = "#DD8452", B3 = "#55A868"), guide = "none") +
  labs(x = NULL, y = "n samples", title = "B: Proteomics coverage by batch\n(0 missing/1221 proteins)") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(size = 13, face = "bold"), axis.text = element_text(size = 11)) +
  ylim(0, max(batch_counts$n) * 1.2)
pB

Panel C: Technical replicate repeatability

replicate_patients <- c("HC14", "P067", "P070")
prot_cols2 <- proteomics_raw %>% select(-linked_id, -pathology, -file_name, -batch) %>% names()

cv_results <- lapply(replicate_patients, function(pid) {
  reps <- proteomics_raw %>% filter(linked_id == pid) %>% select(all_of(prot_cols2))
  cv <- apply(reps, 2, function(x) 100 * sd(x) / mean(x))
  data.frame(patient = pid, cv = cv)
}) %>% bind_rows()

cat("Median absolute %CV per replicate pair:\n")
## Median absolute %CV per replicate pair:
cv_results %>% group_by(patient) %>% summarise(median_cv = median(abs(cv), na.rm = TRUE))
## # A tibble: 3 × 2
##   patient median_cv
##   <chr>       <dbl>
## 1 HC14         51.8
## 2 P067         58.6
## 3 P070         46.9
pC <- ggplot(cv_results, aes(x = patient, y = abs(cv), fill = patient)) +
  geom_boxplot(width = 0.5, outlier.size = 0.4, outlier.alpha = 0.3) +
  scale_fill_manual(values = c(HC14 = "#4C72B0", P067 = "#DD8452", P070 = "#55A868"), guide = "none") +
  labs(x = NULL, y = "Absolute %CV", title = "C: Technical replicate CV\n(3 patients, 1221 proteins)") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(size = 13, face = "bold"), axis.text = element_text(size = 11)) +
  coord_cartesian(ylim = c(0, quantile(abs(cv_results$cv), 0.98, na.rm = TRUE)))
pC

This finding is flagged plainly, not smoothed over: median absolute %CV between the two replicate runs is 51.8% (HC14), 58.6% (P067), and 46.9% (P070) – well above the ~20-30% typically considered acceptable reproducibility for large-scale DIA-MS proteomics, and far above what would be expected for a well-controlled technical replicate. This does not necessarily indicate a fault in this pipeline’s handling of the data (both runs were averaged as documented in 01_Proteomics_Import_QC.Rmd), but it does mean the averaged values for these three specific patients carry more measurement uncertainty than a typical single-run sample, and it is worth raising with the supervisor whether this level of run-to-run variability is consistent with what is expected for this specific DIA-MS acquisition protocol, or whether it warrants closer inspection of these three patients’ raw spectra.

Combined Figure 3

figure3 <- pA / (pB | pC) + plot_layout(heights = c(1.3, 1)) +
  plot_annotation(title = "Figure 3. Missingness, completeness, and repeatability",
                   theme = theme(plot.title = element_text(face = "bold", size = 16)))
figure3

ggsave("Figure3_Feature_Sample_QC.png", figure3, width = 11, height = 9, dpi = 200, bg = "white")

Results text (draft)

Feature-level completeness, sample-level coverage, and technical repeatability were assessed prior to normalization and batch correction (Figure 3; batch correction itself is demonstrated separately in Figure 4). At the feature level, 8 of 37 quantified metabolites were excluded for falling below an 80% detection threshold (Figure 3A), while all 1221 proteins were fully quantified with no missing values. Proteomic batches were approximately balanced in size (53-63 samples per batch; Figure 3B). Technical repeatability, assessed using the three genuine replicate pairs identified during import QC (no pooled QC samples were included in this dataset’s acquisition), showed a median absolute %CV of 46.9-58.6% between replicate runs (Figure 3C) – higher than typically expected, and noted here for discussion with the supervisor regarding acceptable variability for this acquisition protocol.

Session info

sessionInfo()
## R version 4.6.1 (2026-06-24 ucrt)
## Platform: x86_64-w64-mingw32/x64
## Running under: Windows 11 x64 (build 26200)
## 
## Matrix products: default
##   LAPACK version 3.12.1
## 
## locale:
## [1] LC_COLLATE=English_United States.utf8 
## [2] LC_CTYPE=English_United States.utf8   
## [3] LC_MONETARY=English_United States.utf8
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.utf8    
## 
## time zone: Africa/Johannesburg
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] patchwork_1.3.2 ggplot2_4.0.3   dplyr_1.2.1    
## 
## loaded via a namespace (and not attached):
##  [1] sass_0.4.10        utf8_1.2.6         generics_0.1.4     tidyr_1.3.2       
##  [5] stringi_1.8.7      digest_0.6.39      magrittr_2.0.5     evaluate_1.0.5    
##  [9] grid_4.6.1         timechange_0.4.0   RColorBrewer_1.1-3 fastmap_1.2.0     
## [13] cellranger_1.1.0   jsonlite_2.0.0     purrr_1.2.2        scales_1.4.0      
## [17] textshaping_1.0.5  jquerylib_0.1.4    cli_3.6.6          rlang_1.2.0       
## [21] withr_3.0.3        cachem_1.1.0       yaml_2.3.12        otel_0.2.0        
## [25] tools_4.6.1        vctrs_0.7.3        R6_2.6.1           lifecycle_1.0.5   
## [29] lubridate_1.9.5    snakecase_0.11.1   stringr_1.6.0      ragg_1.5.2        
## [33] janitor_2.2.1      pkgconfig_2.0.3    pillar_1.11.1      bslib_0.11.0      
## [37] gtable_0.3.6       glue_1.8.1         systemfonts_1.3.2  xfun_0.59         
## [41] tibble_3.3.1       tidyselect_1.2.1   rstudioapi_0.19.0  knitr_1.51        
## [45] farver_2.1.2       htmltools_0.5.9    rmarkdown_2.31     labeling_0.4.3    
## [49] compiler_4.6.1     S7_0.2.2           readxl_1.5.0