01_AddIDs.R

# ==============================================================================
#' This script reads the data from the mzmine output and transfers the ID of
#' mzmine to method specific ID for the respective project. The number from 
#' mzmine is prefixed with a method prefix and numbers are filled by zeros. This
#' is done for MS1 and MS2 data.
#' 
#' Version: 1.0
#' Created by: PD Dr. Michael Witting
#' Date: 22.01.2025

# ==============================================================================
# read positive mode data
# ==============================================================================
# read MS1 features and format ID ----------------------------------------------
# read csv exported from mzmine
pos_data <- read_csv(paste0(project_path,
                            "/04_mzmine_output/pos_data.csv"))

# select only relevant columns
pos_data_selected <- pos_data %>% 
  select(c("id", "rt", "mz"),
         starts_with("datafile") & ends_with("area")) %>% 
  mutate(id = paste0(toupper(method_prefix),
                     "POSFT",
                     str_pad(id, width = 5, side = "left", pad = "0")))

# rename columns to remove flags
colnames(pos_data_selected) <- pos_data_selected %>% 
  colnames() %>% 
  str_replace("^datafile:", "") %>% 
  str_replace(":area$", "")

# get sample names and create study design from it -----------------------------
pos_sample_names <- pos_data_selected %>% 
  select(-c("id", "rt", "mz")) %>% 
  colnames()

pos_studydesign <- data.frame(samples = pos_sample_names,
                              batch = as.integer(str_extract(str_extract(pos_sample_names, "Btc\\d+"), "\\d+")),
                              injection = as.integer(str_extract(str_extract(pos_sample_names, "_\\d+_"), "\\d+")),
                              type = str_extract(pos_sample_names, "CMTRX|MTRX|blank|Blank|Warmup|RTI|PRCS|STRD|QC"))

pos_studydesign$type[which(is.na(pos_studydesign$type))] <- "sample"
pos_studydesign$type[which(str_detect(pos_studydesign$samples, "_MS2.mzML$"))] <- "MS2"

# normalize type for non-standard/non-NT nomenclature in samples
pos_studydesign$type[which(pos_studydesign$type %in% c("blank", "Blank"))] <- "PRCS"
pos_studydesign$type[which(pos_studydesign$type %in% c("QC"))] <- "CMTRX"

#get sample name
pos_studydesign$sample <- pos_studydesign$samples
pos_studydesign$sample[which(!pos_studydesign$type == "sample")] <- NA_character_
pos_studydesign$sample <- str_replace(pos_studydesign$sample, "^Btc\\d+_\\d+_", "")
pos_studydesign$sample <- str_replace(pos_studydesign$sample, ".mzML$", "")

# add sample group column
pos_studydesign$group <- NA_character_

# order for data output
pos_studydesign <- pos_studydesign %>% 
  arrange(match(type, c("Warmup", "blank", "Blank", "PRCS", "RTI", "STRD", "MTRX", "CMTRX")),
          batch,
          injection)

pos_data_selected <- pos_data_selected %>% 
  relocate(any_of(c("id", "rt", "mz", pos_studydesign$samples)))

# filter Warmup samples away
pos_data_selected <- pos_data_selected %>% 
  select(-which(str_detect(colnames(.),"Warmup")))

# remove Warmup studies
pos_studydesign <- pos_studydesign %>% 
  filter(type != "Warmup") 

# add potential sample weights -------------------------------------------------
if(file.exists(paste0(project_path,
                      "/00_project_view/sampleWeights.tsv"))) {
  
  # read sample weights and add to table
  sampleWeights <- read_tsv(paste0(project_path,
                                   "/00_project_view/sampleWeights.tsv"),
                            col_types = "cn")
  
  pos_studydesign <- full_join(pos_studydesign,
                               sampleWeights,
                               by = join_by(sample == sample))
  
  pos_studydesign$sampleWeight_quant[which(is.na(pos_studydesign$sampleWeight_quant))] <- sampleWeight_target
  pos_studydesign$sampleWeight_norm <- 1
  
} else {
  
  pos_studydesign <- pos_studydesign %>% 
    mutate(sampleWeight_quant = sampleWeight_target,
           sampleWeight_norm = 1)
  
}

# read annotations and clean ---------------------------------------------------
pos_anno <- read_csv(paste0(project_path,
                            "/04_mzmine_output/pos_annotations.csv")) %>% 
  mutate(id = paste0(toupper(method_prefix),
                     "POSFT",
                     str_pad(id, width = 5, side = "left", pad = "0")))
  
# filter based on annotation source and library
pos_anno_ms2_inhouse <- pos_anno %>% 
  filter(method == "spectral_library_annotation") %>% 
  filter(str_detect(compound_name, "_MPC$")) %>% 
  mutate(compound_name = str_replace(compound_name, "_MPC$", ""))

pos_anno_ms2_external <- pos_anno %>% 
  filter(method == "spectral_library_annotation") %>% 
  filter(!str_detect(compound_name, "_MPC$"))

pos_anno_ms2_lipid <- pos_anno %>% 
  filter(method == "matchedlipid")

pos_anno_ms1_inhouse <- pos_anno %>% 
  filter(method == "simple_compound_db_annotation") %>% 
  mutate(compound_name = str_replace(compound_name, "_MPC$", ""))

# TODO add aggregation of results....

# read positive mode MS2 data --------------------------------------------------
pos_spec <- Spectra(file = paste0(project_path,
                                  "/04_mzmine_output/pos_sirius.mgf"),
                    source = MsBackendMgf(),
                    backend = MsBackendDataFrame(),
                    mapping = c(rtime = "RTINSECONDS",
                                acquisitionNum = "SCANS",
                                precursorMz = "PEPMASS",
                                precursorIntensity = "PEPMASSINT",
                                precursorCharge = "CHARGE",
                                msLevel = "MSLEVEL",
                                id = "FEATURE_ID"))

pos_spec$id <- paste0(toupper(method_prefix),
                      "POSFT",
                      str_pad(pos_spec$id, width = 5, side = "left", pad = "0"))

pos_spec$title <- pos_spec$id

# write data to new files ------------------------------------------------------
# MS1 data
pos_data_selected %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/pos_data_id.tsv"))
pos_studydesign %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/pos_studydesign_id.tsv"))
pos_anno_ms2_inhouse %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/pos_anno_ms2_inhouse_id.tsv"))
pos_anno_ms2_external %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/pos_anno_ms2_external_id.tsv"))
pos_anno_ms2_lipid %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/pos_anno_ms2_lipid_id.tsv"))
pos_anno_ms1_inhouse %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/pos_anno_ms1_inhouse_id.tsv"))

# MS2 data
export(pos_spec,
       backend = MsBackendMgf(),
       file = paste0(project_path,
                     "/05_results/01_ids/pos_sirius_id.mgf"),
       mapping = c(rtime = "RTINSECONDS",
                   acquisitionNum = "SCANS",
                   precursorMz = "PEPMASS",
                   precursorIntensity = "PEPMASSINT",
                   precursorCharge = "CHARGE",
                   msLevel = "MSLEVEL",
                   id = "FEATURE_ID",
                   title = "TITLE"))

# ==============================================================================
# read negative mode data
# ==============================================================================
# read MS1 features and format ID ----------------------------------------------
# read csv exported from mzmine
neg_data <- read_csv(paste0(project_path,
                            "/04_mzmine_output/neg_data.csv"))

# select only relevant columns
neg_data_selected <- neg_data %>% 
  select(c("id", "rt", "mz"),
         starts_with("datafile") & ends_with("area")) %>% 
  mutate(id = paste0(toupper(method_prefix),
                     "NEGFT",
                     str_pad(id, width = 5, side = "left", pad = "0")))

# rename columns to remove flags
colnames(neg_data_selected) <- neg_data_selected %>% 
  colnames() %>% 
  str_replace("^datafile:", "") %>% 
  str_replace(":area$", "")

# get sample names and create study design from it -----------------------------
neg_sample_names <- neg_data_selected %>% 
  select(-c("id", "rt", "mz")) %>% 
  colnames()

neg_studydesign <- data.frame(samples = neg_sample_names,
                              batch = as.integer(str_extract(str_extract(neg_sample_names, "Btc\\d+"), "\\d+")),
                              injection = as.integer(str_extract(str_extract(neg_sample_names, "_\\d+_"), "\\d+")),
                              type = str_extract(neg_sample_names, "CMTRX|MTRX|blank|Blank|Warmup|RTI|PRCS|STRD|QC"))

neg_studydesign$type[which(is.na(neg_studydesign$type))] <- "sample"
neg_studydesign$type[which(str_detect(neg_studydesign$samples, "_MS2.mzML$"))] <- "MS2"

# normalize type for non-standard/non-NT nomenclature in samples
neg_studydesign$type[which(neg_studydesign$type %in% c("blank", "Blank"))] <- "PRCS"
neg_studydesign$type[which(neg_studydesign$type %in% c("QC"))] <- "CMTRX"

#get sample name
neg_studydesign$sample <- neg_studydesign$samples
neg_studydesign$sample[which(!neg_studydesign$type == "sample")] <- NA_character_
neg_studydesign$sample <- str_replace(neg_studydesign$sample, "^Btc\\d+_\\d+_", "")
neg_studydesign$sample <- str_replace(neg_studydesign$sample, ".mzML$", "")

# add sample group column
neg_studydesign$group <- NA_character_

# order for data output
neg_studydesign <- neg_studydesign %>% 
  arrange(match(type, c("Warmup", "blank", "Blank", "PRCS", "RTI", "STRD", "MTRX", "CMTRX")),
          batch,
          injection)

neg_data_selected <- neg_data_selected %>% 
  relocate(any_of(c("id", "rt", "mz", neg_studydesign$samples)))

# filter Warmup samples away
neg_data_selected <- neg_data_selected %>% 
  select(-which(str_detect(colnames(.),"Warmup")))

# remove Warmup studies
neg_studydesign <- neg_studydesign %>% 
  filter(type != "Warmup")

# add potential sample weights -------------------------------------------------
if(file.exists(paste0(project_path,
                      "/00_project_view/sampleWeights.tsv"))) {
  
  # read sample weights and add to table
  sampleWeights <- read_tsv(paste0(project_path,
                                   "/00_project_view/sampleWeights.tsv"),
                            col_types = "cn")
  
  neg_studydesign <- full_join(neg_studydesign,
                               sampleWeights,
                               by = join_by(sample == sample))
  
  neg_studydesign$sampleWeight_quant[which(is.na(neg_studydesign$sampleWeight_quant))] <- sampleWeight_target
  neg_studydesign$sampleWeight_norm <- 1
  
} else {
  
  neg_studydesign <- neg_studydesign %>% 
    mutate(sampleWeight_quant = sampleWeight_target,
           sampleWeight_norm = 1)
  
}

# read annotations and clean ---------------------------------------------------
neg_anno <- read_csv(paste0(project_path,
                            "/04_mzmine_output/neg_annotations.csv")) %>% 
mutate(id = paste0(toupper(method_prefix),
                   "NEGFT",
                   str_pad(id, width = 5, side = "left", pad = "0")))

# filter based on annotation source and library
neg_anno_ms2_inhouse <- neg_anno %>% 
  filter(method == "spectral_library_annotation") %>% 
  filter(str_detect(compound_name, "_MPC$")) %>% 
  mutate(compound_name = str_replace(compound_name, "_MPC$", ""))

neg_anno_ms2_external <- neg_anno %>% 
  filter(method == "spectral_library_annotation") %>% 
  filter(!str_detect(compound_name, "_MPC$"))

neg_anno_ms2_lipid <- neg_anno %>% 
  filter(method == "matchedlipid")

neg_anno_ms1_inhouse <- neg_anno %>% 
  filter(method == "simple_compound_db_annotation") %>% 
  mutate(compound_name = str_replace(compound_name, "_MPC$", ""))

# read negative mode MS2 data --------------------------------------------------
neg_spec <- Spectra(file = paste0(project_path,
                                  "/04_mzmine_output/neg_sirius.mgf"),
                    source = MsBackendMgf(),
                    backend = MsBackendDataFrame(),
                    mapping = c(rtime = "RTINSECONDS",
                                acquisitionNum = "SCANS",
                                precursorMz = "PEPMASS",
                                precursorIntensity = "PEPMASSINT",
                                precursorCharge = "CHARGE",
                                msLevel = "MSLEVEL",
                                id = "FEATURE_ID"))

neg_spec$id <- paste0(toupper(method_prefix),
                      "NEGFT",
                      str_pad(neg_spec$id, width = 5, side = "left", pad = "0"))

neg_spec$title <- neg_spec$id

# write data to new files ------------------------------------------------------
# MS1 data
neg_data_selected %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/neg_data_id.tsv"))
neg_studydesign %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/neg_studydesign_id.tsv"))
neg_anno_ms2_inhouse %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/neg_anno_ms2_inhouse_id.tsv"))
neg_anno_ms2_external %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/neg_anno_ms2_external_id.tsv"))
neg_anno_ms2_lipid %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/neg_anno_ms2_lipid_id.tsv"))
neg_anno_ms1_inhouse %>%
  write_tsv(paste0(project_path,
                   "/05_results/01_ids/neg_anno_ms1_inhouse_id.tsv"))

# MS2 data
export(neg_spec,
       backend = MsBackendMgf(),
       file = paste0(project_path,
                     "/05_results/01_ids/neg_sirius_id.mgf"),
       mapping = c(rtime = "RTINSECONDS",
                   acquisitionNum = "SCANS",
                   precursorMz = "PEPMASS",
                   precursorIntensity = "PEPMASSINT",
                   precursorCharge = "CHARGE",
                   msLevel = "MSLEVEL",
                   id = "FEATURE_ID",
                   title = "TITLE"))

# ==============================================================================
# clean up, remove datasets
# ==============================================================================
# positive mode
rm(pos_data, pos_data_selected, pos_studydesign,
   pos_anno, pos_anno_ms2_inhouse, pos_anno_ms2_external,
   pos_anno_ms2_lipid, pos_anno_ms1_inhouse, pos_spec,
   pos_sample_names)

# negative mode
rm(neg_data, neg_data_selected, neg_studydesign,
   neg_anno, neg_anno_ms2_inhouse, neg_anno_ms2_external,
   neg_anno_ms2_lipid, neg_anno_ms1_inhouse, neg_spec,
   neg_sample_names)