2.P640-1_Gene_Analysis.qmd

---
title: "P640-1: Cheetah RNA-seq analysis"
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

Load 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(pullString)
library(mgcv)
library(ComplexHeatmap)
#library(ica)
library(gtools)
library(readr)
#pak::pak("basilinb/HoardR")
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

## Theme, Utilities, and Color Palettes

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, linewidth = 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))
```

## Functions and Plot Utilities



## 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_swab <- "P640-1"
dirBoxBase <-
  file.path("~", "Library", "CloudStorage", "Box-Box")
dirRoot <-
  file.path(
    dirBoxBase,
    "Altman_Lab", # Box version
    # "~", "Documents", "Projects_local", # local version
    projectNumber_swab
  )
dirPlots <- file.path("Figs")
dirresults <- file.path("Results")
dirDataclean <- file.path("Data_clean")
dirDataraw <- file.path("Data_raw")

subdirs_used <- c(
  "Volcano",
  "Contrast_Scatter",
  "Boxplot",
  "Enrichment",
  "Heatmap",
  "Density"
)

dirPlots_swab <- file.path("Figs", projectNumber_swab)
dirresults_swab <- file.path("Results", projectNumber_swab)
dirPlots_volcano_swab <- file.path(dirPlots_swab, "Volcano")
dirPlots_contrast_scatter_swab <- file.path(
  dirPlots_swab,
  "Contrast_Scatter"
)
dirPlots_boxplot_swab <- file.path(
  dirPlots_swab,
  "Boxplot"
)

dirPlots_enrichment_swab <- file.path(
  dirPlots_swab,
  "Enrichment"
)

dirPlots_heatmap_swab <- file.path(
  dirPlots_swab,
  "Heatmap"
)


if (!dir.exists(dirPlots_swab)) {
  dir.create(dirPlots_swab)
}
if (!dir.exists(dirresults_swab)) {
  dir.create(dirresults_swab)
}


if (!dir.exists(dirPlots_contrast_scatter_swab)) {
  dir.create(dirPlots_contrast_scatter_swab)
}
if (!dir.exists(dirPlots_boxplot_swab)) {
  dir.create(dirPlots_boxplot_swab)
}
if (!dir.exists(dirPlots_enrichment_swab)) {
  dir.create(dirPlots_enrichment_swab)
}
if (!dir.exists(dirPlots_heatmap_swab)) {
  dir.create(dirPlots_heatmap_swab)
}
if (!dir.exists(dirresults_swab)) {
  dir.create(dirresults_swab)
}

dirT2plots_swab <- HoardR::init_experiment_dirs(
  parent_dir = dirPlots_swab,
  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_swab <- read_csv(file.path(
  dirDataclean,
  paste(projectNumber_swab, "_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"))
file.tmp <- file.path(
  dirDataclean,
  paste(projectNumber_swab, "_dat_voom_normalized_w_clinical.RDS", sep = "")
)
dat_voom <- readRDS(file = file.tmp)

dat_voom$targets <- dat_voom$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(
    EOS_cat = factor(
      EOS_cat,
      levels = c(
        "EOS_Less_than_150",
        "EOS_between_150_300",
        "EOS_Greater_than_300"
      )
    )
  ) %>%
  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 the ? and Unknown visit
```{r wite_load__swab_data_wo_unknown_w_clinical}
dat_voom_wo_unknown <- HoardR::cache(
  paths = make_path(
    project = projectNumber_swab,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical",
    ext = "rds",
    dir = dirDataclean
  ),
  {
    libs_wo_unknown <- dat_voom$targets %>%
      subset(visit %notin% c("?", "Unknown")) %>%
      pull(libid)
    dat_voom_wo_unknown <- dat_voom[, libs_wo_unknown]

    dat_voom_wo_unknown$targets <- dat_voom_wo_unknown$targets %>%
      mutate(
        visit = ifelse(visit == "Screening/Enrollment", "Screening", visit)
      ) %>%
      mutate(
        visit = factor(
          visit,
          levels = c("Screening", "Cold A Clinic", "Cold B Home", "Visit 1")
        )
      )
    dat_voom_wo_unknown
  },
  overwrite = FALSE
)
```

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

Swab_samples_with_T2 <- dat_voom_wo_unknown$targets %>%
  janitor::clean_names() %>%
  left_join(
    T2_HIgh_Low_classification_by_eos1_meand_1sd_swab %>%
      select(subject_id, gene_exp, t2_1sd)
  )

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


# Differential Gene Expression Analysis

## Visit Comparison Analysis
```{r}
model_visit_wo_unknown_kimma <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    analysis = "visit_median_cv_comparison_wo_unknown_kimma",
    ext = "RDS",
    dir = dirresults
  ),
  {
    model_visit_wo_unknown_kimma <- kimma::kmFit(
      dat = dat_voom_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
    )
  },
  overwrite = FALSE
)


kimma::summarize_kmFit(
  model_visit_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)
```



# Extract and Process Visit Comparison Results
```{r}
#| label = "visit_model_results"
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_visit_median_cv_comparison_wo_unknown_kimma.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  model_visit_wo_unknown_kimma_results <- model_visit_wo_unknown_kimma$lme.contrast %>%
    left_join(
      dat_voom_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_kimma_results, file.tmp)
} else {
  model_visit_wo_unknown_kimma_results <- read_csv(file.tmp)
}
```


# Make a list of the contrast results

```{r}
fdr_cutoff <- 0.05

model_visit_wo_unknown_kimma_results_list <- model_visit_wo_unknown_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
```{r}
upset(
  fromList(model_visit_wo_unknown_kimma_results_list[c(1, 2, 7, 8, 5, 6)]),
  text.scale = 1.5,
  order.by = "freq",
  mainbar.y.label = "Significant genes\nFDR < 0.05",
  keep.order = TRUE,
  sets = names(model_visit_wo_unknown_kimma_results_list[c(1, 2, 7, 8, 5, 6)])
)
```


```{r}
volcano_funcxn <- function(
  limma_output,
  coeff = NULL,
  fdr = 0.01,
  n_top_gene = 10,
  add_genes_int = NULL
) {
  for (package in c('tidyverse', 'ggplot2', 'ggrepel')) {
    if (!require(package, character.only = T, quietly = T)) {
      install.packages(package)
      library(package, character.only = T)
    }
  }

  if (is.null(coeff)) {
    coeff <- limma_output$volcano_comp %>% unique()
  } else {
    coeff <- coeff
  }

  options(ggrepel.max.overlaps = Inf)
  # Read in file generated from extract pvalue function for model of interest
  p_val <- limma_output
  vol_plt_list <- list()
  # If you want to run function for all coefficients

  # running for loop for each given coefficient
  for (i in coeff) {
    pval_coef <- subset(p_val, volcano_comp == i)
    # getting all hgnc symbol for up regulated genes
    up_reg <- head(
      subset(pval_coef, estimate > 0 & FDR < fdr) %>% arrange(FDR),
      n_top_gene
    ) %>%
      rename(delabel = hgnc_symbol) %>%
      select(gene, delabel)
    # getting all hgnc symbol for down regulated genes
    down_reg <- head(
      subset(pval_coef, estimate < 0 & FDR < fdr) %>% arrange(FDR),
      n_top_gene
    ) %>%
      rename(delabel = hgnc_symbol) %>%
      select(gene, delabel)
    add_gene_reg <- head(
      subset(pval_coef, hgnc_symbol %in% add_genes_int) %>% arrange(FDR),
      n_top_gene
    ) %>%
      rename(delabel = hgnc_symbol) %>%
      select(gene, delabel)
    # making the dataframe to plot
    pval_plot <- left_join(pval_coef, down_reg, by = "gene") %>%
      left_join(., up_reg, by = "gene") %>%
      left_join(., add_gene_reg, by = "gene") %>%
      mutate(delabel_final = coalesce(delabel.x, delabel.y, delabel)) %>%
      mutate(FC.group = ifelse(estimate > 0, "up", "down")) %>%
      mutate(
        FC.group = case_when(
          FDR <= fdr ~ FC.group,
          FDR > fdr & !is.na(delabel_final) ~ "Gene of Interest(NS)",
          FDR > fdr ~ "Not significant"
        )
      ) %>%
      #mutate(FC.group = factor(FC.group,levels = c("Not significant","Gene of Interest","down","up"))) %>%
      arrange(desc(FC.group))
    cols <- c(
      "up" = "red",
      "down" = "blue",
      "Gene of Interest(NS)" = "steelblue4",
      "Not significant" = "grey"
    )
    # running ggplot to make volcano plot
    plt_vol <- ggplot(
      data = pval_plot,
      aes(x = estimate, y = -log10(FDR), col = FC.group)
    ) +
      geom_point(size = 3, alpha = 0.7) +
      theme_classic() +
      geom_text_repel(aes(label = delabel_final), color = "black") +
      scale_color_manual(values = cols) +
      geom_vline(xintercept = 0, col = "black", linetype = 2) + #xlim(c(-6.1,6.7)) + ylim(c(0,14.9)) +
      geom_hline(yintercept = -log10(fdr), col = "black", linetype = 2) +
      ggtitle(i)
    vol_plt_list[[i]] <- plt_vol
  }
  return(vol_plt_list)
}
```

# Volcano plot: Screening vs Cold A (top genes)
```{r}
nasal_swab_visit_volc <- volcano_funcxn(
  model_visit_wo_unknown_kimma_results %>%
    mutate(
      volcano_comp = paste(contrast_ref, " <------> ", contrast_lvl, sep = "")
    ),
  fdr = 0.05,
  n_top_gene = 15
)

nasal_swab_visit_volc_FDR0.25 <- volcano_funcxn(
  model_visit_wo_unknown_kimma_results %>%
    mutate(
      volcano_comp = paste(contrast_ref, " <------> ", contrast_lvl, sep = "")
    ),
  fdr = 0.25,
  n_top_gene = 15
)

nasal_swab_visit_volc$`Screening <------> Cold_A_Clinic` +
  theme_classic_custom +
  ggtitle("Screening <------> Cold A Clinic FDR < 0.05")
nasal_swab_visit_volc$`Screening <------> Cold_B_Home` +
  theme_classic_custom +
  ggtitle("Screening <------> Cold B Home FDR < 0.05")
nasal_swab_visit_volc_FDR0.25$`Cold_A_Clinic <------> Cold_B_Home` +
  theme_classic_custom +
  ggtitle("Cold A Clinic <------> Cold B Home FDR < 0.25")


HoardR::save_figures(
  plot = nasal_swab_visit_volc$`Screening <------> Cold_A_Clinic` +
    theme_classic_custom +
    ggtitle("Screening <------> Cold A Clinic FDR < 0.05"),
  path = HoardR::make_path(
    dir = dirPlots_volcano_swab,
    project = projectNumber_swab,
    analysis = "Volcano_plot_ColdA_vs_Screening",
    ext = "pdf"
  ),
  height = 6,
  width = 7,
  overwrite = FALSE
)

HoardR::save_figures(
  plot = nasal_swab_visit_volc$`Screening <------> Cold_B_Home` +
    theme_classic_custom +
    ggtitle("Screening <------> Cold B Home FDR < 0.05"),
  path = HoardR::make_path(
    dir = dirPlots_volcano_swab,
    project = projectNumber_swab,
    analysis = "Volcano_plot_ColdB_vs_Screening",
    ext = "pdf"
  ),
  height = 6,
  width = 7,
  overwrite = FALSE
)

HoardR::save_figures(
  plot = nasal_swab_visit_volc$`Cold_A_Clinic <------> Cold_B_Home` +
    theme_classic_custom +
    ggtitle("Cold A Clinic <------> Cold B Home FDR < 0.25"),
  path = HoardR::make_path(
    dir = dirPlots_volcano_swab,
    project = projectNumber_swab,
    analysis = "Volcano_plot_ColdA_vs_ColdB",
    ext = "pdf"
  ),
  height = 6,
  width = 7,
  overwrite = FALSE
)


HoardR::save_figures(
  plot = nasal_swab_visit_volc$`Screening <------> Cold_A_Clinic` +
    theme_classic_custom +
    ggtitle("Screening <------> Cold A Clinic FDR < 0.05") +
    nasal_swab_visit_volc$`Screening <------> Cold_B_Home` +
    theme_classic_custom +
    ggtitle("Screening <------> Cold B Home FDR < 0.05") +
    nasal_swab_visit_volc_FDR0.25$`Cold_A_Clinic <------> Cold_B_Home` +
    theme_classic_custom +
    ggtitle("Cold A Clinic <------> Cold B Home FDR < 0.25") +
    plot_layout(ncol = 1, guides = "collect"),
  path = HoardR::make_path(
    dir = dirPlots_volcano_swab,
    project = projectNumber_swab,
    analysis = "Volcano_plot_ColdA_vs_ColdB_vs_Screening_combined",
    ext = "pdf"
  ),
  height = 15,
  width = 8,
  overwrite = TRUE
)
```



# Boxplots of Top Genes
```{r}
dat_wo_unknown_w_meta <- dat_voom_wo_unknown$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  pivot_longer(!libid, names_to = "geneName", values_to = "geneExp") %>%
  left_join(dat_voom_wo_unknown$targets) %>%
  left_join(dat_voom_wo_unknown$genes %>% select(geneName, hgnc_symbol))

top_4_cold_A_genes <- model_visit_wo_unknown_kimma_results %>%
  subset(
    hgnc_symbol %in%
      c("IFIH1", "STAT2", "IFIT3", "IRF7", "OAS3") &
      comparison == "Cold_A_Clinic_vs_Screening"
  )

box_plot_top_4_cold_A_genes <- dat_wo_unknown_w_meta %>%
  subset(geneName %in% top_4_cold_A_genes$gene) %>%
  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 +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(~hgnc_symbol, scales = "free", nrow = 1) +
  ylab("Gene Expression (log2)")

ggsave(
  file.path(dirPlots, "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_kimma_results %>%
  subset(comparison == "Cold_A_Clinic_vs_Screening" & estimate <= c(-1)) %>%
  arrange(FDR) %>%
  head(5)

box_plot_top_4_screening_genes <- dat_wo_unknown_w_meta %>%
  subset(geneName %in% top_4_screening_genes$gene) %>%
  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 +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  facet_wrap(~hgnc_symbol, scales = "free", nrow = 1) +
  ylab("Gene Expression (log2)")

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


ggsave(
  file.path(dirPlots, "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
)
```





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

model_visit_wo_unknown_swab_kimma_results_list[[
  "Up in Cold B Home"
]] <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_B_Home_vs_Screening") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
model_visit_wo_unknown_swab_kimma_results_list[[
  "Up in Screening"
]] <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_B_Home_vs_Screening") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()

model_visit_wo_unknown_swab_kimma_results_list[[
  "Up in Cold B Home"
]] <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_B_Home_vs_Cold_A_Clinic") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
model_visit_wo_unknown_swab_kimma_results_list[[
  "Up in Cold A Clinic"
]] <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_B_Home_vs_Cold_A_Clinic") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()

enrichment_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "enrichment_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_list <- list()

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

    enrichment_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = model_visit_wo_unknown_swab_kimma_results_list,
      collection = "H"
    )

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


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

plot_H_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_swab_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_swab_kimma_results_list_GO_H_muppits_swab
)

HoardR::save_figures(
  plot = plot_H_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "H_enrichment_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits",
    ext = "png"
  ),
  height = 5,
  width = 10,
  overwrite = TRUE
)

plot_GO_BP_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_swab_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_swab_kimma_results_list_GO_H_muppits_swab
)

HoardR::save_figures(
  plot = plot_GO_BP_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "GO_BP_enrichment_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits",
    ext = "png"
  ),
  height = 20,
  width = 15,
  overwrite = FALSE
)


plot_GO_BP_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab_top_15_term <- plot_enrichment_comparisons(
  df = enrichment_model_visit_wo_unknown_swab_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_swab_kimma_results_list_GO_H_muppits_swab_top_15_term
)

HoardR::save_figures(
  plot = plot_GO_BP_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab_top_15_term +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "GO_BP_model_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_top_15_term",
    ext = "png"
  ),
  height = 9,
  width = 15,
  overwrite = FALSE
)


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

HoardR::save_figures(
  plot = plot_muppits_model_enrichment_visit_wo_unknown_swab_kimma_results_list_GO_H_muppits_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "muppits_model_enrichment_visit_wo_unknown_swab_kimma_results",
    ext = "png"
  ),
  height = 5,
  width = 10,
  overwrite = TRUE
)
```

## Old Enrichment Analysis Code (can delete later)
```{r}
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_visit_comparison_wo_unknown_Enrichment_Hallmark_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  Enrichment_Visit_model_wo_unknown_H <- SEARchways::BIGprofiler(
    gene_list = model_visit_wo_unknown_kimma_results_list,
    collection = "H"
  )
  write_csv(Enrichment_Visit_model_wo_unknown_H, file.tmp)
} else {
  Enrichment_Visit_model_wo_unknown_H <- read_csv(file.tmp)
}

file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_visit_comparison_wo_unknown_Enrichment_GO_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  Enrichment_Visit_model_wo_unknown_GO <- SEARchways::BIGprofiler(
    model_visit_wo_unknown_kimma_results_list,
    collection = "C5",
    subcollection = "GO:BP",
    ID = "SYMBOL"
  )
  write_csv(Enrichment_Visit_model_wo_unknown_GO, file.tmp)
} else {
  Enrichment_Visit_model_wo_unknown_GO <- read_csv(file.tmp)
}

Enrichment_Visit_model_wo_unknown_GO <- HoardR::cache(
  make_path(
    project = projectNumber_swab,
    analysis = "visit_comparison_wo_unknown_Enrichment_GO",
    ext = "csv",
    dir = dirresults_swab
  ),
  {
    Enrichment_Visit_model_wo_unknown_GO <- SEARchways::BIGprofiler(
      model_visit_wo_unknown_kimma_results_list,
      collection = "C5",
      subcollection = "GO:BP",
      ID = "SYMBOL"
    )
  },
  overwrite = FALSE
)

plot_enrichment_Visit_model_wo_unknown_H <- BIGpicture::plot_enrich(
  Enrichment_Visit_model_wo_unknown_H
)
save_figures(
  plot_enrichment_Visit_model_wo_unknown_H,
  paths = make_path(
    projectNumber_swab,
    "enrichment_Visit_model_wo_unknown_Hallmark",
    c("pdf"),
    dirPlots_swab
  ),
  width = 18,
  height = 6
)
```

# Plot the Enrichment
```{r}
fdr_cutoff <- 0.2

#Hallmark
dat_enrichment_Visit_model_wo_unknown_H <- Enrichment_Visit_model_wo_unknown_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)),
      "Up in Contrast Lvl",
      ifelse(
        sig & str_detect(group, paste0("^Up in ", contrast_ref)),
        "Up in Contrast Ref",
        NA
      )
    )
  )


plot_enrichment_Visit_model_wo_unknown_H <- ggplot(
  dat_enrichment_Visit_model_wo_unknown_H %>% subset(FDR <= fdr_cutoff),
  aes(x = comparison, y = pathway)
) +
  geom_point(
    aes(size = `k/K`, color = Significance),
    shape = 17, # triangle
    alpha = 0.5
  ) +
  scale_color_manual(
    values = c("Up in Contrast Lvl" = "red", "Up in Contrast Ref" = "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, "Enrichment_Visit_wo_unknown_Hallmark.pdf"),
  plot = plot_enrichment_Visit_model_wo_unknown_H,
  height = 5,
  width = 15
)

#GO
fdr_cutoff <- 0.01
dat_enrichment_Visit_model_wo_unknown_GO <- Enrichment_Visit_model_wo_unknown_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)),
      "Up in Contrast Lvl",
      ifelse(
        sig & str_detect(group, paste0("^Up in ", contrast_ref)),
        "Up in Contrast Ref",
        NA
      )
    )
  )

dat_enrichment_Visit_model_wo_unknown_GO_2 <- Enrichment_Visit_model_wo_unknown_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_GO <- ggplot(
  dat_enrichment_Visit_model_wo_unknown_GO %>%
    subset(FDR <= fdr_cutoff) %>%
    slice_min(order_by = FDR, n = 10),
  aes(x = comparison, y = pathway)
) +
  geom_point(
    aes(size = `k/K`, color = Significance),
    shape = 17, # triangle
    alpha = 0.5
  ) +
  scale_color_manual(
    values = c("Up in Contrast Lvl" = "red", "Up in Contrast Ref" = "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, "Enrichment_Visit_wo_unknown_GO.pdf"),
  plot = plot_enrichment_Visit_model_wo_unknown_GO,
  height = 15,
  width = 25
)


plot_enrichment_Visit_model_wo_unknown_GO_coldA_vs_screening <- ggplot(
  dat_enrichment_Visit_model_wo_unknown_GO_2 %>%
    subset(
      comparison == "Contrast Lvl: Cold_A_Clinic\nContrast Ref: Screening"
    ) %>%
    subset(
      pathway !=
        "GOBP_ADAPTIVE_IMMUNE_RESPONSE_BASED_ON_SOMATIC_RECOMBINATION_OF_IMMUNE_RECEPTORS_BUILT_FROM_IMMUNOGLOBULIN_SUPERFAMILY_DOMAINS"
    ) %>%
    subset(FDR <= fdr_cutoff) %>%
    group_by(Significance) %>%
    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, "Enrichment_Visit_coldA_vs_screening_wo_unknown_GO.pdf"),
  plot = plot_enrichment_Visit_model_wo_unknown_GO_coldA_vs_screening,
  height = 10,
  width = 15
)
```

# Box plot of top genes
```{r}
dat_wo_unknown_w_meta <- dat_voom_wo_unknown$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  pivot_longer(!libid, names_to = "geneName", values_to = "geneExp") %>%
  left_join(dat_voom_wo_unknown$targets) %>%
  left_join(dat_voom_wo_unknown$genes %>% select(geneName, hgnc_symbol))

top_2_cold_A_vs_screening_genes <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_A_Clinic_vs_Screening" & estimate > 1) %>%
  arrange(FDR) %>%
  head(2)

box_plot_top_2_cold_A_vs_screening_genes <- dat_wo_unknown_w_meta %>%
  subset(geneName %in% top_2_cold_A_vs_screening_genes$gene) %>%
  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") +
  ylab("Gene Expression (log2)")

ggsave(
  file.path(dirPlots, "box_plot_top_2_cold_A_vs_screening_genes.pdf"),
  plot = box_plot_top_2_cold_A_vs_screening_genes,
  height = 5,
  width = 7
)


top_2_cold_B_vs_screening_genes <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_B_Home_vs_Screening" & estimate > 1) %>%
  arrange(FDR) %>%
  head(2)

box_plot_top_2_cold_B_vs_screening_genes <- dat_wo_unknown_w_meta %>%
  subset(geneName %in% top_2_cold_B_vs_screening_genes$gene) %>%
  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") +
  ylab("Gene Expression (log2)")

ggsave(
  file.path(dirPlots, "box_plot_top_2_cold_B_vs_screening_genes.pdf"),
  plot = box_plot_top_2_cold_B_vs_screening_genes,
  height = 5,
  width = 7
)


top_2_screening_genes <- model_visit_wo_unknown_kimma_results %>%
  subset(comparison == "Cold_A_Clinic_vs_Screening" & estimate < 1) %>%
  arrange(FDR) %>%
  head(2)

box_plot_top_2_screening_genes <- dat_wo_unknown_w_meta %>%
  subset(geneName %in% top_2_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") +
  ylab("Gene Expression (log2)")

ggsave(
  file.path(dirPlots, "box_plot_top_2_screening_genes.pdf"),
  plot = box_plot_top_2_screening_genes,
  height = 5,
  width = 7
)
```


# WGCNA on Significant Genes (FDR 0.2)

```{r}
#| dependson = visit_model_results
file.tmp <- file.path(
  dirDataclean,
  paste(
    projectNumber_swab,
    "_dat_voom_normalized_wo_unknown_subset_Visit_significant_genes(FDR_0.2).RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  sig_genes_fdr_0.2 <- model_visit_wo_unknown_kimma_results %>%
    subset(FDR <= 0.2) %>%
    pull(gene) %>%
    unique()
  dat_voom_wo_unknown_subset_visit_sig_genes_FDR_0.2 <- dat_voom_wo_unknown[
    sig_genes_fdr_0.2,
  ]

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

# Run WGCNA: Module Detection
```{r}
file.tmp <- file.path(
  dirDataclean,
  paste(projectNumber_swab, "_dat_WGCNA_Visit_significant_genes.RDS", sep = "")
)
if (file.exists(file.tmp)) {
  load(file.tmp)
} else {
  fit_WGCNA_visit_sig_genes_FDR_0.2 <- RNAetc::fit_modules(
    dat = dat_voom_wo_unknown_subset_visit_sig_genes_FDR_0.2
  )

  fit_WGCNA_visit_sig_genes_FDR_0.2$top.plot +
    fit_WGCNA_visit_sig_genes_FDR_0.2$connect.plot

  WGCNA_visit_sig_genes_FDR_0.2 <- RNAetc::make_modules(
    fit_WGCNA_visit_sig_genes_FDR_0.2,
    sft_value = 10,
    nThread = 4,
    mods_mean = TRUE,
    minModuleSize = 50
  )
  #Make WGCNA modules EList

  dat.WGCNA_Visit_significant_genes <- new("EList")
  dat.WGCNA_Visit_significant_genes$E <- as.matrix(
    WGCNA_visit_sig_genes_FDR_0.2$mods.mean %>%
      select(-c(module.char)) %>%
      select(
        dat_voom_wo_unknown_subset_visit_sig_genes_FDR_0.2$targets %>%
          arrange(libid) %>%
          pull(libid)
      )
  )
  dat.WGCNA_Visit_significant_genes$targets <- dat_voom_wo_unknown_subset_visit_sig_genes_FDR_0.2$targets %>%
    arrange(libid)
  dat.WGCNA_Visit_significant_genes$modules <- WGCNA_visit_sig_genes_FDR_0.2$mods

  save(dat.WGCNA_Visit_significant_genes, file = file.tmp)
}

dat.WGCNA_Visit_significant_genes$modules %>%
  count(module.char) %>%
  kable(align = "l", caption = "Modules") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```


# Module-level Differential Expression (KIMMA)
```{r}
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "WGCNA_module_visit_comparison_wo_unknown_kimma.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  model_module_visit_wo_unknown_kimma <- kimma::kmFit(
    dat = dat.WGCNA_Visit_significant_genes,
    patientID = "subjectId",
    libraryID = "libid",
    model = "~ visit + (1|subjectId)",
    run_lme = TRUE,
    run_contrast = TRUE,
    use_weights = FALSE,
    metrics = TRUE,
    contrast_var = "visit",
    processors = 6
  )
  saveRDS(model_module_visit_wo_unknown_kimma, file = file.tmp)
} else {
  model_module_visit_wo_unknown_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_module_visit_wo_unknown_kimma$lme.contrast,
  FCgroup = TRUE,
  fdr_cutoff = c("0.01", "0.05")
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```


# Module Contrast Results
```{r}
#| label = "visit_model_results"
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "Module_visit_comparison_wo_unknown_kimma_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  model_module_visit_wo_unknown_kimma_results <- model_module_visit_wo_unknown_kimma$lme.contrast %>%
    mutate(contrast_ref = gsub(" ", "_", contrast_ref)) %>%
    mutate(contrast_lvl = gsub(" ", "_", contrast_lvl)) %>%
    mutate(comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_")) %>%
    select(gene, comparison, estimate, pval, FDR, contrast_ref, contrast_lvl)
  write_csv(model_module_visit_wo_unknown_kimma_results, file.tmp)
} else {
  model_module_visit_wo_unknown_kimma_results <- read_csv(file.tmp)
}
```


# Module Contrast Results: Lists

```{r}
fdr_cutoff <- 0.01

model_module_visit_wo_unknown_kimma_results_list <- model_module_visit_wo_unknown_kimma_results %>%
  filter(FDR <= fdr_cutoff) %>%
  group_by(comparison) %>%
  summarise(
    up_in_lvl = list(unique(gene[estimate > 0])),
    up_in_ref = list(unique(gene[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()
```

# Module Upset Plot
```{r}
upset(
  fromList(model_module_visit_wo_unknown_kimma_results_list),
  text.scale = 1.5,
  order.by = "freq",
  mainbar.y.label = "Significant Modules\nFDR < 0.01",
  keep.order = TRUE,
  sets = names(model_module_visit_wo_unknown_kimma_results_list)
)

source("../../Bioinformatics_R_scripts/overlapGroups.R")
model_module_visit_wo_unknown_kimma_results_list_overlaps <- overlapGroups(
  model_module_visit_wo_unknown_kimma_results_list
)

model_module_visit_wo_unknown_kimma_results_list_overlaps_list <- list()
for (overlap_ in names(
  model_module_visit_wo_unknown_kimma_results_list_overlaps
)) {
  model_module_visit_wo_unknown_kimma_results_list_overlaps_list[[
    overlap_
  ]] <- attr(
    model_module_visit_wo_unknown_kimma_results_list_overlaps,
    "elements"
  )[model_module_visit_wo_unknown_kimma_results_list_overlaps[[overlap_]]]
}
```

# Module Annotation

# Module Enrichment (Hallmark, GO, KEGG)
```{r}
module_gene_list <- list()

for (mod_ in unique(dat.WGCNA_Visit_significant_genes$modules$module.char)) {
  module_gene_list[[mod_]] <- dat.WGCNA_Visit_significant_genes$modules %>%
    subset(module.char == mod_) %>%
    pull(hgnc_symbol)
}

file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_Module_visit_comparison_wo_unknown_Enrichment_Hallmark_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  Enrichment_Module_Visit_model_wo_unknown_H <- SEARchways::BIGprofiler(
    gene_list = module_gene_list,
    collection = "H"
  )
  write_csv(Enrichment_Module_Visit_model_wo_unknown_H, file.tmp)
} else {
  Enrichment_Module_Visit_model_wo_unknown_H <- read_csv(file.tmp)
}

file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_Module_visit_comparison_wo_unknown_Enrichment_GO_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  Enrichment_Module_Visit_model_wo_unknown_GO <- SEARchways::BIGprofiler(
    module_gene_list,
    collection = "C5",
    subcollection = "GO:BP",
    ID = "SYMBOL"
  )
  write_csv(Enrichment_Module_Visit_model_wo_unknown_GO, file.tmp)
} else {
  Enrichment_Module_Visit_model_wo_unknown_GO <- read_csv(file.tmp)
}


file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_Module_visit_comparison_wo_unknown_Enrichment_KEGG_result.csv",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  Enrichment_Module_Visit_model_wo_unknown_KEGG <- SEARchways::BIGprofiler(
    module_gene_list,
    collection = "C2",
    subcollection = "CP:KEGG_LEGACY",
    ID = "SYMBOL"
  )
  write_csv(Enrichment_Module_Visit_model_wo_unknown_KEGG, file.tmp)
} else {
  Enrichment_Module_Visit_model_wo_unknown_KEGG <- read_csv(file.tmp)
}
```

### Cold A: Key Modules
"module_09", "module_14", "module_16"
```{r}
cold_A_mods <- c("module_09", "module_14", "module_16")

dat_WGCNA_Visit_significant_genes_w_meta <- dat.WGCNA_Visit_significant_genes$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  left_join(dat.WGCNA_Visit_significant_genes$targets) %>%
  pivot_longer(
    starts_with("module_"),
    names_to = "Module",
    values_to = "moduleExp"
  )

library(ggbeeswarm)
Module_09_visit_plot <- dat_WGCNA_Visit_significant_genes_w_meta %>%
  subset(Module == "module_09") %>%
  ggplot(aes(x = visit, y = moduleExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec)
ggsave(
  file.path(dirPlots, "Module_visit_significant_gene", "Visit_Module_09.pdf"),
  plot = Module_09_visit_plot,
  height = 5,
  width = 7
)

ggsave(
  file.path(
    dirPlots,
    "Module_visit_significant_gene",
    "Enrichment_Module_09.pdf"
  ),
  plot = BIGpicture::plot_enrich(
    Enrichment_Module_Visit_model_wo_unknown_H %>% subset(group == "09")
  ),
  height = 5,
  width = 7
)
#dir.create(file.path(dirPlots,"Module_visit_significant_gene"))

Module_14_visit_plot <- dat_WGCNA_Visit_significant_genes_w_meta %>%
  subset(Module == "module_14") %>%
  ggplot(aes(x = visit, y = moduleExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec)

ggsave(
  file.path(dirPlots, "Module_visit_significant_gene", "Visit_Module_14.pdf"),
  plot = Module_14_visit_plot,
  height = 5,
  width = 7
)
ggsave(
  file.path(
    dirPlots,
    "Module_visit_significant_gene",
    "Enrichment_Module_14.pdf"
  ),
  plot = BIGpicture::plot_enrich(
    Enrichment_Module_Visit_model_wo_unknown_H %>% subset(group == "14")
  ),
  height = 5,
  width = 7
)

Module_16_visit_plot <- dat_WGCNA_Visit_significant_genes_w_meta %>%
  subset(Module == "module_16") %>%
  ggplot(aes(x = visit, y = moduleExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec)

ggsave(
  file.path(dirPlots, "Module_visit_significant_gene", "Visit_Module_16.pdf"),
  plot = Module_16_visit_plot,
  height = 5,
  width = 7
)
ggsave(
  file.path(
    dirPlots,
    "Module_visit_significant_gene",
    "Enrichment_Module_16.pdf"
  ),
  plot = BIGpicture::plot_enrich(
    Enrichment_Module_Visit_model_wo_unknown_H %>% subset(group == "16")
  ),
  height = 5,
  width = 7
)
```


### Screening: Key Modules
"module_03", "module_04", "module_18"
```{r}
screening_mods <- c("module_03", "module_04", "module_18")

dat_WGCNA_Visit_significant_genes_w_meta <- dat.WGCNA_Visit_significant_genes$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  left_join(dat.WGCNA_Visit_significant_genes$targets) %>%
  pivot_longer(
    starts_with("module_"),
    names_to = "Module",
    values_to = "moduleExp"
  )

library(ggbeeswarm)
Module_03_visit_plot <- dat_WGCNA_Visit_significant_genes_w_meta %>%
  subset(Module == "module_03") %>%
  ggplot(aes(x = visit, y = moduleExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec)
ggsave(
  file.path(dirPlots, "Module_visit_significant_gene", "Visit_Module_03.pdf"),
  plot = Module_03_visit_plot,
  height = 5,
  width = 7
)
ggsave(
  file.path(
    dirPlots,
    "Module_visit_significant_gene",
    "Enrichment_Module_03.pdf"
  ),
  plot = BIGpicture::plot_enrich(
    Enrichment_Module_Visit_model_wo_unknown_H %>% subset(group == "03")
  ),
  height = 5,
  width = 7
)
BIGpicture::plot_enrich(
  Enrichment_Module_Visit_model_wo_unknown_GO %>% subset(group == "03")
)
#dir.create(file.path(dirPlots,"Module_visit_significant_gene"))

Module_04_visit_plot <- dat_WGCNA_Visit_significant_genes_w_meta %>%
  subset(Module == "module_04") %>%
  ggplot(aes(x = visit, y = moduleExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec)
ggsave(
  file.path(dirPlots, "Module_visit_significant_gene", "Visit_Module_04.pdf"),
  plot = Module_04_visit_plot,
  height = 5,
  width = 7
)
ggsave(
  file.path(
    dirPlots,
    "Module_visit_significant_gene",
    "Enrichment_Module_04.pdf"
  ),
  plot = BIGpicture::plot_enrich(
    Enrichment_Module_Visit_model_wo_unknown_H %>% subset(group == "04")
  ),
  height = 5,
  width = 7
)

Module_18_visit_plot <- dat_WGCNA_Visit_significant_genes_w_meta %>%
  subset(Module == "module_18") %>%
  ggplot(aes(x = visit, y = moduleExp, color = visit)) +
  geom_boxplot(outlier.shape = NA) +
  geom_beeswarm(alpha = 0.8, size = 3) +
  theme_classic_custom +
  scale_color_manual(values = visit_col_vec)
ggsave(
  file.path(dirPlots, "Module_visit_significant_gene", "Visit_Module_18.pdf"),
  plot = Module_18_visit_plot,
  height = 5,
  width = 7
)
ggsave(
  file.path(
    dirPlots,
    "Module_visit_significant_gene",
    "Enrichment_Module_18.pdf"
  ),
  plot = BIGpicture::plot_enrich(
    Enrichment_Module_Visit_model_wo_unknown_H %>% subset(group == "18")
  ),
  height = 5,
  width = 7
)
```


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


# Add Viral Load Metadata

```{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))


dat_voom_wo_unknown_swab_w_viral <- HoardR::cache(
  HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_swab,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_w_viral",
    ext = "RDS"
  ),
  {
    dat_voom_wo_unknown_swab_w_viral <- dat_voom_wo_unknown
    dat_voom_wo_unknown_swab_w_viral$targets <- dat_voom_wo_unknown_swab_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")
        )
      )
    dat_voom_wo_unknown_swab_w_viral
  },
  overwrite = FALSE
)
```

Have the cold B visit which fall within 15 days of cold A visit be viral positive or negative
```{r}
coldA_visit_by_virus_status <- dat_voom_wo_unknown_swab_w_viral$targets %>%
  select(subjectId, visit, Visit_start_date, Virus_status) %>%
  subset(visit == "Cold A Clinic") %>%
  select(subjectId, Visit_start_date, Virus_status) %>%
  rename("Visit_start_date_ColdA" = "Visit_start_date")

coldB_visit_by_virus_status <- dat_voom_wo_unknown_swab_w_viral$targets %>%
  select(subjectId, visit, Visit_start_date) %>%
  subset(visit == "Cold B Home") %>%
  left_join(coldA_visit_by_virus_status) %>%
  subset(
    Visit_start_date >= Visit_start_date_ColdA &
      Visit_start_date <= Visit_start_date_ColdA + days(15)
  ) %>%
  select(-c(Visit_start_date_ColdA))

dat_voom_wo_unknown_swab_w_viral$targets <- dat_voom_wo_unknown_swab_w_viral$targets %>%
  left_join(
    coldB_visit_by_virus_status %>%
      select(subjectId, visit, Visit_start_date, Virus_status) %>%
      rename(Virus_status_coldB = Virus_status),
    by = c("subjectId", "visit", "Visit_start_date")
  ) %>%
  mutate(
    Virus_status = case_when(
      visit == "Cold B Home" &
        Virus_status == "No Viral Data" &
        !is.na(Virus_status_coldB) ~ Virus_status_coldB,
      TRUE ~ Virus_status
    )
  ) %>%
  select(-Virus_status_coldB)
```

# KIMMA: Visit × Viral Status Model
```{r}
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_visit_viral_incomplete_comparison_wo_unknown_kimma.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  dat_voom_wo_unknown_swab_w_viral$targets <- dat_voom_wo_unknown_swab_w_viral$targets %>%
    mutate(
      Virus_status = ifelse(Virus_status == "No Viral Data", NA, Virus_status)
    )
  model_visit_viral_wo_unknown_kimma <- kimma::kmFit(
    dat = dat_voom_wo_unknown_swab_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_wo_unknown_kimma, file = file.tmp)
} else {
  model_visit_wo_unknown_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_visit_viral_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_wo_unknown_kimma_visit1_coldA_posvsneg <- model_visit_viral_wo_unknown_kimma$lme.contrast %>%
  subset(
    contrast_ref %in%
      c("Visit 1 Negative", "Cold A Clinic Negative", "Cold B Home Negative") &
      contrast_lvl %in%
        c("Visit 1 Positive", "Cold A Clinic Positive", "Cold B Home Positive")
  ) %>%
  left_join(
    dat_voom_wo_unknown_swab_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  ) %>%
  mutate(comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_")) %>%
  subset(
    comparison %in%
      c(
        "Cold A Clinic Positive_vs_Cold A Clinic Negative",
        "Cold B Home Positive_vs_Cold B Home Negative",
        "Visit 1 Positive_vs_Visit 1 Negative"
      )
  )

kimma::summarize_kmFit(
  model_visit_viral_wo_unknown_kimma_visit1_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 <- model_visit_viral_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_wo_unknown_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  )

Significant_genes_ColdA_posvsneg_list <- list()

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

# Volcano Plot: Cold A Positive vs Negative
```{r}
Volcano_plot_Significant_genes_ColdA_posvsneg <- volcano_funcxn(
  model_visit_viral_wo_unknown_kimma_visit1_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$`Cold A Clinic Negative<----->Cold A Clinic Positive` +
  theme_classic_custom
```

# Enrichment Analysis: Cold A Positive vs Negative
```{r}
enrichment_GO_Significant_genes_ColdA_posvsneg <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_ColdA_posvsneg_list,
  collection = "C5",
  subcollection = "GO:BP"
)

enrichment_H_Significant_genes_ColdA_posvsneg <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_ColdA_posvsneg_list,
  collection = "H"
)

BIGpicture::plot_enrich(
  enrichment_GO_Significant_genes_ColdA_posvsneg,
  fdr_cutoff = c(0.01)
)
BIGpicture::plot_enrich(enrichment_H_Significant_genes_ColdA_posvsneg)
```


# KIMMA: Visit × Viral Status Model (including Screening & Cold B)
```{r}
file.tmp <- file.path(
  dirresults,
  paste(
    projectNumber_swab,
    "_visit_viral_incomplete_comparison_wo_unknown_kimma_with_Screening_and_coldB.RDS",
    sep = ""
  )
)
if (!file.exists(file.tmp)) {
  dat_voom_wo_unknown_swab_w_viral$targets <- dat_voom_wo_unknown_swab_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_swab_wo_unknown_with_Screening_and_coldB_kimma <- kimma::kmFit(
    dat = dat_voom_wo_unknown_swab_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_swab_wo_unknown_with_Screening_and_coldB_kimma,
    file = file.tmp
  )
} else {
  model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma <- readRDS(
    file.tmp
  )
}
kimma::summarize_kmFit(
  model_visit_viral_swab_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_swab_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_swab_wo_unknown_with_Screening_and_coldB_kimma_dat <- model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma$lme.contrast %>%
  left_join(
    dat_voom_wo_unknown_swab_w_viral$genes %>% select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  ) %>%
  mutate(comparison = paste(contrast_lvl, contrast_ref, sep = "_vs_"))
```

# Upset Plot (Screening & Cold B)
```{r}
fdr_cutoff <- 0.05
model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_list <- list()

for (comp_ in unique(
  model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_dat$comparison
)) {
  model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_list[[
    comp_
  ]] <- model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
    subset(FDR <= fdr_cutoff & comparison == comp_) %>%
    pull(hgnc_symbol) %>%
    unique()
}
upset(
  fromList(
    model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_list
  ),
  text.scale = 1.5,
  order.by = "freq",
  mainbar.y.label = "Significant genes\nFDR < 0.05",
  keep.order = TRUE,
  sets = names(
    model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_list
  )
)
```

# Volcano Plots: Screening vs ColdA/ColdB
```{r}
Volcano_plot_Significant_Screening_and_coldB <- volcano_funcxn(
  model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
    subset(
      comparison %in%
        c(
          "Screening_vs_Cold A Clinic_Positive",
          "Screening_vs_Cold B Home_Positive"
        )
    ) %>%
    mutate(
      volcano_comp = ifelse(
        comparison == "Screening_vs_Cold A Clinic_Positive",
        "Cold A Clinic Positive<----->Screening",
        "Cold B Home Positive<----->Screening"
      )
    ),
  fdr = 0.05,
  n_top_gene = 10
)

Volcano_plot_Significant_Screening_and_coldB$`Cold A Clinic Positive<----->Screening` +
  theme_classic_custom
Volcano_plot_Significant_Screening_and_coldB$`Cold B Home Positive<----->Screening` +
  theme_classic_custom
```

# Enrichment Analysis: Screening & Cold B Comparisons
```{r}
Significant_genes_Screening_and_coldB_list <- list()
Significant_genes_Screening_and_coldB_list[[
  "Up in Screening compared to Cold A Clinic"
]] <- model_visit_viral_swab_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[[
  "Up in Cold A Clinic compared to Screening"
]] <- model_visit_viral_swab_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[[
  "Up in Screening compared to Cold B Home Positive"
]] <- model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
  subset(comparison %in% c("Screening_vs_Cold B Home_Positive")) %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_Screening_and_coldB_list[[
  "Up in Cold B Home Positive compared to Screening"
]] <- model_visit_viral_swab_wo_unknown_with_Screening_and_coldB_kimma_dat %>%
  subset(comparison %in% c("Screening_vs_Cold B Home_Positive")) %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()

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

enrichment_H_Significant_genes_Screening_and_coldB <- SEARchways::BIGprofiler(
  gene_list = Significant_genes_Screening_and_coldB_list,
  collection = "H"
)

BIGpicture::plot_enrich(
  enrichment_GO_Significant_Screening_and_coldB,
  fdr_cutoff = c(0.001)
)
BIGpicture::plot_enrich(enrichment_H_Significant_genes_Screening_and_coldB)
```

## Asthma vs Non-Asthma at Screening
### Subset the voom object
```{r}
dat_voom_wo_unknown_swab_w_viral_Screening <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    dir = dirDataclean,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_viral_subset_Screening",
    ext = "RDS"
  ),
  {
    libs_screening_swab <- dat_voom_wo_unknown_swab_w_viral$targets %>%
      subset(visit == "Screening") %>%
      pull(libid)
    dat_voom_wo_unknown_swab_w_viral_Screening <- dat_voom_wo_unknown_swab_w_viral[,
      libs_screening_swab
    ]
    dat_voom_wo_unknown_swab_w_viral_Screening
  },
  overwrite = FALSE
)
```

## Asthma vs Non-Asthma
```{r}
file.tmp <- file.path(
  dirresults,
  paste(projectNumber_swab, "_Astham_vs_NonAsthma_Screening.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  model_Astham_vs_NonAsthma_Screening_swab_kimma <- kimma::kmFit(
    dat = dat_voom_wo_unknown_swab_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_swab_kimma, file = file.tmp)
} else {
  model_Astham_vs_NonAsthma_Screening_swab_kimma <- readRDS(file.tmp)
}

model_Astham_vs_NonAsthma_Screening_swab_kimma <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    dir = dirresults_swab,
    analysis = "model_Astham_vs_NonAsthma_Screening_swab_kimma",
    ext = "rds"
  ),
  {
    model_Astham_vs_NonAsthma_Screening_swab_kimma <- kimma::kmFit(
      dat = dat_voom_wo_unknown_swab_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
    )
  },
  overwrite = FALSE
)

kimma::summarize_kmFit(
  model_Astham_vs_NonAsthma_Screening_swab_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_swab_kimma$lm,
  FCgroup = FALSE,
  p_cutoff = c(0.01, 0.05)
) %>%
  kable(align = "l", caption = "DEGS") %>%
  kable_styling(bootstrap_options = "striped", full_width = FALSE)
```


# EOS Category Comparison
```{r}
file.tmp <- file.path(
  dirresults,
  paste(projectNumber_swab, "_EOS_Category_Screening.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  model_EOS_cat_Screening_kimma <- kimma::kmFit(
    dat = dat_voom_wo_unknown_swab_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)
```

"Data_clean/Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_swab.csv"

# T2 High vs Low (eos1-based classification)
### Add eos1-based T2 status
```{r}
library(janitor)


dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2 <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    dir = dirDataclean,
    analysis = "dat_voom_normalized_wo_unknown_w_clinical_viral_subset_Screening_w_eos1_based_T2",
    ext = "rds"
  ),
  {
    T2_HIgh_Low_classification_swab <- read_csv(
      "Data_clean/Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_swab.csv"
    ) %>%
      clean_names()
    dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2 <- dat_voom_wo_unknown_swab_w_viral_Screening
    dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets <- dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets %>%
      left_join(
        T2_HIgh_Low_classification_swab %>%
          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_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2
  },
  overwrite = FALSE
)
```

### T2 High vs Low
```{r}
model_eos1_based_T2_High_vs_Low_Screening_swab_kimma <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    dir = dirresults_swab,
    analysis = "eos1_based_T2_High_vs_Low_Screening",
    ext = "RDS"
  ),
  {
    model_eos1_based_T2_High_vs_Low_Screening_swab_kimma <- kimma::kmFit(
      dat = dat_voom_swab_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_swab_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_swab_kimma_dat <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    dir = dirresults_swab,
    analysis = "eos1_based_T2_High_vs_Low_Screening_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_based_T2_High_vs_Low_Screening_swab_kimma$lm %>%
      left_join(
        dat_voom_swab_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(
        model,
        gene,
        hgnc_symbol,
        variable,
        comparison,
        df,
        statistic,
        estimate,
        pval,
        FDR
      )
  },
  overwrite = FALSE
)
```



### Volcano Plot: T2 High vs Low
```{r}
Volcano_plot_eos1_based_T2_High_vs_Low_Screening_swab <- volcano_funcxn(
  model_eos1_based_T2_High_vs_Low_Screening_swab_kimma_dat %>%
    subset(comparison == "T2_high_vs_T2_low") %>%
    mutate(volcano_comp = "T2 Low<----->T2 High"),
  fdr = 0.05,
  n_top_gene = 10,
  add_genes_int = c("None")
)

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

HoardR::save_figures(
  plot = Volcano_plot_eos1_based_T2_High_vs_Low_Screening_swab$`T2 Low<----->T2 High` +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_volcano_swab,
    project = projectNumber_swab,
    analysis = "Volcano_plot_eos1_based_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 6,
  width = 6,
  overwrite = TRUE
)
```

### Enrichment Analysis: T2 High vs Low
```{r}
Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab <- list()
Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab[[
  "Up in T2 High"
]] <- model_eos1_based_T2_High_vs_Low_Screening_swab_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_swab[[
  "Up in T2 Low"
]] <- model_eos1_based_T2_High_vs_Low_Screening_swab_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_Screening_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab_GO_H_muppits_list <- list()

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

    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab,
      collection = "H"
    )

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

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_Screening_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab_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_Screening_swab
)

HoardR::save_figures(
  plot = plot_H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "H_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 5,
  width = 8,
  overwrite = FALSE
)

plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab_GO_H_muppits_list$GO_BP,
  comparisons = comparisons_to_plot,
  fdr_cutoff = 0.001
) +
  theme_classic_custom
print(
  plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab
)

HoardR::save_figures(
  plot = plot_GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "GO_BP_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 20,
  width = 15,
  overwrite = FALSE
)


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

HoardR::save_figures(
  plot = plot_muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "muppits_enrichment_Significant_genes_eos1_based_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 5,
  width = 8,
  overwrite = TRUE
)
```

### String Network: T2 High vs Low
```{r}  
dir.create(
  file.path(dirPlots_swab, "String_Network"),
  showWarnings = FALSE,
  recursive = TRUE
)
for (comp_ in names(
  Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab
)) {
  pullString::get_png_network(
    ids = Significant_genes_eos1_based_T2_High_vs_Low_Screening_swab[[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_swab,
      "String_Network",
      paste(
        projectNumber_swab,
        "_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_swab_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening",
    ext = "rds"
  ),
  kimma::kmFit(
    dat = dat_voom_swab_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_swab_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_swab_kimma$lm.contrast,
  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_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_kimma_lm_contrast_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_kimma_lm_contrast_dat",
    ext = "csv"
  ),
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_kimma$lm.contrast %>%
    left_join(
      dat_voom_swab_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_swab <- volcano_funcxn(
  model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab$`T2 Low<----->T2 High` +
  theme_classic_custom +
  Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab$`Non-Asthma<----->T2 High` +
  theme_classic_custom

HoardR::save_figures(
  plot = Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab$`T2 Low<----->T2 High` +
    theme_classic_custom +
    Volcano_plot_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab$`Non-Asthma<----->T2 High` +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_volcano_swab,
    project = projectNumber_swab,
    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_swab <- list()

Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab[[
  "Up in T2 High vs T2 Low"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab[[
  "Up in T2 Low vs T2 High"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab[[
  "Up in T2 High vs Non-Asthma"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab[[
  "Up in Non-Asthma vs T2 High"
]] <- model_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_kimma_lm_contrast_dat %>%
  subset(comparison == "(T2-high)_vs_Non_Asthma") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_GO_H_muppits",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_GO_H_muppits_list <- list()
    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_GO_H_muppits_list[[
      "GO_BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab,
      collection = "H"
    )

    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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"
  )
)

plot_H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "H_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 5,
  width = 8,
  overwrite = FALSE
)

plot_muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    analysis = "muppits_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening",
    ext = "png"
  ),
  height = 5,
  width = 8,
  overwrite = FALSE
)


plot_GO_BP_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_GO_BP_enrichment_Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab +
    theme_classic_custom,
  paths = HoardR::make_path(
    dir = dirPlots_enrichment_swab,
    project = projectNumber_swab,
    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_swab, "String_Network"), showWarnings = FALSE,recursive = TRUE)
for (comp_ in names(
  Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab
)) {
  pullString::get_png_network(
    ids = Significant_genes_eos1_based_Non_Asthma_vs_T2_High_vs_Low_Screening_swab[[
      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_swab,
      "String_Network",
      paste(
        projectNumber_swab,
        "_String_Network_T2_High_vs_Low_Screening_",
        gsub(pattern = " ", replacement = "_", comp_),
        ".png",
        sep = ""
      )
    )
  )
}
```


# T2 High vs Low by Asthma Status
```{r}
model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    analysis = "eos1_based_T2_High_vs_Low_by_asthma_status_Screening",
    ext = "rds",
    dir = dirresults_swab
  ),
  {
    kimma::kmFit(
      dat = dat_voom_swab_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_swab_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_T2_High_vs_Low_by_asthma_status_Screening_swab_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 Non_Asthma)_vs_(T2-high Asthma)",
          "(T2-low Non_Asthma)_vs_(T2-low Asthma)",
          "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)",
          "(T2-high Asthma)_vs_(T2-low Asthma)"
        )
    ),
  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_by_asthma_status_Screening_swab_kimma_lm_dat <- HoardR::cache(
  HoardR::make_path(
    project = projectNumber_swab,
    dir = dirresults_swab,
    analysis = "eos1_based_T2_High_vs_Low_by_asthma_status_Screening_kimma_lm_dat",
    ext = "csv"
  ),
  {
    model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma$lm %>%
      left_join(
        dat_voom_swab_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 == "asthma_statusNon_Asthma",
            "Non_Asthma_vs_Asthma",
            variable
          )
        )
      )
  },
  overwrite = FALSE
)

model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_kimma_lm_contrast_dat <-
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma_lm_contrast_dat <- HoardR::cache(
    HoardR::make_path(
      project = projectNumber_swab,
      dir = dirresults_swab,
      analysis = "eos1_based_T2_High_vs_Low_by_asthma_status_Screening_kimma_lm_contrast_dat",
      ext = "csv"
    ),
    {
      model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma$lm.contrast %>%
        left_join(
          dat_voom_swab_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 Non_Asthma)_vs_(T2-high Asthma)",
              "(T2-low Non_Asthma)_vs_(T2-low Asthma)",
              "(T2-high Non_Asthma)_vs_(T2-low Non_Asthma)",
              "(T2-high Asthma)_vs_(T2-low Asthma)"
            )
        )
    },
    overwrite = FALSE
  )
```

## Volcano Plot: T2 High vs Low by Asthma Status
```{r}
Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lm_swab <- volcano_funcxn(
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma_lm_dat %>%
    subset(comparison %in% c("T2_high_vs_T2_low", "Non_Asthma_vs_Asthma")) %>%
    mutate(
      volcano_comp = ifelse(
        comparison == "T2_high_vs_T2_low",
        "T2 Low<----->T2 High",
        "Non-Asthma<----->Asthma"
      )
    ),
  fdr = 0.05,
  n_top_gene = 10,
  add_genes_int = c("None")
)


volcano_plot_T2_High_vs_Low_T2_Asthma_status_model_swab <- Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lm_swab$`T2 Low<----->T2 High` +
  theme_classic_custom

HoardR::save_figures(
  plot = volcano_plot_T2_High_vs_Low_T2_Asthma_status_model_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_swab,
    analysis = "Volcano_plot_T2_High_vs_Low_T2_Asthma_status_model",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)


volcano_plot_Non_Asthma_vs_Asthma_T2_Asthma_status_model_swab <- Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lm_swab$`Non-Asthma<----->Asthma` +
  theme_classic_custom

HoardR::save_figures(
  plot = volcano_plot_Non_Asthma_vs_Asthma_T2_Asthma_status_model_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_volcano_swab,
    analysis = "Volcano_plot_Non_Asthma_vs_Asthma_T2_Asthma_status_model",
    ext = "png"
  ),
  overwrite = FALSE
)
```

## Volcano Plot: T2 High vs Low by Asthma Status For the contrast results
```{r}
Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_lm_contrast_swab <- volcano_funcxn(
  model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_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",
            ifelse(
              comparison == "(T2-low Non_Asthma)_vs_(T2-low Asthma)",
              "Asthma<----->Non-Asthma in T2 Low",
              comparison
            )
          )
        )
      )
    ),
  fdr = 0.2,
  n_top_gene = 10,
  add_genes_int = c("None")
)


volcano_plot_All_contrast_T2_Asthma_status_model_swab <- Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_lm_contrast_swab$`T2 Low<----->T2 High in Asthma` +
  theme_classic_custom +
  Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_lm_contrast_swab$`T2 Low<----->T2 High in Non-Asthma` +
  theme_classic_custom +
  Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_lm_contrast_swab$`Asthma<----->Non-Asthma in T2 High` +
  theme_classic_custom +
  Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_lm_contrast_swab$`Asthma<----->Non-Asthma in T2 Low` +
  theme_classic_custom +
  plot_layout(ncol = 2, guides = "collect")

HoardR::save_figures(
  plot = volcano_plot_All_contrast_T2_Asthma_status_model_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_volcano_swab,
    analysis = "Volcano_plot_All_contrast_T2_Asthma_status_model",
    ext = "png"
  ),
  height = 9,
  width = 12,
  overwrite = TRUE
)


volcano_plot_Non_Asthma_vs_Asthma_T2_Asthma_status_model_swab <- Volcano_plot_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_lm_swab$`Non-Asthma<----->Asthma` +
  theme_classic_custom

HoardR::save_figures(
  plot = volcano_plot_Non_Asthma_vs_Asthma_T2_Asthma_status_model_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_swab,
    analysis = "Volcano_plot_Non_Asthma_vs_Asthma_T2_Asthma_status_model",
    ext = "png"
  ),
  overwrite = FALSE
)
```

## 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_swab <- list()

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

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab[[
  "Up in Non-Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma_lm_dat %>%
  subset(comparison == "Non_Asthma_vs_Asthma") %>%
  subset(estimate > 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab[[
  "Up in Asthma"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma_lm_dat %>%
  subset(comparison == "Non_Asthma_vs_Asthma") %>%
  subset(estimate < 0 & FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()

Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab[[
  "Interaction T2 High vs Low by Asthma Status"
]] <- model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_kimma_lm_dat %>%
  subset(comparison == "t2_1sdT2-high:asthma_statusNon_Asthma") %>%
  subset(FDR <= 0.2) %>%
  pull(hgnc_symbol) %>%
  unique()


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

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

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

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


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

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

enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    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_swab_GO_H_muppits_list <- list()

    enrichment_GO_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab,
      collection = "C5",
      subcollection = "GO:BP"
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_GO_H_muppits_list[[
      "GO:BP"
    ]] <- enrichment_GO_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab

    enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab,
      collection = "H"
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_GO_H_muppits_list[[
      "H"
    ]] <- enrichment_H_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab

    enrichment_muppits_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_GO_H_muppits_list[[
      "muppits"
    ]] <- enrichment_muppits_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab

    enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_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_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_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_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_based_T2_High_vs_Low_and_Asthma_Non_Asthma_Screening_swab_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_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab <- muppits_mods_expression_w_meta_all_swab %>%
  subset(
    libid %in%
      dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets$libid
  ) %>%
  select(libid, geneSet, geneExp) %>%
  left_join(
    dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2$targets
  )

colors_T2_w_Non_Asthma <- make_color_palette(
  variables = dat_voom_swab_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_swab %>%
  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_swab <- muppits_mods_expression_w_meta_wo_unknown_Screening_swab %>%
  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_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_boxplot_swab,
    analysis = "Boxplot_Non_Asthma_T2_High_Low_muppits_modules_of_interest_Screening_swab",
    ext = "png"
  ),
  height = 7,
  width = 12,
  overwrite = TRUE
)
```



# EOS1 Expression (continuous)
## Run Model: EOS1 Expression as a continuous variable
```{r}
model_eos1_exp_Screening_swab_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_Screening",
    ext = "RDS"
  ),
  kimma::kmFit(
    dat = dat_voom_swab_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_swab_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_swab_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_Screening_kimma_lm_dat",
    ext = "csv"
  ),
  model_eos1_exp_Screening_swab_kimma$lm %>%
    left_join(
      dat_voom_swab_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)),
  overwrite = FALSE
)
```

### Volcano Plot: EOS1 Expression Correlations
```{r}
Volcano_plot_eos1_exp_Screening_swab <- volcano_funcxn(
  model_eos1_exp_Screening_swab_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_swab$EOS1_exp + theme_classic_custom

HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_Screening_swab$EOS1_exp +
    theme_classic_custom +
    labs(title = "Negative Correlation <-----> Positive Correlation to EOS1"),
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_volcano_swab,
    analysis = "Volcano_plot_EOS1_expression_Screening",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```

### Enrichment Analysis: EOS1 Expression
```{r}
Significant_genes_eos1_exp_Screening_swab <- list()
Significant_genes_eos1_exp_Screening_swab[[
  "Positive Correlation to EOS1"
]] <- model_eos1_exp_Screening_swab_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate > 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_Screening_swab[[
  "Negative Correlation to EOS1"
]] <- model_eos1_exp_Screening_swab_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate < 0 & FDR <= 0.05) %>%
  pull(hgnc_symbol) %>%
  unique()


enrichment_Significant_genes_eos1_exp_Screening_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "Enrichment_GO_Hallmark_muppits_Significant_genes_eos1_exp_Screening",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_Screening_swab_GO_H_muppits_list <- list()
    enrichment_Significant_genes_eos1_exp_Screening_swab_GO_H_muppits_list[[
      "GO:BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_swab,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_exp_Screening_swab_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_swab,
      collection = "H",
    )

    enrichment_Significant_genes_eos1_exp_Screening_swab_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_swab,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_exp_Screening_swab_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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_H_Significant_genes_eos1_exp_Screening_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab_top_15_term <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_swab_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.2
) +
  theme_classic_custom

print(
  plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_swab_top_15_term
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_swab_top_15_term,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_exp_Screening_top_15_term",
    ext = "png"
  ),
  height = 9,
  width = 10,
  overwrite = FALSE
)

plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_Screening",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = TRUE
)

plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_swab +
  plot_enrichment_H_Significant_genes_eos1_exp_Screening_swab +
  plot_layout(ncol = 2, guides = "collect")

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_swab +
    plot_enrichment_H_Significant_genes_eos1_exp_Screening_swab +
    plot_layout(ncol = 2, guides = "collect"),
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_H_and_muppits_combined_Significant_genes_eos1_exp_Screening",
    ext = "png"
  ),
  height = 6,
  width = 14,
  overwrite = FALSE
)
```


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

```{r}
model_eos1_exp_interation_Asthma_status_Screening_swab_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_interation_Asthma_status_Screening_kimma",
    ext = "rds"
  ),
  {
    kimma::kmFit(
      dat = dat_voom_swab_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_swab_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_swab_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_interation_Asthma_status_Screening_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_interation_Asthma_status_Screening_swab_kimma$lm %>%
      left_join(
        dat_voom_swab_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_swab <- volcano_funcxn(
  model_eos1_exp_interation_Asthma_status_Screening_swab_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_swab$EOS1_exp +
  theme_classic_custom
HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_swab$EOS1_exp +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_volcano_swab,
    analysis = "Volcano_plot_EOS1_exp_interation_Asthma_status_Screening_swab",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_swab$`Non-Asthma vs Asthma` +
  theme_classic_custom
HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_swab$`Non-Asthma vs Asthma` +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_volcano_swab,
    analysis = "Volcano_plot_Non_Asthma_vs_Asthma_interation_Asthma_status_Screening_swab",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_swab$`gene_exp:asthma_statusNon_Asthma` +
  theme_classic_custom
HoardR::save_figures(
  plot = Volcano_plot_eos1_exp_interation_Asthma_status_Screening_Screening_swab$`gene_exp:asthma_statusNon_Asthma` +
    theme_classic_custom,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_volcano_swab,
    analysis = "Volcano_plot_EOS1_exp_Asthma_Interaction_term_interation_Asthma_status_Screening_swab",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```

### EOS1 Expression Enrichment Analysis

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

Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab[[
  "Up in Non-Asthma vs Asthma"
]] <- model_eos1_exp_interation_Asthma_status_Screening_swab_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_swab[[
  "Up in Asthma vs Non-Asthma"
]] <- model_eos1_exp_interation_Asthma_status_Screening_swab_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_swab[[
  "EOS1 Exp Interaction:Asthma"
]] <- model_eos1_exp_interation_Asthma_status_Screening_swab_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_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "Enrichment_GO_Hallmark_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_GO_H_muppits_list <- list()
    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_GO_H_muppits_list[[
      "GO:BP"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab,
      collection = "C5",
      subcollection = "GO:BP"
    )

    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_GO_H_muppits_list[[
      "H"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab,
      collection = "H",
    )

    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_H_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    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_swab$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1 (Interaction)" = Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab$`Negative Correlation to EOS1`,
    "Positive Correlation to EOS1" = Significant_genes_eos1_exp_Screening_swab$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1" = Significant_genes_eos1_exp_Screening_swab$`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_swab$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1 (Interaction)" = Significant_genes_eos1_exp_interation_Asthma_status_Screening_swab$`Negative Correlation to EOS1`,
    "Positive Correlation to EOS1" = Significant_genes_eos1_exp_Screening_swab$`Positive Correlation to EOS1`,
    "Negative Correlation to EOS1" = Significant_genes_eos1_exp_Screening_swab$`Negative Correlation to EOS1`
  ))
)
```


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

#### Asthma only Model

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

model_eos1_exp_Screening_Asthma_only_swab_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_Screening_Asthma_only_kimma",
    ext = "rds"
  ),
  {
    model_eos1_exp_Screening_Asthma_only_swab_kimma <- kimma::kmFit(
      dat = dat_voom_swab_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_swab_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_swab_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_Screening_Asthma_only_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_Screening_Asthma_only_swab_kimma$lm %>%
      left_join(
        dat_voom_swab_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_swab <- volcano_funcxn(
  model_eos1_exp_Screening_Asthma_only_swab_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_swab$EOS1_exp +
  theme_classic_custom

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

### EOS1 Expression Enrichment Analysis

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

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


enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_GO_H_muppits_list <- list()

    enrichment_GO_Significant_genes_eos1_exp_Screening_Asthma_only_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Asthma_only_swab,
      collection = "C5",
      subcollection = "GO:BP"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_GO_H_muppits_list[[
      "GO:BP"
    ]] <- enrichment_GO_Significant_genes_eos1_exp_Screening_Asthma_only_swab
    enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Asthma_only_swab,
      collection = "H"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_GO_H_muppits_list[[
      "H"
    ]] <- enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_swab

    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_GO_H_muppits_list[[
      "muppits"
    ]] <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Asthma_only_swab,
      db = muppits_mods_for_enrichment
    )
    enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_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_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_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_swab)

HoardR::save_figures(
  plot = plot_enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_H_Significant_genes_eos1_exp_Screening_Asthma_only_swab",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)

plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_GO_BP_Significant_genes_eos1_exp_Screening_Asthma_only_swab",
    ext = "png"
  ),
  height = 28,
  width = 20,
  overwrite = TRUE
)


plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Asthma_only_swab_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_swab
)

HoardR::save_figures(
  plot = plot_enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_Screening_Asthma_only_swab",
    ext = "png"
  ),
  height = 5,
  width = 7,
  overwrite = FALSE
)
```


#### Non_Asthma only Model

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


model_eos1_exp_Screening_Non_Asthma_only_swab_kimma <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_Screening_Non_Asthma_only_kimma",
    ext = "rds"
  ),
  {
    kimma::kmFit(
      dat = dat_voom_swab_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_swab_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_swab_kimma_dat <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "eos1_exp_Screening_Non_Asthma_only_kimma_dat",
    ext = "csv"
  ),
  {
    model_eos1_exp_Screening_Non_Asthma_only_swab_kimma$lm %>%
      left_join(
        dat_voom_swab_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_swab <- volcano_funcxn(
  model_eos1_exp_Screening_Non_Asthma_only_swab_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_swab$EOS1_exp +
  theme_classic_custom

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

### EOS1 Expression Enrichment Analysis

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

Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab[[
  "Positive Correlation to EOS1"
]] <- model_eos1_exp_Screening_Non_Asthma_only_swab_kimma_dat %>%
  subset(comparison == "EOS1_exp") %>%
  subset(estimate > 0 & FDR <= 0.1) %>%
  pull(hgnc_symbol) %>%
  unique()
Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab[[
  "Negative Correlation to EOS1"
]] <- model_eos1_exp_Screening_Non_Asthma_only_swab_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_swab_GO_H_muppits_list <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirresults_swab,
    project = projectNumber_swab,
    analysis = "enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_GO_H_muppits_list",
    ext = "rds"
  ),
  {
    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_GO_H_muppits_list <- list()
    enrichment_GO_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab,
      collection = "C5",
      subcollection = "GO:BP"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_GO_H_muppits_list[[
      "GO:BP"
    ]] <- enrichment_GO_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab
    enrichment_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab <- SEARchways::BIGprofiler(
      gene_list = Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab,
      collection = "H"
    )
    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_GO_H_muppits_list[[
      "H"
    ]] <- enrichment_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab

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

    enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_GO_H_muppits_list
  },
  overwrite = FALSE
)


enrichment_plot_GO_BP_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_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_swab
)

HoardR::save_figures(
  plot = enrichment_plot_GO_BP_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_GO_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab",
    ext = "png"
  ),
  height = 20,
  width = 30,
  overwrite = FALSE
)

enrichment_plot_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_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_swab
)

HoardR::save_figures(
  plot = enrichment_plot_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_H_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab",
    ext = "png"
  ),
  height = 6,
  width = 10,
  overwrite = FALSE
)


enrichment_plot_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab <- plot_enrichment_comparisons(
  df = enrichment_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab_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_swab
)

HoardR::save_figures(
  plot = enrichment_plot_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab,
  paths = HoardR::make_path(
    project = projectNumber_swab,
    dir = dirPlots_enrichment_swab,
    analysis = "Enrichment_muppits_Significant_genes_eos1_exp_Screening_Non_Asthma_only_swab",
    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_swab_kimma_lm_contrast_plot <- make_contrast_scatterplot(
#   dat = model_eos1_based_T2_High_vs_Low_by_asthma_status_Screening_swab_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_swab_kimma_lm_contrast_plot$plot

scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_eos1_continuous_model_Screening_swab_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_swab_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_swab_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_swab,
    dir = dirPlots_swab,
    analysis = "scatterplot_eos1_based_Non_Asthma_vs_Asthma_in_eos1_continuous_model_Screening_swab_kimma_lm_plot",
    ext = "png"
  ),
  height = 12,
  width = 18,
  overwrite = TRUE
)
```


# Blood Eosinophil Expression (continuous)
```{r}
file.tmp <- file.path(
  dirresults,
  paste(projectNumber_swab, "_blood_eos_cont_Screening.RDS", sep = "")
)
if (!file.exists(file.tmp)) {
  model_blood_eos_cont_Screening_swab_kimma <- kimma::kmFit(
    dat = dat_voom_swab_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
  )
  saveRDS(model_blood_eos_cont_Screening_swab_kimma, file = file.tmp)
} else {
  model_blood_eos_cont_Screening_swab_kimma <- readRDS(file.tmp)
}
kimma::summarize_kmFit(
  model_blood_eos_cont_Screening_swab_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_swab_kimma_dat <- model_blood_eos_cont_Screening_swab_kimma$lm %>%
  left_join(
    dat_voom_swab_wo_unknown_w_viral_Screening_w_eos1_based_T2$genes %>%
      select(geneName, hgnc_symbol),
    by = c("gene" = "geneName")
  ) %>%
  mutate(comparison = variable)
```

# Volcano Plot: Blood Eosinophil Associations
```{r}
Volcano_plot_blood_eos_cont_Screening_swab <- volcano_funcxn(
  model_blood_eos_cont_Screening_swab_kimma_dat %>%
    subset(comparison == "log10_EOS") %>%
    mutate(volcano_comp = "log10_EOS"),
  fdr = 0.05,
  n_top_gene = 10,
  add_genes_int = c("None")
)

Volcano_plot_blood_eos_cont_Screening_swab$log10_EOS + theme_classic_custom
```




# Module Annotation with madRich
```{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,
  paste(projectNumber_swab, "_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,
  paste(projectNumber_swab, "_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,
  paste(
    projectNumber_swab,
    "_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,
  paste(
    projectNumber_swab,
    "_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)
```


# T2 Marker Heatmaps and Plots

## Color Palettes for T2 Analyses
```{r}
color_EOS <-
  c(
    paletteer::paletteer_d("khroma::highcontrast") %>% as.character(),
    "grey"
  ) %>%
  setNames(
    dat_voom_wo_unknown$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_wo_unknown$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_wo_unknown$targets %>%
      dplyr::arrange(Allergy_Sensitization_cat) %>%
      dplyr::pull(Allergy_Sensitization_cat) %>%
      unique()
  )
```

## Boxplots and Violin Plots for T2 Marker Genes
```{r}
dir.create(file.path(dirPlots, "T2_Plots"))
dirPlotsT2 <- file.path(dirPlots, "T2_Plots")

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

library(ggbeeswarm)

dat_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)


Violin_plot_screening_3_genes <- dat_wo_unknown_w_meta %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening", "Visit 1")) %>%
  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, scales = "free")


# Violin plot split by the 3 EOS category
Violin_plot_screening_3_genes_EOS_cat <- dat_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")

#Make a density histogram
density_plot_screening_3_genes <- ggplot(
  dat_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")

boxplot_screening_3_genes_EOS_cat <- dat_wo_unknown_w_meta %>%
  subset(!is.na(EOS_cat)) %>%
  subset(
    geneName %in% c("ENSG00000133110", "ENSG00000197632", "ENSG00000016490")
  ) %>%
  subset(visit %in% c("Screening")) %>%
  ggplot(aes(x = visit, y = geneExp)) +
  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 = "red"
  ) +

  facet_wrap(~hgnc_symbol, scales = "free")

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

ggsave(
  file.path(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_boxplot_density_plot_Screening.pdf"
  ),
  plot = boxplot_screening_3_genes_EOS_cat +
    density_plot_screening_3_genes +
    plot_layout(ncol = 1, guides = "collect"),
  height = 10,
  width = 10
)


Violin_plot_screening_3_genes_EOS_cat <- dat_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(
    dirPlotsT2,
    "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_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(
    dirPlotsT2,
    "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_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(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_boxplot_Screening_by_Allergy_Sensitization_Category.pdf"
  ),
  plot = Violin_plot_screening_3_genes_Allergy_Sensitization_cat,
  height = 5,
  width = 8
)
```

```{r}
postn_serbinb2_clca1_w_clinical_screening <- dat_voom_wo_unknown$E %>%
  t() %>%
  as.data.frame() %>%
  rownames_to_column("libid") %>%
  select(libid, ENSG00000133110, ENSG00000197632, ENSG00000016490) %>%
  left_join(
    dat_voom_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_mat <- postn_serbinb2_clca1_w_clinical_screening %>%
  rownames_to_column("t") %>%
  select(-c(t)) %>%
  column_to_rownames("subjectId") %>%
  as.matrix()

postn_serbinb2_clca1_w_clinical_screening_scale <- scale(
  postn_serbinb2_clca1_w_clinical_screening_mat
)


ranks_data <- postn_serbinb2_clca1_w_clinical_screening %>%
  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 <- rownames(ranks_data)

Rank_data_meta <- postn_serbinb2_clca1_w_clinical_screening %>%
  mutate(did = subjectId) %>%
  rownames_to_column("c") %>%
  select(-c("c")) %>%
  column_to_rownames("did")
Rank_data_meta <- Rank_data_meta[rn, ]

Rank_data_meta_mat <- Rank_data_meta %>% select(-c(subjectId)) %>% as.matrix()


pdf(
  file.path(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Ranked_heatmap_Screening.pdf"
  ),
  height = 15,
  width = 15
)
ComplexHeatmap::Heatmap(
  postn_serbinb2_clca1_w_clinical_screening_scale[rn, ],
  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_mat[rn, ], i, j)
      ),
      x,
      y,
      gp = gpar(fontsize = 10)
    )
  }
)
invisible(dev.off())
```


# Correlation Heatmaps: T2 Markers & Clinical Measures
```{r}
cor_matrix_T2 <- dat_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(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_heatmap_Screening.pdf"
  ),
  height = 15,
  width = 15
)
pheatmap(
  cor_matrix_T2,
  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",
)
invisible(dev.off())


cor_matrix_T2_log10 <- dat_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(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_heatmap_Screening_log10_EOS_FENO.pdf"
  ),
  height = 8,
  width = 10
)
pheatmap(
  cor_matrix_T2_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())
```

# Regression Plots: T2 Markers vs Clinical Measures
```{r}
correlation_regression_T2 <- dat_wo_unknown_w_meta %>%
  subset(visit == "Screening") %>%
  select(
    subjectId,
    hgnc_symbol,
    geneExp,
    log10_EOS,
    log10_FENO,
    `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(
    subjectId,
    log10_EOS,
    log10_FENO,
    `Eosinophil Count (V0)`,
    `FeNO (V0)`,
    `Allergen Sensitizations per AST (V0)`,
    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 <- cor_matrix_T2 %>%
  as.data.frame() %>%
  rownames_to_column("var") %>%
  pivot_longer(-var, names_to = "var2", values_to = "correlation")
# make all combinations (pairs)
pairs_T2 <- 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 <- pairs_T2 %>%
  mutate(x_rank = match(x, vars_T2), y_rank = match(y, vars_T2)) %>%
  filter(x_rank < y_rank) %>%
  left_join(corr_long_T2, by = c("x" = "var", "y" = "var2"))
# generate plots for each pair
correlation_regression_T2_plots <- apply(pairs_T2, 1, function(row) {
  xvar <- row[["x"]]
  yvar <- row[["y"]]
  r_val <- round(as.numeric(row[["correlation"]]), 2)
  ggplot(correlation_regression_T2, 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_combined_plot <- wrap_plots(
  correlation_regression_T2_plots,
  ncol = 5
)

correlation_regression_T2_combined_plot

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


vars_T2_log10 <- c(
  "log10_EOS",
  "log10_FENO",
  "Allergy_Sensitization_per_AST",
  "CLCA1",
  "POSTN",
  "SERPINB2"
)
corr_long_T2_log10 <- cor_matrix_T2_log10 %>%
  as.data.frame() %>%
  rownames_to_column("var") %>%
  pivot_longer(-var, names_to = "var2", values_to = "correlation")
# make all combinations (pairs)
pairs_T2_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_log10 <- pairs_T2_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_log10, by = c("x" = "var", "y" = "var2"))
# generate plots for each pair
correlation_regression_T2_log10_plots <- apply(
  pairs_T2_log10,
  1,
  function(row) {
    xvar <- row[["x"]]
    yvar <- row[["y"]]
    r_val <- round(as.numeric(row[["correlation"]]), 2)
    ggplot(correlation_regression_T2, 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_log10_combined_plot <- wrap_plots(
  correlation_regression_T2_log10_plots,
  ncol = 5
)

correlation_regression_T2_log10_combined_plot

ggsave(
  file.path(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_log10_regression_Screening.pdf"
  ),
  plot = correlation_regression_T2_log10_combined_plot,
  width = 17,
  height = 10
)
```

# Regression Plots by EOS Category
```{r}
correlation_regression_T2_EOS_cat <- dat_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 <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  eso_3_gene_scatter_plot_list_eos_together[[paste(
    "EOS vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_EOS_cat,
    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_combined <- wrap_plots(
  eso_3_gene_scatter_plot_list_eos_together,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_eos_3_gene_by_cat_Screening.pdf"
  ),
  plot = eso_3_gene_scatter_plot_list_eos_together_combined,
  width = 10,
  height = 10
)


eso_3_gene_scatter_plot_list_log10eos_together <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  eso_3_gene_scatter_plot_list_log10eos_together[[paste(
    "EOS vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_EOS_cat,
    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_log10eos_together_combined <- wrap_plots(
  eso_3_gene_scatter_plot_list_log10eos_together,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlotsT2,
    "POST_CLCA1_SERPINB2_genes_w_Clinical_Correlation_log10eos_3_gene_by_cat_Screening.pdf"
  ),
  plot = eso_3_gene_scatter_plot_list_log10eos_together_combined,
  width = 10,
  height = 10
)
```


# Regression Plots by FENO Category
```{r}
correlation_regression_T2_FENO_cat <- dat_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 <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO[[paste(
    "FENO vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_FENO_cat,
    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_combined <- wrap_plots(
  FENO_3_gene_scatter_plot_list_FENO_together_log10_FENO,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlotsT2,
    "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_combined,
  width = 10,
  height = 10
)
```


# Regression Plots by Allergy Sensitization Category
```{r}
correlation_regression_T2_Allergy_Sensitization_cat <- dat_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 <- list()

for (gene_ in c("CLCA1", "POSTN", "SERPINB2")) {
  Allergy_Sensitization_3_gene_scatter_plot_list_together[[paste(
    "Allergy_Sensitization vs ",
    gene_
  )]] <- ggplot(
    correlation_regression_T2_Allergy_Sensitization_cat,
    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_combined <- wrap_plots(
  Allergy_Sensitization_3_gene_scatter_plot_list_together,
  ncol = 1,
  guides = "collect"
)
ggsave(
  file.path(
    dirPlotsT2,
    "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_combined,
  width = 10,
  height = 10
)
```


# MUPPITS Modules: Correlations with EOS Expression
## 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_swab <- 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_swab_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_swab_w_meta <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_swab,
    analysis = "Muppits_expression_data",
    ext = "csv"
  ),
  {
    Muppits_modules_mean_expression_swab %>%
      pivot_longer(!geneSet, names_to = "libid", values_to = "geneExp") %>%
      left_join(dat_voom_swab_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_swab_w_meta <- Muppits_modules_mean_expression_swab_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_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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 swab")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  theme_classic_custom

ggplot(
  Muppits_modules_mean_expression_swab_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 swab")
  ) +
  stat_cor(method = "pearson") +
  theme_classic_custom


plot_unma11_vs_EOS_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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 swab")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  theme_classic_custom


plot_eos1_vs_EOS_swab +
  plot_unma11_vs_EOS_swab +
  plot_layout(ncol = 1, guides = "collect")


HoardR::save_figures(
  plot = plot_eos1_vs_EOS_swab +
    plot_unma11_vs_EOS_swab +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Contrast_Scatter,
    project = projectNumber_swab,
    analysis = "EOS_vs_Muppits_EOS1_and_unma11",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)


plot_eos1_vs_EOS_swab_asthma_status <- ggplot(
  Muppits_modules_mean_expression_swab_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 swab")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  facet_grid(~asthma_status) +
  theme_classic_custom


plot_unma11_vs_EOS_swab_asthma_status <- ggplot(
  Muppits_modules_mean_expression_swab_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 swab")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_EOS) +
  facet_grid(~asthma_status) +
  theme_classic_custom


plot_eos1_vs_EOS_swab_asthma_status +
  plot_unma11_vs_EOS_swab_asthma_status +
  plot_layout(ncol = 1, guides = "collect")


HoardR::save_figures(
  plot = plot_eos1_vs_EOS_swab_asthma_status +
    plot_unma11_vs_EOS_swab_asthma_status +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Contrast_Scatter,
    project = projectNumber_swab,
    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_swab <- Muppits_modules_mean_expression_swab_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 swab") +
  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_swab <- Muppits_modules_mean_expression_swab_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 swab"
  ) +
  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_swab +
  boxplot_plot_screening_muppits_unma11_EOS_cat_swab +
  plot_layout(ncol = 2, guides = "collect")


HoardR::save_figures(
  plot = boxplot_plot_screening_muppits_eos1_EOS_cat_swab +
    boxplot_plot_screening_muppits_unma11_EOS_cat_swab +
    plot_layout(ncol = 2, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Boxplot,
    project = projectNumber_swab,
    analysis = "EOS_cat_Muppits_eos1_and_unma11",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)


boxplot_plot_screening_muppits_eos1_EOS_cat_swab_asthma_status <- Muppits_modules_mean_expression_swab_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 swab") +
  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_swab_asthma_status <- Muppits_modules_mean_expression_swab_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 swab"
  ) +
  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_swab_asthma_status +
  boxplot_plot_screening_muppits_unma11_EOS_cat_swab_asthma_status +
  plot_layout(ncol = 2, guides = "collect")


HoardR::save_figures(
  plot = boxplot_plot_screening_muppits_eos1_EOS_cat_swab_asthma_status +
    boxplot_plot_screening_muppits_unma11_EOS_cat_swab_asthma_status +
    plot_layout(ncol = 2, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Boxplot,
    project = projectNumber_swab,
    analysis = "EOS_cat_Muppits_eos1_and_unma11_by_asthma_status",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 13,
  height = 7
)


boxplot_plot_screening_muppits_eos1_swab <- Muppits_modules_mean_expression_swab_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 swab") +
  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_swab,
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Boxplot,
    project = projectNumber_swab,
    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_swab <- ggplot(
  Muppits_modules_mean_expression_swab_w_meta %>%
    subset(!is.na(EOS_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", 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_swab <- ggplot(
  Muppits_modules_mean_expression_swab_w_meta %>%
    subset(!is.na(EOS_cat)) %>%

    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_swab +
    density_plot_screening_muppits_eos1_unma11_eos_cat_swab +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Density,
    project = projectNumber_swab,
    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_swab_asthma_status <- ggplot(
  Muppits_modules_mean_expression_swab_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_swab_asthma_status <- ggplot(
  Muppits_modules_mean_expression_swab_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_swab_asthma_status +
    density_plot_screening_muppits_eos1_unma11_eos_cat_swab_asthma_status +
    plot_layout(ncol = 1, guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Density,
    project = projectNumber_swab,
    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_swab <- ggplot(
  Muppits_modules_mean_expression_swab_w_meta %>%
    subset(!is.na(EOS_cat)) %>%

    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 swab") +
  theme(legend.position = "none")


density_plot_screening_muppits_eos1_eos_cat_swab /
  (boxplot_plot_screening_muppits_eos1_EOS_cat_swab +
    xlab("EOS1 Expression") +
    theme(legend.position = "none") +
    plot_eos1_vs_EOS_swab +
    xlab("Blood eos counts (log10)")) +
  plot_layout(guides = "collect")

HoardR::save_figures(
  plot = density_plot_screening_muppits_eos1_eos_cat_swab /
    (boxplot_plot_screening_muppits_eos1_EOS_cat_swab +
      xlab("EOS1 Expression") +
      theme(legend.position = "none") +
      plot_eos1_vs_EOS_swab +
      xlab("Blood eos counts (log10)")) +
    plot_layout(guides = "collect"),
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots,
    project = projectNumber_swab,
    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_swab_w_meta_eos1_screening <- Muppits_modules_mean_expression_swab_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_swab <- Muppits_modules_mean_expression_swab_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_swab_w_meta_eos1_screening_classified <- Muppits_modules_mean_expression_swab_w_meta_eos1_screening %>%
  #left_join(stats_per_group_eos1_swab, by = "EOS_cat") %>%
  mutate(
    T2_mean = ifelse(
      geneExp <= stats_per_group_eos1_swab$mean_val,
      "T2-low",
      "T2-high"
    ),
    T2_1SD = ifelse(
      geneExp <= stats_per_group_eos1_swab$mean_plus_1SD,
      "T2-low",
      "T2-high"
    ),
    T2_2SD = ifelse(
      geneExp <= stats_per_group_eos1_swab$mean_plus_2SD,
      "T2-low",
      "T2-high"
    ),
    T2_1SE = ifelse(
      geneExp <= stats_per_group_eos1_swab$mean_plus_1SE,
      "T2-low",
      "T2-high"
    ),
    T2_2SE = ifelse(
      geneExp <= stats_per_group_eos1_swab$mean_plus_2SE,
      "T2-low",
      "T2-high"
    ),
  )
Muppits_modules_mean_expression_swab_w_meta_eos1_screening_classified <- HoardR::cache(
  paths = HoardR::make_path(
    dir = dirDataclean,
    project = projectNumber_swab,
    analysis = "Donor_classification_T2_high_low_based_on_eos1_exp_in_EOS_less_than_150_nasal_swab",
    ext = "csv"
  ),
  {
    Muppits_modules_mean_expression_swab_w_meta_eos1_screening %>%
      #left_join(stats_per_group_eos1_swab, by = "EOS_cat") %>%
      mutate(
        T2_mean = ifelse(
          geneExp <= stats_per_group_eos1_swab$mean_val,
          "T2-low",
          "T2-high"
        ),
        T2_1SD = ifelse(
          geneExp <= stats_per_group_eos1_swab$mean_plus_1SD,
          "T2-low",
          "T2-high"
        ),
        T2_2SD = ifelse(
          geneExp <= stats_per_group_eos1_swab$mean_plus_2SD,
          "T2-low",
          "T2-high"
        ),
        T2_1SE = ifelse(
          geneExp <= stats_per_group_eos1_swab$mean_plus_1SE,
          "T2-low",
          "T2-high"
        ),
        T2_2SE = ifelse(
          geneExp <= stats_per_group_eos1_swab$mean_plus_2SE,
          "T2-low",
          "T2-high"
        ),
      )
  },
  overwrite = FALSE
)


# write_csv(
#   Muppits_modules_mean_expression_swab_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_swab.csv"
# )

# Count samples per EOS group for each cutpoint
Muppits_modules_mean_expression_swab_w_meta_eos1_screening_classified_summary_table <- Muppits_modules_mean_expression_swab_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_swab <- stats_per_group_eos1_swab %>%
  select(mean_val, mean_plus_1SD, mean_plus_2SD, mean_plus_1SE, mean_plus_2SE)


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

# Density plot with vertical lines for each EOS group
density_plot_screening_muppits_eos1_unma11_eos_cat_swab_split_T2_high_low <- ggplot(
  Muppits_modules_mean_expression_swab_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_swab,
    aes(xintercept = mean_val),
    color = "black",
    linetype = "solid",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_swab,
    aes(xintercept = mean_plus_1SD),
    color = "blue",
    linetype = "dashed",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_swab,
    aes(xintercept = mean_plus_2SD),
    color = "red",
    linetype = "dotted",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_swab,
    aes(xintercept = mean_plus_1SE),
    color = "darkgreen",
    linetype = "longdash",
    linewidth = 1
  ) +
  geom_vline(
    data = stats_per_group_eos1_swab,
    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_swab_split_T2_high_low,
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Density,
    project = projectNumber_swab,
    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_swab_split_T2_high_low_only_mean_1sd <- ggplot(
  Muppits_modules_mean_expression_swab_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_swab,
    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_swab_split_T2_high_low_only_mean_1sd,
  paths = HoardR::make_path(
    dir = dirT2plots_swab$plots_Density,
    project = projectNumber_swab,
    analysis = "Density_plot_eos1_expression_by_EOS_cat_only_mean_1sd_cutpoint",
    ext = "pdf"
  ),
  overwrite = FALSE,
  width = 10,
  height = 7
)
```

# plot log10FENO vs eos1 mod and unma11

```{r}
plot_eos1_vs_FENO_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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 = "eos1 Expression",
    title = paste0("log10_FENO", " vs ", "eos1", " Nasal Swab")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_FENO) +
  theme_classic_custom


plot_unma11_vs_FENO_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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 Swab")
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_FENO) +
  theme_classic_custom


plot_eos1_vs_FENO_swab +
  plot_unma11_vs_FENO_swab +
  plot_layout(ncol = 1, guides = "collect")


ggsave(
  file.path(dirPlotsT2, "FENO_vs_Muppits_eos1_and_unma11.pdf"),
  plot = plot_eos1_vs_FENO_swab +
    plot_unma11_vs_FENO_swab +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)


boxplot_plot_screening_muppits_eos1_FENO_cat_swab <- Muppits_modules_mean_expression_swab_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 Swab") +
  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_swab <- Muppits_modules_mean_expression_swab_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 Swab"
  ) +
  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_swab +
  boxplot_plot_screening_muppits_unma11_FENO_cat_swab +
  plot_layout(ncol = 2, guides = "collect")


ggsave(
  file.path(dirPlotsT2, "box_plot_FENO_cat_Muppits_eos1_and_unma11.pdf"),
  plot = boxplot_plot_screening_muppits_eos1_FENO_cat_swab +
    boxplot_plot_screening_muppits_unma11_FENO_cat_swab +
    plot_layout(ncol = 2, guides = "collect"),
  width = 10,
  height = 7
)

#Make the density plot just eos1 expression coloered by category
#Make a density histogram
density_plot_screening_muppits_eos1_unma11_FENO_cat_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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(
    dirPlotsT2,
    "density_plot_screening_muppits_eos1_unma11_FENO_cat.pdf"
  ),
  plot = density_plot_screening_muppits_eos1_unma11 +
    density_plot_screening_muppits_eos1_unma11_FENO_cat_swab +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)
```


# plot log10Allergy_Sensitization vs eos1 mod and unma11

```{r}
plot_eos1_vs_Allergy_Sensitization_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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 Swab"
    )
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_Allergy_sens) +
  theme_classic_custom


plot_unma11_vs_Allergy_Sensitization_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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 Swab"
    )
  ) +
  stat_cor(method = "pearson") +
  scale_color_manual(values = color_Allergy_sens) +
  theme_classic_custom


plot_eos1_vs_Allergy_Sensitization_swab +
  plot_unma11_vs_Allergy_Sensitization_swab +
  plot_layout(ncol = 1, guides = "collect")


ggsave(
  file.path(dirPlotsT2, "Allergy_Sensitization_vs_Muppits_eos1_and_unma11.pdf"),
  plot = plot_eos1_vs_Allergy_Sensitization_swab +
    plot_unma11_vs_Allergy_Sensitization_swab +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)


boxplot_plot_screening_muppits_eos1_Allergy_Sensitization_cat_swab <- Muppits_modules_mean_expression_swab_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 Swab"
  ) +
  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_swab <- Muppits_modules_mean_expression_swab_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 Swab"
  ) +
  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_swab +
  boxplot_plot_screening_muppits_unma11_Allergy_Sensitization_cat_swab +
  plot_layout(ncol = 2, guides = "collect")


ggsave(
  file.path(
    dirPlotsT2,
    "box_plot_Allergy_Sensitization_cat_Muppits_eos1_and_unma11.pdf"
  ),
  plot = boxplot_plot_screening_muppits_eos1_Allergy_Sensitization_cat_swab +
    boxplot_plot_screening_muppits_unma11_Allergy_Sensitization_cat_swab +
    plot_layout(ncol = 2, guides = "collect"),
  width = 10,
  height = 7
)


density_plot_screening_muppits_eos1_unma11_Allergy_Sensitization_cat_swab <- ggplot(
  Muppits_modules_mean_expression_swab_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(
    dirPlotsT2,
    "density_plot_screening_muppits_eos1_unma11_Allergy_Sensitization_cat.pdf"
  ),
  plot = density_plot_screening_muppits_eos1_unma11 +
    density_plot_screening_muppits_eos1_unma11_Allergy_Sensitization_cat_swab +
    plot_layout(ncol = 1, guides = "collect"),
  width = 10,
  height = 7
)
```


# Make a correlation heatmap with lower triange shwoing the correlation and upper triangle showing scatterplot of correlation
```{r}
df_num_swab <- Muppits_modules_mean_expression_swab_w_meta %>%
  subset(geneSet %in% c("eos1", "unma11") & visit == "Screening") %>%
  select(
    subjectId,
    visit,
    geneSet,
    geneExp,
    log10_EOS,
    log10_FENO,
    `Allergen Sensitizations per AST (V0)`
  ) %>%
  distinct() %>%
  pivot_wider(names_from = geneSet, values_from = geneExp) %>%
  select(-c(visit, subjectId)) %>%
  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_swab,
  lower = list(continuous = my_cor_tile),
  upper = list(continuous = my_scatter_lm),
  diag = list(continuous = wrap("densityDiag")),
  axisLabels = "show"
)


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


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

P640_1_muppits_module_coherence <-
  calculate_module_coherence(
    mods = muppits_mods %>%
      rename("geneName" = "ensemblID", "module.char" = "geneSet"),
    dat = dat_voom_wo_unknown,
    mods_title = "MUPPITS Modules",
    dat_title = "P640-1 Swab",
    r_cutoff = 0.3,
    p_cutoff = 0.05
  )
```