E_Proteomics_PathwayAnalysisResults.Rmd

---
title: "Visualization and Interpretation of Pathway Analysis Results"
output: html_document
date: "2025-07-02"
---

# Load packages

```{r message = FALSE, warning = FALSE}
# Clears global environment
rm(list = ls(all.names = TRUE)) 

# If needed, install and load packages
library(this.path) # For file path
library(openxlsx) # For data import
library(readxl) # For data import
library(tidyverse) # For data organization
library(ggrepel) # for ggplot text labels
library(cowplot) # for plot themes
library(janitor) # for data cleaning
library(scales) # for graphing
library(rstatix) # for correlations

# Redefine function that is masked
select <- dplyr::select

# Set working directory
setwd(this.dir())
```

# Import data

Import IPA results as well as proteomics, demographics, and sputum data. For the demographic and sputum data, filter to only participants with proteomics data. For the sputum data, remove 6 hours post exposure. Then, keep only participants with sputum data in the proteomics and demographic dataframes. 
```{r}
# IPA results
ipa_all <- read_excel("1_InputData/IPA_Pathways_Aggregate_10.01.2024.xls")
ipa_fem <- read_excel("1_InputData/IPA_Pathways_Females_10.01.2024.xls")
ipa_male <- read_excel("1_InputData/IPA_Pathways_Males_10.01.2024.xls")

# Proteomics data
prot <- read.xlsx("2_ProcessedData/ProcessedData_Proteomics_Norm_Filt_Log2_07.29.2025.xlsx")

# Demographics data
demo <- read.xlsx("2_ProcessedData/ProcessedData_Demographics_07.29.2025.xlsx") %>%
  filter(Proteomics_Run == "Yes") %>%
  clean_names()

# Sputum data
sputum <- read.xlsx("2_ProcessedData/ProcessedData_SputumCellsAndCytokines_07.29.2025.xlsx") %>% 
  filter(Subject_ID_Original %in% demo$subject_id_original) %>% 
  clean_names() %>%
  filter(!visit == "Post_6hrs")

# Filter demographics data to only the participants who have sputum data
demo_filtered <- demo %>% 
  filter(subject_id_original %in% sputum$subject_id_original)
```

# Write out IPA results as supplemental tables

```{r}
write.xlsx(ipa_all, "3_OutputTables/SupplementalTableS6_Proteomics_PathwayAnalysisAggregate.xlsx")
write.xlsx(ipa_male, "3_OutputTables/SupplementalTableS7_Proteomics_PathwayAnalysisMale.xlsx")
write.xlsx(ipa_fem, "3_OutputTables/SupplementalTableS8_Proteomics_PathwayAnalysisFemale.xlsx")
```


# Create bar graph

Clean the data to prepare for graphing: 
```{r}
ipa_all <- clean_names(ipa_all)
ipa_fem <- clean_names(ipa_fem)
ipa_male <- clean_names(ipa_male)

# I will represent this graph using log2 p-values, therefore I will multiply the Neg_Log_P.value by -1.
ipa_all$log_pval <- -1 * ipa_all$log_p_value
ipa_fem$log_pval <- -1 * ipa_fem$log_p_value
ipa_male$log_pval <- -1 * ipa_male$log_p_value

colnames(ipa_all)[which(names(ipa_all) == "log_p_value")] <- "neg_log_p_value"
colnames(ipa_fem)[which(names(ipa_fem) == "log_p_value")] <- "neg_log_p_value"
colnames(ipa_male)[which(names(ipa_male) == "log_p_value")] <- "neg_log_p_value"

# Selecting for biologically relevant pathways.
all_10 <- ipa_all[c(1, 5, 8, 20, 57, 64, 78, 85, 191),]
fem_10 <- ipa_fem[c(1, 5, 8, 20, 57, 64, 78, 85, 191),]
male_10 <- ipa_male[c(1, 5, 8, 20, 57, 64, 78, 85, 191),]

# Add a col identifier before combining dfs
all_10$group <- 'All'
fem_10$group <- 'F'
male_10$group <- 'M'
 
data <- bind_rows(all_10, fem_10)
data <- bind_rows(data, male_10)
```

```{r}
# Create plot.
comp_plot <- data %>% 
  ggplot(aes(x = reorder(ingenuity_canonical_pathways, -log_pval), y = z_score, fill = group, color = group)) +
  geom_col(position = position_dodge(0.8), linewidth = 1.1, width = 0.65) +
  
  scale_fill_manual(labels = c("All", "Females", "Males"), values = c("#a475ff",  "#c261c2", "#0086f9")) +
  scale_color_manual(values = c("#3c00b0", "#702b70", "#005197"), guide = "none") +
  
  guides(fill = guide_legend(reverse=TRUE)) +
  
  scale_x_discrete(position = "top", labels = label_wrap(width = 40)) +
  scale_y_continuous(expand = c(0, 0), limits = c(-12, 5)) +
  
  geom_hline(yintercept = 0) +
  geom_hline(yintercept = 2, linetype = 2, color = "gray25") +
  geom_hline(yintercept = -2, linetype = 2, color = "gray25") +
  
  labs(title = "Canonical Pathways", x = "", y = "Z-Score", fill = "Group") +
  
  theme(axis.line = element_line(colour = "black"),
        axis.ticks = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border = element_blank(),
        plot.margin = unit(c(1, 1, 1, 1), "cm"),
        
        plot.title = element_text(hjust = 0.5, face = "bold", size = 27),
        axis.title.y = element_blank(), 
        axis.text.y = element_text(colour = "black", size = 22),
        axis.title.x = element_text(colour = "black", size = 25, vjust = -0.5), 
        axis.text.x = element_text(colour = "black", size = 25),
        
        legend.position = c(0.1, 0.5),
        legend.title = element_text(size = 25),
        legend.text = element_text(size = 25),
        legend.key.size = unit(8, 'mm'),
        legend.box.background = element_rect(color = "gray10", linewidth = 0.75)) +
  
  coord_flip() 

comp_plot

png("4_OutputFigures/PathwaysBarplot.png", width = 18, height = 10, units = "in", res = 1200)
comp_plot
invisible(dev.off())
```