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)

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

dim(clinical_clean)
## [1] 233  54
dim(proteomics_clean)
## [1]  175 1223
dim(metabolomics_clean)
## [1] 81 39
batch_info <- readxl::read_excel("NMR_Metabolomics_Data__main_.xlsx", sheet = "Sheet2") %>%
  janitor::clean_names()

metabolomics_clean <- metabolomics_clean %>%
  left_join(batch_info %>% select(nmr_id, id, batch), by = "nmr_id")

metabolomics_clean %>% select(nmr_id, id, batch, group) %>% slice_head(n = 5)
## # A tibble: 5 × 4
##   nmr_id          id    batch   group
##   <chr>           <chr> <chr>   <chr>
## 1 CRS-20221027-01 P41   Batch 1 RPC  
## 2 CRS-20221027-02 P64   Batch 1 RPC  
## 3 CRS-20221027-03 P34   Batch 1 RPC  
## 4 CRS-20221027-04 P66   Batch 1 MPC  
## 5 CRS-20221027-05 P51   Batch 1 RPC
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)))
}

clinical_std <- clinical_clean %>%
  mutate(id_std = standardize_id(linked_id)) %>%
  select(patient_id, linked_id, id_std, pathology)

proteomics_std <- proteomics_clean %>%
  mutate(id_std = standardize_id(linked_id)) %>%
  select(linked_id, id_std, pathology)

metabolomics_std <- metabolomics_clean %>%
  mutate(id_std = standardize_id(id)) %>%
  select(nmr_id, id, id_std, group)

# quick look at each, side by side
clinical_std %>% slice_head(n = 5)
## # A tibble: 5 × 4
##   patient_id linked_id id_std pathology
##   <chr>      <chr>     <chr>  <chr>    
## 1 RPC1       P001      P001   RPC      
## 2 RPC2       P002      P002   RPC      
## 3 RPC3       P003      P003   RPC      
## 4 RPC4       P004      P004   RPC      
## 5 RPC5       P005      P005   RPC
proteomics_std %>% slice_head(n = 5)
## # A tibble: 5 × 3
##   linked_id id_std pathology
##   <chr>     <chr>  <chr>    
## 1 BP001     BP001  Bmass    
## 2 BP003     BP003  Bmass    
## 3 BP004     BP004  Bmass    
## 4 BP006     BP006  Bmass    
## 5 BP010     BP010  Bmass
metabolomics_std %>% slice_head(n = 5)
## # A tibble: 5 × 4
##   nmr_id          id    id_std group
##   <chr>           <chr> <chr>  <chr>
## 1 CRS-20221027-01 P41   P041   RPC  
## 2 CRS-20221027-02 P64   P064   RPC  
## 3 CRS-20221027-03 P34   P034   RPC  
## 4 CRS-20221027-04 P66   P066   MPC  
## 5 CRS-20221027-05 P51   P051   RPC
clinical_ids <- clinical_std$id_std
proteomics_ids <- proteomics_std$id_std
metabolomics_ids <- metabolomics_std$id_std

overlap_summary <- tibble::tibble(
  comparison = c(
    "Clinical only (total)",
    "Proteomics only (total)",
    "Metabolomics only (total)",
    "Clinical + Proteomics",
    "Clinical + Metabolomics",
    "Proteomics + Metabolomics",
    "All three"
  ),
  n_patients = c(
    length(unique(clinical_ids)),
    length(unique(proteomics_ids)),
    length(unique(metabolomics_ids)),
    length(intersect(clinical_ids, proteomics_ids)),
    length(intersect(clinical_ids, metabolomics_ids)),
    length(intersect(proteomics_ids, metabolomics_ids)),
    length(Reduce(intersect, list(clinical_ids, proteomics_ids, metabolomics_ids)))
  )
)

overlap_summary
## # A tibble: 7 × 2
##   comparison                n_patients
##   <chr>                          <int>
## 1 Clinical only (total)            233
## 2 Proteomics only (total)          175
## 3 Metabolomics only (total)         81
## 4 Clinical + Proteomics            175
## 5 Clinical + Metabolomics           81
## 6 Proteomics + Metabolomics         64
## 7 All three                         64
matched_cohort_ids <- Reduce(intersect, list(clinical_ids, proteomics_ids, metabolomics_ids))

saveRDS(matched_cohort_ids, "matched_cohort_ids.rds")
length(matched_cohort_ids)
## [1] 64
matched_clinical <- clinical_std %>%
  filter(id_std %in% matched_cohort_ids) %>%
  rename(pathology_clinical = pathology)

matched_proteomics <- proteomics_std %>%
  filter(id_std %in% matched_cohort_ids) %>%
  rename(pathology_proteomics = pathology)

matched_metabolomics <- metabolomics_std %>%
  filter(id_std %in% matched_cohort_ids) %>%
  rename(group_metabolomics = group)

# confirm all three pathology/group labels for the same 64 patients, side by side
matched_clinical %>%
  select(id_std, pathology_clinical) %>%
  left_join(matched_proteomics %>% select(id_std, pathology_proteomics), by = "id_std") %>%
  left_join(matched_metabolomics %>% select(id_std, group_metabolomics), by = "id_std") %>%
  count(pathology_clinical, pathology_proteomics, group_metabolomics)
## # A tibble: 5 × 4
##   pathology_clinical pathology_proteomics group_metabolomics     n
##   <chr>              <chr>                <chr>              <int>
## 1 LAPC               PDAC                 LAPC                  10
## 2 MPC                PDAC                 MPC                    4
## 3 RPC                PDAC                 LAPC                   2
## 4 RPC                PDAC                 MPC                    5
## 5 RPC                PDAC                 RPC                   43
mismatches <- matched_clinical %>%
  select(id_std, pathology_clinical) %>%
  left_join(matched_proteomics %>% select(id_std, pathology_proteomics), by = "id_std") %>%
  left_join(matched_metabolomics %>% select(id_std, group_metabolomics), by = "id_std") %>%
  filter(pathology_clinical != group_metabolomics)

mismatches
## # A tibble: 7 × 4
##   id_std pathology_clinical pathology_proteomics group_metabolomics
##   <chr>  <chr>              <chr>                <chr>             
## 1 P011   RPC                PDAC                 MPC               
## 2 P019   RPC                PDAC                 LAPC              
## 3 P036   RPC                PDAC                 MPC               
## 4 P046   RPC                PDAC                 MPC               
## 5 P074   RPC                PDAC                 LAPC              
## 6 P075   RPC                PDAC                 MPC               
## 7 P081   RPC                PDAC                 MPC
saveRDS(mismatches, "pathology_mismatches.rds")
write.csv(mismatches, "pathology_mismatches.csv", row.names = FALSE)

Update following re-verification (see Appendix): re-checking this step against the source spreadsheets identified three additional patients (P035, P064, P076) whose clinical pathology value was inconsistent with metabolomics in a direction that PDAC cannot follow biologically – metastatic-to-resectable or locally-advanced-to-resectable – since disease stage cannot regress. These three could not be explained by progression and were confirmed as clinical data-entry errors; the clinical pathology values for P035, P064, and P076 have been corrected to match metabolomics.

Seven of 64 matched patients still show disagreement between clinical pathology and metabolomics group classification (all agree with proteomics, which uses a single “PDAC” category for all cancer patients): P011, P019, P036, P046, P074, P075, and P081. All seven follow the direction consistent with disease progression (RPC to LAPC and/or MPC), which PDAC can plausibly follow between clinical staging and metabolomics sample collection. These are retained as documented, unresolved discordances rather than corrected outright, since directional plausibility alone does not confirm progression actually occurred – confirming this would require the clinical staging date and metabolomics sample collection date for each patient, which have not yet been cross-checked. Flagged to the supervisor for confirmation before these patients are used in any pathology-stratified analysis.

Table 3: Dataset structure and participant overlap

table3 <- tibble::tibble(
  Dataset = c("Clinical", "Proteomic", "Metabolomic", "Proteomic + Metabolomic"),
  `Initial participants` = c(
    length(unique(clinical_ids)),
    length(unique(proteomics_ids)),
    length(unique(metabolomics_ids)),
    NA
  ),
  `Matched participants (all 3 sources)` = c(
    length(matched_cohort_ids),
    length(matched_cohort_ids),
    length(matched_cohort_ids),
    length(matched_cohort_ids)
  ),
  `Key features` = c(
    "Demographic, clinical, and pathology variables",
    paste0(ncol(proteomics_clean) - 2, " proteins (DIA-MS)"),
    paste0(ncol(metabolomics_clean %>% select(-nmr_id, -id, -batch, -group)),
           " metabolites (1H-NMR); 29 retained after 80% detection filtering"),
    "Matched multi-omics cohort"
  )
)

knitr::kable(table3, caption = "Table 3. Structure and participant overlap across the clinical, proteomic and metabolomic datasets")
Table 3. Structure and participant overlap across the clinical, proteomic and metabolomic datasets
Dataset Initial participants Matched participants (all 3 sources) Key features
Clinical 233 64 Demographic, clinical, and pathology variables
Proteomic 175 64 1221 proteins (DIA-MS)
Metabolomic 81 64 37 metabolites (1H-NMR); 29 retained after 80% detection filtering
Proteomic + Metabolomic NA 64 Matched multi-omics cohort

This table establishes the technical basis for the 64-patient matched multi-omics cohort used throughout the remainder of Objective 1 and as the input population for Objective 2. It is distinct from the clinical characterization in Tables 1 and 2 (see 06_Clinical_Characterization.Rmd), which describe the full recruited clinical cohort (185 and 102 patients respectively) rather than the subset with complete multi-omics data – readers should note that Table 1’s 102 PDAC patients are not all represented in the 64-patient harmonized dataset; only those with matched proteomic and metabolomic data are.

Figure 2: Data integration workflow

The diagram below traces the same steps computed above – ID harmonization, duplicate checking, sample correspondence, and clinical-omics linkage – laid out as a workflow, with the actual participant counts and processing decisions carried through from this dataset rather than a generic template.

library(ggplot2)
library(grid)

boxes <- data.frame(
  id = 1:9,
  x = c(6, 1.8, 6, 10.2, 6, 6, 6, 6, 6),
  y = c(11.3, 9.2, 9.2, 9.2, 7.3, 6.1, 4.9, 3.7, 2.1),
  w = c(5.0, 3.0, 3.0, 3.0, 8.6, 8.6, 8.6, 8.6, 4.8),
  h = c(0.9, 1.1, 1.1, 1.1, 1.0, 1.0, 1.0, 1.0, 0.95),
  label = c(
    paste0("Parent PDAC study\n(n = ", length(unique(clinical_ids)), " clinical records)"),
    paste0("Clinical dataset\n(n = ", length(unique(clinical_ids)), ")"),
    paste0("Proteomics: DIA-MS\n(n = ", length(unique(proteomics_ids)), ")"),
    paste0("Metabolomics: \u00b9H-NMR\n(n = ", length(unique(metabolomics_ids)), ")"),
    "Patient ID harmonization\nunpadded / zero-padded / prefix -> single standardized format",
    "Duplicate checking\n3 proteomics patients with replicate runs identified & averaged",
    "Sample correspondence\nID overlap checked across all three sources",
    "Clinical-omics linkage\nintersection of clinical, proteomics & metabolomics IDs",
    paste0("Integrated PDAC dataset\n(n = ", length(matched_cohort_ids), " matched multi-omics cohort)")
  ),
  fill = c("#2C3E50","#4C72B0","#DD8452","#55A868","#8172B2","#8172B2","#8172B2","#8172B2","#C44E52")
)

arrows <- data.frame(
  x = c(6, 1.8, 6, 10.2, 6, 6, 6, 6),
  y = c(10.85,8.65,8.65,8.65, 6.8, 5.6, 4.4, 3.2),
  xend=c(6, 6, 6, 6, 6, 6, 6, 6),
  yend=c(9.75,7.85,7.85,7.85, 6.6, 5.4, 4.2, 2.575)
)

figure2 <- ggplot() +
  geom_rect(data=boxes, aes(xmin=x-w/2, xmax=x+w/2, ymin=y-h/2, ymax=y+h/2, fill=fill),
            color="grey20", linewidth=0.4) +
  scale_fill_identity() +
  geom_text(data=boxes, aes(x=x, y=y, label=label), color="white", size=3.0, lineheight=1.0, fontface="bold") +
  geom_segment(data=arrows, aes(x=x,y=y,xend=xend,yend=yend),
               arrow=arrow(length=unit(0.12,"inches"), type="closed"), linewidth=0.5, color="grey30") +
  coord_cartesian(xlim=c(-1,13), ylim=c(1.4,11.9)) +
  theme_void() +
  labs(title="Figure 2. Data integration workflow: clinical, proteomic, and metabolomic datasets") +
  theme(plot.title=element_text(size=11, face="bold", hjust=0.5, margin=margin(b=10)))

figure2

ggsave("Figure2_Data_Integration_Workflow.png", figure2, width=10.5, height=10.5, dpi=200, bg="white")