NK fatty acid validation

library(tidyverse)

## ── Attaching packages ────────────────────────────────── tidyverse 1.3.0 ──

## ✔ ggplot2 3.3.0     ✔ purrr   0.3.3
## ✔ tibble  2.1.3     ✔ dplyr   0.8.5
## ✔ tidyr   1.0.2     ✔ stringr 1.4.0
## ✔ readr   1.3.1     ✔ forcats 0.4.0

## ── Conflicts ───────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()

library(viridis)

## Loading required package: viridisLite

library(fgsea)

## Loading required package: Rcpp

library(cowplot)

## 
## ********************************************************

## Note: As of version 1.0.0, cowplot does not change the

##   default ggplot2 theme anymore. To recover the previous

##   behavior, execute:
##   theme_set(theme_cowplot())

## ********************************************************

FIG_fgsea_OUT_PATH <- "plots/validations/SchulteSchrepping/nk_fatty_acid/nk_fatty_acid_cohorts_1and2_severe-mild.pdf"
dir.create(dirname(FIG_fgsea_OUT_PATH), recursive = TRUE)

GENESETS_IN_PATH <- "input/genesets/kegg_go_btm_reactome_foointerferon.rds"

genesets <- readRDS(GENESETS_IN_PATH)

COHORT1_PATH <- "data/externalData/SchulteSchrepping/differential_expression/sample_groups/t0_covid_only/results/severity/fgsea_tables/severe-mild/id.celltype/15_NK cells--model@severity--coef@severe-mild--fgsea.tsv"
COHORT2_PATH <- "data/externalData/SchulteSchrepping/differential_expression_cohort2/sample_groups/t0_covid_only/results/severity/fgsea_tables/severe-mild/cluster_labels_res.0.4/NK cells--model@severity--coef@severe-mild--fgsea.tsv"

files <- list(cohort1 = COHORT1_PATH, cohort2 = COHORT2_PATH)

fgsea_list <- lapply(files, function(path){
  read_tsv(path,
           col_types = cols(
                       pathway = col_character(),
                       pval = col_double(),
                       padj = col_double(),
                       ES = col_double(),
                       NES = col_double(),
                       nMoreExtreme = col_double(),
                       size = col_double(),
                       leadingEdge = col_character()
                       ))
})
#names(fgsea_list) <- sapply(strsplit(files, "\\/"), `[[`, 7)


COHORT1_TOPTAB_PATH <- "data/externalData/SchulteSchrepping/differential_expression/sample_groups/t0_covid_only/results/severity/toptables/severe-mild/id.celltype/15_NK cells--model@severity--coef@severe-mild--toptab.tsv"
COHORT2_TOPTAB_PATH <- "data/externalData/SchulteSchrepping/differential_expression_cohort2/sample_groups/t0_covid_only/results/severity/toptables/severe-mild/cluster_labels_res.0.4/NK cells--model@severity--coef@severe-mild--toptab.tsv"

toptab_files <- list(cohort1 = COHORT1_TOPTAB_PATH, cohort2 = COHORT2_TOPTAB_PATH)

toptab_list <- lapply(toptab_files, function(path){
                              read_tsv(path)

})

## Parsed with column specification:
## cols(
##   gene = col_character(),
##   logFC = col_double(),
##   AveExpr = col_double(),
##   t = col_double(),
##   P.Value = col_double(),
##   adj.P.Val = col_double(),
##   B = col_double()
## )
## Parsed with column specification:
## cols(
##   gene = col_character(),
##   logFC = col_double(),
##   AveExpr = col_double(),
##   t = col_double(),
##   P.Value = col_double(),
##   adj.P.Val = col_double(),
##   B = col_double()
## )

keep_pathways <- list("reactome_Fatty acid metabolism")

plot_list <- list()
for(cohort in c("cohort1", "cohort2")){
  for(path in unique(keep_pathways)){
    pval <- fgsea_list[[cohort]] %>% filter(pathway == path) %>% pull(pval) %>%
            round(5)
    ranks <- toptab_list[[cohort] ]%>%
            select(gene, t) %>%
            deframe()
    plot_list[[paste(cohort, path, sep = "_")]] <- 
            plotEnrichment(stats = ranks, pathway = genesets[[path]]) +
            ylim(c(-.6, .6)) +
            ggtitle(paste(cohort, "Severe - Mild", path, sep = "\n")) +
             annotate("text",  x=Inf, y = Inf, label = paste("p =", pval), vjust=1, hjust=1)
  }

}
#combined_dat %>% filter(pathway == path & coef == coeff.) %>% as.data.frame()
p <- plot_grid(plotlist = plot_list, ncol = 2)
print(p)

#ggsave(plot = p, filename = FIG_fgsea_OUT_PATH, width = 8, height =4)