library(cowplot)
library(ggsignif)
library(scales)
library(tidyverse)

# Sourcing analysis scripts
source("analysis_formatting.R")
source("analysis_processing.R")
source("analysis_plotting.R")

#==============================================================================>
# Reading data

if (! "data.RData" %in% list.files() ) {
  
  samples <- c(
    "C" = "Control",
    "S1" = "ANP",
    "S3" = "ANP + A7"
  )

  d <- read_file(
      './data.csv', 
      sample = "(?<=\\[).*?(?=\\])",
      charge = "(pos)|(neg)",
      full = ".*"
    ) %>% 
    separate(sample, c("sample", "day"), sep="-") %>%
    mutate(
      day = recode(day, "4" = "1", "5" = "2", "6" = "3"),
      sample = toupper(sample),
      sample = as.character(sample),
      id = paste(sample, day, sep = "-")
    ) %>%
    filter(sample != "S2") %>%
    mutate(sample = factor(samples[sample], levels = samples))

  # Adding grouping sets 
  selection <- c("Cer", "DG", "LPC", "LPE", "PC", "PE", "PG", "PI", "SM", "TG")

  d <- filter(d, class %in% selection)

  # Checking proper parsing
  for ( column in c("sample", "day", "charge", "set") ) {
    print(column)
    x <- d[[column]]
    print(unique(x))
    print(sum(is.na(x)))
    print(unique(d$full[is.na(x)]))
  }

  # Filtering sn and pq (and deuterated standards)
  d <- d %>%
    group_by(line, sample) %>%
    filter_quality(sn.min = 2, pq.min = 0.8, n.min = 2) %>%
    filter(! str_detect(fatty.acid, "[D]"))

  # Calculating maximum number of FA chains observed per class
  d.chains <- d %>%
    group_by(class) %>%
    summarize(max.chains = max(str_count(fatty.acid, "_")))

  save(d, d.chains, file = "data.RData")
} else {
  load(file = "data.RData")
}

#==============================================================================
# By sample

d %>%
  group_by(sample, id, class) %>%
  summarize(norm.area = sum(norm.area, na.rm = TRUE)) %>%
  pivot_wider(names_from = "class", values_from = "norm.area") %>%
  write.table(
    file = 'replicate_sums.csv', sep = ",", 
    quote = FALSE, row.names = FALSE
  ) 


#==============================================================================
# Unique species

d.areas <- d %>%
  group_by(sample, class, fatty.acid, line) %>%
  average_and_sum("norm.area") %>%
  group_by(class, fatty.acid) %>%
  filter_number(n.min = 2) %>%
  filter(
    ( (sample == "Control") & (! is.na(norm.area)) ) | 
    (sample != "Control")
  ) %>%
  group_by(class) %>%
  filter_number(n.min = 6)

f <- function(x, y) {log2(x/y)}

d.control <- d.areas %>%
  filter(sample == "Control") %>%
  filter(!is.na(norm.area)) %>%
  select(class, fatty.acid, control.norm.area = norm.area)

d.stats <- d.areas %>%
  filter(sample != "Control") %>%
  right_join(d.control, by = c("class", "fatty.acid")) %>%
  mutate(norm.area = log2(norm.area/control.norm.area)) %>%
  group_by(class, sample) %>%
  do(test_wilcox(.$norm.area, .$control.norm.area)) %>%
  ungroup() %>%
  mutate(
    reference = factor("Control", levels = levels(sample)), 
    label = test_significance(p)
  )


#------------------------------------------------------------------------------
# Writing tables

d.temp <- d.areas %>%
  pivot_wider(names_from = "sample", values_from = "norm.area")

write.table(
  d.temp,
  file = 'areas_all.csv', sep = ",", 
  quote = FALSE, row.names = FALSE
)

d.temp <- d.areas %>%
  filter(class %in% c("DG", "TG")) %>%
  pivot_wider(names_from = "sample", values_from = "norm.area")

write.table(
  d.temp,
  file = 'areas_dgtg.csv', sep = ",", 
  quote = FALSE, row.names = FALSE
)

d.stats %>%
  select(class, sample, fold.change, p, label) %>%
  mutate(p = sprintf("%.4f", as.numeric(p)),
         fold.change = sprintf("%.2f", as.numeric(fold.change))) %>%
  write.table(file = 'significance.csv', sep = ",", quote = FALSE,
              row.names = FALSE)


#------------------------------------------------------------------------------
# Box plots

# All areas
p.areas <- plot_boxes(
    d.areas, x = "class", y = "norm.area", fill = "sample", 
    log10 = TRUE, alpha = 0.5
  ) +
  ylab("Normalized area") +
  xlab("Lipid class") +
  theme(legend.position = c(0.1, 0.85), legend.title = element_blank())

ggsave('boxplots_all.tiff', width = 10, height = 5)


p.sig <- add_comparisons(
    p.areas, d.stats, x = "class", x.split = "sample", 
    y = "reference.max", tick.height = .1
  ) +
  ylab("Normalized area") +
  xlab("Lipid class") +
  theme(legend.position = c(0.1, 0.85), legend.title = element_blank())

ggsave('boxplots_all_sig.tiff', width = 10, height = 5)
ggsave('boxplots_all_linear.tiff', p.areas + scale_y_continuous(), width = 10, height = 5)


# Just DG/TG
p.areas <- plot_boxes(
    filter(d.areas, class %in% c("DG", "TG")), 
    x = "class", y = "norm.area", fill = "sample", 
    log10 = TRUE, alpha = 0.5
  ) +
  ylab("Normalized area") +
  xlab("Lipid class") +
  theme(legend.position = "bottom")

p.sig <- add_comparisons(
    p.areas, filter(d.stats, class %in% c("DG", "TG")), 
    x = "class", x.split = "sample", 
    y = "reference.max", tick.height = .1
  ) +
  ylab("Normalized area") +
  xlab("Lipid class") +
  theme(legend.position = "bottom", legend.title = element_blank())

ggsave('boxplots_dgtg.tiff', width = 5, height = 5)
ggsave('boxplots_dgtg_linear.tiff', p.sig + scale_y_continuous(), width = 5, height = 5)

# Redoing significance with full pairing
d.stats <- d.areas %>%
  filter(sample != "Control") %>%
  right_join(d.control, by = c("class", "fatty.acid")) %>%
  mutate(norm.area = log2(norm.area/control.norm.area)) %>%
  group_by(class, fatty.acid) %>%
  filter_number(n.min = 2) %>%
  group_by(class, sample) %>%
  do(test_wilcox(.$norm.area, .$control.norm.area)) %>%
  ungroup() %>%
  mutate(
    reference = factor("Control", levels = levels(sample)), 
    label = test_significance(p)
  )

p.sig <- add_comparisons(
    p.areas, filter(d.stats, class %in% c("DG", "TG")), 
    x = "class", x.split = "sample", 
    y = "reference.max", tick.height = .1
  ) +
  ylab("Normalized area") +
  xlab("Sample") +
  theme(legend.position = "bottom", legend.title = element_blank())

ggsave('boxplots_dgtg_full_paired.tiff', width = 5, height = 5)
ggsave('boxplots_dgtg_full_paired_linear.tiff', p.areas + scale_y_continuous(), width = 5, height = 5)

#==============================================================================
# By sample

d %>%
  filter(fatty.acid %in% unique(d.areas$fatty.acid)) %>%
  group_by(sample, id, class) %>%
  summarize(norm.area = sum(norm.area, na.rm = TRUE)) %>%
  pivot_wider(names_from = "class", values_from = "norm.area") %>%
  write.table(
    file = 'replicate_sums.csv', sep = ",", 
    quote = FALSE, row.names = FALSE
  ) 

#==============================================================================>
# Individual fatty acids

d.individual <- d.areas %>%
  group_by(class, fatty.acid) %>%
  filter_number(n.min = 2) %>%
  left_join(d.chains, by = c("class")) %>%
  filter(str_count(fatty.acid, "_") == max.chains) %>%
  separate(fatty.acid, c("fa1", "fa2", "fa3"), sep = "_") %>%
  pivot_longer(
    starts_with("fa"), names_to = "chain", values_to = "fatty.acid"
  ) %>%
  filter(! is.na(fatty.acid)) %>%
  group_by(sample, class, fatty.acid) %>%
  summarize(norm.area = sum(norm.area)) %>%
  ungroup() %>%
  separate(fatty.acid, c("first", "second"), sep = ":", remove = FALSE) %>%
  mutate(
    length = as.numeric(str_extract(first, "\\d+$")), 
    saturation = as.numeric(str_extract(second, "\\d+"))
  ) %>%
  arrange(length, saturation) %>%
  mutate(
    fatty.acid = as.character(fatty.acid),
    fatty.acid = factor(fatty.acid, levels = unique(fatty.acid))
  )

# Stats
d.control <- d.individual %>%
  filter(sample == "Control") %>%
  filter(!is.na(norm.area)) %>%
  select(class, fatty.acid, control.norm.area = norm.area)

d.individual %>%
  filter(sample != "Control") %>%
  right_join(d.control, by = c("class", "fatty.acid")) %>%
  mutate(norm.area = log2(norm.area/control.norm.area)) %>%
  group_by(class, sample) %>%
  do(test_wilcox(.$norm.area, .$control.norm.area)) %>%
  ungroup() %>%
  mutate(
    reference = factor("Control", levels = levels(sample)), 
    label = test_significance(p)
  ) %>%
  write.table(
    file = 'fatty_acid_chains_significance.csv', sep = ",", quote = FALSE,
    row.names = FALSE
  )


# DG
p <- plot_points(
  filter(d.individual, class %in% c("DG", "TG")), 
  x = "fatty.acid", y = "norm.area", 
  colour = "sample",  size = 2, log10 = TRUE) +
  ylab("Normalized area") +
  xlab("Fatty acid chain") +
  scale_colour_brewer("Sample", palette = "Dark2") +
  facet_wrap(~ class, scale = "free_y", ncol = 1) +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 9)) +
  coord_flip()

ggsave('fatty_acid_chains_dgtg.pdf', width = 10, height = 8)
ggsave('fatty_acid_chains_dgtg.tiff', width = 10, height = 8)

# Overall
p <- plot_points(
  group_by(d.individual, sample, fatty.acid) %>% summarize(norm.area = sum(norm.area)), 
  x = "fatty.acid", y = "norm.area", 
  colour = "sample",  size = 2, log10 = TRUE) +
  ylab("Normalized area") +
  xlab("Fatty acid chain") +
  scale_colour_brewer("Sample", palette = "Dark2") +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 9)) +
  coord_flip()

ggsave('fatty_acid_chains_overall.pdf', width = 10, height = 8)
ggsave('fatty_acid_chains_overall.tiff', width = 10, height = 8)


#==============================================================================>
# Sums

f_sum <- function(x) {
  x <- str_split(x, "_")
  carbons <- map(x, ~ sum(as.numeric(str_extract(.x, "\\d+(?=:)"))))
  saturations <- map(x, ~ sum(as.numeric(str_extract(.x, "(?<=:)\\d+"))))
  paste(carbons, saturations, sep = ":")
}


d.sums <- d.areas %>%
  group_by(class, fatty.acid) %>%
  filter_number(n.min = 2) %>%
  mutate(fatty.acid = f_sum(fatty.acid)) %>%
  ungroup() %>%
  separate(fatty.acid, c("first", "second"), sep = ":", remove = FALSE) %>%
  mutate(
    length = as.numeric(str_extract(first, "\\d+$")), 
    saturation = as.numeric(str_extract(second, "\\d+"))
  ) %>%
  arrange(length, saturation) %>%
  mutate(
    fatty.acid = as.character(fatty.acid),
    fatty.acid = factor(fatty.acid, levels = unique(fatty.acid))
  )

# Stats
d.control <- d.sums %>%
  filter(sample == "Control") %>%
  filter(!is.na(norm.area)) %>%
  select(class, fatty.acid, control.norm.area = norm.area)

d.sums %>%
  filter(sample != "Control") %>%
  right_join(d.control, by = c("class", "fatty.acid")) %>%
  mutate(norm.area = log2(norm.area/control.norm.area)) %>%
  group_by(class, sample) %>%
  do(test_wilcox(.$norm.area, .$control.norm.area)) %>%
  ungroup() %>%
  mutate(
    reference = factor("Control", levels = levels(sample)), 
    label = test_significance(p)
  ) %>%
  write.table(
    file = 'fatty_acid_total_significance.csv', sep = ",", quote = FALSE,
    row.names = FALSE
  )


# DG
p <- plot_points(
  filter(d.sums, class %in% c("DG", "TG")), 
  x = "fatty.acid", y = "norm.area", 
  colour = "sample",  size = 2, log10 = TRUE) +
  ylab("Normalized area") +
  xlab("Total carbons:double bonds") +
  scale_colour_brewer("Sample", palette = "Dark2") +
  facet_wrap(~ class, scale = "free_y", ncol = 1) +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 9)) +
  coord_flip()

ggsave('fatty_acid_total_dgtg.pdf', width = 10, height = 8)
ggsave('fatty_acid_total_dgtg.tiff', width = 10, height = 8)


# Overall
p <- plot_points(
  group_by(d.sums, sample, fatty.acid) %>% summarize(norm.area = sum(norm.area)), 
  x = "fatty.acid", y = "norm.area", 
  colour = "sample",  size = 2, log10 = TRUE) +
  ylab("Normalized area") +
  xlab("Total carbons:double bonds") +
  scale_colour_brewer("Sample", palette = "Dark2") +
  theme(legend.position = "bottom",
        axis.text.y = element_text(size = 9)) +
  coord_flip()

ggsave('fatty_acid_total_overall.pdf', width = 10, height = 12)
ggsave('fatty_acid_total_overall.tiff', width = 10, height = 12)