---
title: "BMD vs POD"
author: "Ly Ly Pham"
date: "October 10, 2018"
output: html_document
---

Data used in BMDS 

```{r setup, include=FALSE}

library(RMySQL)
library(data.table)
library(magrittr)
library(readxl)

```

```{r}

sessionInfo()

```


The U.S. EPA's toxicity database (ToxRefDB) version 2

MySQL Quarry from ToxRefdB version 2_0

```{r eval=FALSE}

con <- dbConnect(drv = RMySQL::MySQL(),  group="toxrefdb_2_0")

toxRefv2 <- dbGetQuery(con, "SELECT 
          chemical.chemical_id,
          chemical.dsstox_gsid,
          chemical.casrn,
          chemical.preferred_name,
          study.study_id,
          study.chemical_id,
          study.dsstox_gsid,
          study.processed,
          study.batch_name,
          study.data_entry_status,
          study.data_entry_level,
          study.data_usability,
          study.study_type,
          study.study_year,
          study.study_source,
          study.species,
          study.strain_group,
          study.admin_route,
          study.admin_method,
          study.substance_purity,
          study.dose_start,
          study.dose_start_unit,
          study.dose_end,
          study.dose_end_unit,
          endpoint.endpoint_category,
          endpoint.endpoint_type,
          endpoint.endpoint_target,
          endpoint.endpoint_id,
          tg_effect.life_stage,
          tg_effect.tg_effect_id,
          effect.effect_id,
          effect.effect_desc,
          tg.sex,
          tg.generation,
          dose.dose_level,
          dtg.dose_adjusted,
          dtg.dose_adjusted_unit,
          dtg_effect.treatment_related,
          dtg_effect.critical_effect,
          tested_status,
          reported_status
      FROM 
          (((((((((chemical INNER JOIN study ON chemical.chemical_id=study.chemical_id)
              LEFT JOIN dose ON dose.study_id=study.study_id)
              LEFT JOIN tg ON tg.study_id=study.study_id)
                  LEFT JOIN dtg ON tg.tg_id=dtg.tg_id AND dose.dose_id=dtg.dose_id)
                      LEFT JOIN tg_effect ON tg.tg_id=tg_effect.tg_id)
                          LEFT JOIN dtg_effect ON tg_effect.tg_effect_id=dtg_effect.tg_effect_id AND dtg.dtg_id=dtg_effect.dtg_id)
                              LEFT JOIN effect ON effect.effect_id=tg_effect.effect_id)
                                  LEFT JOIN endpoint ON endpoint.endpoint_id=effect.endpoint_id)
                                   LEFT JOIN obs ON obs.study_id=study.study_id AND obs.endpoint_id=endpoint.endpoint_id)") %>% 
  data.table() 

bmd_table<- dbGetQuery(con, "SELECT * FROM bmd_models") %>% data.table()
```

```{r include=FALSE}

load("C:/Users/lpham/Desktop/Projects/ToxRefDB/BMDvsPOD/toxRefv2.RData")
load("C:/Users/lpham/Desktop/Projects/ToxRefDB/BMDvsPOD/bmd_table.RData")

```


```{r warning=FALSE}

toxRefv2[ , dataset_id:= paste(study_id, endpoint_id, tg_effect_id, sex, sep = "_")]

toxRefv2.2 <- toxRefv2[ !is.na(endpoint_id)]

## Convert NA to zero for fields: treatment_related and critical_effect 
toxRefv2.2[ is.na(treatment_related), treatment_related:=0]
toxRefv2.2[ is.na(critical_effect), critical_effect:=0]

# LEL at the Effect level
toxRefv2.2[ , tr := max(treatment_related, na.rm = TRUE), by = c("dataset_id") ]
toxRefv2.2[ tr == 1, effect.lel := min(dose_adjusted[treatment_related == 1]), by = c("dataset_id") ]

# LOAEL at the effect level
toxRefv2.2[ , ce:= max(critical_effect, na.rm = TRUE), by = c("dataset_id")]
toxRefv2.2[ ce ==1 , effect.loael:= min(dose_adjusted[critical_effect== 1]), by = c("dataset_id") ]
toxRefv2.2[ is.na(effect.loael), effect.loael:=effect.lel] #when no loael is id, lel is used


# NEL at the effect level
toxRefv2.2[ tr == 0 , effect.nel:=max(dose_adjusted[treatment_related==0]), by = c("dataset_id") ]
toxRefv2.2[ tr == 1 , effect.nel:=max(dose_adjusted[treatment_related==0]), by = c("dataset_id")]


# NOAEL at the effect level
#toxRefv2.2[ , ce.dose:= max(critical_effect), by = c("study_id", "endpoint_id", "tg_effect_id", "dose_level")]
# if there was a critical effect, at what dose level 
toxRefv2.2[ critical_effect==1 , loael.dose.level:=min(dose_level), by = c("dataset_id")]

#tr2[ ce.dose == 0 , loael.dose.level:=0, by = c("study_id", "generation")]
#tr2[ is.na(loael.dose.level), loael.dose.level:=0]
toxRefv2.2[ , loael.dose.level:=max(loael.dose.level), by = c("dataset_id")]

# for critical effects seens at one dose level but not at a higher dose level
toxRefv2.2[ dose_level>=loael.dose.level & loael.dose.level>0, ce.dose.correct:=1]
toxRefv2.2[ is.na(ce.dose.correct), ce.dose.correct:=0]

toxRefv2.2[ce.dose.correct==0 , effect.noael:=max(dose_adjusted), by = c("dataset_id")]
toxRefv2.2[is.na(effect.noael), effect.noael:=0, by = c("dataset_id")]
toxRefv2.2[ , effect.noael:=max(effect.noael), by = c("dataset_id")]

```


Get only the data from bmd_table

```{r}

trv2 <- toxRefv2.2[ ,list(dataset_id, effect.nel, effect.noael, effect.lel, effect.loael)] %>% unique()

bmd.rec <- bmd_table[ recommended=="1"]

```


```{r}

tr.bmd <- merge(bmd.rec, trv2, by="dataset_id")

```

Comparing BMDL with NOAEL 

* N = NOAEL/NEL > BMDL 
* B = NOAEL/NEL < BMDL 
* E = NOAEL/NEL = BMDL 

```{r}

tr.bmd[ effect.noael > BMDL, noael.bmdl:="N"]
tr.bmd[ BMDL > effect.noael, noael.bmdl:="B"]
tr.bmd[ effect.noael == BMDL, noael.bmdl:="E"]

tr.bmd[ effect.nel > BMDL, nel.bmdl:="N"]
tr.bmd[ BMDL > effect.nel, nel.bmdl:="B"]
tr.bmd[ effect.nel == BMDL, nel.bmdl:="E"]

tr.bmd5 <- tr.bmd[ bmr %in% c(1, 5)]
tr.bmd10 <- tr.bmd[ bmr %in% c(1, 10)]

```

```{r}

table(tr.bmd5$nel.bmdl)

table(tr.bmd5$noael.bmdl)

table(tr.bmd10$nel.bmdl)

table(tr.bmd10$noael.bmdl)


```


Comparing BMD with LEL 

* L = LEL/LOAEL > BMD 
* B = LEL/LOAEL < BMD 
* E = LEL/LOAEL = BMD 

```{r}

tr.bmd[ effect.lel > BMD, lel.bmd:="L"]
tr.bmd[ BMD > effect.lel, lel.bmd:="B"]
tr.bmd[ effect.lel == BMD, lel.bmd:="E"]

tr.bmd[ effect.loael > BMD, loael.bmd:="L"]
tr.bmd[ BMD > effect.loael, loael.bmd:="B"]
tr.bmd[ effect.loael == BMD, loael.bmd:="E"]

tr.bmd5 <- tr.bmd[ bmr %in% c(1, 5)]
tr.bmd10 <- tr.bmd[ bmr %in% c(1, 10)]

```

```{r}

table(tr.bmd5$lel.bmd)

table(tr.bmd5$loael.bmd)

table(tr.bmd10$lel.bmd)

table(tr.bmd10$loael.bmd)


```















