#!/usr/bin/env r

#-------------------------
## Load required libraries
#-------------------------
library(tidyverse)
library(vroom)

#-------------------
## Working directory
#-------------------
setwd("/users/project/precisiontox_008279_no_backup/production_metabolomics_data/")

#-----------
## Load data
#-----------
#orgs <- c("CElegans", "Daphnia", "Drosophila", "Danio", "Xenopus", "HepG2")

# Annotated metabolites
annotations_dir <- "annotations/"
annotations <-
  list.files(annotations_dir) %>%
  map(\(x) vroom(
    paste0(annotations_dir, x),
    col_types = "ccccccccclciccc",
    progress = FALSE
  )) %>%
  list_rbind()

# Results HILIC metabolomics
metabolomics_dir <- "HILIC/"
metabolomics <-
  list.files(metabolomics_dir) %>%
  map(\(x) vroom(
    paste0(metabolomics_dir, x),
    col_types = "cccccccdddd",
    progress = FALSE
  )) %>%
  list_rbind() %>%
  filter(compound_id != "000") %>%
  left_join(
    annotations %>%
      filter(assay_type == "HILIC") %>%
      select(metabolite_id, biosystem_id, exposure_batch_id) %>%
      mutate(annotated = TRUE) %>%
      distinct()
  )

# Results lipidomics
lipidomics_dir <- "LIPIDS/"
lipidomics <-
  list.files(lipidomics_dir) %>%
  map(\(x) vroom(
    paste0(lipidomics_dir, x),
    col_types = "cccccccdddd",
    progress = FALSE
  )) %>%
  list_rbind() %>%
  filter(compound_id != "000") %>%
  left_join(
    annotations %>%
      filter(assay_type == "LIPIDS") %>%
      select(metabolite_id, biosystem_id, exposure_batch_id) %>%
      mutate(annotated = TRUE) %>%
      distinct()
  )


#-------------------------------------------------------
## Save logFC data for all tested compounds/metabolites
#-------------------------------------------------------
metabolomics %>%
  select(metabolite_id:compound_id, ratio) %>%
  arrange(biosystem_id, compound_id, metabolite_id, exposure_batch_id) %>%
  unite("id", compound_id, timepoint_level, sep = "", remove = TRUE) %>%
  pivot_wider(values_from = "ratio", names_from = "id") %>%
  vroom_write("logFC_tables/logFC_metabolomics_unannotated.tsv")

lipidomics %>%
  select(metabolite_id:compound_id, ratio) %>%
  arrange(biosystem_id, compound_id, metabolite_id, exposure_batch_id) %>%
  unite("id", compound_id, timepoint_level, sep = "", remove = TRUE) %>%
  pivot_wider(values_from = "ratio", names_from = "id") %>%
  vroom_write("logFC_tables/logFC_lipidomics_unannotated.tsv")

#--------------------------------------------------------------
## Save logFC data for uniquely annotated compounds/metabolites
#--------------------------------------------------------------
annotated_tested <-
  metabolomics %>%
  filter(annotated) %>%
  select(metabolite_id, exposure_batch_id, biosystem_id) %>%
  distinct() %>%
  mutate(assay_type = "HILIC") %>%
  bind_rows(
    lipidomics %>%
      filter(annotated) %>%
      select(metabolite_id, exposure_batch_id, biosystem_id) %>%
      distinct() %>%
      mutate(assay_type = "LIPIDS")
  )

# The ones with unique compound name
unique_compounds <-
  annotated_tested %>%
  left_join(annotations) %>%
  select(
    metabolite_id, compound, assay_type, biosystem_id, exposure_batch_id
  ) %>%
  #select(-metabolite_id) %>%
  group_by(compound, assay_type, biosystem_id, exposure_batch_id) %>%
  summarise(N = n()) %>%
  filter(N == 1) %>%
  ungroup() %>%
  select(-N)

# The ones with unique metabolite id
unique_metabolite_ids <-
  annotated_tested %>%
  left_join(annotations) %>%
  select(
    metabolite_id, compound, assay_type, biosystem_id, exposure_batch_id
  ) %>%
  #select(-compound) %>%
  group_by(metabolite_id, assay_type, biosystem_id, exposure_batch_id) %>%
  summarise(N = n()) %>%
  filter(N == 1) %>%
  ungroup() %>%
  select(-N)

# The intersection of the previous two
annotated_unique <-
  annotations %>%
  inner_join(unique_metabolite_ids) %>%
  inner_join(unique_compounds) %>%
  select(
    metabolite_id, compound, assay_type, biosystem_id, exposure_batch_id
  )

# How many metabolites/compounds per species
annotated_unique %>%
  select(biosystem_id, assay_type, compound) %>%
  distinct() %>%
  group_by(biosystem_id, assay_type) %>%
  summarise(N = n())

# How many metabolites/compounds per combination of species
annotated_unique %>%
  select(biosystem_id, assay_type, compound) %>%
  distinct() %>%
  group_by(assay_type, compound) %>%
  summarise(
    Norgs = n(),
    orgs = paste(biosystem_id, collapse = ",")
  ) %>%
  #pull(Norgs) %>% table
  filter(Norgs > 1) %>%
  group_by(assay_type, orgs, Norgs) %>%
  summarise(Ncompounds = n()) %>%
  arrange(desc(Norgs), desc(Ncompounds)) %>%
  print(n = Inf)

# Save the tables in the correct format
annotated_unique %>%
  filter(assay_type == "HILIC") %>%
  select(-assay_type) %>%
  left_join(metabolomics) %>%
  unite(
    "id", biosystem_id, exposure_batch_id, compound_id, dose_code,
    timepoint_level, sep = "", remove = FALSE
  ) %>%
  arrange(biosystem_id, compound_id, compound, exposure_batch_id) %>%
  select(compound, id, ratio) %>%
  pivot_wider(values_from = "ratio", names_from = "id") %>%
  vroom_write("logFC_tables/logFC_metabolomics_uniquely_annotated.tsv")

annotated_unique %>%
  filter(assay_type == "LIPIDS") %>%
  select(-assay_type) %>%
  left_join(lipidomics) %>%
  unite(
    "id", biosystem_id, exposure_batch_id, compound_id, dose_code,
    timepoint_level, sep = "", remove = FALSE
  ) %>%
  arrange(biosystem_id, compound_id, compound, exposure_batch_id) %>%
  select(compound, id, ratio) %>%
  pivot_wider(values_from = "ratio", names_from = "id") %>%
  vroom_write("logFC_tables/logFC_lipidomics_uniquely_annotated.tsv")