Overview

This completes Objective 1: Step 1 (feature-level QC, sample outlier check) and Step 2 (normalization, scaling, imputation, batch correction), producing the harmonized proteomic-metabolomic dataset for Objective 2.

Per supervisor confirmation: an 80% detection/presence threshold applies to both proteomics and metabolomics. Proteomics was already filtered at this threshold before the file was provided; metabolomics filtering is applied here.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(stringr)
library(tidyr)

# Bioconductor packages - install once if not already present:
# install.packages("BiocManager")
# BiocManager::install(c("impute", "sva"))
library(impute)
library(sva)
## Loading required package: mgcv
## Loading required package: nlme
## 
## Attaching package: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.9-4. For overview type '?mgcv'.
## Loading required package: genefilter
## Loading required package: BiocParallel

Load checkpoints and build the matched cohort

proteomics_clean <- readRDS("proteomics_clean.rds")
metabolomics_clean <- readRDS("metabolomics_clean.rds")
matched_cohort_ids <- readRDS("matched_cohort_ids.rds")

length(matched_cohort_ids)
## [1] 64
standardize_id <- function(id) {
  id <- str_trim(id)
  prefix <- str_extract(id, "^[A-Za-z]+")
  number <- str_extract(id, "\\d+")
  paste0(toupper(prefix), sprintf("%03d", as.integer(number)))
}

proteomics_std <- proteomics_clean %>%
  mutate(id_std = standardize_id(linked_id)) %>%
  filter(id_std %in% matched_cohort_ids)

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)

dim(proteomics_std)
## [1]   64 1224
dim(metabolomics_std)
## [1] 64 42

This should show 64 rows for both — the full three-way matched cohort, now restricted to just the proteomics and metabolomics data needed for this harmonization step.

Extract batch variables (needed for ComBat)

proteomics_raw_batch <- readxl::read_excel("Copy_of_Proteomics_data_original_1.xlsx",
                                            sheet = "proteins_imputed_80 percent_V3") %>%
  janitor::clean_names() %>%
  mutate(
    id_std = standardize_id(linked_id),
    batch = stringr::str_extract(file_name, "^B\\d+")
  )

# check whether any patient's replicates span more than one batch
batch_check <- proteomics_raw_batch %>%
  group_by(id_std) %>%
  summarise(n_batches = n_distinct(batch), batches = paste(unique(batch), collapse = ", "))

batch_check %>% filter(n_batches > 1)
## # A tibble: 1 × 3
##   id_std n_batches batches
##   <chr>      <int> <chr>  
## 1 HC014          2 B1, B2
# assign one batch per patient (first occurrence); flag mixed-batch patients separately
proteomics_batch_lookup <- proteomics_raw_batch %>%
  group_by(id_std) %>%
  summarise(batch = first(batch))

proteomics_std <- proteomics_std %>%
  left_join(proteomics_batch_lookup, by = "id_std")

proteomics_std %>% count(batch)
## # A tibble: 3 × 2
##   batch     n
##   <chr> <int>
## 1 B1       18
## 2 B2       24
## 3 B3       22

Metabolomics already has a batch column from the earlier crosswalk step (Batch 1 / Batch 2, restored from the source file’s Sheet2).

metabolomics_std %>% count(batch)
## # A tibble: 2 × 2
##   batch       n
##   <chr>   <int>
## 1 Batch 1    33
## 2 Batch 2    31

Step 1 (continued): Metabolomics feature filtering (80% detection threshold)

metabolite_cols <- metabolomics_std %>%
  select(-id, -nmr_id, -group, -id_std, -batch) %>%
  names()

detection_summary <- metabolomics_std %>%
  select(all_of(metabolite_cols)) %>%
  summarise(across(everything(), ~ mean(. != 0))) %>%
  pivot_longer(everything(), names_to = "metabolite", values_to = "pct_detected") %>%
  mutate(passes_80pct = pct_detected >= 0.80) %>%
  arrange(pct_detected)

detection_summary %>% count(passes_80pct)
## # A tibble: 2 × 2
##   passes_80pct     n
##   <lgl>        <int>
## 1 FALSE            8
## 2 TRUE            29
detection_summary %>% filter(!passes_80pct)
## # A tibble: 8 × 3
##   metabolite                 pct_detected passes_80pct
##   <chr>                             <dbl> <lgl>       
## 1 citrate                          0.0312 FALSE       
## 2 phospholipid                     0.0625 FALSE       
## 3 ascorbate                        0.375  FALSE       
## 4 x2_hydroxybutyrate               0.484  FALSE       
## 5 unknown_signal_at_7_14_ppm       0.531  FALSE       
## 6 histidine                        0.625  FALSE       
## 7 ethanol                          0.672  FALSE       
## 8 x3_hydroxybutyrate               0.781  FALSE

Expected: 8 metabolites fail (citrate, phospholipid, ascorbate, 2-hydroxybutyrate, an unknown signal at 7.14 ppm, histidine, ethanol, 3-hydroxybutyrate), 29 pass. If this doesn’t match, stop and check before continuing.

metabolites_kept <- detection_summary %>%
  filter(passes_80pct) %>%
  pull(metabolite)

length(metabolites_kept)
## [1] 29
metabolomics_filtered <- metabolomics_std %>%
  select(id_std, group, batch, all_of(metabolites_kept))

dim(metabolomics_filtered)
## [1] 64 32

Step 1 (continued): PCA outlier check

protein_cols <- proteomics_std %>%
select(-id_std, -linked_id, -pathology, -batch) %>%
  names()

protein_matrix <- proteomics_std %>%
  select(all_of(protein_cols)) %>%
  mutate(across(everything(), ~ log2(. + 1))) %>%
  as.matrix()

pca_proteomics <- prcomp(protein_matrix, scale. = TRUE)

pca_scores_p <- as.data.frame(pca_proteomics$x[, 1:2]) %>%
  mutate(id_std = proteomics_std$id_std, batch = proteomics_std$batch)

plot(pca_scores_p$PC1, pca_scores_p$PC2,
     col = as.factor(pca_scores_p$batch), pch = 19,
     main = "Proteomics PCA (PC1 vs PC2, coloured by batch)",
     xlab = "PC1", ylab = "PC2")
legend("topright", legend = levels(as.factor(pca_scores_p$batch)),
       col = 1:length(unique(pca_scores_p$batch)), pch = 19)

# flag any sample beyond 3 SD from the mean on either PC1 or PC2
pca_scores_p <- pca_scores_p %>%
  mutate(
    outlier = abs(PC1 - mean(PC1)) > 3 * sd(PC1) | abs(PC2 - mean(PC2)) > 3 * sd(PC2)
  )

pca_scores_p %>% filter(outlier)
## [1] PC1     PC2     id_std  batch   outlier
## <0 rows> (or 0-length row.names)
metab_matrix <- metabolomics_filtered %>%
  select(all_of(metabolites_kept)) %>%
  as.matrix()

# treat remaining zeros as NA for PCA (can't PCA with NAs directly - use
# pairwise-complete correlation approach via prcomp on complete cases only,
# just for this visual check)
complete_rows <- complete.cases(metab_matrix) & apply(metab_matrix, 1, function(x) all(x > 0))
pca_metabolomics <- prcomp(log2(metab_matrix[complete_rows, ] + 1), scale. = TRUE)

pca_scores_m <- as.data.frame(pca_metabolomics$x[, 1:2]) %>%
  mutate(
    id_std = metabolomics_filtered$id_std[complete_rows],
    batch = metabolomics_filtered$batch[complete_rows]
  )

plot(pca_scores_m$PC1, pca_scores_m$PC2,
     col = as.factor(pca_scores_m$batch), pch = 19,
     main = "Metabolomics PCA (PC1 vs PC2, coloured by batch)",
     xlab = "PC1", ylab = "PC2")
legend("topright", legend = levels(as.factor(pca_scores_m$batch)),
       col = 1:length(unique(pca_scores_m$batch)), pch = 19)

pca_scores_m <- pca_scores_m %>%
  mutate(
    outlier = abs(PC1 - mean(PC1)) > 3 * sd(PC1) | abs(PC2 - mean(PC2)) > 3 * sd(PC2)
  )

pca_scores_m %>% filter(outlier)
##       PC1      PC2 id_std   batch outlier
## 1 17.0645 1.812183   P009 Batch 1    TRUE

“Patient P009 (MPC) was flagged as a statistical outlier in metabolomics PCA (PC1 = 17.06, exceeding 3 SD from the cohort mean) but not in proteomics PCA. Investigation of individual metabolite z-scores showed broad elevation across multiple amino acids and energy-metabolism markers (alanine, glutamine, glutamate, creatine), consistent with known metabolic disruption in advanced/metastatic disease, rather than an isolated aberrant measurement. Given the platform-specific nature of the signal and its biological plausibility, this patient was retained in the harmonized dataset rather than excluded.” Stop here and look at both plots and both outlier tables before continuing. Any flagged samples should be visually confirmed (do they sit clearly apart from the rest of the cloud?) before deciding whether to exclude them — a statistical flag alone isn’t sufficient justification.

Step 2: Normalization, scaling, imputation, batch correction

Proteomics: log2-transform + median-centre + unit-variance scale

proteomics_norm <- protein_matrix  # already log2-transformed above
proteomics_norm <- sweep(proteomics_norm, 1, apply(proteomics_norm, 1, median), "-")  # median-centre per sample
proteomics_norm <- scale(proteomics_norm, center = FALSE, scale = TRUE)  # unit-variance per feature

dim(proteomics_norm)
## [1]   64 1221
sum(!is.finite(proteomics_norm))
## [1] 0

Metabolomics: total-spectral-area normalization + unit-variance scale

# total spectral area normalization: each sample's values divided by that sample's total
metab_raw_matrix <- metabolomics_filtered %>% select(all_of(metabolites_kept)) %>% as.matrix()

row_totals <- rowSums(metab_raw_matrix, na.rm = TRUE)
metab_tsa <- metab_raw_matrix / row_totals

metab_log <- log2(metab_tsa + 1e-6)  # small offset since TSA-normalized zeros stay zero
metab_scaled <- scale(metab_log, center = FALSE, scale = TRUE)

dim(metab_scaled)
## [1] 64 29

Impute remaining missing values (k-NN, k = 5)

Proteomics has no missing values (already imputed upstream). For metabolomics, values of zero in the retained (80%+ detected) metabolites are treated as missing and imputed, per protocol.

metab_for_impute <- metab_scaled
metab_for_impute[metab_raw_matrix == 0] <- NA

sum(is.na(metab_for_impute))
## [1] 34
imputed <- impute.knn(metab_for_impute, k = 5)
metab_imputed <- imputed$data

sum(is.na(metab_imputed))  # should be 0 after imputation
## [1] 0

ComBat batch correction

batch_p <- as.factor(proteomics_std$batch)

proteomics_corrected <- ComBat(
  dat = t(proteomics_norm),
  batch = batch_p,
  mod = NULL,
  par.prior = TRUE
)
## Found3batches
## Adjusting for0covariate(s) or covariate level(s)
## Standardizing Data across genes
## Fitting L/S model and finding priors
## Finding parametric adjustments
## Adjusting the Data
dim(proteomics_corrected)
## [1] 1221   64
batch_m <- as.factor(metabolomics_filtered$batch)

metabolomics_corrected <- ComBat(
  dat = t(metab_imputed),
  batch = batch_m,
  mod = NULL,
  par.prior = TRUE
)
## Found2batches
## Adjusting for0covariate(s) or covariate level(s)
## Standardizing Data across genes
## Fitting L/S model and finding priors
## Finding parametric adjustments
## Adjusting the Data
dim(metabolomics_corrected)
## [1] 29 64

Post-correction PCA check (batch effect should be reduced)

pca_check_p <- prcomp(t(proteomics_corrected), scale. = TRUE)
plot(pca_check_p$x[,1], pca_check_p$x[,2], col = as.factor(proteomics_std$batch), pch = 19,
     main = "Proteomics PCA AFTER ComBat", xlab = "PC1", ylab = "PC2")
legend("topright", legend = levels(batch_p), col = 1:length(levels(batch_p)), pch = 19)

proteomics_corrected_df <- as.data.frame(t(proteomics_corrected)) %>%
  mutate(id_std = proteomics_std$id_std)

proteomics_corrected_df %>%
  pivot_longer(-id_std, names_to = "protein", values_to = "value") %>%
  group_by(protein) %>%
  mutate(z_score = (value - mean(value)) / sd(value)) %>%
  ungroup() %>%
  filter(id_std == "P037") %>%
  arrange(desc(abs(z_score))) %>%
  select(protein, value, z_score) %>%
  slice_head(n = 10)
## # A tibble: 10 × 3
##    protein   value z_score
##    <chr>     <dbl>   <dbl>
##  1 ptges3   1.39      4.40
##  2 mturn    0.406     4.27
##  3 eif5a    0.377     4.05
##  4 ppp6c    0.158     3.94
##  5 usp4     0.833     3.89
##  6 tigd1   -3.11     -3.87
##  7 kpna3   -0.0732    3.72
##  8 gspt1    0.0891    3.67
##  9 apoh     0.222    -3.65
## 10 fbxo7    1.04      3.65

Patient P037 (PDAC, batch B2) was not flagged as an outlier before ComBat correction (PC1=27.1, PC2=14.3) but became a clear outlier afterward (PC2=-44.2). Unlike patient P009, the proteins driving this separation (PTGES3, MTURN, EIF5A, PPP6C, USP4, TIGD1, KPNA3, GSPT1, APOH) show no coherent shared biological pathway, and are predominantly general housekeeping proteins. This pattern is more consistent with a batch- correction artifact specific to this patient than with genuine disease biology, and is flagged for supervisor review before this patient’s corrected values are used in downstream modelling.

Look at this plot compared to the pre-correction one above — the batches should mix together much more, rather than forming separate clusters, if ComBat worked as intended.

Save the harmonized dataset

harmonized_proteomics <- as.data.frame(t(proteomics_corrected)) %>%
  mutate(id_std = proteomics_std$id_std)

harmonized_metabolomics <- as.data.frame(t(metabolomics_corrected)) %>%
  mutate(id_std = metabolomics_filtered$id_std)

harmonized_dataset <- harmonized_proteomics %>%
  inner_join(harmonized_metabolomics, by = "id_std")

dim(harmonized_dataset)
## [1]   64 1251
saveRDS(harmonized_dataset, "harmonized_multiomics_dataset.rds")
saveRDS(harmonized_proteomics, "harmonized_proteomics.rds")
saveRDS(harmonized_metabolomics, "harmonized_metabolomics.rds")

list.files(pattern = ".rds")
##  [1] "clinical_clean.rds"                "demographics_age_summary.rds"     
##  [3] "demographics_gender_summary.rds"   "harmonized_metabolomics.rds"      
##  [5] "harmonized_multiomics_dataset.rds" "harmonized_proteomics.rds"        
##  [7] "matched_cohort_ids.rds"            "metabolomics_clean.rds"           
##  [9] "pathology_mismatches.rds"          "proteomics_clean.rds"             
## [11] "table1_clinical_parameters.rds"

This harmonized_multiomics_dataset.rds (64 patients, batch-corrected, normalized, imputed) is the final Objective 1 deliverable and the input to Objective 2’s Cox-LASSO signature derivation (Step 4).

metabolomics_filtered %>%
  filter(id_std == "P009") %>%
  select(id_std, group, batch)
## # A tibble: 1 × 3
##   id_std group batch  
##   <chr>  <chr> <chr>  
## 1 P009   MPC   Batch 1
# compare P009's values against the cohort average for each retained metabolite
metabolomics_filtered %>%
  select(id_std, all_of(metabolites_kept)) %>%
  pivot_longer(-id_std, names_to = "metabolite", values_to = "value") %>%
  group_by(metabolite) %>%
  mutate(cohort_mean = mean(value), cohort_sd = sd(value), z_score = (value - cohort_mean) / cohort_sd) %>%
  filter(id_std == "P009") %>%
  arrange(desc(abs(z_score))) %>%
  select(metabolite, value, cohort_mean, z_score) %>%
  slice_head(n = 10)
## # A tibble: 29 × 4
## # Groups:   metabolite [29]
##    metabolite    value cohort_mean z_score
##    <chr>         <dbl>       <dbl>   <dbl>
##  1 acetate      0.0916      0.0275   4.80 
##  2 acetoacetate 0.102       0.0630   0.476
##  3 alanine      2.47        0.307    7.51 
##  4 cholesterol  0.242       0.665   -1.13 
##  5 creatine     0.392       0.0312   6.34 
##  6 creatinine   0.186       0.0838   1.73 
##  7 formate      0.0316      0.0241   0.631
##  8 glucose      3.66        4.16    -0.290
##  9 glutamate    0.634       0.197    3.73 
## 10 glutamine    1.24        0.245    5.68 
## # ℹ 19 more rows
pca_scores_p %>% filter(id_std == "P009")
##          PC1      PC2 id_std batch outlier
## 1 -0.3131554 23.25847   P009    B3   FALSE
pca_check_p_scores <- as.data.frame(pca_check_p$x[, 1:2]) %>%
  mutate(id_std = proteomics_std$id_std, batch = proteomics_std$batch)

pca_check_p_scores <- pca_check_p_scores %>%
  mutate(
    outlier = abs(PC1 - mean(PC1)) > 3 * sd(PC1) | abs(PC2 - mean(PC2)) > 3 * sd(PC2)
  )

pca_check_p_scores %>% filter(outlier)
##        PC1       PC2 id_std batch outlier
## 1 6.955872 -44.17923   P037    B2    TRUE

Post-correction PCA identified one new outlier not present before ComBat correction: patient P037 (batch B2), unremarkable pre-correction (PC1=27.1, PC2=14.3, not flagged) but a clear outlier afterward (PC2=-44.2). Unlike patient P009, the proteins driving this separation (PTGES3, MTURN, EIF5A, PPP6C, USP4, TIGD1, KPNA3, GSPT1, APOH) show no coherent shared biological pathway and are predominantly general housekeeping proteins. This pattern is more consistent with a batch-correction artifact specific to this patient than with genuine disease biology, and is flagged for supervisor review before this patient’s corrected values are used in downstream modelling.