This script builds Table 4 (structural quality and completeness across the three datasets) and Figure 3A (a patient-by-variable missingness map for the clinical dataset), using the raw source files and cleaned checkpoints from the earlier import/QC steps.
library(dplyr)
library(readxl)
library(janitor)
library(ggplot2)
clinical_clean <- readRDS("clinical_clean.rds")
proteomics_raw <- read_excel("Copy_of_Proteomics_data_original_1.xlsx",
sheet = "proteins_imputed_80 percent_V3") %>% clean_names()
proteomics_clean <- readRDS("proteomics_clean.rds")
metabolomics_raw <- read_excel("NMR_Metabolomics_Data__main_.xlsx", sheet = "Sheet1") %>% clean_names()
metabolomics_clean <- readRDS("metabolomics_clean.rds")
cat("Clinical: ", nrow(clinical_clean), "rows\n")
## Clinical: 233 rows
cat("Proteomics: raw =", nrow(proteomics_raw), " final =", nrow(proteomics_clean), "\n")
## Proteomics: raw = 178 final = 175
cat("Metabolomics: raw =", nrow(metabolomics_raw), " final =", nrow(metabolomics_clean), "\n")
## Metabolomics: raw = 81 final = 81
proteomics_raw %>% count(linked_id) %>% filter(n > 1)
## # A tibble: 3 × 2
## linked_id n
## <chr> <int>
## 1 HC14 2
## 2 P067 2
## 3 P070 2
# already established in 01_Proteomics_Import_QC.Rmd: zero missing values
protein_cols <- proteomics_clean %>% select(-linked_id, -pathology)
cat("Missing values across", ncol(protein_cols), "proteins:", sum(is.na(protein_cols)), "\n")
## Missing values across 1221 proteins: 0
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 = round(100 * n_zero / nrow(metabolomics_raw), 1))
cat("Excluded (>20% below detection):\n")
## Excluded (>20% below detection):
zero_summary %>% filter(pct_zero > 20) %>% arrange(desc(pct_zero))
## # A tibble: 8 × 3
## metabolite n_zero pct_zero
## <chr> <int> <dbl>
## 1 citrate 79 97.5
## 2 phospholipid 77 95.1
## 3 ascorbate 53 65.4
## 4 x2_hydroxybutyrate 39 48.1
## 5 unknown_signal_at_7_14_ppm 38 46.9
## 6 histidine 30 37
## 7 ethanol 25 30.9
## 8 x3_hydroxybutyrate 20 24.7
retained <- zero_summary %>% filter(pct_zero <= 20)
cat("\nRetained:", nrow(retained), "metabolites\n")
##
## Retained: 29 metabolites
cat("Mean % below detection among retained:", round(mean(retained$pct_zero), 1), "\n")
## Mean % below detection among retained: 1.7
Structural feature issue, metabolomics: one
retained-adjacent metabolite, unknown_signal_at_7_14_ppm,
is an NMR peak that could not be assigned to a known metabolite
identity. It also happens to fail the 80% detection threshold (46.9%
below detection) and is excluded on that basis already; it is flagged
here separately as a structural/identity issue distinct from its
missingness.
cont_vars <- c(Age = "age", WCC = "white_cell_count", Haemoglobin = "haemaglobin",
Creatinine = "creatinine_umol_l", CRP = "c_reactive_protein_mg_l",
`Total protein` = "total_protein_g_l", Albumin = "albumin_g_l",
Tbil = "total_bilirubin_tbil_umol_l", Dbil = "conjugated_bilirubin_dbil_umol_l",
ALT = "alanine_transaminase_alt_u_l", AST = "aspartate_transaminase_ast_u_l",
ALP = "alkaline_phosphatase_alp_u_l", GGT = "gamma_glutamyl_transferase_ggt_u_l")
miss_counts <- sapply(cont_vars, function(v) sum(is.na(clinical_clean[[v]])))
miss_pct <- round(100 * miss_counts / nrow(clinical_clean), 1)
data.frame(variable = names(cont_vars), n_missing = miss_counts, pct_missing = miss_pct) %>%
arrange(desc(pct_missing))
## variable n_missing pct_missing
## Total protein Total protein 87 37.3
## Albumin Albumin 78 33.5
## ALT ALT 78 33.5
## AST AST 78 33.5
## Creatinine Creatinine 71 30.5
## WCC WCC 70 30.0
## CRP CRP 70 30.0
## Haemoglobin Haemoglobin 69 29.6
## Tbil Tbil 57 24.5
## Dbil Dbil 57 24.5
## ALP ALP 57 24.5
## GGT GGT 57 24.5
## Age Age 15 6.4
cat("\nOverall mean % missing across panel:", round(mean(miss_pct), 1), "\n")
##
## Overall mean % missing across panel: 27.9
cat("Range:", round(min(miss_pct),1), "-", round(max(miss_pct),1), "\n")
## Range: 6.4 - 37.3
table4 <- tibble::tibble(
`QC characteristic` = c("Initial samples", "Duplicate records/technical replicates",
"Missing identifiers", "Missing clinical data",
"Missing molecular measurements", "Features with structural issues",
"Final samples assessed"),
Clinical = c("233", "0", "0",
paste0(round(mean(miss_pct),1), "% average (range ",
round(min(miss_pct),1), "-", round(max(miss_pct),1), "%)"),
"-", "-", "233"),
Proteomic = c("178", "3 patients (2 runs each, averaged)", "0", "-",
"0/1221 (fully quantified)", "0 identified", "175"),
Metabolomic = c("81", "0", "0", "-",
"8/37 excluded (<80% detection); 29 retained",
"1 (unassigned peak: unknown signal at 7.14 ppm)", "81")
)
knitr::kable(table4, caption = "Table 4. Structural quality and completeness assessment of the datasets")
| QC characteristic | Clinical | Proteomic | Metabolomic |
|---|---|---|---|
| Initial samples | 233 | 178 | 81 |
| Duplicate records/technical replicates | 0 | 3 patients (2 runs each, averaged) | 0 |
| Missing identifiers | 0 | 0 | 0 |
| Missing clinical data | 27.9% average (range 6.4-37.3%) | - | - |
| Missing molecular measurements | - | 0/1221 (fully quantified) | 8/37 excluded (<80% detection); 29 retained |
| Features with structural issues | - | 0 identified | 1 (unassigned peak: unknown signal at 7.14 ppm) |
| Final samples assessed | 233 | 175 | 81 |
cc <- clinical_clean %>% mutate(grp = ifelse(pathology == "HC", "HC (n=42)", "Non-HC (n=191)"))
miss_long <- lapply(names(cont_vars), function(nm) {
v <- cont_vars[[nm]]
cc %>% group_by(grp) %>%
summarise(pct_missing = 100 * mean(is.na(.data[[v]])), n = n()) %>%
mutate(variable = nm)
}) %>% bind_rows()
var_order <- miss_long %>% group_by(variable) %>% summarise(m = mean(pct_missing)) %>%
arrange(desc(m)) %>% pull(variable)
miss_long$variable <- factor(miss_long$variable, levels = var_order)
figure3a <- ggplot(miss_long, aes(x = variable, y = pct_missing, fill = grp)) +
geom_col(position = position_dodge(width = 0.75), width = 0.65) +
scale_fill_manual(values = c("HC (n=42)" = "#B0B0B0", "Non-HC (n=191)" = "#4C72B0"), name = NULL) +
labs(x = NULL, y = "% missing", title = "Figure 3A. Missingness per clinical variable, HC vs. non-HC") +
theme_minimal(base_size = 11) +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
plot.title = element_text(face = "bold", size = 12),
legend.position = "top")
figure3a
ggsave("Figure3A_Missingness_Barchart.png", figure3a, width = 9, height = 5.5, dpi = 200, bg = "white")
Missingness is split by healthy control (HC) status rather than shown per patient, since the pattern is driven almost entirely by this one subgroup: HC patients (n=42) are missing essentially the entire core laboratory panel (at or near 100% for every variable except age), while the rest of the cohort (n=191) shows a modest, fairly consistent missingness rate across variables (5-24%). This reflects a data-collection design difference – healthy controls were not given the full diagnostic laboratory workup administered to symptomatic patients – rather than a data quality problem. This is a “missing not at random” pattern tied to patient subgroup rather than missing-completely-at-random, which should be considered before any imputation strategy is applied to these variables; in particular, imputing HC values from the observed distribution of symptomatic patients would not be appropriate.
Structural quality and completeness were assessed separately for each of the three datasets prior to harmonization (Table 4). No duplicate patient identifiers were found in the clinical or metabolomic datasets; three proteomic patients (HC14, P067, P070) had two technical replicate runs each, which were averaged prior to further analysis. Proteomic data were fully quantified with no missing values across 1221 proteins. Metabolomic data required detection-threshold filtering: 8 of 37 quantified metabolites were excluded for falling below an 80% detection threshold, including one NMR peak that could not be assigned a metabolite identity; the remaining 29 metabolites had a mean 1.7% rate of below-detection values. Clinical data showed the highest overall missingness (mean 27.9% across the core laboratory panel, ranging from 6.4% for age to 37.3% for the most incomplete variable), attributable in part to a structured subgroup of patients missing the full laboratory panel rather than missingness scattered at random across the cohort (Figure 3A).
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] ggplot2_4.0.3 janitor_2.2.1 readxl_1.5.0 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] pkgconfig_2.0.3 pillar_1.11.1 bslib_0.11.0 gtable_0.3.6
## [37] glue_1.8.1 systemfonts_1.3.2 xfun_0.59 tibble_3.3.1
## [41] tidyselect_1.2.1 rstudioapi_0.19.0 knitr_1.51 farver_2.1.2
## [45] htmltools_0.5.9 rmarkdown_2.31 labeling_0.4.3 compiler_4.6.1
## [49] S7_0.2.2