Kim Lab: IDD v IAD (Added GSEA)

1 Libraries Setup

library("hpgltools")
library("ggplot2")
library("patchwork")
library("plotly")
library("RColorBrewer")
library("DT")
library("dplyr")
library("gprofiler2")
library("DESeq2")
source("helper_functions.R")

1.1 Changelog

• hard-coded the year/month for downloading annotations (atb)
• removed Theresa’s notes as per Najib’s request
• minor formatting changes while I read

1.2 Updates

• CL versus Saline after adjusting for diet effects
• Added IDD_Saline - IAD_Saline / IDD_CL - IAD_CL contrast
  
• Added IDD_Saline-IDD_CL / IAD_Saline-IAD_CL contrast
  

2 M. musculus

2.1 Annotations

I am using mm38_100.

mm_annot <- load_biomart_annotations(species="mmusculus", year="2023", month="07")

## Using mart: ENSEMBL_MART_ENSEMBL from host: jul2023.archive.ensembl.org.

## Successfully connected to the mmusculus_gene_ensembl database.

## Finished downloading ensembl gene annotations.

## Finished downloading ensembl structure annotations.

## symbol columns is null, pattern matching 'symbol'.

## Including symbols, there are 57550 vs the 149547 gene annotations.

## Not dropping haplotype chromosome annotations, set drop_haplotypes = TRUE if this is bad.

## Saving annotations to mmusculus_biomart_annotations.rda.

## Finished save().

mm_annot <- mm_annot[["annotation"]]
mm_annot[["txid"]] <- paste0(mm_annot[["ensembl_transcript_id"]], ".", mm_annot[["version"]])
rownames(mm_annot) <- make.names(mm_annot[["ensembl_gene_id"]], unique=TRUE)
tx_gene_map <- mm_annot[, c("txid", "ensembl_gene_id")]

So, I now have a table of mouse annotations.

2.2 Hisat2 expressionset

The sample sheet is called ‘Experimental_design_Kim_v1.0.xlsx’ and has a column to point to the names of the count tables to load. Here I combine the metadata, count data, and annotations.

hisat_annot <- mm_annot
##rownames(hisat_annot) <- paste0("gene:", rownames(hisat_annot))
mm38_hisat <- create_expt("sample_sheets/Experimental_design_Kim_v1.0.xlsx",
                          gene_info=hisat_annot)

## Reading the sample metadata.

## Did not find the column: sampleid.

## Setting the ID column to the first column.

## The sample definitions comprises: 16 rows(samples) and 23 columns(metadata fields).

## Matched 25583 annotations and counts.

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## Saving the expressionset to 'expt.rda'.

## The final expressionset has 25760 features and 16 samples.

color_condition <- mm38_hisat$initial_metadata$condition
color_condition <- factor(color_condition, levels = c("IDD-Saline", "IDD-CL", "IAD-Saline", "IAD-CL"))

color_values <- c("#F8766D", "#FBADA7",  "#7CAE00", "#D3EE8F")#, "#00BFC4",  "#9AD2D3",  )

3 Summary Plots

plot_libsize(mm38_hisat)

plot_nonzero(mm38_hisat)

##  [1] "5214_S1"   "3740_S3"   "3750_S4"   "3772_S5"   "3774_S6"   "3775_S7"  
##  [7] "3766_S9"   "3767_S10"  "3742_S12"  "3747_S13"  "3741b_S14"

4 Normalize Expression

First we will filter out the genes which have low counts across all samples. Then I will normalize without background correction (it was determined there is no need for adjusting for background noise after evaluating the SVA correction).

mm38_filt <- normalize_expt(mm38_hisat, filter=TRUE)

## Removing 13318 low-count genes (12442 remaining).

mm38_norm <- normalize_expt(mm38_filt, convert="cpm", norm = "quant", transform="log2")

## transform_counts: Found 2 values equal to 0, adding 1 to the matrix.

mm38_nb <- normalize_expt(mm38_hisat, filter = TRUE, convert = "cpm",
                          transform = "log2", batch = "svaseq")

## Removing 13318 low-count genes (12442 remaining).

## Setting 128 low elements to zero.

## transform_counts: Found 128 values equal to 0, adding 1 to the matrix.

5 PCA

pca_norm <- plot_pca(mm38_norm, plot_labels = FALSE)
pca_norm$plot +
    scale_color_manual(values = color_values,
                    aesthetics = "color") +
    scale_fill_manual(values = color_values,
                    aesthetics = "fill") +
  guides(fill = "none")

plot_pca(mm38_nb)

6 Tables for Number of Samples

How many samples do we have currently for each group?

table(pData(mm38_hisat)[,c( "diet", "treatment")])

##      treatment
## diet  CL Saline (Vehicle)
##   IAD  4                4
##   IDD  5                3

7 Differential Expression Analysis

In the conditions column, we have IDD-Saline, IDD-CL, IAD-Saline, and IAD-CL. For the differential expression analysis, we will perform the following contrasts:

Normal Pairwise Contrasts:
IDD-CL vs IDD-Saline
IAD-CL vs IAD-Saline
IDD-CL vs IAD-CL
IDD-Saline vs IAD-Saline

Contrasts after adjustment:
CL versus Saline after adjusting for diet
IDD versus IAD after adjusting for CL/Saline treatment

I will conduct these contrasts below and provide a volcano plot (both interactive and static), as well as a searchable table of the significant DE results for each contrast.

mm_de_normal <- all_pairwise(mm38_filt, do_ebseq = FALSE,
                             model_batch = "svaseq", parallel = FALSE)

## 
## IDD-Saline     IDD-CL IAD-Saline     IAD-CL 
##          3          5          4          4

## Removing 0 low-count genes (12442 remaining).

## Setting 128 low elements to zero.

## transform_counts: Found 128 values equal to 0, adding 1 to the matrix.

## Starting basic pairwise comparison.

## Basic step 0/3: Normalizing data.

## Basic step 0/3: Converting data.

## Basic step 0/3: Transforming data.

## Basic step 1/3: Creating mean and variance tables.

## Basic step 2/3: Performing 10 comparisons.

## Basic step 3/3: Creating faux DE Tables.

## Basic: Returning tables.

## The data should be suitable for EdgeR/DESeq/EBSeq.
## If they freak out, check the state of the count table
## and ensure that it is in integer counts.

## converting counts to integer mode

## gene-wise dispersion estimates

## mean-dispersion relationship

## final dispersion estimates

## Starting edgeR pairwise comparisons.

## The data should be suitable for EdgeR/DESeq/EBSeq.
## If they freak out, check the state of the count table
## and ensure that it is in integer counts.

## EdgeR step 1/9: Importing and normalizing data.

## EdgeR step 2/9: Estimating the common dispersion.

## EdgeR step 3/9: Estimating dispersion across genes.

## EdgeR step 4/9: Estimating GLM Common dispersion.

## EdgeR step 5/9: Estimating GLM Trended dispersion.

## EdgeR step 6/9: Estimating GLM Tagged dispersion.

## EdgeR step 7/9: Running glmFit, switch to glmQLFit by changing the argument 'edger_test'.

## EdgeR step 8/9: Making pairwise contrasts.

## libsize was not specified, this parameter has profound effects on limma's result.

## Using the libsize from expt$best_libsize.

## Limma step 1/6: choosing model.

## Limma step 2/6: running limma::voom(), switch with the argument 'which_voom'.

## Using normalize.method = quantile for voom.

## Limma step 3/6: running lmFit with method: ls.

## Limma step 4/6: making and fitting contrasts with no intercept. (~ 0 + factors)

## Finished make_pairwise_contrasts.

## Limma step 5/6: Running eBayes with robust = FALSE and trend = FALSE.

## Limma step 6/6: Writing limma outputs.

## Limma step 6/6: 1/6: Creating table: IADSaline_vs_IADCL.  Adjust = BH

## Limma step 6/6: 2/6: Creating table: IDDCL_vs_IADCL.  Adjust = BH

## Limma step 6/6: 3/6: Creating table: IDDSaline_vs_IADCL.  Adjust = BH

## Limma step 6/6: 4/6: Creating table: IDDCL_vs_IADSaline.  Adjust = BH

## Limma step 6/6: 5/6: Creating table: IDDSaline_vs_IADSaline.  Adjust = BH

## Limma step 6/6: 6/6: Creating table: IDDSaline_vs_IDDCL.  Adjust = BH

## Limma step 6/6: 1/4: Creating table: IADCL.  Adjust = BH

## Limma step 6/6: 2/4: Creating table: IADSaline.  Adjust = BH

## Limma step 6/6: 3/4: Creating table: IDDCL.  Adjust = BH

## Limma step 6/6: 4/4: Creating table: IDDSaline.  Adjust = BH

## Starting noiseq pairwise comparisons.

## The data should be suitable for EdgeR/DESeq/EBSeq.
## If they freak out, check the state of the count table
## and ensure that it is in integer counts.

mm_de_normal

## A pairwise differential expression with results from: basic, deseq, edger, limma, noiseq.

## This used a surrogate/batch estimate from: svaseq.

## The primary analysis performed 10 comparisons.

keepers <- list("IDDCL_vs_IADCL" = c("IDDCL", "IADCL"),
                "IDDSaline_vs_IADSaline" = c("IDDSaline", "IADSaline"),
                "IDDSaline_vs_IDDCL" = c("IDDSaline", "IDDCL"),
                "IADSaline_vs_IADCL" = c("IADSaline", "IADCL"))
mm_de_tables <- combine_de_tables(mm_de_normal, excel="excel/DE_20221003.xlsx",
                                  keepers = keepers)
mm_de_tables

## A set of combined differential expression results.

##                    table deseq_sigup deseq_sigdown edger_sigup edger_sigdown
## 1         IDDCL_vs_IADCL         132           115         142           117
## 2 IDDSaline_vs_IADSaline          48            55          53            58
## 3     IDDSaline_vs_IDDCL           4            27           5            40
## 4     IADSaline_vs_IADCL          12             5          12             7
##   limma_sigup limma_sigdown
## 1          84            71
## 2          15             7
## 3           0             0
## 4           1             0

## Plot describing unique/shared genes in a differential expression table.

mm_de_sig <- extract_significant_genes(mm_de_tables,
                                       excel="excel/DE_sig_20221003.xlsx",
                                       according_to = "deseq")
mm_de_sig

## A set of genes deemed significant according to deseq.

## The parameters defining significant were:

## LFC cutoff: 1 adj P cutoff: 0.05

##                        deseq_up deseq_down
## IDDCL_vs_IADCL              132        115
## IDDSaline_vs_IADSaline       48         55
## IDDSaline_vs_IDDCL            4         27
## IADSaline_vs_IADCL           12          5

res_tbls <- c()
for (t in seq_along(keepers)) {
  tbl <- names(keepers)[t]
  deseq_tbl <- mm_de_normal$deseq$all_tables[[tbl]]
  res_tbls[[tbl]] <- merge(deseq_tbl, hisat_annot, by = "row.names")
  res_tbls[[tbl]] <- set_sig_limma(res_tbls[[tbl]],
                                   factors = c(paste(gsub("\\_.*","",tbl), "\nEnriched"),
                                               paste(sub('.*\\_', '', tbl), "\nEnriched")))
  res_tbls[[paste0(tbl, "volc_plotly")]] <- volc_plot(res_tbls[[tbl]],
                                                      title = paste("Volcano Plot:", tbl),
                                                      type = "plotly",
                                                      tbl = "limma",
                                                      gene_name_col = "external_gene_name")
  res_tbls[[paste0(tbl, "volc_ggplot")]] <- volc_plot(res_tbls[[tbl]],
                                                      title = paste("Volcano Plot:", tbl),
                                                      tbl = "limma")
}

7.1 IDD: CL vs Saline

Volcano Plot

res_tbls$IDDSaline_vs_IDDCLvolc_plotly

Table of Significant DE results

This table is significant results with adjusted p-value < 0.05 and log2FC greater than .5.

res_tbls$IDDSaline_vs_IDDCL %>%
  filter(abs(logFC) > .5 & adj.P.Val < 0.05) %>%
  arrange(desc(logFC)) %>%
  select("ensembl_gene_id", "mgi_symbol", "description", "logFC", "adj.P.Val", "Significance") %>%
  DT::datatable(rownames = FALSE)

GSEA

IDD_genes <- res_tbls$IDDSaline_vs_IDDCL %>%
  filter(abs(logFC) >= .5 & adj.P.Val <= 0.05) %>%
    arrange(desc(abs(logFC)))
IDD_gost <- gost(query = IDD_genes$ensembl_gene_id,
                organism = "mmusculus", ordered_query = TRUE,
                evcodes = TRUE)

#IDD Saline Enriched
IDD_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  datatable(rownames = FALSE, caption = "GO: IDD Significant GSEA")

IDD_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "REAC") %>%
  datatable(rownames = FALSE, caption = "REAC: IDD Significant GSEA")

IDD_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "KEGG") %>%
  datatable(rownames = FALSE, caption = "KEGG: IDD Significant GSEA")

##### IDD Enriched
gostplot(IDD_gost, capped = FALSE, interactive = TRUE)

7.2 IAD: CL vs Saline

Volcano Plot

res_tbls$IADSaline_vs_IADCLvolc_plotly

Table of Significant DE results

res_tbls$IADSaline_vs_IADCL %>%
  filter(abs(logFC) > .5 & adj.P.Val < 0.05) %>%
  arrange(desc(logFC)) %>%
  select("ensembl_gene_id", "mgi_symbol", "description", "logFC", "adj.P.Val", "Significance") %>%
  DT::datatable(rownames = FALSE)

GSEA

IAD_genes <- res_tbls$IADSaline_vs_IADCL %>%
  filter(abs(logFC) >= .5 & adj.P.Val <= 0.05) %>%
    arrange(desc(abs(logFC)))
IAD_gost <- gost(query = IAD_genes$ensembl_gene_id,
                organism = "mmusculus", ordered_query = TRUE,
                evcodes = TRUE)

#IAD Saline Enriched
IAD_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  datatable(rownames = FALSE, caption = "GO: IAD Significant GSEA")

##### IAD Enriched
gostplot(IAD_gost, capped = FALSE, interactive = TRUE)

7.3 CL: IDD vs IAD

Volcano Plot

res_tbls$IDDCL_vs_IADCLvolc_plotly

Table of Significant DE results

res_tbls$IDDCL_vs_IADCL %>%
  filter(abs(logFC) > .5 & adj.P.Val < 0.05) %>%
  arrange(desc(logFC)) %>%
  select("ensembl_gene_id", "mgi_symbol", "description", "logFC", "adj.P.Val", "Significance") %>%
  DT::datatable(rownames = FALSE)

GSEA

CL_genes <- res_tbls$IDDCL_vs_IADCL %>%
  filter(abs(logFC) >= .5 & adj.P.Val <= 0.05) %>%
    arrange(desc(abs(logFC)))
CL_gost <- gost(query = CL_genes$ensembl_gene_id,
                organism = "mmusculus", ordered_query = TRUE,
                evcodes = TRUE)

#CL Saline Enriched
CL_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  datatable(rownames = FALSE, caption = "GO: CL Significant GSEA")

CL_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "REAC") %>%
  datatable(rownames = FALSE, caption = "REAC: CL Significant GSEA")

CL_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "KEGG") %>%
  datatable(rownames = FALSE, caption = "KEGG: CL Significant GSEA")

##### CL Enriched
gostplot(CL_gost, capped = FALSE, interactive = TRUE)

7.4 Saline: IDD vs IAD

Volcano Plot

res_tbls$IDDSaline_vs_IADSalinevolc_plotly

Table of Significant DE results

res_tbls$IDDSaline_vs_IADSaline %>%
  filter(abs(logFC) > 1 & adj.P.Val < 0.05) %>%
  arrange(desc(logFC)) %>%
  select("ensembl_gene_id", "mgi_symbol", "description", "logFC", "adj.P.Val", "Significance") %>%
  DT::datatable(rownames = FALSE)

GSEA

Saline_genes <- res_tbls$IDDSaline_vs_IADSaline %>%
  filter(abs(logFC) >= .5 & adj.P.Val <= 0.05) %>%
    arrange(desc(abs(logFC)))
Saline_gost <- gost(query = Saline_genes$ensembl_gene_id,
                organism = "mmusculus", ordered_query = TRUE,
                evcodes = TRUE)

#CL Saline Enriched
Saline_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  datatable(rownames = FALSE, caption = "GO: Saline Significant GSEA")

Saline_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "REAC") %>%
  datatable(rownames = FALSE, caption = "REAC: Saline Significant GSEA")

Saline_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "KEGG") %>%
  datatable(rownames = FALSE, caption = "KEGG: Saline Significant GSEA")

##### CL Enriched
gostplot(Saline_gost, capped = FALSE, interactive = TRUE)

8 Comparison Between Contrasts

8.1 CL versus Saline after adjusting for diet

countdata <- exprs(mm38_filt)
coldata <- colData(mm38_filt)
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData = coldata,
                                 design = ~ diet + treatment)

## converting counts to integer mode

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

dds <- DESeq(ddsMat)

## estimating size factors
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## final dispersion estimates

## fitting model and testing

res <- results(dds)
res <- merge(as.data.frame(res), hisat_annot, by = 0)
res <- res[!is.na(res$padj),]
res <- set_sig(res, factors = c("Saline \n Enriched", "CL \n Enriched"))
volc_plot(res_tbl = res, type = "plotly", title = "Saline versus CL after adjusting for Diet")

res %>%
  filter(padj < 0.05 & abs(log2FoldChange) >= 1) %>%
  select("ensembl_gene_id", "mgi_symbol", "description", "log2FoldChange", "padj", "Significance") %>%
  arrange(desc(log2FoldChange)) %>%
  datatable(rownames = FALSE)

GSEA

adj_genes <- res %>%
  filter(padj < 0.05 & abs(log2FoldChange) >= 1) %>%
    arrange(desc(abs(log2FoldChange)))
adj_gost <- gost(query = adj_genes$ensembl_gene_id,
                organism = "mmusculus", ordered_query = TRUE,
                evcodes = TRUE)

adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  datatable(rownames = FALSE, caption = "GO: Adjusted for Diet Significant GSEA")

adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "REAC") %>%
  datatable(rownames = FALSE, caption = "REAC: Adjusted for Diet Significant GSEA")

adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "KEGG") %>%
  datatable(rownames = FALSE, caption = "KEGG: Adjusted for Diet Significant GSEA")

##### CL Enriched
gostplot(adj_gost, capped = FALSE, interactive = TRUE)

library(openxlsx)
res_df <- res %>%
  filter(padj < 0.05 & abs(log2FoldChange) >= .5) %>%
    arrange(desc(log2FoldChange)) %>%
    select(-Row.names, -lfcSE, -stat, -pvalue, -ensembl_gene_id, -version, -transcript_version, -mgi_symbol)
openxlsx::write.xlsx(res_df, file = "excel/adjdiet_DE.xlsx")

8.2 IDD versus IAD after adjusting for treatment

countdata <- exprs(mm38_filt)
coldata <- colData(mm38_filt)
ddsMat <- DESeqDataSetFromMatrix(countData = countdata,
                                 colData = coldata,
                                 design = ~ treatment + diet)

## converting counts to integer mode

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

dds <- DESeq(ddsMat)

## estimating size factors
##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## estimating dispersions

## gene-wise dispersion estimates

## mean-dispersion relationship

##   Note: levels of factors in the design contain characters other than
##   letters, numbers, '_' and '.'. It is recommended (but not required) to use
##   only letters, numbers, and delimiters '_' or '.', as these are safe characters
##   for column names in R. [This is a message, not a warning or an error]

## final dispersion estimates

## fitting model and testing

res <- results(dds)
res <- merge(as.data.frame(res), hisat_annot, by = 0)
res <- res[!is.na(res$padj),]
res <- set_sig(res, factors = c("IDD \n Enriched", "IAD \n Enriched"))
volc_plot(res_tbl = res, type = "plotly", title = "IDD versus IAD after adjusting for Treatment")

res %>%
  filter(padj < 0.05 & abs(log2FoldChange) >= 1) %>%
  select("ensembl_gene_id", "mgi_symbol", "description", "log2FoldChange", "padj", "Significance") %>%
  arrange(desc(log2FoldChange)) %>%
  datatable(rownames = FALSE)

GSEA

adj_genes <- res %>%
  filter(padj < 0.05 & abs(log2FoldChange) >= 1) %>%
    arrange(desc(abs(log2FoldChange)))
adj_gost <- gost(query = adj_genes$ensembl_gene_id,
                organism = "mmusculus",
                ordered_query = TRUE,
                evcodes = TRUE)

adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  datatable(rownames = FALSE, caption = "GO: Diet after Adjusted for Treatment Significant GSEA")

adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "REAC") %>%
  datatable(rownames = FALSE, caption = "REAC: Diet after Adjusted for Treatment Significant GSEA")

adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "KEGG") %>%
  datatable(rownames = FALSE, caption = "KEGG: Diet after Adjusted for Treatment Significant GSEA")

##### CL Enriched
gostplot(adj_gost, capped = FALSE, interactive = TRUE)

gseaplot_GO_diet_adjtreat <- adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "GO:BP" | source == "GO:MF" | source == "GO:CC") %>%
  head(n = 10) %>%
    ggplot(aes(x = recall, y = reorder(term_name, recall), fill = p_value)) +
  geom_bar(stat = "identity")+
  scale_fill_continuous(low = "blue", high = "red") +
  theme_bw()+
  ylab("") +
  xlab("GSEA Score")

gseaplot_REAC_diet_adjtreat <- adj_gost$result %>%
  select(term_name, p_value, term_size, intersection_size, recall, source, intersection) %>%
  arrange(desc(recall)) %>%
  filter(source == "REAC") %>%
  head(n = 10) %>%
    ggplot(aes(x = recall, y = reorder(term_name, recall), fill = p_value)) +
  geom_bar(stat = "identity")+
  scale_fill_continuous(low = "blue", high = "red") +
  theme_bw()+
  ylab("") +
  xlab("GSEA Score")

gseaplot_GO_diet_adjtreat

gseaplot_REAC_diet_adjtreat

library(openxlsx)
res_df <- res %>%
  filter(padj < 0.05 & abs(log2FoldChange) >= .5) %>%
    arrange(desc(log2FoldChange)) %>%
    select(-Row.names, -lfcSE, -stat, -pvalue, -ensembl_gene_id, -version, -transcript_version, -hgnc_symbol)
openxlsx::write.xlsx(res_df, file = "excel/adjtreatment_DE.xlsx")