2 Tidy PCB Data

2.1 Libraries

library(mixOmics)
library(tidyverse)

2.2 Metadata

# My metadata
meta_me <- 
  read.csv("Input Files/metadata_edit.csv")

# Wendy's metadata
meta_wendy <- 
  read.csv("Input Files/Hg_virology_2023.csv") %>% 
  select(Tag.ID., year, field.site, IAV, IAV.serology, sex, molt) %>% 
  setNames(c("Sample", "year", "location", "iav", "iavser", "sex", "molt.stage")) %>% 
  filter(Sample %in% meta_me$Sample) %>% 
   mutate(iav = ifelse(iav == "N", "neg", 
                      ifelse(iav == "NEG", "neg", 
                             ifelse(iav == "P", "pos", "NA"))),
         iavser = ifelse(iavser == "N", "neg", 
                      ifelse(iavser == "P", "pos", "NA")))

# Compare metadata
comp <- 
  merge(meta_me, meta_wendy, by="Sample")

unique(comp$iav.x == comp$iav.y)

## [1] TRUE

unique(comp$iavser.x == comp$iavser.y)

## [1] TRUE

unique(comp$year.x == comp$year.y)

## [1] TRUE

unique(comp$location.x == comp$location.y)

## [1]  TRUE FALSE

unique(comp$sex.x == comp$sex.y)

## [1]  TRUE FALSE    NA

unique(comp$molt.stage.x == comp$molt.stage.y)

## [1]  TRUE    NA FALSE

loc <- 
  comp %>% 
  select(Sample, year.x, location.x, location.y) %>% 
  mutate(check = ifelse(location.x==location.y, "y", "nope"))

sex <-
  comp %>% 
  select(Sample, year.x, sex.x, sex.y) %>% 
  mutate(check = ifelse(sex.x==sex.y, "y", "nope"))

molt <- 
  comp %>% 
  select(Sample, year.x, molt.stage.x, molt.stage.y) %>% 
  mutate(check = ifelse(molt.stage.x == molt.stage.y, "y", "nope"))

extra <-
  meta_me %>% 
  filter(Sample == "204-268") %>% 
  select(Sample, iav, iavser, year, location, sex, molt.stage)

meta_tidy <-
  comp %>% 
  select(Sample, iav.x, iavser.x, year.x, location.x, sex.x, molt.stage.x) %>%
  setNames(c("Sample", "iav", "iavser", "year", "location", "sex", "molt.stage")) %>% 
  mutate(molt.stage = ifelse(row_number() == 136, "III", molt.stage)) %>% 
  rbind(extra)

write.csv(meta_tidy, file="Input Files/metadata_tidy.csv", row.names=FALSE)

2.3 Max LOD

2.3.1 Read in Raw Data & Pre-process

pcb1 <- read.csv("Input Files/Contaminants/160134_PCB_blood_spot_card.csv")
pcb1 <- pcb1[10:31,c(3,6:40)]
colnames(pcb1) <- pcb1[1,]
pcb1$`Field ID`[22] <- paste("LOD")
pcb1$`Field ID`[17] <- paste("204-268")
pcb1$analysis.year <- "2016"

pcb2 <- read.csv("Input Files/Contaminants/170101_PCB_blood_spot_card.csv")
pcb2 <- pcb2[10:66,c(2,5:25)]
colnames(pcb2) <- pcb2[1,]
pcb2$`Field ID`[57] <- paste("LOD")
pcb2$analysis.year <- "2017"


pcb3 <- read.csv("Input Files/Contaminants/190137_PCBs.csv")
pcb3 <- pcb3[10:31,c(2,5:25)]
colnames(pcb3) <- pcb3[1,]
pcb3$`Field ID`[22] <- paste("LOD")
pcb3$analysis.year <- "2019"


pcb4 <- read.csv("Input Files/Contaminants/220264_UMaine PCBs.csv")
pcb4 <- pcb4[10:44,c(2,5:25)]
colnames(pcb4) <- pcb4[1,]
pcb4$`Field ID`[35] <- paste("LOD")
pcb4$analysis.year <- "2022"


pcb5 <- read.csv("Input Files/Contaminants/230045_UMaine PCB.csv")
pcb5 <- pcb5[10:14,c(2,5:25)]
colnames(pcb5) <- pcb5[1,]
pcb5$`Field ID`[5] <- paste("LOD")
pcb5$`Field ID` <- gsub("Hg","", pcb5$`Field ID`)
pcb5$analysis.year <- "2023"

2.3.2 Combine data

# Combine datasets
data <- 
  bind_rows(pcb1,pcb2,pcb3,pcb4,pcb5)

# Get rid of PCBs not measured in all datasets (has NAs)
data <- 
  data %>% 
  select_if(~ !any(is.na(.)))

# Create dataframe with all LODs
LOD <- 
  data %>% 
  filter(`Field ID` == "LOD") %>% 
  select(!analysis.year)

2.3.3 Exclude samples

• 271: PCBs not actually quantified
• 1348: duplicate PCB quantifications, not sure which to use
• 2016 Tag: Not sure which animal this is

# Exclude samples: 271, 1348 x2, 2016 Tag, all header rows, all LODs
data <- filter(data, 
               !`Field ID`=='271' 
               & !`Field ID`=="1348"
               & !`Field ID`=="2016 Tag"
               & !`Field ID`=="LOD"
               & !`Field ID`=="Field ID")

2.3.4 Shorten PCB names to IUPAC number

iupac <- as.data.frame(colnames(data)) %>% 
  setNames("pcb") %>% 
  filter(!pcb=="analysis.year")

for (i in 1:nrow(iupac)) {
  name <- iupac[i,]
  {if(grepl('\\(',name)) {
    newname <- gsub("[\\(\\)]", "", regmatches(name, gregexpr("\\(.*?\\)", name))[[1]])
  }
    else{newname<-name}
  }
  iupac[i,] <- newname
}
iupac$pcb <- paste("pcb", iupac$pcb, sep="")

# Rename column names in both dataframes
colnames(data)[2:17] <- iupac$pcb[2:17]
colnames(LOD)[2:17] <- iupac$pcb[2:17]

2.3.5 Remove PCBs below max LOD across datasets

# Change "ND" to 0's
for (i in 1:nrow(data)) {
  for (j in 1:ncol(data)) {
    {if (data[i,j] == "ND") {
      data[i,j] <- 0
    }
      else{next}
    }
  }
}

# Get rid of Field ID column
data <- 
  data %>% 
  column_to_rownames(var="Field ID")

LOD <- as.data.frame(LOD[,-1])

# Make PCB & LOD concentrations numeric                                         
data[,1:16] <- sapply(data[,1:16], as.numeric)
LOD[,1:16] <- sapply(LOD[,1:16], as.numeric)

# Find max LOD for each PCB
LOD["max",] <- apply(LOD[,1:16], 2, max)
LOD

##     pcb18 pcb44 pcb66 pcb87 pcb105 pcb118 pcb128 pcb138 pcb153 pcb180 pcb183
## 1    13.3  26.7  13.3  26.7   26.7   13.3   53.3   53.3   26.7  266.7  266.7
## 2    25.0  25.0  50.0  75.0   25.0   75.0   75.0  100.0   50.0  100.0   75.0
## 3    25.0  25.0  50.0  75.0   25.0   75.0   75.0  100.0   50.0  100.0   75.0
## 4     3.2   0.7   0.7   1.0    1.3    3.3    3.5    3.2    1.9    3.0    2.8
## 5     3.2   0.7   0.7   1.0    1.3    3.3    3.5    3.2    1.9    3.0    2.8
## max  25.0  26.7  50.0  75.0   26.7   75.0   75.0  100.0   50.0  266.7  266.7
##     pcb184 pcb187 pcb195 pcb206 pcb209
## 1    133.3  266.7   66.7   26.7   66.7
## 2    100.0  100.0   25.0   25.0   10.0
## 3    100.0  100.0   25.0   25.0   10.0
## 4      1.3    1.1    1.9    2.8    2.3
## 5      1.3    1.1    1.9    2.8    2.3
## max  133.3  266.7   66.7   26.7   66.7

# Raw sumPCB per sample
data_raw <- 
  data %>% 
  mutate(sum=rowSums(data[,1:16])) %>% 
  rownames_to_column(var="sample") %>% 
  select(sample, sum)

# Get rid of PCB values below the highest LOD for that PCB
for (i in 1:nrow(data)) {
  for (j in 1:16) {
    value <- data[i,j]
    {if (value < LOD[6,j]) {
      data[i,j] <- 0
    }
      else{next}
      }
  }
}

LOD <-
  LOD %>% 
  t() %>% 
  data.frame()

write.csv(LOD, "Output Files/lod_all.csv")

2.3.6 Create variables for analysis

# Create variables: sumpcb, pcbbin, co/noncoplanar, groups 1-3
data$sumpcb <- rowSums(data[,1:16])
data$pcbbin <- ifelse(data$sumpcb > 0, 1, 0)
data <- data %>% mutate (coplanar = pcb105+pcb118)
data <- data %>% mutate (noncoplanar = pcb18+pcb44+pcb66+pcb87+pcb128+pcb138+pcb153+pcb184+pcb195+pcb206+pcb209)
data$cobin <- ifelse(data$coplanar > 0, 1, 0)
data$noncobin <- ifelse(data$noncoplanar > 0, 1, 0)

2.3.7 Write out cleaned data for future use

write.csv(x=data, file = "Output Files/cleaned_pcb_maxLOD.csv")

2.4 PCB Dataset A

2.4.1 Tidy

Exclude sample 271 (see Max LOD for reason)

# Dataset specific LOD
lod1 <- 
  pcb1 %>% 
  filter(`Field ID`=="LOD") %>% 
  select(!`Field ID` &
         !analysis.year) %>% 
  t() %>% 
  data.frame()

# Exclude sample 271, all header rows, all LODs
pcb1 <- filter(pcb1, !`Field ID`=='271' 
               & !`Field ID`=="LOD"
               & !`Field ID`=="Field ID")

2.4.2 Shorten PCB names to IUPAC number

iupac1 <- 
  as.data.frame(colnames(pcb1)) %>% 
  setNames("pcb")

for (i in 1:nrow(iupac1)) {
  name <- iupac1[i,]
  {if(grepl('\\(',name)) {
    newname <- gsub("[\\(\\)]", "", regmatches(name, gregexpr("\\(.*?\\)", name))[[1]])
  }
    else{newname<-name}
  }
  iupac1[i,] <- newname
}

iupac1 <- 
  iupac1 %>% 
  mutate(pcbname = paste("pcb", pcb, sep="."),
         pcbname = gsub("pcb.PCB.47.PCB.48.PCB.49.PCB.52", "pcb.47.48.49.52", pcbname),
         pcbname = gsub("pcb.PCB.28.PCB.31", "pcb.28.31", pcbname),
         pcbname = gsub("pcb.PCB.90.PCB.101", "pcb.90.101", pcbname),
         pcbname = gsub("pcb.PCB.170.PCB.190", "pcb.170.190", pcbname))


colnames(pcb1)[2:36] <- iupac1$pcbname[2:36]
rownames(lod1) <- iupac1$pcb[2:36]

2.4.3 Remove PCBs below LOD

# Change "ND" to 0's
pcb1 <-
  pcb1 %>% 
  mutate(across(everything(), function(x) gsub("ND", "0", x)))

# Get rid of Field ID column
pcb1 <- 
  pcb1 %>% 
  column_to_rownames(var="Field ID")

# Get rid of PCB values below the LOD for that PCB
for (i in 1:nrow(pcb1)) {
  for (j in 1:35) {
    value <- pcb1[i,j]
    {if (value < lod1[j,1]) {
      pcb1[i,j] <- 0
    }
      else{next}
      }
  }
}

write.csv(lod1, "Output Files/lod_A.csv")

2.4.4 Create variables for analysis

# Create variables: sumpcb, pcbbin, co/noncoplanar
pcb1[,1:35] <- sapply(pcb1[,1:35], as.numeric)

pcb1$sumpcb <- rowSums(pcb1[,1:35])
pcb1$pcbbin <- ifelse(pcb1$sumpcb > 0, 1, 0)

pcb1 <- pcb1 %>% mutate(coplanar = `pcb.105` + `pcb.118` + `pcb.156` + `pcb.157`)
pcb1 <- pcb1 %>% mutate(noncoplanar = sumpcb - coplanar)

pcb1$cobin <- ifelse(pcb1$coplanar > 0, 1, 0)
pcb1$noncobin <- ifelse(pcb1$noncoplanar > 0, 1, 0)

2.4.5 Write out cleaned data for future use

write.csv(x=pcb1, file = "Output Files/cleaned_pcb_A.csv")

2.5 PCB Dataset B

2.5.1 Exclude sample 1348, 2016 Tag (see Max LOD for reason)

# Combine pcb2, pcb3
pcbB <- bind_rows(pcb2,pcb3)

# Dataset specific LOD
lod2 <- 
  pcbB %>% 
  filter(`Field ID`=="LOD") %>% 
  select(!`Field ID` &
         !analysis.year) %>% 
  t() %>% 
  data.frame()

# Exclude sample 1348, 2016 Tag, all header rows, all LODs
pcbB <- filter(pcbB, !`Field ID`=="1348"
               & !`Field ID`=="2016 Tag"
               & !`Field ID`=="LOD"
               & !`Field ID`=="Field ID")

2.5.2 Shorten PCB names to IUPAC number

iupac1 <- 
  as.data.frame(colnames(pcbB)) %>% 
  setNames("pcb") %>% 
  filter(!pcb == "analysis.year")

for (i in 1:nrow(iupac1)) {
  name <- iupac1[i,]
  {if(grepl('\\(',name)) {
    newname <- gsub("[\\(\\)]", "", regmatches(name, gregexpr("\\(.*?\\)", name))[[1]])
  }
    else{newname<-name}
  }
  iupac1[i,] <- newname
}
iupac1$pcb <- paste("pcb", iupac1$pcb, sep="")
colnames(pcbB)[2:22] <- iupac1$pcb[2:22]

2.5.3 Remove PCBs below LOD

# Change "ND" to 0's
for (i in 1:nrow(pcbB)) {
  for (j in 1:ncol(pcbB)) {
    {if (pcbB[i,j] == "ND") {
      pcbB[i,j] <- 0
    }
      else{next}
    }
  }
}

# Get rid of Field ID column
pcbB <- 
  pcbB %>% 
  column_to_rownames(var="Field ID")


# Get rid of PCB values below the LOD for that PCB
for (i in 1:nrow(pcbB)) {
  for (j in 1:21) {
    value <- pcbB[i,j]
    {if (value < lod2[j,1]) {
      pcbB[i,j] <- 0
    }
      else{next}
      }
  }
}

write.csv(lod2, "Output Files/lod_B.csv")

2.5.4 Create variables for analysis

# Create variables: sumpcb, pcbbin, co/noncoplanar
pcbB[,1:21] <- sapply(pcbB[,1:21], as.numeric)

pcbB$sumpcb <- rowSums(pcbB[,1:21])
pcbB$pcbbin <- ifelse(pcbB$sumpcb > 0, 1, 0)

pcbB <- pcbB %>% mutate(coplanar = pcb105 + pcb118)
pcbB <- pcbB %>% mutate(noncoplanar = sumpcb - coplanar)

pcbB$cobin <- ifelse(pcbB$coplanar > 0, 1, 0)
pcbB$noncobin <- ifelse(pcbB$noncoplanar > 0, 1, 0)

2.5.5 Write out cleaned data for future use

write.csv(x=pcbB, file = "Output Files/cleaned_pcb_B.csv")

2.6 PCB Dataset C

2.6.1 Exclude header/LOD rows

# Combine pcb4, pcb5
pcbC <- bind_rows(pcb4,pcb5)

# Dataset specific LOD
lod3 <- 
  pcbC %>% 
  filter(`Field ID`=="LOD") %>% 
  select(!`Field ID` &
         !analysis.year) %>% 
  t() %>% 
  data.frame()

# Exclude all header rows, all LODs
pcbC <- filter(pcbC, !`Field ID`=="LOD" 
               & !`Field ID`=="Field ID")

2.6.2 Shorten PCB names to IUPAC number

iupac1 <- 
  as.data.frame(colnames(pcbC)) %>% 
  setNames("pcb") %>% 
  filter(!pcb == "analysis.year")

for (i in 1:nrow(iupac1)) {
  name <- iupac1[i,]
  {if(grepl('\\(',name)) {
    newname <- gsub("[\\(\\)]", "", regmatches(name, gregexpr("\\(.*?\\)", name))[[1]])
  }
    else{newname<-name}
  }
  iupac1[i,] <- newname
}
iupac1$pcb <- paste("pcb", iupac1$pcb, sep="")
colnames(pcbC)[2:22] <- iupac1$pcb[2:22]

2.6.3 Remove PCBs below max LOD across datasets

# Change "ND" to 0's
for (i in 1:nrow(pcbC)) {
  for (j in 1:ncol(pcbC)) {
    {if (pcbC[i,j] == "ND") {
      pcbC[i,j] <- 0
    }
      else{next}
    }
  }
}

# Get rid of Field ID column
pcbC <-
  pcbC %>% 
  column_to_rownames(var="Field ID")

# Get rid of PCB values below the LOD for that PCB
for (i in 1:nrow(pcbC)) {
  for (j in 1:21) {
    value <- pcbC[i,j]
    {if (value < lod3[j,1]) {
      pcbC[i,j] <- 0
    }
      else{next}
      }
  }
}

write.csv(lod3, "Output Files/lod_C.csv")

2.6.4 Create variables for analysis

# Create variables: sumpcb, pcbbin, co/noncoplanar
pcbC[,1:21] <- sapply(pcbC[,1:21], as.numeric)

pcbC$sumpcb <- rowSums(pcbC[,1:21])
pcbC$pcbbin <- ifelse(pcbC$sumpcb > 0, 1, 0)

pcbC <- pcbC %>% mutate(coplanar = pcb105 + pcb118)
pcbC <- pcbC %>% mutate(noncoplanar = sumpcb - coplanar)

pcbC$cobin <- ifelse(pcbC$coplanar > 0, 1, 0)
pcbC$noncobin <- ifelse(pcbC$noncoplanar > 0, 1, 0)

2.6.5 Write out cleaned data for future use

write.csv(x=pcbC, file = "Output Files/cleaned_pcb_C.csv")

2.7 Whole Blood PCBs

2.7.1 Read in data and pre-process

whole <- 
  read.csv("Input Files/Contaminants/170101_PCB_whole blood.csv") %>% 
  select(2, 5:25) %>% 
  slice(c(-1:-10, -43)) 
colnames(whole) <- whole[1,]
whole <- whole[-1,]

whole_lod <- 
  whole[31,] %>% 
  t() %>% 
  data.frame() %>% 
  rownames_to_column(var="pcb") %>% 
  slice(-1) %>% 
  rename("LOD" = "X32")

whole <- 
  whole %>% 
  slice(-31)

2.7.2 Shorten PCB names to IUPAC

iupac1 <- as.data.frame(colnames(whole))
colnames(iupac1) <- "pcb"
for (i in 1:nrow(iupac1)) {
  name <- iupac1[i,]
  {if(grepl('\\(',name)) {
    newname <- gsub("[\\(\\)]", "", regmatches(name, gregexpr("\\(.*?\\)", name))[[1]])
  }
    else{newname<-name}
  }
  iupac1[i,] <- newname
}
iupac1$pcb <- paste("pcb", iupac1$pcb, sep="")
colnames(whole)[2:22] <- iupac1$pcb[2:22]

2.7.3 Change “ND” to 0’s

for (i in 1:nrow(whole)) {
  for (j in 1:ncol(whole)) {
    {if (whole[i,j] == "ND") {
      whole[i,j] <- 0
    }
      else{next}
    }
  }
}

2.7.4 Configure LOD data for BSC validation

whole_lod$pcb <- colnames(whole)[2:22]

write.csv(whole_lod, "Output Files/170101_LOD.csv")

2.7.5 Export for BSC analyses

whole <- 
  whole %>% 
  column_to_rownames(var="Field ID")

write.csv(whole, "Output Files/cleaned_pcb2_wholeblood.csv")