Objective 1 — Step 10: Normalization and Batch Correction (Figure 4)
Judy Ngobeni
2026-09-17
Overview
This script builds Figure 4, demonstrating the effect of ComBat batch correction directly – before vs. after – for both proteomics and metabolomics side by side, using PCA and sample intensity distributions as the two lines of evidence. This is distinct from Figure 3, which covers missingness, completeness, and repeatability rather than the correction step itself.
library(dplyr)
library(ggplot2)
library(patchwork)standardize_id <- function(id) {
id <- stringr::str_trim(id)
prefix <- stringr::str_extract(id, "^[A-Za-z]+")
number <- stringr::str_extract(id, "\\d+")
paste0(toupper(prefix), sprintf("%03d", as.integer(number)))
}
matched_cohort_ids <- readRDS("matched_cohort_ids.rds")Proteomics: before and after batch correction
proteomics_clean <- readRDS("proteomics_clean.rds")
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+"), id_std = standardize_id(linked_id))
proteomics_std <- proteomics_clean %>%
mutate(id_std = standardize_id(linked_id)) %>%
filter(id_std %in% matched_cohort_ids) %>%
left_join(proteomics_raw %>% select(id_std, batch) %>% distinct(id_std, .keep_all = TRUE), by = "id_std")
pcols <- proteomics_std %>% select(-id_std, -linked_id, -pathology, -batch) %>% names()
pmat_pre <- proteomics_std %>% select(all_of(pcols)) %>% mutate(across(everything(), ~log2(. + 1))) %>% as.matrix()
pca_p_pre <- prcomp(pmat_pre, scale. = TRUE)
scores_p_pre <- as.data.frame(pca_p_pre$x[, 1:2]) %>% mutate(batch = proteomics_std$batch)
harmonized_proteomics <- readRDS("harmonized_proteomics.rds")
hp <- harmonized_proteomics %>% left_join(proteomics_std %>% select(id_std, batch), by = "id_std")
pmat_post <- hp %>% select(-id_std, -batch) %>% as.matrix()
pca_p_post <- prcomp(pmat_post, scale. = TRUE)
scores_p_post <- as.data.frame(pca_p_post$x[, 1:2]) %>% mutate(batch = hp$batch)
prot_intensity_pre <- proteomics_std %>% rowwise() %>%
mutate(med = median(c_across(all_of(pcols)), na.rm = TRUE)) %>% ungroup() %>% select(batch, med)
prot_intensity_post <- hp %>% rowwise() %>%
mutate(med = median(c_across(-c(id_std, batch)), na.rm = TRUE)) %>% ungroup() %>% select(batch, med)Metabolomics: before and after batch correction
metabolomics_clean <- readRDS("metabolomics_clean.rds")
batch_info <- readxl::read_excel("NMR_Metabolomics_Data__main_.xlsx", sheet = "Sheet2") %>% janitor::clean_names()
metabolomics_std <- metabolomics_clean %>%
left_join(batch_info %>% select(nmr_id, id, batch), by = "nmr_id") %>%
mutate(id_std = standardize_id(id)) %>%
filter(id_std %in% matched_cohort_ids)
mcols <- metabolomics_std %>% select(-nmr_id, -id, -batch, -group, -id_std) %>% names()
mmat_pre <- metabolomics_std %>% select(all_of(mcols)) %>% as.matrix()
pca_m_pre <- prcomp(mmat_pre, scale. = TRUE)
scores_m_pre <- as.data.frame(pca_m_pre$x[, 1:2]) %>% mutate(batch = metabolomics_std$batch)
harmonized_metabolomics <- readRDS("harmonized_metabolomics.rds")
hm <- harmonized_metabolomics %>% left_join(metabolomics_std %>% select(id_std, batch), by = "id_std")
mmat_post <- hm %>% select(-id_std, -batch) %>% as.matrix()
pca_m_post <- prcomp(mmat_post, scale. = TRUE)
scores_m_post <- as.data.frame(pca_m_post$x[, 1:2]) %>% mutate(batch = hm$batch)
metab_intensity_pre <- metabolomics_std %>% rowwise() %>%
mutate(med = median(c_across(all_of(mcols)), na.rm = TRUE)) %>% ungroup() %>% select(batch, med)
metab_intensity_post <- hm %>% rowwise() %>%
mutate(med = median(c_across(-c(id_std, batch)), na.rm = TRUE)) %>% ungroup() %>% select(batch, med)Combined Figure 4
prot_colors <- c(B1 = "#4C72B0", B2 = "#DD8452", B3 = "#55A868")
metab_colors <- c(`Batch 1` = "#4C72B0", `Batch 2` = "#DD8452")
mk_pca <- function(df, colors, title) {
ggplot(df, aes(x = PC1, y = PC2, color = batch)) +
geom_point(size = 1.3, alpha = 0.85) +
scale_color_manual(values = colors) +
labs(title = title, color = "Batch") +
theme_minimal(base_size = 8) +
theme(plot.title = element_text(size = 8.5, face = "bold"), legend.position = "right",
legend.text = element_text(size = 6), legend.title = element_text(size = 7))
}
mk_box <- function(df, colors, title) {
ggplot(df, aes(x = batch, y = med, fill = batch)) +
geom_boxplot(width = 0.5, outlier.size = 0.5) +
scale_fill_manual(values = colors, guide = "none") +
labs(x = NULL, y = "Median value", title = title) +
theme_minimal(base_size = 8) + theme(plot.title = element_text(size = 8.5, face = "bold"))
}
pA <- mk_pca(scores_p_pre, prot_colors, "A: Proteomics PCA - before correction")
pB <- mk_pca(scores_p_post, prot_colors, "B: Proteomics PCA - after correction")
pC <- mk_box(prot_intensity_pre, prot_colors, "C: Proteomics intensity - before")
pD <- mk_box(prot_intensity_post, prot_colors, "D: Proteomics intensity - after")
pE <- mk_pca(scores_m_pre, metab_colors, "E: Metabolomics PCA - before correction")
pF <- mk_pca(scores_m_post, metab_colors, "F: Metabolomics PCA - after correction")
pG <- mk_box(metab_intensity_pre, metab_colors, "G: Metabolomics intensity - before")
pH <- mk_box(metab_intensity_post, metab_colors, "H: Metabolomics intensity - after")
figure4 <- (pA | pB | pC | pD) / (pE | pF | pG | pH) +
plot_annotation(
title = "Figure 4. Normalization and batch correction: before vs. after (proteomics and metabolomics)",
theme = theme(plot.title = element_text(face = "bold", size = 13))
)
figure4
ggsave("Figure4_Batch_Correction.png", figure4, width = 15, height = 7.5, dpi = 200, bg = "white")Results text (draft)
The effect of ComBat batch correction was assessed directly for both platforms (Figure 4). In proteomics, sample intensity distributions showed a clear pre-correction batch effect, with batch 2 running higher than batches 1 and 3 (Figure 4C); this difference was substantially reduced following correction, with all three batches converging to a comparable, near-zero-centred distribution (Figure 4D). Proteomic PCA showed batches already reasonably intermixed on the first two components both before and after correction (Figure 4A-B), indicating that the batch effect visible in the intensity distributions was not the dominant source of variance captured by PC1/PC2. In metabolomics, the two acquisition batches showed only a modest difference in intensity distribution before correction (Figure 4G), consistent with less pronounced batch drift than proteomics, and remained similarly comparable after correction (Figure 4H); PCA showed intermixed batches both before and after (Figure 4E-F).
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] gtable_0.3.6 jsonlite_2.0.0 compiler_4.6.1 tidyselect_1.2.1
## [5] stringr_1.6.0 snakecase_0.11.1 jquerylib_0.1.4 textshaping_1.0.5
## [9] systemfonts_1.3.2 scales_1.4.0 readxl_1.5.0 yaml_2.3.12
## [13] fastmap_1.2.0 R6_2.6.1 labeling_0.4.3 generics_0.1.4
## [17] knitr_1.51 tibble_3.3.1 janitor_2.2.1 lubridate_1.9.5
## [21] bslib_0.11.0 pillar_1.11.1 RColorBrewer_1.1-3 rlang_1.2.0
## [25] stringi_1.8.7 cachem_1.1.0 xfun_0.59 sass_0.4.10
## [29] S7_0.2.2 otel_0.2.0 timechange_0.4.0 cli_3.6.6
## [33] withr_3.0.3 magrittr_2.0.5 digest_0.6.39 grid_4.6.1
## [37] rstudioapi_0.19.0 lifecycle_1.0.5 vctrs_0.7.3 evaluate_1.0.5
## [41] glue_1.8.1 cellranger_1.1.0 farver_2.1.2 ragg_1.5.2
## [45] rmarkdown_2.31 tools_4.6.1 pkgconfig_2.0.3 htmltools_0.5.9