GeoTox Package Data • GeoTox

This package includes example data geo_tox_data. Below is a description of the data and example code for how it was gathered.

NOTE: FIPS codes can change. Since data is being pulled from various sources, ensure that the FIPS values can be used to connect data across these sources. For example, in 2022 Connecticut began the process of going from 8 legacy counties to 9 planning regions.

library(dplyr)
library(sf)
library(tidyr)
library(readr)
library(stringr)
library(purrr)
library(readxl)
library(httk)
library(httr2)

geo_tox_data <- list()

Chemical data

Exposure data

Download modeled exposure data from AirToxScreen. Results from AirToxScreen 2019 for a subset of chemicals in North Carolina counties are included in the package data as geo_tox_data$exposure.

NOTE: The 2020 release does not currently provide a single file for the exposure national concentration summaries. The 2019 data can be found following the “Previous air toxics assessments” link.

filename <- "2019_Toxics_Exposure_Concentrations.xlsx"
tmp <- tempfile(filename)
download.file(
  paste0("https://www.epa.gov/system/files/documents/2022-12/", filename),
  tmp
)
exposure <- read_xlsx(tmp)

# Normalization function
min_max_norm = function(x) {
  min_x <- min(x, na.rm = TRUE)
  max_x <- max(x, na.rm = TRUE)
  if (min_x == max_x) {
    rep(0, length(x))
  } else {
    (x - min_x) / (max_x - min_x)
  }
}

geo_tox_data$exposure <- exposure |> 
  # North Carolina counties
  filter(State == "NC", !grepl("0$", FIPS)) |> 
  # Aggregate chemicals by county
  summarize(across(-c(State:Tract), c(mean, sd)), .by = FIPS) |> 
  pivot_longer(-FIPS, names_to = "chemical") |> 
  mutate(stat = if_else(grepl("_1$", chemical), "mean", "sd"),
         chemical = gsub('.{2}$', '', chemical)) |> 
  pivot_wider(names_from = stat) |> 
  # Normalize concentrations
  mutate(norm = min_max_norm(mean), .by = chemical)

Get CASRN and PREFERRED_NAME from CompTox Dashboard

# Copy/paste the chemical names into the batch search field
# https://comptox.epa.gov/dashboard/batch-search
cat(geo_tox_data$exposure |> distinct(chemical) |> pull(), sep = "\n")
# Export results with "CAS-RN" identifiers as a csv file, then process in R

exposure_casrn <- read_csv("CCD-Batch-Search.csv",
                           show_col_types = FALSE) |> 
  filter(DTXSID != "N/A") |> 
  # Prioritize results based on FOUND_BY status
  arrange(INPUT,
          grepl("Approved Name", FOUND_BY),
          grepl("^Synonym", FOUND_BY)) |> 
  # Keep one result per INPUT
  group_by(INPUT) |> 
  slice(1) |> 
  ungroup()

# Update exposure data with CompTox Dashboard data
geo_tox_data$exposure <- geo_tox_data$exposure |> 
  inner_join(exposure_casrn, by = join_by(chemical == INPUT)) |> 
  select(FIPS, casn = CASRN, chnm = PREFERRED_NAME, mean, sd, norm)

ICE cHTS data

Use ICE cHTS data to identify active chemicals for a given set of assays.

get_cHTS_hits <- function(assays = NULL, chemids = NULL) {
  
  if (is.null(assays) & is.null(chemids)) {
    stop("Must provide at least one of 'assays' or 'chemids'")
  }
  
  # Format query parameters
  req_params <- list()
  
  if (!is.null(assays)) {
    if (!is.list(assays)) assays <- as.list(assays)
    req_params$assays <- assays
  }
  
  if (!is.null(chemids)) {
    if (!is.list(chemids)) chemids <- as.list(chemids)
    req_params$chemids <- chemids
  }
  
  # Query ICE API
  resp <- request("https://ice.ntp.niehs.nih.gov/api/v1/search") |> 
    req_body_json(req_params) |> 
    req_perform()
  
  if (resp$status_code != 200) {
    stop("Failed to retrieve data from ICE API")
  }
  
  # Return active chemicals
  result <- resp |> resp_body_json() |> pluck("endPoints")
  
  fields <- c("assay", "casrn", "dtxsid", "substanceName",
              "endpoint", "value")
  
  map(fields, \(x) map_chr(result, x)) |> 
    set_names(fields) |>
    bind_cols() |> 
    filter(endpoint == "Call", value == "Active") |> 
    select(-c(endpoint, value)) |> 
    distinct()
}

assays <- c("APR_HepG2_p53Act_1h_dn",
            "APR_HepG2_p53Act_1h_up",
            "APR_HepG2_p53Act_24h_dn",
            "APR_HepG2_p53Act_24h_up",
            "APR_HepG2_p53Act_72h_dn",
            "APR_HepG2_p53Act_72h_up",
            "ATG_p53_CIS_up",
            "TOX21_DT40",
            "TOX21_DT40_100",
            "TOX21_DT40_657",
            "TOX21_ELG1_LUC_Agonist",
            "TOX21_H2AX_HTRF_CHO_Agonist_ratio",
            "TOX21_p53_BLA_p1_ratio",
            "TOX21_p53_BLA_p2_ratio",
            "TOX21_p53_BLA_p3_ratio",
            "TOX21_p53_BLA_p4_ratio",
            "TOX21_p53_BLA_p5_ratio")

chemids <- unique(geo_tox_data$exposure$casn)
  
cHTS_hits_API <- get_cHTS_hits(assays = assays, chemids = chemids)

Dose-response data

Use the ICE API to retrieve dose-response data for selected assays and chemicals.

get_ICE_dose_resp <- function(assays = NULL, chemids = NULL) {
  
  if (is.null(assays) & is.null(chemids)) {
    stop("Must provide at least one of 'assays' or 'chemids'")
  }
  
  # Format query parameters
  req_params <- list()
  
  if (!is.null(assays)) {
    if (!is.list(assays)) assays <- as.list(assays)
    req_params$assays <- assays
  }
  
  if (!is.null(chemids)) {
    if (!is.list(chemids)) chemids <- as.list(chemids)
    req_params$chemids <- chemids
  }
  
  # Query ICE API
  resp <- request("https://ice.ntp.niehs.nih.gov/api/v1/curves") |> 
    req_body_json(req_params) |> 
    req_perform()
  
  if (resp$status_code != 200) {
    stop("Failed to retrieve data from ICE API")
  }
  
  # Return dose-response data
  result <- resp |> resp_body_json() |> pluck("curves")
  
  map(result, function(x) {
    tibble(
      endp = x[["assay"]],
      casn = x[["casrn"]],
      call = x[["dsstoxsid"]],
      chnm = x[["substance"]],
      call = x[["call"]],
      logc = map_dbl(x[["concentrationResponses"]], "concentration") |> log10(),
      resp = map_dbl(x[["concentrationResponses"]], "response")
    )
  }) |> 
    bind_rows()
}

assays <- unique(cHTS_hits_API$assay)

chemids <- intersect(cHTS_hits_API$casrn, geo_tox_data$exposure$casn)

dose_response <- get_ICE_dose_resp(assays = assays, chemids = chemids)

# Only keep active calls for assay/chemical combinations
geo_tox_data$dose_response <- dose_response |> 
  filter(call == "Active") |> 
  select(-call)

# Update dose-response data with CompTox Dashboard data
geo_tox_data$dose_response <- geo_tox_data$dose_response |>
  inner_join(exposure_casrn, by = join_by(casn == CASRN)) |>
  select(endp, casn, chnm = PREFERRED_NAME, logc, resp)

Population data

Age

Download age data from the U.S. Census Bureau by searching for “County Population by Characteristics”. A subset of data for North Carolina from 2019 is included in the package data as geo_tox_data$age.

# Data for North Carolina
url <- paste0("https://www2.census.gov/programs-surveys/popest/datasets/",
              "2010-2019/counties/asrh/cc-est2019-alldata-37.csv")
age <- read_csv(url, show_col_types = FALSE)

geo_tox_data$age <- age |>
  # 7/1/2019 population estimate
  filter(YEAR == 12) |>
  # Create FIPS
  mutate(FIPS = str_c(STATE, COUNTY)) |>
  # Keep selected columns
  select(FIPS, AGEGRP, TOT_POP)

Obesity

Follow the “Data Portal” link from CDC PLACES and search for “places county data”. Go to the desired dataset webpage, for example 2020 county data, and download the data by selecting Actions → API → Download file. A subset of data for North Carolina is included in the package data as geo_tox_data$obesity.

places <- read_csv("PLACES__County_Data__GIS_Friendly_Format___2020_release.csv",
                   show_col_types = FALSE)

# Convert confidence interval to standard deviation
extract_SD <- function(x) {
  range <- as.numeric(str_split_1(str_sub(x, 2, -2), ","))
  diff(range) / 3.92
}

geo_tox_data$obesity <- places |>
  # North Carolina Counties
  filter(StateAbbr == "NC") |>
  # Select obesity data
  select(FIPS = CountyFIPS, OBESITY_CrudePrev, OBESITY_Crude95CI) |>
  # Change confidence interval to standard deviation
  rowwise() |>
  mutate(OBESITY_SD = extract_SD(OBESITY_Crude95CI)) |>
  ungroup() |>
  select(-OBESITY_Crude95CI)

Steady-state plasma concentration (Css)

Use httk to generate Css values for combinations of age group and weight status for each chemical. The generation of these values is a time-intensive step, so one approach is to generate populations of Css values initially and then sample them later.

set.seed(2345)
n_samples <- 500

# Get CASN for which httk simulation is possible. Try using load_dawson2021,
# load_sipes2017, or load_pradeep2020 to increase availability.
load_sipes2017()
casn <- intersect(unique(geo_tox_data$dose_response$casn),
                  get_cheminfo(suppress.messages = TRUE))

# Define population demographics for httk simulation
pop_demo <- cross_join(
  tibble(age_group = list(c(0, 2), c(3, 5), c(6, 10), c(11, 15),
                          c(16, 20), c(21, 30), c(31, 40), c(41, 50),
                          c(51, 60), c(61, 70), c(71, 100))),
  tibble(weight = c("Normal", "Obese"))) |>
  # Create column of lower age_group values
  rowwise() |>
  mutate(age_min = age_group[1]) |>
  ungroup()

# Create wrapper function around httk steps
simulate_css <- function(chem.cas, agelim_years, weight_category,
                         samples, verbose = TRUE) {
  
  if (verbose) {
    cat(chem.cas,
        paste0("(", paste(agelim_years, collapse = ", "), ")"),
        weight_category,
        "\n")
  }
  
  httkpop <- list(method = "vi",
                  gendernum = NULL,
                  agelim_years = agelim_years,
                  agelim_months = NULL,
                  weight_category = weight_category,
                  reths = c(
                    "Mexican American",
                    "Other Hispanic",
                    "Non-Hispanic White",
                    "Non-Hispanic Black",
                    "Other"
                  ))
  
  css <- try(
    suppressWarnings({
      mcs <- create_mc_samples(chem.cas = chem.cas,
                               samples = samples,
                               httkpop.generate.arg.list = httkpop,
                               suppress.messages = TRUE)
      
      calc_analytic_css(chem.cas = chem.cas,
                        parameters = mcs,
                        model = "3compartmentss",
                        suppress.messages = TRUE)
    }),
    silent = TRUE
  )
  
  # Return
  if (is(css, "try-error")) {
    warning(paste0("simulate_css failed to generate data for CASN ", chem.cas))
    list(NA)
  } else {
    list(css)
  }
}

# Simulate Css values
simulated_css <- map(casn, function(chem.cas) {
  pop_demo |>
    rowwise() |>
    mutate(
      css = simulate_css(.env$chem.cas, age_group, weight, .env$n_samples)
    ) |>
    ungroup()
})
simulated_css <- setNames(simulated_css, casn)

# Remove CASN that failed simulate_css
casn_keep <- map_lgl(simulated_css, function(df) {
  !(length(df$css[[1]]) == 1 && is.na(df$css[[1]]))
})
simulated_css <- simulated_css[casn_keep]

# Get median Css values for each age_group
simulated_css <- map(
  simulated_css,
  function(cas_df) {
    cas_df |>
      nest(.by = age_group) |>
      mutate(
        age_median_css = map_dbl(data, function(df) median(unlist(df$css),
                                                           na.rm = TRUE))
      ) |>
      unnest(data)
  }
)

# Get median Css values for each weight
simulated_css <- map(
  simulated_css,
  function(cas_df) {
    cas_df |>
      nest(.by = weight) |>
      mutate(
        weight_median_css = map_dbl(data, function(df) median(unlist(df$css),
                                                              na.rm = TRUE))
      ) |>
      unnest(data) |>
      arrange(age_min, weight)
  }
)

geo_tox_data$simulated_css <- simulated_css

Prune data

Retain only those chemicals found in exposure, dose-response and Css datasets.

casn_keep <- names(geo_tox_data$simulated_css)

geo_tox_data$exposure <- geo_tox_data$exposure |> 
  filter(casn %in% casn_keep)

geo_tox_data$dose_response <- geo_tox_data$dose_response |> 
  filter(casn %in% casn_keep)

County/State boundaries

Download cartographic boundary files for counties and states from the U.S. Census Bureau. The geometry data for North Carolina counties and the state are included in the package data as geo_tox_data$boundaries.

county <- st_read("cb_2019_us_county_5m/cb_2019_us_county_5m.shp")
state <- st_read("cb_2019_us_state_5m/cb_2019_us_state_5m.shp")

geo_tox_data$boundaries <- list(
  county = county |>
    filter(STATEFP == 37) |>
    select(FIPS = GEOID, geometry),
  state = state |>
    filter(STATEFP == 37) |>
    select(geometry)
)