6 PCB ~ Cytokines

6.1 Libraries

library(mixOmics)
library(tidyverse)
library(ggfortify)
library(vegan)
library(car)
library(EnvStats)
library(gridExtra)
library(ggpubr)
library(mlr3)
library(mlr3tuning)
library(rpart)
library(rpart.plot)
library(grid)
library(ggmosaic)
library(kableExtra)

6.2 Data

pcbmax <- 
  read.csv("Output Files/cleaned_pcb_maxLOD.csv") %>% 
  rename("Sample" = c(1)) %>% 
  mutate(pcbbin = ifelse(pcbbin=="0", "PCB Absent", "PCB Present"))

meta <- 
  read.csv("Input Files/metadata_tidy.csv")

data <- 
  read.csv("../hg-cyto/Output Files/cleaned_cytokine_all.csv") %>% 
  filter(Sample %in% pcbmax$Sample) %>% 
  merge(., pcbmax, by="Sample") %>% 
  merge(., meta, by = "Sample")

databin <- 
  read.csv("../hg-cyto/Output Files/cleaned_cytokine_bin_all.csv") %>% 
  filter(Sample %in% pcbmax$Sample) %>%  
  merge(., pcbmax, by="Sample") %>% 
  merge(., meta, by = "Sample")

uncleaned_cyto <-
  read.csv("../hg-cyto/Output Files/uncleaned_cyto.csv") %>% 
  select(!X) %>% 
  filter(Sample %in% data$Sample)
write.csv(uncleaned_cyto, "Output Files/uncleaned_cyto.csv")

6.3 PCB - IAV interaction

Is there a difference in the proportion of IAV+ pups across PCB groups? Should IAV be considered as a confounding factor?
** Not sig

table(data$iav, data$pcbbin)

##      
##       PCB Absent PCB Present
##   neg         37          11
##   pos         32           7

chisq.test(data$iav, data$pcbbin)

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  data$iav and data$pcbbin
## X-squared = 0.091686, df = 1, p-value = 0.762

6.4 Cytokine Exploratory

6.4.1 Cytokine detectionss

detect <- 
  data %>% 
  pivot_longer(., cols=c(2:14), names_to="cytokine", values_to="conc") %>% 
  group_by(cytokine) %>% 
  summarize(n = sum(conc > 0)) %>% 
  arrange(desc(n))

detect %>% 
  kable() %>% 
  kable_styling("basic")

cytokine	n
IL.18	87
IL.7	58
IFNg	49
IL.2	42
IL.10	28
KC.like	22
IL.15	16
IL.8	13
IL.6	7
GM.CSF	1
IP.10	1
MCP.1	0
TNFa	0

6.4.2 Remove cytokines with low detection rates

detect <- 
  detect %>% 
  filter(n < 2)

data <- 
  data %>% 
  select(!detect$cytokine)

databin <-
  databin %>% 
  select(!detect$cytokine)

6.4.3 PCA

6.4.3.1 PCA & data

cyto_pca <- 
  prcomp(data[,2:10], scale = TRUE, center = TRUE)

cyto_pca_data <-
  data.frame(
    x = cyto_pca$x[,1],
    y = cyto_pca$x[,2],
    pcbbin = factor(data$pcbbin),
    analysis.year = factor(data$analysis.year.x),
    iav = data$iav,
    iavser = data$iavser,
    location = data$location,
    year = factor(data$year), 
    sex = data$sex,
    molt.stage = data$molt.stage
  )

cyto_pca_labelx <-
  paste("PC1 (", 
        round(abs(summary(cyto_pca)$importance[2,1]*100), digits=0),
        "%)",
        sep="")
cyto_pca_labely <-
  paste("PC2 (", 
        round(abs(summary(cyto_pca)$importance[2,2]*100), digits=0),
        "%)",
        sep="")

6.4.3.2 PCBS

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=pcbbin)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("PCBs",
                     values=c("#9d9d9d", "#008ba2")) +
  theme_bw() +
  theme(panel.grid=element_blank())

6.4.3.3 Analysis Year

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=analysis.year)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("Analysis Year",
                     values=c("#008ba2","#9d9d9d")) +
  theme_bw() +
  theme(panel.grid=element_blank())

6.4.3.4 IAV

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=iav)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("IAV",
                     values=c("#9d9d9d", "#008ba2")) +
  theme_bw() +
  theme(panel.grid=element_blank())

6.4.3.5 IAV Serology

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=iavser)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("IAV Serology",
                     values=c("#9d9d9d", "#008ba2")) +
  theme_bw() +
  theme(panel.grid=element_blank())

6.4.3.6 Location

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=location)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("Location",
                     values=c("#9d9d9d", "#008ba2", "black")) +
  theme_bw() +
  theme(panel.grid=element_blank())

6.4.3.7 Year

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=year)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("Year",
                     values=c("#9d9d9d", "#008ba2", "orchid1", "darkseagreen", 
                              "slategray2", "magenta4", "black")) +
  theme_bw() +
  theme(panel.grid=element_blank())

## Too few points to calculate an ellipse

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

## Warning: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).

6.4.3.8 Sex

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=sex)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("Sex",
                     values=c("#9d9d9d", "#008ba2")) +
  theme_bw() +
  theme(panel.grid=element_blank())

## Too few points to calculate an ellipse

## Warning: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).

6.4.3.9 Molt Stage

ggplot(cyto_pca_data, 
       aes(x=x, 
           y=y,
           col=molt.stage)) +
  geom_point() +
  stat_ellipse() +
  labs(x=cyto_pca_labelx,
       y = cyto_pca_labely) +
  scale_color_manual("Molt Stage",
                     values=c("#9d9d9d", "#008ba2", "black")) +
  theme_bw() +
  theme(panel.grid=element_blank())

## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure

6.5 PERMANOVA

6.5.1 Cytokine Presence/Absence

6.5.1.1 Create similarity matrix (Sorensen)

Dissimilarity = 1 - Sorensen

databin_dist <-
  databin %>% 
  select(2:10) %>% 
  betadiver(., method=11)

plot(databin_dist)

# Dissimilarity summary stats
mean(1-databin_dist); min(1-databin_dist); max(1-databin_dist); median(1-databin_dist)

## [1] 0.4049534

## [1] 0

## [1] 0.8

## [1] 0.4285714

6.5.1.2 Run PERMANOVA

R2 = partial R2 = factor (pcbs) is explaining at least R2 * 100% of variation in variables (cytokines)
** Not sig

adonis2((1-databin_dist) ~ pcbbin, 
        data=databin, 
        permutations=4999)

## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 4999
## 
## adonis2(formula = (1 - databin_dist) ~ pcbbin, data = databin, permutations = 4999)
##          Df SumOfSqs      R2      F Pr(>F)  
## Model     1   0.2656 0.03076 2.6973 0.0582 .
## Residual 85   8.3690 0.96924                
## Total    86   8.6346 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

6.5.1.3 Homogeneity of group dispersions

** Not sig

disper <- 
  betadisper((1-databin_dist), group=databin$pcbbin, type="centroid")

anova(disper)

## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq   Mean Sq F value Pr(>F)
## Groups     1 0.02255 0.0225523  2.7397 0.1016
## Residuals 85 0.69969 0.0082317

boxplot(disper)

6.5.1.4 Visualization

Principal coordinates analysis

# Use 1 - Sorensen's for dissimilarity 
databin_pco <- 
  cmdscale((1-databin_dist), eig=TRUE)

plot(databin_pco$points, type="n", 
     cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordiellipse(ord = databin_pco, 
            groups = factor(data$pcbbin), 
            display = "sites", 
            col = c("grey", "darkseagreen"),
            lwd = 2,
            label = TRUE)

plot(databin_pco$points, type="n", 
     cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordispider(ord = databin_pco,
           groups = factor(data$pcbbin),
           display = "sites",
           col = c("grey", "darkseagreen"),
           lwd=2,
           label=TRUE)

6.5.2 Cytokine Concentration

6.5.2.1 Run PERMANOVA

** Not sig

data_dist <-
  data %>% 
  select(2:10) %>% 
  vegdist(., method="bray")

adonis2(data_dist ~ pcbbin, 
        data=data, 
        permuations=4999)

## Permutation test for adonis under reduced model
## Permutation: free
## Number of permutations: 999
## 
## adonis2(formula = data_dist ~ pcbbin, data = data, permuations = 4999)
##          Df SumOfSqs      R2     F Pr(>F)
## Model     1   0.2141 0.01425 1.229   0.26
## Residual 85  14.8055 0.98575             
## Total    86  15.0196 1.00000

6.5.2.2 Homogeneity of group dispersions

** Not sig

disper <- 
  betadisper(data_dist, group=data$pcbbin, type="centroid")

anova(disper)

## Analysis of Variance Table
## 
## Response: Distances
##           Df  Sum Sq  Mean Sq F value Pr(>F)
## Groups     1 0.02607 0.026067  1.0693  0.304
## Residuals 85 2.07218 0.024379

boxplot(disper)

6.5.2.3 Visualization

# Use Bray-Curtis dissimilarity matrix
data_pco <- 
  cmdscale(data_dist, eig=TRUE)

plot(data_pco$points, type="n", 
     cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordiellipse(ord = data_pco, 
            groups = factor(data$pcbbin), 
            display = "sites", 
            col = c("grey", "darkseagreen"),
            lwd = 2,
            label = TRUE)

plot(data_pco$points, type="n", 
     cex.lab=1.5, cex.axis=1.1, cex.sub=1.1)
ordispider(ord = data_pco,
           groups = factor(data$pcbbin),
           display = "sites",
           col = c("grey", "darkseagreen"),
           lwd=2,
           label=TRUE)

6.6 PLS-DA

6.6.1 Data

pls_x <-
  data[,2:10]

pls_y <- 
  data %>% 
  mutate(pcbbin=as.factor(pcbbin)) %>% 
  select(pcbbin)

6.6.2 PLSDA

plsda_model <- 
  plsda(pls_x,
        pls_y$pcbbin,
        scale=TRUE)

6.6.3 Visualization

6.6.3.1 Set x and y labels

plsda_labelx <-
  paste("LV1 (", 
        round(abs(plsda_model$prop_expl_var$X[1]*100), digits=1),
        "%)",
        sep="")
plsda_labely <-
  paste("LV2 (", 
        round(abs(plsda_model$prop_expl_var$X[2]*100), digits=1),
        "%)",
        sep="")

6.6.3.2 Create plot data frame

plsda_ggplot <- 
  data.frame(x = plsda_model$variates$X[,1], 
             y = plsda_model$variates$X[,2],
             pcb = plsda_model$Y)

6.6.3.3 Plot

pcb_plsda_ggplot <- 
  ggplot(plsda_ggplot, 
         aes(x=x, 
             y=y,
             col=pcb)) +
  geom_point(size = 2) +
  stat_ellipse(linewidth = 1.5) +
  labs(x = plsda_labelx,
       y = plsda_labely) +
  scale_color_manual("PCB Status",
                     values=c("#b4b4b4","#7eaaac")) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "inside",
        legend.position.inside = c(0.12, 0.9),
        text = element_text(size = 16))

pcb_plsda_ggplot

ggsave("Figures/pcb-cyto_plsda.jpeg", pcb_plsda_ggplot, 
       width=8, height=6, units="in")

6.6.3.4 Plot Presentation

pcb_plsda_pres <- 
  pcb_plsda_ggplot +
  theme(legend.position = "bottom",
        text = element_text(size = 18))
ggsave("Figures/pcb-cyto_plsda_presentation.jpeg", pcb_plsda_pres, 
       width=8, height=6, units="in")

6.6.3.5 Production Plots

pcb_plsda_pres_prod <- 
 ggplot(plsda_ggplot, 
         aes(x=x, 
             y=y,
             col=pcb)) +
  geom_point(size = 0.05) +
  stat_ellipse(linewidth = 0.2) +
  labs(x = plsda_labelx,
       y = plsda_labely) +
  scale_color_manual("PCB Status",
                     values=c("#b4b4b4","#7eaaac")) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "inside",
        legend.position.inside = c(0.12, 0.90),
        text = element_text(size = 5),
        legend.background = element_blank(),
        legend.text = element_text(size = 3),
        legend.title = element_text(size = 4, margin = margin(b = 1)),
        legend.key.size = unit(0.15, "cm"))
pcb_plsda_pres_prod

ggsave("Figures/Figure_4.jpeg", pcb_plsda_pres_prod, 
       width = 2500, height = 1875, units="px", dpi = 1000)

6.6.4 Cytokine Contributions

plsda_model$loadings$X %>% 
  data.frame() %>% 
  select(1:2) %>% 
  arrange(comp1) %>% 
  kable() %>% 
  kable_styling("basic")

	comp1	comp2
KC.like	-0.4577447	0.1665368
IL.2	-0.4049218	0.0832380
IL.7	-0.2650913	0.2187258
IL.6	-0.2520067	0.1077794
IFNg	-0.1864299	-0.0175947
IL.18	0.0162000	0.4373388
IL.10	0.1605509	0.8412679
IL.15	0.1765294	-0.0080581
IL.8	0.6330637	0.0806724

biplot(plsda_model, ind.names=FALSE, legend.title="PCB Status")

plotVar(plsda_model)

6.7 CART

6.7.1 Prepare data

data_cart <- 
  data %>% 
  select(2:10, pcbbin) %>% 
  mutate(pcbbin=factor(pcbbin))

6.7.2 Cross Validation

# Creating task and learner
task <- as_task_classif(pcbbin ~ ., 
                        data = data_cart)
        
task <- task$set_col_roles(cols="pcbbin",
                           add_to="stratum")

# min pcb group = 18
learner <- lrn("classif.rpart",
               predict_type = "prob",
               maxdepth = to_tune(2, 5),
               minbucket = to_tune(1, 18),
               minsplit = to_tune(1, 30)
               )

# Define tuning instance - info ~ tuning process
instance <- ti(task = task,
               learner = learner,
               resampling = rsmp("cv", folds = 10),
               measures = msr("classif.ce"),
               terminator = trm("none")
               )

# Define how to tune the model
tuner <- tnr("grid_search", 
             batch_size = 10
             )

# Trigger the tuning process
#tuner$optimize(instance)

# optimal values: 
  # maxdepth = 3
  # minbucket = 7
  # minsplit = 7
  # classif.ce = 0.1855159

6.7.3 Run model with optimized parameters

cart_model <- rpart(formula=pcbbin ~ .,
                    data=data_cart,
                    method="class",
                    maxdepth=3,
                    minbucket=7,
                    minsplit=7)

# plot optimized tree
rpart.plot(cart_model)

cart_model$variable.importance 

##      IL.7      IL.8     IL.18     IL.10     IL.15   KC.like      IFNg      IL.6 
## 5.0186666 4.0368368 2.3390381 1.8603130 1.6147347 1.3507454 1.0096066 0.9354804 
##      IL.2 
## 0.5457439

  # IL-7, IL-8 are top cytokines

6.7.4 Manuscript Plot

Variable importance: the sum of the goodness of split measures for each split for which it was the primary variable, plus goodness * (adjusted agreement) for all splits in which it was a surrogate. In the printout these are scaled to sum to 100 and the rounded values are shown, omitting any variable whose proportion is less than 1%. (https://cran.r-project.org/web/packages/rpart/vignettes/longintro.pdf)

6.7.4.1 Variable Importance

varimp <- 
  data.frame(imp=cart_model$variable.importance) %>% 
  rownames_to_column() %>% 
  rename("variable" = rowname) %>% 
  arrange(imp) %>%
  mutate(variable=gsub("\\.","-",variable),
         variable = forcats::fct_inorder(variable))

varimp_plot <- 
  ggplot(varimp) + 
  geom_segment(aes(x = variable, 
                   y = 0, 
                   xend = variable, 
                   yend = imp), 
               linewidth = 0.5, 
               alpha = 0.7) +
  geom_point(aes(x = variable, 
                 y = imp), 
             size = 1, 
             show.legend = F,
             col="black") +
  labs(x="Cytokine", 
       y="Variable Importance") +
  coord_flip() +
  theme_classic() +
  theme(text=element_text(size=10, color="black"))

6.7.4.2 Data Prep

databin_plot <-
  databin %>% 
  mutate(across(2:10, \(x) gsub("1", "Present", x)),
         across(2:10, \(x) gsub("0", "Absent", x)),
         across(2:10, \(x) factor(x, levels=c("Present", "Absent"))))

6.7.4.3 Proportions

prop.table(table(databin$pcbbin, databin$IL.7), margin=1)*100

##              
##                      0        1
##   PCB Absent  39.13043 60.86957
##   PCB Present 11.11111 88.88889

prop.table(table(databin$pcbbin, databin$IL.8), margin=1)*100

##              
##                      0        1
##   PCB Absent  89.85507 10.14493
##   PCB Present 66.66667 33.33333

il7_prop <- 
  ggplot(data=databin_plot) +
  geom_mosaic(aes(x=product(IL.7,pcbbin), fill=IL.7)) +
  scale_fill_manual(values=c("Present" = "gray60", "Absent" = "lightgray"),
                    name="Cytokine") +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-7") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="bottom",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))
il7_prop

cyto_leg <- get_legend(il7_prop)

il7_prop <-
  il7_prop +
  theme(legend.position="none")

il8_prop <- 
  ggplot(data=databin_plot) +
  geom_mosaic(aes(x=product(IL.8,pcbbin), fill=IL.8)) +
  scale_fill_manual(values=c("Present" = "gray60", "Absent" = "lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-8") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(), 
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))
il8_prop

6.7.4.4 Concentrations

il7 <- 
  data %>% 
  select(pcbbin, IL.7) %>% 
  filter(!IL.7==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.7))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL) +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
         panel.border=element_rect(color="black", fill=NA, linewidth=0.5))
il7

il8 <- 
  data %>% 
  select(pcbbin, IL.8) %>% 
  filter(!IL.8==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.8))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL) +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
         panel.border=element_rect(color="black", fill=NA, linewidth=0.5))
il8

6.7.4.5 Combine together

prop_grid <- grid.arrange(il7_prop, il8_prop, ncol=2,
                          left=textGrob("Proportion of Samples",
                                        rot=90, gp=gpar(fontsize=10)))

prop_grid <- grid.arrange(prop_grid, cyto_leg,
                          ncol=1, 
                          heights=c(3, 0.3))

conc_grid <- grid.arrange(il7, il8, ncol=2, 
                          left=textGrob("log(Concentration (pg/mL))", 
                                        rot=90, gp=gpar(fontsize=10)))

cyto_grid <- grid.arrange(prop_grid, conc_grid,
                          ncol=1, 
                          heights=c(4,3.5),
                          bottom=textGrob("PCB Absent/Present", gp=gpar(fontsize=10)))

jpeg("Figures/cytoimportance.jpeg", width=6, height=7, units="in", res=300)
final_grid <- 
  grid.arrange(varimp_plot, cyto_grid, 
               ncol=1, heights=c(2,4)) 
grid.polygon(x=c(0, 0, 1, 1),
             y=c(0, 0.67, 0.67, 0),
             gp=gpar(fill=NA))
grid.rect(x=unit(0.5, "npc"), y = unit(0.958, "npc"),
          width=unit(0.948, "npc"), height=unit(0.05, "npc"), 
          gp=gpar(lwd=1.5, col="gray30", fill=NA))
dev.off()

## png 
##   2

6.8 Supplemental Figures

6.8.1 Cytokine Presence/Absence

kc_prop <-
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(KC.like,pcbbin), fill=KC.like)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray"),
                    name="Cytokine") +
  scale_y_continuous(labels = scales::percent) +
  labs(title="KC-like") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il6_prop <-
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.6,pcbbin), fill=IL.6)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray"),
                    name="Cytokine") +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-6") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

ifng_prop <-
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IFNg,pcbbin), fill=IFNg)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray"),
                    name="Cytokine") +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IFNg") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il2_prop <- 
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.2,pcbbin), fill=IL.2)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-2") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il7_prop <- 
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.7,pcbbin), fill=IL.7)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-7") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il8_prop <-
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.8,pcbbin), fill=IL.8)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-8") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il10_prop <- 
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.10,pcbbin), fill=IL.10)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-10") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il15_prop <- 
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.15,pcbbin), fill=IL.15)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-15") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

il18_prop <-
  databin_plot %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(IL.18,pcbbin), fill=IL.18)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray")) +
  scale_y_continuous(labels = scales::percent) +
  labs(title="IL-18") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(),
        axis.title.y = element_blank(),
        legend.position="none",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

cyto_prop_grid <-
  grid.arrange(ifng_prop, il2_prop, il6_prop,
               il7_prop, il8_prop, il15_prop, 
               kc_prop, il10_prop, il18_prop,
               left=textGrob("Proportion of Samples", rot=90),
               bottom=textGrob("PCB Status"))

cyto_prop_grid <- 
  grid.arrange(cyto_prop_grid, cyto_leg, 
               ncol=1, 
               heights=c(9,0.5))

ggsave("Figures/supplemental_cyto_prop_grid.jpeg", cyto_prop_grid, 
       width=8, height=8, units="in")

6.8.2 Cytokine Concentration - Log

ifng <- 
  data %>% 
  select(pcbbin, IFNg) %>% 
  filter(!IFNg==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IFNg))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IFNg") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il2 <- 
  data %>% 
  select(pcbbin, IL.2) %>% 
  filter(!IL.2==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.2))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-2") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il6 <- 
  data %>% 
  select(pcbbin, IL.6) %>% 
  filter(!IL.6==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.6))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-6") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il7 <- 
  data %>% 
  select(pcbbin, IL.7) %>% 
  filter(!IL.7==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.7))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-7") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il8 <- 
  data %>% 
  select(pcbbin, IL.8) %>% 
  filter(!IL.8==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.8))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-8") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il15 <- 
  data %>% 
  select(pcbbin, IL.15) %>% 
  filter(!IL.15==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.15))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-15") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

kc <- 
  data %>% 
  select(pcbbin, KC.like) %>% 
  filter(!KC.like==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(KC.like))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "KC-like") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il10 <- 
  data %>% 
  select(pcbbin, IL.10) %>% 
  filter(!IL.10==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.10))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-10") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

il18 <- 
  data %>% 
  select(pcbbin, IL.18) %>% 
  filter(!IL.18==0) %>% 
  ggplot(., aes(x=pcbbin, y=log(IL.18))) +
  geom_boxplot(color="black", fill="gray60", lwd=0.25,
               outliers=FALSE) +
  geom_point(size=0.75, position=position_jitter(width=0.2)) +
  labs(x=NULL, y=NULL, title = "IL-18") +
  stat_n_text(size=3) +
  theme_classic() +
  theme(text = element_text(size=10),
        panel.border=element_rect(color="black", fill=NA, linewidth=0.5),
        plot.title = element_text(hjust=0.5))

cyto_conc_grid <-
  grid.arrange(ifng, il2, il6,
               il7, il8, il15, 
               kc, il10, il18,
               ncol=3,
               left=textGrob("log(Cytokine Concentration (pg/mL))", rot=90),
               bottom=textGrob("PCBs"))

ggsave("Figures/supplemental_cyto_conc_grid.jpeg", cyto_conc_grid, 
       width=8, height=8, units="in")