3 BSC - WB Validation

3.1 Load required libraries

library(tidyverse)
library(ggpubr)
library(ggmosaic)
library(gridExtra)
library(kableExtra)
library(knitr)
library(coin)
library(ggrepel)

3.2 Load & configure data

lod <- 
  read.csv("Output Files/170101_LOD.csv", row.names=1)

whole <- 
  read.csv("Output Files/cleaned_pcb2_wholeblood.csv") %>% 
  rename("id" = "X") %>% 
  pivot_longer(!id, names_to="pcb", values_to="wholeblood") %>% 
  merge(., lod, by="pcb")
  

spot <- 
  read.csv("Output Files/cleaned_pcb_B.csv") %>% 
  rename("id" = "X") %>% 
  select(1:22) %>% 
  pivot_longer(!id, names_to="pcb", values_to="bloodspotcard") %>% 
  merge(., lod, by="pcb")

data <- 
  merge(whole, spot, by = c("id", "pcb", "LOD")) %>% 
  mutate(wlod=ifelse(wholeblood>LOD, "above", "below"),
         slod=ifelse(bloodspotcard>LOD, "above", "below"))

meta <- 
  read.csv("Input Files/metadata_tidy.csv") %>% 
  filter(Sample %in% data$id)

3.3 Data Summary

3.3.1 Metadata

table(meta$year)

## 
## 2017 
##   30

table(meta$sex)

## 
##  F  M 
## 13 12

table(meta$location)

## 
## Monomoy 
##      30

table(meta$molt.stage)

## 
##  II III  IV   V 
##   1   9   8  11

table(meta$iav)

## 
## neg pos 
##  21   9

table(meta$iavser)

## 
## neg pos 
##  27   3

3.3.2 Number of Samples

sample_summary <- 
  data %>% 
  group_by(id) %>% 
  summarize(wbsum = sum(wholeblood),
            bscsum = sum(bloodspotcard)) %>% 
  mutate(wbpres = ifelse(wbsum > 0, "Present", "Absent"),
         bscpres = ifelse(bscsum > 0, "Present", "Absent"))

table(sample_summary$wbpres)

## 
##  Absent Present 
##       9      21

table(sample_summary$bscpres)

## 
##  Absent Present 
##      10      20

3.3.3 ~PCB Congeners

3.3.3.1 Configure Data

data %>% 
  select(!ends_with(c("lod", "LOD"))) %>% 
  mutate(pcb=factor(pcb, levels=lod$pcb)) %>% 
  group_by(pcb) %>% 
  summarize(`WB SumPCB`=sum(wholeblood),
            `WB Count`=sum(wholeblood>0),
            `BSC SumPCB`=sum(bloodspotcard), 
            `BSC Count`=sum(bloodspotcard>0)) %>% 
  kable() %>% 
  kable_styling("basic")

pcb	WB SumPCB	WB Count	BSC SumPCB	BSC Count
pcb18	118.7	1	0.0	0
pcb28	65.9	1	86.5	1
pcb44	148.7	4	384.0	6
pcb49	378.9	7	280.8	5
pcb52	404.2	6	472.7	7
pcb66	234.9	3	72.1	1
pcb87	287.3	2	0.0	0
pcb101	111.9	1	0.0	0
pcb105	75.7	2	88.7	2
pcb118	0.0	0	0.0	0
pcb128	0.0	0	0.0	0
pcb138	0.0	0	0.0	0
pcb153	74.7	1	140.3	2
pcb170	244.5	2	87.2	1
pcb180	0.0	0	0.0	0
pcb183	0.0	0	0.0	0
pcb184	0.0	0	0.0	0
pcb187	107.3	1	131.2	1
pcb195	0.0	0	0.0	0
pcb206	0.0	0	0.0	0
pcb209	119.8	2	204.4	2

congener_numdet <- 
  data %>% 
  select(!ends_with(c("lod", "LOD"))) %>% 
  mutate(pcb=factor(pcb, levels=lod$pcb)) %>% 
  group_by(pcb) %>% 
  summarize(wb_count = sum(wholeblood>0),
            bsc_count = sum(bloodspotcard>0)) %>% 
  mutate(wb_prop = wb_count/30, 
         bsc_prop = bsc_count/30) %>% 
  pivot_longer(c(wb_prop, bsc_prop), 
               names_to = "method", values_to = "proportion") 

3.3.3.2 Plot

congener_prop <- 
  ggplot(congener_numdet, aes(x = pcb, y = proportion, fill = method)) +
  geom_col(position = "dodge") +
  scale_fill_manual(values = c("#b4b4b4", "#7eaaac"), 
                    labels = c("DBS", "WB"),
                    name = "Method") +
  labs(y = "Proportion of Samples", 
       x = "PCB Congener") +
  theme_bw() +
  theme(panel.grid= element_blank(),
       axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "inside", 
        legend.position.inside = c(0.9, 0.86),
       text = element_text(size = 12))
congener_prop

congener_meanconc <- 
  data %>% 
  mutate(pcb=factor(pcb, levels=lod$pcb)) %>% 
  pivot_longer(cols=c("wholeblood", "bloodspotcard"), 
               names_to="method", values_to="concentration") %>% 
  filter(!concentration == "0") %>% 
  group_by(pcb, method) %>% 
  summarize(mean = mean(concentration),
            sd = sd(concentration)) %>% 
  ggplot(., aes(x = pcb, y = mean, fill = method)) +
  geom_bar(position="dodge", stat="identity") +
  scale_fill_manual(name="Sample Matrix", 
                    labels=c("DBS", "WB"), values=c("#b4b4b4", "#7eaaac")) +
  labs(y="Mean Concentration (ng/g wet weight)") +
  theme_bw() +
  theme(panel.grid= element_blank(),
        axis.title.x = element_blank(), 
        axis.text.x=element_text(angle=45, hjust=1),
        legend.position="none")

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

congener_meanconc

3.3.4 ~Individual Seals

3.3.4.1 Configure data

whole_summary_sample <- 
  data %>% 
  select(!ends_with(c("lod", "LOD", "card"))) %>% 
  group_by(id) %>% 
  summarize(whole_sumpcb=sum(wholeblood),
            whole_counts=sum(wholeblood>0))
sum(whole_summary_sample$whole_sumpcb>0)

## [1] 21

spot_summary_sample <- 
  data %>% 
  select(!ends_with(c("lod", "LOD", "blood"))) %>% 
  group_by(id) %>% 
  summarize(spot_sumpcb=sum(bloodspotcard),
            spot_counts=sum(bloodspotcard>0))
sum(spot_summary_sample$spot_counts>0)

## [1] 20

3.3.4.2 Plot

seal_summary_plot <- 
  data %>% 
  mutate(id=factor(id)) %>% 
  pivot_longer(cols=c("wholeblood", "bloodspotcard"), 
               names_to="method", values_to="concentration") %>% 
  group_by(pcb, id, method) %>% 
  summarize(mean = mean(concentration), 
            sd = sd(concentration)) %>% 
  ggplot(., aes(x=id, y=mean, fill=method)) +
  geom_bar(position="dodge", stat="identity") +
  geom_errorbar(aes(ymin = mean - sd, ymax = mean + sd)) +
  scale_fill_manual(name="Sample Matrix", 
                    labels=c("DBS", "WB"), values=c("#b4b4b4", "#7eaaac")) +
  labs(x="Grey Seal Pup ID", 
       y="Mean Concentration (ng/g wet weight)") +
  theme_classic() +
  theme(axis.text.x=element_text(angle=45, hjust=1),
        legend.position="bottom")

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

seal_summary_plot

ggsave("Figures/bsc_validation_seal_suppl.jpg", seal_summary_plot, 
       width=8, height=6, units="in")

3.4 t-test

3.4.1 Number of Detections

No sig difference in number of detections per congener ** but not normally distributed - count data, can’t use t-test?

hist(congener_numdet$wb_count)

hist(congener_numdet$bsc_count)

shapiro.test(congener_numdet$wb_count)

## 
##  Shapiro-Wilk normality test
## 
## data:  congener_numdet$wb_count
## W = 0.76763, p-value = 9.68e-07

shapiro.test(congener_numdet$bsc_count)

## 
##  Shapiro-Wilk normality test
## 
## data:  congener_numdet$bsc_count
## W = 0.67404, p-value = 2.196e-08

congener_wilcox <-
  congener_numdet %>% 
  select(pcb, bsc_count, wb_count) %>% 
  pivot_longer(!pcb, names_to = "method", values_to = "count") %>% 
  mutate(method = factor(method))
  
wilcox_test(data = congener_wilcox,
            count ~ method,
            paired = TRUE)

## 
##  Asymptotic Wilcoxon-Mann-Whitney Test
## 
## data:  count by method (bsc_count, wb_count)
## Z = -1.1551, p-value = 0.2481
## alternative hypothesis: true mu is not equal to 0

cor.test(congener_numdet$wb_count, congener_numdet$bsc_count)

## 
##  Pearson's product-moment correlation
## 
## data:  congener_numdet$wb_count and congener_numdet$bsc_count
## t = 11.921, df = 40, p-value = 9.631e-15
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.7921168 0.9360038
## sample estimates:
##       cor 
## 0.8833761

validation_detections <- 
  ggplot(congener_numdet, aes(x=wb_count, y = bsc_count)) +
  geom_point() +
  geom_smooth(method="lm", color="#7eaaac") +
  labs(x = "WB PCB Detections", 
       y = "BSC PCB Detections") +
  theme_classic() +
  theme(text=element_text(size=12))
validation_detections

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

mean(congener_numdet$wb_count)

## [1] 1.571429

sd(congener_numdet$wb_count)

## [1] 1.964848

min(congener_numdet$wb_count)

## [1] 0

max(congener_numdet$wb_count)

## [1] 7

mean(congener_numdet$bsc_count)

## [1] 1.333333

sd(congener_numdet$bsc_count)

## [1] 2.079712

min(congener_numdet$bsc_count)

## [1] 0

max(congener_numdet$bsc_count)

## [1] 7

3.5 PCB Concentrations

data_above <- 
  data %>% 
  filter(wlod=="above" & slod=="above")

cor.test(data_above$wholeblood, data_above$bloodspotcard)

## 
##  Pearson's product-moment correlation
## 
## data:  data_above$wholeblood and data_above$bloodspotcard
## t = 4.0278, df = 13, p-value = 0.001435
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.3765319 0.9100127
## sample estimates:
##      cor 
## 0.745081

validation_plot <- 
  ggplot(data_above, aes(x=wholeblood, y = bloodspotcard)) +
  geom_point(size = 2) +
  geom_smooth(method="lm", color="#7eaaac") +
  stat_cor() +
  stat_regline_equation(label.y = 135) +
  geom_text_repel(data = data_above, aes(label = pcb), size = 2.5) +
  scale_y_continuous(breaks = c(40, 60, 80, 100, 120, 140)) +
  scale_x_continuous(breaks = c(40, 60, 80, 100, 120, 140)) +
  labs(x = "WB PCB (ng/g wet weight)", 
       y = "DBS PCB (ng/g wet weight)") +
  theme_classic() +
  theme(text=element_text(size=12))
validation_plot

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

## Create Figure for Manuscript

manuscript <- 
  grid.arrange(congener_prop, validation_plot, 
               ncol = 2, widths = c(2,1.5))

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

ggsave("Figures/bsc_validation_manuscript.jpg", manuscript, 
       width = 10, height = 4, units = "in")

3.6 Create production-ready figures

congener_prop_prod <-
  congener_prop +
  theme(text = element_text(size = 9),
        legend.position.inside = c(0.9, 0.86),
        legend.key.size = unit(0.25, "cm"))
  
validation_plot_prod <-
  ggplot(data_above, aes(x=wholeblood, y = bloodspotcard)) +
  geom_point(size = 0.75) +
  geom_smooth(method="lm", color="#7eaaac", size = 0.5) +
  stat_cor(size = 3, label.y = 130) +
  stat_regline_equation(label.y = 140, size = 3) +
  geom_text_repel(data = data_above, aes(label = pcb), size = 2, 
                  point.size = 0.05) +
  scale_y_continuous(limits = c(20, 140), breaks = c(40, 60, 80, 100, 120, 140)) +
  scale_x_continuous(limits = c(20, 140), breaks = c(40, 60, 80, 100, 120, 140)) +
  labs(x = "WB PCB (ng/g wet weight)", 
       y = "DBS PCB (ng/g wet weight)") +
  theme_classic() +
  theme(text=element_text(size=9))

## Warning: Using `size` aesthetic for lines was deprecated
## in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to
## see where this warning was generated.

production <- 
  grid.arrange(congener_prop_prod, validation_plot_prod, 
               ncol = 2, widths = c(2,1.5))

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

ggsave("Figures/Figure_1.jpg", production, 
       width = 7500, height = 3000, units = "px", dpi = 1000)