4.P639-1_Gene_Analysis.qmd

---
title: "P639-1: Cheetah RNA-seq analysis Nasal Lavage"
author: "Basilin Benson, basilinb@uw.edu"
date: "version `r format(Sys.time(),'%B %d, %Y')`"
output: 
  html_document:
    toc: yes
    embed-resources: true
    toc_float:
      collapsed: no
  pdf_document:
    toc: yes
editor_options:
  chunk_output_type: console
header-includes:
  \usepackage{float}
---

```{r}
#| include = FALSE
knitr::opts_chunk$set(
  fig.height = 8.5,
  fig.width = 8.5,
  echo = FALSE,
  message = FALSE,
  warning = FALSE
)
```

# Background

The purpose of this workflow is to complete basic data cleaning of metadata and RNA-seq libraries. This includes 1) removing low coverage libraries, 2) filtering rare genes, 3) removing outlying libraries and duplicates, and 4) normalizing for RNA composition.

# Setup and Configuration

## Load Required Packages

```{r}
#| label = "setup",
#| message = FALSE,
#| warning = FALSE
##Data and figures
#library(MOFA2)
library(tidyverse)
library(BIGpicture)
#library(RNAetc)
library(ggdendro)
#Addtl figure customization
library(patchwork)
library(ggrepel)
library(scales)
#Correlation
library(psych)
library(corrplot)
library(cowplot)
#Linear modeling
library(limma)
library(kimma)
library(UpSetR)
#Rmd
library(kableExtra)
library(gginnards)
#rrvgo
library(rrvgo)
library(mgcv)
library(ComplexHeatmap)
#library(ica)
library(gtools)
library(readr)
library(smplot2)
library(ggpubr)
library(HoardR)
#db.format <- msigdbr::msigdbr("human", "C5")
#db.format_new <- db.format %>% select(gs_name,gs_exact_source) %>% rename(pathway = gs_name) %>% distinct()

`%notin%` <- Negate(`%in%`)
select <- dplyr::select
rename <- dplyr::rename
filter <- dplyr::filter
distinct <- dplyr::distinct
mutate <- dplyr::mutate
options(ggrepel.max.overlaps = Inf)
theme_classic_custom <- theme_bw() +
  theme(
    text = element_text(size = 10, color = "black"),
    axis.text = element_text(size = 10, color = "black"),
    panel.border = element_rect(color = "black", fill = NA, size = 1),
    axis.ticks = element_line(color = "black")
  )

visit_col_vec <- c(
  "Screening" = "#15134BFF",
  "Visit 1" = "#F56455FF",
  "Cold A Clinic" = "#87C785FF",
  "Cold B Home" = "#572F30FF"
)

#stim_col_vec <- c("Screening"="black","Visit 1" = "steelblue2","Cold A Clinic" = "goldenrod3","Cold B Home"="darkolivegreen","10Day"="darkred","No_Treatment" = "#E69F00","IL-13" = "#0072B2","No_Treatment_No_Infection" = "#CC79A7","IL-13_RV16" = "#fc4e2a","No_Treatment_RV16"="aquamarine","T2_High" = "darkred","T2_Low"="black","Severe_EX" = "darkred","No_Severe_EX"="black")

#IL_13_treatment_col <- c("No_Treatment_RV16" = "#00cb9e", "IL-13_No_Infection" = "orange", "IL-13_RV16" = "#7d25a9")
#IL_13_treatment_studyGroup_col <- group_colors <- c(
#  "No_Treatment_RV16_Healthy"     = "#66f2d1",  # lighter teal
#  "No_Treatment_RV16_Asthma"      = "#00cb9e",  # original teal (darker)
#
#   "IL-13_No_Infection_Healthy"    = "#ffd27f",  # light orange
#   "IL-13_No_Infection_Asthma"     = "orange",   # original orange (darker)
#
#   "IL-13_RV16_Healthy"            = "#b985d6",  # light purple
#   "IL-13_RV16_Asthma"             = "#7d25a9"   # original purple (darker)
# )
# break_time_w_control <- c(0,2,4,7,10)
# time_label_w_control <- c("Pre-Infection (Day 0)","Post-Infection (Day 2)","Post-Infection (Day 4)","Post-Infection (Day 7)","Post-Infection (Day 10)")
# shape_time <- c("48hr"=16,"96hr"=17,"7Day"=15,"10Day"=18)
# break_viral_load <- c(0,1,2,3,4,5)
# label_viral_load <- c("0",expression("10"^1),expression("10"^2),expression("10"^3),expression("10"^4),expression("10"^5))
```

## Load Visualization Functions

```{r}
# Scatterplot function for comparing two contrasts
source(
  "https://raw.githubusercontent.com/basilinb/Bioinformatics_R_scripts/refs/heads/main/make_contrast_scatterplot.R"
)
# Volcano plot function
source(
  "https://raw.githubusercontent.com/basilinb/Bioinformatics_R_scripts/refs/heads/main/Volcano_fxn.R"
)
# Plot enrichment comparisons function
source(
  "https://raw.githubusercontent.com/basilinb/Bioinformatics_R_scripts/refs/heads/main/plot_enrichment_comparisons.R"
)

# Function for color selector
source(
  "https://raw.githubusercontent.com/basilinb/Bioinformatics_R_scripts/refs/heads/main/make_color_palette.R"
)
```




## Directory Structure Configuration

```{r}
projectNumber_lavage <- "P639-1"
dirBoxBase <-
  file.path("~", "Library", "CloudStorage", "Box-Box")
dirRoot <-
  file.path(
    dirBoxBase,
    "Altman_Lab", # Box version
    # "~", "Documents", "Projects_local", # local version
    projectNumber_lavage
  )
dirPlots_lavage <- file.path("Figs", projectNumber_lavage)
subdirs_used <- c(
  "Volcano",
  "Contrast_Scatter",
  "Boxplot",
  "Enrichment",
  "Heatmap",
  "Density"
)
dirPlots_volcano_lavage <- file.path(dirPlots_lavage, "Volcano")
dirPlots_contrast_scatter_lavage <- file.path(
  dirPlots_lavage,
  "Contrast_Scatter"
)
dirPlots_boxplot_lavage <- file.path(
  dirPlots_lavage,
  "Boxplot"
)

dirPlots_enrichment_lavage <- file.path(
  dirPlots_lavage,
  "Enrichment"
)

dirPlots_heatmap_lavage <- file.path(
  dirPlots_lavage,
  "Heatmap"
)

dirresults_lavage <- file.path("Results", projectNumber_lavage)
dirDataclean <- file.path("Data_clean")
dirDataraw <- file.path("Data_raw")


if (!dir.exists(dirPlots_lavage)) {
  dir.create(dirPlots_lavage)
}
if (!dir.exists(dirPlots_volcano_lavage)) {
  dir.create(dirPlots_volcano_lavage)
}
if (!dir.exists(dirPlots_contrast_scatter_lavage)) {
  dir.create(dirPlots_contrast_scatter_lavage)
}
if (!dir.exists(dirPlots_boxplot_lavage)) {
  dir.create(dirPlots_boxplot_lavage)
}
if (!dir.exists(dirPlots_enrichment_lavage)) {
  dir.create(dirPlots_enrichment_lavage)
}
if (!dir.exists(dirPlots_heatmap_lavage)) {
  dir.create(dirPlots_heatmap_lavage)
}
if (!dir.exists(dirresults_lavage)) {
  dir.create(dirresults_lavage)
}

dirT2plots_lavage <- HoardR::init_experiment_dirs(
  parent_dir = dirPlots_lavage,
  experiment_name = "T2_Classification",
  plot_subdirs = subdirs_used
)
```

## Load the Muppits modules
```{r}
muppits_annotations <- read_csv("Data_raw/Muppits_data/module_annotation.csv")

muppits_mods <- read_csv(file.path(
  dirDataraw,
  "MUPPITS_geneSetsByAssignedCell.csv"
)) %>%
  left_join(muppits_annotations, by = c("geneSet" = "mod")) %>%
  mutate(
    annot = case_when(
      is.na(annot) ~ "Unannotated",
      TRUE ~ annot
    )
  ) %>%
  mutate(gs_name = paste0(annot, " (", geneSet, ")"))

muppits_mods_for_enrichment <- muppits_mods %>%
  select(gs_name, symbol) %>%
  rename(gene_symbol = symbol)


muppits_mods_expression_w_meta_all_lavage <- read_csv(file.path(
  dirDataclean,
  paste(projectNumber_lavage, "_Muppits_expression_data.csv", sep = "")
))
```

# Data Loading and Preprocessing

## Load Raw Data and Clinical Metadata
```{r}
asthma_status <- readxl::read_xlsx(file.path(dirDataraw, "Asthma_status.xlsx"))
clinical_data_list <- rio::import_list(file.path(dirDataraw, "Cold_Data.xlsx"))

file.tmp <- file.path(
  dirDataclean,
  paste(projectNumber_lavage, "_dat_voom_normalized_w_clinical.RDS", sep = "")
)
dat_voom_lavage <- readRDS(file = file.tmp)

dat_voom_lavage$targets <- dat_voom_lavage$targets %>%
  mutate(
    EOS_cat = case_when(
      `Eosinophil Count (V0)` < 150 ~ "EOS_Less_than_150",
      `Eosinophil Count (V0)` >= 150 &
        `Eosinophil Count (V0)` < 300 ~ "EOS_between_150_300",
      `Eosinophil Count (V0)` >= 300 ~ "EOS_Greater_than_300",
      TRUE ~ NA_character_ # safety fallback
    )
  ) %>%
  mutate(
    FENO_cat = case_when(
      `FeNO (V0)` < 20 ~ "FENO_Less_than_20",
      `FeNO (V0)` >= 20 & `FeNO (V0)` < 35 ~ "FENO_between_20_35",
      `FeNO (V0)` >= 35 ~ "FENO_Greater_than_35",
      TRUE ~ NA_character_ # safety fallback
    )
  ) %>%
  mutate(
    Allergy_Sensitization_cat = case_when(
      `Allergen Sensitizations per AST (V0)` <=
        1 ~ "Allergy_Sensitization_Less_than_1",
      `Allergen Sensitizations per AST (V0)` >= 2 &
        `Allergen Sensitizations per AST (V0)` <
          5 ~ "Allergy_Sensitization_between_2_5",
      `Allergen Sensitizations per AST (V0)` >=
        5 ~ "Allergy_Sensitization_Greater_than_5",
      TRUE ~ NA_character_ # safety fallback
    )
  ) %>%
  mutate(
    EOS_cat = factor(
      EOS_cat,
      levels = c(
        "EOS_Less_than_150",
        "EOS_between_150_300",
        "EOS_Greater_than_300"
      )
    )
  ) %>%
  mutate(
    FENO_cat = factor(
      FENO_cat,
      levels = c(
        "FENO_Less_than_20",
        "FENO_between_20_35",
        "FENO_Greater_than_35"
      )
    )
  ) %>%
  mutate(
    Allergy_Sensitization_cat = factor(
      Allergy_Sensitization_cat,
      levels = c(
        "Allergy_Sensitization_Less_than_1",
        "Allergy_Sensitization_between_2_5",
        "Allergy_Sensitization_Greater_than_5"
      )
    )
  ) %>%
  mutate(
    log10_EOS = log10(`Eosinophil Count (V0)` + 1),
    log10_FENO = log10(`FeNO (V0)` + 1)
  ) %>%
  mutate(
    asthma_status = ifelse(
      subjectId %in% unique(asthma_status$subjectID),
      "Non_Asthma",
      "Asthma"
    )
  )
```



## Data Cleaning and Variable Creation

### Remove Unknown Visits and Clean Visit Labels

```{r}
file.tmp <- file.path(
  dirDataclean,
  paste(
    projectNumber_lavage,
    "_dat_voom_normalized_wo_unknown_w_clinical.rds",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  libs_wo_unknown <- dat_voom_lavage$targets %>%
    subset(visit %notin% c("Unknown")) %>%
    pull(libid)
  dat_voom_lavage_wo_unknown <- dat_voom_lavage[, libs_wo_unknown]

  dat_voom_lavage_wo_unknown$targets <- dat_voom_lavage_wo_unknown$targets %>%
    mutate(
      visit = ifelse(visit == "Screening/Enrollment", "Screening", visit)
    ) %>%
    mutate(
      visit = factor(visit, levels = c("Screening", "Cold A Clinic", "Visit 1"))
    )

  saveRDS(dat_voom_lavage_wo_unknown, file.tmp)
} else {
  dat_voom_lavage_wo_unknown <- readRDS(file.tmp)
}
```

## Count Visit by T2
```{r}
T2_HIgh_Low_classification_by_eos1_meand_1sd_lavage <- read_csv(
  HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_lavage",
    ext = "csv"
  )
) %>%
  janitor::clean_names()

Lavage_samples_with_T2 <- dat_voom_lavage_wo_unknown$targets %>%
  janitor::clean_names() %>%
  left_join(
    T2_HIgh_Low_classification_by_eos1_meand_1sd_lavage %>%
      select(subject_id, gene_exp, t2_1sd)
  )

Lavage_samples_with_T2 %>% dplyr::count(visit, t2_1sd)
```

# Visit Differential Gene Expression Analysis

## Visit Comparison Analysis

### KIMMA Modeling for Visit Effects
```{r}
file.tmp <- file.path(
  dirresults_lavage,
  paste(
    projectNumber_lavage,
    "_visit_comparison_wo_unknown_kimma.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  model_visit_wo_unknown_lavage_kimma <- kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ visit + median_cv_coverage + (1|subjectId)",
    run_lme = TRUE,
    run_contrast = TRUE,
    use_weights = TRUE,
    metrics = TRUE,
    contrast_var = "visit",
    processors = 6
  )
  saveRDS(model_visit_wo_unknown_lavage_kimma, file = file.tmp)
} else {
  model_visit_wo_unknown_lavage_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_visit_wo_unknown_lavage_kimma$lme.contrast,
  FCgroup = TRUE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```



### Extract and Process Visit Comparison Results

```{r}
#| label = "visit_model_results"
file.tmp <- file.path(
  dirresults_lavage,
  paste(
    projectNumber_lavage,
    "_visit_comparison_wo_unknown_kimma_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  model_visit_wo_unknown_lavage_kimma_results <- model_visit_wo_unknown_lavage_kimma$lme.contrast %>%
    left_join(
      dat_voom_lavage_wo_unknown$genes %>% select(geneName, hgnc_symbol),
      by = c("gene" = "geneName")
    ) %>%
    mutate(contrast_ref = gsub(" ", "_", contrast_ref)) %>%
    mutate(contrast_lvl = gsub(" ", "_", contrast_lvl)) %>%
    mutate(comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_")) %>%
    select(
      gene,
      hgnc_symbol,
      comparison,
      estimate,
      pval,
      FDR,
      contrast_ref,
      contrast_lvl
    )
  write_csv(model_visit_wo_unknown_lavage_kimma_results, file.tmp)
} else {
  model_visit_wo_unknown_lavage_kimma_results <- read_csv(file.tmp)
}
```


## Visualization of Visit Comparison Results

### Volcano Plot for Cold A vs Screening Comparison

```{r}
nasal_lavage_visit_volc <- volcano_funcxn(
  model_visit_wo_unknown_lavage_kimma_results %>%
    mutate(
      volcano_comp = paste(contrast_ref, " <------> ", contrast_lvl, sep = "")
    ),
  fdr = 0.05,
  n_top_gene = 15
)


ggsave(
  file.path(dirPlots_lavage, "Volcano_plot_ColdA_vs_Screening.pdf"),
  plot = nasal_lavage_visit_volc$`Screening <------> Cold_A_Clinic`,
  height = 7,
  width = 10
)
```

### Generate Gene Lists for Enrichment Analysis

```{r}
fdr_cutoff <- 0.05

model_visit_wo_unknown_lavage_kimma_results_list <- model_visit_wo_unknown_lavage_kimma_results %>%
  filter(FDR <= fdr_cutoff) %>%
  group_by(comparison) %>%
  summarise(
    up_in_lvl = list(unique(hgnc_symbol[estimate > 0])),
    up_in_ref = list(unique(hgnc_symbol[estimate < 0])),
    lvl = first(contrast_lvl),
    ref = first(contrast_ref)
  ) %>%
  mutate(
    res = pmap(
      list(up_in_lvl, up_in_ref, lvl, ref, comparison),
      ~ set_names(
        list(..1, ..2),
        c(
          paste0("Up in ", ..3, " (", ..5, ")"),
          paste0("Up in ", ..4, " (", ..5, ")")
        )
      )
    )
  ) %>%
  pull(res) %>%
  flatten()
```

### UpSet Plot of Differential Genes

```{r}
upset(
  fromList(model_visit_wo_unknown_lavage_kimma_results_list),
  text.scale = 1.5,
  order.by = "freq",
  mainbar.y.label = "Significant genes\nFDR < 0.05",
  keep.order = TRUE,
  sets = names(model_visit_wo_unknown_lavage_kimma_results_list)
)
```

## Pathway Enrichment Analysis

### Run Enrichment Analysis

```{r}
model_visit_wo_unknown_lavage_kimma_results_list <- list()
model_visit_wo_unknown_lavage_kimma_results_list[[
  "Up in Cold A Clinic"
]] <- model_visit_wo_unknown_lavage_kimma_results %>%
  subset(comparison == "Cold_A_Clinic_vs_Screening") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
model_visit_wo_unknown_lavage_kimma_results_list[[
  "Up in Screening"
]] <- model_visit_wo_unknown_lavage_kimma_results %>%
  subset(comparison == "Cold_A_Clinic_vs_Screening") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list <- list()

    enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list[[
      "GO_BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = model_visit_wo_unknown_lavage_kimma_results_list,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = model_visit_wo_unknown_lavage_kimma_results_list,
      collection = "H"
    )

    enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = model_visit_wo_unknown_lavage_kimma_results_list,
      db = muppits_mods_for_enrichment
    )
    enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list
  },
  overwrite = FALSE
)


comparisons_to_plot <- list(
  "Group 1:Cold A Clinic vs Group 2:Screening " = c(
    "Up in Cold A Clinic",
    "Up in Screening"
  )
)

plot_H_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list$H,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  plot_H_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage
)

HoardR::save_figures(
  plot = plot_H_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "H_enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits",
    ext = "png"
  ),
  height = 4.5,
  width = 8,
  overwrite = FALSE
)

plot_GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list$GO_BP,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  plot_GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage
)

HoardR::save_figures(
  plot = plot_GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "GO_BP_enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits",
    ext = "png"
  ),
  height = 20,
  width = 15,
  overwrite = FALSE
)


plot_GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage_top_15_term <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list$GO_BP %>%
    group_by(group) %>%
    slice_min(order_by = FDR, n = 15) %>%
    slice_max(order_by = `k/K`, n = 15),
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  plot_GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage_top_15_term
)

HoardR::save_figures(
  plot = plot_GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage_top_15_term +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "GO_BP_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_top_15_term",
    ext = "png"
  ),
  height = 9,
  width = 10,
  overwrite = TRUE
)


plot_muppits_model_enrichment_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_list$muppits,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.05
) +
  theme_classic_custom
print(
  plot_muppits_model_enrichment_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage
)

HoardR::save_figures(
  plot = plot_muppits_model_enrichment_visit_wo_unknown_lavage_kimma_results_list_GO_H_muppits_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "muppits_model_enrichment_visit_wo_unknown_lavage_kimma_results",
    ext = "png"
  ),
  height = 4.5,
  width = 8,
  overwrite = TRUE
)
```

### Visualize Enrichment Results
```{r}
fdr_cutoff <- 0.2

#Hallmark
dat_enrichment_Visit_model_wo_unknown_lavage_H <- Enrichment_Visit_model_wo_unknown_lavage_H %>%
  mutate(
    sig = FDR <= fdr_cutoff,
    # extract text inside parentheses
    inside_parens = str_extract(group, "(?<=\\().*?(?=\\))"),
    # split at "_vs_"
    contrast_pair = str_split(inside_parens, "_vs_"),
    # assign contrast_lvl and contrast_ref correctly
    contrast_lvl = map_chr(contrast_pair, 1),
    contrast_ref = map_chr(contrast_pair, 2)
  ) %>%
  select(-inside_parens, -contrast_pair) %>%
  mutate(
    comparison = paste(
      "Contrast Lvl: ",
      contrast_lvl,
      "\nContrast Ref: ",
      contrast_ref,
      sep = ""
    )
  ) %>%
  group_by(comparison) %>%
  mutate(
    Significance = ifelse(
      sig & str_detect(group, paste0("^Up in ", contrast_lvl)),
      paste("Up in ", contrast_lvl, sep = ""),
      ifelse(
        sig & str_detect(group, paste0("^Up in ", contrast_ref)),
        paste("Down in ", contrast_lvl, sep = ""),
        NA
      )
    )
  )


plot_enrichment_Visit_model_wo_unknown_lavage_H <- ggplot(
  dat_enrichment_Visit_model_wo_unknown_lavage_H %>%
    subset(
      FDR <= fdr_cutoff &
        group %in%
          c(
            "Up in Cold_A_Clinic (Cold_A_Clinic_vs_Screening)",
            "Up in Screening (Cold_A_Clinic_vs_Screening)"
          )
    ),
  aes(x = comparison, y = pathway)
) +
  geom_point(
    aes(
      size = `k/K`,
      color = Significance,
      fill = Significance,
      shape = Significance
    ), # triangle
    alpha = 0.5
  ) +
  scale_shape_manual(
    values = c("Up in Cold_A_Clinic" = 24, "Down in Cold_A_Clinic" = 25)
  ) +
  scale_color_manual(
    values = c("Up in Cold_A_Clinic" = "red", "Down in Cold_A_Clinic" = "blue")
  ) +
  scale_fill_manual(
    values = c("Up in Cold_A_Clinic" = "red", "Down in Cold_A_Clinic" = "blue")
  ) +

  scale_size_continuous(range = c(1, 8)) + # adjust triangle sizes
  theme_classic_custom +
  theme(
    axis.title.x = element_blank()
  ) +
  labs(y = "Pathway", size = "k/K")

ggsave(
  file.path(dirPlots_lavage, "Enrichment_Visit_wo_unknown_lavage_Hallmark.pdf"),
  plot = plot_enrichment_Visit_model_wo_unknown_lavage_H,
  height = 5,
  width = 10
)

#GO
fdr_cutoff <- 0.01
dat_enrichment_Visit_model_wo_unknown_lavage_GO <- Enrichment_Visit_model_wo_unknown_lavage_GO %>%
  mutate(
    sig = FDR <= fdr_cutoff,
    # extract text inside parentheses
    inside_parens = str_extract(group, "(?<=\\().*?(?=\\))"),
    # split at "_vs_"
    contrast_pair = str_split(inside_parens, "_vs_"),
    # assign contrast_lvl and contrast_ref correctly
    contrast_lvl = map_chr(contrast_pair, 1),
    contrast_ref = map_chr(contrast_pair, 2)
  ) %>%
  select(-inside_parens, -contrast_pair) %>%
  mutate(
    comparison = paste(
      "Contrast Lvl: ",
      contrast_lvl,
      "\nContrast Ref: ",
      contrast_ref,
      sep = ""
    )
  ) %>%
  group_by(comparison) %>%
  mutate(
    Significance = ifelse(
      sig & str_detect(group, paste0("^Up in ", contrast_lvl)),
      paste("Up in ", contrast_lvl, sep = ""),
      ifelse(
        sig & str_detect(group, paste0("^Up in ", contrast_ref)),
        paste("Down in ", contrast_lvl, sep = ""),
        NA
      )
    )
  )


plot_enrichment_Visit_model_wo_unknown_lavage_GO <- ggplot(
  dat_enrichment_Visit_model_wo_unknown_lavage_GO %>%
    subset(
      FDR <= fdr_cutoff &
        group %in%
          c(
            "Up in Cold_A_Clinic (Cold_A_Clinic_vs_Screening)",
            "Up in Screening (Cold_A_Clinic_vs_Screening)"
          )
    ) %>%
    group_by(group) %>%
    slice_min(order_by = FDR, n = 10),
  aes(x = comparison, y = pathway)
) +
  geom_point(
    aes(
      size = `k/K`,
      color = Significance,
      fill = Significance,
      shape = Significance
    ),
    alpha = 0.5
  ) +
  scale_shape_manual(
    values = c("Up in Cold_A_Clinic" = 24, "Down in Cold_A_Clinic" = 25)
  ) +
  scale_color_manual(
    values = c("Up in Cold_A_Clinic" = "red", "Down in Cold_A_Clinic" = "blue")
  ) +
  scale_fill_manual(
    values = c("Up in Cold_A_Clinic" = "red", "Down in Cold_A_Clinic" = "blue")
  ) +
  scale_size_continuous(range = c(2, 8)) + # adjust triangle sizes
  theme_classic_custom +
  theme(
    axis.title.x = element_blank()
  ) +
  labs(y = "Pathway", size = "k/K")
ggsave(
  file.path(dirPlots_lavage, "Enrichment_Visit_wo_unknown_lavage_GO.pdf"),
  plot = plot_enrichment_Visit_model_wo_unknown_lavage_GO,
  height = 7,
  width = 9
)
```



### Generate Top Gene Expression Plots

```{r}
dat_lavage_wo_unknown_w_meta <- dat_voom_lavage_wo_unknown$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  pivot_longer(!libid, names_to = "geneName", values_to = "geneExp") %>%
  left_join(dat_voom_lavage_wo_unknown$targets) %>%
  left_join(dat_voom_lavage_wo_unknown$genes %>% select(geneName, hgnc_symbol))

top_4_cold_A_genes <- model_visit_wo_unknown_lavage_kimma_results %>%
  subset(
    hgnc_symbol %in%
      c("OAS3", "SOCS1", "IFI6", "IRF7", "CXCL10") &
      comparison == "Cold_A_Clinic_vs_Screening"
  )

box_plot_top_4_cold_A_genes <- dat_lavage_wo_unknown_w_meta %>%
  subset(geneName %in% top_4_cold_A_genes$gene) %>%
  subset(visit %in% c("Screening", "Cold A Clinic")) %>%
  ggplot(aes(x = visit, y = geneExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +

  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "red"
  ) +
  scale_color_manual(values = visit_col_vec) +
  theme_classic_custom +
  facet_wrap(~hgnc_symbol, scales = "free", nrow = 1) +
  ylab("Gene Expression (log2)")

ggsave(
  file.path(dirPlots_lavage, "box_plot_top_4_cold_A_genes.pdf"),
  plot = box_plot_top_4_cold_A_genes,
  height = 5,
  width = 12
)


top_4_screening_genes <- model_visit_wo_unknown_lavage_kimma_results %>%
  subset(
    comparison == "Cold_A_Clinic_vs_Screening" &
      hgnc_symbol %in% c("CXCL1", "FOS", "TMSB10", "ERAP1", "PKM")
  ) %>%
  arrange(FDR)

box_plot_top_4_screening_genes <- dat_lavage_wo_unknown_w_meta %>%
  subset(geneName %in% top_4_screening_genes$gene) %>%
  subset(visit %in% c("Screening", "Cold A Clinic")) %>%
  ggplot(aes(x = visit, y = geneExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +

  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "red"
  ) +
  scale_color_manual(values = visit_col_vec) +
  theme_classic_custom +
  facet_wrap(~hgnc_symbol, scales = "free", nrow = 1) +
  ylab("Gene Expression (log2)")

ggsave(
  file.path(dirPlots_lavage, "box_plot_top_4_screening_genes.pdf"),
  plot = box_plot_top_4_screening_genes,
  height = 5,
  width = 12
)


ggsave(
  file.path(
    dirPlots_lavage,
    "box_plot_top_4_screening_and_cold_A_genes_combined.pdf"
  ),
  plot = box_plot_top_4_cold_A_genes +
    box_plot_top_4_screening_genes +
    plot_layout(ncol = 1, guides = "collect"),
  height = 10,
  width = 12
)
```

### Generate STRING Network Diagrams

```{r}
dir.create(file.path(dirPlots_lavage, "String_Network"))
library(pullString)
for (mod_ in names(model_visit_wo_unknown_lavage_kimma_results_list)) {
  get_png_network(
    ids = model_visit_wo_unknown_lavage_kimma_results_list[[mod_]],
    png_filename = paste(
      file.path(dirPlots_lavage, "String_Network"),
      "/",
      mod_,
      ".png",
      sep = ""
    )
  )
}
```

# Viral Load Analysis

## Load and Process Viral Data

```{r}
viral_data_list <- rio::import_list(file.path(
  dirDataraw,
  "CHEETAH Nasal Sample Virology_working 101425.xlsx"
))

viral_data_positive_rna <- viral_data_list$Sheet1 %>%
  janitor::clean_names() %>%
  mutate(
    collection_date = suppressWarnings(as.numeric(collection_date)),
    collection_date = as.Date(collection_date, origin = "1899-12-30")
  ) %>%
  subset(!is.na(virus_1)) %>%
  mutate(Visit_start_date = collection_date) %>%
  rename("subjectId" = "subject_id") %>%
  mutate(barcode = as.character(barcode))


file.tmp <- file.path(
  dirDataclean,
  paste(
    projectNumber_lavage,
    "_dat_voom_normalized_wo_unknown_w_clinical_w_viral.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  dat_voom_lavage_wo_unknown_w_viral <- dat_voom_lavage_wo_unknown
  dat_voom_lavage_wo_unknown_w_viral$targets <- dat_voom_lavage_wo_unknown_w_viral$targets %>%
    mutate(barcode = as.character(barcode)) %>%
    rename("Visit_start_date" = "dateCollected") %>%
    left_join(
      viral_data_positive_rna %>%
        select(subjectId, Visit_start_date, visit, virus_1, virus_quantity_ct)
    ) %>%
    mutate(
      Virus_status = ifelse(
        is.na(virus_1),
        "No Viral Data",
        ifelse(virus_1 == "Negative", "Negative", "Positive")
      )
    )

  saveRDS(dat_voom_lavage_wo_unknown_w_viral, file.tmp)
} else {
  dat_voom_lavage_wo_unknown_w_viral <- readRDS(file.tmp)
}
```



## Viral Status Analysis at Cold A Clinic Visit

### KIMMA Modeling for Visit-Viral Status Interactions

```{r}
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_lavage,
    "_visit_viral_incomplete_comparison_wo_unknown_kimma.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  dat_voom_lavage_wo_unknown_w_viral$targets <- dat_voom_lavage_wo_unknown_w_viral$targets %>%
    mutate(
      Virus_status = ifelse(Virus_status == "No Viral Data", NA, Virus_status)
    )
  model_visit_viral_lavage_wo_unknown_kimma <- kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ visit*Virus_status + (1|subjectId)",
    run_lme = TRUE,
    run_contrast = TRUE,
    use_weights = TRUE,
    metrics = TRUE,
    contrast_var = "visit:Virus_status",
    processors = 6
  )
  saveRDS(model_visit_viral_lavage_wo_unknown_kimma, file = file.tmp)
} else {
  model_visit_viral_lavage_wo_unknown_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_visit_viral_lavage_wo_unknown_kimma$lme.contrast,
  FCgroup = TRUE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

model_visit_viral_lavage_wo_unknown_kimma_coldA_posvsneg <- model_visit_viral_lavage_wo_unknown_kimma$lme.contrast %>%
  subset(
    contrast_ref %in%
      c("Cold A Clinic Negative") &
      contrast_lvl %in% c("Cold A Clinic Positive")
  ) %>%
  left_join(
    dat_voom_wo_unknown_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  ) %>%
  mutate(comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_"))

kimma::summarize_kmFit(
  model_visit_viral_lavage_wo_unknown_kimma_coldA_posvsneg,
  FCgroup = TRUE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

Significant_genes_ColdA_posvsneg_lavage <- model_visit_viral_lavage_wo_unknown_kimma$lme.contrast %>%
  subset(
    contrast_ref == "Cold A Clinic Negative" &
      contrast_lvl == "Cold A Clinic Positive" &
      FDR <= 0.05
  ) %>%
  left_join(
    dat_voom_lavage_wo_unknown_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  )


Significant_genes_lavage_ColdA_posvsneg_list <- list()

Significant_genes_lavage_ColdA_posvsneg_list[[
  "Up in Viral Positive (Cold A)"
]] <- Significant_genes_ColdA_posvsneg_lavage %>%
  subset(estimate > 0) %>%
  pull(hgnc_symbol)
Significant_genes_lavage_ColdA_posvsneg_list[[
  "Down in Viral Positive (Cold A)"
]] <- Significant_genes_ColdA_posvsneg_lavage %>%
  subset(estimate < 0) %>%
  pull(hgnc_symbol)
```

### Viral Positive vs Negative Volcano Plot

```{r}
Volcano_plot_Significant_genes_ColdA_posvsneg_lavage <- volcano_funcxn(
  model_visit_viral_lavage_wo_unknown_kimma_coldA_posvsneg %>%
    subset(
      comparison %in% c("Cold A Clinic Positive_vs_Cold A Clinic Negative")
    ) %>%
    mutate(
      volcano_comp = "Cold A Clinic Negative<----->Cold A Clinic Positive"
    ),
  fdr = 0.05,
  n_top_gene = 10
)

Volcano_plot_Significant_genes_ColdA_posvsneg_lavage$`Cold A Clinic Negative<----->Cold A Clinic Positive` +
  theme_classic_custom
```

### Viral Status Enrichment Analysis

```{r}
enrichment_GO_Significant_genes_ColdA_posvsneg_lavage <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_lavage_ColdA_posvsneg_list,
  collection = "C5",
  subcollection = "GO:BP"
)

enrichment_H_Significant_genes_ColdA_posvsneg_lavage <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_lavage_ColdA_posvsneg_list,
  collection = "H"
)

BIGpicture::plot_enrich(
  enrichment_GO_Significant_genes_ColdA_posvsneg_lavage,
  fdr_cutoff = c(0.01)
)
BIGpicture::plot_enrich(enrichment_H_Significant_genes_ColdA_posvsneg_lavage)
```

## Extended Viral Analysis Including Screening Visit

### KIMMA Modeling with Visit-by-Viral Status Design
```{r}
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_lavage,
    "_visit_viral_incomplete_comparison_wo_unknown_kimma_with_Screening_and_coldB.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  dat_voom_lavage_wo_unknown_w_viral$targets <- dat_voom_lavage_wo_unknown_w_viral$targets %>%
    mutate(
      Virus_status = ifelse(Virus_status == "No Viral Data", NA, Virus_status)
    ) %>%
    mutate(
      visit_by_virus_status = ifelse(
        visit == "Screening",
        "Screening",
        ifelse(is.na(Virus_status), NA, paste(visit, Virus_status, sep = "_"))
      )
    )

  model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma <- kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ visit_by_virus_status + (1|subjectId)",
    run_lme = TRUE,
    run_contrast = TRUE,
    use_weights = TRUE,
    metrics = TRUE,
    contrast_var = "visit_by_virus_status",
    processors = 6
  )
  saveRDS(
    model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma,
    file = file.tmp
  )
} else {
  model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma <- readRDS(
    file.tmp
  )
}
kimma::summarize_kmFit(
  model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma$lme.contrast,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

kimma::summarize_kmFit(
  model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma$lme.contrast,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05)
) %>%
  select(-c(variable)) %>%
  filter(!is.na(contrast_ref), !is.na(contrast_lvl)) %>%
  pivot_wider(
    names_from = contrast_lvl,
    values_from = `fdr<0.05`,
    values_fill = 0
  ) %>%
  column_to_rownames("contrast_ref") %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma_dat <- model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma$lme.contrast %>%
  left_join(
    dat_voom_lavage_wo_unknown_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  ) %>%
  mutate(comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_"))
```



#### Make volcano plot
```{r}
Volcano_plot_Significant_Screening_and_coldB_lavage <- volcano_funcxn(
  model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
    subset(comparison %in% c("Screening_vs_Cold A Clinic_Positive")) %>%
    mutate(volcano_comp = "Cold A Clinic Positive<----->Screening"),
  fdr = 0.05,
  n_top_gene = 10
)

Volcano_plot_Significant_Screening_and_coldB_lavage$`Cold A Clinic Positive<----->Screening` +
  theme_classic_custom
```

#### Run Enrichment
```{r}
Significant_genes_Screening_and_coldB_list_lavage <- list()
Significant_genes_Screening_and_coldB_list_lavage[[
  "Up in Screening compared to Cold A Clinic"
]] <- model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
  subset(comparison %in% c("Screening_vs_Cold A Clinic_Positive")) %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_Screening_and_coldB_list_lavage[[
  "Up in Cold A Clinic compared to Screening"
]] <- model_visit_viral_lavage_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
  subset(comparison %in% c("Screening_vs_Cold A Clinic_Positive")) %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_GO_Significant_Screening_and_coldB_lavage <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_Screening_and_coldB_list_lavage,
  collection = "C5",
  subcollection = "GO:BP"
)

enrichment_H_Significant_genes_Screening_and_coldB_lavage <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_Screening_and_coldB_list_lavage,
  collection = "H"
)

BIGpicture::plot_enrich(
  enrichment_GO_Significant_Screening_and_coldB_lavage,
  fdr_cutoff = c(0.001)
)
BIGpicture::plot_enrich(
  enrichment_H_Significant_genes_Screening_and_coldB_lavage
)
```



# Screening Visit Analysis

## Subset Data to Screening Visit Only

```{r}
dat_voom_lavage_wo_unknown_w_viral_Screening <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_viral_subset_Screening",
    ext = "RDS"
  ),
  {
    libs_screening_lavage <- dat_voom_lavage_wo_unknown_w_viral$targets %>%
      subset(visit == "Screening") %>%
      pull(libid)
    dat_voom_lavage_wo_unknown_w_viral_Screening <- dat_voom_lavage_wo_unknown_w_viral[,
      libs_screening_lavage
    ]
    dat_voom_lavage_wo_unknown_w_viral_Screening
  },
  overwrite = FALSE
)
```

## Asthma vs Non-Asthma Comparison at Screening

### KIMMA Modeling for Asthma Status

```{r}
file.tmp <- file.path(
  dirresults,
  paste(projectNumber_lavage, "_Astham_vs_NonAsthma_Screening.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  model_Astham_vs_NonAsthma_Screening_kimma <- kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral_Screening,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ asthma_status",
    run_lme = FALSE,
    run_lm = TRUE,
    run_contrast = FALSE,
    use_weights = TRUE,
    metrics = TRUE,
    processors = 6
  )
  saveRDS(model_Astham_vs_NonAsthma_Screening_kimma, file = file.tmp)
} else {
  model_Astham_vs_NonAsthma_Screening_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_Astham_vs_NonAsthma_Screening_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.2, 0.25, 0.5)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

kimma::summarize_kmFit(
  model_Astham_vs_NonAsthma_Screening_kimma$lm,
  FCgroup = FALSE,
  p_cutoff = c(0.01, 0.05)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```

## Eosinophil Category Analysis at Screening

### KIMMA Modeling for Eosinophil Categories
```{r}
file.tmp <- file.path(
  dirresults,
  paste(projectNumber_lavage, "_EOS_Category_Screening.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  model_EOS_cat_Screening_kimma <- kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral_Screening,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ EOS_cat",
    run_lme = FALSE,
    run_lm = TRUE,
    run_contrast = FALSE,
    use_weights = TRUE,
    metrics = TRUE,
    processors = 6
  )
  saveRDS(model_EOS_cat_Screening_kimma, file = file.tmp)
} else {
  model_EOS_cat_Screening_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_EOS_cat_Screening_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

kimma::summarize_kmFit(
  model_EOS_cat_Screening_kimma$lm,
  FCgroup = FALSE,
  p_cutoff = c(0.01, 0.05)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

model_EOS_cat_Screening_kimma_sig_genes_EOS_300 <- model_EOS_cat_Screening_kimma$lm %>%
  subset(FDR <= 0.05 & variable == "EOS_catEOS_Greater_than_300") %>%
  left_join(
    dat_voom_lavage_wo_unknown_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  )
write_csv(
  model_EOS_cat_Screening_kimma_sig_genes_EOS_300,
  file.path(
    dirresults,
    paste(
      projectNumber_lavage,
      "_EOS_Category_Screening_sig_genes_EOS_300.csv",
      sep = ""
    )
  )
)

eos1_genes <- muppits_mods %>%
  subset(geneSet == "eos1") %>%
  pull(symbol) %>%
  unique()

intersect(
  model_EOS_cat_Screening_kimma_sig_genes_EOS_300$hgnc_symbol,
  eos1_genes
)
```




# T2 Endotype Analysis

## Data Preparation for T2 High vs Low Analysis

### Load T2 Classification Based on EOS1 Module Expression

```{r}
dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2 <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_viral_subset_Screening_w_eos1_based_T2",
    ext = "RDS"
  ),
  {
    T2_HIgh_Low_classification_lavage <- read_csv(
      HoardR::make_path(
        dir = dirDataclean,
        project = projectNumber_lavage,
        analysis = "Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_lavage",
        ext = "csv"
      )
    ) %>%
      janitor::clean_names()
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2 <- dat_voom_lavage_wo_unknown_w_viral_Screening
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets <- dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      left_join(
        T2_HIgh_Low_classification_lavage %>%
          select(subject_id, gene_exp, t2_mean, t2_1sd),
        by = c("subjectId" = "subject_id")
      ) %>%
      mutate(t2_1sd = factor(t2_1sd, levels = c("T2-low", "T2-high"))) %>%
      mutate(
        T2_w_Non_Asthma = ifelse(
          as.character(asthma_status) == "Non_Asthma",
          "Non_Asthma",
          as.character(t2_1sd)
        )
      ) %>%
      mutate(
        T2_w_Non_Asthma = factor(
          T2_w_Non_Asthma,
          levels = c("Non_Asthma", "T2-low", "T2-high")
        )
      )
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2
  },
  overwrite = FALSE
)
```

## T2 High vs Low Differential Expression Analysis

### KIMMA Modeling for T2 High vs Low (DONT run this because did not add Median CV coverage)
```{r}
model_eos1_based_T2_High_vs_Low_Screening_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_Screening",
    ext = "RDS"
  ),
  kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ t2_1sd",
    run_lme = FALSE,
    run_lm = TRUE,
    run_contrast = FALSE,
    use_weights = TRUE,
    metrics = TRUE,
    processors = 6
  ),
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_based_T2_High_vs_Low_Screening_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_based_T2_High_vs_Low_Screening_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_Screening_kimma_dat",
    ext = "csv"
  ),
  model_eos1_based_T2_High_vs_Low_Screening_kimma$lm %>%
    left_join(
      dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
        select(geneName, hgnc_symbol),
      by = c("gene" = "geneName")
    ) %>%
    mutate(
      comparison = ifelse(
        variable == "t2_1sdT2-high",
        "T2_high_vs_T2_low",
        variable
      )
    ) %>%
    select(
      gene,
      hgnc_symbol,
      variable,
      comparison,
      df,
      statistic,
      estimate,
      pval,
      FDR
    ),
  overwrite = FALSE
)
```

### KIMMA Modeling for T2 High vs Low with Median cv coverage
```{r}
model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_w_median_cv_Screening",
    ext = "RDS"
  ),
  kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ t2_1sd + median_cv_coverage",
    run_lme = FALSE,
    run_lm = TRUE,
    run_contrast = FALSE,
    use_weights = TRUE,
    metrics = TRUE,
    processors = 6
  ),
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma_dat",
    ext = "csv"
  ),
  model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma$lm %>%
    left_join(
      dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
        select(geneName, hgnc_symbol),
      by = c("gene" = "geneName")
    ) %>%
    mutate(
      comparison = ifelse(
        variable == "t2_1sdT2-high",
        "T2_high_vs_T2_low",
        variable
      )
    ) %>%
    select(
      gene,
      hgnc_symbol,
      variable,
      comparison,
      df,
      statistic,
      estimate,
      pval,
      FDR
    ),
  overwrite = FALSE
)
```

### T2 High vs Low Volcano Plot

```{r}
Volcano_plot_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage <- volcano_funcxn(
  model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma_dat %>%
    subset(comparison == "T2_high_vs_T2_low") %>%
    mutate(volcano_comp = "T2 Low<----->T2 High"),
  fdr = 0.05,
  n_top_gene = 15,
  add_genes_int = c("None")
)

Volcano_plot_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage$`T2 Low<----->T2 High` +
  theme_classic_custom

HoardR::save_figures(
  plot = Volcano_plot_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage$`T2 Low<----->T2 High` +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_volcano_lavage,
    project = projectNumber_lavage,
    analysis = "Volcano_plot_eos1_based_T2_High_vs_Low_w_median_cv_Screening",
    ext = "png"
  ),
  height = 6,
  width = 7,
  overwrite = FALSE
)
```

### T2 High vs Low Enrichment Analysis

```{r}
Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage <- list()
Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage[[
  "Up in T2 High"
]] <- model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma_dat %>%
  subset(comparison == "T2_high_vs_T2_low") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage[[
  "Up in T2 Low"
]] <- model_eos1_based_T2_High_vs_Low_w_median_cv_Screening_kimma_dat %>%
  subset(comparison == "T2_high_vs_T2_low") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list <- list()

    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list[[
      "GO_BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage,
      collection = "H"
    )

    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list
  },
  overwrite = FALSE
)


BIGpicture::plot_enrich(
  enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$H,
  fdr_cutoff = c(0.2)
)
BIGpicture::plot_enrich(
  enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$muppits,
  fdr_cutoff = c(0.2)
)
BIGpicture::plot_enrich(
  enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$GO_BP,
  fdr_cutoff = c(0.01)
)

comparisons_to_plot <- list(
  "Group 1:T2 High vs Group 2:T2 Low " = c(
    "Up in T2 High",
    "Up in T2 Low"
  )
)

plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$H,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage
)

HoardR::save_figures(
  plot = plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening",
    ext = "png"
  ),
  height = 4.5,
  width = 8,
  overwrite = FALSE
)

plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$GO_BP,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage
)

HoardR::save_figures(
  plot = plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening",
    ext = "png"
  ),
  height = 20,
  width = 15,
  overwrite = FALSE
)


plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_top_15_term <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$GO_BP %>%
    group_by(group) %>%
    slice_min(order_by = FDR, n = 15) %>%
    slice_max(order_by = `k/K`, n = 15),
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_top_15_term
)

HoardR::save_figures(
  plot = plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_top_15_term +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_top_15_term",
    ext = "png"
  ),
  height = 9,
  width = 10,
  overwrite = TRUE
)


plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage_GO_H_muppits_list$muppits,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.001
) +
  theme_classic_custom
print(
  plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage
)

HoardR::save_figures(
  plot = plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening",
    ext = "png"
  ),
  height = 4.5,
  width = 8,
  overwrite = TRUE
)

HoardR::save_figures(
  plot = plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage +
    theme_classic_custom +
    plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening_lavage +
    theme_classic_custom +
    plot_layout(nrow = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "muppits_H_combined_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_w_median_cv_Screening",
    ext = "png"
  ),
  height = 5,
  width = 11,
  overwrite = TRUE
)
```




### STRING Network Analysis for T2 High vs Low

```{r}  
#dir.create(file.path(dirPlots_lavage, "String_Network"), showWarnings = FALSE,recursive = TRUE)
for (comp_ in names(
  Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage
)) {
  pullString::get_png_network(
    ids = Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage[[comp_]],
    flat_node_design = 1, # disable 3D bubble design
    block_structure_pics_in_bubbles = 1, # removes 3d protein structure pictures inside the bubble
    png_filename = file.path(
      dirPlots_lavage,
      "String_Network",
      paste(
        projectNumber_lavage,
        "_String_Network_T2_High_vs_Low_Screening_",
        gsub(pattern = " ", replacement = "_", comp_),
        ".png",
        sep = ""
      )
    )
  )
}
```

## Non-Asthma vs T2 High vs Low Differential Expression Analysis
### KIMMA Modeling for Non-Asthma vs T2 High vs Low
```{r}
model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening",
    ext = "RDS"
  ),
  kimma::kmFit(
    dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ T2_w_Non_Asthma",
    run_lme = FALSE,
    run_lm = TRUE,
    run_contrast = TRUE,
    use_weights = TRUE,
    metrics = TRUE,
    processors = 6
  ),
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

kimma::summarize_kmFit(
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma$lm.contrast,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2, 0.25, 0.3)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

kimma::summarize_kmFit(
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma$lm.contrast,
  FCgroup = FALSE,
  p_cutoff = c(0.01, 0.05, 0.2, 0.25, 0.3)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat",
    ext = "csv"
  ),
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma$lm.contrast %>%
    left_join(
      dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
        select(geneName, hgnc_symbol),
      by = c("gene" = "geneName")
    ) %>%
    mutate(
      comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_")
    ) %>%
    select(
      gene,
      hgnc_symbol,
      variable,
      contrast_lvl,
      contrast_ref,
      comparison,
      df,
      statistic,
      estimate,
      pval,
      FDR
    ),
  overwrite = FALSE
)
```



### Non-Asthma vs T2 High vs Low Volcano Plot

```{r}
Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage <- volcano_funcxn(
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
    mutate(
      volcano_comp = ifelse(
        comparison == "(T2-high)_vs_(T2-low)",
        "T2 Low<----->T2 High",
        ifelse(
          comparison == "(T2-low)_vs_Non_Asthma",
          "Non-Asthma<----->T2 Low",
          "Non-Asthma<----->T2 High"
        )
      )
    ),
  fdr = 0.2,
  n_top_gene = 10,
  add_genes_int = c("None")
)

Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage$`T2 Low<----->T2 High` +
  theme_classic_custom +
  Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage$`Non-Asthma<----->T2 High` +
  theme_classic_custom

HoardR::save_figures(
  plot = Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage$`T2 Low<----->T2 High` +
    theme_classic_custom +
    Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage$`Non-Asthma<----->T2 High` +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_volcano_lavage,
    project = projectNumber_lavage,
    analysis = "Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 6,
  width = 12,
  overwrite = FALSE
)
```

### Non-Asthma vs T2 High vs Low  Enrichment Analysis

```{r}
Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage <- list()

Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
  "Up in T2 High vs T2 Low"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high)_vs_(T2-low)") %>%
  subset(estimate > 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
  "Up in T2 Low vs T2 High"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high)_vs_(T2-low)") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
  "Up in T2 High vs Non-Asthma"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high)_vs_Non_Asthma") %>%
  subset(estimate > 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
  "Up in Non-Asthma vs T2 High"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high)_vs_Non_Asthma") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()


Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
  "Up in T2 Low vs Non-Asthma (pval)"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-low)_vs_Non_Asthma") %>%
  subset(estimate > 0 & pval <= 0.01) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
  "Up in Non-Asthma vs T2 Low (pval)"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-low)_vs_Non_Asthma") %>%
  subset(estimate < 0 & pval <= 0.01) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list <- list()
    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list[[
      "GO_BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage,
      collection = "H"
    )

    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list
  },
  overwrite = TRUE
)


comparisons_to_plot <- list(
  "Group 1:T2 High vs Group 2:T2 Low " = c(
    "Up in T2 High vs T2 Low",
    "Up in T2 Low vs T2 High"
  ),
  "Group 1:T2 High vs Group 2:Non-Asthma " = c(
    "Up in T2 High vs Non-Asthma",
    "Up in Non-Asthma vs T2 High"
  ),
  "Group 1:T2 Low vs Group 2:Non-Asthma \n (pvalue) " = c(
    "Up in T2 Low vs Non-Asthma (pval)",
    "Up in Non-Asthma vs T2 Low (pval)"
  )
)

plot_H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list$H,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  plot_H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage
)

HoardR::save_figures(
  plot = plot_H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_with_pval_cutoff_Non_Asthma_vs_T2_Low_w_pvalue_Non_Asthma_vs_T2_Low",
    ext = "png"
  ),
  height = 5,
  width = 12,
  overwrite = FALSE
)

plot_muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list$muppits,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  plot_muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage
)

HoardR::save_figures(
  plot = plot_muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_w_pvalue_Non_Asthma_vs_T2_Low",
    ext = "png"
  ),
  height = 5,
  width = 12,
  overwrite = FALSE
)


plot_GO_BP_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage_GO_H_muppits_list$GO_BP,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  plot_GO_BP_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage
)

HoardR::save_figures(
  plot = plot_GO_BP_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_lavage,
    project = projectNumber_lavage,
    analysis = "GO_BP_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 30,
  width = 20,
  overwrite = FALSE
)
```

### STRING Network Analysis for Non-Asthma vs T2 High vs Low

```{r}  
#dir.create(file.path(dirPlots_lavage, "String_Network"), showWarnings = FALSE,recursive = TRUE)
for (comp_ in names(
  Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage
)) {
  pullString::get_png_network(
    ids = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_lavage[[
      comp_
    ]],
    flat_node_design = 1, # disable 3D bubble design
    block_structure_pics_in_bubbles = 1, # removes 3d protein structure pictures inside the bubble
    png_filename = file.path(
      dirPlots_lavage,
      "String_Network",
      paste(
        projectNumber_lavage,
        "_String_Network_T2_High_vs_Low_Screening_",
        gsub(pattern = " ", replacement = "_", comp_),
        ".png",
        sep = ""
      )
    )
  )
}
```






## T2 High vs Low by Asthma Status Interaction Analysis

### KIMMA Modeling with T2-Asthma Interaction

```{r}
model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_by_asthma_status_Screening",
    ext = "RDS"
  ),
  {
    model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_kimma <- kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ t2_1sd*asthma_status",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = TRUE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm",
    ext = "csv"
  ),
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma$lm %>%
    left_join(
      dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
        select(geneName, hgnc_symbol),
      by = c("gene" = "geneName")
    ) %>%
    mutate(
      comparison = ifelse(
        variable == "t2_1sdT2-high",
        "T2_high_vs_T2_low",
        ifelse(
          variable == "t2_1sdT2-high:asthma_statusNon_Asthma",
          "Interaction_T2_high_vs_T2_low_by_asthma_status",
          ifelse(
            variable == "asthma_statusNon_Asthma",
            "NonAsthma_vs_Asthma",
            variable
          )
        )
      )
    ) %>%
    select(
      model,
      gene,
      hgnc_symbol,
      variable,
      comparison,
      df,
      statistic,
      estimate,
      pval,
      FDR
    ),
  overwrite = FALSE
)


kimma::summarize_kmFit(
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma$lm.contrast %>%
    mutate(
      comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_")
    ) %>%
    subset(
      comparison %in%
        c(
          "Non_Asthma_vs_Asthma",
          "(T2-high)_vs_(T2-low)",
          "(T2-high Asthma)_vs_(T2-low Asthma)",
          "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)",
          "(T2-high Non_Asthma)_vs_(T2-high Asthma)",
          "(T2-low Non_Asthma)_vs_(T2-low Asthma)"
        )
    ),
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2, 0.25)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast",
    ext = "csv"
  ),
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma$lm.contrast %>%
    left_join(
      dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
        select(geneName, hgnc_symbol),
      by = c("gene" = "geneName")
    ) %>%
    mutate(
      comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_")
    ) %>%
    subset(
      comparison %in%
        c(
          "(T2-high Asthma)_vs_(T2-low Asthma)",
          "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)",
          "(T2-high Non_Asthma)_vs_(T2-high Asthma)",
          "(T2-low Non_Asthma)_vs_(T2-low Asthma)"
        )
    ) %>%
    select(
      model,
      gene,
      hgnc_symbol,
      variable,
      comparison,
      df,
      statistic,
      estimate,
      pval,
      FDR
    ),
  overwrite = FALSE
)
```

#### Volcano Plots for Interaction Analysis
```{r}
Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_lavage <- volcano_funcxn(
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat %>%
    subset(comparison %in% c("T2_high_vs_T2_low", "NonAsthma_vs_Asthma")) %>%
    mutate(
      volcano_comp = ifelse(
        comparison == "T2_high_vs_T2_low",
        "T2 Low<----->T2 High",
        "Non-Asthma<----->Asthma"
      )
    ),
  fdr = 0.25,
  n_top_gene = 10,
  add_genes_int = c("None")
)

volcano_plot_T2_High_vs_Low_T2_Asthma_status_model_lavage <- Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_lavage$`T2 Low<----->T2 High` +
  theme_classic_custom
print(volcano_plot_T2_High_vs_Low_T2_Asthma_status_model_lavage)

HoardR::save_figures(
  plot = volcano_plot_T2_High_vs_Low_T2_Asthma_status_model_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_T2_High_vs_Low_T2_Asthma_status_model",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

volcano_plot_NonAsthma_vs_Asthma_T2_High_vs_Low_model_lavage <- Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_lavage$`Non-Asthma<----->Asthma` +
  theme_classic_custom

print(volcano_plot_NonAsthma_vs_Asthma_T2_High_vs_Low_model_lavage)

HoardR::save_figures(
  plot = volcano_plot_NonAsthma_vs_Asthma_T2_High_vs_Low_model_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_NonAsthma_vs_Asthma_T2_High_vs_Low_model",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```

##### Make volcano plot for lm contrasts
```{r}
Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_contrast_lavage <- volcano_funcxn(
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
    mutate(
      volcano_comp = ifelse(
        comparison == "(T2-high Asthma)_vs_(T2-low Asthma)",
        "T2 Low<----->T2 High in Asthma",
        ifelse(
          comparison == "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)",
          "T2 Low<----->T2 High in Non-Asthma",
          ifelse(
            comparison == "(T2-high Non_Asthma)_vs_(T2-high Asthma)",
            "Asthma<----->Non-Asthma in T2 High",
            "Asthma<----->Non-Asthma in T2 Low"
          )
        )
      )
    ),
  fdr = 0.2,
  n_top_gene = 10,
  add_genes_int = c("None")
)

volcano_plot_T2_High_vs_Low_in_Asthma_T2_Asthma_status_model_lavage <- Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_contrast_lavage$`T2 Low<----->T2 High in Asthma` +
  theme_classic_custom

volcano_plot_T2_High_vs_Low_in_Non_Asthma_T2_Asthma_status_model_lavage <- Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_contrast_lavage$`T2 Low<----->T2 High in Non-Asthma` +
  theme_classic_custom

HoardR::save_figures(
  plot = volcano_plot_T2_High_vs_Low_in_Asthma_T2_Asthma_status_model_lavage +
    volcano_plot_T2_High_vs_Low_in_Non_Asthma_T2_Asthma_status_model_lavage +
    plot_layout(guides = "collect", ncol = 2),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_T2_High_vs_Low_in_Asthma_and_Non_Asthma_contrast_T2_Asthma_status_model",
    ext = "png"
  ),
  height = 5,
  width = 10,
  overwrite = FALSE
)

volcano_plot_Non_Asthma_vs_Asthma_in_T2_High_T2_Asthma_status_model_lavage <- Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_contrast_lavage$`Asthma<----->Non-Asthma in T2 High` +
  theme_classic_custom

volcano_plot_Non_Asthma_vs_Asthma_in_T2_Low_T2_Asthma_status_model_lavage <- Volcano_plot_eos1_based_T2_High_vs_Low_Screening_lm_contrast_lavage$`Asthma<----->Non-Asthma in T2 Low` +
  theme_classic_custom

HoardR::save_figures(
  plot = volcano_plot_Non_Asthma_vs_Asthma_in_T2_High_T2_Asthma_status_model_lavage +
    volcano_plot_Non_Asthma_vs_Asthma_in_T2_Low_T2_Asthma_status_model_lavage +
    plot_layout(guides = "collect", ncol = 2),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_Non_Asthma_vs_Asthma_in_T2_High_and_T2_Low_contrast_T2_Asthma_status_model",
    ext = "png"
  ),
  height = 5,
  width = 10,
  overwrite = FALSE
)
```

#### Make a scatterplot of the estimates of T2 high vs low in asthma vs non-asthma and Non-Asthma vs Asthma in T2 high and T2 low

##### Scatterplot of T2 High vs Low in Asthma vs Non-Asthma
```{r}
scatterplot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_plot <- make_contrast_scatterplot(
  dat = model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat,
  comp_x = "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)",
  comp_y = "(T2-high Asthma)_vs_(T2-low Asthma)",
  xlab = "Estimate of T2 High vs Low in Non-Asthma",
  ylab = "Estimate of T2 High vs Low in Asthma",
  group_x = "Non-Asthma",
  group_y = "Asthma",
  fdr_cutoff = 0.05,
  top_n = 10
)

HoardR::save_figures(
  plot = scatterplot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_plot$plot +
    scale_color_manual(
      values = c(
        "Significant in both" = "purple",
        "Significant in Asthma only" = "darkred",
        "Significant in Non-Asthma only" = "orange",
        "Not significant" = "grey"
      )
    ),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_contrast_scatter_lavage,
    analysis = "Scatter_plot_T2_High_vs_Low_in_Asthma_vs_Non_Asthma_eos1_based_T2_Asthma_model",
    ext = "png"
  ),
  height = 7,
  width = 10,
  overwrite = FALSE
)
```

##### Scatterplot of Non-Asthma vs Asthma in T2 High vs T2 Low
```{r}
scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_T2_High_and_T2_Low_status_Screening_lavage_kimma_lm_contrast_plot <- make_contrast_scatterplot(
  dat = model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat,
  comp_x = "(T2-high Non_Asthma)_vs_(T2-high Asthma)",
  comp_y = "(T2-low Non_Asthma)_vs_(T2-low Asthma)",
  xlab = "Estimate of Non_Asthma vs Asthma in T2 High",
  ylab = "Estimate of Non_Asthma vs Asthma in T2 Low",
  group_x = "T2-High",
  group_y = "T2-Low",
  fdr_cutoff = 0.25,
  top_n = 10
)

HoardR::save_figures(
  plot = scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_T2_High_and_T2_Low_status_Screening_lavage_kimma_lm_contrast_plot$plot +
    scale_color_manual(
      values = c(
        "Significant in both" = "purple",
        "Significant in T2-High only" = "darkred",
        "Significant in T2-Low only" = "orange",
        "Not significant" = "grey"
      )
    ),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_contrast_scatter_lavage,
    analysis = "Scatter_plot_Asthma_vs_Non_Asthma_in_T2_High_and_T2_Low_eos1_based_T2_Asthma_model",
    ext = "png"
  ),
  height = 7,
  width = 10,
  overwrite = FALSE
)
```

#### Run Enrichment
```{r}
Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage <- list()
Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage[[
  "Up in T2 High"
]] <- model_eos1_based_T2_High_vs_Low_Screening_kimma_dat %>%
  subset(comparison == "T2_high_vs_T2_low") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage[[
  "Up in T2 Low"
]] <- model_eos1_based_T2_High_vs_Low_Screening_kimma_dat %>%
  subset(comparison == "T2_high_vs_T2_low") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_GO_Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage,
  collection = "C5",
  subcollection = "GO:BP"
)

enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage,
  collection = "H"
)

BIGpicture::plot_enrich(
  enrichment_GO_Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage,
  fdr_cutoff = c(0.01)
)
BIGpicture::plot_enrich(
  enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_Screening_lavage
)
```


#### Run Enrichment with T2 High vs Low interaction with Asthma and Non-Asthma 
```{r}
Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- list()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in T2 High"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat %>%
  subset(comparison == "T2_high_vs_T2_low") %>%
  subset(estimate > 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in T2 Low"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat %>%
  subset(comparison == "T2_high_vs_T2_low") %>%
  subset(estimate < 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in Non-Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat %>%
  subset(comparison == "NonAsthma_vs_Asthma") %>%
  subset(estimate > 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat %>%
  subset(comparison == "NonAsthma_vs_Asthma") %>%
  subset(estimate < 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Interaction T2 High vs Low by Asthma Status"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_dat %>%
  subset(comparison == "Interaction_T2_high_vs_T2_low_by_asthma_status") %>%
  subset(FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()


Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in T2 High in Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high Asthma)_vs_(T2-low Asthma)") %>%
  subset(estimate > 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in T2 Low in Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high Asthma)_vs_(T2-low Asthma)") %>%
  subset(estimate < 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in T2 High in Non-Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)") %>%
  subset(estimate > 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in T2 Low in Non-Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)") %>%
  subset(estimate < 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in Non-Asthma in T2 Low"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-low Non_Asthma)_vs_(T2-low Asthma)") %>%
  subset(estimate > 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in Asthma in T2 Low"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-low Non_Asthma)_vs_(T2-low Asthma)") %>%
  subset(estimate < 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()


Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in Non-Asthma in T2 High"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high Non_Asthma)_vs_(T2-high Asthma)") %>%
  subset(estimate > 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage[[
  "Up in Asthma in T2 High"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high Non_Asthma)_vs_(T2-high Asthma)") %>%
  subset(estimate < 0 & FDR <= 0.25) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "Enrichment_GO_Hallmark_muppits_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list <- list()

    enrichment_GO_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list[[
      "GO:BP"
    ]] <- enrichment_GO_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage

    enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage,
      collection = "H"
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list[[
      "H"
    ]] <- enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage

    enrichment_muppits_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- enrichment_muppits_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage

    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list
  },
  overwrite = FALSE
)
comparisons_to_plot <- list(
  "Group 1:T2 High vs Group 2:T2 Low \n(Non-Asthma)" = c(
    "Up in T2 High in Non-Asthma",
    "Up in T2 Low in Non-Asthma"
  ),
  "Group 1:T2 High vs Group 2:T2 Low \n(Asthma)" = c(
    "Up in T2 High in Asthma",
    "Up in T2 Low in Asthma"
  ),

  "Group 1:Asthma vs Group 2:Non-Asthma \n(T2 High)" = c(
    "Up in Asthma in T2 High",
    "Up in Non-Asthma in T2 High"
  ),
  "Group 1:Asthma vs Group 2:Non-Asthma \n(T2 Low)" = c(
    "Up in Asthma in T2 Low",
    "Up in Non-Asthma in T2 Low"
  )
)
plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list$H,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom

HoardR::save_figures(
  plot = plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_interaction_model_Screening",
    ext = "png"
  ),
  height = 7,
  width = 14,
  overwrite = FALSE
)

plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list$`GO:BP`,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.001
) +
  theme_classic_custom

HoardR::save_figures(
  plot = plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_interaction_model_Screening",
    ext = "png"
  ),
  height = 35,
  width = 30,
  overwrite = FALSE
)

plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage_GO_H_muppits_list$muppits,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom

HoardR::save_figures(
  plot = plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_muppits_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_interaction_model_Screening",
    ext = "png"
  ),
  height = 7,
  width = 13,
  overwrite = FALSE
)
```


# Plot the muppits modules by T2 High vs Low in Asthma vs Non-Asthma in screening only
```{r}
muppits_mods_expression_w_meta_wo_unknown_Screening_lavage <- muppits_mods_expression_w_meta_all_lavage %>%
  subset(
    libid %in%
      dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets$libid
  ) %>%
  select(libid, geneSet, geneExp) %>%
  left_join(
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets
  )

colors_T2_w_Non_Asthma <- make_color_palette(
  variables = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
    dplyr::arrange(desc(T2_w_Non_Asthma)) %>%
    dplyr::pull(T2_w_Non_Asthma),
  registry_key = "T2_w_Non_Asthma"
)

muppits_mods_expression_w_meta_wo_unknown_Screening_lavage %>%
  ggplot(aes(x = T2_w_Non_Asthma, y = geneExp, color = T2_w_Non_Asthma)) +
  ggbeeswarm::geom_beeswarm(width = 0.5) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) +
  facet_wrap(~geneSet) +
  scale_color_manual(values = colors_T2_w_Non_Asthma) +
  theme_classic_custom

muppits_mods_of_int <- c(
  "eos1",
  "eos5",
  "neut1",
  "unma23",
  "epi1",
  "unma11",
  "unma13",
  "unma24"
)


boxplot_Non_Asthma_T2_High_Low_muppits_mod_of_interest_Screening_lavage <- muppits_mods_expression_w_meta_wo_unknown_Screening_lavage %>%
  subset(geneSet %in% muppits_mods_of_int) %>%
  ggplot(aes(x = T2_w_Non_Asthma, y = geneExp, color = T2_w_Non_Asthma)) +
  ggbeeswarm::geom_beeswarm(dodge.width = 0.5) +
  geom_boxplot(alpha = 0.5, outlier.shape = NA) +
  facet_wrap(~geneSet, nrow = 2) +
  scale_color_manual(values = colors_T2_w_Non_Asthma) +
  theme_classic_custom

HoardR::save_figures(
  plot = boxplot_Non_Asthma_T2_High_Low_muppits_mod_of_interest_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_boxplot_lavage,
    analysis = "Boxplot_Non_Asthma_T2_High_Low_muppits_modules_of_interest_Screening_lavage",
    ext = "png"
  ),
  height = 7,
  width = 12,
  overwrite = TRUE
)
```

# EOS1 Continuous Variable Analyses

## EOS1 Module Expression Analysis

### KIMMA Modeling with Continuous EOS1 Expression

```{r}
model_eos1_exp_Screening_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_Screening_kimma",
    ext = "RDS"
  ),
  {
    model_eos1_exp_Screening_lavage_kimma <- kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ gene_exp",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = FALSE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_exp_Screening_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  mutate(
    variable = ifelse(
      variable == "gene_exp",
      "EOS1_exp",
      "total (nonredundant)"
    )
  ) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_exp_Screening_lavage_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_Screening_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_Screening_lavage_kimma$lm %>%
      left_join(
        dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
          select(geneName, hgnc_symbol),
        by = c("gene" = "geneName")
      ) %>%
      mutate(
        comparison = ifelse(variable == "gene_exp", "EOS1_exp", variable)
      ) %>%
      select(
        model,
        gene,
        hgnc_symbol,
        variable,
        comparison,
        df,
        statistic,
        estimate,
        pval,
        FDR
      )
  },
  overwrite = FALSE
)


# with median cv
model_eos1_exp_median_cv_Screening_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_median_cv_Screening_kimma",
    ext = "RDS"
  ),
  {
    kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ gene_exp + median_cv_coverage",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = FALSE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_exp_median_cv_Screening_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  mutate(
    variable = ifelse(
      variable == "gene_exp",
      "EOS1_exp",
      "total (nonredundant)"
    )
  ) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)

model_wo_median_cv_genes <- model_eos1_exp_Screening_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(FDR <= 0.05) %>%
  pull(gene)

model_w_median_cv_genes <- model_eos1_exp_median_cv_Screening_lavage_kimma$lm %>%
  subset(variable == "gene_exp") %>%
  subset(FDR <= 0.05) %>%
  pull(gene)

# make a venn diagram of the significant genes in the model with and without median cv
venn_data <- list(
  "Without Median CV" = model_wo_median_cv_genes,
  "With Median CV" = model_w_median_cv_genes
)
venn_plot <- ggvenn::ggvenn(
  venn_data,
  fill_color = c("skyblue", "salmon"),
  stroke_size = 0.5,
  set_name_size = 4
) +
  theme_classic_custom
```

### EOS1 Expression Volcano Plot

```{r}
Volcano_plot_eos1_exp_Screening_lavage <- volcano_funcxn(
  model_eos1_exp_Screening_lavage_kimma_dat %>%
    subset(comparison == "EOS1_exp") %>%
    mutate(volcano_comp = "EOS1_exp"),
  fdr = 0.05,
  n_top_gene = 15,
  add_genes_int = c("None")
)

Volcano_plot_eos1_exp_Screening_lavage$EOS1_exp + theme_classic_custom


HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_Screening_lavage$EOS1_exp +
    theme_classic_custom +
    labs(title = "Negative Correlation <-----> Positive Correlation to EOS1"),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_EOS1_exp_Screening_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = TRUE
)
```

### EOS1 Expression Enrichment Analysis

```{r}
Significant_genes_eos1_exp_Screening_lavage <- list()
Significant_genes_eos1_exp_Screening_lavage[[
  "Positive Correlation to EOS1"
]] <- model_eos1_exp_Screening_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_Screening_lavage[[
  "Negative Correlation to EOS1"
]] <- model_eos1_exp_Screening_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "Enrichment_GO_Hallmark_muppits_Significant_genes_eos1_exp_Screening",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list <- list()
    enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list[[
      "GO:BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_lavage,
      collection = "H",
    )

    enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list
  },
  overwrite = FALSE
)

enrichment_comparison_plot <- list(
  "Group 1: Positive Correlation vs Group 2: Negative Correlation\nto EOS1" = c(
    "Positive Correlation to EOS1",
    "Negative Correlation to EOS1"
  )
)
plot_enrichment_H_Significant_genes_eos1_exp_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list$H,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.05
) +
  theme_classic_custom
print(plot_enrichment_H_Significant_genes_eos1_exp_Screening_lavage)

HoardR::save_figures(
  plot = plot_enrichment_H_Significant_genes_eos1_exp_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_Significant_genes_eos1_exp_Screening",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list$`GO:BP`,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.001
) +
  theme_classic_custom

print(plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_lavage)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_exp_Screening",
    ext = "png"
  ),
  height = 35,
  width = 25,
  overwrite = FALSE
)


plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_lavage_top_15_term <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list$`GO:BP` %>%
    group_by(group) %>%
    slice_min(order_by = FDR, n = 15) %>%
    slice_max(order_by = `k/K`, n = 15),
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom

print(
  plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_lavage_top_15_term
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_lavage_top_15_term,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_exp_Screening_top_15_term",
    ext = "png"
  ),
  height = 9,
  width = 10,
  overwrite = TRUE
)

plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_lavage_GO_H_muppits_list$muppits,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.05
) +
  theme_classic_custom

print(plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_lavage)

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_Screening",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_lavage +
  theme_classic_custom +
  plot_enrichment_H_Significant_genes_eos1_exp_Screening_lavage +
  theme_classic_custom +
  plot_layout(ncol = 2, guides = "collect")

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_lavage +
    theme_classic_custom +
    plot_enrichment_H_Significant_genes_eos1_exp_Screening_lavage +
    theme_classic_custom +
    plot_layout(ncol = 2, guides = "collect"),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_and_muppits_combined_Significant_genes_eos1_exp_Screening",
    ext = "png"
  ),
  height = 5,
  width = 14,
  overwrite = FALSE
)
```


### KIMMA Modeling with Continuous EOS1 Expression Interaction with Asthma status

```{r}
model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_interation_Asthma_status_Screening_kimma",
    ext = "rds"
  ),
  {
    kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ gene_exp * asthma_status",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = FALSE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  mutate(
    variable = ifelse(
      variable == "gene_exp",
      "EOS1_exp",
      as.character(variable)
    )
  ) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_interation_Asthma_status_Screening_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma$lm %>%
      left_join(
        dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
          select(geneName, hgnc_symbol),
        by = c("gene" = "geneName")
      ) %>%
      mutate(
        comparison = ifelse(
          variable == "gene_exp",
          "EOS1_exp",
          as.character(variable)
        )
      ) %>%
      select(
        model,
        gene,
        hgnc_symbol,
        variable,
        comparison,
        df,
        statistic,
        estimate,
        pval,
        FDR
      )
  },
  overwrite = FALSE
)
```

### EOS1 Expression in the Asthma Interaction Model Volcano Plot

```{r}
Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage <- volcano_funcxn(
  model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat %>%
    subset(comparison != "(Intercept)") %>%
    mutate(
      volcano_comp = ifelse(
        comparison == "asthma_statusNon_Asthma",
        "Non-Asthma vs Asthma",
        as.character(comparison)
      )
    ),
  fdr = 0.2,
  n_top_gene = 10,
  add_genes_int = c("None")
)

Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage$EOS1_exp +
  theme_classic_custom
HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage$EOS1_exp +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_EOS1_exp_interation_Asthma_status_Screening_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage$`Non-Asthma vs Asthma` +
  theme_classic_custom
HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage$`Non-Asthma vs Asthma` +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_Non_Asthma_vs_Asthma_interation_Asthma_status_Screening_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage$`gene_exp:asthma_statusNon_Asthma` +
  theme_classic_custom
HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_lavage$`gene_exp:asthma_statusNon_Asthma` +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_EOS1_exp_Asthma_Interaction_term_interation_Asthma_status_Screening_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```

### EOS1 Expression Enrichment Analysis

```{r}
Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage <- list()
Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage[[
  "Positive Correlation to EOS1"
]] <- model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate > 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage[[
  "Negative Correlation to EOS1"
]] <- model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage[[
  "Up in Non-Asthma vs Asthma"
]] <- model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat %>%
  subset(comparison == "asthma_statusNon_Asthma") %>%
  subset(estimate > 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage[[
  "Up in Asthma vs Non-Asthma"
]] <- model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat %>%
  subset(comparison == "asthma_statusNon_Asthma") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()


Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage[[
  "EOS1 Exp Interaction:Asthma"
]] <- model_eos1_exp_interation_Asthma_status_Screening_lavage_kimma_dat %>%
  subset(comparison == "gene_exp:asthma_statusNon_Asthma") %>%
  subset(FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "Enrichment_GO_Hallmark_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list <- list()
    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list[[
      "GO:BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage,
      collection = "H",
    )

    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list
  },
  overwrite = FALSE
)

enrichment_comparison_plot <- list(
  "EOS1 Exp" = c(
    "Positive Correlation to EOS1",
    "Negative Correlation to EOS1"
  )
)
plot_enrichment_H_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list$H,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  plot_enrichment_H_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage
)

HoardR::save_figures(
  plot = plot_enrichment_H_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_Significant_genes_eos1_exp_interation_Asthma_status_Screening",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

plot_enrichment_GO_BP_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list$`GO:BP`,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.001
) +
  theme_classic_custom

print(
  plot_enrichment_GO_BP_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_exp_interation_Asthma_status_Screening",
    ext = "png"
  ),
  height = 35,
  width = 25,
  overwrite = FALSE
)

plot_enrichment_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage_GO_H_muppits_list$muppits,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom

print(
  plot_enrichment_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage
)

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```

#### Compare the signal with and without Asthma interaction term
```{r}
upset(
  fromList(list(
    "Positive Correlation to EOS1 (Interaction)" = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1 (Interaction)" = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage$`Negative Correlation to EOS1`,
    "Positive Correlation to EOS1" = Significant_genes_eos1_exp_Screening_lavage$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1" = Significant_genes_eos1_exp_Screening_lavage$`Negative Correlation to EOS1`
  )),
  text.scale = 1.5,
  order.by = "freq",
  mainbar.y.label = "Significant genes\nFDR < 0.2",
  keep.order = TRUE,
  sets = names(list(
    "Positive Correlation to EOS1 (Interaction)" = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1 (Interaction)" = Significant_genes_eos1_exp_interation_Asthma_status_Screening_lavage$`Negative Correlation to EOS1`,
    "Positive Correlation to EOS1" = Significant_genes_eos1_exp_Screening_lavage$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1" = Significant_genes_eos1_exp_Screening_lavage$`Negative Correlation to EOS1`
  ))
)
```


### KIMMA Modeling with Continuous EOS1 Expression split by Asthma Status

#### Asthma only Model

```{r}
dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Asthma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_viral_subset_Screening_w_eos1_based_T2_Asthma",
    ext = "rds"
  ),
  {
    libs_asthma <- dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      subset(asthma_status == "Asthma") %>%
      pull(libid)
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Asthma <- dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2[,
      libs_asthma
    ]
  },
  overwrite = FALSE
)

model_eos1_exp_Screening_Asthma_only_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_Screening_Asthma_only_kimma",
    ext = "rds"
  ),
  {
    model_eos1_exp_Screening_Asthma_only_lavage_kimma <- kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Asthma,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ gene_exp",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = FALSE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_exp_Screening_Asthma_only_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.25)
) %>%
  mutate(
    variable = ifelse(
      variable == "gene_exp",
      "EOS1_exp",
      "total (nonredundant)"
    )
  ) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_exp_Screening_Asthma_only_lavage_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_Screening_Asthma_only_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_Screening_Asthma_only_lavage_kimma$lm %>%
      left_join(
        dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Asthma$genes %>%
          select(geneName, hgnc_symbol),
        by = c("gene" = "geneName")
      ) %>%
      mutate(comparison = ifelse(variable == "gene_exp", "EOS1_exp", variable))
  },
  overwrite = FALSE
)
```

##### EOS1 Expression Volcano Plot

```{r}
Volcano_plot_eos1_exp_Screening_Asthma_only_lavage <- volcano_funcxn(
  model_eos1_exp_Screening_Asthma_only_lavage_kimma_dat %>%
    subset(comparison == "EOS1_exp") %>%
    mutate(volcano_comp = "EOS1_exp"),
  fdr = 0.05,
  n_top_gene = 10,
  add_genes_int = c("None")
)

Volcano_plot_eos1_exp_Screening_Asthma_only_lavage$EOS1_exp +
  theme_classic_custom

HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_Screening_Asthma_only_lavage$EOS1_exp +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_eos1_exp_Screening_Asthma_only_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```

### EOS1 Expression Enrichment Analysis

```{r}
Significant_genes_eos1_exp_Screening_Asthma_only_lavage <- list()

Significant_genes_eos1_exp_Screening_Asthma_only_lavage[[
  "Positive Correlation to EOS1"
]] <- model_eos1_exp_Screening_Asthma_only_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_Screening_Asthma_only_lavage[[
  "Negative Correlation to EOS1"
]] <- model_eos1_exp_Screening_Asthma_only_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list <- list()

    enrichment_GO_Significant_genes_eos1_exp_Screening_Asthma_only_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Asthma_only_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list[[
      "GO:BP"
    ]] <- enrichment_GO_Significant_genes_eos1_exp_Screening_Asthma_only_lavage
    enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Asthma_only_lavage,
      collection = "H"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list[[
      "H"
    ]] <- enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_lavage

    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Asthma_only_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list
  },
  overwrite = FALSE
)


enrichment_comparison_plot <- list(
  "EOS1 Exp" = c(
    "Positive Correlation to EOS1",
    "Negative Correlation to EOS1"
  )
)

plot_enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list$H,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(plot_enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_lavage)

HoardR::save_figures(
  plot = plot_enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list$`GO:BP`,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_lavage
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_lavage",
    ext = "png"
  ),
  height = 28,
  width = 20,
  overwrite = TRUE
)


plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_lavage_GO_H_muppits_list$muppits,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_lavage
)

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_lavage",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```


#### Non_Asthma only Model

```{r}
dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Non_Asthma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_viral_subset_Screening_w_eos1_based_T2_Non_Asthma",
    ext = "rds"
  ),
  {
    libs_Non_Asthma <- dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      subset(asthma_status == "Non_Asthma") %>%
      pull(libid)
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Non_Asthma <- dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2[,
      libs_Non_Asthma
    ]
  },
  overwrite = TRUE
)


model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_Screening_Non_Asthma_only_kimma",
    ext = "rds"
  ),
  {
    kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Non_Asthma,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ gene_exp",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = FALSE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.1, 0.2, 0.25)
) %>%
  mutate(
    variable = ifelse(
      variable == "gene_exp",
      "EOS1_exp",
      "total (nonredundant)"
    )
  ) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "eos1_exp_Screening_Non_Asthma_only_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma$lm %>%
      left_join(
        dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2_Non_Asthma$genes %>%
          select(geneName, hgnc_symbol),
        by = c("gene" = "geneName")
      ) %>%
      mutate(
        comparison = ifelse(variable == "gene_exp", "EOS1_exp", variable)
      ) %>%
      select(
        model,
        gene,
        hgnc_symbol,
        variable,
        comparison,
        df,
        statistic,
        estimate,
        pval,
        FDR
      )
  },
  overwrite = FALSE
)
```

##### EOS1 Expression Volcano Plot

```{r}
Volcano_plot_eos1_exp_Screening_Non_Asthma_only_lavage <- volcano_funcxn(
  model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma_dat %>%
    subset(comparison == "EOS1_exp") %>%
    mutate(volcano_comp = "EOS1_exp"),
  fdr = 0.1,
  n_top_gene = 10,
  add_genes_int = c("None")
)

Volcano_plot_eos1_exp_Screening_Non_Asthma_only_lavage$EOS1_exp +
  theme_classic_custom

HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_Screening_Non_Asthma_only_lavage$EOS1_exp +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_eos1_exp_Screening_Non_Asthma_only_lavage",
    ext = "png"
  ),
  height = 6,
  width = 8,
  overwrite = FALSE
)
```

### EOS1 Expression Enrichment Analysis

```{r}
Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage <- list()

Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage[[
  "Positive Correlation to EOS1"
]] <- model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate > 0 & FDR <= 0.1) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage[[
  "Negative Correlation to EOS1"
]] <- model_eos1_exp_Screening_Non_Asthma_only_lavage_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate < 0 & FDR <= 0.1) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list <- list()
    enrichment_GO_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list[[
      "GO:BP"
    ]] <- enrichment_GO_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage
    enrichment_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage,
      collection = "H"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list[[
      "H"
    ]] <- enrichment_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage

    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage,
      db = muppits_mods_for_enrichment
    )

    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list
  },
  overwrite = FALSE
)


enrichment_plot_GO_BP_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list$`GO:BP`,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.01
) +
  theme_classic_custom
print(
  enrichment_plot_GO_BP_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage
)

HoardR::save_figures(
  plot = enrichment_plot_GO_BP_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage",
    ext = "png"
  ),
  height = 20,
  width = 30,
  overwrite = FALSE
)

enrichment_plot_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list$H,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  enrichment_plot_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage
)

HoardR::save_figures(
  plot = enrichment_plot_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage",
    ext = "png"
  ),
  height = 6,
  width = 10,
  overwrite = FALSE
)


enrichment_plot_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage_GO_H_muppits_list$muppits,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  enrichment_plot_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage
)

HoardR::save_figures(
  plot = enrichment_plot_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_lavage",
    ext = "png"
  ),
  height = 6,
  width = 10,
  overwrite = FALSE
)
```

## Make a scatter plot of EOS1 expression model between Asthma and Non-Asthma
```{r}
model_eos1_exp_Screening_Asthma_and_Non_Asthma_combined_kimma_dat <- model_eos1_exp_Screening_Asthma_only_kimma_dat %>%
  subset(comparison != "(Intercept)") %>%
  mutate(comparison = "EOS1_exp_Asthma") %>%
  rbind(
    model_eos1_exp_Screening_Non_Asthma_only_kimma_dat %>%
      subset(comparison != "(Intercept)") %>%
      mutate(comparison = "EOS1_exp_Non_Asthma")
  )

# scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_T2_High_and_T2_Low_status_Screening_lavage_kimma_lm_contrast_plot <- make_contrast_scatterplot(
#   dat = model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lavage_kimma_lm_contrast_dat,
#   comp_x = "(T2-high Non_Asthma)_vs_(T2-high Asthma)",
#   comp_y = "(T2-low Non_Asthma)_vs_(T2-low Asthma)",
#   xlab = "Estimate of Non_Asthma vs Asthma in T2 High",
#   ylab = "Estimate of Non_Asthma vs Asthma in T2 Low",
#   group_x = "T2-High",
#   group_y = "T2-Low",
#   fdr_cutoff = 0.25,
#   top_n = 10
# )
# scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_T2_High_and_T2_Low_status_Screening_lavage_kimma_lm_contrast_plot$plot

scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_eos1_continuous_model_Screening_lavage_kimma_lm_plot <- make_contrast_scatterplot(
  dat = model_eos1_exp_Screening_Asthma_and_Non_Asthma_combined_kimma_dat,
  comp_x = "EOS1_exp_Asthma",
  comp_y = "EOS1_exp_Non_Asthma",
  xlab = "Negative Correlation to EOS1 <----> Positive Correlation to EOS1 \nin Asthma",
  ylab = "Negative Correlation to EOS1 <----> Positive Correlation to EOS1 \nin Non-Asthma",
  group_x = "Asthma",
  group_y = "Non-Asthma",
  fdr_cutoff = 0.2,
  top_n = 25
)

scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_eos1_continuous_model_Screening_lavage_kimma_lm_plot$plot +
  scale_color_manual(
    values = c(
      "Significant in both" = "purple",
      "Significant in Asthma only" = "darkred",
      "Significant in Non-Asthma only" = "orange",
      "Not significant" = "grey"
    )
  ) +
  ggtitle(
    "Scatter Plot of EOS1 Expression Correlation in Asthma vs Non-Asthma FDR <= 0.2"
  )

HoardR::save_figures(
  plot = scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_eos1_continuous_model_Screening_lavage_kimma_lm_plot$plot +
    scale_color_manual(
      values = c(
        "Significant in both" = "purple",
        "Significant in Asthma only" = "darkred",
        "Significant in Non-Asthma only" = "orange",
        "Not significant" = "grey"
      )
    ) +
    ggtitle(
      "Scatter Plot of EOS1 Expression Correlation in Asthma vs Non-Asthma FDR <= 0.2"
    ),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_lavage,
    analysis = "scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_eos1_continuous_model_Screening_lavage_kimma_lm_plot",
    ext = "png"
  ),
  height = 12,
  width = 18,
  overwrite = TRUE
)
```

## Blood Eosinophil Count Analysis

### KIMMA Modeling with Continuous Blood Eosinophil Count

```{r}
model_blood_eos_cont_Screening_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "blood_eos_cont_Screening",
    ext = "RDS"
  ),
  {
    kimma::kmFit(
      dat = dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2,
      patientID = "subjectId",
      libraryID = "libid",
      model = "~ log10_EOS",
      run_lme = FALSE,
      run_lm = TRUE,
      run_contrast = FALSE,
      use_weights = TRUE,
      metrics = TRUE,
      processors = 6
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_blood_eos_cont_Screening_kimma$lm,
  FCgroup = FALSE,
  fdr_cutoff = c(0.05, 0.2)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)


model_blood_eos_cont_Screening_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "blood_eos_cont_Screening_kimma_dat",
    ext = "csv"
  ),
  {
    model_blood_eos_cont_Screening_kimma$lm %>%
      left_join(
        dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
          select(geneName, hgnc_symbol),
        by = c("gene" = "geneName")
      ) %>%
      mutate(comparison = variable) %>%
      select(
        gene,
        hgnc_symbol,
        variable,
        comparison,
        statistic,
        estimate,
        pval,
        FDR,
        df,
        model
      )
  },
  overwrite = TRUE
)
```

### Blood Eosinophil Count Volcano Plot

```{r}
Volcano_plot_blood_eos_cont_Screening_lavage <- volcano_funcxn(
  model_blood_eos_cont_Screening_kimma_dat %>%
    subset(comparison == "log10_EOS") %>%
    mutate(volcano_comp = "log10_EOS"),
  fdr = 0.2,
  n_top_gene = 15,
  add_genes_int = c("None")
)

HoardR::save_figures(
  plot = Volcano_plot_blood_eos_cont_Screening_lavage$log10_EOS +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_volcano_lavage,
    analysis = "Volcano_plot_blood_eos_cont_Screening_lavage",
    ext = "png"
  ),
  height = 6,
  width = 8,
  overwrite = TRUE
)
Volcano_plot_blood_eos_cont_Screening_lavage$log10_EOS + theme_classic_custom
```

### Blood Eosinophil Count Enrichment Analysis

```{r}
Significant_genes_blood_eos_cont_Screening_lavage <- list()
Significant_genes_blood_eos_cont_Screening_lavage[[
  "Positive Correlation to Blood eos"
]] <- model_blood_eos_cont_Screening_kimma_dat %>%
  subset(comparison == "log10_EOS") %>%
  subset(estimate > 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_blood_eos_cont_Screening_lavage[[
  "Negative Correlation to Blood eos"
]] <- model_blood_eos_cont_Screening_kimma_dat %>%
  subset(comparison == "log10_EOS") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_lavage,
    project = projectNumber_lavage,
    analysis = "enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list <- list()
    enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list[[
      "GO_BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_blood_eos_cont_Screening_lavage,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_blood_eos_cont_Screening_lavage,
      collection = "H"
    )
    enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_blood_eos_cont_Screening_lavage,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list
  },
  overwrite = TRUE
)

enrichment_comparison_plot <- list(
  "Group 1:Positive Correlation vs Group 2:Negative Correlation \n Blood Eosinophil Count" = c(
    "Positive Correlation to Blood eos",
    "Negative Correlation to Blood eos"
  )
)

enrichment_plot_GO_BP_Significant_genes_blood_eos_cont_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list$GO_BP,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.05
) +
  theme_classic_custom
print(
  enrichment_plot_GO_BP_Significant_genes_blood_eos_cont_Screening_lavage
)

HoardR::save_figures(
  plot = enrichment_plot_GO_BP_Significant_genes_blood_eos_cont_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_GO_Significant_genes_blood_eos_cont_Screening_lavage",
    ext = "png"
  ),
  height = 10,
  width = 10,
  overwrite = FALSE
)

enrichment_plot_H_Significant_genes_blood_eos_cont_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list$H,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  enrichment_plot_H_Significant_genes_blood_eos_cont_Screening_lavage
)

HoardR::save_figures(
  plot = enrichment_plot_H_Significant_genes_blood_eos_cont_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_H_Significant_genes_blood_eos_cont_Screening_lavage",
    ext = "png"
  ),
  height = 6,
  width = 10,
  overwrite = FALSE
)


enrichment_plot_muppits_Significant_genes_blood_eos_cont_Screening_lavage <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_blood_eos_cont_Screening_lavage_GO_H_muppits_list$muppits,
  comparisons = enrichment_comparison_plot,
  fdr_cutoff = 0.2
) +
  theme_classic_custom
print(
  enrichment_plot_muppits_Significant_genes_blood_eos_cont_Screening_lavage
)

HoardR::save_figures(
  plot = enrichment_plot_muppits_Significant_genes_blood_eos_cont_Screening_lavage,
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Enrichment_muppits_Significant_genes_blood_eos_cont_Screening_lavage",
    ext = "png"
  ),
  height = 6,
  width = 10,
  overwrite = FALSE
)
```

#### Combined plot for volcano and muppits enrichment for blood eosinophil count continuous model
```{r}
Volcano_plot_blood_eos_cont_Screening_lavage$log10_EOS +
  theme_classic_custom +
  labs(
    title = "Negative Correlation <-----> Positive Correlation to Blood Eosinophil Count"
  ) +
  enrichment_plot_muppits_Significant_genes_blood_eos_cont_Screening_lavage +
  plot_layout(ncol = 1, heights = c(1, 1.2), widths = c(1.5, 1))

HoardR::save_figures(
  plot = Volcano_plot_blood_eos_cont_Screening_lavage$log10_EOS +
    theme_classic_custom +
    labs(
      title = "Negative Correlation <-----> Positive Correlation to Blood Eosinophil Count"
    ) +
    enrichment_plot_muppits_Significant_genes_blood_eos_cont_Screening_lavage +
    plot_layout(ncol = 1, heights = c(1, 1.2)),
  paths = HoardR::make_path(
    project = projectNumber_lavage,
    dir = dirPlots_enrichment_lavage,
    analysis = "Volcano_and_muppits_Enrichment_Significant_genes_blood_eos_cont_Screening_lavage",
    ext = "png"
  ),
  height = 8,
  width = 10,
  overwrite = TRUE
)
```

# Advanced Enrichment and Module Analysis

## madRich Module Annotation Analysis

```{r}
library(madRich)

madRich_module_df <- select(
  Enrichment_Module_Visit_model_wo_unknown_GO,
  colnames(Enrichment_Module_Visit_model_wo_unknown_H)
) %>%
  rbind(Enrichment_Module_Visit_model_wo_unknown_H) %>%
  rbind(Enrichment_Module_Visit_model_wo_unknown_KEGG)

file.tmp <- file.path(
  dirresults_lavage,
  paste(projectNumber_lavage, "_madRich_enrichment_height.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  madrich_height <- chooseCutHeight(
    df = madRich_module_df,
    enrich_method = "hypergeometric",
    ID = "SYMBOL",
    collections = c("C5", "H", "C2"),
    subcollections = c("C5" = "GO:BP", "C2" = "CP:KEGG_LEGACY"),
    hclust_heights = c(0.1, 0.3, 0.5, 0.7, 0.9, 1),
    fdr_cutoff = 0.2
  )
  saveRDS(madrich_height, file.tmp)
} else {
  madrich_height <- readRDS(file.tmp)
}
file.tmp <- file.path(
  dirresults_lavage,
  paste(projectNumber_lavage, "_madRich_enrichment_clustersets.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  madrich_res <- clusterSets(
    df = madRich_module_df,
    enrich_method = "hypergeometric",
    ID = "SYMBOL",
    collections = c("C5", "H", "C2"),
    subcollections = c("C5" = "GO:BP", "C2" = "CP:KEGG_LEGACY"),
    hclust_height = c(0.9),
    fdr_cutoff = 0.2
  )
  saveRDS(madrich_res, file.tmp)
} else {
  madrich_res <- readRDS(file.tmp)
}


file.tmp <- file.path(
  dirresults_lavage,
  paste(
    projectNumber_lavage,
    "_madRich_enrichment_height_module_3_4_18.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  madRich_module_df_3_4_18 <- madRich_module_df %>%
    subset(group %in% c("03", "04", "18"))
  madrich_height_3_4_18 <- chooseCutHeight(
    df = madRich_module_df_3_4_18,
    enrich_method = "hypergeometric",
    ID = "SYMBOL",
    collections = c("C5", "H", "C2"),
    subcollections = c("C5" = "GO:BP", "C2" = "CP:KEGG_LEGACY"),
    hclust_heights = c(0.3, 0.5, 0.9, 1),
    fdr_cutoff = 0.2
  )
  saveRDS(madrich_height_3_4_18, file.tmp)
} else {
  madrich_height_3_4_18 <- readRDS(file.tmp)
}
file.tmp <- file.path(
  dirresults_lavage,
  paste(
    projectNumber_lavage,
    "_madRich_enrichment_clustersets_3_4_18.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  madrich_res_3_4_18 <- clusterSets(
    df = madRich_module_df_3_4_18,
    enrich_method = "hypergeometric",
    ID = "SYMBOL",
    collections = c("C5", "H", "C2"),
    subcollections = c("C5" = "GO:BP", "C2" = "CP:KEGG_LEGACY"),
    hclust_height = c(0.9),
    fdr_cutoff = 0.2
  )
  saveRDS(madrich_res_3_4_18, file.tmp)
} else {
  madrich_res_3_4_18 <- readRDS(file.tmp)
}

madrich_res_3 <- clusterSets(
  df = madRich_module_df_3_4_18,
  enrich_method = "hypergeometric",
  ID = "SYMBOL",
  collections = c("C5", "H", "C2"),
  subcollections = c("C5" = "GO:BP", "C2" = "CP:KEGG_LEGACY"),
  hclust_height = c(0.9),
  group_name = "03",
  fdr_cutoff = 0.2
)
madRich::clusterScatter(cluster_result = madrich_res_3)
madRich::clusterTreemap(cluster_result = madrich_res_3)
```

# Clinical Visualization and T2 Gene Analysis

## Color Scheme Configuration for Clinical Variables
```{r}
color_EOS <-
  c(
    paletteer::paletteer_d("khroma::highcontrast") %>% as.character(),
    "grey"
  ) %>%
  setNames(
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      dplyr::arrange(EOS_cat) %>%
      dplyr::pull(EOS_cat) %>%
      unique()
  )


color_FENO <-
  c(paletteer::paletteer_d("unikn::pal_signal") %>% as.character(), "grey") %>%
  setNames(
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      dplyr::arrange(FENO_cat) %>%
      dplyr::pull(FENO_cat) %>%
      unique()
  )

color_Allergy_sens <-
  c(
    paletteer::paletteer_d("nbapalettes::suns_retro") %>% as.character(),
    "grey"
  ) %>%
  setNames(
    dat_voom_lavage_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      dplyr::arrange(Allergy_Sensitization_cat) %>%
      dplyr::pull(Allergy_Sensitization_cat) %>%
      unique()
  )
```



## T2 Signature Gene Expression Analysis

### Create Plots Directory and Data Preparation

```{r}
dir.create(file.path(dirPlots_lavage, "T2_Plots"))
dirPlots_lavageT2 <- file.path(dirPlots_lavage, "T2_Plots")

dat_lavage_wo_unknown_w_meta <- dat_voom_lavage_wo_unknown$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  pivot_longer(!libid, names_to = "geneName", values_to = "geneExp") %>%
  left_join(dat_voom_lavage_wo_unknown$targets) %>%
  left_join(dat_voom_lavage_wo_unknown$genes %>% select(geneName, hgnc_symbol))

library(ggbeeswarm)

dat_lavage_wo_unknown_w_meta %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  ggplot(aes(x = visit, y = geneExp, color = visit)) +
  geom_violin() +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec) +
  facet_wrap(~hgnc_symbol)


box_plot_screening_3_genes <- dat_lavage_wo_unknown_w_meta %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = visit, y = geneExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +

  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "red"
  ) +
  scale_color_manual(values = visit_col_vec) +
  facet_wrap(~hgnc_symbol, scales = "free")


box_plot_screening_3_genes <- dat_lavage_wo_unknown_w_meta %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = visit, y = geneExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +

  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "red"
  ) +
  scale_color_manual(values = visit_col_vec) +
  facet_wrap(~hgnc_symbol, scales = "free")


# Violin plot split by the 3 EOS category
Boxplot_screening_3_genes_EOS_cat <- dat_lavage_wo_unknown_w_meta %>%
  subset(!is.na(EOS_cat)) %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = EOS_cat, y = geneExp, color = EOS_cat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_EOS) +
  facet_wrap(~hgnc_symbol, scales = "free")

Boxplot_screening_3_genes_FENO_cat <- dat_lavage_wo_unknown_w_meta %>%
  subset(!is.na(FENO_cat)) %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = FENO_cat, y = geneExp, color = FENO_cat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_FENO) +
  facet_wrap(~hgnc_symbol, scales = "free")

#Make a density histogram
density_plot_screening_3_genes <- ggplot(
  dat_lavage_wo_unknown_w_meta %>%
    subset(
      geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
    ) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp, color = visit, fill = visit)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25) +
  scale_color_manual(values = visit_col_vec) +
  scale_fill_manual(values = visit_col_vec) +
  theme_classic_custom +
  facet_wrap(~hgnc_symbol, scales = "free")

box_plot_screening_3_genes +
  density_plot_screening_3_genes +
  Boxplot_screening_3_genes_EOS_cat +
  plot_layout(ncol = 1, guides = "collect")

ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_Violin_density_plot_Screening_visit1.pdf"
  ),
  plot = box_plot_screening_3_genes +
    density_plot_screening_3_genes +
    Boxplot_screening_3_genes_EOS_cat +
    plot_layout(ncol = 1, guides = "collect"),
  height = 15,
  width = 10
)


Violin_plot_screening_3_genes_EOS_cat <- dat_lavage_wo_unknown_w_meta %>%
  subset(!is.na(EOS_cat)) %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = EOS_cat, y = geneExp, color = EOS_cat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_EOS) +
  facet_wrap(~hgnc_symbol, scales = "free")

ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_boxplot_Screening_by_EOS_Category.pdf"
  ),
  plot = Violin_plot_screening_3_genes_EOS_cat,
  height = 5,
  width = 8
)


Violin_plot_screening_3_genes_FENO_cat <- dat_lavage_wo_unknown_w_meta %>%
  subset(!is.na(FENO_cat)) %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = FENO_cat, y = geneExp, color = FENO_cat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_FENO) +
  facet_wrap(~hgnc_symbol, scales = "free")

ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_boxplot_Screening_by_FENO_Category.pdf"
  ),
  plot = Violin_plot_screening_3_genes_FENO_cat,
  height = 5,
  width = 8
)


Violin_plot_screening_3_genes_Allergy_Sensitization_cat <- dat_lavage_wo_unknown_w_meta %>%
  subset(!is.na(Allergy_Sensitization_cat)) %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(
    x = Allergy_Sensitization_cat,
    y = geneExp,
    color = Allergy_Sensitization_cat
  )) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_Allergy_sens) +
  facet_wrap(~hgnc_symbol, scales = "free")

ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_boxplot_Screening_by_Allergy_Sensitization_Category.pdf"
  ),
  plot = Violin_plot_screening_3_genes_Allergy_Sensitization_cat,
  height = 5,
  width = 8
)
```

## Clinical Correlation Analysis

### Prepare T2 Gene-Clinical Data Matrix

```{r}
postn_serbinb2_clca1_w_clinical_screening_lavage <- dat_voom_lavage_wo_unknown$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  select(libid, ENSG00000133110, ENSG00000197632, ENSG00000016490) %>%
  left_join(
    dat_voom_lavage_wo_unknown$targets %>%
      select(
        libid,
        visit,
        subjectId,
        `Eosinophil Count (V0)`,
        `FeNO (V0)`,
        `Allergen Sensitizations per AST (V0)`
      )
  ) %>%
  subset(visit == "Screening") %>%
  rename(
    "POSTN" = "ENSG00000133110",
    "SERPINB2" = "ENSG00000197632",
    "CLCA1" = "ENSG00000016490"
  ) %>%
  select(
    subjectId,
    POSTN,
    SERPINB2,
    CLCA1,
    `Eosinophil Count (V0)`,
    `FeNO (V0)`,
    `Allergen Sensitizations per AST (V0)`
  )
postn_serbinb2_clca1_w_clinical_screening_lavage_mat <- postn_serbinb2_clca1_w_clinical_screening_lavage %>%
  rownames_to_column("t") %>%
  select(-c(t)) %>%
  column_to_rownames("subjectId") %>%
  as.matrix()

postn_serbinb2_clca1_w_clinical_screening_lavage_scale <- scale(
  postn_serbinb2_clca1_w_clinical_screening_lavage_mat
)
```

### Create Ranked Data for Clinical Ordering

```{r}
ranks_data_lavage <- postn_serbinb2_clca1_w_clinical_screening_lavage %>%
  mutate(across(
    c(POSTN:`Allergen Sensitizations per AST (V0)`),
    dense_rank
  )) %>%
  mutate(Rank_Means = rowMeans(.[, 2:7], na.rm = TRUE)) %>%
  arrange(Rank_Means) %>%
  select(subjectId, POSTN:`Allergen Sensitizations per AST (V0)`) %>%
  column_to_rownames("subjectId")
rn_lavage <- rownames(ranks_data_lavage)

Rank_data_lavage_meta <- postn_serbinb2_clca1_w_clinical_screening_lavage %>%
  mutate(did = subjectId) %>%
  rownames_to_column("c") %>%
  select(-c("c")) %>%
  column_to_rownames("did")
Rank_data_lavage_meta <- Rank_data_lavage_meta[rn_lavage, ]

Rank_data_meta_lavage_mat <- Rank_data_lavage_meta %>%
  select(-c(subjectId)) %>%
  as.matrix()


pdf(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Ranked_heatmap_Screening.pdf"
  ),
  height = 15,
  width = 15
)
ComplexHeatmap::Heatmap(
  postn_serbinb2_clca1_w_clinical_screening_lavage_scale[rn_lavage, ],
  name = "GeneExp",
  row_gap = unit(0, "mm"),
  column_gap = unit(0, "mm"),
  column_names_rot = 45,
  border = TRUE, #col = color_fun,
  cluster_rows = FALSE,
  cluster_columns = FALSE,
  layer_fun = function(j, i, x, y, width, height, fill) {
    # since grid.text can also be vectorized
    grid.text(
      sprintf(
        "%.1f",
        pindex(
          postn_serbinb2_clca1_w_clinical_screening_lavage_mat[rn_lavage, ],
          i,
          j
        )
      ),
      x,
      y,
      gp = gpar(fontsize = 10)
    )
  }
)
invisible(dev.off())
```


# Make a correlation heatmap of EOS, FENO, SERPINB2, CLCA1 and POSTN
```{r}
cor_matrix_T2_lavage <- dat_lavage_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(
    hgnc_symbol,
    geneExp,
    `Eosinophil Count (V0)`,
    `FeNO (V0)`,
    `Allergen Sensitizations per AST (V0)`,
    `Allergen Sensitizations per sIgE (V0)`
  ) %>%
  subset(hgnc_symbol %in% c("SERPINB2", "CLCA1", "POSTN")) %>%
  pivot_wider(names_from = hgnc_symbol, values_from = geneExp) %>%
  select(
    `Eosinophil Count (V0)`,
    `FeNO (V0)`,
    `Allergen Sensitizations per AST (V0)`,
    CLCA1,
    POSTN,
    SERPINB2
  ) %>%
  rename(
    "EOS" = "Eosinophil Count (V0)",
    "FENO" = "FeNO (V0)",
    "Allergy_Sensitization_per_AST" = "Allergen Sensitizations per AST (V0)"
  ) %>%
  cor(use = "pairwise.complete.obs")
colors_grad_heatmap <- colorRampPalette(c("blue", "white", "red"))(100)
breaks_heatmap <- seq(-1, 1, length.out = 101)


# Plot correlation heatmap
pdf(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_heatmap_Screening.pdf"
  ),
  height = 10,
  width = 15
)
pheatmap(
  cor_matrix_T2_lavage,
  color = colors_grad_heatmap,
  breaks = breaks_heatmap,
  display_numbers = TRUE, # <-- show numbers in cells
  number_format = "%.3f",
  clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean",
  main = "Correlation Heatmap",
  fontsize_number = 10
)
invisible(dev.off())


cor_matrix_T2_lavage_log10 <- dat_lavage_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(
    hgnc_symbol,
    geneExp,
    `Allergen Sensitizations per AST (V0)`,
    `Allergen Sensitizations per sIgE (V0)`,
    log10_EOS,
    log10_FENO
  ) %>%
  subset(hgnc_symbol %in% c("SERPINB2", "CLCA1", "POSTN")) %>%
  pivot_wider(names_from = hgnc_symbol, values_from = geneExp) %>%
  select(
    `Allergen Sensitizations per AST (V0)`,
    log10_EOS,
    log10_FENO,
    CLCA1,
    POSTN,
    SERPINB2
  ) %>%
  rename(
    "Allergy_Sensitization_per_AST" = "Allergen Sensitizations per AST (V0)"
  ) %>%
  cor(use = "pairwise.complete.obs")
colors_grad_heatmap <- colorRampPalette(c("blue", "white", "red"))(100)
breaks_heatmap <- seq(-1, 1, length.out = 101)


# Plot correlation heatmap
pdf(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_heatmap_Screening_log10_EOS_FENO.pdf"
  ),
  height = 8,
  width = 10
)
pheatmap(
  cor_matrix_T2_lavage_log10,
  color = colors_grad_heatmap,
  breaks = breaks_heatmap,
  display_numbers = TRUE, # <-- show numbers in cells
  number_format = "%.3f",
  clustering_distance_rows = "euclidean",
  clustering_distance_cols = "euclidean",
  main = "Correlation Heatmap",
  fontsize_number = 10
)
invisible(dev.off())
```

# Make regression plots for each of the comparisons
```{r}
correlation_regression_T2_lavage <- dat_lavage_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(
    subjectId,
    hgnc_symbol,
    geneExp,
    `Eosinophil Count (V0)`,
    `FeNO (V0)`,
    `Allergen Sensitizations per AST (V0)`,
    `Allergen Sensitizations per sIgE (V0)`,
    log10_EOS,
    log10_FENO
  ) %>%
  subset(hgnc_symbol %in% c("SERPINB2", "CLCA1", "POSTN")) %>%
  pivot_wider(names_from = hgnc_symbol, values_from = geneExp) %>%
  select(
    subjectId,
    `Eosinophil Count (V0)`,
    `FeNO (V0)`,
    `Allergen Sensitizations per AST (V0)`,
    log10_EOS,
    log10_FENO,
    CLCA1,
    POSTN,
    SERPINB2
  ) %>%
  distinct() %>%
  rename(
    "EOS" = "Eosinophil Count (V0)",
    "FENO" = "FeNO (V0)",
    "Allergy_Sensitization_per_AST" = "Allergen Sensitizations per AST (V0)"
  )

vars_T2 <- c(
  "EOS",
  "FENO",
  "Allergy_Sensitization_per_AST",
  "CLCA1",
  "POSTN",
  "SERPINB2"
)
corr_long_T2_lavage <- cor_matrix_T2_lavage %>%
  as.data.frame() %>%
  rownames_to_column("var") %>%
  pivot_longer(-var, names_to = "var2", values_to = "correlation")
# make all combinations (pairs)
pairs_T2_lavage <- expand.grid(
  x = vars_T2,
  y = vars_T2,
  stringsAsFactors = FALSE
) %>%
  filter(x != y) # remove same-vs-same

# optional: if you only want each pair once (not mirrored), keep lower triangle:
pairs_T2_lavage <- pairs_T2_lavage %>%
  mutate(x_rank = match(x, vars_T2), y_rank = match(y, vars_T2)) %>%
  filter(x_rank < y_rank) %>%
  left_join(corr_long_T2_lavage, by = c("x" = "var", "y" = "var2"))
# generate plots for each pair
correlation_regression_T2_lavage_plots <- apply(
  pairs_T2_lavage,
  1,
  function(row) {
    xvar <- row[["x"]]
    yvar <- row[["y"]]
    r_val <- round(as.numeric(row[["correlation"]]), 2)
    ggplot(correlation_regression_T2_lavage, aes_string(x = xvar, y = yvar)) +
      geom_point() +
      geom_smooth(method = "lm", se = TRUE, color = "red") +
      labs(
        x = xvar,
        y = yvar,
        title = paste0(yvar, " vs ", xvar, "\n(r = ", r_val, ")")
      ) +
      theme_classic_custom +
      theme(base_size = 10)
  }
)

# combine with patchwork
correlation_regression_T2_lavage_combined_plot <- wrap_plots(
  correlation_regression_T2_lavage_plots,
  ncol = 5
)

correlation_regression_T2_lavage_combined_plot

ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_regreesion_Screening.pdf"
  ),
  plot = correlation_regression_T2_combined_plot,
  width = 17,
  height = 10
)


vars_T2_log10 <- c(
  "Allergy_Sensitization_per_AST",
  "CLCA1",
  "POSTN",
  "SERPINB2",
  "log10_EOS",
  "log10_FENO"
)
corr_long_T2_lavage_log10 <- cor_matrix_T2_lavage_log10 %>%
  as.data.frame() %>%
  rownames_to_column("var") %>%
  pivot_longer(-var, names_to = "var2", values_to = "correlation")
# make all combinations (pairs)
pairs_T2_lavage_log10 <- expand.grid(
  x = vars_T2_log10,
  y = vars_T2_log10,
  stringsAsFactors = FALSE
) %>%
  filter(x != y) # remove same-vs-same

# optional: if you only want each pair once (not mirrored), keep lower triangle:
pairs_T2_lavage_log10 <- pairs_T2_lavage_log10 %>%
  mutate(x_rank = match(x, vars_T2_log10), y_rank = match(y, vars_T2_log10)) %>%
  filter(x_rank < y_rank) %>%
  left_join(corr_long_T2_lavage_log10, by = c("x" = "var", "y" = "var2"))
# generate plots for each pair
correlation_regression_T2_lavage_log10_plots <- apply(
  pairs_T2_lavage_log10,
  1,
  function(row) {
    xvar <- row[["x"]]
    yvar <- row[["y"]]
    r_val <- round(as.numeric(row[["correlation"]]), 2)
    ggplot(correlation_regression_T2_lavage, aes_string(x = xvar, y = yvar)) +
      geom_point() +
      geom_smooth(method = "lm", se = TRUE, color = "red") +
      labs(
        x = xvar,
        y = yvar,
        title = paste0(yvar, " vs ", xvar, "\n(r = ", r_val, ")")
      ) +
      theme_classic_custom +
      theme(base_size = 10)
  }
)

# combine with patchwork
correlation_regression_T2_lavage_log10_combined_plot <- wrap_plots(
  correlation_regression_T2_lavage_log10_plots,
  ncol = 5
)

correlation_regression_T2_lavage_log10_combined_plot

ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_log10_EOS_FENO_regression_Screening.pdf"
  ),
  plot = correlation_regression_T2_lavage_log10_combined_plot,
  width = 17,
  height = 10
)
```

# Make regression plots for each of the comparisons split by the 3 EOS categories
```{r}
correlation_regression_T2_EOS_cat_lavage <- dat_lavage_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(
    subjectId,
    hgnc_symbol,
    geneExp,
    log10_EOS,
    `Eosinophil Count (V0)`,
    EOS_cat
  ) %>%
  subset(hgnc_symbol %in% c("SERPINB2", "CLCA1", "POSTN")) %>%
  pivot_wider(names_from = hgnc_symbol, values_from = geneExp) %>%
  select(
    subjectId,
    log10_EOS,
    `Eosinophil Count (V0)`,
    EOS_cat,
    CLCA1,
    POSTN,
    SERPINB2
  ) %>%
  distinct() %>%
  rename("EOS" = "Eosinophil Count (V0)") %>%
  subset(!is.na(EOS_cat))


eso_3_gene_scatter_plot_list_eos_together_lavage <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  eso_3_gene_scatter_plot_list_eos_together_lavage[[paste(
    "EOS vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_EOS_cat_lavage,
    aes_string(x = "EOS", y = gene_, color = "EOS_cat")
  ) +
    geom_point(size = 3) +
    geom_smooth(method = "lm", se = TRUE) +
    labs(x = "EOS", y = gene_, title = paste0("EOS", " vs ", gene_)) +
    stat_cor(method = "pearson") +
    scale_color_manual(values = color_EOS) +
    theme_classic_custom
}

eso_3_gene_scatter_plot_list_eos_together_lavage_combined <- wrap_plots(
  eso_3_gene_scatter_plot_list_eos_together_lavage,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_eos_3_gene_by_cat_Screening.pdf"
  ),
  plot = eso_3_gene_scatter_plot_list_eos_together_lavage_combined,
  width = 10,
  height = 10
)


eso_3_gene_scatter_plot_list_eos_together_log10_EOS_lavage <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  eso_3_gene_scatter_plot_list_eos_together_log10_EOS_lavage[[paste(
    "EOS vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_EOS_cat_lavage,
    aes_string(x = "log10_EOS", y = gene_, color = "EOS_cat")
  ) +
    geom_point(size = 3) +
    geom_smooth(method = "lm", se = TRUE) +
    labs(
      x = "log10_EOS",
      y = gene_,
      title = paste0("log10_EOS", " vs ", gene_)
    ) +
    stat_cor(method = "pearson") +
    scale_color_manual(values = color_EOS) +
    theme_classic_custom
}

eso_3_gene_scatter_plot_list_eos_together_log10_EOS_lavage_combined <- wrap_plots(
  eso_3_gene_scatter_plot_list_eos_together_log10_EOS_lavage,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_eos_3_gene_by_cat_Screening_log10_EOS.pdf"
  ),
  plot = eso_3_gene_scatter_plot_list_eos_together_log10_EOS_lavage_combined,
  width = 10,
  height = 10
)
```


# Make regression plots for each of the comparisons split by the 3 FENO categories
```{r}
correlation_regression_T2_FENO_cat_lavage <- dat_lavage_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(subjectId, hgnc_symbol, geneExp, log10_FENO, FENO_cat) %>%
  subset(hgnc_symbol %in% c("SERPINB2", "CLCA1", "POSTN")) %>%
  pivot_wider(names_from = hgnc_symbol, values_from = geneExp) %>%
  select(subjectId, log10_FENO, FENO_cat, CLCA1, POSTN, SERPINB2) %>%
  distinct() %>%
  #rename("FENO"="FENOinophil Count (V0)") %>%
  subset(!is.na(FENO_cat))


FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO_lavage <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO_lavage[[paste(
    "FENO vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_FENO_cat_lavage,
    aes_string(x = "log10_FENO", y = gene_, color = "FENO_cat")
  ) +
    geom_point(size = 3) +
    geom_smooth(method = "lm", se = TRUE) +
    labs(
      x = "log10_FENO",
      y = gene_,
      title = paste0("log10_FENO", " vs ", gene_)
    ) +
    stat_cor(method = "pearson") +
    scale_color_manual(values = color_FENO) +
    theme_classic_custom
}

FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO_lavage_combined <- wrap_plots(
  FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO_lavage,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_FENO_3_gene_by_cat_Screening_log10_FENO.pdf"
  ),
  plot = FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO_lavage_combined,
  width = 10,
  height = 10
)
```


# Make regression plots for each of the comparisons split by the 3 Allergy_Sensitization categories
```{r}
correlation_regression_T2_Allergy_Sensitization_cat_lavage <- dat_lavage_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(
    subjectId,
    hgnc_symbol,
    geneExp,
    `Allergen Sensitizations per AST (V0)`,
    Allergy_Sensitization_cat
  ) %>%
  subset(hgnc_symbol %in% c("SERPINB2", "CLCA1", "POSTN")) %>%
  pivot_wider(names_from = hgnc_symbol, values_from = geneExp) %>%
  select(
    subjectId,
    `Allergen Sensitizations per AST (V0)`,
    Allergy_Sensitization_cat,
    CLCA1,
    POSTN,
    SERPINB2
  ) %>%
  distinct() %>%
  #rename("Allergy_Sensitization"="Allergy_Sensitizationinophil Count (V0)") %>%
  subset(!is.na(Allergy_Sensitization_cat))


Allergy_Sensitization_3_gene_scatter_plot_list_together_lavage <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  Allergy_Sensitization_3_gene_scatter_plot_list_together_lavage[[paste(
    "Allergy_Sensitization vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_Allergy_Sensitization_cat_lavage,
    aes_string(
      x = "`Allergen Sensitizations per AST (V0)`",
      y = gene_,
      color = "Allergy_Sensitization_cat"
    )
  ) +
    geom_point(size = 3) +
    geom_smooth(method = "lm", se = TRUE) +
    labs(
      x = "`Allergen Sensitizations per AST (V0)`",
      y = gene_,
      title = paste0("`Allergen Sensitizations per AST (V0)`", " vs ", gene_)
    ) +
    stat_cor(method = "pearson") +
    scale_color_manual(values = color_Allergy_sens) +
    theme_classic_custom
}

Allergy_Sensitization_3_gene_scatter_plot_list_together_lavage_combined <- wrap_plots(
  Allergy_Sensitization_3_gene_scatter_plot_list_together_lavage,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlots_lavageT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_Allergy_Sensitization_3_gene_by_cat_Screening_log10_Allergy_Sensitization.pdf"
  ),
  plot = Allergy_Sensitization_3_gene_scatter_plot_list_together_lavage_combined,
  width = 10,
  height = 10
)
```



# Look at the correlation of MUPPITS eos1 and unma11 modules to EOS expression
```{r}
muppits_mods <- read_csv(file.path(
  dirDataraw,
  "MUPPITS_geneSetsByAssignedCell.csv"
))

Muppits_modules_mean_expression_lavage <- muppits_mods %>%
  rename("geneName" = "ensemblID") %>%
  #Combine count and module data
  dplyr::select(geneName, geneSet) %>%
  dplyr::left_join(
    tibble::rownames_to_column(
      as.data.frame(dat_voom_lavage_wo_unknown$E),
      "geneName"
    ),
    by = "geneName"
  ) %>%
  #Calculate mean by module
  dplyr::group_by(geneSet) %>%
  dplyr::summarise_if(is.numeric, mean, na.rm = TRUE)


Muppits_modules_mean_expression_lavage_w_meta <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "Muppits_expression_data",
    ext = "csv"
  ),
  {
    Muppits_modules_mean_expression_lavage %>%
      pivot_longer(!geneSet, names_to = "libid", values_to = "geneExp") %>%
      left_join(dat_voom_lavage_wo_unknown$targets) %>%
      mutate(
        EOS_cat = factor(
          EOS_cat,
          levels = c(
            "EOS_Less_than_150",
            "EOS_between_150_300",
            "EOS_Greater_than_300"
          )
        ),
        FENO_cat = factor(
          FENO_cat,
          levels = c(
            "FENO_Less_than_20",
            "FENO_between_20_35",
            "FENO_Greater_than_35"
          )
        ),
        Allergy_Sensitization_cat = factor(
          Allergy_Sensitization_cat,
          levels = c(
            "Allergy_Sensitization_Less_than_1",
            "Allergy_Sensitization_between_2_5",
            "Allergy_Sensitization_Greater_than_5"
          )
        )
      )
  },
  overwrite = FALSE
)

# Muppits_modules_mean_expression_lavage_w_meta <- Muppits_modules_mean_expression_lavage_w_meta %>%
#   mutate(
#     EOS_cat = factor(
#       EOS_cat,
#       levels = c(
#         "EOS_Less_than_150",
#         "EOS_between_150_300",
#         "EOS_Greater_than_300"
#       )
#     ),
#     FENO_cat = factor(
#       FENO_cat,
#       levels = c(
#         "FENO_Less_than_20",
#         "FENO_between_20_35",
#         "FENO_Greater_than_35"
#       )
#     ),
#     Allergy_Sensitization_cat = factor(
#       Allergy_Sensitization_cat,
#       levels = c(
#         "Allergy_Sensitization_Less_than_1",
#         "Allergy_Sensitization_between_2_5",
#         "Allergy_Sensitization_Greater_than_5"
#       )
#     )
#   )
```

# plot log10EOS vs eos1 mod and unma11

```{r}
plot_eos1_vs_EOS_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "eos1" & !is.na(EOS_cat) & visit == "Screening"),
  aes_string(x = "log10_EOS", y = "geneExp", color = "EOS_cat")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "log10_EOS",
    y = "EOS1 Expression",
    title = paste0("log10_EOS", " vs ", "EOS1", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  theme_classic_custom

ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "eos1" & !is.na(EOS_cat) & visit == "Screening"),
  aes_string(x = "log10_EOS", y = "geneExp")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "log10_EOS",
    y = "EOS1 Expression",
    title = paste0("log10_EOS", " vs ", "EOS1", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  theme_classic_custom


plot_unma11_vs_EOS_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "unma11" & !is.na(EOS_cat) & visit == "Screening"),
  aes_string(x = "log10_EOS", y = "geneExp", color = "EOS_cat")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "log10_EOS",
    y = "unma11 Expression",
    title = paste0("log10_EOS", " vs ", "unma11", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  theme_classic_custom


plot_eos1_vs_EOS_lavage +
  plot_unma11_vs_EOS_lavage +
  plot_layout(ncol = 1, guides = "collect")


HoardR::save_figures(
  plot = plot_eos1_vs_EOS_lavage +
    plot_unma11_vs_EOS_lavage +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Contrast_Scatter,
    project = projectNumber_lavage,
    analysis = "EOS_vs_Muppits_EOS1_and_unma11",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)


plot_eos1_vs_EOS_lavage_asthma_status <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "eos1" & !is.na(EOS_cat) & visit == "Screening") %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ),
  aes_string(x = "log10_EOS", y = "geneExp", color = "EOS_cat")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    x = "log10_EOS",
    y = "EOS1 Expression",
    title = paste0("log10_EOS", " vs ", "EOS1", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  facet_grid(~asthma_status) +
  theme_classic_custom


plot_unma11_vs_EOS_lavage_asthma_status <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "unma11" & !is.na(EOS_cat) & visit == "Screening") %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ),
  aes_string(x = "log10_EOS", y = "geneExp", color = "EOS_cat")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = FALSE) +
  labs(
    x = "log10_EOS",
    y = "unma11 Expression",
    title = paste0("log10_EOS", " vs ", "unma11", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  facet_grid(~asthma_status) +
  theme_classic_custom


plot_eos1_vs_EOS_lavage_asthma_status +
  plot_unma11_vs_EOS_lavage_asthma_status +
  plot_layout(ncol = 1, guides = "collect")


HoardR::save_figures(
  plot = plot_eos1_vs_EOS_lavage_asthma_status +
    plot_unma11_vs_EOS_lavage_asthma_status +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Contrast_Scatter,
    project = projectNumber_lavage,
    analysis = "EOS_vs_Muppits_EOS1_and_unma11_by_asthma_status",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 14,
  height = 7
)


boxplot_plot_screening_muppits_eos1_EOS_cat_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "eos1" & !is.na(EOS_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = EOS_cat, y = geneExp, color = EOS_cat)) +
  geom_boxplot(outlier.shape = NA) +
  ggbeeswarm::geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "EOS Category", y = "eos1 Expression", title = "eos1 Nasal Lavage") +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_EOS)


boxplot_plot_screening_muppits_unma11_EOS_cat_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "unma11" & !is.na(EOS_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = EOS_cat, y = geneExp, color = EOS_cat)) +
  geom_boxplot(outlier.shape = NA) +
  ggbeeswarm::geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "EOS Category",
    y = "unma11 Expression",
    title = "unma11 Nasal Lavage"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_EOS)


boxplot_plot_screening_muppits_eos1_EOS_cat_lavage +
  boxplot_plot_screening_muppits_unma11_EOS_cat_lavage +
  plot_layout(ncol = 2, guides = "collect")


HoardR::save_figures(
  plot = boxplot_plot_screening_muppits_eos1_EOS_cat_lavage +
    boxplot_plot_screening_muppits_unma11_EOS_cat_lavage +
    plot_layout(ncol = 2, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Boxplot,
    project = projectNumber_lavage,
    analysis = "EOS_cat_Muppits_eos1_and_unma11",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)


boxplot_plot_screening_muppits_eos1_EOS_cat_lavage_asthma_status <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "eos1" & !is.na(EOS_cat)) %>%
  mutate(
    asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = EOS_cat, y = geneExp, color = EOS_cat)) +
  geom_boxplot(outlier.shape = NA) +
  ggbeeswarm::geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "EOS Category", y = "eos1 Expression", title = "eos1 Nasal Lavage") +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  facet_grid(~asthma_status) +
  scale_color_manual(values = color_EOS)


boxplot_plot_screening_muppits_unma11_EOS_cat_lavage_asthma_status <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "unma11" & !is.na(EOS_cat)) %>%
  mutate(
    asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = EOS_cat, y = geneExp, color = EOS_cat)) +
  geom_boxplot(outlier.shape = NA) +
  ggbeeswarm::geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "EOS Category",
    y = "unma11 Expression",
    title = "unma11 Nasal Lavage"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  facet_grid(~asthma_status) +
  scale_color_manual(values = color_EOS)


boxplot_plot_screening_muppits_eos1_EOS_cat_lavage_asthma_status +
  boxplot_plot_screening_muppits_unma11_EOS_cat_lavage_asthma_status +
  plot_layout(ncol = 2, guides = "collect")


HoardR::save_figures(
  plot = boxplot_plot_screening_muppits_eos1_EOS_cat_lavage_asthma_status +
    boxplot_plot_screening_muppits_unma11_EOS_cat_lavage_asthma_status +
    plot_layout(ncol = 2, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Boxplot,
    project = projectNumber_lavage,
    analysis = "EOS_cat_Muppits_eos1_and_unma11_by_asthma_status",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 13,
  height = 7
)


boxplot_plot_screening_muppits_eos1_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet %in% c("eos1", "unma11") & !is.na(EOS_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = visit, y = geneExp)) +
  geom_boxplot(outlier.shape = NA) +
  ggbeeswarm::geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(x = "Visit", y = "Mean Expression", title = "Nasal Lavage") +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "red"
  ) +
  facet_wrap(~geneSet, nrow = 1)

HoardR::save_figures(
  plot = boxplot_plot_screening_muppits_eos1_lavage,
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Boxplot,
    project = projectNumber_lavage,
    analysis = "Muppits_eos1_and_unma11_boxplot",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 5,
  height = 5
)

# make the desnity plot
#Make the density plot just eos1 expression
##Make a density histogram
density_plot_screening_muppits_eos1_unma11_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(EOS_cat)) %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25, color = "black", fill = "grey") +
  theme_classic_custom +
  facet_wrap(~geneSet, scales = "free", nrow = 1)

#Make the density plot just eos1 expression coloered by category
#Make a density histogram colored by EOS category
density_plot_screening_muppits_eos1_unma11_eos_cat_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(EOS_cat)) %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp, color = EOS_cat, fill = EOS_cat)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25) +
  scale_color_manual(values = color_EOS) +
  scale_fill_manual(values = color_EOS) +
  theme_classic_custom +
  facet_wrap(~geneSet, scales = "free", nrow = 1)


HoardR::save_figures(
  plot = density_plot_screening_muppits_eos1_unma11_lavage +
    density_plot_screening_muppits_eos1_unma11_eos_cat_lavage +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Density,
    project = projectNumber_lavage,
    analysis = "Muppits_eos1_and_unma11_density_by_EOS_cat",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 14,
  height = 7
)


##Make a density histogram split by asthma status
density_plot_screening_muppits_eos1_unma11_lavage_asthma_status <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(EOS_cat)) %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25, color = "black", fill = "grey") +
  theme_classic_custom +
  facet_wrap(geneSet ~ asthma_status, scales = "free", nrow = 1)

#Make the density plot just eos1 expression coloered by category
#Make a density histogram split by asthma status and colored by EOS category
density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_asthma_status <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(EOS_cat)) %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp, color = EOS_cat, fill = EOS_cat)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25) +
  scale_color_manual(values = color_EOS) +
  scale_fill_manual(values = color_EOS) +
  theme_classic_custom +
  facet_wrap(geneSet ~ asthma_status, scales = "free", nrow = 1)


HoardR::save_figures(
  plot = density_plot_screening_muppits_eos1_unma11_lavage_asthma_status +
    density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_asthma_status +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Density,
    project = projectNumber_lavage,
    analysis = "Muppits_eos1_and_unma11_density_by_EOS_cat_and_asthma_status",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 14,
  height = 7
)
```

## Make a combined plot just for eos1 with the boxplot, scatter plot, and density plot together
```{r}
density_plot_screening_muppits_eos1_eos_cat_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(EOS_cat)) %>%
    mutate(
      asthma_status = factor(asthma_status, levels = c("Non_Asthma", "Asthma"))
    ) %>%
    subset(geneSet %in% c("eos1")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp, color = EOS_cat, fill = EOS_cat)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25) +
  scale_color_manual(values = color_EOS) +
  scale_fill_manual(values = color_EOS) +
  theme_classic_custom +
  xlab("EOS1 Expression") +
  ggtitle("EOS1 Nasal Lavage") +
  theme(legend.position = "none")


density_plot_screening_muppits_eos1_eos_cat_lavage /
  (boxplot_plot_screening_muppits_eos1_EOS_cat_lavage +
    xlab("EOS1 Expression") +
    theme(legend.position = "none") +
    plot_eos1_vs_EOS_lavage +
    xlab("Blood eos counts (log10)")) +
  plot_layout(guides = "collect")

HoardR::save_figures(
  plot = density_plot_screening_muppits_eos1_eos_cat_lavage /
    (boxplot_plot_screening_muppits_eos1_EOS_cat_lavage +
      xlab("EOS1 Expression") +
      theme(legend.position = "none") +
      plot_eos1_vs_EOS_lavage +
      xlab("Blood eos counts (log10)")) +
    plot_layout(guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots,
    project = projectNumber_lavage,
    analysis = "EOS1_combined_density_boxplot_scatter",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 8
)
```

# Segregrate the donors based on eos1 in low cat
```{r}
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening <- Muppits_modules_mean_expression_lavage_w_meta %>%
  #subset(!is.na(EOS_cat)) %>%
  subset(geneSet %in% c("eos1")) %>%
  subset(visit %in% c("Screening")) %>%
  select(subjectId, geneExp, geneSet, EOS_cat) %>%
  distinct()


# Compute mean and SD per EOS group
stats_per_group_eos1_lavage <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening %>%
  subset(EOS_cat == "EOS_Less_than_150") %>%
  #group_by(EOS_cat) %>%
  summarise(
    n = n(),
    mean_val = mean(geneExp),
    sd_val = sd(geneExp),
    se_val = sd_val / sqrt(n),
    .groups = "drop"
  ) %>%
  mutate(
    mean_plus_1SD = mean_val + sd_val,
    mean_plus_2SD = mean_val + 2 * sd_val,
    mean_plus_1SE = mean_val + se_val,
    mean_plus_2SE = mean_val + 2 * se_val
  )

# Join stats back to df
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening %>%
  #left_join(stats_per_group_eos1_lavage, by = "EOS_cat") %>%
  mutate(
    T2_mean = ifelse(
      geneExp <= stats_per_group_eos1_lavage$mean_val,
      "T2-low",
      "T2-high"
    ),
    T2_1SD = ifelse(
      geneExp <= stats_per_group_eos1_lavage$mean_plus_1SD,
      "T2-low",
      "T2-high"
    ),
    T2_2SD = ifelse(
      geneExp <= stats_per_group_eos1_lavage$mean_plus_2SD,
      "T2-low",
      "T2-high"
    ),
    T2_1SE = ifelse(
      geneExp <= stats_per_group_eos1_lavage$mean_plus_1SE,
      "T2-low",
      "T2-high"
    ),
    T2_2SE = ifelse(
      geneExp <= stats_per_group_eos1_lavage$mean_plus_2SE,
      "T2-low",
      "T2-high"
    ),
  )
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_lavage,
    analysis = "Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_lavage",
    ext = "csv"
  ),
  {
    Muppits_modules_mean_expression_lavage_w_meta_eos1_screening %>%
      #left_join(stats_per_group_eos1_lavage, by = "EOS_cat") %>%
      mutate(
        T2_mean = ifelse(
          geneExp <= stats_per_group_eos1_lavage$mean_val,
          "T2-low",
          "T2-high"
        ),
        T2_1SD = ifelse(
          geneExp <= stats_per_group_eos1_lavage$mean_plus_1SD,
          "T2-low",
          "T2-high"
        ),
        T2_2SD = ifelse(
          geneExp <= stats_per_group_eos1_lavage$mean_plus_2SD,
          "T2-low",
          "T2-high"
        ),
        T2_1SE = ifelse(
          geneExp <= stats_per_group_eos1_lavage$mean_plus_1SE,
          "T2-low",
          "T2-high"
        ),
        T2_2SE = ifelse(
          geneExp <= stats_per_group_eos1_lavage$mean_plus_2SE,
          "T2-low",
          "T2-high"
        ),
      )
  },
  overwrite = FALSE
)


# write_csv(
#   Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified %>%
#     select(subjectId, geneExp, EOS_cat, T2_mean, T2_1SD),
#   "Data_clean/Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_lavage.csv"
# )

# Count samples per EOS group for each cutpoint
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified %>%
  pivot_longer(
    cols = starts_with("T2_"),
    names_to = "cutpoint",
    values_to = "T2_status"
  ) %>%
  group_by(EOS_cat, cutpoint, T2_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = T2_status, values_from = n, values_fill = 0)


cutpoints_values_lavage <- stats_per_group_eos1_lavage %>%
  select(mean_val, mean_plus_1SD, mean_plus_2SD, mean_plus_1SE, mean_plus_2SE)


# Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table %>%
#   left_join(cutpoints_values, by = "EOS_cat")
#
# Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table

# Density plot with vertical lines for each EOS group
density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified %>%
    subset(!is.na(EOS_cat)),
  aes(x = geneExp, fill = EOS_cat, color = EOS_cat)
) +
  geom_density(alpha = 0.25, lwd = 1) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_val),
    color = "black",
    linetype = "solid",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_1SD),
    color = "blue",
    linetype = "dashed",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_2SD),
    color = "red",
    linetype = "dotted",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_1SE),
    color = "darkgreen",
    linetype = "longdash",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_2SE),
    color = "darkorchid4",
    linetype = "dotdash",
    linewidth = 1
  ) +
  #facet_wrap(~EOS_cat, scales = "fixed") +
  labs(
    x = "Gene Expression (eos1)",
    y = "Density",
    title = "Density of eos1 Expression by EOS Category",
    subtitle = "Vertical lines: black=mean, blue=mean+1SD, red=mean+2SD, green=mean+1SE,purple=mean+2SE"
  ) +
  scale_color_manual(values = color_EOS) +
  scale_fill_manual(values = color_EOS) +
  theme_classic_custom


HoardR::save_figures(
  plot = density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low,
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Density,
    project = projectNumber_lavage,
    analysis = "Density_plot_eos1_expression_by_EOS_cat_with_cutpoints",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)


density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low_only_mean_1sd <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified %>%
    subset(!is.na(EOS_cat)),
  aes(x = geneExp, fill = EOS_cat, color = EOS_cat)
) +
  geom_density(alpha = 0.25, lwd = 1) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_1SD),
    color = "blue",
    linetype = "dashed",
    linewidth = 1
  ) +
  #facet_wrap(~EOS_cat, scales = "fixed") +
  labs(
    x = "Gene Expression (eos1)",
    y = "Density",
    title = "Density of eos1 Expression by EOS Category",
    subtitle = "Vertical lines: blue=mean+1SD"
  ) +
  scale_color_manual(values = color_EOS) +
  scale_fill_manual(values = color_EOS) +
  theme_classic_custom

HoardR::save_figures(
  plot = density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low_only_mean_1sd,
  paths = HoardR::make_path(
    dir = dirT2plots_lavage$plots_Density,
    project = projectNumber_lavage,
    analysis = "Density_plot_eos1_expression_by_EOS_cat_only_mean_1sd_cutpoint",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)
```


# Segregrate the donors based on eos1 in low cat and Asthma
```{r}
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_asthma <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(!is.na(EOS_cat)) %>%
  subset(geneSet %in% c("eos1")) %>%
  subset(visit %in% c("Screening")) %>%
  select(subjectId, geneExp, geneSet, EOS_cat, asthma_status) %>%
  distinct()


# Compute mean and SD per EOS group
stats_per_group_eos1_lavage_asthma <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_asthma %>%
  subset(EOS_cat == "EOS_Less_than_150" & asthma_status == "Asthma") %>%
  #group_by(EOS_cat) %>%
  summarise(
    n = n(),
    mean_val = mean(geneExp),
    sd_val = sd(geneExp),
    se_val = sd_val / sqrt(n),
    .groups = "drop"
  ) %>%
  mutate(
    mean_plus_1SD = mean_val + sd_val,
    mean_plus_2SD = mean_val + 2 * sd_val,
    mean_plus_1SE = mean_val + se_val,
    mean_plus_2SE = mean_val + 2 * se_val
  )

# Join stats back to df
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_asthma <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_asthma %>%
  #left_join(stats_per_group_eos1_lavage, by = "EOS_cat") %>%
  mutate(
    T2_mean = ifelse(
      geneExp <= stats_per_group_eos1_lavage_asthma$mean_val,
      "T2-low",
      "T2-high"
    ),
    T2_1SD = ifelse(
      geneExp <= stats_per_group_eos1_lavage_asthma$mean_plus_1SD,
      "T2-low",
      "T2-high"
    ),
    T2_2SD = ifelse(
      geneExp <= stats_per_group_eos1_lavage_asthma$mean_plus_2SD,
      "T2-low",
      "T2-high"
    ),
    T2_1SE = ifelse(
      geneExp <= stats_per_group_eos1_lavage_asthma$mean_plus_1SE,
      "T2-low",
      "T2-high"
    ),
    T2_2SE = ifelse(
      geneExp <= stats_per_group_eos1_lavage_asthma$mean_plus_2SE,
      "T2-low",
      "T2-high"
    ),
  )

# Count samples per EOS group for each cutpoint
Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_asthma_summary_table <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_asthma %>%
  pivot_longer(
    cols = starts_with("T2_"),
    names_to = "cutpoint",
    values_to = "T2_status"
  ) %>%
  group_by(EOS_cat, asthma_status, cutpoint, T2_status) %>%
  summarise(n = n(), .groups = "drop") %>%
  pivot_wider(names_from = T2_status, values_from = n, values_fill = 0)

cutpoints_values_lavage_asthma <- stats_per_group_eos1_lavage_asthma %>%
  select(mean_val, mean_plus_1SD, mean_plus_2SD, mean_plus_1SE, mean_plus_2SE)


# Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table <- Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table %>%
#   left_join(cutpoints_values, by = "EOS_cat")
#
# Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_summary_table

# Density plot with vertical lines for each EOS group
density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low_asthma <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta_eos1_screening_classified_asthma,
  aes(x = geneExp, fill = EOS_cat, color = EOS_cat)
) +
  geom_density(alpha = 0.25, lwd = 1) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_val),
    color = "black",
    linetype = "solid",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_1SD),
    color = "blue",
    linetype = "dashed",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_2SD),
    color = "red",
    linetype = "dotted",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_1SE),
    color = "darkgreen",
    linetype = "longdash",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_lavage,
    aes(xintercept = mean_plus_2SE),
    color = "darkorchid4",
    linetype = "dotdash",
    linewidth = 1
  ) +
  facet_wrap(~asthma_status, scales = "fixed") +
  labs(
    x = "Gene Expression (eos1)",
    y = "Density",
    title = "Density of eos1 Expression by EOS Category",
    subtitle = "Vertical lines: black=mean, blue=mean+1SD, red=mean+2SD, green=mean+1SE,purple=mean+2SE"
  ) +
  scale_color_manual(values = color_EOS) +
  scale_fill_manual(values = color_EOS) +
  theme_classic_custom

ggsave(
  file.path(
    dirPlots_lavageT2,
    "density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low_asthma.pdf"
  ),
  plot = density_plot_screening_muppits_eos1_unma11_eos_cat_lavage_split_T2_high_low_asthma,
  width = 10,
  height = 7
)
```


# Make a correlation heatmap with lower triange shwoing the correlation and upper triangle showing scatterplot of correlation
```{r}
df_num <- Muppits_modules_mean_expression_lavage_w_meta %>%
  select(
    visit,
    geneSet,
    geneExp,
    log10_EOS,
    log10_FENO,
    `Allergen Sensitizations per AST (V0)`
  ) %>%
  subset(geneSet %in% c("eos1", "unma11") & visit == "Screening") %>%
  pivot_wider(names_from = geneSet, values_from = geneExp) %>%
  select(-c(visit)) %>%
  rename(Allergen_Sensitization = `Allergen Sensitizations per AST (V0)`)


library(GGally)
library(dplyr)

# Example: keep only numeric columns

p_to_stars <- function(p) {
  if (is.na(p)) {
    return("")
  }
  if (p < 0.001) {
    return("***")
  }
  if (p < 0.01) {
    return("**")
  }
  if (p < 0.05) {
    return("*")
  }
  return("")
}

# Lower panel: colored tile (fill = correlation) + text = r (and stars)
my_cor_tile <- function(data, mapping, digits = 2, size = 5, ...) {
  x <- eval_data_col(data, mapping$x)
  y <- eval_data_col(data, mapping$y)

  # correlation & p-value
  ct <- tryCatch(cor.test(x, y), error = function(e) NULL)
  if (!is.null(ct)) {
    r <- as.numeric(ct$estimate)
    p <- ct$p.value
  } else {
    r <- NA_real_
    p <- NA_real_
  }

  lbl <- if (is.na(r)) {
    ""
  } else {
    paste0(formatC(r, format = "f", digits = digits), "\n", p_to_stars(p))
  }

  # create a single-tile panel that fills the cell
  ggplot(data.frame(x = 1, y = 1, r = r), aes(x = x, y = y, fill = r)) +
    geom_tile(width = 1.0, height = 1.0) +
    geom_text(aes(label = lbl), size = size) +
    scale_fill_gradient2(
      low = "blue",
      mid = "white",
      high = "red",
      limits = c(-1, 1),
      na.value = "grey90"
    ) +
    scale_x_continuous(expand = c(0, 0), limits = c(0.5, 1.5)) +
    scale_y_continuous(expand = c(0, 0), limits = c(0.5, 1.5)) +
    theme_void() +
    theme(legend.position = "none")
}

# Upper panel: scatter + regression line (with CI)
my_scatter_lm <- function(data, mapping, ...) {
  ggplot(data = data, mapping = mapping) +
    geom_point(alpha = 0.7, size = 1.5) +
    geom_smooth(method = "lm", se = TRUE, size = 0.6) +
    theme_minimal()
}

# Build the ggpairs plot
p <- ggpairs(
  df_num,
  lower = list(continuous = my_cor_tile),
  upper = list(continuous = my_scatter_lm),
  diag = list(continuous = wrap("densityDiag")),
  axisLabels = "show"
)


pdf(
  file.path(
    dirPlots_lavageT2,
    "Correlation_heatmap_with_scatterplot_between_muppits_and_clinical.pdf"
  ),
  height = 15,
  width = 15
)
print(p)
invisible(dev.off())
```


# plot log10FENO vs eos1 mod and unma11

```{r}
plot_eos1_vs_FENO_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "eos1" & !is.na(FENO_cat) & visit == "Screening"),
  aes_string(x = "log10_FENO", y = "geneExp", color = "FENO_cat")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "log10_FENO",
    y = "FENO1 Expression",
    title = paste0("log10_FENO", " vs ", "FENO1", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_FENO) +
  theme_classic_custom


plot_unma11_vs_FENO_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(geneSet == "unma11" & !is.na(FENO_cat) & visit == "Screening"),
  aes_string(x = "log10_FENO", y = "geneExp", color = "FENO_cat")
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "log10_FENO",
    y = "unma11 Expression",
    title = paste0("log10_FENO", " vs ", "unma11", " Nasal Lavage")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_FENO) +
  theme_classic_custom


plot_eos1_vs_FENO_lavage +
  plot_unma11_vs_FENO_lavage +
  plot_layout(ncol = 1, guides = "collect")


ggsave(
  file.path(dirPlots_lavageT2, "FENO_vs_Muppits_FENO1_and_unma11.pdf"),
  plot = plot_eos1_vs_FENO_lavage +
    plot_unma11_vs_FENO_lavage +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)


boxplot_plot_screening_muppits_eos1_FENO_cat_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "eos1" & !is.na(FENO_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = FENO_cat, y = geneExp, color = FENO_cat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "FENO Category",
    y = "eos1 Expression",
    title = "eos1 Nasal Lavage"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_FENO)


boxplot_plot_screening_muppits_unma11_FENO_cat_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "unma11" & !is.na(FENO_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = FENO_cat, y = geneExp, color = FENO_cat)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "FENO Category",
    y = "unma11 Expression",
    title = "unma11 Nasal Lavage"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_FENO)


boxplot_plot_screening_muppits_eos1_FENO_cat_lavage +
  boxplot_plot_screening_muppits_unma11_FENO_cat_lavage +
  plot_layout(ncol = 2, guides = "collect")


ggsave(
  file.path(dirPlots_lavageT2, "box_plot_FENO_cat_Muppits_eos1_and_unma11.pdf"),
  plot = boxplot_plot_screening_muppits_eos1_FENO_cat_lavage +
    boxplot_plot_screening_muppits_unma11_FENO_cat_lavage +
    plot_layout(ncol = 2, guides = "collect"),
  width = 10,
  height = 7
)


# make the desnity plot
#Make the density plot just eos1 expression
#Make a density histogram
density_plot_screening_muppits_eos1_unma11_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(FENO_cat)) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25, color = "black", fill = "grey") +
  theme_classic_custom +
  facet_wrap(~geneSet, scales = "free")

#Make the density plot just eos1 expression coloered by category
#Make a density histogram
density_plot_screening_muppits_eos1_unma11_FENO_cat_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(FENO_cat)) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp, color = FENO_cat, fill = FENO_cat)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25) +
  scale_color_manual(values = color_FENO) +
  scale_fill_manual(values = color_FENO) +
  theme_classic_custom +
  facet_wrap(~geneSet, scales = "free")


ggsave(
  file.path(
    dirPlots_lavageT2,
    "density_plot_screening_muppits_eos1_unma11_FENO_cat.pdf"
  ),
  plot = density_plot_screening_muppits_eos1_unma11_lavage +
    density_plot_screening_muppits_eos1_unma11_FENO_cat_lavage +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)
```


# plot log10Allergy_Sensitization vs eos1 mod and unma11

```{r}
plot_eos1_vs_Allergy_Sensitization_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(
      geneSet == "eos1" &
        !is.na(Allergy_Sensitization_cat) &
        visit == "Screening"
    ),
  aes_string(
    x = "`Allergen Sensitizations per AST (V0)`",
    y = "geneExp",
    color = "Allergy_Sensitization_cat"
  )
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "Allergen Sensitizations per AST (V0)",
    y = "eos1 Expression",
    title = paste0(
      "Allergen Sensitizations per AST (V0)",
      " vs ",
      "eos1",
      " Nasal Lavage"
    )
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_Allergy_sens) +
  theme_classic_custom


plot_unma11_vs_Allergy_Sensitization_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(
      geneSet == "unma11" &
        !is.na(Allergy_Sensitization_cat) &
        visit == "Screening"
    ),
  aes_string(
    x = "`Allergen Sensitizations per AST (V0)`",
    y = "geneExp",
    color = "Allergy_Sensitization_cat"
  )
) +
  geom_point(size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(
    x = "Allergen Sensitizations per AST (V0)",
    y = "unma11 Expression",
    title = paste0(
      "Allergen Sensitizations per AST (V0)",
      " vs ",
      "unma11",
      " Nasal Lavage"
    )
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_Allergy_sens) +
  theme_classic_custom


plot_eos1_vs_Allergy_Sensitization_lavage +
  plot_unma11_vs_Allergy_Sensitization_lavage +
  plot_layout(ncol = 1, guides = "collect")


ggsave(
  file.path(
    dirPlots_lavageT2,
    "Allergy_Sensitization_vs_Muppits_eos1_and_unma11.pdf"
  ),
  plot = plot_eos1_vs_Allergy_Sensitization_lavage +
    plot_unma11_vs_Allergy_Sensitization_lavage +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)


boxplot_plot_screening_muppits_eos1_Allergy_Sensitization_cat_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "eos1" & !is.na(Allergy_Sensitization_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(
    x = Allergy_Sensitization_cat,
    y = geneExp,
    color = Allergy_Sensitization_cat
  )) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "Allergy_Sensitization Category",
    y = "eos1 Expression",
    title = "eos1 Nasal Lavage"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_Allergy_sens)


boxplot_plot_screening_muppits_unma11_Allergy_Sensitization_cat_lavage <- Muppits_modules_mean_expression_lavage_w_meta %>%
  subset(geneSet == "unma11" & !is.na(Allergy_Sensitization_cat)) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(
    x = Allergy_Sensitization_cat,
    y = geneExp,
    color = Allergy_Sensitization_cat
  )) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(
    x = "Allergy_Sensitization Category",
    y = "unma11 Expression",
    title = "unma11 Nasal Lavage"
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    shape = 18,
    size = 3,
    color = "black"
  ) +
  scale_color_manual(values = color_Allergy_sens)


boxplot_plot_screening_muppits_eos1_Allergy_Sensitization_cat_lavage +
  boxplot_plot_screening_muppits_unma11_Allergy_Sensitization_cat_lavage +
  plot_layout(ncol = 2, guides = "collect")


ggsave(
  file.path(
    dirPlots_lavageT2,
    "box_plot_Allergy_Sensitization_cat_Muppits_eos1_and_unma11.pdf"
  ),
  plot = boxplot_plot_screening_muppits_eos1_Allergy_Sensitization_cat_lavage +
    boxplot_plot_screening_muppits_unma11_Allergy_Sensitization_cat_lavage +
    plot_layout(ncol = 2, guides = "collect"),
  width = 10,
  height = 7
)


#Make a density histogram
density_plot_screening_muppits_eos1_unma11_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(Allergy_Sensitization_cat)) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(x = geneExp)
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25, color = "black", fill = "grey") +
  theme_classic_custom +
  facet_wrap(~geneSet, scales = "free")

#Make the density plot just eos1 expression coloered by category
#Make a density histogram
density_plot_screening_muppits_eos1_unma11_Allergy_Sensitization_cat_lavage <- ggplot(
  Muppits_modules_mean_expression_lavage_w_meta %>%
    subset(!is.na(Allergy_Sensitization_cat)) %>%
    subset(geneSet %in% c("eos1", "unma11")) %>%
    subset(visit %in% c("Screening")),
  aes(
    x = geneExp,
    color = Allergy_Sensitization_cat,
    fill = Allergy_Sensitization_cat
  )
) +
  #geom_histogram(aes(y = ..density..),alpha=0.5) +
  geom_density(lwd = 1, alpha = 0.25) +
  scale_color_manual(values = color_Allergy_sens) +
  scale_fill_manual(values = color_Allergy_sens) +
  theme_classic_custom +
  facet_wrap(~geneSet, scales = "free")


ggsave(
  file.path(
    dirPlots_lavageT2,
    "density_plot_screening_muppits_eos1_unma11_Allergy_Sensitization_cat.pdf"
  ),
  plot = density_plot_screening_muppits_eos1_unma11_lavage +
    density_plot_screening_muppits_eos1_unma11_Allergy_Sensitization_cat_lavage +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)
```


#Muppits Module decoherence in Lavage
```{r}
source("../../Bioinformatics_R_scripts/calculate_module_coherence.R")

P639_1_muppits_module_coherence <-
  calculate_module_coherence(
    mods = muppits_mods %>%
      rename("geneName" = "ensemblID", "module.char" = "geneSet"),
    dat = dat_voom_lavage_wo_unknown,
    mods_title = "MUPPITS Modules",
    dat_title = "P639-1 Lavage",
    r_cutoff = 0.3,
    p_cutoff = 0.05
  )
```