#!/usr/bin/env r

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

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

#-----------
## Load data
#-----------

# Annotated metabolites
annotations_dir <- "annotations/"
annotations <-
  list.files(annotations_dir) %>%
  str_subset(".tsv.gz$") %>%
  map(\(x) vroom(
    paste0(annotations_dir, x),
    col_types = "cccccccclcicc",
    progress = FALSE
  ) %>% mutate(file_name = x)) %>%
  list_rbind()

annotations <-
  annotations %>%
  mutate(
    biosystem_id =
      str_extract(file_name, "_._") %>%
      str_replace_all("_", "")
  ) %>%
  select(-file_name) %>%
  select(biosystem_id, assay_type, everything())

#d_dmagna_annot <- read_tsv("annotations/annotations_D_Dmagna.tsv") %>% mutate(biosystem_id = "D")
#
#write_tsv(d_dmagna_annot, "annotations/taggedannotations/annotations_D_Dmagna.tsv")
#
#E_Celegans_annot <- read_tsv("annotations/annotations_E_Celegans.tsv") %>% mutate(biosystem_id = "E")
#
#write_tsv(E_Celegans_annot, "annotations/taggedannotations/annotations_E_Celegans.tsv")
#
#F_Dmelanogaster_annot <- read_tsv("annotations/annotations_F_Dmelanogaster.tsv") %>% mutate(biosystem_id = "F")
#
#write_tsv(F_Dmelanogaster_annot, "annotations/taggedannotations/annotations_F_Dmelanogaster.tsv")
#
#G_HepG2_rerun_annot <- read_tsv("annotations/annotations_G_HepG2_rerun.tsv") %>% mutate(biosystem_id = "G")
#
#write_tsv(G_HepG2_rerun_annot, "annotations/taggedannotations/annotations_G_HepG2_rerun.tsv")
#
#G_HepG2_annot <- read_tsv("annotations/annotations_G_HepG2.tsv") %>% mutate(biosystem_id = "G")
#
#write_tsv(G_HepG2_annot, "annotations/taggedannotations/annotations_G_HepG2.tsv")
#
#M_Dmelanogaster_annot <- read_tsv("annotations/annotations_M_Dmelanogaster.tsv") %>% mutate(biosystem_id = "M")
#
#write_tsv(M_Dmelanogaster_annot , "annotations/taggedannotations/annotations_M_Dmelanogaster.tsv")
#
#R_Danio_annot <- read_tsv("annotations/annotations_R_Danio.tsv") %>% mutate(biosystem_id = "R")
#
#write_tsv(R_Danio_annot , "annotations/taggedannotations/annotations_R_Danio.tsv")
#
#X_Xenopus_annot <- read_tsv("annotations/annotations_X_Xenopus.tsv") %>% mutate(biosystem_id = "X")
#
#write_tsv(X_Xenopus_annot , "annotations/taggedannotations/annotations_X_Xenopus.tsv")
#
#
#
#annotations_dir <- "annotations/taggedannotations/"
#
#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 = "ccccccccddddc",
    progress = FALSE
  )) %>%
  list_rbind() %>%
  left_join(
    annotations %>%
      filter(assay_type == "HILIC") %>%
      select(metabolite_id, biosystem_id) %>%
      mutate(annotated = TRUE) %>%
      distinct()
  )

#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) %>%
#
#      mutate(annotated = TRUE) %>%
#
#      distinct()
#
#  )



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

#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) %>%
#
#      mutate(annotated = TRUE) %>%
#
#      distinct()
#
#  )



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

# no exposure_batch_id or dose_code?
#unannotated_tested <- metabolomics %>%
#  select(metabolite_id,compound_id, timepoint_level, biosystem_id,ratio,padj) %>%
#  distinct() %>%
#  mutate(assay_type = "HILIC") %>%
#  bind_rows(
#    lipidomics %>%
#      select(metabolite_id, compound_id, timepoint_level, biosystem_id,ratio,padj) %>%
#      distinct() %>%
#      mutate(assay_type = "LIPIDS")
#  )
#
#
#
#unannotated_tested %>%
#  filter(assay_type == "HILIC") %>%
#  vroom_write("logFC_tables/logFC_metabolomics_unannotated.tsv")
#
#
#
#unannotated_tested %>%
#  filter(assay_type == "LIPIDS") %>%
#  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, assay_type) %>%
  distinct() %>%
  bind_rows(
    lipidomics %>%
      filter(annotated) %>%
      select(metabolite_id, exposure_batch_id, biosystem_id, assay_type) %>%
      distinct()
  )



# The ones with unique compound name
compounds_matched <-
  annotated_tested %>%
  left_join(annotations) %>%
  select(
    metabolite_id, compound, assay_type, biosystem_id,exposure_batch_id
  )

unique_compounds <- compounds_matched%>%
  #select(-metabolite_id) %>%
  group_by(compound, metabolite_id, assay_type, biosystem_id, exposure_batch_id) %>%
  summarise(N = n()) %>%
  filter(N == 1) %>%
  ungroup() %>%
  select(-N)





# The ones with unique metabolite id
pure_annotations <- annotations %>%
  select(metabolite_id, compound) %>%
  distinct() %>%
  group_by(metabolite_id) %>%
  filter(n() == 1) %>%
  ungroup()

unique_metabolite_ids <- annotated_tested %>%
  inner_join(pure_annotations, by = "metabolite_id")



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



# Save the tables in the correct format
annotated_unique %>%
  filter(assay_type == "HILIC") %>%
  select(-assay_type) %>%
  left_join(metabolomics) %>%
  vroom_write("logFC_tables/logFC_metabolomics_uniquely_annotated.tsv")



annotated_unique %>%
  filter(assay_type == "LIPIDS") %>%
  select(-assay_type) %>%
  left_join(lipidomics) %>%
  vroom_write("logFC_tables/logFC_lipidomics_uniquely_annotated.tsv")