library(readxl)
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(janitor)
## 
## Attaching package: 'janitor'
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
metabolomics_raw <- read_excel("NMR_Metabolomics_Data__main_.xlsx", sheet = "Sheet1") %>%
  clean_names()

dim(metabolomics_raw)
## [1] 81 39
metabolomics_raw <- metabolomics_raw %>%
  mutate(across(-c(nmr_id, group), as.numeric))

# sanity check: did every value convert cleanly, with no accidental NAs introduced?
sum(is.na(metabolomics_raw %>% select(-nmr_id, -group)))
## [1] 0
metabolomics_raw %>% tabyl(group)
##  group  n   percent
##   LAPC 15 0.1851852
##    MPC  9 0.1111111
##    RPC 57 0.7037037
metabolomics_raw %>%
  select(-nmr_id, -group) %>%
  summarise(across(everything(), ~ sum(is.na(.)))) %>%
  rowSums()
## [1] 0
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)) %>%
  arrange(desc(n_zero))
## # A tibble: 37 × 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
##  9 leucine                        12     14.8
## 10 creatine                        8      9.9
## # ℹ 27 more rows
metabolomics_raw %>%
  count(nmr_id) %>%
  filter(n > 1)
## # A tibble: 0 × 2
## # ℹ 2 variables: nmr_id <chr>, n <int>
metabolite_values <- metabolomics_raw %>% select(-nmr_id, -group)

range(metabolite_values, na.rm = TRUE)
## [1]   0.0000 278.3874
batch_info <- read_excel("NMR_Metabolomics_Data__main_.xlsx", sheet = "Sheet2") %>%
  clean_names()

metabolomics_raw %>%
  left_join(batch_info, by = "nmr_id") %>%
  select(nmr_id, group = group.x, id, batch, group_check = group.y) %>%
  slice_head(n = 5)
## # A tibble: 5 × 5
##   nmr_id          group id    batch   group_check
##   <chr>           <chr> <chr> <chr>   <chr>      
## 1 CRS-20221027-01 RPC   P41   Batch 1 RPC        
## 2 CRS-20221027-02 RPC   P64   Batch 1 RPC        
## 3 CRS-20221027-03 RPC   P34   Batch 1 RPC        
## 4 CRS-20221027-04 MPC   P66   Batch 1 MPC        
## 5 CRS-20221027-05 RPC   P51   Batch 1 RPC

Two items require supervisor input before this data can be used in modelling or joined via the crosswalk:

  1. Citrate and phospholipid are recorded as zero in 97.5% and 95.1% of patients respectively, with only 2-4 non-zero readings each across 81 patients. Please confirm whether this reflects expected below-detection- limit behaviour for this sample type, or a measurement/processing issue — this determines whether these two metabolites should be excluded from modelling or retained with this limitation noted.

  2. This metabolomics cohort uses only three diagnosis groups (RPC, LAPC, MPC — cancer subtypes only, no healthy controls or benign groups), which differs from both the clinical data (six groups: BBP, RPC, HC, LAPC, MPC, CP) and the proteomics data (three groups, but Bmass/HC/PDAC — a different three). All three classification schemes will need to be reconciled before the crosswalk join.

saveRDS(metabolomics_raw, "metabolomics_clean.rds")
list.files(pattern = ".rds")
## [1] "clinical_clean.rds"     "metabolomics_clean.rds" "proteomics_clean.rds"
metabolomics_raw %>%
  left_join(batch_info, by = "nmr_id") %>%
  filter(group.x != group.y) %>%
  select(nmr_id, group.x, group.y)
## # A tibble: 0 × 3
## # ℹ 3 variables: nmr_id <chr>, group.x <chr>, group.y <chr>
saveRDS(metabolomics_raw, "metabolomics_clean.rds")
list.files(pattern = ".rds")
## [1] "clinical_clean.rds"     "metabolomics_clean.rds" "proteomics_clean.rds"