Extended Proteomic and Metabolomic Profiling

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

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; this needs a separate data pull, agreed with supervisors first.
• 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 the proteomic/metabolomic PDAC-vs-BBP comparison called for in 5.5.2, which still 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

Reading this table: any row with p < 0.05 differs significantly across at least one pair of the six groups (the test does not say which pair — follow up with pairwise Wilcoxon tests if a specific comparison, e.g. RPC vs BBP, is needed for the thesis text). Liver-function markers and inflammatory markers are expected to track disease group most closely, since biliary obstruction and systemic inflammation both increase with more advanced disease.

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

5 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

6 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

7 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	p	rho	fdr
serpine1	lipid_alpha_ch2	0	0.6237	0.0013
serpine1	protein_nh	0	-0.5636	0.0166
serpine1	glycorol_phospholipid	0	0.5612	0.0166
serpine1	lipid_ch2	0	0.5469	0.0225
galnt15	glutamine	0	0.5452	0.0225
itih4	lipid_alpha_ch2	0	0.5400	0.0244
gsn	lipid_alpha_ch2	0	-0.5312	0.0318
prdx4	glyc_b	0	0.5281	0.0323
ctsd	protein_nh	0	-0.5200	0.0358
itih4	glycorol_phospholipid	0	0.5191	0.0358
prdx4	glyc_a	0	0.5190	0.0358
acta1	lipid_alpha_ch2	0	-0.5148	0.0398
iglv2_8	glutamine	0	0.5125	0.0406
ywhag	lipid_alpha_ch2	0	0.5089	0.0443
pcbp2	lactate	0	-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")
}

7.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.

8 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")

ann <- data.frame(pathology = clin64$pathology, row.names = clin64$id_std)
ann_colors <- list(pathology = pal)

pheatmap(t(scale(X_all)), show_colnames = FALSE, show_rownames = FALSE,
         annotation_col = ann, annotation_colors = ann_colors,
         main = "Hierarchical clustering: all 64 patients x 1,250 features")

Reading this section: if patients do not visibly separate by pathology colour in either the PCA or the clustering, 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.

9 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%

10 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),
  "extended_profiling_results.rds"
)

11 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