4_Dysregulated_Pathways_labeled.Rmd

---
title: "4_Dysregulated Pathways"
author: "Christian Ayala"
date: "4/14/2021"
output: html_document
---

This Notebook will explore the KEGG (and Metabolika) Pathways that are dysregulated

# 1. Importing Libraries

```{r libraries, message=FALSE, warning=FALSE}
library(ggpubr)
library(ggsci)
library(ggrepel)
library(pheatmap)
library(RColorBrewer)
library(scales)
library(cowplot)
library(readxl)
library(KEGGREST)
library(tidyverse)
source('functions_dysreg_path.R')
```

# 2. Import Data

```{r}
# Flag for labeled / unlabeled data, set TRUE or FALSE

label = TRUE
```

```{r set_path, message=FALSE, warning=FALSE}
# set path variables
project_dir <- getwd()
project_name <- 'Bog_labeled_all'
figures_dir <- file.path(project_dir, paste0(project_name, '_output_figures'))
tables_dir <- file.path(project_dir,  paste0(project_name, '_output_tables'))
```

This analysis uses the tables generated during the Differential Analysis script.

```{r message=FALSE, warning=FALSE}
# Load compounds table

# If the flag was set before, nothing needs to be changed here, the correct file will be automatically used
if(label == TRUE){
  compounds_table_file <- file.path(tables_dir, 'labeled_compounds_table_in_all.csv')
}else{
  compounds_table_file <- file.path(tables_dir, 'gap_filled_compounds_table.csv')
}
compounds_table <- read_csv(compounds_table_file)

# unlabeled compounds file

unlabeled_compounds_table <- read_csv(file.path(tables_dir, 'unlabeled_compounds_table.csv'))

# Load normalized abundances

norm_auc_table_file <- file.path(tables_dir, 'normalized_transformed_auc_table.csv')
norm.matrix <- read_csv(norm_auc_table_file)
norm.matrix <- column_to_rownames(norm.matrix, var = 'X1')

# Import metadata and fix names
metadata_file <- file.path(tables_dir, 'fixed_metadata.csv')
metadata <- read_csv(metadata_file)%>% 
  filter(Label == 'Labeled' | SampleID == 'BNC')

# Import table(s) obtained during the Differential analysis
table_file <- file.path(tables_dir, 'Diff_expressed_T1.csv')
T1.diff_table <- read_csv(table_file)
table_file <- file.path(tables_dir, 'Diff_expressed_T2.csv')
T2.diff_table <- read_csv(table_file)
table_file <- file.path(tables_dir, 'Diff_expressed_T3.csv')
T3.diff_table <- read_csv(table_file)

# Import diff tables of NMR data
# Import table(s) obtained during the Differential analysis
table_file <- file.path(tables_dir, 'Diff_expressed_NMR_T1.csv')
T1.NMR.diff_table <- read_csv(table_file)
table_file <- file.path(tables_dir, 'Diff_expressed_NMR_T2.csv')
T2.NMR.diff_table <- read_csv(table_file)
table_file <- file.path(tables_dir, 'Diff_expressed_NMR_T3.csv')
T3.NMR.diff_table <- read_csv(table_file)

nmr.matrix <- read_csv(file.path(tables_dir, 'NMR_matrix.csv')) %>% 
  column_to_rownames(var = 'X1')

nmr.matrix[is.na(nmr.matrix)] <- 0

```

From Compound Discover export a **2-level** table with: 

- *KEGG Pathways* in the first level with the columns: PathwayID, Pathway Name
- *Compounds* in the second level (Only way I found to be able to export the **KEGG Compounds IDs**) with the columns: Checked, Molecular Weight, Name, Formula, KEGG Compound Ids, KEGG Compound Names

```{r message=FALSE, warning=FALSE}
kegg_table_file <- file.path(tables_dir, 'KEGG Pathways_Bog_labeled_all.xlsx')
kegg_table <- read_xlsx(kegg_table_file)

# Fix column names
colnames(kegg_table) <- c('PathwayID', 'Pathway_Name', 'Name', 'Formula', 'Molecular Weight', 'KEGG_ID', 'KEGG_Name')

# Fix 2-level table into regular table

kegg_table <- kegg_table %>% 
  filter(Pathway_Name != 'Checked') %>% 
  mutate(Pathway_Name = ifelse(Pathway_Name == 'FALSE', NA, Pathway_Name)) %>% 
  fill(PathwayID, Pathway_Name, .direction = 'down') %>% 
  filter(!is.na(`Molecular Weight`)) %>% 
  select(-`Molecular Weight`)

comp_kegg <- left_join(compounds_table, kegg_table, by = c('Name', 'Formula')) %>% 
  distinct()

table_file <- file.path(tables_dir, 'Compounds_with_KEGG.csv')
write_csv(comp_kegg, table_file)
```

# 3. Analyze KEGG Pathways annotated at each time point

Merge KEGG annotation with Differential analysis data

```{r}
# This step includes filtering for significant DE metabolites (pval.adj < 0.05)
T1.kegg <- T1.diff_table %>%
  select(FeatureID, Name, Formula, log2FC, pval, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct() %>% 
  filter(pval < 0.05) %>% 
  filter(!is.na(Pathway_Name)) %>% 
  mutate(Time = 'T1')

T2.kegg <- T2.diff_table %>%
  select(FeatureID, Name, Formula, log2FC, pval, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct() %>% 
  filter(pval < 0.05)%>% 
  filter(!is.na(Pathway_Name)) %>% 
  mutate(Time = 'T2')

T3.kegg <- T3.diff_table %>%
  select(FeatureID, Name, Formula, log2FC, pval, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct() %>% 
  filter(pval < 0.05)%>% 
  filter(!is.na(Pathway_Name)) %>% 
  mutate(Time = 'T3')
```

Get not significant KEGGS for KEGG mapper

```{r}

T1.kegg_ns <- T1.diff_table %>%
  select(FeatureID, Name, Formula, log2FC, pval, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct() %>% 
  filter(pval >= 0.05) %>% 
  filter(!is.na(Pathway_Name)) %>% 
  mutate(Time = 'T1')

T2.kegg_ns <- T2.diff_table %>%
  select(FeatureID, Name, Formula, log2FC, pval, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct() %>% 
  filter(pval >= 0.05) %>% 
  filter(!is.na(Pathway_Name)) %>% 
  mutate(Time = 'T2')

T3.kegg_ns <- T3.diff_table %>%
  select(FeatureID, Name, Formula, log2FC, pval, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct() %>% 
  filter(pval >= 0.05) %>% 
  filter(!is.na(Pathway_Name)) %>% 
  mutate(Time = 'T3')

```


## 3.1 Plot number of KEGG pathways

### 3.1.1 Number of dysregulated metabolites with KEGG annotations

```{r}
# Extract significant features of each comparison based on adjusted pvalue
sig_features <- c(T1.diff_table$FeatureID[T1.diff_table$pval < 0.05],
                  T2.diff_table$FeatureID[T2.diff_table$pval < 0.05],
                  T3.diff_table$FeatureID[T3.diff_table$pval < 0.05])

sig_features <- unique(sig_features)

kegg_annotated <- compounds_table %>% 
  filter(FeatureID %in% sig_features) %>% 
  mutate(kegg_annot = ifelse(!is.na(`KEGG Pathways`), 'Annotated', 'No annotation')) %>% 
  group_by(SampleID) %>% 
  count(kegg_annot) %>% 
  filter(SampleID %in% metadata$SampleID)

kegg_annotated_plot <- plot_col(kegg_annotated, SampleID, n, kegg_annot) +
  scale_fill_jama() +
  labs(title = 'Number of dysregulated metabolites annotated in KEGG')

kegg_annotated_plot

figure_file <- file.path(figures_dir, 'dysregulated_annotated.png')
ggsave(figure_file, kegg_annotated_plot, dpi = 300)

```


```{r}
# Set colors for plot
my_colors = c('Not present in control' = '#d8365e', 
              'Upregulated' = '#d34849',
              'Downregulated' = '#005193',
              'Only present in control' = '#4b70de')

T1_counts.kegg <- T1.kegg %>% 
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>% 
  mutate(Time = 'T1')


T1_counts.kegg.plot <- plot_col(T1_counts.kegg, n, Pathway_Name, Comment, my_colors = my_colors) +
  facet_wrap(~Comment) +
  labs(title = 'Dysregulated KEGG Pathways at T1',
       y = 'Pathway Name') +
  scale_x_continuous(breaks = pretty_breaks(3))
T1_counts.kegg.plot

figure_file <- file.path(figures_dir, 'Dysreg_KEGG_T1.png')
ggsave(figure_file, T1_counts.kegg.plot, dpi = 300)

T2_counts.kegg <- T2.kegg %>% 
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>% 
  mutate(Time = 'T2')

T2_counts.kegg.plot <- plot_col(T2_counts.kegg, n, Pathway_Name, Comment, my_colors = my_colors) +
  facet_wrap(~Comment) +
  labs(title = 'Dysregulated KEGG Pathways at T2',
       y = 'Pathway Name') +
  scale_x_continuous(breaks = pretty_breaks(1))
T2_counts.kegg.plot

figure_file <- file.path(figures_dir, 'Dysreg_KEGG_T2.png')
ggsave(figure_file, T2_counts.kegg.plot, dpi = 300)

T3_counts.kegg <- T3.kegg %>% 
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>% 
  mutate(Time = 'T3')

T3_counts.kegg.plot <- plot_col(T3_counts.kegg, n, Pathway_Name, Comment, my_colors = my_colors) +
  facet_wrap(~Comment) +
  labs(title = 'Dysregulated KEGG Pathways at T3',
       y = 'Pathway Name') +
  scale_x_continuous(breaks = pretty_breaks(2))
T3_counts.kegg.plot

figure_file <- file.path(figures_dir, 'Dysreg_KEGG_T3.png')
ggsave(figure_file, T3_counts.kegg.plot, dpi = 300)

```

## 3.2 General visualization of dysregulated KEGG pathways

```{r fig.height=10}

lcms_kegg <- rbind(T1_counts.kegg, T2_counts.kegg, T3_counts.kegg) %>% 
  mutate(Source = paste(Time, Comment, sep = '_')) %>% 
  select(-Comment, -Time)

lcms_kegg <- pivot_wider(lcms_kegg, names_from = Source, values_from = n)

lcms_kegg[is.na(lcms_kegg)] <- 0

lcms_kegg <- pivot_longer(lcms_kegg, cols = 2:ncol(lcms_kegg), names_to = 'Source', values_to = 'n')

lcms_kegg <- plyr::ddply(lcms_kegg, plyr::.(Source), transform, rescale = rescale(n))

write_csv(lcms_kegg, file.path(tables_dir, 'all_kegg.csv'))

table_file <- file.path('data', 'only_pathways.csv')
pathways <- read_csv(table_file)

lcms_kegg_filt <- lcms_kegg[lcms_kegg$Pathway_Name %in% pathways$Pathway_Name,]

kegg_heatmap <- ggplot(lcms_kegg_filt,
                       aes(x = Source,
                           y = Pathway_Name,
                           fill = n)) +
  geom_tile(color = 'white') +
  coord_fixed() +
  scale_fill_distiller(palette = 'Blues', direction = 1) +
  theme_bw() +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        axis.title.x = element_blank()) +
  labs(title = 'Dysregulated KEGG Pathways',
       y = 'Pathway Name',
       fill = 'Counts')
kegg_heatmap

figure_file <- file.path(figures_dir, 'KEGG_tile_heatmap.png')
ggsave(figure_file, kegg_heatmap, height = 5)
```

## 3.3 KEGG Pathways from NMR data

Use the KEGGREST API to find the *KEGG Pathways* associated with the *KEGG Compound ID* from each metabolite identified with NMR

```{r message=FALSE, warning=FALSE, results='hide'}

T1.NMR.diff_table <- filter(T1.NMR.diff_table, !is.na(`KEGG Compound ID`))

T1_NMR.kegg <- tibble(Pathway_Name = NA, KEGG_ID = NA, Name = NA, Comment = NA, Time = 'T1', .rows = 0)
for(i in 1:length(T1.NMR.diff_table$`KEGG Compound ID`)){
  k_id <- T1.NMR.diff_table$`KEGG Compound ID`[i]
  query <- keggGet(k_id)
  pathway <- query[[1]]$PATHWAY
  if(is.null(pathway)){
    pathway <- NA
  }
  t_tibble <- tibble(Pathway_Name = pathway,
                     KEGG_ID = k_id,
                     Name = T1.NMR.diff_table$Name[i],
                     Comment = T1.NMR.diff_table$Comment[i],
                     Time = 'T1')
  T1_NMR.kegg <- rbind(T1_NMR.kegg, t_tibble)
  print(i)
}

T2.NMR.diff_table <- filter(T2.NMR.diff_table, !is.na(`KEGG Compound ID`))

T2_NMR.kegg <- tibble(Pathway_Name = NA, KEGG_ID = NA, Name = NA, Comment = NA, Time = 'T2', .rows = 0)
for(i in 1:length(T2.NMR.diff_table$`KEGG Compound ID`)){
  k_id <- T2.NMR.diff_table$`KEGG Compound ID`[i]
  query <- keggGet(k_id)
  pathway <- query[[1]]$PATHWAY
  if(is.null(pathway)){
    pathway <- NA
  }
  t_tibble <- tibble(Pathway_Name = pathway,
                     KEGG_ID = k_id,
                     Name = T2.NMR.diff_table$Name[i],
                     Comment = T2.NMR.diff_table$Comment[i],
                     Time = 'T2')
  T2_NMR.kegg <- rbind(T2_NMR.kegg, t_tibble)
  print(i)
}

T3.NMR.diff_table <- filter(T3.NMR.diff_table, !is.na(`KEGG Compound ID`))

T3_NMR.kegg <- tibble(Pathway_Name = NA, KEGG_ID = NA, Name = NA, Comment = NA, Time = 'T3', .rows = 0)
for(i in 1:length(T3.NMR.diff_table$`KEGG Compound ID`)){
  k_id <- T3.NMR.diff_table$`KEGG Compound ID`[i]
  query <- keggGet(k_id)
  pathway <- query[[1]]$PATHWAY
  if(is.null(pathway)){
    pathway <- NA
  }
  t_tibble <- tibble(Pathway_Name = pathway,
                     KEGG_ID = k_id,
                     Name = T3.NMR.diff_table$Name[i],
                     Comment = T3.NMR.diff_table$Comment[i],
                     Time = 'T3')
  T3_NMR.kegg <- rbind(T3_NMR.kegg, t_tibble)
  print(i)
}

unlabeled_nmr <- read_csv(file.path(tables_dir, 'unlabeled_NMR.csv')) %>% 
  mutate(Comment = 'Unlabeled') %>% 
  select(-SampleID, -Label, -Comp, -concentration) %>% 
  distinct()

unlabeled_nmr <- filter(unlabeled_nmr, !is.na(`KEGG Compound ID`))

unlabeled_nmr_kegg <- tibble(Pathway_Name = NA, KEGG_ID = NA, Name = NA, Comment = NA, Time = NA, .rows = 0)
for(i in 1:length(unlabeled_nmr$`KEGG Compound ID`)){
  k_id <- unlabeled_nmr$`KEGG Compound ID`[i]
  time <- unlabeled_nmr$Time[i]
  query <- keggGet(k_id)
  pathway <- query[[1]]$PATHWAY
  if(is.null(pathway)){
    pathway <- NA
  }
  t_tibble <- tibble(Pathway_Name = pathway,
                     KEGG_ID = k_id,
                     Name = T3.NMR.diff_table$Name[i],
                     Comment = 'Unlabeled',
                     Time = time)
  unlabeled_nmr_kegg <- rbind(unlabeled_nmr_kegg, t_tibble)
  print(i)
}
  
```

Plotting the Heatmap of the KEGG Pathways associated with NMR data

```{r fig.height=10}

T1_NMR_counts.kegg <- T1_NMR.kegg %>% 
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>% 
  filter(n > 2) %>%
  mutate(Time = 'T1')

T2_NMR_counts.kegg <- T2_NMR.kegg %>% 
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>% 
  filter(n > 2) %>%
  mutate(Time = 'T2')

T3_NMR_counts.kegg <- T3_NMR.kegg %>% 
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>% 
  filter(n > 2) %>%
  mutate(Time = 'T3')

nmr_kegg <- rbind(T1_NMR_counts.kegg, T2_NMR_counts.kegg, T3_NMR_counts.kegg) %>% 
  mutate(Source = paste(Time, Comment, sep = '_')) %>% 
  select(-Comment, -Time)

nmr_kegg <- pivot_wider(nmr_kegg, names_from = Source, values_from = n)

nmr_kegg[is.na(nmr_kegg)] <- 0

nmr_kegg <- pivot_longer(nmr_kegg, cols = 2:ncol(nmr_kegg), names_to = 'Source', values_to = 'n')

nmr_kegg <- plyr::ddply(nmr_kegg, plyr::.(Source), transform, rescale = rescale(n))

kegg_nmr_heatmap <- ggplot(nmr_kegg,
                       aes(x = Source,
                           y = Pathway_Name,
                           fill = n)) +
  geom_tile(color = 'white') +
  coord_fixed() +
  scale_fill_distiller(palette = 'Blues', direction = 1) +
  theme_bw() +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5),
        axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
        axis.title.x = element_blank()) +
  labs(title = 'Dysregulated KEGG Pathways',
       y = 'Pathway Name',
       fill = 'Scaled counts')
kegg_nmr_heatmap

figure_file <- file.path(figures_dir, 'KEGG_tile_heatmap_with_NMR.png')
ggsave(figure_file, kegg_nmr_heatmap, height = 5)
```

## 3.4 KEGG Pathways from the sum of LC-MS2 and NMR data

```{r fig.height=10}
# Join both datasets together
T1_all_counts <- T1.kegg %>% 
  select(Pathway_Name, KEGG_ID, Name, Comment, Time) %>% 
  rbind(T1_NMR.kegg) %>%
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>%
  mutate(Time = 'T1')

T2_all_counts <- T2.kegg %>% 
  select(Pathway_Name, KEGG_ID, Name, Comment, Time) %>% 
  rbind(T2_NMR.kegg) %>%
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>%
  mutate(Time = 'T2')

T3_all_counts <- T3.kegg %>% 
  select(Pathway_Name, KEGG_ID, Name, Comment, Time) %>% 
  rbind(T3_NMR.kegg) %>%
  group_by(Comment) %>% 
  count(Pathway_Name) %>% 
  ungroup() %>%
  mutate(Time = 'T3')

all_kegg <- rbind(T1_all_counts, T2_all_counts, T3_all_counts) %>% 
  mutate(Source = paste(Time, Comment, sep = '_')) %>% 
  select(-Comment, -Time)

all_kegg <- pivot_wider(all_kegg, names_from = Source, values_from = n) %>% 
  filter(!is.na(Pathway_Name)) %>% 
  column_to_rownames(var = 'Pathway_Name')

all_kegg[is.na(all_kegg)] <- 0

all_kegg <- all_kegg %>% 
  rownames_to_column(var = 'Pathway_Name')

all_kegg <- all_kegg %>% 
  pivot_longer(cols = 2:ncol(all_kegg), names_to = 'Source', values_to = 'n') %>% 
  mutate(dysreg = gsub('.*_', '', Source),
         Time = gsub('_.*', '', Source))
  

all_kegg <- plyr::ddply(all_kegg, plyr::.(Source), transform, rescale = rescale(n))

table_file <- file.path(tables_dir, 'KEGG_Pathways_NMR_LC.csv')
write_csv(all_kegg, table_file)

filt_kegg <- all_kegg[all_kegg$Pathway_Name %in% pathways$Pathway_Name,]
filt_kegg$Source <- factor(filt_kegg$Source, 
                           levels = c('T1_Only present in control', 'T1_Downregulated', 'T1_Upregulated',
                                      'T1_Not present in control', 
                                      'T2_Only present in control', 'T2_Downregulated', 'T2_Upregulated',
                                      'T2_Not present in control',
                                      'T3_Only present in control', 'T3_Downregulated', 'T3_Upregulated',
                                      'T3_Not present in control'))

my_colors = c('Not present in control' = '#d8365e', 
              'Upregulated' = '#d34849',
              'Downregulated' = '#005193',
              'Only present in control' = '#4b70de')

kegg_all_heatmap <- ggplot(filt_kegg) +
  geom_tile(aes(x = dysreg,
                y = Pathway_Name,
                fill = rescale,
                color = dysreg),
            size = .85) +
  coord_fixed() +
  scale_fill_distiller(palette = 'Blues', direction = 1) +
  scale_color_manual(values = my_colors) +
  theme_bw() +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank()) +
  labs(title = 'KEGG Pathways (LC_MS/MS and NMR data)',
       y = 'Pathway Name',
       fill = 'Scaled counts') +
  guides(color = FALSE) +
  facet_wrap(~ Time)
kegg_all_heatmap

figure_file <- file.path(figures_dir, 'KEGG_tile_heatmap_all_data.png')
ggsave(figure_file, kegg_all_heatmap, height = 10)


### Individual heatmaps for each time point
kegg_all_heatmap_t1 <- filt_kegg %>% 
  filter(Time == 'T1') %>% 
  ggplot() +
  geom_tile(aes(x = Source,
                y = Pathway_Name,
                fill = rescale,
                color = dysreg),
            size = .85) +
  coord_fixed(ratio = 0.5) +
  scale_fill_gradient(low = "snow1", high = '#EFC000FF')  +
  scale_color_manual(values = my_colors) +
  theme_bw() +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank()) +
  labs(title = 'KEGG Pathways at T1',
       y = 'Pathway Name',
       fill = 'Scaled counts') +
  guides(color = FALSE)
kegg_all_heatmap_t1

figure_file <- file.path(figures_dir, 'KEGG_tile_heatmap_T1.png')
ggsave(figure_file, kegg_all_heatmap_t1, height = 10)

kegg_all_heatmap_t2 <- filt_kegg %>% 
  filter(Time == 'T2') %>% 
  ggplot() +
  geom_tile(aes(x = Source,
                y = Pathway_Name,
                fill = rescale,
                color = dysreg),
            size = .85) +
  coord_fixed(ratio = 0.5) +
  scale_fill_gradient(low = 'snow1', high = '#868686FF')  +
  scale_color_manual(values = my_colors) +
  theme_bw() +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank()) +
  labs(title = 'KEGG Pathways at T2',
       y = 'Pathway Name',
       fill = 'Scaled counts') +
  guides(color = FALSE)
kegg_all_heatmap_t2

figure_file <- file.path(figures_dir, 'KEGG_tile_heatmap_T2.png')
ggsave(figure_file, kegg_all_heatmap_t2, height = 10)

kegg_all_heatmap_t3 <- filt_kegg %>% 
  filter(Time == 'T3') %>% 
  ggplot() +
  geom_tile(aes(x = Source,
                y = Pathway_Name,
                fill = rescale,
                color = dysreg),
            size = .85) +
  coord_fixed(ratio = 0.5) +
  scale_fill_gradient(low = 'snow1', high = '#CD534CFF')  +
  scale_color_manual(values = my_colors) +
  theme_bw() +
  theme(plot.title = element_text(face = 'bold', hjust = 0.5),
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        axis.title.x = element_blank()) +
  labs(title = 'KEGG Pathways at T3',
       y = 'Pathway Name',
       fill = 'Scaled counts') +
  guides(color = FALSE)
kegg_all_heatmap_t3

figure_file <- file.path(figures_dir, 'KEGG_tile_heatmap_T3.png')
ggsave(figure_file, kegg_all_heatmap_t3, height = 10)


```

# 4. Heatmaps of each Dysregulated KEGG Pathway

```{r}
col_annot <- column_to_rownames(metadata, var = 'SampleID') 
mapcolor <- colorRampPalette(brewer.pal(11, 'RdYlBu'))(100)[100:1]
annot_colors <- list(
  Time = c(T0 ='#0073C2FF', T1 = '#EFC000FF', T2 = '#868686FF', T3 = '#CD534CFF'))

diff_summary <- rbind(T1.diff_table %>% 
                        select(FeatureID, Name, Comment, pval) %>% 
                        distinct() %>% 
                        mutate(Timepoint = 'T1'),
                      T1.diff_table %>% 
                        select(FeatureID, Name, Comment, pval) %>% 
                        distinct() %>% 
                        mutate(Timepoint = 'T2'),
                      T1.diff_table %>% 
                        select(FeatureID, Name, Comment, pval) %>% 
                        distinct() %>% 
                        mutate(Timepoint = 'T3'))

nmr_summary <- rbind(T1_NMR.kegg, T2_NMR.kegg, T3_NMR.kegg) %>% 
  mutate(FeatureID = Name,
         name4plot = Name,
         pval = 1) %>% 
  select(FeatureID, Name, Pathway_Name, name4plot, pval, Time) %>% 
  rename(Timepoint = Time)

sig_pathways <- filt_kegg %>% 
  pull(Pathway_Name)

sig_pathways <- unique(sig_pathways)

```

```{r}

col_annot <- column_to_rownames(metadata, var = 'SampleID') %>% 
  select(-Label)
mapcolor <- colorRampPalette(brewer.pal(11, 'RdYlBu'))(100)[100:1]

figure_file <- file.path(figures_dir, 'All_features_heatmap.pdf')

annot_colors <- list(
  Time = c(T0 ='#0073C2FF', T1 = '#EFC000FF', T2 = '#868686FF', T3 = '#CD534CFF'),
  Comp = c(litter = 'green4', peat_and_litter = 'chocolate4'),
  Origin = c(LC_MS2 = '#0086cb', NMR = '#d23936'))

dev.off()
for(path in sig_pathways){
  path_table_t <- comp_kegg %>% 
    select(FeatureID, Name, Pathway_Name, name4plot) %>% 
    left_join(diff_summary, by = c('FeatureID', 'Name')) %>% 
    select(-Comment) %>% 
    distinct() 
  
  path_table <- rbind(path_table_t, nmr_summary)%>% 
    filter(Pathway_Name == path)
  
  path.matrix_t <- norm.matrix[rownames(norm.matrix) %in% path_table$FeatureID,] %>% 
    select(all_of(metadata$SampleID))
  
  path.matrix_t2 <- nmr.matrix[rownames(nmr.matrix) %in% path_table$FeatureID,]
  
  path.matrix <- rbind(path.matrix_t, path.matrix_t2) %>% 
    select(-BNC) %>% 
    relocate(BNE)
  
  rownames(path.matrix) <- path_table$name4plot[match(rownames(path.matrix), path_table$FeatureID)]
  
  sig.matrix <- path_table %>% 
    left_join(metadata %>% 
                select(Time, SampleID), 
              by = c('Timepoint' = 'Time')) %>% 
    select(-Timepoint) %>% 
    pivot_wider(names_from = 'SampleID', values_from = 'pval') %>% 
    select(-Name, -name4plot, -Pathway_Name) %>%
    column_to_rownames(var = 'FeatureID')
  
  sig.matrix <- as.data.frame(ifelse(sig.matrix < 0.05, '*', '')) %>% 
    mutate(BNC = '',
           BNE = '') %>% 
    select(all_of(metadata$SampleID)) %>% 
    select(-BNC)
  rownames(sig.matrix) <- path_table$name4plot[match(rownames(sig.matrix), path_table$FeatureID)]
  sig.matrix <- sig.matrix[rownames(path.matrix),]
  
  if(nrow(path.matrix) > 2){
    clust = TRUE
  } else{
    clust = FALSE
  }
  
  figure_file <- file.path(figures_dir, paste0('heatmap_', str_replace(path, '/', '-'), '.pdf'))
  pdf(figure_file, width = 8, height = 4.5)
  pheatmap(path.matrix,
           cluster_rows = clust,
           cluster_cols = FALSE,
           clustering_distance_rows = 'correlation',
           clustering_distance_cols = 'correlation',
           scale = 'row',
           annotation_col = col_annot,
           annotation_colors = annot_colors,
           color = mapcolor,
           fontsize_row = 8,
           fontsize_number = 15,
           main = path,
           display_numbers = sig.matrix
  )
  dev.off()
  print(paste0(path, ' has finished'))
}
```

# 5. Extract KEGG IDs

```{r}

lc_kegg_ids <- rbind(T1.kegg, T2.kegg, T3.kegg) %>% 
  select(KEGG_ID, Time, Comment, Pathway_Name) %>% 
  #mutate(Source = paste(Time, Comment, sep = '_')) %>% 
  #select(-Comment, - Time) %>% 
  separate_rows(KEGG_ID, sep = '; ') %>% 
  group_by(Time, Comment) %>% 
  distinct() %>% 
  ungroup()

lc_kegg_ns <- rbind(T1.kegg_ns, T2.kegg_ns, T3.kegg_ns) %>% 
  select(KEGG_ID, Time, Comment, Pathway_Name) %>% 
  #mutate(Source = paste(Time, Comment, sep = '_')) %>% 
  #select(-Comment, - Time) %>% 
  separate_rows(KEGG_ID, sep = '; ') %>% 
  mutate(Comment = paste0(Comment, ' - ns')) %>% 
  group_by(Time, Comment) %>% 
  distinct() %>% 
  ungroup()

lc_kegg_unlabeled <- unlabeled_compounds_table %>%
  mutate(Comment = 'Unlabeled') %>% 
  select(FeatureID, Name, Formula, Time, Comment) %>% 
  left_join(kegg_table, by = c('Name', 'Formula')) %>% 
  distinct()  %>% 
  select(KEGG_ID, Time, Comment, Pathway_Name) %>% 
  separate_rows(KEGG_ID, sep = '; ')

  
nmr_kegg_ids <- rbind(T1_NMR.kegg, T2_NMR.kegg, T3_NMR.kegg) %>% 
  filter(!is.na(KEGG_ID)) %>% 
  select(KEGG_ID, Time, Comment, Pathway_Name) %>% 
  #mutate(Source = paste(Time, Comment, sep = '_')) %>% 
  #select(-Comment, - Time) %>% 
  mutate(Comment = paste0(Comment)) %>% 
  group_by(Time, Comment) %>% 
  distinct() %>% 
  ungroup()

unlabeled_nmr_ids <- unlabeled_nmr_kegg %>%
  separate_rows(KEGG_ID, sep = '; ') %>% 
  select(KEGG_ID, Time, Comment, Pathway_Name) %>% 
  group_by(Time, Comment) %>% 
  distinct() %>% 
  ungroup()

for(path in sig_pathways){
  all_kegg_ids <- rbind(lc_kegg_ids, lc_kegg_ns, lc_kegg_unlabeled, nmr_kegg_ids, unlabeled_nmr_ids) %>% 
    filter(Pathway_Name == path) %>%
    select(-Pathway_Name) %>% 
    distinct() %>% 
    mutate(Color = case_when(Comment == 'Not present in control' ~ 'red',
                             Comment == 'Upregulated' ~ 'orange',
                             Comment == 'Downregulated' ~ 'blue',
                             Comment == 'Only present in control' ~ 'purple',
                             # non significant keggs
                             Comment == 'Not present in control - ns' ~ 'pink',
                             Comment == 'Upregulated - ns' ~ 'pink',
                             Comment == 'Downregulated - ns' ~ '#0099ff',
                             Comment == 'Only present in control - ns' ~ '#0099ff',
                             # unlabeled samples
                             Comment == 'Unlabeled' ~ '#66cdaa')) %>% 
    # mutate(Color = case_when(Source == 'T1_Not present in control' ~ '#FF3300',
    #                          Source == 'T1_Upregulated' ~ '#FF6600',
    #                          Source == 'T1_Downregulated' ~ '#FF9900',
    #                          Source == 'T1_Only present in control' ~ '#FFCC00',
    #                          Source == 'T2_Not present in control' ~ '#000066',
    #                          Source == 'T2_Upregulated' ~ '#000099',
    #                          Source == 'T2_Downregulated' ~ '#3333CC',
    #                          Source == 'T2_Only present in control' ~ '#6969FF',
    #                          Source == 'T3_Not present in control' ~ '#005C2E',
    #                          Source == 'T3_Upregulated' ~ '#3329A3',
    #                          Source == 'T3_Downregulated' ~ '#66CC99',
    #                          Source == 'T3_Only present in control' ~ '#99FFCC')) %>% 
    mutate(for_kegg_mapper = paste(KEGG_ID, Color, sep = ' '))
  
  table_file <- file.path(tables_dir, paste0('KEGG_ids_for_mapper_', str_replace(path, '/', '-'), '.csv'))
  write_csv(all_kegg_ids, table_file)
  print(paste0(path, ' has finished'))
}

```