1 Overview

This file covers the sections that were not yet analysed, using the full 64-patient harmonized cohort (all PDAC, no controls — a BBP comparison is a separate data question, not covered here). It follows the reporting style used in the group’s own AP and GBC papers: an ordinal Spearman correlation against disease severity (here, pathology group), plus a Kruskal–Wallis test across the three groups, both reported with p-value and FDR side by side, as in those papers’ tables.

Pathology groups are ordered by clinical severity: RPC (resectable) = 1, LAPC (locally advanced) = 2, MPC (metastatic) = 3. This mirrors the AP paper’s HC/MAP/MSAP/SAP ordinal coding.

Not covered here, and why:

PDAC vs BBP — no BBP samples exist in the 64-patient matched omics cohort; a related clinical-level comparison is reported separately in the thesis (Section 5.5.2), using the full 233-patient clinical file.
Pathway-level integration — needs new tooling (clusterProfiler / ReactomePA, protein-to-gene mapping) not yet set up; scoped separately.
KODAMA — the unsupervised method used in the group’s AP/GBC papers is a specific package not used elsewhere in this project. Standard PCA and hierarchical clustering are used instead below, which show the same thing (whether molecular phenotype tracks pathology group). Confirm with supervisors whether KODAMA itself is expected for consistency with prior group papers.

2 Load data

harm      <- readRDS("harmonized_multiomics_dataset.rds")
prot      <- readRDS("harmonized_proteomics.rds")
metab     <- readRDS("harmonized_metabolomics.rds")
clin      <- readRDS("clinical_clean.rds")

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)))
}

# All 64 patients, with pathology, age and gender
clin64 <- clin %>%
  mutate(id_std = standardize_id(linked_id)) %>%
  filter(id_std %in% harm$id_std) %>%
  distinct(id_std, .keep_all = TRUE) %>%
  mutate(stage_ord = case_when(
    pathology == "RPC"  ~ 1L,
    pathology == "LAPC" ~ 2L,
    pathology == "MPC"  ~ 3L,
    TRUE ~ NA_integer_
  )) %>%
  select(id_std, pathology, stage_ord, age, gender)

stopifnot(nrow(clin64) == 64, sum(is.na(clin64$stage_ord)) == 0)
table(clin64$pathology)

## 
## LAPC  MPC  RPC 
##   10    4   50

prot_mat  <- prot  %>% column_to_rownames("id_std") %>% as.matrix()
metab_mat <- metab %>% column_to_rownames("id_std") %>% as.matrix()

# align row order to clin64
prot_mat  <- prot_mat[clin64$id_std, ]
metab_mat <- metab_mat[clin64$id_std, ]

dim(prot_mat); dim(metab_mat)

## [1]   64 1221

## [1] 64 29

3 Helper functions

# Per-feature Kruskal-Wallis across 3 groups + Spearman correlation with
# ordinal severity, both with BH-adjusted FDR. Mirrors the AP paper's Table 2.
group_association_table <- function(X, group3, stage_ord) {
  kw <- apply(X, 2, function(v) {
    tryCatch(kruskal.test(v, group3)$p.value, error = function(e) NA_real_)
  })
  sp <- apply(X, 2, function(v) {
    tryCatch({
      r <- suppressWarnings(cor.test(v, stage_ord, method = "spearman"))
      c(rho = unname(r$estimate), p = r$p.value)
    }, error = function(e) c(rho = NA_real_, p = NA_real_))
  })
  tibble(
    feature   = colnames(X),
    kw_p      = kw,
    kw_fdr    = p.adjust(kw, method = "BH"),
    spearman_rho = sp["rho", ],
    spearman_p   = sp["p", ],
    spearman_fdr = p.adjust(sp["p", ], method = "BH")
  ) %>% arrange(spearman_p)
}

# Per-feature association with a continuous clinical variable (age)
clinical_cor_table <- function(X, clinvar) {
  res <- apply(X, 2, function(v) {
    tryCatch({
      r <- suppressWarnings(cor.test(v, clinvar, method = "spearman"))
      c(rho = unname(r$estimate), p = r$p.value)
    }, error = function(e) c(rho = NA_real_, p = NA_real_))
  })
  tibble(feature = colnames(X), rho = res["rho", ], p = res["p", ]) %>%
    mutate(fdr = p.adjust(p, method = "BH")) %>%
    arrange(p)
}

# Per-feature association with a two-level clinical variable (gender)
clinical_wilcox_table <- function(X, clinvar) {
  p <- apply(X, 2, function(v) {
    tryCatch(wilcox.test(v ~ clinvar)$p.value, error = function(e) NA_real_)
  })
  tibble(feature = colnames(X), p = p) %>%
    mutate(fdr = p.adjust(p, method = "BH")) %>%
    arrange(p)
}

4 Clinical parameters across all diagnosis groups (Table 1)

This uses the full 233-patient clinical file, not the 64-patient omics cohort, because clinical_clean.rds contains three groups with no omics data attached: BBP (benign biliary pathology), CP (chronic pancreatitis) and HC (healthy control). This answers a clinical baseline question — how routine lab parameters differ across all six diagnostic groups — and is a different, narrower thing than a proteomic/metabolomic PDAC-vs-BBP comparison, which cannot be done because BBP patients were never run through DIA-MS or NMR.

group_order <- c("HC", "BBP", "CP", "RPC", "LAPC", "MPC")

clin_all <- clin %>%
  filter(pathology %in% group_order) %>%
  mutate(pathology = factor(pathology, levels = group_order))

table(clin_all$pathology)

## 
##   HC  BBP   CP  RPC LAPC  MPC 
##   42   83    6   79   13   10

lab_vars <- c(
  white_cell_count                   = "WCC",
  haemaglobin                        = "Haemoglobin",
  creatinine_umol_l                  = "Creatinine",
  c_reactive_protein_mg_l            = "CRP",
  total_protein_g_l                  = "Total Protein",
  albumin_g_l                        = "Albumin",
  total_bilirubin_tbil_umol_l        = "Tbil",
  conjugated_bilirubin_dbil_umol_l   = "Dbil",
  alanine_transaminase_alt_u_l       = "ALT",
  aspartate_transaminase_ast_u_l     = "AST",
  alkaline_phosphatase_alp_u_l       = "ALP",
  gamma_glutamyl_transferase_ggt_u_l = "GGT"
)

stopifnot(all(names(lab_vars) %in% names(clin_all)))

median_iqr <- function(x) {
  x <- x[!is.na(x)]
  if (length(x) == 0) return("-")
  sprintf("%.1f [%.1f\u2013%.1f]", median(x), quantile(x, 0.25), quantile(x, 0.75))
}

table1_row <- function(var, label) {
  x <- clin_all[[var]]
  n_total <- sum(!is.na(x))
  by_grp <- sapply(group_order, function(g) median_iqr(x[clin_all$pathology == g]))
  kw_p <- tryCatch(kruskal.test(x ~ clin_all$pathology)$p.value, error = function(e) NA_real_)
  as_tibble_row(c(list(Feature = label, N = n_total), as.list(by_grp),
                  list(p = kw_p)))
}

table1 <- bind_rows(lapply(names(lab_vars), function(v) table1_row(v, lab_vars[[v]])))
table1 <- table1 %>% mutate(p = ifelse(p < 0.01, "<0.01", sprintf("%.2f", p)))

knitr::kable(table1)

Feature	N	HC	BBP	CP	RPC	LAPC	MPC	p
WCC	163	-	8.1 [6.5–9.8]	6.8 [6.2–8.2]	8.9 [7.0–13.2]	10.8 [7.1–11.4]	13.4 [9.1–15.5]	0.09
Haemoglobin	164	-	11.9 [10.7–13.7]	12.7 [9.6–14.6]	10.4 [9.2–12.2]	9.9 [7.9–11.6]	11.1 [9.7–12.2]	<0.01
Creatinine	162	-	69.5 [61.0–82.8]	90.0 [64.2–130.0]	72.0 [58.0–97.2]	63.0 [56.0–82.0]	71.0 [57.0–78.0]	0.78
CRP	163	-	20.0 [8.0–73.0]	70.0 [20.0–162.8]	56.0 [19.0–135.0]	94.0 [39.0–136.0]	113.0 [50.8–222.8]	<0.01
Total Protein	146	-	72.5 [67.2–78.8]	66.0 [61.2–71.5]	64.0 [57.5–68.0]	65.0 [55.5–67.5]	72.0 [69.0–74.2]	<0.01
Albumin	155	-	39.5 [35.2–42.8]	36.5 [35.2–38.5]	31.0 [27.0–35.0]	27.0 [25.0–30.0]	36.0 [31.0–38.0]	<0.01
Tbil	176	-	30.0 [10.0–92.0]	5.0 [5.0–8.8]	160.0 [65.0–268.0]	100.0 [43.0–185.0]	60.0 [40.0–335.0]	<0.01
Dbil	176	-	23.0 [4.0–73.5]	2.0 [2.0–2.8]	141.0 [37.0–246.0]	90.0 [31.0–151.0]	51.0 [38.0–315.0]	<0.01
ALT	155	-	47.5 [17.0–107.5]	18.0 [15.0–22.5]	92.0 [41.0–151.0]	64.0 [17.0–101.0]	48.0 [30.0–362.0]	0.01
AST	155	-	53.0 [28.0–90.0]	28.5 [22.5–34.5]	104.0 [59.0–162.0]	106.0 [51.0–118.0]	123.0 [57.0–268.0]	<0.01
ALP	176	-	200.0 [111.0–403.5]	74.0 [70.5–104.5]	615.0 [302.0–1045.0]	386.0 [237.0–789.0]	397.0 [302.0–1645.0]	<0.01
GGT	176	-	171.0 [80.5–500.5]	61.5 [53.5–246.5]	502.0 [198.0–1155.0]	313.0 [151.0–683.0]	839.0 [400.0–1021.0]	<0.01

4.1 Pairwise follow-up: which groups differ

Kruskal-Wallis only shows that some pair of the six groups differs, not which one. This runs pairwise Wilcoxon tests for every parameter above, focusing on two questions: does each PDAC stage differ from BBP (benign disease), and do the three PDAC stages differ from each other. FDR correction is applied within each parameter, across its full set of pairwise comparisons.

pairwise_for_var <- function(var, label) {
  x <- clin_all[[var]]
  g <- clin_all$pathology
  keep <- !is.na(x) & !is.na(g)
  pw <- pairwise.wilcox.test(x[keep], droplevels(g[keep]),
                              p.adjust.method = "BH")
  as.data.frame(as.table(pw$p.value)) %>%
    filter(!is.na(Freq)) %>%
    transmute(Feature = label,
              group1 = as.character(Var1), group2 = as.character(Var2),
              p_fdr = Freq) %>%
    arrange(p_fdr)
}

pairwise_all <- bind_rows(lapply(names(lab_vars), function(v) {
  pairwise_for_var(v, lab_vars[[v]])
}))

# Focus 1: each PDAC stage vs BBP
bbp_comparisons <- pairwise_all %>%
  filter((group1 == "BBP" & group2 %in% c("RPC", "LAPC", "MPC")) |
         (group2 == "BBP" & group1 %in% c("RPC", "LAPC", "MPC"))) %>%
  filter(p_fdr < 0.05) %>%
  arrange(Feature, p_fdr)

cat("Significant PDAC-stage vs BBP comparisons (FDR < 0.05):", nrow(bbp_comparisons), "\n")

## Significant PDAC-stage vs BBP comparisons (FDR < 0.05): 20

knitr::kable(bbp_comparisons, digits = 4)

Feature	group1	group2	p_fdr
ALP	RPC	BBP	0.0000
ALP	MPC	BBP	0.0337
AST	RPC	BBP	0.0027
Albumin	RPC	BBP	0.0000
Albumin	LAPC	BBP	0.0003
CRP	RPC	BBP	0.0291
CRP	LAPC	BBP	0.0434
CRP	MPC	BBP	0.0434
Dbil	RPC	BBP	0.0000
Dbil	LAPC	BBP	0.0246
Dbil	MPC	BBP	0.0299
GGT	RPC	BBP	0.0002
GGT	MPC	BBP	0.0381
Haemoglobin	RPC	BBP	0.0016
Haemoglobin	LAPC	BBP	0.0055
Tbil	RPC	BBP	0.0000
Tbil	LAPC	BBP	0.0413
Tbil	MPC	BBP	0.0413
Total Protein	RPC	BBP	0.0000
Total Protein	LAPC	BBP	0.0130

# Focus 2: the three PDAC stages against each other
stage_comparisons <- pairwise_all %>%
  filter(group1 %in% c("RPC", "LAPC", "MPC") & group2 %in% c("RPC", "LAPC", "MPC")) %>%
  filter(p_fdr < 0.05) %>%
  arrange(Feature, p_fdr)

cat("Significant PDAC-stage vs PDAC-stage comparisons (FDR < 0.05):", nrow(stage_comparisons), "\n")

## Significant PDAC-stage vs PDAC-stage comparisons (FDR < 0.05): 2

knitr::kable(stage_comparisons, digits = 4)

Feature	group1	group2	p_fdr
Total Protein	MPC	RPC	0.0130
Total Protein	MPC	LAPC	0.0131

5 5.5.3 Proteomic differences across LAPC, MPC and RPC

prot_stage <- group_association_table(prot_mat, clin64$pathology, clin64$stage_ord)
sig_prot_stage <- prot_stage %>% filter(kw_fdr < 0.10 | spearman_fdr < 0.10)
nrow(sig_prot_stage)

## [1] 0

knitr::kable(head(sig_prot_stage, 15), digits = 4)

feature	kw_p	kw_fdr	spearman_rho	spearman_p	spearman_fdr

if (nrow(sig_prot_stage) == 0) {
  cat("**No proteins passed FDR < 0.10 across pathology groups in this cohort.**",
      "This is a legitimate result to report as-is: with 50 RPC, 10 LAPC and",
      "4 MPC patients, the MPC group is small and statistical power for a",
      "3-group comparison is limited. Report the top nominal (unadjusted",
      "p < 0.05) associations descriptively, with this caveat stated",
      "explicitly, rather than implying a validated finding.\n")
}

No proteins passed FDR < 0.10 across pathology groups in this cohort. This is a legitimate result to report as-is: with 50 RPC, 10 LAPC and 4 MPC patients, the MPC group is small and statistical power for a 3-group comparison is limited. Report the top nominal (unadjusted p < 0.05) associations descriptively, with this caveat stated explicitly, rather than implying a validated finding.

prot_stage %>% filter(kw_p < 0.05 | spearman_p < 0.05) %>%
  head(15) %>% knitr::kable(digits = 4)

feature	kw_p	kw_fdr	spearman_rho	spearman_p	spearman_fdr
hspa1a	0.0082	0.999	0.3746	0.0023	0.9964
igf2	0.0230	0.999	0.3335	0.0071	0.9964
sptb	0.0332	0.999	-0.3219	0.0095	0.9964
ppp6c	0.0365	0.999	-0.3216	0.0096	0.9964
galnt15	0.0595	0.999	0.2948	0.0181	0.9964
prr4	0.0669	0.999	0.2926	0.0190	0.9964
kpnb1	0.0734	0.999	-0.2877	0.0212	0.9964
fam114a2	0.0669	0.999	-0.2854	0.0223	0.9964
igfbp1	0.0846	0.999	0.2796	0.0253	0.9964
septin7	0.0678	0.999	0.2747	0.0280	0.9964
ppp6r3	0.0983	0.999	-0.2711	0.0303	0.9964
hspa4l	0.0922	0.999	0.2669	0.0330	0.9964
ahnak	0.0967	0.999	0.2645	0.0347	0.9964
coro1b	0.0360	0.999	0.2641	0.0350	0.9964
psmd2	0.1151	0.999	-0.2609	0.0373	0.9964

6 5.6.2 Metabolomic differences across LAPC, MPC and RPC

metab_stage <- group_association_table(metab_mat, clin64$pathology, clin64$stage_ord)
sig_metab_stage <- metab_stage %>% filter(kw_fdr < 0.10 | spearman_fdr < 0.10)
nrow(sig_metab_stage)

## [1] 0

knitr::kable(metab_stage %>% head(15), digits = 4)

feature	kw_p	kw_fdr	spearman_rho	spearman_p	spearman_fdr
glucose	0.1839	0.8614	0.2279	0.0702	0.8213
mannose	0.2218	0.8614	0.2159	0.0866	0.8213
pyruvate	0.2519	0.8614	0.2090	0.0974	0.8213
glycine	0.2937	0.8614	0.1932	0.1261	0.8213
isoleucine	0.3564	0.8614	0.1805	0.1534	0.8213
lactate	0.3431	0.8614	0.1721	0.1739	0.8213
acetoacetate	0.2707	0.8614	0.1345	0.2895	0.8213
leucine	0.2272	0.8614	0.1343	0.2899	0.8213
glyc_b	0.2615	0.8614	0.1273	0.3162	0.8213
cholesterol	0.5688	0.9104	0.1241	0.3286	0.8213
protein_nh	0.3515	0.8614	-0.1165	0.3591	0.8213
lipid_beta_ch2	0.6088	0.9104	0.1115	0.3806	0.8213
glyc_a	0.5843	0.9104	0.1018	0.4235	0.8213
glutamine	0.7227	0.9104	0.0979	0.4417	0.8213
tyrosine	0.5970	0.9104	0.0923	0.4684	0.8213

7 5.5.4 / 5.6.3 Clinical and pathological associations

Age (continuous, Spearman) and gender (two-level, Wilcoxon) are tested against every protein and metabolite. This uses what the clinical file actually has consistently for all 64 patients; if other variables (e.g. smoking, comorbidities) should be included, they can be added the same way.

prot_age  <- clinical_cor_table(prot_mat, clin64$age)
metab_age <- clinical_cor_table(metab_mat, clin64$age)

cat("Proteins with age FDR < 0.10:", sum(prot_age$fdr < 0.10, na.rm = TRUE), "\n")

## Proteins with age FDR < 0.10: 0

cat("Metabolites with age FDR < 0.10:", sum(metab_age$fdr < 0.10, na.rm = TRUE), "\n")

## Metabolites with age FDR < 0.10: 0

knitr::kable(head(prot_age, 10), digits = 4)

feature	rho	p	fdr
acly	0.3816	0.0024	0.8952
ighv5_51	0.3810	0.0025	0.8952
hspb1	-0.3777	0.0027	0.8952
iglv3_19	0.3707	0.0033	0.8952
tagln2	-0.3638	0.0039	0.8952
iglv3_25	0.3556	0.0049	0.8952
ppia	-0.3540	0.0051	0.8952
ilk	-0.3373	0.0079	0.9379
iglv2_18	0.3359	0.0081	0.9379
lrg1	0.3305	0.0093	0.9379

knitr::kable(head(metab_age, 10), digits = 4)

feature	rho	p	fdr
creatinine	0.3441	0.0066	0.1920
glycorol_phospholipid	0.2433	0.0588	0.7913
glucose	-0.2246	0.0819	0.7913
phenylalanine	0.1983	0.1255	0.8064
lipid_alpha_ch2	0.1916	0.1390	0.8064
creatine	-0.1721	0.1847	0.8870
glycine	-0.1481	0.2547	0.8870
glutamate	-0.1367	0.2934	0.8870
leucine	-0.1254	0.3357	0.8870
glyc_b	-0.1220	0.3489	0.8870

gender_f <- factor(clin64$gender)
table(gender_f)

## gender_f
## female   male 
##     20     41

prot_gender  <- clinical_wilcox_table(prot_mat, gender_f)
metab_gender <- clinical_wilcox_table(metab_mat, gender_f)

cat("Proteins with gender FDR < 0.10:", sum(prot_gender$fdr < 0.10, na.rm = TRUE), "\n")

## Proteins with gender FDR < 0.10: 0

cat("Metabolites with gender FDR < 0.10:", sum(metab_gender$fdr < 0.10, na.rm = TRUE), "\n")

## Metabolites with gender FDR < 0.10: 0

knitr::kable(head(prot_gender, 10), digits = 4)

feature	p	fdr
mapk14	0.0013	0.9859
psmd6	0.0029	0.9859
ighv1_18	0.0046	0.9859
iglv3_25	0.0054	0.9859
pcdh12	0.0060	0.9859
hpse	0.0073	0.9859
iglv3_9	0.0076	0.9859
aco2	0.0122	0.9859
eif3b	0.0122	0.9859
ppm1f	0.0122	0.9859

knitr::kable(head(metab_gender, 10), digits = 4)

feature	p	fdr
creatine	0.0411	0.8302
lipid_ch_ch2_ch	0.0661	0.8302
unsaturated_lipid_ch_ch	0.0989	0.8302
phenylalanine	0.1233	0.8302
formate	0.1804	0.8302
tyrosine	0.2678	0.8302
protein_nh	0.3097	0.8302
lipid_ch3	0.3097	0.8302
valine	0.3322	0.8302
glucose	0.3477	0.8302

8 5.7.2 Protein-metabolite associations

Every protein is correlated (Spearman) against every metabolite: 1,221 × 29 = 35,409 pairs. FDR is applied across the full set of pairwise tests.

cor_mat <- cor(prot_mat, metab_mat, method = "spearman", use = "pairwise.complete.obs")
dim(cor_mat)

## [1] 1221   29

# p-values via cor.test, pairwise (slower but exact)
n <- nrow(prot_mat)
t_stat <- cor_mat * sqrt((n - 2) / (1 - cor_mat^2))
p_mat  <- 2 * pt(-abs(t_stat), df = n - 2)

pm_long <- as_tibble(as.table(p_mat), .name_repair = "minimal")
colnames(pm_long) <- c("protein", "metabolite", "p")
rho_long <- as_tibble(as.table(cor_mat), .name_repair = "minimal")
colnames(rho_long) <- c("protein", "metabolite", "rho")

pm_results <- pm_long %>%
  left_join(rho_long, by = c("protein", "metabolite")) %>%
  mutate(fdr = p.adjust(p, method = "BH")) %>%
  arrange(p)

cat("Pairs with FDR < 0.10:", sum(pm_results$fdr < 0.10, na.rm = TRUE), "\n")

## Pairs with FDR < 0.10: 232

knitr::kable(head(pm_results, 15), digits = 4)

protein	metabolite	rho	fdr
serpine1	lipid_alpha_ch2	0.6237	0.0013
serpine1	protein_nh	-0.5636	0.0166
serpine1	glycorol_phospholipid	0.5612	0.0166
serpine1	lipid_ch2	0.5469	0.0225
galnt15	glutamine	0.5452	0.0225
itih4	lipid_alpha_ch2	0.5400	0.0244
gsn	lipid_alpha_ch2	-0.5312	0.0318
prdx4	glyc_b	0.5281	0.0323
ctsd	protein_nh	-0.5200	0.0358
itih4	glycorol_phospholipid	0.5191	0.0358
prdx4	glyc_a	0.5190	0.0358
acta1	lipid_alpha_ch2	-0.5148	0.0398
iglv2_8	glutamine	0.5125	0.0406
ywhag	lipid_alpha_ch2	0.5089	0.0443
pcbp2	lactate	-0.5036	0.0471

top_metab_by_hits <- pm_results %>% filter(fdr < 0.10) %>%
  count(metabolite, sort = TRUE) %>% pull(metabolite)

if (length(top_metab_by_hits) >= 2) {
  top_prot_by_hits <- pm_results %>% filter(fdr < 0.10) %>%
    count(protein, sort = TRUE) %>% slice_head(n = 40) %>% pull(protein)
  pheatmap(cor_mat[top_prot_by_hits, top_metab_by_hits, drop = FALSE],
           main = "Spearman rho: top proteins x metabolites with FDR < 0.10",
           fontsize_row = 6, fontsize_col = 8)
} else {
  cat("Fewer than 2 metabolites had any FDR < 0.10 hit; heatmap skipped.",
      "The strongest nominal pairs are shown in the table above instead.\n")
}

8.1 SERPINE1: a closer look

SERPINE1 appears in 4 of the top 15 protein-metabolite pairs. This checks whether it is a genuinely broad hub across metabolites, or driven by one or two outlying patients.

serpine1_cor <- pm_results %>% filter(protein == "serpine1") %>% arrange(p)
cat("SERPINE1 pairs with FDR < 0.10:", sum(serpine1_cor$fdr < 0.10, na.rm = TRUE),
    "of 29 metabolites tested\n")

## SERPINE1 pairs with FDR < 0.10: 4 of 29 metabolites tested

knitr::kable(serpine1_cor, digits = 4)

protein	metabolite	p	rho	fdr
serpine1	lipid_alpha_ch2	0.0000	0.6237	0.0013
serpine1	protein_nh	0.0000	-0.5636	0.0166
serpine1	glycorol_phospholipid	0.0000	0.5612	0.0166
serpine1	lipid_ch2	0.0000	0.5469	0.0225
serpine1	unsaturated_lipid_ch_ch	0.0192	0.2921	0.2515
serpine1	creatinine	0.0201	0.2899	0.2542
serpine1	alanine	0.0300	-0.2715	0.2867
serpine1	glycine	0.0476	-0.2486	0.3366
serpine1	cholesterol	0.0590	-0.2374	0.3649
serpine1	glutamate	0.1149	-0.1990	0.4680
serpine1	mannose	0.1772	0.1708	0.5475
serpine1	glutamine	0.2247	0.1539	0.5994
serpine1	glucose	0.2275	-0.1530	0.6024
serpine1	lipid_beta_ch2	0.2375	0.1498	0.6110
serpine1	pyruvate	0.2437	0.1478	0.6167
serpine1	phenylalanine	0.2604	0.1428	0.6316
serpine1	acetate	0.2929	-0.1335	0.6589
serpine1	lactate	0.4364	0.0990	0.7622
serpine1	acetoacetate	0.5293	0.0801	0.8156
serpine1	isoleucine	0.6062	0.0657	0.8553
serpine1	lipid_ch_ch2_ch	0.6424	0.0592	0.8715
serpine1	formate	0.7068	-0.0479	0.8987
serpine1	lipid_ch3	0.7127	0.0469	0.9007
serpine1	glyc_b	0.7714	-0.0370	0.9250
serpine1	leucine	0.8235	-0.0284	0.9435
serpine1	creatine	0.8395	-0.0258	0.9483
serpine1	glyc_a	0.9230	-0.0123	0.9762
serpine1	valine	0.9305	-0.0111	0.9783
serpine1	tyrosine	0.9573	0.0068	0.9878

pal <- c(RPC = "#1D9E75", LAPC = "#EF9F27", MPC = "#D85A30")
top_pair <- serpine1_cor$metabolite[1]
plot(prot_mat[, "serpine1"], metab_mat[, top_pair],
     xlab = "SERPINE1 (standardized)", ylab = paste0(top_pair, " (standardized)"),
     pch = 19, col = pal[clin64$pathology],
     main = sprintf("SERPINE1 vs %s (rho = %.2f, FDR = %.4f)",
                     top_pair, serpine1_cor$rho[1], serpine1_cor$fdr[1]))
legend("topleft", legend = names(pal), col = pal, pch = 19, bty = "n")

Reading this: if the strongest pair’s scatter looks like a smooth trend across most patients, that supports a real biological association. If it looks driven by one or two extreme points sitting apart from the rest, the correlation is fragile and should be reported as such, the same way SPTAN1’s influence from P009 was checked in Objective 2.

9 5.7.3 Integrated molecular patterns

A joint PCA on the combined 1,250-feature harmonized matrix (all 64 patients), coloured by pathology group. This is the standard-methods equivalent of the papers’ KODAMA unsupervised analysis: it asks the same question — does the molecular profile separate by pathology group at all?

X_all <- harm %>% column_to_rownames("id_std") %>% as.matrix()
X_all <- X_all[clin64$id_std, ]

pca <- prcomp(X_all, scale. = FALSE)  # data already scaled in Objective 1
var_exp <- summary(pca)$importance[2, 1:2] * 100
var_exp

##    PC1    PC2 
## 34.178  9.439

pal <- c(RPC = "#1D9E75", LAPC = "#EF9F27", MPC = "#D85A30")
plot(pca$x[, 1], pca$x[, 2], col = pal[clin64$pathology], pch = 19,
     xlab = sprintf("PC1 (%.1f%%)", var_exp[1]),
     ylab = sprintf("PC2 (%.1f%%)", var_exp[2]),
     main = "Joint proteomic-metabolomic PCA, coloured by pathology")
legend("topright", legend = names(pal), col = pal, pch = 19, bty = "n")

Reading this section: if patients do not visibly separate by pathology colour on the PCA, that is itself the honest result — it says the dominant variation in the omics profile is not primarily driven by resectability stage, which is worth stating plainly rather than searching for a cluster structure that isn’t there.

10 Results summary

summary_tbl <- tibble(
  Section = c("5.5.3 Proteins x stage", "5.6.2 Metabolites x stage",
              "5.5.4 Proteins x age", "5.5.4 Proteins x gender",
              "5.6.3 Metabolites x age", "5.6.3 Metabolites x gender",
              "5.7.2 Protein-metabolite pairs", "5.7.3 Joint PCA (PC1+PC2 variance)"),
  Result = c(
    sprintf("%d of %d proteins FDR<0.10 (KW or Spearman)", nrow(sig_prot_stage), ncol(prot_mat)),
    sprintf("%d of %d metabolites FDR<0.10", nrow(sig_metab_stage), ncol(metab_mat)),
    sprintf("%d of %d FDR<0.10", sum(prot_age$fdr < 0.10, na.rm = TRUE), ncol(prot_mat)),
    sprintf("%d of %d FDR<0.10", sum(prot_gender$fdr < 0.10, na.rm = TRUE), ncol(prot_mat)),
    sprintf("%d of %d FDR<0.10", sum(metab_age$fdr < 0.10, na.rm = TRUE), ncol(metab_mat)),
    sprintf("%d of %d FDR<0.10", sum(metab_gender$fdr < 0.10, na.rm = TRUE), ncol(metab_mat)),
    sprintf("%d of %d pairs FDR<0.10", sum(pm_results$fdr < 0.10, na.rm = TRUE), nrow(pm_results)),
    sprintf("%.1f%% + %.1f%%", var_exp[1], var_exp[2])
  )
)
knitr::kable(summary_tbl)

Section	Result
5.5.3 Proteins x stage	0 of 1221 proteins FDR<0.10 (KW or Spearman)
5.6.2 Metabolites x stage	0 of 29 metabolites FDR<0.10
5.5.4 Proteins x age	0 of 1221 FDR<0.10
5.5.4 Proteins x gender	0 of 1221 FDR<0.10
5.6.3 Metabolites x age	0 of 29 FDR<0.10
5.6.3 Metabolites x gender	0 of 29 FDR<0.10
5.7.2 Protein-metabolite pairs	232 of 35409 pairs FDR<0.10
5.7.3 Joint PCA (PC1+PC2 variance)	34.2% + 9.4%

11 Save results

saveRDS(
  list(prot_stage = prot_stage, metab_stage = metab_stage,
       prot_age = prot_age, prot_gender = prot_gender,
       metab_age = metab_age, metab_gender = metab_gender,
       pm_results = pm_results, pca = pca, clin64 = clin64,
       table1 = table1, bbp_comparisons = bbp_comparisons,
       stage_comparisons = stage_comparisons),
  "extended_profiling_results.rds"
)

12 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] pheatmap_1.0.13 tibble_3.3.1    stringr_1.6.0   dplyr_1.2.1    
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.7.3        cli_3.6.6          knitr_1.51         rlang_1.2.0       
##  [5] xfun_0.59          stringi_1.8.7      otel_0.2.0         generics_0.1.4    
##  [9] jsonlite_2.0.0     glue_1.8.1         htmltools_0.5.9    sass_0.4.10       
## [13] scales_1.4.0       rmarkdown_2.31     grid_4.6.1         evaluate_1.0.5    
## [17] jquerylib_0.1.4    fastmap_1.2.0      yaml_2.3.12        lifecycle_1.0.5   
## [21] compiler_4.6.1     codetools_0.2-20   RColorBrewer_1.1-3 pkgconfig_2.0.3   
## [25] rstudioapi_0.19.0  farver_2.1.2       digest_0.6.39      R6_2.6.1          
## [29] tidyselect_1.2.1   pillar_1.11.1      magrittr_2.0.5     bslib_0.11.0      
## [33] withr_3.0.3        gtable_0.3.6       tools_4.6.1        cachem_1.1.0

Extended Proteomic and Metabolomic Profiling

Sections 5.5.3–5.5.4, 5.6.2–5.6.3, 5.7.2–5.7.3

Judy Ngobeni

2026-09-23