---
title: "Stress type–specific small extracellular vesicle signatures_R"
output: html_document
date: "2025-07-05"
editor_options: 
  chunk_output_type: console
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

## preparation
```{r preparation}
# set working directory
setwd("~")

# load library
library(readxl)
library(writexl)
library(haven)
library(tidyverse)
library(rstatix)
library(psych)
library(afex)
library(emmeans)
library(ggplot2)
library(ggpubr)
library(ggtext)
library(smplot2)
library(randomForest)
library(datasets)
library(caret)
library(hrbrthemes)
library(ROCR)
library(rpart)
library(rpart.plot)
library(reshape2)
library(yardstick)

# color design
library(paletteer)
library(ggthemes)
library(pals)

# import EV data
data_prep <- read_excel("EVs_data.xlsx")

# scientific function
scientific_10 = function(x) {
  ifelse(
    x==0, "0",
    parse(text = sub("e\\+*", " %*% 10^", scales::scientific_format()(x)))
  )
} 

```

# prepare data
```{r prepare data}

# create a data frame
data_prep <- data.frame(data_prep)

# remove unnecessary variables
data_prep$reference <- NULL
data_prep$HMGB1 <- NULL # HMGB1

# drop low standard variables
data_prep$C3b.C3i <- NULL
data_prep$CD273 <- NULL

# note missing values as NA
data_prep[data_prep == -99] <- NA

# correct the order of factor levels for the facet grid
data_prep$stress = factor(data_prep$stress, levels = c("TSST", "Sport"))

# change format to long
data_long <- pivot_longer(data_prep, cols = "CD41":"CD206", names_to = "marker", values_to = "values")

# correct the order of all markers
marker_order <- c("CD9", "CD13", "CD14", "CD16", "CD24", "CD41", "CD44", "CD61", "CD63", "CD66b", "CD81", "CD82", "CD90", "CD100", "CD171", "CD206", "CD227", "CD235a", "CTLA.4", "HLA.ABC", "HLA.DR", "PD.L1", "PS")

data_long$marker = factor(data_long$marker, levels = marker_order)

# remove zero-values
data_long <- data_long %>%
  filter(values != 0)

```

# outliers
```{r outliers}

# check normal distribution for each marker
norm_distr <- data_long %>%
  group_by(marker) %>%
  shapiro_test(values)

norm_distr %>%
  filter(p <= .05) # all 23 markers deviate from a normal distribution

# check normal distribution for each marker, stress and time point combination
norm_distr_comb <- data_long %>%
  group_by(stress, time_point, marker) %>%
  shapiro_test(values)

norm_distr_comb %>%
  filter(p <= .05) # 195 combinations deviate from a normal distribution

# calculate the median
data_median <- data_prep %>%
  group_by(stress, time_point) %>%
  summarize(across(marker_order, ~ median(.x, na.rm = TRUE)))

# calculate the median absolute deviation (MAD) for non-normally distributed values
data_mad <- data_prep %>%
  group_by(stress, time_point) %>%
  summarize(across(marker_order, ~ mad(.x, constant = 1, na.rm = TRUE)))

# change format of the median and MAD data frame to long
data_median_long <- pivot_longer(data_median, cols = "CD9":"PS", names_to = "marker", values_to = "median")
data_mad_long <- pivot_longer(data_mad, cols = "CD9":"PS", names_to = "marker", values_to = "mad")

# add median and MAD info to the data_long dataframe
data_median_long <- data_median_long %>%
  mutate(merger = paste0(stress, time_point, marker)) # merger variable for the median data frame

data_mad_long <- data_mad_long %>%
  mutate(merger = paste0(stress, time_point, marker)) # merger variable for the mad data frame

data_long <- data_long %>%
  mutate(merger = paste0(stress, time_point, marker)) # merger variable for the data_long data frame

data_long <- merge(data_long, data_median_long %>% select(merger, median), by = "merger") # add median info
data_long <- merge(data_long, data_mad_long %>% select(merger, mad), by = "merger") # add mad info

# remove unnecessary variables
data_long$merger <- NULL
data_long$stress.x <- NULL
data_long$stress.y <- NULL

# order of the variables
data_long <- data_long %>%
  relocate(stress, .after = time_point)

# create a M ± 3*MAD variable
data_long <- data_long %>%
  mutate(median_mad_plus_3 = median + 3 * mad)

data_long <- data_long %>%
  mutate(median_mad_minus_3 = median - 3 * mad)

# remove values ± 3 MAD
data_long <- data_long %>%
  filter(values >= median_mad_minus_3 & values <= median_mad_plus_3)

# remove merger variable
data_long$merger <- NULL

```

# normal distribution
```{r normal distribution}

# check the assumption of normality again
norm_distr_mad <- data_long %>%
  group_by(stress, time_point, marker) %>%
  shapiro_test(values)

count <- norm_distr_mad %>%
  filter(p <= .05) # 146 combinations deviate  from a normal distribution

table(count$marker) # remove when proceeding; CD44 shows 0 entries and therefore normal distribution

# visualize distribution of CD44
ggdensity(data_long %>%
         select(time_point, stress, marker, values) %>%
         filter(marker == "CD44"),
         "values", ggtheme = theme_minimal()) +
  facet_grid(stress ~ time_point, labeller = "label_both") +
  labs(title = "CD44 Density Plots by Time Point and Condition",
       x = "Values",
       y = "Density")

# transform non-normally distributed variables using ln
data_long <- data_long %>%
  mutate(logvalues = values)

data_long_log <- data_long %>%
  filter(!marker %in% c("CD9", "CD44", "CD171",
                        "CD61", "CD81", "CD100")) # filter all markers that need to be ln-transformed

data_long_log <- data_long_log %>%
  mutate(logvalues = log(logvalues)) # ln-transformation

data_long_other <- data_long %>%
  filter(marker %in% c("CD61", "CD81", "CD100")) # filter all markers that need another transformation than ln

data_long_other <- data_long_other %>%
  mutate(logvalues = sqrt(logvalues)) # sqrt-transformation

data_long <- data_long %>%
  filter(marker %in% c("CD9", "CD44", "CD171")) # filter all markers that do not need to be ln-transformed

data_long <- rbind(data_long_log, data_long)
data_long <- rbind(data_long_other, data_long)

# order of the variables
data_long <- data_long %>%
  relocate(logvalues, .after = values)

# check the assumption of normality again
norm_distr_log <- data_long %>%
  group_by(stress, time_point, marker) %>%
  shapiro_test(logvalues)

count_log <- norm_distr_log %>%
  filter(p <= .05) # 33 combinations deviate from a normal distribution

table(count_log$marker)

```

# incomplete pairs
```{r incomplete pairs}

# identify incomplete pairs
comb_complete <- expand.grid(
  VPN = 1:20,             # 20 participants
  stress = c("TSST", "Sport"),            # 2 conditions
  time_point = 2:6,               # 5 time points
  marker = marker_order  # 23 markers
) # full data sample template

comb_incomplete <- comb_complete %>%
  anti_join(data_long, by = c("VPN", "stress", "time_point", "marker")) # identify missing values

comb_incomplete <- comb_incomplete %>%
  mutate(idf = paste0(VPN, stress, marker)) # identifier variable

data_long <- data_long %>%
  mutate(idf = paste0(VPN, stress, marker)) # identifier variable

idf_var <- comb_incomplete$idf

# remove incomplete pairs
data_long_complete <- data_long %>%
  filter(!idf %in% idf_var)

# baseline incomplete
comp_baseline_incomplete <- expand.grid(
  VPN = 1:20,             # 20 participants
  stress = c("TSST", "Sport"),            # 2 conditions
  time_point = 2,               # baseline time point
  marker = marker_order  # 23 markers
) # full data sample template

data_baseline <- data_long %>%
  filter(time_point == 2) # baseline data

baseline_incomplete <- comp_baseline_incomplete %>%
  anti_join(data_long, by = c("VPN", "stress", "time_point", "marker")) # identify missing values

baseline_incomplete <- baseline_incomplete %>%
  mutate(idf = paste0(VPN, marker, time_point)) # identifier variable

data_baseline <- data_baseline %>%
  mutate(idf = paste0(VPN, marker, time_point)) # identifier variable


idf_var <- baseline_incomplete$idf

# remove incomplete pairs
data_long_baseline <- data_baseline %>%
  filter(!idf %in% idf_var)

# check results
data_long_baseline %>%
  group_by(marker, stress) %>%
  summarise(n_samples = sum(!is.na(values)), .groups = "drop")

```


# time point -2 min
```{r time point -2 min}

# statistical analysis
comp_baseline <- data_long_baseline %>%
  group_by(marker) %>%
  pairwise_t_test(
    logvalues ~ stress, paired = TRUE,
    comparisons = list(c("Sport", "TSST")),
    p.adjust.method = "bonferroni"
  ) # or pairwise t-test

# add y. position
y_positions_baseline <- data_long_baseline %>%
  group_by(marker, stress) %>%
  summarise(mean_value = mean(values, na.rm = TRUE), .groups = "drop") %>%
  mutate(y.position = mean_value * 1.5)  # Adjust for spacing

# Merge the y-position into the pairwise comparison results
comp_baseline <- comp_baseline %>%
  left_join(y_positions_baseline, by = c("marker"))

# filter one stress type
comp_baseline <- comp_baseline %>%
  filter(stress == "Sport")

# compute effect sizes
effect_sizes_baseline <- data_long_baseline %>%
  group_by(marker) %>%
  cohens_d(logvalues ~ stress, paired = TRUE, 
           comparisons = list(c("Sport", "TSST")),
           ci = TRUE)

# merge effect sizes to existing data frame
comp_baseline <- comp_baseline %>%
  left_join(effect_sizes_baseline %>%
              select(marker, group1, group2, effsize, magnitude, conf.low, conf.high),
            by = c("marker", "group1", "group2"))

# relocate variables
comp_baseline <- relocate(comp_baseline, effsize, .after = p.adj.signif)
comp_baseline <- relocate(comp_baseline, magnitude, .after = effsize)
comp_baseline <- relocate(comp_baseline, conf.low, .after = magnitude)
comp_baseline <- relocate(comp_baseline, conf.high, .after = conf.low)

# plot
ggplot(mapping = aes(x = as.factor(stress), y = values, fill = stress)) +
  sm_raincloud(data = data_long_complete %>% filter(time_point == 2), position = position_nudge(x = +0.15),
               point.params = list(size = 3, shape = 21,
                                   color = 'transparent', alpha = 0.4,
                                   position = sdamr::position_jitternudge(nudge.x = 0.06, 
                                                                   jitter.width = 0.06)),
               boxplot.params = list(outlier.shape = NA)) +
  facet_wrap(~marker, scales = "free") +
  scale_fill_manual(values=c("#679289", "#0C2342"),
                    labels = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress"))) +
  scale_y_continuous(label = scientific_10) +
  stat_pvalue_manual(comp_baseline, label = "p.adj.signif", 
                     tip.length = 0, 
                     bracket.size = 0.75, bracket.shorten = 0.4 ,size = 8, 
                     hide.ns = TRUE, color = "black") +
  labs(x = "", y = "objects/mL") +
  theme_minimal() +
  theme(legend.position = "bottom", legend.title = element_blank(),
        text = element_text(size = 18),
        axis.text.x=element_blank(),
        panel.spacing = unit(2, "lines")) +
  theme(
      strip.text.x = element_text(
        size = 18, face = "bold"
        )
      )
ggsave("sEV_baseline.tiff", width = 5*3, height = 4*3, device='tiff', dpi=300)

```

# rmANOVA
```{r rmANOVA}

# test for significant differences over time
rmanova_EVs <- data_long_complete %>%
  group_by(marker) %>%
  anova_test(dv = logvalues, wid = VPN, within = c(stress, time_point), effect.size = "pes")
rmanova_EVs_output <- get_anova_table(rmanova_EVs)

# filter for significant interaction
rmanova_EVs_interac <- rmanova_EVs_output %>%
  filter(Effect == "stress:time_point" & p < .05)

# significant sEV names
EVs_sign <- as.character(rmanova_EVs_interac$marker)

# filter dataframe for significant sEVs
data_sign <- data_long_complete %>%
  filter(marker %in% EVs_sign)

```

# pairwise comparisons
```{r pairwise comparisons}

# calculate pairwise comparisons
comp_pairwise <- data_sign %>%
  group_by(marker, stress) %>%
  pairwise_t_test(
    logvalues ~ time_point, paired = TRUE,
    comparisons = list(c("2", "3"), c("2", "4"), c("2", "5"), c("2", "6")),
    p.adjust.method = "bonferroni"
  )

# add y. position
y_positions <- data_sign %>%
  group_by(marker, stress) %>%
  summarise(max_value = max(values, na.rm = TRUE), .groups = "drop") %>%
  mutate(y.position = max_value * 1.1)  # Adjust for spacing

# Merge the y-position into the pairwise comparison results
comp_pairwise <- comp_pairwise %>%
  left_join(y_positions, by = c("marker", "stress"))

# compute effect sizes
effect_sizes_pairwise <- data_sign %>%
  group_by(marker, stress) %>%
  cohens_d(logvalues ~ time_point, paired = TRUE, 
           comparisons = list(c("2", "3"), c("2", "4"), c("2", "5"), c("2", "6")),
           ci = TRUE)

# merge effect sizes to existing data frame
comp_pairwise <- comp_pairwise %>%
  left_join(effect_sizes_pairwise %>%
              select(marker, stress, group1, group2, effsize, magnitude, conf.low, conf.high),
            by = c("marker", "stress", "group1", "group2"))

# relocate variables
comp_pairwise <- relocate(comp_pairwise, effsize, .after = p.adj.signif)
comp_pairwise <- relocate(comp_pairwise, magnitude, .after = effsize)
comp_pairwise <- relocate(comp_pairwise, conf.low, .after = magnitude)
comp_pairwise <- relocate(comp_pairwise, conf.high, .after = conf.low)

# export Excel file
write_xlsx(comp_pairwise, "sEV_comp_pairwise.xlsx")


```

# plot pairwise comparisons
```{r plot pairwise comparisons}

p1 <- ggplot(data = data_sign %>% filter(marker == "CD9"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "CD9"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD9**<br>(<i>n</i> = 19, <i>n</i> = 14)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p2 <- ggplot(data = data_sign %>% filter(marker == "CD13"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "CD13"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD13**<br>(<i>n</i> = 13, <i>n</i> = 12)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p3 <- ggplot(data = data_sign %>% filter(marker == "CD14"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "CD14"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD14**<br>(<i>n</i> = 8, <i>n</i> = 10)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p4 <- ggplot(data = data_sign %>% filter(marker == "CD16"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "CD16"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD16**<br>(<i>n</i> = 17, <i>n</i> = 13)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p5 <- ggplot(data = data_sign %>% filter(marker == "CD41"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "CD41"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD41**<br>(<i>n</i> = 12, <i>n</i> = 12)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p6 <- ggplot(data = data_sign %>% filter(marker == "CD44"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "CD44"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD44**<br>(<i>n</i> = 13, <i>n</i> = 11)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p7 <- ggplot(data = data_sign %>% filter(marker == "CD63"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD63**<br>(<i>n</i> = 11, <i>n</i> = 10)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

p8 <- ggplot(data = data_sign %>% filter(marker == "HLA.DR"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(size = 1, position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean, size = 2.5,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = 0.25, color = "black",
               position = position_dodge(.7), size = 2.5) +
  stat_pvalue_manual(comp_pairwise %>%
                       filter(marker == "HLA.DR"), label = "p.adj.signif", tip.length = 0.02, 
                     bracket.size = 1.5, size = 12, step.increase = 0.075,
                     hide.ns = TRUE, color = "black") +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**HLA.DR**<br>(<i>n</i> = 11, <i>n</i> = 14)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        text = element_text(size = 30),
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      size = 26, color = "#404040", face = "bold"
    )
  )

# Arrange the plots in 2 rows and 4 columns with labels A-H
ggarrange(p1, p2, p3, p4, p5, p6, p7, p8, 
          nrow = 2, ncol = 4, 
          labels = LETTERS[1:8],
          font.label = list(size = 64) # ,common.legend = TRUE, legend = "right"
          )
ggsave("sEV_significant.tiff", width = 4*15.4, height = 1.4*14, device='tiff', dpi=300, limitsize = FALSE)

```

# plot other sEVs
```{r plot other sEVs}

# filter nonsignificant sEVs
EVs_nonsign <- data_long_complete %>%
  select(marker) %>%
  filter(!(marker %in% EVs_sign))

EVs_nonsign <- unique(EVs_nonsign) # remove duplicates
EVs_nonsign <- as.character(EVs_nonsign$marker) # nonsignificant sEVs

# calculate n per marker
print(
  data_long_complete %>%
  filter(marker %in% EVs_nonsign) %>%
  group_by(marker, stress) %>%
  summarise(n_samples = (sum(!is.na(values))/5), .groups = "drop"),
  n = 30)


p9 <- ggplot(data = data_long_complete %>% filter(marker == "CD24"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD24**<br>(<i>n</i> = 14, <i>n</i> = 12)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p10 <- ggplot(data = data_long_complete %>% filter(marker == "CD61"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD61**<br>(<i>n</i> = 13, <i>n</i> = 13)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p11 <- ggplot(data = data_long_complete %>% filter(marker == "CD66b"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD66b**<br>(<i>n</i> = 10, <i>n</i> = 10)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p12 <- ggplot(data = data_long_complete %>% filter(marker == "CD81"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD81**<br>(<i>n</i> = 15, <i>n</i> = 17)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p13 <- ggplot(data = data_long_complete %>% filter(marker == "CD82"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD82**<br>(<i>n</i> = 15, <i>n</i> = 14)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p14 <- ggplot(data = data_long_complete %>% filter(marker == "CD90"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD90**<br>(<i>n</i> = 14, <i>n</i> = 13)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p15 <- ggplot(data = data_long_complete %>% filter(marker == "CD100"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD100**<br>(<i>n</i> = 11, <i>n</i> = 8)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p16 <- ggplot(data = data_long_complete %>% filter(marker == "CD171"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD171**<br>(<i>n</i> = 15, <i>n</i> = 14)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p17 <- ggplot(data = data_long_complete %>% filter(marker == "CD206"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD206**<br>(<i>n</i> = 11, <i>n</i> = 10)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p18 <- ggplot(data = data_long_complete %>% filter(marker == "CD227"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD227**<br>(<i>n</i> = 9, <i>n</i> = 9)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p19 <- ggplot(data = data_long_complete %>% filter(marker == "CD235a"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CD235a**<br>(<i>n</i> = 13, <i>n</i> = 13)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p20 <- ggplot(data = data_long_complete %>% filter(marker == "CTLA.4"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**CTLA.4**<br>(<i>n</i> = 15, <i>n</i> = 14)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p21 <- ggplot(data = data_long_complete %>% filter(marker == "HLA.ABC"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**HLA.ABC**<br>(<i>n</i> = 11, <i>n</i> = 12)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p22 <- ggplot(data = data_long_complete %>% filter(marker == "PD.L1"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**PD.L1**<br>(<i>n</i> = 10, <i>n</i> = 9)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

p23 <- ggplot(data = data_long_complete %>% filter(marker == "PS"), 
       aes(x = as.factor(time_point), y = values, group = as.factor(VPN), colour = as.factor(VPN))) + 
  geom_line(position = position_dodge(.05)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "point", fun = mean, shape = 18,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), colour = "black",
               geom = "line", fun = mean,
               position = position_dodge(.7)) +
  stat_summary(aes(group = as.factor(stress)), 
               fun.data = mean_se, geom = "errorbar", width = .25 ,color = "black",
               position = position_dodge(.7)) +
  labs(x = NULL, y = "objects/ml", colour = NULL) +
  scale_x_discrete(labels = c("2" = "-2 Min", "3" = "+2 Min", "4" = "+15 Min", "5" = "+30 Min", "6" = "+40 Min")) +
  scale_y_continuous(label = scientific_10) +
  ggtitle("**PS**<br>(<i>n</i> = 14, <i>n</i> = 14)") +  # Use Markdown-style formatting
  facet_wrap(vars(stress), 
             labeller = as_labeller(c("TSST" = "psychosocial stress", "Sport" = "physical stress")), scales = "fixed") +
  theme_minimal() +
  theme(legend.position = "none",
        plot.title = element_markdown(hjust = 0.5)) +
  theme(
    strip.text.x = element_text(
      color = "#404040", face = "bold"
    )
  )

ggarrange(p9, p10, p11, p12, p13, p14, p15, p16, p17, p18, p19, p20, p21, p22, p23, 
          nrow = 5, ncol = 3, 
          labels = LETTERS[1:15],
          font.label = list(size = 24) # ,common.legend = TRUE, legend = "right"
)
ggsave("sEV_nonsignificant.tiff", width = 12*1.45 , height = 4*3.2, device='tiff', dpi=300, limitsize = FALSE)

```

# rpart classification
```{r rpart classification}

# Change data format to wide
data_rpart <- pivot_wider(data_long %>% dplyr::select(VPN, time_point, stress, marker, logvalues), 
                                  names_from = marker, values_from = logvalues) # logvalues used, NO imputation following

# rename TSST and Sport
data_rpart$stress <- plyr::revalue(data_rpart$stress, c("TSST" = "psychosocial", "Sport" = "physical"))

# no imputation

## Split Data Into Training And Testing Sets
set.seed(2025)  # For reproducibility

# Create an index of 10 random numbers to split the data frame
rpart_index <- sample(1:20, 10, replace = FALSE)

# Split the data into four blocks of non duplicate VPNs
train_1 <- data_rpart %>% filter(VPN %in% rpart_index, stress == "physical")
train_2 <- data_rpart %>% filter(!(VPN %in% rpart_index), stress == "psychosocial")

test_1 <- data_rpart %>% filter(!(VPN %in% rpart_index), stress == "physical")
test_2 <- data_rpart %>% filter(VPN %in% rpart_index, stress == "psychosocial")

# # Split data into 50% training and 50% testing
train <- rbind(train_1, train_2)
test <- rbind(test_1, test_2)

# perform rpart time series analysis
rpart_model <- rpart(stress ~ CD9 + CD13 + CD14 + CD16 + CD24 + 
                         CD41 + CD44 + CD61 + CD63 + CD66b + CD81 + 
                         CD82 + CD90 + CD100 + CD171 + CD206 + CD227 + 
                         CD235a + CTLA.4 + HLA.ABC + HLA.DR + PD.L1 + PS, 
                       data = train, method = "class",
                  control = rpart.control(cp = 0.05, xval = 10))

# print results summary
print(rpart_model)
summary(rpart_model)

# prediction and confusion matrix - train data
table_rpart <- data.frame(value = train$stress, pred = predict(rpart_model, train, type = "class"))
conf_mat(table_rpart, truth = value, estimate = pred)

# metrics - train data
class_metrics_rpart <- metric_set(accuracy, sensitivity, specificity)
class_metrics_rpart(table_rpart, truth = value, estimate = pred)

# prediction and confusion matrix - test data
table_rpart2 <- data.frame(value = test$stress, pred = predict(rpart_model, test, type = "class"))
conf_mat(table_rpart2, truth = value, estimate = pred)

# metrics - test data
class_metrics_rpart2 <- metric_set(accuracy, sensitivity, specificity)
class_metrics_rpart2(table_rpart2, truth = value, estimate = pred)


```

# rpart plots
```{r rpart plots}

## Decision Tree
tiff("sEV_rpart_Tree.tiff", units = "in", width = 3, height = 3, res = 400)
rpart.plot(rpart_model, extra = 104, under = TRUE, box.palette = list("white", "white")) # extra = 106, 104 or 100?
dev.off()

## Variable Importance Plot
# Get variable importance from the model
var_imp <- data.frame(Variable = names(rpart_model$variable.importance),
                      Importance = rpart_model$variable.importance)

# Plot variable importance
ggplot(data = var_imp, aes(x = Importance, y = reorder(Variable, Importance))) +
  geom_segment( aes(xend=0, yend=Variable)) +
  geom_point( size=4, color="#2C3E50") +
  theme_bw() +
  ylab("") +
  xlab("Variable Importance")
ggsave("sEV_rpart_VarImportance.tiff", width = 3.5, height = 3, device='tiff', dpi=400)

## Class probability heatmap
# Get class probabilities
pred_probs <- predict(rpart_model, type = "prob")
pred_probs_df <- as.data.frame(pred_probs)

# Convert to long format for ggplot
pred_probs_long <- melt(pred_probs_df)

# Plot as a heatmap
ggplot(pred_probs_long, aes(x = variable, y = value)) +
  geom_tile(aes(fill = value)) +
  scale_fill_gradient(low = "#FF5252", high = "#4CAF50") +
  theme_minimal() +
  theme(legend.title=element_blank()) + # legend.title=element_blank() or legend.position = "none"
  labs(x = "Predicted", y = "Probability")
ggsave("sEV_rpart_Probability.png", width = 3.5, height = 3, device='png', dpi=400)


## ROC Curve
# calculate AUC
# Predict probabilities on training set
train_preds <- predict(rpart_model, train, type = "prob") %>%
  as.data.frame() %>%
  mutate(Class = train$stress) %>%
  rename(prob = 1)  # Assuming 1st column contains probability for the positive class

pr_train <- prediction(train_preds$prob, train$stress)
auc_train <- performance(pr_train, measure = "auc")@y.values[[1]] 
auc_train # AUC

# Predict probabilities on testing set
test_preds <- predict(rpart_model, test, type = "prob") %>%
  as.data.frame() %>%
  mutate(Class = test$stress) %>%
  rename(prob = 1)

pr_test <- prediction(test_preds$prob, test$stress)
auc_test <- performance(pr_test, measure = "auc")@y.values[[1]] 
auc_test # AUC

# calculate sensitivity and specificity
perf_rpart_train <- performance(pr_train, "sens", "spec")
perf_rpart_test <- performance(pr_test, "sens", "spec")

# plot ROC curve
tiff("sEV_rpart_ROC.tiff", units = "in", width = 3.5, height = 3, res = 400)
plot(perf_rpart_train,
     avg= "threshold",
     col= "#4CABAF",
     lwd= 3)
plot(perf_rpart_test,
     avg= "threshold",
     col="#FFD3B6",
     lwd= 3,
     add = TRUE)
legend("bottomleft",
  legend = c("training", "testing"), 
  col = c("#4CABAF", "#FFD3B6"), 
  lty = 1,
  lwd = 3,
  bty = "n", 
  cex = 0.75, 
  text.col = "black", 
  horiz = F,
  inset = c(0.05, 0.05))
legend("bottomleft", title = "    AUC",
  legend = c("0.78", "0.76"), 
  bty = "n", 
  cex = 0.75, 
  text.col = "black", 
  horiz = F,
  inset = c(0.4, 0.05))
dev.off()


## Confusion Matrix
# create data frame for plotting training data
crplot <- conf_mat(table_rpart, truth = value, estimate = pred)
crplot <- crplot$table
crplot <- as.matrix(crplot)
# crplot <- crplot/10 # Accuracy

# rename factor levels
crplot <- as.data.frame(crplot) %>%
  dplyr::rename(Predicted = Prediction, Actual = Truth)

# plot confusion matrix
ggplot(crplot, aes(x = Predicted, y = Actual, fill = Freq)) +
  geom_tile(color = "white") +
  geom_text(aes(label = Freq), size = 6, color = "black") +
  scale_fill_gradient2(low = "red", mid = "white", high = "#4CABAF", midpoint = 0) +
  labs(title = "Overall Accuracy: 78%") +
  theme_minimal() +
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5, face = "bold"))
ggsave("sEV_rpart_Matrix_Train.tiff", width = 3.5, height = 3, device='tiff', dpi=400)


# create data frame for plotting test data
crplot2 <- conf_mat(table_rpart2, truth = value, estimate = pred)
crplot2 <- crplot2$table
crplot2 <- as.matrix(crplot2)
# crplot2 <- crplot/10 # Accuracy

# rename factor levels
crplot2 <- as.data.frame(crplot2) %>%
  dplyr::rename(Predicted = Prediction, Actual = Truth)

# plot confusion matrix
ggplot(crplot2, aes(x = Predicted, y = Actual, fill = Freq)) +
  geom_tile(color = "white") +
  geom_text(aes(label = Freq), size = 6, color = "black") +
  scale_fill_gradient2(low = "red", mid = "white", high = "#FFD3B6", midpoint = 0) +
  labs(title = "Overall Accuracy: 76%") +
  theme_minimal() +
  theme(legend.position="none",
        plot.title = element_text(hjust = 0.5, face = "bold"))
ggsave("sEV_rpart_Matrix_Test.tiff", width = 3.5, height = 3, device='tiff', dpi=400)


```