9 PCB Dataset B

2017 & 2019 data sets, 72 pups, file name = 170101, 190137

9.1 Libraries

library(tidyverse)
library(gridExtra)
library(ggpubr)
library(car)
library(arm)
library(ggmosaic)
library(EnvStats)

9.2 Data

pcb <- 
  read.csv("Output Files/cleaned_pcb_B.csv") %>% 
  rename("Sample" = 1) 

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

data <- 
  merge (pcb, meta) %>% 
  mutate(year = as.numeric(year),
         molt.stage = as.factor(molt.stage),
         iav = as.factor(iav),
         location = as.factor(location),
         sex = as.factor(sex))

lod <- 
  read.csv("Output Files/lod_B.csv")

9.3 Exploratory

9.3.1 Sample Description

table(data$year)

## 
## 2017 2018 2019 
##   51    1   20

table(data$sex)

## 
##  F  M 
## 38 28

table(data$location)

## 
## Great Point     Monomoy    Muskeget 
##           2          46          24

table(data$molt.stage)

## 
##  II III  IV   V 
##   1   9  14  31

table(data$iav)

## 
## neg pos 
##  46  26

table(data$pcbbin)

## 
##  0  1 
## 34 38

9.3.2 Descriptive Stats

9.3.2.1 Mean Median, Range

data %>% 
  filter(sumpcb>0) %>% 
  summarise(mean=mean(sumpcb), 
            se=sd(sumpcb)/(sqrt(length(sumpcb[sumpcb>0]))),
            median=median(sumpcb),
            max=max(sumpcb),
            min=min(sumpcb))

##       mean       se median   max  min
## 1 112.9237 10.14842  93.25 321.4 27.9

datalong <- 
  data %>% 
  pivot_longer(2:22, names_to = "pcb", values_to = "conc")

9.3.2.2 Congener Total

sumcongener <- 
  data.frame(colSums(data[,2:22])) %>% 
  mutate(congener=factor(row.names(.), levels=row.names(.))) %>% 
  rename(sum=c(1)) 

sumcongenerplot <- 
  ggplot(sumcongener, aes(x=congener, y=sum)) +
  geom_bar(stat="identity") + 
  labs(x="PCB Congener", y="Sum PCBs (ng/g wet weight)") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=90, hjust=1), 
        text=element_text(size=12))
sumcongenerplot

ggsave("Figures/sumpcb_B.jpg", sumcongenerplot, height=5, width=10, units="in")

9.3.2.3 Congener Mean

meancongener <- 
  data %>% 
  dplyr::select(2:22) %>% 
  pivot_longer(everything(), names_to="pcb", values_to="conc") %>%
  filter(!conc==0) %>% 
  group_by(pcb) %>% 
  summarize(meanconc=mean(conc),
            se=sd(conc)/(sqrt(length(conc)))) 

meancongenerplot <-
  ggplot(meancongener, aes(x=pcb, y=meanconc)) +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin = meanconc - se, ymax = meanconc + se), width = 0.2) +
  labs(x="PCB Congener", y="Mean PCBs (ng/g wet weight)") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=90, hjust=1), 
        text=element_text(size=12))
meancongenerplot

ggsave("Figures/meanpcb_B.jpg", meancongenerplot, height=5, width=8, units="in")

9.3.2.4 Year Total

yearsumcongener <-
  datalong %>% 
  group_by(year, pcb) %>% 
  summarize(sum = sum(conc))

## `summarise()` has grouped output by 'year'. You
## can override using the `.groups` argument.

yearsumcongener

## # A tibble: 63 × 3
## # Groups:   year [3]
##     year pcb      sum
##    <dbl> <chr>  <dbl>
##  1  2017 pcb101   0  
##  2  2017 pcb105 150. 
##  3  2017 pcb118   0  
##  4  2017 pcb128   0  
##  5  2017 pcb138   0  
##  6  2017 pcb153 194. 
##  7  2017 pcb170  87.2
##  8  2017 pcb18    0  
##  9  2017 pcb180   0  
## 10  2017 pcb183   0  
## # ℹ 53 more rows

ggplot(yearsumcongener, aes(x=pcb, y=sum)) +
       geom_bar(stat="identity") + 
       labs(x="PCB Congener", y="sum PCBs (ng/g wet weight)") +
       facet_wrap(~year, ncol = 4) +
       theme(axis.text.x = element_text(angle=90))

yearsum <- 
  yearsumcongener %>% 
  group_by(year) %>% 
  summarise(sum=sum(sum)) 

yearsumplot <-
  ggplot(yearsum, aes(x=factor(year), y=sum)) +
  geom_bar(stat="identity") + 
  labs(x="Year", y="sum PCBs (ng/g wet weight)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90),
        panel.grid = element_blank())
yearsumplot

ggsave("Figures/yearsumpcb_B.jpeg", yearsumplot, height=5, width=7, units="in")

9.3.2.5 Year Mean

yearmean <-
  datalong %>% 
  filter(!conc==0) %>% 
  group_by(year) %>% 
  summarize(mean = mean(conc),
            se = sd(conc)/(sqrt(length(conc)))) 
yearmean

## # A tibble: 3 × 3
##    year  mean    se
##   <dbl> <dbl> <dbl>
## 1  2017  73.8  4.55
## 2  2018  80.9 28.9 
## 3  2019  88.1  7.70

yearmeanplot <- 
  ggplot(yearmean, aes(x=factor(year), y=mean)) +
  geom_bar(stat="identity") + 
  geom_errorbar(aes(ymin = mean-se, ymax = mean+se), 
                width = 0.2) +
  labs(x="PCB Congener", y="Mean PCBs (ng/g wet weight)") +
  theme_bw() +
  theme(axis.text.x = element_text(angle=90),
        panel.grid = element_blank())
yearmeanplot

ggsave("Figures/yearmeanpcb_B.jpeg", yearmeanplot, height=5, width=7, units="in")

9.3.2.6 What congeners are present in each animal (separated by year/iav status)?

Additive/synergistic effects - is # congeners associated with IAV?

ggplot(data=datalong, aes(x=Sample, y=conc, fill=pcb)) +
  geom_col() +
  labs(x="Seal ID", y="PCB Concentration") +
  facet_wrap(~iav, drop=TRUE, scales="free", ncol=1) +
  theme_bw() +
  theme(panel.grid=element_blank(), 
        axis.text.x=element_text(angle=90, vjust=0.5, hjust=1))

9.3.2.7 How many seals were each congener detected in?

datalong %>%
  group_by(pcb) %>%
  summarise(count=sum(conc>0)) %>% 
  ggplot(data=., aes(x=pcb, y=count)) +
  geom_col() +
  labs(x="Congener", y="Number of Seals") + 
  geom_text(aes(label=count, vjust=-0.75)) +
  theme_bw() +
  theme(panel.grid=element_blank(), 
        axis.text.x=element_text(angle=90, vjust=0.5, hjust=1))

9.3.3 Table of Mean, SE, Ranges for each congener

Gives some errors becauause some PCBs have no detections

summetrics <- 
  data.frame(
    LOD = "NA",
    sampsize = table(data$pcbbin)[2],
    sampperc = (table(data$pcbbin)[2]/nrow(data))*100,
    mean = mean(data$sumpcb[data$sumpcb>0]),
    se = sd(data$sumpcb)/(sqrt(length(data$sumpcb[data$sumpcb>0]))),
    median = median(data$sumpcb[data$sumpcb>0]),
    min = min(data$sumpcb[data$sumpcb>0]),
    max = max(data$sumpcb[data$sumpcb>0])
  )

sumpcb <- 
  data.frame(
    X = "sumPCB",
    X1 = "NA",
    X2 = "NA")

lod <- 
  lod %>% 
  rbind(sumpcb)

pcbtable <-
  data.frame(matrix(nrow=21, ncol=0)) %>% 
  mutate(LOD = lod$X1[1:21],
         sampsize=sapply(data[,2:22], function(x) length(x[x>0])),
         sampperc=sapply(data[,2:22], function(x) (length(x[x>0])/127)*100),
         mean=sapply(data[,2:22], function(x) {mean(x[x>0])}),
         se=sapply(data[,2:22], function(x) sd(x)/(sqrt(length(x[x>0])))),
         median=sapply(data[,2:22], function(x) {median(x[x>0])}),
         min=sapply(data[,2:22], function(x) {min(x[x>0])}),
         max=sapply(data[,2:22], max)) %>% 
  rbind(summetrics) %>%
  mutate(across(everything(), as.numeric)) %>% 
  sapply(., function(x) round(x, digits=2)) %>% 
  data.frame() %>% 
  add_column(pcb = lod$X, .before=1) %>% 
  mutate(sampnum=paste0(sampsize, " (", sampperc, "%)"),
         meanse=paste0(mean, " ± ", se),
         medrange=paste0(median, " (", min, " - ", max, ")"))

write.csv(pcbtable, "Output Files/pcb_summarystats_B.csv", row.names=FALSE)

9.4 What influences PCBs?

9.4.1 Data

# Remove pcbs with 0 detections
data_tidy <- 
  pcb %>% 
  column_to_rownames(var = "Sample") %>% 
  dplyr::select(1:21) %>% 
  select_if(function (x) (sum(x) > 0)) %>% 
  rownames_to_column(var = "Sample") %>% 
  merge(., data[,c(1,23:35)], by = "Sample")

# Subset only pups with PCBs
dataP <- 
  data_tidy %>% 
  filter(sumpcb > 0) %>% 
  mutate(logpcb=log(sumpcb))

9.4.2 Year

9.4.2.1 Presence/Absence

Significant - declining presence of PCBs from 2017 to 2019 (only 1 pup in 2018)

# Logistic regression
yearmodel <- 
  glm(pcbbin ~ year, data=data_tidy, family=binomial) 

Anova(yearmodel)

## Analysis of Deviance Table (Type II tests)
## 
## Response: pcbbin
##      LR Chisq Df Pr(>Chisq)  
## year    5.269  1    0.02171 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Make predictions from model
logit_preds <- data.frame(stats::predict(yearmodel, type="link", se.fit=TRUE))
logit_preds$lwr <- logit_preds$fit + 1.96 * logit_preds$se.fit
logit_preds$upr <- logit_preds$fit - 1.96 * logit_preds$se.fit
real_preds <- apply(logit_preds, 2, invlogit)

# Combine with original data & plot
data.frame(data_tidy, real_preds) %>% 
  ggplot(., aes(x = year, y = fit)) +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.25) +
  geom_line() +
  labs(x="Year", y="Probality of PCBs") +
  theme_bw() +
  theme(panel.grid = element_blank(),
        text=element_text(size=13))

yearprop <- 
  data_tidy %>% 
  mutate(pcbbin=ifelse(pcbbin=="1", "Present", "Absent"),
         pcbbin=factor(pcbbin, levels=c("Present", "Absent"))) %>% 
  ggplot(data=.) +
  geom_mosaic(aes(x=product(pcbbin, year), fill=pcbbin)) +
  scale_fill_manual(values=c("Present"="gray60", "Absent"="lightgray"), name="PCB Status") +
  scale_y_continuous(labels = scales::percent) +
  labs(y="Proportion of Samples") +
  theme_classic() +
  theme(legend.position = "inside",
        legend.position.inside = c(0.83, 0.15), 
        axis.title.x=element_blank(),
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
yearprop

9.4.2.2 Concentration

Not significant

# Pearson's Correlation
cor.test(dataP$year, dataP$sumpcb, method="pearson")

## 
##  Pearson's product-moment correlation
## 
## data:  dataP$year and dataP$sumpcb
## t = 1.6467, df = 36, p-value = 0.1083
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  -0.06009896  0.53876742
## sample estimates:
##       cor 
## 0.2646694

# Plot
yearconc <-
  ggplot(dataP, aes(x=year, y=sumpcb)) +
  geom_point(size=0.5, position=position_jitter(width=0.2)) +
  geom_smooth(method = "lm", colour="black") +
  labs(y="log(ΣPCB ng/g wet weight)") +
  scale_x_continuous(breaks = seq(2017,2019, by=1)) +
  stat_n_text() +
  theme_classic() +
  theme(axis.title.x=element_blank(), 
        text=element_text(size=13),
        axis.text.x=element_text(angle=45, hjust=1))
yearconc

## `geom_smooth()` using formula = 'y ~ x'

9.4.2.3 Combine

yearplots <- grid.arrange(yearprop, yearconc, ncol=2)

## `geom_smooth()` using formula = 'y ~ x'

ggsave("Figures/year-pcbbin-sumpcb_B.jpeg", yearplots, width=12, height=6, units="in")

9.4.3 Sex

9.4.3.1 Presence/Absence

Not significant

tryCatch({
  chisq.test(table(data_tidy$sex, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$sex, data_tidy$pcbbin))
})

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(data_tidy$sex, data_tidy$pcbbin)
## X-squared = 0.55827, df = 1, p-value = 0.455

9.4.3.2 Concentration

Not significant

# Concentration
t.test(sumpcb ~ sex, data=dataP)

## 
##  Welch Two Sample t-test
## 
## data:  sumpcb by sex
## t = -1.2827, df = 16.812, p-value = 0.217
## alternative hypothesis: true difference in means between group F and group M is not equal to 0
## 95 percent confidence interval:
##  -86.95215  21.23549
## sample estimates:
## mean in group F mean in group M 
##        100.8667        133.7250

9.4.4 Location

9.4.4.1 Presence/Absence

Not significant

tryCatch({
  chisq.test(table(data_tidy$location, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$location, data_tidy$pcbbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data_tidy$location, data_tidy$pcbbin)
## p-value = 0.1966
## alternative hypothesis: two.sided

9.4.4.2 Concentration

Not significant

t.test(sumpcb ~ location, data=dataP)

## 
##  Welch Two Sample t-test
## 
## data:  sumpcb by location
## t = 0.15452, df = 18.659, p-value = 0.8789
## alternative hypothesis: true difference in means between group Monomoy and group Muskeget is not equal to 0
## 95 percent confidence interval:
##  -43.97707  50.97841
## sample estimates:
##  mean in group Monomoy mean in group Muskeget 
##               113.9370               110.4364

9.4.5 Molt Stage

9.4.5.1 Presence/Absence

Not significant

tryCatch({
  chisq.test(table(data_tidy$molt.stage, data_tidy$pcbbin))
}, warning = function(w) {
  fisher.test(table(data_tidy$molt.stage, data_tidy$pcbbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data_tidy$molt.stage, data_tidy$pcbbin)
## p-value = 0.8681
## alternative hypothesis: two.sided

9.4.5.2 Concentration

Not significant

summary(aov(sumpcb ~ molt.stage, data=dataP))

##             Df Sum Sq Mean Sq F value Pr(>F)
## molt.stage   3   5535    1845   0.428  0.735
## Residuals   24 103441    4310               
## 10 observations deleted due to missingness

9.5 PCB Dataset B - Cytokines

n = 37 pups, all measured in 2022

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

9.7 Data

pcb <- 
  read.csv("Output Files/cleaned_pcb_B.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% pcb$Sample) %>% 
  merge(., pcb, by="Sample") %>% 
  merge(., meta, by = "Sample")

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

9.8 Sample Description

table(data$year)

## 
## 2017 2018 2019 
##   24    1   12

table(data$sex)

## 
##  F  M 
## 19 17

table(data$location)

## 
## Great Point     Monomoy    Muskeget 
##           2          17          18

table(data$molt.stage)

## 
## III  IV   V 
##   3  10  13

table(data$iav)

## 
## neg pos 
##  12  25

table(data$pcbbin)

## 
##  PCB Absent PCB Present 
##          19          18

table(data$analysis.year.x)

## 
## 2022 
##   37

9.9 PCB - IAV interaction

Not significant

table(data$iav, data$pcbbin)

##      
##       PCB Absent PCB Present
##   neg          6           6
##   pos         13          12

tryCatch({
  chisq.test(table(data$iav, data$pcbbin))
}, warning = function(w) {
  fisher.test(table(data$iav, data$pcbbin))
})

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(data$iav, data$pcbbin)
## X-squared = 2.6315e-31, df = 1, p-value = 1

9.10 Cytokine Exploratory

9.10.1 Cytokine detections

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	37
IL.7	26
IFNg	22
IL.2	20
KC.like	9
IL.10	2
IL.6	2
IL.15	1
IP.10	1
GM.CSF	0
IL.8	0
MCP.1	0
TNFa	0

9.10.2 Remove cytokines with low detection rates

Filtering out: IL-10, IL-6, IL-15, IP-10, GM-CSF, IL-8, MCP-1, TNFa

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

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

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

9.10.3 PCA

9.10.3.1 PCA & data

cyto_pca <- 
  prcomp(data[,2:6], scale = 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="")

9.10.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("black", "#008ba2")) +
  theme_bw() +
  theme(panel.grid=element_blank())

9.10.3.3 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())

9.10.3.4 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())

9.10.3.5 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())

## Too few points to calculate an ellipse

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

9.10.3.6 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: Removed 1 row containing missing values or
## values outside the scale range (`geom_path()`).

9.10.3.7 Sex

Maybe some separation here?

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()`).

9.10.3.8 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())

## Too few points to calculate an ellipse

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

9.11 PERMANOVA

9.11.1 Cytokine Presence/Absence

9.11.1.1 Create similarity matrix (Sorensen)

Dissimilarity = 1 - Sorensen

databin_dist <-
  databin %>% 
  select(2:6) %>% 
  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.3068044

## [1] 0

## [1] 0.6666667

## [1] 0.3333333

9.11.1.2 Run PERMANOVA

Not significant

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.11627 0.04902 1.8041 0.2146
## Residual 35  2.25572 0.95098              
## Total    36  2.37199 1.00000

9.11.1.3 Homogeneity of group dispersions

Not significant

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.00851 0.0085092  0.8197 0.3715
## Residuals 35 0.36332 0.0103806

boxplot(disper)

9.11.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)

9.11.2 Cytokine Concentration

9.11.2.1 Run PERMANOVA

Significant

data_dist <-
  data %>% 
  select(2:6) %>% 
  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.21877 0.07125 2.6852  0.029 *
## Residual 35  2.85154 0.92875                
## Total    36  3.07032 1.00000                
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

9.11.2.2 Homogeneity of group dispersions

Not significant

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.00879 0.0087902  0.5512 0.4628
## Residuals 35 0.55813 0.0159467

boxplot(disper)

9.11.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)

9.12 PLS-DA

9.12.1 Data

pls_x <-
  data[,2:6]

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

9.12.2 PLSDA

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

9.12.3 Visualization

9.12.3.1 Set x and y labels

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

9.12.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)

9.12.3.3 Plot

pcb_plsda_ggplot <- 
  ggplot(plsda_ggplot, 
       aes(x=x, 
           y=y,
           col=pcb)) +
  geom_point() +
  stat_ellipse() +
  labs(x = plsda_labelx,
       y = plsda_labely) +
  scale_color_manual("PCB Status",
                     values=c("black", "#9d9d9d")) +
  theme_bw() +
  theme(panel.grid = element_blank(),
        legend.position = "bottom")

pcb_plsda_ggplot

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

9.12.4 Cytokine Contributions

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

	comp1	comp2
IFNg	-0.1645339	-0.2727963
KC.like	-0.0346947	-0.3662218
IL.2	0.2126143	-0.4767680
IL.7	0.3285638	0.7202198
IL.18	0.9047463	-0.2131654

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

plotVar(plsda_model)

9.13 CART

9.13.1 Prepare data

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

9.13.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 = 2
  # minbucket = 11
  # minsplit = 14
  # classif.ce = 0.25

9.13.3 Run model with optimized parameters

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

# plot optimized tree
rpart.plot(cart_model)

cart_model$variable.importance 

##     IL.18      IL.7      IL.2 
## 5.0400579 1.2600145 0.6300072

  # IL-18 is top cytokines

9.13.4 Plot

9.13.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))

varimpplot <- 
  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"))

9.13.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"))))

9.13.4.3 Proportions

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

##              
##                 1
##   PCB Absent  100
##   PCB Present 100

il18pres <- 
  ggplot(data=databin_plot) +
  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(y = "Proportion of Samples", title="IL-18") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title.x = element_blank(), 
        legend.position="bottom",
        text = element_text(size=10),
        plot.title = element_text(hjust = 0.5))

9.13.4.4 Concentrations

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="log(Concentration (pg/mL))", 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))

9.13.4.5 Combine together

cyto <- grid.arrange(il18pres, il18, ncol=2)

all <- grid.arrange(varimpplot, cyto,
                    ncol=1, 
                    heights=c(1, 2))

ggsave("Figures/cytoimportance_B.jpeg", all, width = 6, height = 6, units="in")

9.14 Supplemental Figures

9.14.1 Cytokine Presence/Absence

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

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")) +
  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))

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, il7_prop,
               kc_prop, il18_prop,
               ncol=3,
               left=textGrob("Proportion of Samples", rot=90),
               bottom=textGrob("Influenza A Virus Infection Status"))

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

9.14.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))

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

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

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, il7,
               kc, il18,
               ncol=3,
               left=textGrob("log(Cytokine Concentration (pg/mL))", rot=90),
               bottom=textGrob("PCBs"))

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

9.15 PCB Dataset B - IAV

n = 72 pups

9.16 Load required libraries

library(tidyverse)
library(car)
library(ggmosaic)
library(gridExtra)
library(EnvStats)
library(rstatix)

9.17 Data

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

pcb <- 
  read.csv("Output Files/cleaned_pcb_B.csv") %>% 
  dplyr::rename("Sample" = 1)

viro <- 
  read.csv("Input Files/Hg_virology_2023.csv", sep=',', strip.white=TRUE) %>% 
  select(Tag.ID., year, IAV, IAV.Nasal, IAV.Conj, IAV.Rectal) %>% 
  setNames(c("Sample", "year", "IAV", "Nasal", "Conj", "Rectal"))

viro_two <-
  viro %>% 
  filter(Sample == "204" | 
         Sample == "268")

viro_tidy <- 
  viro %>% 
  filter(Sample %in% pcb$Sample) %>% 
  rbind(., viro_two)

9.18 Merge data & configure variables

data <- 
  merge(meta, pcb, by="Sample") %>% 
  select(Sample, iav, sumpcb, pcbbin, year, location, sex, molt.stage) %>% 
  mutate(iavbin = as.numeric(ifelse(iav=="neg", "0", "1")))

swabdata <- 
  merge(viro_tidy, pcb, by="Sample")

9.19 IAV ~ Variables

9.19.1 Year

Not significant

yeariav <- glm(iavbin ~ year, data=data, family=binomial) 
Anova(yeariav)

## Analysis of Deviance Table (Type II tests)
## 
## Response: iavbin
##      LR Chisq Df Pr(>Chisq)
## year  0.75586  1     0.3846

9.19.2 Sex

Not significant

tryCatch({
  chisq.test(table(data$sex, data$iavbin))
}, warning = function(w) {
  fisher.test(table(data$sex, data$iavbin))
})

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(data$sex, data$iavbin)
## X-squared = 0, df = 1, p-value = 1

9.19.3 Location

Not significant

tryCatch({
  chisq.test(table(data$location, data$iavbin))
}, warning = function(w) {
  fisher.test(table(data$location, data$iavbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data$location, data$iavbin)
## p-value = 0.05626
## alternative hypothesis: two.sided

9.19.4 Molt Stage

Not significant

tryCatch({
  chisq.test(table(data$molt.stage, data$iavbin))
}, warning = function(w) {
  fisher.test(table(data$molt.stage, data$iavbin))
})

## 
##  Fisher's Exact Test for Count Data
## 
## data:  table(data$molt.stage, data$iavbin)
## p-value = 0.1247
## alternative hypothesis: two.sided

9.20 IAV ~ PCB p/a

Not significant

tryCatch({
  chisq.test(table(data$iav, data$pcbbin))
}, warning = function(w) {
  fisher.test(table(data$iav, data$pcbbin))
})

## 
##  Pearson's Chi-squared test with Yates' continuity correction
## 
## data:  table(data$iav, data$pcbbin)
## X-squared = 0.011928, df = 1, p-value = 0.913

9.21 IAV ~ PCB conc

Not significant

# Fit model
model <- glm(iavbin ~ sumpcb, data=data, family=binomial)
Anova(model)

## Analysis of Deviance Table (Type II tests)
## 
## Response: iavbin
##        LR Chisq Df Pr(>Chisq)
## sumpcb  0.82336  1     0.3642

9.22 Figures

9.22.1 Presence/Absence

IAV+ pups less likely to have PCBs present

iavmosaic <- 
  data %>% 
  mutate(iav = ifelse(iav == "pos", "IAV+", "IAV-"),
         iav = factor(iav, levels=c("IAV+", "IAV-")),
         pcbbin = ifelse(pcbbin == "1", "Present", "Absent"),
         pcbbin = factor(pcbbin, levels=c("Present", "Absent"))) %>%
  ggplot(data=.) +
  geom_mosaic(aes(x=product(pcbbin, iav), fill=pcbbin)) +
  scale_fill_manual(values = c("Absent"="lightgray", "Present"="gray60"), 
                    name = "PCB Status") +
  scale_y_continuous(labels = scales::percent) +
  labs(y = "Proportion of Samples") +  
  theme_classic() +
  theme(legend.position=c(0.83, 0.15),
        axis.title.x = element_blank(),
        text = element_text(size=13),
        axis.text = element_text(size=13))

9.22.2 Concentration

…and groups have roughly equal concentrations of PCBs

iavbox <- 
  data %>% 
  mutate(iav=ifelse(iav=="pos", "IAV+", "IAV-"),
         iav=factor(iav, levels=c("IAV+", "IAV-"))) %>% 
  ggplot(., aes(x=iav, y=log(sumpcb))) +
  geom_boxplot () +
  geom_point(position = position_jitter(width = 0.2)) +
  labs(y="log(ΣPCB Concentration)") +
  stat_n_text(y.pos=3) +
  theme_classic() +
  theme(axis.title.x=element_blank(), 
        text=element_text(size=13),
        axis.text = element_text(size=13))

9.22.3 Combine

pcbiav <- grid.arrange(iavmosaic, iavbox, ncol=2)

## Warning: Removed 34 rows containing non-finite outside
## the scale range (`stat_boxplot()`).

## Warning: Removed 34 rows containing non-finite outside
## the scale range (`stat_n_text()`).

ggsave("Figures/pcb-iav_B.jpeg", pcbiav, width=10, height=5, units="in")