Nothing
inst/doc/Vi_Kapraun2022.R
In httk: High-Throughput Toxicokinetics
## ----configure_knitr, eval = TRUE---------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
options(rmarkdown.html_vignette.check_title = FALSE)
## ----clear_memory, eval = TRUE------------------------------------------------
rm(list=ls())
## ----runchunks, eval = TRUE---------------------------------------------------
# Set whether or not the following chunks will be executed (run):
execute.vignette <- FALSE
## ----load_packages, eval = execute.vignette-----------------------------------
# #
# #
# # Setup
# #
# #
# #library(readxl)
# library(ggplot2)
# library(httk)
# library(scales)
# library(gridExtra)
# library(cowplot)
## ----load_code, eval = execute.vignette---------------------------------------
# scientific_10 <- function(x) {
# out <- gsub("1e", "10^", scientific_format()(x))
# out <- gsub("\\+","",out)
# out <- gsub("10\\^01","10",out)
# out <- parse(text=gsub("10\\^00","1",out))
# }
#
# RMSE <- function(x)
# {
# mean(x$residuals^2)^(1/2)
# }
## ----chem_numbers, eval=FALSE-------------------------------------------------
# length(get_cheminfo(model="fetal_pbtk",suppress.messages=TRUE))
## ----load_aylward_data, eval = execute.vignette-------------------------------
# #MFdata <- read_excel("Aylward-MatFet.xlsx")
# MFdata <- httk::aylward2014
#
# cat(paste("summarized data from over 100 studies covering ",
# length(unique(MFdata$DTXSID)[!(unique(MFdata$DTXSID)%in%c("","-"))]),
# " unique chemicals structures\n",sep=""))
#
# # We don't want to exclude the volatiles and PFAS at this point:
# MFdata.httk <- subset(MFdata,DTXSID %in% get_cheminfo(
# info="DTXSID",
# suppress.messages=TRUE))
# MFdata.httk[MFdata.httk$Chemical.Category=="bromodiphenylethers",
# "Chemical.Category"] <- "Bromodiphenylethers"
# MFdata.httk[MFdata.httk$Chemical.Category=="organochlorine Pesticides",
# "Chemical.Category"] <- "Organochlorine Pesticides"
# MFdata.httk[MFdata.httk$Chemical.Category=="polyaromatic hydrocarbons",
# "Chemical.Category"] <- "Polyaromatic Hydrocarbons"
#
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
# "infant.maternal.conc...Central.tendency..calculate.j.k..or.report.paired.result."] <-
# "MFratio"
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
# "PREFERRED_NAME"] <-
# "Chemical"
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
# "details.on.matrix.comparison...e.g...cord.blood.lipid..maternal.serum.lipid..or.cord.blood.wet.weight..maternal.whole.blood.wet.weight"] <-
# "Matrix"
#
# # Format the columns:
# MFdata.httk$MFratio <- as.numeric(MFdata.httk$MFratio)
# MFdata.httk$Chemical <- as.factor(MFdata.httk$Chemical)
# MFdata.httk$Matrix <- as.factor(MFdata.httk$Chemical)
# MFdata.httk$Chemical.Category <- as.factor(MFdata.httk$Chemical.Category)
#
# colnames(MFdata.httk)[15] <- "infant"
# colnames(MFdata.httk)[16] <- "maternal"
# colnames(MFdata.httk)[17] <- "obs.units"
#
# MFdata.httk$infant <- suppressWarnings(as.numeric(MFdata.httk$infant))
# MFdata.httk$maternal <- suppressWarnings(as.numeric(MFdata.httk$maternal))
# MFdata.httk$AVERAGE_MASS <- as.numeric(MFdata.httk$AVERAGE_MASS)
## ----process_aylward_data, eval = execute.vignette----------------------------
# convert1.units <- c("ng/ml","ng/mL","ug/L","ug/l","ng/mL serum","ng/g",
# "ng/g wet wt.","ppb")
#
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] / # ng/ml = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] / # ng/ml = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"obs.units"] <- "uM"
#
# convert2.units <- c("mg/L","ppm")
#
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] * 1000 / # mg/L = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"]* 1000 / # mg/L = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"obs.units"] <- "uM"
## ----make_aylward_preds, eval = execute.vignette------------------------------
# # Simulate starting from the 13th week of gestation until full term (40 weeks):
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
#
# # For each chemical with maternal-fetal ratio data and HTTK in vitro data:
# for (this.id in unique(MFdata.httk$DTXSID))
# {
# print(this.id)
# p <- parameterize_steadystate(dtxsid=this.id,
# suppress.messages=TRUE)
# # skip chemicals where Fup is 0:
# if (p$Funbound.plasma>1e-4)
# {
#
# # Load the chemical-specifc paramaters:
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =TRUE,
# suppress.messages=TRUE)
# # Skip chemicals where the 95% credible interval for Fup spans more than 0.1 to
# # 0.9 (that is, Fup is effectively unknown):
# if (!is.na(p$Funbound.plasma.dist))
# {
# if (as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][3])>0.9 &
# as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][2])<0.11)
# {
# skip <- TRUE
# } else skip <- FALSE
# } else skip <- FALSE
#
# if (!skip)
# {
# # Solve the PBTK equations for the full simulation time assuming 1 total daily
# # dose of 1 mg/kg/day spread out over three evenly spaces exposures:
# out <- solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# suppress.messages=TRUE)
# # Identify the concentrations from the final (279th) day:
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# # Average over the final day:
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred"] <- mean(out[last.row,"Cplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred"] <- mean(out[last.row,"Cfplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred"] <-
# mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
#
# # Repeat the analysis without the adjustment to fetal Fup:
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =FALSE,
# suppress.messages = TRUE)
# out <- solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# maxsteps=1e7,
# suppress.messages = TRUE)
#
# last.row <- which(out[,"time"]>279) # The whole final day
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred.nofup"] <- mean(out[last.row,"Cplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred.nofup"] <- mean(out[last.row,"Cfplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred.nofup"] <-
# mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# }
# }
# }
#
# # Something is wrong with cotinine:
# MFdata.httk <- subset(MFdata.httk,Chemical!="Cotinine")
#
# max.chem <- MFdata.httk[which(MFdata.httk$MFratio==max(MFdata.httk$MFratio)),]
# min.chem <- MFdata.httk[which(MFdata.httk$MFratio==min(MFdata.httk$MFratio)),]
# cat(paste("The minimum observed ratio was ",
# signif(min.chem[,"MFratio"],2),
# " for ",
# min.chem[,"Chemical"],
# " and the maximum was ",
# signif(max.chem[,"MFratio"],2),
# " for ",
# max.chem[,"Chemical"],
# ".\n",sep=""))
#
#
## ----cleanup_repeated_aylward_figure, eval = execute.vignette-----------------
# MFdata.main <- NULL
# MFdata.outliers <- NULL
# for (this.id in unique(MFdata.httk$DTXSID))
# {
# this.subset <- subset(MFdata.httk,DTXSID==this.id)
# this.row <- this.subset[1,]
# this.row$N.obs <- dim(this.subset)[1]
# this.row$MFratio <- median(this.subset$MFratio)
# this.row$MFratio.Q25 <- quantile(this.subset$MFratio,0.25)
# this.row$MFratio.Q75 <- quantile(this.subset$MFratio,0.75)
# MFdata.main <- rbind(MFdata.main,this.row)
# this.subset <- subset(this.subset,
# MFratio<this.row$MFratio.Q25 |
# MFratio>this.row$MFratio.Q75)
# MFdata.outliers <- rbind(MFdata.outliers,this.subset)
# }
## ----aylward_figure_no_fup, eval = execute.vignette---------------------------
# FigCa <- ggplot(data=MFdata.main) +
# geom_segment(color="grey",aes(
# x=MFratio.pred.nofup,
# y=MFratio.Q25,
# xend=MFratio.pred.nofup,
# yend=MFratio.Q75))+
# geom_point(aes(
# x=MFratio.pred.nofup,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=4) +
# scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
# geom_point(data=MFdata.outliers,aes(
# x=MFratio.pred.nofup,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=2) +
# xlim(0,2) +
# ylim(0,3) +
# # geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
# # scale_y_log10(label=scientific_10) +
# #,limits=c(10^-7,100)) +
# # scale_x_log10(label=scientific_10) +
# # ,limits=c(10^-7,100)) +
# # annotation_logticks() +
# geom_abline(slope=1, intercept=0) +
# # geom_abline(slope=1, intercept=1,linetype="dashed") +
# # geom_abline(slope=1, intercept=-1,linetype="dashed") +
# ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) +
# xlab("Generic PBTK Predicted Ratio") +
# # scale_color_brewer(palette="Set2") +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=14))+
# theme(legend.text=element_text(size=10))+
# guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+
# guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))
#
# print(FigCa)
## ----aylward_analysis_text, eval = execute.vignette---------------------------
#
# cat(paste("In Figure 4 we compare predictions made with our high-throughput \
# human gestational PBTK model with experimental observations on a per chemical \
# basis wherever we had both in vitro HTTK data and in vivo observations (",
# length(unique(MFdata.main$DTXSID)),
# " chemicals).\n",sep=""))
#
#
# repeats <- subset(MFdata.main,N.obs>1)
#
# cat(paste("Multiple observations were available for ",
# dim(repeats)[1],
# " of the chemicals,\n",sep=""))
#
#
# max.chem <- as.data.frame(repeats[which(repeats$MFratio==max(repeats$MFratio)),])
# min.chem <- as.data.frame(repeats[which(repeats$MFratio==min(repeats$MFratio)),])
#
# cat(paste("However, among the chemicals with repeated observations, the median \
# observations only ranged from ",
# signif(min.chem[,"MFratio"],2),
# " for ",
# min.chem[,"Chemical"],
# " to ",
# signif(max.chem[,"MFratio"],2),
# " for ",
# max.chem[,"Chemical"],
# ".\n",sep=""))
#
# max.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==max(MFdata.main$MFratio.pred,na.rm=TRUE)),])
# min.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==min(MFdata.main$MFratio.pred,na.rm=TRUE)),])
#
# cat(paste("The predictions for all chemicals ranged from ",
# signif(min.chem[,"MFratio.pred"],2),
# " for ",
# min.chem[,"Chemical"],
# " to ",
# signif(max.chem[,"MFratio.pred"],2),
# " for ",
# max.chem[,"Chemical"],
# ".\n",sep=""))
#
#
#
# fit1 <- lm(data=MFdata.main,MFratio~MFratio.pred.nofup)
# summary(fit1)
# RMSE(fit1)
#
# fit1sub <- lm(data=subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred.nofup)
# summary(fit1sub)
# RMSE(fit1sub)
#
## ----aylward_figure_with_fup, eval = execute.vignette-------------------------
# FigCb <- ggplot(data=MFdata.main) +
# geom_segment(color="grey",aes(
# x=MFratio.pred,
# y=MFratio.Q25,
# xend=MFratio.pred,
# yend=MFratio.Q75))+
# geom_point(aes(
# x=MFratio.pred,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=3) +
# scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
# geom_point(data=MFdata.outliers,aes(
# x=MFratio.pred,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=1) +
# xlim(0,2) +
# ylim(0,3) +
# # geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
# # scale_y_log10(label=scientific_10) +
# #,limits=c(10^-7,100)) +
# # scale_x_log10(label=scientific_10) +
# # ,limits=c(10^-7,100)) +
# # annotation_logticks() +
# geom_abline(slope=1, intercept=0) +
# # geom_abline(slope=1, intercept=1,linetype="dashed") +
# # geom_abline(slope=1, intercept=-1,linetype="dashed") +
# ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) +
# xlab(expression(paste(italic("In vitro")," Predicted Ratio"))) +
# # scale_color_brewer(palette="Set2") +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=14))+
# theme(legend.text=element_text(size=10))+
# guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+
# guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))
#
# print(FigCb)
#
## ----voc_analysis, eval = execute.vignette------------------------------------
# # Mean logHenry's law constant for PAH's:
# mean(subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFdata.main,Chemical.Category=="Polyaromatic Hydrocarbons")$DTXSID)$logHenry)
#
# nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$DTXSID
# fluoros <- chem.physical_and_invitro.data$DTXSID[regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))!=-1]
#
# fit2 <- lm(data=MFdata.main,MFratio~MFratio.pred)
# summary(fit2)
# RMSE(fit2)
#
# fit2sub <- lm(data=subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred)
# summary(fit2sub)
#
# RMSE(fit2sub)
#
# cat(paste("When volatile and fluorinated chemicals are omitted only ",
# dim(subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))))[1],
# " evaluation chemicals remain, but the R2 is ",
# signif(summary(fit1sub)$adj.r.squared,2),
# " and the RMSE is ",
# signif(RMSE(fit1sub),2),
# " without the correction. When the fetal plasma fraction unbound correction is used, the predictivity decreases, slightly: R2 is ",
# signif(summary(fit2sub)$adj.r.squared,2),
# " and the RMSE is ",
# signif(RMSE(fit2sub),2),
# " for the non-volatile, non-fluorinated chemicals.\n",
# sep=""))
#
#
#
# cat(paste("We compare the RMSE for our predictions to the standard deviation \
# of the observations ",
# signif(sd(MFdata.main$MFratio)[1],2),
# " (",
# signif(sd(subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))$MFratio),2),
# " for non PAH or fluorinated compounds).\n",sep=""))
#
# cat(paste("The average standard deviation for chemicals with repeated observations was ",
# signif(sd(subset(MFdata.main,N.obs>1)$MFratio)[1],2),
# " (",
# signif(sd(subset(MFdata.main,
# N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))$MFratio),2),
# " for non PAH or fluorinated compounds).\n",sep=""))
#
# fit3 <- lm(data=repeats,MFratio~MFratio.pred.nofup)
# summary(fit3)
# RMSE(fit3)
#
# fit3sub <- lm(data=subset(MFdata.main, N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred.nofup)
# summary(fit3sub)
#
#
# fit4 <- lm(data=subset(MFdata.main,N.obs > 1),MFratio~MFratio.pred)
# summary(fit4)
# RMSE(fit4)
#
# fit4sub <- lm(data=subset(MFdata.main, N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred)
# summary(fit4sub)
#
# cat(paste("Overall, we observed relatively poor correlation (R2 = ",
# signif(summary(fit1)$adj.r.squared,2),
# ", RMSE = ",
# signif(RMSE(fit1),2),
# ") without our fetal fup correction that was unchanged with that assumption (R2 = ",
# signif(summary(fit2)$adj.r.squared,2),
# ", RMSE = ",
# signif(RMSE(fit2),2),
# ").\n",sep=""))
#
# repeats <-subset(MFdata.main,N.obs > 1)
# cat(paste("The RMSE of the predictions for the ",
# dim(subset(repeats,!(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))))[1],
# " non-PAH and non-fluorinated compounds with repeated observations is ",
# signif(RMSE(fit4sub),2),
# " with the fup correction and ",
# signif(RMSE(fit3sub),2),
# " without.\n",sep=""))
#
# cat(paste("These values are close to the standard deviation of the mean but the p-value for the chemicals with repeated observations is ",
# signif(pf(
# summary(fit4sub)$fstatistic[1],
# summary(fit4sub)$fstatistic[2],
# summary(fit4sub)$fstatistic[3]),2),
# " indicating some value for the predictive model, albeit for only ",
# dim(subset(repeats,!(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))))[1],
# " chemicals",sep=""))
#
#
#
# Fig4.table <- data.frame()
# Fig4.table["Number of Chemicals","All Fig 4"] <- length(unique(MFdata.httk$DTXSID))
# Fig4.table["Number of Observations","All Fig 4"] <- dim(MFdata.httk)[1]
# Fig4.table["Observed Mean (Min - Max)","All Fig 4"] <- paste(
# signif(mean(MFdata.httk$MFratio),2)," (",
# signif(min(MFdata.httk$MFratio),2)," - ",
# signif(max(MFdata.httk$MFratio),2),")",sep="")
# Fig4.table["Observed Standard Deviation","All Fig 4"] <- signif(sd(MFdata.httk$MFratio),2)
# Fig4.table["Predicted Mean (Min - Max)","All Fig 4"] <- paste(
# signif(mean(MFdata.main$MFratio.pred,na.rm=TRUE),2)," (",
# signif(min(MFdata.main$MFratio.pred,na.rm=TRUE),2)," - ",
# signif(max(MFdata.main$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","All Fig 4"] <- signif(RMSE(fit2),2)
# Fig4.table["R-squared","All Fig 4"] <- signif(summary(fit2)[["adj.r.squared"]],2)
# Fig4.table["p-value","All Fig 4"] <- signif(summary(fit2)[["coefficients"]]["MFratio.pred",4],2)
#
# MFdata.sub1 <- subset(MFdata.httk,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))
#
# Fig4.table["Number of Chemicals","No PFAS/PAH"] <- length(unique(MFdata.sub1$DTXSID))
# Fig4.table["Number of Observations","No PFAS/PAH"] <- dim(MFdata.sub1)[1]
# Fig4.table["Observed Mean (Min - Max)","No PFAS/PAH"] <- paste(
# signif(mean(MFdata.sub1$MFratio,na.rm=TRUE),2)," (",
# signif(min(MFdata.sub1$MFratio,na.rm=TRUE),2)," - ",
# signif(max(MFdata.sub1$MFratio,na.rm=TRUE),2),")",sep="")
# Fig4.table["Observed Standard Deviation","No PFAS/PAH"] <- signif(sd(MFdata.sub1$MFratio),2)
#
# MFdata.sub2 <- subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))
#
# Fig4.table["Predicted Mean (Min - Max)","No PFAS/PAH"] <- paste(
# signif(mean(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," (",
# signif(min(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," - ",
# signif(max(MFdata.sub2$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","No PFAS/PAH"] <- signif(RMSE(fit2sub),2)
# Fig4.table["R-squared","No PFAS/PAH"] <- signif(summary(fit2sub)[["adj.r.squared"]],2)
# Fig4.table["p-value","No PFAS/PAH"] <- signif(summary(fit2sub)[["coefficients"]]["MFratio.pred",4],2)
#
#
# MFdata.sub3 <- subset(MFdata.main,N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))
#
# Fig4.table["Number of Chemicals","Replicates Only"] <- length(unique(MFdata.sub3$DTXSID))
# Fig4.table["Number of Observations","Replicates Only"] <- dim(MFdata.sub3)[1]
# Fig4.table["Observed Mean (Min - Max)","Replicates Only"] <- paste(
# signif(mean(MFdata.sub3$MFratio),2)," (",
# signif(min(MFdata.sub3$MFratio),2)," - ",
# signif(max(MFdata.sub3$MFratio),2),")",sep="")
# Fig4.table["Observed Standard Deviation","Replicates Only"] <- signif(sd(MFdata.sub3$MFratio),2)
#
# Fig4.table["Predicted Mean (Min - Max)","Replicates Only"] <- paste(
# signif(mean(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," (",
# signif(min(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," - ",
# signif(max(MFdata.sub3$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","Replicates Only"] <- signif(RMSE(fit4sub),2)
# Fig4.table["R-squared","Replicates Only"] <- signif(summary(fit4sub)[["adj.r.squared"]],2)
# Fig4.table["p-value","Replicates Only"] <- signif(summary(fit4sub)[["coefficients"]]["MFratio.pred",4],2)
#
# write.csv(Fig4.table[1:6,],file="Fig4Table.csv")
## ----dallmann2018_data, eval = execute.vignette-------------------------------
# #TKstats <- as.data.frame(read_excel("MaternalFetalAUCData.xlsx"))
# TKstats <- httk::pregnonpregaucs
#
#
# ids <- unique(TKstats$DTXSID)
# cat(paste(sum(ids %in% get_cheminfo(
# model="fetal_pbtk",
# info="dtxsid",
# suppress.messages=TRUE)),
# "chemicals for which the model fetal_pbtk can run that are in the Dallmann data set."))
#
#
# TKstats.Cmax <- subset(TKstats,DTXSID!="" & Parameter=="Cmax")
# TKstats <- subset(TKstats,DTXSID!="" & Parameter %in% c("AUCinf","AUClast"))
#
# ids <- unique(TKstats$DTXSID)
# cat(paste(sum(ids %in% get_cheminfo(
# model="fetal_pbtk",
# info="dtxsid",
# suppress.messages=TRUE)),
# "chemicals for which the model fetal_pbtk can run that have AUCs in the Dallmann data set."))
#
#
# # Assumed body weight (kg) from Kapraun 2019
# BW.nonpreg <- 61.103
#
# #TKstats$Dose <- TKstats$Dose/BW
# #TKstats$Dose.Units <- "mg/kg"
# colnames(TKstats)[colnames(TKstats)=="Observed Pregnant"] <- "Observed.Pregnant"
# colnames(TKstats)[colnames(TKstats)=="Observed Pregnant5"] <- "Observed.Pregnant5"
# colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant"] <- "Observed.Non.Pregnant"
# colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant2"] <- "Observed.Non.Pregnant2"
# colnames(TKstats)[colnames(TKstats)=="Dose Route"] <- "Dose.Route"
## ----dallmann2018_preds, eval = execute.vignette------------------------------
# for (this.id in unique(TKstats$DTXSID))
# {
# if (any(regexpr("ng",TKstats[TKstats$DTXSID==this.id,"Dose Units"])!=-1))
# {
# }
# if (this.id %in% get_cheminfo(
# info="DTXSID",
# model="fetal_pbtk",
# suppress.messages=TRUE))
# {
# this.subset <- subset(TKstats,DTXSID==this.id)
# p <- parameterize_pbtk(
# dtxsid=this.id,
# suppress.messages=TRUE)
# p$hematocrit <- 0.39412 # Kapraun 2019 (unitless)
# p$Rblood2plasma <- calc_rblood2plasma(
# parameters=p,
# suppress.messages=TRUE)
# p$BW <- BW.nonpreg
# p$Qcardiacc <- 301.78 / p$BW^(3/4) # Kapraun 2019 (L/h/kg^3/4)
# for (this.route in unique(this.subset$Dose.Route))
# {
# this.route.subset <- subset(this.subset, Dose.Route==this.route)
# if (this.route.subset[1,"Gestational.Age.Weeks"] > 12)
# {
# this.dose <- this.route.subset$Dose
# # Non-pregnant PBPK:
# out.nonpreg <- solve_pbtk(
# parameters=p,
# times = seq(0, this.route.subset[1,"NonPreg.Duration.Days"],
# this.route.subset[1,"NonPreg.Duration.Days"]/100),
# dose=this.dose/BW.nonpreg, # mg/kg
# daily.dose=NULL,
# iv.dose=(this.route=="iv"),
# output.units="uM",
# suppress.messages=TRUE)
# # Pregnant PBPK:
# BW.preg <- suppressWarnings(calc_maternal_bw(
# week=this.route.subset[1,"Gestational.Age.Weeks"]))
# out.preg <- solve_fetal_pbtk(
# dtxsid=this.id,
# times = seq(
# this.route.subset[1,"Gestational.Age.Weeks"]*7,
# this.route.subset[1,"Gestational.Age.Weeks"]*7 +
# this.route.subset[1,"Preg.Duration.Days"],
# this.route.subset[1,"Preg.Duration.Days"]/100),
# dose=this.dose/BW.preg, # mg/kg
# iv.dose=(this.route=="iv"),
# daily.dose=NULL,
# output.units="uM",
# suppress.messages=TRUE)
# if (any(regexpr("AUC",this.subset$Parameter)!=-1))
# {
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("AUC",TKstats$Parameter)!=-1,
# "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"AUC"])*24 #uMdays->uMh
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("AUC",TKstats$Parameter)!=-1,
# "Predicted.Pregnant.httk"] <- max(out.preg[,"AUC"])*24 # uM days -> uM h
# }
# if (any(regexpr("Cmax",this.subset$Parameter)!=-1))
# {
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("Cmax",TKstats$Parameter)!=-1,
# "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"Cplasma"])
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("Cmax",TKstats$Parameter)!=-1,
# "Predicted.Pregnant.httk"] <- max(out.preg[,"Cfplasma"])
# }
# }
# }
# }
# }
#
# TKstats$Ratio.obs <- TKstats$Observed.Pregnant / TKstats$Observed.Non.Pregnant
# TKstats$Ratio.httk <- TKstats$Predicted.Pregnant.httk / TKstats$Predicted.Non.Pregnant.httk
# TKstats <- subset(TKstats, !is.na(TKstats$Ratio.httk))
#
# write.csv(TKstats,file="Table-Dallmann2018.csv",row.names=FALSE)
## ----dallmann2018_figure, eval = execute.vignette-----------------------------
# FigEa <- ggplot(data=TKstats) +
# geom_point(aes(
# y=Observed.Non.Pregnant2,
# x=Predicted.Non.Pregnant.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# scale_shape_manual(values=c(22,24))+
# xlab("httk Predicted (uM*h)") +
# ylab("Observed") +
# scale_x_log10(#limits=c(10^-3,10^3),
# label=scientific_10) +
# scale_y_log10(#limits=c(10^-3,10^3),
# label=scientific_10) +
# ggtitle("Non-Pregnant") +
# theme_bw() +
# theme(legend.position="none")+
# theme( text = element_text(size=12)) +
# annotate("text", x = 0.1, y = 1000, label = "A", fontface =2)
#
# #print(FigEa)
# cat(paste(
# "For the Dallman et al. (2018) data we observe an average-fold error (AFE)\n\ of",
# signif(mean(TKstats$Predicted.Non.Pregnant.httk/TKstats$Observed.Non.Pregnant2),2),
# "and a RMSE (log10-scale) of",
# signif((mean((log10(TKstats$Predicted.Non.Pregnant.httk)-log10(TKstats$Observed.Non.Pregnant2))^2))^(1/2),2),
# "for non-pregnant woman.\n"))
#
# FigEb <- ggplot(data=TKstats) +
# geom_point(aes(
# y=Observed.Pregnant5,
# x=Predicted.Pregnant.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# scale_shape_manual(values=c(22,24))+
# xlab("httk Predicted (uM*h)") +
# ylab("Observed") +
# scale_x_log10(#limits=c(10^-5,10^3),
# label=scientific_10) +
# scale_y_log10(#limits=c(10^-5,10^3),
# label=scientific_10) +
# ggtitle("Pregnant")+
# theme_bw() +
# theme(legend.position="none")+
# theme( text = element_text(size=12)) +
# annotate("text", x = 0.1, y = 300, label = "B", fontface =2)
#
# #print(FigEb)
# cat(paste(
# "We observe an AFE of",
# signif(mean(TKstats$Predicted.Pregnant.httk/TKstats$Observed.Pregnant5),2),
# "and a RMSE (log10-scale) of",
# signif((mean((log10(TKstats$Predicted.Pregnant.httk)-log10(TKstats$Observed.Pregnant5))^2))^(1/2),2),
# "for pregnant woman.\n"))
#
#
#
# FigEc <- ggplot(data=TKstats) +
# geom_point(aes(
# x=Predicted.Non.Pregnant.httk,
# y=Predicted.Pregnant.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# scale_shape_manual(values=c(22,24),name="Route")+
# ylab("httk Pregnant (uM*h)") +
# xlab("httk Non-Pregnant (uM*h)") +
# scale_x_log10(#limits=c(10^-1,10^2),
# label=scientific_10) +
# scale_y_log10(#limits=c(10^-1,10^2),
# label=scientific_10) +
# ggtitle("Model Comparison") +
# theme_bw() +
# theme(legend.position="left")+
# theme( text = element_text(size=12))+
# theme(legend.text=element_text(size=10))+
# guides(fill=guide_legend(ncol=2,override.aes=list(shape=21)),alpha=guide_legend(ncol=2),shape=guide_legend(ncol=2))
# print(FigEc)
# FigEc <- get_legend(FigEc)
#
#
# FigEd <- ggplot(data=subset(TKstats,!is.na(Ratio.httk))) +
# geom_point(aes(
# y=Ratio.obs,
# x=Ratio.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# scale_shape_manual(values=c(22,24))+
# xlab("httk Predicted") +
# ylab("Observed") +
# scale_x_continuous(limits=c(0.25,3)) +
# scale_y_continuous(limits=c(0.25,3)) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# ggtitle("Ratio") +
# theme_bw() +
# theme(legend.position="none")+
# theme( text = element_text(size=12)) +
# annotate("text", x = 0.5, y = 2.75, label = "C", fontface =2)
#
# dev.new()
# grid.arrange(FigEa,FigEb,FigEc,FigEd,nrow=2)
#
# write.csv(subset(TKstats,
# !is.na(Ratio.httk)),
# file="DallmanTable.txt")
#
## ----curley1969_data, eval = execute.vignette---------------------------------
# Curley <- as.data.frame(read_excel("Curley1969.xlsx"))
# dim(Curley)
# Curley.compounds <- Curley[1,4:13]
# Curley <- Curley[4:47,]
# colnames(Curley)[1] <- "Tissue"
# colnames(Curley)[2] <- "N"
# colnames(Curley)[3] <- "Stat"
## ----curley1969_calcpcs, eval = execute.vignette------------------------------
# Curley.pcs <- NULL
# cord.blood <- subset(Curley, Tissue == "Cord Blood")
# suppressWarnings(
# for (this.tissue in unique(Curley$Tissue))
# if (this.tissue != "Cord Blood")
# {
# this.subset <- subset(Curley, Tissue == this.tissue)
# for (this.chemical in colnames(Curley)[4:13])
# {
# if (!is.na(as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical])) &
# !is.na(as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])))
# {
# this.row <- data.frame(
# Compound = Curley.compounds[,this.chemical],
# DTXSID = this.chemical,
# Tissue = this.tissue,
# PC = as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical]) /
# as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
# )
# Curley.pcs <- rbind(Curley.pcs,this.row)
# } else if (!is.na((as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]))) &
# !is.na((as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical]))))
# {
# this.row <- data.frame(
# Compound = Curley.compounds[,this.chemical],
# DTXSID = this.chemical,
# Tissue = this.tissue,
# PC = as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]) /
# as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
# )
# Curley.pcs <- rbind(Curley.pcs,this.row)
# }
# }
# }
# )
# Curley.pcs[Curley.pcs$Tissue=="Lungs","Tissue"] <- "Lung"
#
# pearce2017 <- read_excel("PC_Data.xlsx")
# pearce2017 <- subset(pearce2017,DTXSID%in%Curley.pcs$DTXSID)
# any(pearce2017$DTXSID%in%Curley.pcs$DTXSID)
# print("None of the Curley chems were in the Pearce 2017 PC predictor evaluation.")
## ----csanady2002_pcs, eval = execute.vignette---------------------------------
# csanadybpa <- read_excel("Csanady2002.xlsx",sheet="Table 2")
# csanadybpa$Compound <- "Bisphenol A"
# csanadybpa$DTXSID <- "DTXSID7020182"
# csanadybpa <- csanadybpa[,c("Compound","DTXSID","Tissue","Mean")]
# colnames(csanadybpa) <- colnames(Curley.pcs)
#
# csanadydaid <- read_excel("Csanady2002.xlsx",sheet="Table 3",skip=1)
# csanadydaid$Compound <- "Daidzein"
# csanadydaid$DTXSID <- "DTXSID9022310"
# csanadydaid <- csanadydaid[,c("Compound","DTXSID","Tissue","Mean...6")]
# colnames(csanadydaid) <- colnames(Curley.pcs)
#
# Curley.pcs <- rbind(Curley.pcs,csanadybpa,csanadydaid)
# Curley.pcs <- subset(Curley.pcs,Tissue!="Mammary gland")
## ----weijs2013_loadPCs, eval = execute.vignette-------------------------------
# weijstab3 <- read_excel("Weijs2013.xlsx",sheet="Table3",skip=1)
# weijstab3 <- subset(weijstab3, !is.na(Compound) & !is.na(Tissue))
# weijstab4 <- read_excel("Weijs2013.xlsx",sheet="Table4",skip=1)
# weijstab4 <- subset(weijstab4, !is.na(Compound) & !is.na(Tissue))
#
# for (this.compound in unique(weijstab3$Compound))
# for (this.tissue in unique(weijstab3$Tissue))
# {
# Curley.pcs[
# Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
# "Weijs2013"] <- weijstab3[weijstab3$Compound==this.compound &
# weijstab3$Tissue==this.tissue,"value"]
#
# }
#
# for (this.compound in unique(weijstab4$Compound))
# for (this.tissue in unique(weijstab4$Tissue))
# {
# Curley.pcs[
# Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
# "Weijs2013"] <- weijstab4[weijstab4$Compound==this.compound &
# weijstab4$Tissue==this.tissue,"value"]
#
# }
#
# write.csv(Curley.pcs,"PCs-table.csv",row.names=FALSE)
#
## ----curley1969_makepreds, eval = execute.vignette----------------------------
# Curley.pcs <- httk::fetalpcs
# dtxsid.list <- get_cheminfo(
# info="DTXSID",
# model="fetal_pbtk",
# suppress.messages=TRUE)
#
# suppressWarnings(load_sipes2017())
# for (this.chemical in unique(Curley.pcs$DTXSID))
# if (this.chemical %in% dtxsid.list)
# {
# this.subset <- subset(Curley.pcs,DTXSID==this.chemical)
# p <- parameterize_fetal_pbtk(dtxsid=this.chemical,
# fetal_fup_adjustment = FALSE,
# suppress.messages=TRUE,
# minimum.Funbound.plasma = 1e-16)
# fetal.blood.pH <- 7.26
# Fup <- p$Fraction_unbound_plasma_fetus
# fetal_schmitt_parms <- parameterize_schmitt(
# dtxsid=this.chemical,
# suppress.messages=TRUE,
# minimum.Funbound.plasma = 1e-16)
# fetal_schmitt_parms$plasma.pH <- fetal.blood.pH
# fetal_schmitt_parms$Funbound.plasma <- Fup
# Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Fup"] <- Fup
# # httk gives tissue:fup (unbound chemical in plasma) PC's:
# fetal_pcs <- predict_partitioning_schmitt(
# parameters = fetal_schmitt_parms,
# suppress.messages=TRUE,
# model="fetal_pbtk",
# minimum.Funbound.plasma = 1e-16)
# fetal_pcs.nocal <- predict_partitioning_schmitt(
# parameters = fetal_schmitt_parms,
# regression=FALSE,
# suppress.messages=TRUE,
# model="fetal_pbtk",
# minimum.Funbound.plasma = 1e-16)
# out <- solve_fetal_pbtk(
# dtxsid = this.chemical,
# fetal_fup_adjustment =FALSE,
# suppress.messages=TRUE,
# minimum.Funbound.plasma = 1e-16)
# Rb2p <- out[dim(out)[1],"Rfblood2plasma"]
# Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Rb2p"] <- Rb2p
# # Convert to tissue:blood PC's:
# for (this.tissue in this.subset$Tissue)
# if (tolower(this.tissue) %in%
# unique(subset(tissue.data,Species=="Human")$Tissue))
# {
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred.t2up"] <-
# fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal.t2up"] <-
# fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred"] <-
# fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal"] <-
# fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
# } else {
# print(this.tissue)
# }
# } else print(this.chemical)
# reset_httk()
#
# cat(paste(
# "For the two placental observations (",
# signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"PC"],2),
# " for Bisphenol A and ",
# signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"PC"],2),
# " for Diadzen) the predictions were ",
# signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"HTTK.pred"],2),
# " and ",
# signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"HTTK.pred"],2),
# ", respectively.\n",sep=""))
#
#
## ----curley1969_figure_compounds, eval = execute.vignette---------------------
# FigFa <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred.nocal))) +
# geom_point(size=2,aes(
# y=Weijs2013,
# x=HTTK.pred,
# shape=Compound,
# color=Compound)) +
# geom_point(size=5,aes(
# y=PC,
# x=HTTK.pred,
# shape=Compound,
# color=Compound)) +
# geom_abline(slope=1, intercept=0, size=2) +
# geom_abline(slope=1, intercept=1, linetype=3, size=2) +
# geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
# scale_shape_manual(values=rep(c(15,16,17,18),4))+
# xlab("Predicted Tissue:Blood Partition Coefificent") +
# ylab("Observed") +
# scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
# scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=28))+
# theme(legend.text=element_text(size=12))+
# theme(legend.title = element_blank())+
# guides(shape=guide_legend(nrow=3,byrow=TRUE))
#
# print(FigFa)
#
# fitFa <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred))
# RMSE(fitFa)
# fitFb <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred.nocal))
# RMSE(fitFb)
## ----curley1969_figure_tissues, eval = execute.vignette-----------------------
# FigFb <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred))) +
# geom_point(size=2,aes(
# y=Weijs2013,
# x=HTTK.pred,
# shape=Tissue,
# color=Tissue)) +
# geom_point(size=5, aes(
# y=PC,
# x=HTTK.pred,
# shape=Tissue,
# color=Tissue)) +
# geom_abline(slope=1, intercept=0, size=2) +
# geom_abline(slope=1, intercept=1, linetype=3, size=2) +
# geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
# scale_shape_manual(values=rep(c(15,16,17,18),4))+
# xlab("Predicted Tissue:Blood Partition Coefificent") +
# ylab("Observed") +
# scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
# scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=28))+
# theme(legend.text=element_text(size=20))+
# theme(legend.title = element_blank())+
# guides(shape=guide_legend(nrow=3,byrow=TRUE))
#
# print(FigFb)
## ----read_pksim_pcs, eval = execute.vignette----------------------------------
# pksim.pcs <- as.data.frame(read_excel("PartitionCoefficients.xlsx"))
# dim(pksim.pcs)
# for (this.id in unique(pksim.pcs$DTXSID))
# {
# httk.PCs <- predict_partitioning_schmitt(dtxsid=this.id,
# suppress.messages = TRUE)
# p <-
# parameterize_fetal_pbtk(dtxsid=this.id,
# suppress.messages = TRUE)
# httk.PCs[["Kplacenta2pu"]] <- p[["Kplacenta2pu"]]
# httk.PCs[["Fup"]] <- p[["Funbound.plasma"]]
# lapply(httk.PCs,as.numeric)
# this.subset <- subset(pksim.pcs,DTXSID==this.id)
# for (this.tissue in unique(this.subset$Tissue))
# {
# if (any(regexpr(tolower(this.tissue),names(httk.PCs))!=-1))
# {
# pksim.pcs[pksim.pcs$DTXSID==this.id & pksim.pcs$Tissue==this.tissue,
# "HTTK.pred"] <-
# httk.PCs[regexpr(tolower(this.tissue),names(httk.PCs))!=-1][[1]]*
# httk.PCs[["Fup"]][[1]]
# }
# }
# }
# colnames(pksim.pcs)[colnames(pksim.pcs) ==
# "Tissue-to-plasma partition coefficient"] <- "PKSim.pred"
# colnames(pksim.pcs)[colnames(pksim.pcs) ==
# "Method for calculating tissue-to-plasma partition coefficients"] <-
# "Method"
# pksim.pcs[,"Method"] <- as.factor(pksim.pcs[,"Method"])
# pksim.pcs[,"Drug"] <- as.factor(pksim.pcs[,"Drug"])
# pksim.pcs[,"Tissue"] <- as.factor(pksim.pcs[,"Tissue"])
#
#
# pksim.pcs[,"PKSim.pred"] <- as.numeric(
# pksim.pcs[,"PKSim.pred"])
#
## ----compare_pksim_pcs, eval = execute.vignette-------------------------------
# pksim.pcs <- httk::pksim.pcs
#
# pksim.fit1 <- lm(data=pksim.pcs,
# log10(PKSim.pred)~log10(HTTK.pred))
# summary(pksim.fit1)
#
# pksim.fit2 <- lm(data=pksim.pcs,
# log10(PKSim.pred)~log10(HTTK.pred)+
# Drug+Tissue+Method)
# summary(pksim.fit2)
#
## ----wang2018_loaddata, eval = execute.vignette-------------------------------
# #Wangchems <- read_excel("Wang2018.xlsx",sheet=3,skip=2)
# Wangchems <- httk::wang2018
# maternal.list <- Wangchems$CASRN[Wangchems$CASRN %in%
# get_cheminfo(model="fetal_pbtk",
# suppress.messages=TRUE)]
# nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$CAS
# nonfluoros <- chem.physical_and_invitro.data$CAS[
# regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))==-1]
# maternal.list <- maternal.list[maternal.list %in% intersect(nonvols,nonfluoros)]
## ----wang2018_makepreds, eval = execute.vignette------------------------------
#
# pred.table1 <- subset(get_cheminfo(
# info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
# model="fetal_pbtk"),
# CAS %in% maternal.list,
# suppress.messages=TRUE)
# pred.table1$Compound <- gsub("\"","",pred.table1$Compound)
#
# for (this.cas in maternal.list)
# {
# Fup <- subset(get_cheminfo(
# info="all",
# suppress.messages=TRUE),
# CAS==this.cas)$Human.Funbound.plasma
# # Check if Fup is provided as a distribution:
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
# {
# skip <- TRUE
# } else skip <- FALSE
# Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
# Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
# # These results are actually from a Bayesian posterior, but we can approximate that
# # they are normally distributed about the median to estimate a standard deviation:
# Fup.sigma <- (Fup.upper - Fup)/1.96
# # If it's not a distribution use defaults (see Wambaugh et al. 2019)
# } else if (as.numeric(Fup)== 0)
# {
# skip <- TRUE
# } else {
# skip <- FALSE
# Fup <- as.numeric(Fup)
# Fup.sigma <- Fup*0.4
# Fup.upper <- min(1,(1+1.96*0.4)*Fup)
# }
#
# # Only run the simulation if we have the necessary parameters:
# if (!skip)
# {
# print(this.cas)
# out <- solve_fetal_pbtk(
# chem.cas=this.cas,
# dose=0,
# daily.dose=1,
# doses.per.day=3,
# fetal_fup_adjustenment=TRUE,
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# pred.table1[pred.table1$CAS==this.cas,"fup"] <- signif(Fup,3)
# pred.table1[pred.table1$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
# pred.table1[pred.table1$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
# pred.table1[pred.table1$CAS==this.cas,"Ratio.pred"] <-
# signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
# }
# }
## ----brain_uncertainty_prop_table, eval = execute.vignette--------------------
# PC.table <- pred.table1
# colnames(PC.table)[colnames(PC.table)=="Ratio.pred"] <- "R.plasma.FtoM"
# pred.table1$Uncertainty <- "Predicted F:M Plasma Ratio"
## ----Rfetmat_uncertainty, eval = execute.vignette-----------------------------
# pred.table3 <- pred.table1
# pred.table3$Uncertainty <- "Plasma Error (Fig. 4)"
# empirical.error <- RMSE(fit2sub)
# for (this.cas in maternal.list)
# {
#
# pred.table3[pred.table3$CAS==this.cas,"Ratio.pred"] <- signif(
# 1/((1/pred.table1[pred.table3$CAS==this.cas,"Ratio.pred"])*
# (1-1.96*empirical.error)), 3)
# }
# # Update final table for paper:
# PC.table$RMSE <- signif(RMSE(fit2sub),3)
# PC.table$R.plasma.FtoM.upper <- pred.table3$Ratio.pred
## ----Rfetbrain, eval = execute.vignette---------------------------------------
# pred.table4 <- pred.table1
# pred.table4$Uncertainty <- "Fetal Brain Partitioning"
# for (this.cas in maternal.list)
# {
# print(this.cas)
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=TRUE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
# fup <- pred.table4[pred.table4$CAS==this.cas,"fup"]
# pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"] <- signif(
# PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]*
# Kbrain2pu * fup, 3)
# PC.table[PC.table$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
# PC.table[PC.table$CAS==this.cas,"R.brain.FtoM"] <-
# pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"]
# }
## ----Rfetbrain_uncertainty, eval = execute.vignette---------------------------
#
# pred.table5 <- pred.table1
# pred.table5$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
# for (this.cas in maternal.list)
# {
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=TRUE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
#
# fup <- pred.table5[pred.table5$CAS==this.cas,"fup"]
# sigma.fup <- pred.table5[pred.table5$CAS==this.cas,"fup.sigma"]
# Rmatfet <- 1/PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]
# Rbrain2matblood <- Kbrain2pu * fup / Rmatfet
#
# # From Pearce et al. (2017) PC paper:
# sigma.logKbrain <- 0.647
# Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)
#
# error.Rmatfet <- PC.table[PC.table$CAS==this.cas,"RMSE"]
# sigma.Rbrain2matblood <- Rbrain2matblood *
# (log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2)
# Rbrain2matblood.upper <- Rbrain2matblood *
# (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2))
# pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"] <-
# signif(Rbrain2matblood.upper,3)
# PC.table[PC.table$CAS==this.cas,"sigma.logKbrain"] <-
# signif(sigma.logKbrain, 3)
# PC.table[PC.table$CAS==this.cas,"Kbrain2pu.upper"] <-
# signif(Kbrain2pu.upper, 3)
# PC.table[PC.table$CAS==this.cas,"sigma.Rbrain2matblood"] <-
# signif(sigma.Rbrain2matblood, 3)
# PC.table[PC.table$CAS==this.cas,"R.brain.FtoM.upper"] <-
# pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"]
# }
## ----Wang_Uncertainty_Propagation, eval = execute.vignette--------------------
#
# pred.levels <- pred.table5$Compound[order(pred.table5$Ratio.pred)]
#
# pred.table <- rbind(
# pred.table1,
# # pred.table2,
# pred.table3,
# pred.table4,
# pred.table5)
# pred.table$Compound <- factor(pred.table$Compound,
# levels = pred.levels)
#
# pred.table$Uncertainty <- factor(pred.table$Uncertainty,
# levels = c(pred.table1[1,"Uncertainty"],
# # pred.table2[1,"Uncertainty"],
# pred.table3[1,"Uncertainty"],
# pred.table4[1,"Uncertainty"],
# pred.table5[1,"Uncertainty"]))
## ----wang2018_figure, eval = execute.vignette---------------------------------
# #Wang 2018 confirmed 6 chemical identities:
# confirmed.chemicals <- c(
# "2,4-Di-tert-butylphenol",
# "2,4-Dinitrophenol",
# "Pyrocatechol",
# "2'-Hydroxyacetophenone",
# "3,5-Di-tert-butylsalicylic acid",
# "4-Hydroxycoumarin"
# )
# confirmed.chemicals <- c(
# "96-76-4",
# "19715-19-6",
# "51-28-5",
# "120-80-9",
# "118-93-4",
# "1076-38-6")
#
#
# FigG <- ggplot(data=pred.table) +
# geom_point(aes(
# x=Ratio.pred,
# y=Compound,
# color = Uncertainty,
# shape = Uncertainty),
# size=3) +
# scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
# scale_x_log10(limits=c(10^-2,10^3),label=scientific_10)+
# ylab(expression(paste(
# "Chemicals Found in Maternal Plasma by Wang ",italic("et al.")," (2018)"))) +
# xlab("Predicted Ratio to Maternal Plasma") +
# theme_bw() +
# # theme(legend.position="bottom")+
# theme( text = element_text(size=14))+
# theme(legend.text=element_text(size=10))#+
# # guides(color=guide_legend(nrow=4,byrow=TRUE))+
# #guides(shape=guide_legend(nrow=4,byrow=TRUE))
# #+
# # theme(legend.justification = c(0, 0), legend.position = c(0, 0))
#
# print(FigG)
## ----wang2018_details, eval = execute.vignette--------------------------------
#
# for (this.col in 7:14)
# PC.table[,this.col] <- signif(PC.table[,this.col],3)
#
# PC.table <- PC.table[order(PC.table$R.brain.FtoM.upper,decreasing=TRUE),]
#
# for (this.row in 1:dim(PC.table)[1])
# {
# out <- calc_ionization(
# pH=7.26,
# pKa_Donor=PC.table[this.row,"pKa_Donor"],
# pKa_Accept=PC.table[this.row,"pKa_Accept"])
# if (out$fraction_neutral>0.9) PC.table[this.row,"Charge_726"] <- "Neutral"
# else if (out$fraction_positive>0.1) PC.table[this.row,"Charge_726"] <-
# paste(signif(out$fraction_positive*100,2),"% Positive",sep="")
# else if (out$fraction_negative>0.1) PC.table[this.row,"Charge_726"] <-
# paste(signif(out$fraction_negative*100,2),"% Negative",sep="")
# else if (out$fraction_zwitter>0.1) PC.table[this.row,"Charge_726"] <-
# paste(signif(out$fraction_zwitter*100,2),"% Zwitterion",sep="")
# }
#
# PC.table <- PC.table[,c(
# "Compound",
# "CAS",
# "DTXSID",
# "logP",
# "Charge_726",
# "R.plasma.FtoM",
# "RMSE",
# "R.plasma.FtoM.upper",
# "Kbrain2pu",
# "fup",
# "R.brain.FtoM",
# "Kbrain2pu.upper",
# "R.brain.FtoM.upper")]
#
# write.csv(PC.table,
# file="WangTable.txt",
# row.names=FALSE)
## ----impact_of_model, eval = execute.vignette---------------------------------
# MFbrainfit <- lm(R.brain.FtoM.upper~Kbrain2pu.upper,data=PC.table)
# summary(MFbrainfit)
#
# cat(paste("As expected, the predicted fetal brain to maternal blood ratio was strongly correlated with the predicted brain partition coefficient (R2 = ",
# signif(summary(MFbrainfit)$adj.r.square,2),
# ", p-value ",
# signif(summary(MFbrainfit)$coefficients[2,4],2),
# "). However, the PBTK simulation impacted the plasma and tissue concentrations such that ",
# dim(subset(PC.table, rank(R.brain.FtoM.upper) < rank(Kbrain2pu.upper)))[1],
# " chemicals were ranked higher than they would have been using partitioning alone.",
# sep=""))
## ----wang2018_makepreds_nofup, eval = execute.vignette------------------------
#
# pred.table1.nofup <- subset(get_cheminfo(
# info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
# model="fetal_pbtk",
# suppress.messages=TRUE),
# CAS %in% maternal.list)
# pred.table1.nofup$Compound <- gsub("\"","",pred.table1.nofup$Compound)
#
# for (this.cas in maternal.list)
# {
# Fup <- subset(get_cheminfo(
# info="all",
# suppress.messages=TRUE),
# CAS==this.cas)$Human.Funbound.plasma
# # Check if Fup is provided as a distribution:
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
# {
# skip <- TRUE
# } else skip <- FALSE
# Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
# Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
# # These results are actually from a Bayesian posterior, but we can approximate that
# # they are normally distributed about the median to estimate a standard deviation:
# Fup.sigma <- (Fup.upper - Fup)/1.96
# # If it's not a distribution use defaults (see Wambaugh et al. 2019)
# } else if (as.numeric(Fup)== 0)
# {
# skip <- TRUE
# } else {
# skip <- FALSE
# Fup <- as.numeric(Fup)
# Fup.sigma <- Fup*0.4
# Fup.upper <- min(1,(1+1.96*0.4)*Fup)
# }
#
# # Only run the simulation if we have the necessary parameters:
# if (!skip)
# {
# out <- solve_fetal_pbtk(
# chem.cas=this.cas,
# dose=0,
# daily.dose=1,
# doses.per.day=3,
# fetal_fup_adjustenment=FALSE,
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup"] <- signif(Fup,3)
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"Ratio.pred"] <-
# signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
# }
# }
## ----brain_uncertainty_prop_table_nofup, eval = execute.vignette--------------
# PC.table.nofup <- pred.table1.nofup
# colnames(PC.table.nofup)[colnames(PC.table.nofup)=="Ratio.pred"] <- "R.plasma.FtoM"
# pred.table1.nofup$Uncertainty <- "Predicted F:M Plasma Ratio"
## ----Rfetmat_uncertainty_nofup, eval = execute.vignette-----------------------
# pred.table3.nofup <- pred.table1.nofup
# pred.table3.nofup$Uncertainty <- "Plasma Error (Fig. 4)"
# empirical.error <- RMSE(fit2sub)
# for (this.cas in maternal.list)
# {
#
# pred.table3.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
# 1/((1/pred.table1.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"])*
# (1-1.96*empirical.error)), 3)
# }
# # Update final table for paper:
# PC.table.nofup$RMSE <- signif(RMSE(fit2sub),3)
# PC.table.nofup$R.plasma.FtoM.upper <- pred.table3.nofup$Ratio.pred
## ----Rfetbrain_nofup, eval = execute.vignette---------------------------------
# pred.table4.nofup <- pred.table1.nofup
# pred.table4.nofup$Uncertainty <- "Fetal Brain Partitioning"
# for (this.cas in maternal.list)
# {
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=FALSE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
# fup <- pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"fup"]
# pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
# PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]*
# Kbrain2pu * fup, 3)
# PC.table.nofup[PC.table.nofup$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
# PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.brain.FtoM"] <-
# pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"]
# }
## ----Rfetbrain_uncertainty_nofup, eval = execute.vignette---------------------
#
# pred.table5.nofup <- pred.table1.nofup
# pred.table5.nofup$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
# for (this.cas in maternal.list)
# {
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=FALSE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
#
# fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup"]
# sigma.fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup.sigma"]
# Rmatfet <- 1/PC.table[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]
# Rbrain2matblood <- Kbrain2pu * fup / Rmatfet
#
# # From Pearce et al. (2017) PC paper:
# sigma.logKbrain <- 0.647
# Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)
#
# error.Rmatfet <- PC.table.nofup [PC.table.nofup $CAS==this.cas,"RMSE"]
# sigma.Rbrain2matblood <- Rbrain2matblood *
# (log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2)
# Rbrain2matblood.upper <- Rbrain2matblood *
# (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2))
# pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"] <-
# signif(Rbrain2matblood.upper,3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.logKbrain"] <-
# signif(sigma.logKbrain, 3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"Kbrain2pu.upper"] <-
# signif(Kbrain2pu.upper, 3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.Rbrain2matblood"] <-
# signif(sigma.Rbrain2matblood, 3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"R.brain.FtoM.upper"] <-
# pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"]
# }
## ----compare_fup_correction_case_study,eval=FALSE-----------------------------
# case.study.fup.correct.table <- merge(PC.table,PC.table.nofup,by=c("Compound","CAS","DTXSID"))
#
# MFbrainfit.fup <- lm(R.brain.FtoM.upper.x~R.brain.FtoM.upper.y,
# data=case.study.fup.correct.table)
# summary(MFbrainfit.fup)
#
# cat(paste("The predictions for fetal brain to maternal blood ratio with or without the fetal plasma binding correction (Eqn. 1) were strongly correlated (R2 = ",
# signif(summary(MFbrainfit.fup)$adj.r.square,2),
# "). There were ",
# dim(subset(case.study.fup.correct.table,
# rank(R.brain.FtoM.upper.x) > rank(R.brain.FtoM.upper.y)))[1],
# " chemicals that were ranked higher with the correction than without, with an average increase of ",
# signif(mean(
# case.study.fup.correct.table$R.brain.FtoM.upper.y /
# case.study.fup.correct.table$R.brain.FtoM.upper.x),2),
# " times when the plasma binding change was included.\n",
# sep=""))
#
## ----fup_table, eval = execute.vignette---------------------------------------
# fup.table <- NULL
# all.chems <- get_cheminfo(
# model="fetal_pbtk",
# info="all",
# suppress.messages=TRUE)
# # Get rid of median fup 0:
# all.chems <- subset(all.chems,
# as.numeric(unlist(lapply(strsplit(
# all.chems$Human.Funbound.plasma,","),function(x) x[[1]])))!=0)
# for (this.chem in all.chems[,"CAS"])
# {
# temp <- parameterize_fetal_pbtk(chem.cas=this.chem,suppress.messages = TRUE)
# state <- calc_ionization(
# pH=7.26,
# pKa_Donor=temp$pKa_Donor,
# pKa_Accept=temp$pKa_Accept)
# if (state$fraction_positive > 0.5) this.charge <- "Positive"
# else if (state$fraction_negative > 0.5) this.charge <- "Negative"
# else this.charge <- "Neutral"
# this.row <- data.frame(DTXSID=all.chems[all.chems[,"CAS"]==this.chem,"DTXSID"],
# Compound=all.chems[all.chems[,"CAS"]==this.chem,"Compound"],
# CAS=this.chem,
# Fup.Mat.Pred = temp$Funbound.plasma,
# Fup.Neo.Pred = temp$Fraction_unbound_plasma_fetus,
# Charge = this.charge
# )
# fup.table <- rbind(fup.table,this.row)
# }
#
# fup.table[fup.table$Charge=="Positive","Charge"] <- paste("Positive (n=",
# sum(fup.table$Charge=="Positive"),
# ")",sep="")
# fup.table[fup.table$Charge=="Negative","Charge"] <- paste("Negative (n=",
# sum(fup.table$Charge=="Negative"),
# ")",sep="")
# fup.table[fup.table$Charge=="Neutral","Charge"] <- paste("Neutral (n=",
# sum(fup.table$Charge=="Neutral"),
# ")",sep="")
## ----fup_figure, eval = execute.vignette--------------------------------------
# FigA <- ggplot(data=fup.table) +
# geom_point(alpha=0.25, aes(
# x=Fup.Mat.Pred,
# y=Fup.Neo.Pred,
# shape=Charge,
# color=Charge),
# size=3) +
# geom_abline(slope=1, intercept=0) +
# ylab(expression(paste("Adjusted Neonate ",f[up]))) +
# xlab(expression(paste(italic("In vitro")," Measured Adult ",f[up]))) +
# scale_x_log10(label=scientific_10) +
# scale_y_log10(label=scientific_10) +
# theme_bw() +
# theme( text = element_text(size=14))
#
# print(FigA)
## ----MFratio_allhttk_preds, eval = execute.vignette---------------------------
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
#
# MFratio.pred <- NULL
# all.chems <- get_cheminfo(
# model="fetal_pbtk",
# info=c("Compound","DTXSID","CAS","Funbound.plasma"),
# suppress.messages=TRUE)
# for (this.cas in all.chems$CAS)
# if ((this.cas %in% nonvols) &
# !(this.cas %in% fluoros))
# {
# this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
# Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
# {
# skip <- TRUE
# } else skip <- FALSE
# } else if (Fup== 0)
# {
# skip <- TRUE
# } else skip <- FALSE
#
# if (!skip)
# {
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =TRUE,
# suppress.messages=TRUE)
# out <- solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# new.row <- data.frame(
# Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
# DTXSID = this.id,
# Mat.pred = mean(out[last.row,"Cplasma"]),
# Fet.pred = mean(out[last.row,"Cfplasma"]),
# MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# )
# MFratio.pred <- rbind(MFratio.pred,new.row)
# }
# }
## ----MFratio_allhttk_figure, eval = execute.vignette--------------------------
# FigD <- ggplot(data=MFratio.pred)+
# geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+
# xlab("Maternal:Fetal Plasma Concentration Ratio") +
# ylab("Number of chemicals")+
# theme_bw()+
# theme( text = element_text(size=14))
#
#
# print(FigD)
## ----MFratio_allhttk_stats, eval = execute.vignette---------------------------
# max.chem <- MFratio.pred[which(
# MFratio.pred$MFratio.pred==max(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
# min.chem <- MFratio.pred[which(
# MFratio.pred$MFratio.pred==min(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
# cat(paste("In Figure X we examine the ratios predicted for the ",
# dim(MFratio.pred)[1],
# " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
# sep=""))
#
#
# cat(paste("We observe a median ratio of ",
# signif(median(MFratio.pred$MFratio.pred,na.rm=TRUE),3),
# ", with values ranging from ",
# signif(min.chem[,"MFratio.pred"],3),
# " for ",
# min.chem[,"DTXSID"],
# " to ",
# signif(max.chem[,"MFratio.pred"],3),
# " for ",
# max.chem[,"DTXSID"],
# ".\n",sep=""))
#
# # Check out phys-chem > 1.6, < 1:
# highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred>1.6)$DTXSID)
#
# cat(paste("There are",
# dim(highratio)[1],
# "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))
#
# # all highly bound
# highratio$Compound
# suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=TRUE)))
#
#
# lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred<0.9)$DTXSID)
# # No obvious pattern
## ----MFratio_allhttk_preds_nofup, eval = execute.vignette---------------------
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
#
# MFratio.pred.nofup <- NULL
# all.chems <- get_cheminfo(
# model="fetal_pbtk",
# info=c("Compound","DTXSID","CAS","Funbound.plasma"),
# suppress.messages=TRUE)
# for (this.cas in all.chems$CAS)
# if ((this.cas %in% nonvols) &
# !(this.cas %in% fluoros))
# {
# this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
# Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
# {
# skip <- TRUE
# } else skip <- FALSE
# } else if (Fup== 0)
# {
# skip <- TRUE
# } else skip <- FALSE
#
# if (!skip)
# {
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =FALSE,
# suppress.messages=TRUE))
# out <- suppressWarnings(solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# new.row <- data.frame(
# Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
# DTXSID = this.id,
# Mat.pred = mean(out[last.row,"Cplasma"]),
# Fet.pred = mean(out[last.row,"Cfplasma"]),
# MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# )
# MFratio.pred.nofup <- rbind(MFratio.pred.nofup,new.row)
# }
# }
## ----MFratio_allhttk_figure_nofup, eval = execute.vignette--------------------
# FigSupD <- ggplot(data=MFratio.pred.nofup)+
# geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+
# xlab("Maternal:Fetal Plasma Concentration Ratio (No Fup Corr.)") +
# ylab("Number of chemicals")+
# theme_bw()+
# theme( text = element_text(size=14))
#
#
# print(FigSupD)
## ----MFratio_allhttk_stats_nofup, eval = execute.vignette---------------------
# max.chem <- MFratio.pred.nofup[which(
# MFratio.pred.nofup$MFratio.pred==max(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
# min.chem <- MFratio.pred.nofup[which(
# MFratio.pred.nofup$MFratio.pred==min(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
# cat(paste("In Figure X we examine the ratios predicted for the ",
# dim(MFratio.pred)[1],
# " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
# sep=""))
#
#
# cat(paste("We observe a median ratio of ",
# signif(median(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE),3),
# ", with values ranging from ",
# signif(min.chem[,"MFratio.pred"],3),
# " for ",
# min.chem[,"DTXSID"],
# " to ",
# signif(max.chem[,"MFratio.pred"],3),
# " for ",
# max.chem[,"DTXSID"],
# ".\n",sep=""))
#
# # Check out phys-chem > 1.6, < 1:
# highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred>1.6)$DTXSID)
#
# cat(paste("There are",
# dim(highratio)[1],
# "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))
#
# # all highly bound
# highratio$Compound
# suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=2)))
#
#
# lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred<0.9)$DTXSID)
# # No obvious pattern
## ----create_distributable_data, eval = execute.vignette-----------------------
# aylward2014 <-MFdata
# pregnonpregaucs <- TKstats
# fetalpcs <- Curley.pcs
# wang2018 <- Wangchems
#
# save(aylward2014,pregnonpregaucs,fetalpcs,wang2018,pksim.pcs,
# file="Kapraun2022Vignette.RData",version=2)
Try the httk package in your browser
Any scripts or data that you put into this service are public.
httk documentation built on Sept. 11, 2024, 9:32 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
## ----configure_knitr, eval = TRUE---------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = '#>')
options(rmarkdown.html_vignette.check_title = FALSE)
## ----clear_memory, eval = TRUE------------------------------------------------
rm(list=ls())
## ----runchunks, eval = TRUE---------------------------------------------------
# Set whether or not the following chunks will be executed (run):
execute.vignette <- FALSE
## ----load_packages, eval = execute.vignette-----------------------------------
# #
# #
# # Setup
# #
# #
# #library(readxl)
# library(ggplot2)
# library(httk)
# library(scales)
# library(gridExtra)
# library(cowplot)
## ----load_code, eval = execute.vignette---------------------------------------
# scientific_10 <- function(x) {
# out <- gsub("1e", "10^", scientific_format()(x))
# out <- gsub("\\+","",out)
# out <- gsub("10\\^01","10",out)
# out <- parse(text=gsub("10\\^00","1",out))
# }
#
# RMSE <- function(x)
# {
# mean(x$residuals^2)^(1/2)
# }
## ----chem_numbers, eval=FALSE-------------------------------------------------
# length(get_cheminfo(model="fetal_pbtk",suppress.messages=TRUE))
## ----load_aylward_data, eval = execute.vignette-------------------------------
# #MFdata <- read_excel("Aylward-MatFet.xlsx")
# MFdata <- httk::aylward2014
#
# cat(paste("summarized data from over 100 studies covering ",
# length(unique(MFdata$DTXSID)[!(unique(MFdata$DTXSID)%in%c("","-"))]),
# " unique chemicals structures\n",sep=""))
#
# # We don't want to exclude the volatiles and PFAS at this point:
# MFdata.httk <- subset(MFdata,DTXSID %in% get_cheminfo(
# info="DTXSID",
# suppress.messages=TRUE))
# MFdata.httk[MFdata.httk$Chemical.Category=="bromodiphenylethers",
# "Chemical.Category"] <- "Bromodiphenylethers"
# MFdata.httk[MFdata.httk$Chemical.Category=="organochlorine Pesticides",
# "Chemical.Category"] <- "Organochlorine Pesticides"
# MFdata.httk[MFdata.httk$Chemical.Category=="polyaromatic hydrocarbons",
# "Chemical.Category"] <- "Polyaromatic Hydrocarbons"
#
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
# "infant.maternal.conc...Central.tendency..calculate.j.k..or.report.paired.result."] <-
# "MFratio"
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
# "PREFERRED_NAME"] <-
# "Chemical"
# colnames(MFdata.httk)[colnames(MFdata.httk) ==
# "details.on.matrix.comparison...e.g...cord.blood.lipid..maternal.serum.lipid..or.cord.blood.wet.weight..maternal.whole.blood.wet.weight"] <-
# "Matrix"
#
# # Format the columns:
# MFdata.httk$MFratio <- as.numeric(MFdata.httk$MFratio)
# MFdata.httk$Chemical <- as.factor(MFdata.httk$Chemical)
# MFdata.httk$Matrix <- as.factor(MFdata.httk$Chemical)
# MFdata.httk$Chemical.Category <- as.factor(MFdata.httk$Chemical.Category)
#
# colnames(MFdata.httk)[15] <- "infant"
# colnames(MFdata.httk)[16] <- "maternal"
# colnames(MFdata.httk)[17] <- "obs.units"
#
# MFdata.httk$infant <- suppressWarnings(as.numeric(MFdata.httk$infant))
# MFdata.httk$maternal <- suppressWarnings(as.numeric(MFdata.httk$maternal))
# MFdata.httk$AVERAGE_MASS <- as.numeric(MFdata.httk$AVERAGE_MASS)
## ----process_aylward_data, eval = execute.vignette----------------------------
# convert1.units <- c("ng/ml","ng/mL","ug/L","ug/l","ng/mL serum","ng/g",
# "ng/g wet wt.","ppb")
#
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"infant"] / # ng/ml = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"maternal"] / # ng/ml = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert1.units,"obs.units"] <- "uM"
#
# convert2.units <- c("mg/L","ppm")
#
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"infant"] * 1000 / # mg/L = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"] <-
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"maternal"]* 1000 / # mg/L = ug/L
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"AVERAGE_MASS"] # ug/L -> uM
# MFdata.httk[MFdata.httk$obs.units%in%convert2.units,"obs.units"] <- "uM"
## ----make_aylward_preds, eval = execute.vignette------------------------------
# # Simulate starting from the 13th week of gestation until full term (40 weeks):
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
#
# # For each chemical with maternal-fetal ratio data and HTTK in vitro data:
# for (this.id in unique(MFdata.httk$DTXSID))
# {
# print(this.id)
# p <- parameterize_steadystate(dtxsid=this.id,
# suppress.messages=TRUE)
# # skip chemicals where Fup is 0:
# if (p$Funbound.plasma>1e-4)
# {
#
# # Load the chemical-specifc paramaters:
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =TRUE,
# suppress.messages=TRUE)
# # Skip chemicals where the 95% credible interval for Fup spans more than 0.1 to
# # 0.9 (that is, Fup is effectively unknown):
# if (!is.na(p$Funbound.plasma.dist))
# {
# if (as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][3])>0.9 &
# as.numeric(strsplit(p$Funbound.plasma.dist,",")[[1]][2])<0.11)
# {
# skip <- TRUE
# } else skip <- FALSE
# } else skip <- FALSE
#
# if (!skip)
# {
# # Solve the PBTK equations for the full simulation time assuming 1 total daily
# # dose of 1 mg/kg/day spread out over three evenly spaces exposures:
# out <- solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# suppress.messages=TRUE)
# # Identify the concentrations from the final (279th) day:
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# # Average over the final day:
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred"] <- mean(out[last.row,"Cplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred"] <- mean(out[last.row,"Cfplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred"] <-
# mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
#
# # Repeat the analysis without the adjustment to fetal Fup:
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =FALSE,
# suppress.messages = TRUE)
# out <- solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# maxsteps=1e7,
# suppress.messages = TRUE)
#
# last.row <- which(out[,"time"]>279) # The whole final day
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Mat.pred.nofup"] <- mean(out[last.row,"Cplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"Fet.pred.nofup"] <- mean(out[last.row,"Cfplasma"])
# MFdata.httk[MFdata.httk$DTXSID==this.id,"MFratio.pred.nofup"] <-
# mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# }
# }
# }
#
# # Something is wrong with cotinine:
# MFdata.httk <- subset(MFdata.httk,Chemical!="Cotinine")
#
# max.chem <- MFdata.httk[which(MFdata.httk$MFratio==max(MFdata.httk$MFratio)),]
# min.chem <- MFdata.httk[which(MFdata.httk$MFratio==min(MFdata.httk$MFratio)),]
# cat(paste("The minimum observed ratio was ",
# signif(min.chem[,"MFratio"],2),
# " for ",
# min.chem[,"Chemical"],
# " and the maximum was ",
# signif(max.chem[,"MFratio"],2),
# " for ",
# max.chem[,"Chemical"],
# ".\n",sep=""))
#
#
## ----cleanup_repeated_aylward_figure, eval = execute.vignette-----------------
# MFdata.main <- NULL
# MFdata.outliers <- NULL
# for (this.id in unique(MFdata.httk$DTXSID))
# {
# this.subset <- subset(MFdata.httk,DTXSID==this.id)
# this.row <- this.subset[1,]
# this.row$N.obs <- dim(this.subset)[1]
# this.row$MFratio <- median(this.subset$MFratio)
# this.row$MFratio.Q25 <- quantile(this.subset$MFratio,0.25)
# this.row$MFratio.Q75 <- quantile(this.subset$MFratio,0.75)
# MFdata.main <- rbind(MFdata.main,this.row)
# this.subset <- subset(this.subset,
# MFratio<this.row$MFratio.Q25 |
# MFratio>this.row$MFratio.Q75)
# MFdata.outliers <- rbind(MFdata.outliers,this.subset)
# }
## ----aylward_figure_no_fup, eval = execute.vignette---------------------------
# FigCa <- ggplot(data=MFdata.main) +
# geom_segment(color="grey",aes(
# x=MFratio.pred.nofup,
# y=MFratio.Q25,
# xend=MFratio.pred.nofup,
# yend=MFratio.Q75))+
# geom_point(aes(
# x=MFratio.pred.nofup,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=4) +
# scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
# geom_point(data=MFdata.outliers,aes(
# x=MFratio.pred.nofup,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=2) +
# xlim(0,2) +
# ylim(0,3) +
# # geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
# # scale_y_log10(label=scientific_10) +
# #,limits=c(10^-7,100)) +
# # scale_x_log10(label=scientific_10) +
# # ,limits=c(10^-7,100)) +
# # annotation_logticks() +
# geom_abline(slope=1, intercept=0) +
# # geom_abline(slope=1, intercept=1,linetype="dashed") +
# # geom_abline(slope=1, intercept=-1,linetype="dashed") +
# ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) +
# xlab("Generic PBTK Predicted Ratio") +
# # scale_color_brewer(palette="Set2") +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=14))+
# theme(legend.text=element_text(size=10))+
# guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+
# guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))
#
# print(FigCa)
## ----aylward_analysis_text, eval = execute.vignette---------------------------
#
# cat(paste("In Figure 4 we compare predictions made with our high-throughput \
# human gestational PBTK model with experimental observations on a per chemical \
# basis wherever we had both in vitro HTTK data and in vivo observations (",
# length(unique(MFdata.main$DTXSID)),
# " chemicals).\n",sep=""))
#
#
# repeats <- subset(MFdata.main,N.obs>1)
#
# cat(paste("Multiple observations were available for ",
# dim(repeats)[1],
# " of the chemicals,\n",sep=""))
#
#
# max.chem <- as.data.frame(repeats[which(repeats$MFratio==max(repeats$MFratio)),])
# min.chem <- as.data.frame(repeats[which(repeats$MFratio==min(repeats$MFratio)),])
#
# cat(paste("However, among the chemicals with repeated observations, the median \
# observations only ranged from ",
# signif(min.chem[,"MFratio"],2),
# " for ",
# min.chem[,"Chemical"],
# " to ",
# signif(max.chem[,"MFratio"],2),
# " for ",
# max.chem[,"Chemical"],
# ".\n",sep=""))
#
# max.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==max(MFdata.main$MFratio.pred,na.rm=TRUE)),])
# min.chem <- as.data.frame(MFdata.main[which(MFdata.main$MFratio.pred==min(MFdata.main$MFratio.pred,na.rm=TRUE)),])
#
# cat(paste("The predictions for all chemicals ranged from ",
# signif(min.chem[,"MFratio.pred"],2),
# " for ",
# min.chem[,"Chemical"],
# " to ",
# signif(max.chem[,"MFratio.pred"],2),
# " for ",
# max.chem[,"Chemical"],
# ".\n",sep=""))
#
#
#
# fit1 <- lm(data=MFdata.main,MFratio~MFratio.pred.nofup)
# summary(fit1)
# RMSE(fit1)
#
# fit1sub <- lm(data=subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred.nofup)
# summary(fit1sub)
# RMSE(fit1sub)
#
## ----aylward_figure_with_fup, eval = execute.vignette-------------------------
# FigCb <- ggplot(data=MFdata.main) +
# geom_segment(color="grey",aes(
# x=MFratio.pred,
# y=MFratio.Q25,
# xend=MFratio.pred,
# yend=MFratio.Q75))+
# geom_point(aes(
# x=MFratio.pred,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=3) +
# scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
# geom_point(data=MFdata.outliers,aes(
# x=MFratio.pred,
# y=MFratio,
# shape=Chemical.Category,
# color=Chemical.Category),
# size=1) +
# xlim(0,2) +
# ylim(0,3) +
# # geom_text(aes(x=AUC,y=Critical.concentration,label=Compound.abbrev,color=Chemical)) +
# # scale_y_log10(label=scientific_10) +
# #,limits=c(10^-7,100)) +
# # scale_x_log10(label=scientific_10) +
# # ,limits=c(10^-7,100)) +
# # annotation_logticks() +
# geom_abline(slope=1, intercept=0) +
# # geom_abline(slope=1, intercept=1,linetype="dashed") +
# # geom_abline(slope=1, intercept=-1,linetype="dashed") +
# ylab(expression(paste(italic("In vivo")," Mat:Fet Plasma Ratio"))) +
# xlab(expression(paste(italic("In vitro")," Predicted Ratio"))) +
# # scale_color_brewer(palette="Set2") +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=14))+
# theme(legend.text=element_text(size=10))+
# guides(color=guide_legend(title="Class",nrow=3,byrow=TRUE))+
# guides(shape=guide_legend(title="Class",nrow=3,byrow=TRUE))
#
# print(FigCb)
#
## ----voc_analysis, eval = execute.vignette------------------------------------
# # Mean logHenry's law constant for PAH's:
# mean(subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFdata.main,Chemical.Category=="Polyaromatic Hydrocarbons")$DTXSID)$logHenry)
#
# nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$DTXSID
# fluoros <- chem.physical_and_invitro.data$DTXSID[regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))!=-1]
#
# fit2 <- lm(data=MFdata.main,MFratio~MFratio.pred)
# summary(fit2)
# RMSE(fit2)
#
# fit2sub <- lm(data=subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred)
# summary(fit2sub)
#
# RMSE(fit2sub)
#
# cat(paste("When volatile and fluorinated chemicals are omitted only ",
# dim(subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))))[1],
# " evaluation chemicals remain, but the R2 is ",
# signif(summary(fit1sub)$adj.r.squared,2),
# " and the RMSE is ",
# signif(RMSE(fit1sub),2),
# " without the correction. When the fetal plasma fraction unbound correction is used, the predictivity decreases, slightly: R2 is ",
# signif(summary(fit2sub)$adj.r.squared,2),
# " and the RMSE is ",
# signif(RMSE(fit2sub),2),
# " for the non-volatile, non-fluorinated chemicals.\n",
# sep=""))
#
#
#
# cat(paste("We compare the RMSE for our predictions to the standard deviation \
# of the observations ",
# signif(sd(MFdata.main$MFratio)[1],2),
# " (",
# signif(sd(subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))$MFratio),2),
# " for non PAH or fluorinated compounds).\n",sep=""))
#
# cat(paste("The average standard deviation for chemicals with repeated observations was ",
# signif(sd(subset(MFdata.main,N.obs>1)$MFratio)[1],2),
# " (",
# signif(sd(subset(MFdata.main,
# N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))$MFratio),2),
# " for non PAH or fluorinated compounds).\n",sep=""))
#
# fit3 <- lm(data=repeats,MFratio~MFratio.pred.nofup)
# summary(fit3)
# RMSE(fit3)
#
# fit3sub <- lm(data=subset(MFdata.main, N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred.nofup)
# summary(fit3sub)
#
#
# fit4 <- lm(data=subset(MFdata.main,N.obs > 1),MFratio~MFratio.pred)
# summary(fit4)
# RMSE(fit4)
#
# fit4sub <- lm(data=subset(MFdata.main, N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))),
# MFratio~MFratio.pred)
# summary(fit4sub)
#
# cat(paste("Overall, we observed relatively poor correlation (R2 = ",
# signif(summary(fit1)$adj.r.squared,2),
# ", RMSE = ",
# signif(RMSE(fit1),2),
# ") without our fetal fup correction that was unchanged with that assumption (R2 = ",
# signif(summary(fit2)$adj.r.squared,2),
# ", RMSE = ",
# signif(RMSE(fit2),2),
# ").\n",sep=""))
#
# repeats <-subset(MFdata.main,N.obs > 1)
# cat(paste("The RMSE of the predictions for the ",
# dim(subset(repeats,!(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))))[1],
# " non-PAH and non-fluorinated compounds with repeated observations is ",
# signif(RMSE(fit4sub),2),
# " with the fup correction and ",
# signif(RMSE(fit3sub),2),
# " without.\n",sep=""))
#
# cat(paste("These values are close to the standard deviation of the mean but the p-value for the chemicals with repeated observations is ",
# signif(pf(
# summary(fit4sub)$fstatistic[1],
# summary(fit4sub)$fstatistic[2],
# summary(fit4sub)$fstatistic[3]),2),
# " indicating some value for the predictive model, albeit for only ",
# dim(subset(repeats,!(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons"))))[1],
# " chemicals",sep=""))
#
#
#
# Fig4.table <- data.frame()
# Fig4.table["Number of Chemicals","All Fig 4"] <- length(unique(MFdata.httk$DTXSID))
# Fig4.table["Number of Observations","All Fig 4"] <- dim(MFdata.httk)[1]
# Fig4.table["Observed Mean (Min - Max)","All Fig 4"] <- paste(
# signif(mean(MFdata.httk$MFratio),2)," (",
# signif(min(MFdata.httk$MFratio),2)," - ",
# signif(max(MFdata.httk$MFratio),2),")",sep="")
# Fig4.table["Observed Standard Deviation","All Fig 4"] <- signif(sd(MFdata.httk$MFratio),2)
# Fig4.table["Predicted Mean (Min - Max)","All Fig 4"] <- paste(
# signif(mean(MFdata.main$MFratio.pred,na.rm=TRUE),2)," (",
# signif(min(MFdata.main$MFratio.pred,na.rm=TRUE),2)," - ",
# signif(max(MFdata.main$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","All Fig 4"] <- signif(RMSE(fit2),2)
# Fig4.table["R-squared","All Fig 4"] <- signif(summary(fit2)[["adj.r.squared"]],2)
# Fig4.table["p-value","All Fig 4"] <- signif(summary(fit2)[["coefficients"]]["MFratio.pred",4],2)
#
# MFdata.sub1 <- subset(MFdata.httk,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))
#
# Fig4.table["Number of Chemicals","No PFAS/PAH"] <- length(unique(MFdata.sub1$DTXSID))
# Fig4.table["Number of Observations","No PFAS/PAH"] <- dim(MFdata.sub1)[1]
# Fig4.table["Observed Mean (Min - Max)","No PFAS/PAH"] <- paste(
# signif(mean(MFdata.sub1$MFratio,na.rm=TRUE),2)," (",
# signif(min(MFdata.sub1$MFratio,na.rm=TRUE),2)," - ",
# signif(max(MFdata.sub1$MFratio,na.rm=TRUE),2),")",sep="")
# Fig4.table["Observed Standard Deviation","No PFAS/PAH"] <- signif(sd(MFdata.sub1$MFratio),2)
#
# MFdata.sub2 <- subset(MFdata.main,
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))
#
# Fig4.table["Predicted Mean (Min - Max)","No PFAS/PAH"] <- paste(
# signif(mean(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," (",
# signif(min(MFdata.sub2$MFratio.pred,na.rm=TRUE),2)," - ",
# signif(max(MFdata.sub2$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","No PFAS/PAH"] <- signif(RMSE(fit2sub),2)
# Fig4.table["R-squared","No PFAS/PAH"] <- signif(summary(fit2sub)[["adj.r.squared"]],2)
# Fig4.table["p-value","No PFAS/PAH"] <- signif(summary(fit2sub)[["coefficients"]]["MFratio.pred",4],2)
#
#
# MFdata.sub3 <- subset(MFdata.main,N.obs > 1 &
# !(Chemical.Category %in% c(
# "Fluorinated compounds",
# "Polyaromatic Hydrocarbons")))
#
# Fig4.table["Number of Chemicals","Replicates Only"] <- length(unique(MFdata.sub3$DTXSID))
# Fig4.table["Number of Observations","Replicates Only"] <- dim(MFdata.sub3)[1]
# Fig4.table["Observed Mean (Min - Max)","Replicates Only"] <- paste(
# signif(mean(MFdata.sub3$MFratio),2)," (",
# signif(min(MFdata.sub3$MFratio),2)," - ",
# signif(max(MFdata.sub3$MFratio),2),")",sep="")
# Fig4.table["Observed Standard Deviation","Replicates Only"] <- signif(sd(MFdata.sub3$MFratio),2)
#
# Fig4.table["Predicted Mean (Min - Max)","Replicates Only"] <- paste(
# signif(mean(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," (",
# signif(min(MFdata.sub3$MFratio.pred,na.rm=TRUE),2)," - ",
# signif(max(MFdata.sub3$MFratio.pred,na.rm=TRUE),2),")",sep="")
# Fig4.table["RMSE","Replicates Only"] <- signif(RMSE(fit4sub),2)
# Fig4.table["R-squared","Replicates Only"] <- signif(summary(fit4sub)[["adj.r.squared"]],2)
# Fig4.table["p-value","Replicates Only"] <- signif(summary(fit4sub)[["coefficients"]]["MFratio.pred",4],2)
#
# write.csv(Fig4.table[1:6,],file="Fig4Table.csv")
## ----dallmann2018_data, eval = execute.vignette-------------------------------
# #TKstats <- as.data.frame(read_excel("MaternalFetalAUCData.xlsx"))
# TKstats <- httk::pregnonpregaucs
#
#
# ids <- unique(TKstats$DTXSID)
# cat(paste(sum(ids %in% get_cheminfo(
# model="fetal_pbtk",
# info="dtxsid",
# suppress.messages=TRUE)),
# "chemicals for which the model fetal_pbtk can run that are in the Dallmann data set."))
#
#
# TKstats.Cmax <- subset(TKstats,DTXSID!="" & Parameter=="Cmax")
# TKstats <- subset(TKstats,DTXSID!="" & Parameter %in% c("AUCinf","AUClast"))
#
# ids <- unique(TKstats$DTXSID)
# cat(paste(sum(ids %in% get_cheminfo(
# model="fetal_pbtk",
# info="dtxsid",
# suppress.messages=TRUE)),
# "chemicals for which the model fetal_pbtk can run that have AUCs in the Dallmann data set."))
#
#
# # Assumed body weight (kg) from Kapraun 2019
# BW.nonpreg <- 61.103
#
# #TKstats$Dose <- TKstats$Dose/BW
# #TKstats$Dose.Units <- "mg/kg"
# colnames(TKstats)[colnames(TKstats)=="Observed Pregnant"] <- "Observed.Pregnant"
# colnames(TKstats)[colnames(TKstats)=="Observed Pregnant5"] <- "Observed.Pregnant5"
# colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant"] <- "Observed.Non.Pregnant"
# colnames(TKstats)[colnames(TKstats)=="Observed Non-Pregnant2"] <- "Observed.Non.Pregnant2"
# colnames(TKstats)[colnames(TKstats)=="Dose Route"] <- "Dose.Route"
## ----dallmann2018_preds, eval = execute.vignette------------------------------
# for (this.id in unique(TKstats$DTXSID))
# {
# if (any(regexpr("ng",TKstats[TKstats$DTXSID==this.id,"Dose Units"])!=-1))
# {
# }
# if (this.id %in% get_cheminfo(
# info="DTXSID",
# model="fetal_pbtk",
# suppress.messages=TRUE))
# {
# this.subset <- subset(TKstats,DTXSID==this.id)
# p <- parameterize_pbtk(
# dtxsid=this.id,
# suppress.messages=TRUE)
# p$hematocrit <- 0.39412 # Kapraun 2019 (unitless)
# p$Rblood2plasma <- calc_rblood2plasma(
# parameters=p,
# suppress.messages=TRUE)
# p$BW <- BW.nonpreg
# p$Qcardiacc <- 301.78 / p$BW^(3/4) # Kapraun 2019 (L/h/kg^3/4)
# for (this.route in unique(this.subset$Dose.Route))
# {
# this.route.subset <- subset(this.subset, Dose.Route==this.route)
# if (this.route.subset[1,"Gestational.Age.Weeks"] > 12)
# {
# this.dose <- this.route.subset$Dose
# # Non-pregnant PBPK:
# out.nonpreg <- solve_pbtk(
# parameters=p,
# times = seq(0, this.route.subset[1,"NonPreg.Duration.Days"],
# this.route.subset[1,"NonPreg.Duration.Days"]/100),
# dose=this.dose/BW.nonpreg, # mg/kg
# daily.dose=NULL,
# iv.dose=(this.route=="iv"),
# output.units="uM",
# suppress.messages=TRUE)
# # Pregnant PBPK:
# BW.preg <- suppressWarnings(calc_maternal_bw(
# week=this.route.subset[1,"Gestational.Age.Weeks"]))
# out.preg <- solve_fetal_pbtk(
# dtxsid=this.id,
# times = seq(
# this.route.subset[1,"Gestational.Age.Weeks"]*7,
# this.route.subset[1,"Gestational.Age.Weeks"]*7 +
# this.route.subset[1,"Preg.Duration.Days"],
# this.route.subset[1,"Preg.Duration.Days"]/100),
# dose=this.dose/BW.preg, # mg/kg
# iv.dose=(this.route=="iv"),
# daily.dose=NULL,
# output.units="uM",
# suppress.messages=TRUE)
# if (any(regexpr("AUC",this.subset$Parameter)!=-1))
# {
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("AUC",TKstats$Parameter)!=-1,
# "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"AUC"])*24 #uMdays->uMh
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("AUC",TKstats$Parameter)!=-1,
# "Predicted.Pregnant.httk"] <- max(out.preg[,"AUC"])*24 # uM days -> uM h
# }
# if (any(regexpr("Cmax",this.subset$Parameter)!=-1))
# {
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("Cmax",TKstats$Parameter)!=-1,
# "Predicted.Non.Pregnant.httk"] <- max(out.nonpreg[,"Cplasma"])
# TKstats[TKstats$DTXSID==this.id &
# TKstats$Dose.Route == this.route &
# regexpr("Cmax",TKstats$Parameter)!=-1,
# "Predicted.Pregnant.httk"] <- max(out.preg[,"Cfplasma"])
# }
# }
# }
# }
# }
#
# TKstats$Ratio.obs <- TKstats$Observed.Pregnant / TKstats$Observed.Non.Pregnant
# TKstats$Ratio.httk <- TKstats$Predicted.Pregnant.httk / TKstats$Predicted.Non.Pregnant.httk
# TKstats <- subset(TKstats, !is.na(TKstats$Ratio.httk))
#
# write.csv(TKstats,file="Table-Dallmann2018.csv",row.names=FALSE)
## ----dallmann2018_figure, eval = execute.vignette-----------------------------
# FigEa <- ggplot(data=TKstats) +
# geom_point(aes(
# y=Observed.Non.Pregnant2,
# x=Predicted.Non.Pregnant.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# scale_shape_manual(values=c(22,24))+
# xlab("httk Predicted (uM*h)") +
# ylab("Observed") +
# scale_x_log10(#limits=c(10^-3,10^3),
# label=scientific_10) +
# scale_y_log10(#limits=c(10^-3,10^3),
# label=scientific_10) +
# ggtitle("Non-Pregnant") +
# theme_bw() +
# theme(legend.position="none")+
# theme( text = element_text(size=12)) +
# annotate("text", x = 0.1, y = 1000, label = "A", fontface =2)
#
# #print(FigEa)
# cat(paste(
# "For the Dallman et al. (2018) data we observe an average-fold error (AFE)\n\ of",
# signif(mean(TKstats$Predicted.Non.Pregnant.httk/TKstats$Observed.Non.Pregnant2),2),
# "and a RMSE (log10-scale) of",
# signif((mean((log10(TKstats$Predicted.Non.Pregnant.httk)-log10(TKstats$Observed.Non.Pregnant2))^2))^(1/2),2),
# "for non-pregnant woman.\n"))
#
# FigEb <- ggplot(data=TKstats) +
# geom_point(aes(
# y=Observed.Pregnant5,
# x=Predicted.Pregnant.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# scale_shape_manual(values=c(22,24))+
# xlab("httk Predicted (uM*h)") +
# ylab("Observed") +
# scale_x_log10(#limits=c(10^-5,10^3),
# label=scientific_10) +
# scale_y_log10(#limits=c(10^-5,10^3),
# label=scientific_10) +
# ggtitle("Pregnant")+
# theme_bw() +
# theme(legend.position="none")+
# theme( text = element_text(size=12)) +
# annotate("text", x = 0.1, y = 300, label = "B", fontface =2)
#
# #print(FigEb)
# cat(paste(
# "We observe an AFE of",
# signif(mean(TKstats$Predicted.Pregnant.httk/TKstats$Observed.Pregnant5),2),
# "and a RMSE (log10-scale) of",
# signif((mean((log10(TKstats$Predicted.Pregnant.httk)-log10(TKstats$Observed.Pregnant5))^2))^(1/2),2),
# "for pregnant woman.\n"))
#
#
#
# FigEc <- ggplot(data=TKstats) +
# geom_point(aes(
# x=Predicted.Non.Pregnant.httk,
# y=Predicted.Pregnant.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# scale_shape_manual(values=c(22,24),name="Route")+
# ylab("httk Pregnant (uM*h)") +
# xlab("httk Non-Pregnant (uM*h)") +
# scale_x_log10(#limits=c(10^-1,10^2),
# label=scientific_10) +
# scale_y_log10(#limits=c(10^-1,10^2),
# label=scientific_10) +
# ggtitle("Model Comparison") +
# theme_bw() +
# theme(legend.position="left")+
# theme( text = element_text(size=12))+
# theme(legend.text=element_text(size=10))+
# guides(fill=guide_legend(ncol=2,override.aes=list(shape=21)),alpha=guide_legend(ncol=2),shape=guide_legend(ncol=2))
# print(FigEc)
# FigEc <- get_legend(FigEc)
#
#
# FigEd <- ggplot(data=subset(TKstats,!is.na(Ratio.httk))) +
# geom_point(aes(
# y=Ratio.obs,
# x=Ratio.httk,
# shape=Dose.Route,
# alpha=Drug,
# fill=Drug),
# size=5) +
# scale_shape_manual(values=c(22,24))+
# xlab("httk Predicted") +
# ylab("Observed") +
# scale_x_continuous(limits=c(0.25,3)) +
# scale_y_continuous(limits=c(0.25,3)) +
# geom_abline(slope=1, intercept=0) +
# geom_abline(slope=1, intercept=1, linetype=3) +
# geom_abline(slope=1, intercept=-1, linetype=3) +
# ggtitle("Ratio") +
# theme_bw() +
# theme(legend.position="none")+
# theme( text = element_text(size=12)) +
# annotate("text", x = 0.5, y = 2.75, label = "C", fontface =2)
#
# dev.new()
# grid.arrange(FigEa,FigEb,FigEc,FigEd,nrow=2)
#
# write.csv(subset(TKstats,
# !is.na(Ratio.httk)),
# file="DallmanTable.txt")
#
## ----curley1969_data, eval = execute.vignette---------------------------------
# Curley <- as.data.frame(read_excel("Curley1969.xlsx"))
# dim(Curley)
# Curley.compounds <- Curley[1,4:13]
# Curley <- Curley[4:47,]
# colnames(Curley)[1] <- "Tissue"
# colnames(Curley)[2] <- "N"
# colnames(Curley)[3] <- "Stat"
## ----curley1969_calcpcs, eval = execute.vignette------------------------------
# Curley.pcs <- NULL
# cord.blood <- subset(Curley, Tissue == "Cord Blood")
# suppressWarnings(
# for (this.tissue in unique(Curley$Tissue))
# if (this.tissue != "Cord Blood")
# {
# this.subset <- subset(Curley, Tissue == this.tissue)
# for (this.chemical in colnames(Curley)[4:13])
# {
# if (!is.na(as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical])) &
# !is.na(as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])))
# {
# this.row <- data.frame(
# Compound = Curley.compounds[,this.chemical],
# DTXSID = this.chemical,
# Tissue = this.tissue,
# PC = as.numeric(subset(this.subset,Stat=="Mean")[,this.chemical]) /
# as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
# )
# Curley.pcs <- rbind(Curley.pcs,this.row)
# } else if (!is.na((as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]))) &
# !is.na((as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical]))))
# {
# this.row <- data.frame(
# Compound = Curley.compounds[,this.chemical],
# DTXSID = this.chemical,
# Tissue = this.tissue,
# PC = as.numeric(subset(this.subset,Stat=="Range")[,this.chemical]) /
# as.numeric(subset(cord.blood,Stat=="Mean")[,this.chemical])
# )
# Curley.pcs <- rbind(Curley.pcs,this.row)
# }
# }
# }
# )
# Curley.pcs[Curley.pcs$Tissue=="Lungs","Tissue"] <- "Lung"
#
# pearce2017 <- read_excel("PC_Data.xlsx")
# pearce2017 <- subset(pearce2017,DTXSID%in%Curley.pcs$DTXSID)
# any(pearce2017$DTXSID%in%Curley.pcs$DTXSID)
# print("None of the Curley chems were in the Pearce 2017 PC predictor evaluation.")
## ----csanady2002_pcs, eval = execute.vignette---------------------------------
# csanadybpa <- read_excel("Csanady2002.xlsx",sheet="Table 2")
# csanadybpa$Compound <- "Bisphenol A"
# csanadybpa$DTXSID <- "DTXSID7020182"
# csanadybpa <- csanadybpa[,c("Compound","DTXSID","Tissue","Mean")]
# colnames(csanadybpa) <- colnames(Curley.pcs)
#
# csanadydaid <- read_excel("Csanady2002.xlsx",sheet="Table 3",skip=1)
# csanadydaid$Compound <- "Daidzein"
# csanadydaid$DTXSID <- "DTXSID9022310"
# csanadydaid <- csanadydaid[,c("Compound","DTXSID","Tissue","Mean...6")]
# colnames(csanadydaid) <- colnames(Curley.pcs)
#
# Curley.pcs <- rbind(Curley.pcs,csanadybpa,csanadydaid)
# Curley.pcs <- subset(Curley.pcs,Tissue!="Mammary gland")
## ----weijs2013_loadPCs, eval = execute.vignette-------------------------------
# weijstab3 <- read_excel("Weijs2013.xlsx",sheet="Table3",skip=1)
# weijstab3 <- subset(weijstab3, !is.na(Compound) & !is.na(Tissue))
# weijstab4 <- read_excel("Weijs2013.xlsx",sheet="Table4",skip=1)
# weijstab4 <- subset(weijstab4, !is.na(Compound) & !is.na(Tissue))
#
# for (this.compound in unique(weijstab3$Compound))
# for (this.tissue in unique(weijstab3$Tissue))
# {
# Curley.pcs[
# Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
# "Weijs2013"] <- weijstab3[weijstab3$Compound==this.compound &
# weijstab3$Tissue==this.tissue,"value"]
#
# }
#
# for (this.compound in unique(weijstab4$Compound))
# for (this.tissue in unique(weijstab4$Tissue))
# {
# Curley.pcs[
# Curley.pcs$DTXSID==this.compound & Curley.pcs$Tissue==this.tissue,
# "Weijs2013"] <- weijstab4[weijstab4$Compound==this.compound &
# weijstab4$Tissue==this.tissue,"value"]
#
# }
#
# write.csv(Curley.pcs,"PCs-table.csv",row.names=FALSE)
#
## ----curley1969_makepreds, eval = execute.vignette----------------------------
# Curley.pcs <- httk::fetalpcs
# dtxsid.list <- get_cheminfo(
# info="DTXSID",
# model="fetal_pbtk",
# suppress.messages=TRUE)
#
# suppressWarnings(load_sipes2017())
# for (this.chemical in unique(Curley.pcs$DTXSID))
# if (this.chemical %in% dtxsid.list)
# {
# this.subset <- subset(Curley.pcs,DTXSID==this.chemical)
# p <- parameterize_fetal_pbtk(dtxsid=this.chemical,
# fetal_fup_adjustment = FALSE,
# suppress.messages=TRUE,
# minimum.Funbound.plasma = 1e-16)
# fetal.blood.pH <- 7.26
# Fup <- p$Fraction_unbound_plasma_fetus
# fetal_schmitt_parms <- parameterize_schmitt(
# dtxsid=this.chemical,
# suppress.messages=TRUE,
# minimum.Funbound.plasma = 1e-16)
# fetal_schmitt_parms$plasma.pH <- fetal.blood.pH
# fetal_schmitt_parms$Funbound.plasma <- Fup
# Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Fup"] <- Fup
# # httk gives tissue:fup (unbound chemical in plasma) PC's:
# fetal_pcs <- predict_partitioning_schmitt(
# parameters = fetal_schmitt_parms,
# suppress.messages=TRUE,
# model="fetal_pbtk",
# minimum.Funbound.plasma = 1e-16)
# fetal_pcs.nocal <- predict_partitioning_schmitt(
# parameters = fetal_schmitt_parms,
# regression=FALSE,
# suppress.messages=TRUE,
# model="fetal_pbtk",
# minimum.Funbound.plasma = 1e-16)
# out <- solve_fetal_pbtk(
# dtxsid = this.chemical,
# fetal_fup_adjustment =FALSE,
# suppress.messages=TRUE,
# minimum.Funbound.plasma = 1e-16)
# Rb2p <- out[dim(out)[1],"Rfblood2plasma"]
# Curley.pcs[Curley.pcs$DTXSID==this.chemical,"Rb2p"] <- Rb2p
# # Convert to tissue:blood PC's:
# for (this.tissue in this.subset$Tissue)
# if (tolower(this.tissue) %in%
# unique(subset(tissue.data,Species=="Human")$Tissue))
# {
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred.t2up"] <-
# fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal.t2up"] <-
# fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred"] <-
# fetal_pcs[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
# Curley.pcs[Curley.pcs$DTXSID==this.chemical &
# Curley.pcs$Tissue == this.tissue, "HTTK.pred.nocal"] <-
# fetal_pcs.nocal[[paste("K",tolower(this.tissue),"2pu",sep="")]]*Fup/Rb2p
# } else {
# print(this.tissue)
# }
# } else print(this.chemical)
# reset_httk()
#
# cat(paste(
# "For the two placental observations (",
# signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"PC"],2),
# " for Bisphenol A and ",
# signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"PC"],2),
# " for Diadzen) the predictions were ",
# signif(subset(Curley.pcs,Compound=="Bisphenol A" & Tissue=="Placenta")[,"HTTK.pred"],2),
# " and ",
# signif(subset(Curley.pcs,Compound=="Daidzein" & Tissue=="Placenta")[,"HTTK.pred"],2),
# ", respectively.\n",sep=""))
#
#
## ----curley1969_figure_compounds, eval = execute.vignette---------------------
# FigFa <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred.nocal))) +
# geom_point(size=2,aes(
# y=Weijs2013,
# x=HTTK.pred,
# shape=Compound,
# color=Compound)) +
# geom_point(size=5,aes(
# y=PC,
# x=HTTK.pred,
# shape=Compound,
# color=Compound)) +
# geom_abline(slope=1, intercept=0, size=2) +
# geom_abline(slope=1, intercept=1, linetype=3, size=2) +
# geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
# scale_shape_manual(values=rep(c(15,16,17,18),4))+
# xlab("Predicted Tissue:Blood Partition Coefificent") +
# ylab("Observed") +
# scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
# scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=28))+
# theme(legend.text=element_text(size=12))+
# theme(legend.title = element_blank())+
# guides(shape=guide_legend(nrow=3,byrow=TRUE))
#
# print(FigFa)
#
# fitFa <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred))
# RMSE(fitFa)
# fitFb <- lm(data=Curley.pcs,log10(PC)~log10(HTTK.pred.nocal))
# RMSE(fitFb)
## ----curley1969_figure_tissues, eval = execute.vignette-----------------------
# FigFb <- ggplot(data=subset(Curley.pcs,!is.na(HTTK.pred))) +
# geom_point(size=2,aes(
# y=Weijs2013,
# x=HTTK.pred,
# shape=Tissue,
# color=Tissue)) +
# geom_point(size=5, aes(
# y=PC,
# x=HTTK.pred,
# shape=Tissue,
# color=Tissue)) +
# geom_abline(slope=1, intercept=0, size=2) +
# geom_abline(slope=1, intercept=1, linetype=3, size=2) +
# geom_abline(slope=1, intercept=-1, linetype=3, size=2) +
# scale_shape_manual(values=rep(c(15,16,17,18),4))+
# xlab("Predicted Tissue:Blood Partition Coefificent") +
# ylab("Observed") +
# scale_x_log10(label=scientific_10,limits=c(10^-1,10^2)) +
# scale_y_log10(label=scientific_10,limits=c(10^-1,10^4)) +
# theme_bw() +
# theme(legend.position="bottom")+
# theme( text = element_text(size=28))+
# theme(legend.text=element_text(size=20))+
# theme(legend.title = element_blank())+
# guides(shape=guide_legend(nrow=3,byrow=TRUE))
#
# print(FigFb)
## ----read_pksim_pcs, eval = execute.vignette----------------------------------
# pksim.pcs <- as.data.frame(read_excel("PartitionCoefficients.xlsx"))
# dim(pksim.pcs)
# for (this.id in unique(pksim.pcs$DTXSID))
# {
# httk.PCs <- predict_partitioning_schmitt(dtxsid=this.id,
# suppress.messages = TRUE)
# p <-
# parameterize_fetal_pbtk(dtxsid=this.id,
# suppress.messages = TRUE)
# httk.PCs[["Kplacenta2pu"]] <- p[["Kplacenta2pu"]]
# httk.PCs[["Fup"]] <- p[["Funbound.plasma"]]
# lapply(httk.PCs,as.numeric)
# this.subset <- subset(pksim.pcs,DTXSID==this.id)
# for (this.tissue in unique(this.subset$Tissue))
# {
# if (any(regexpr(tolower(this.tissue),names(httk.PCs))!=-1))
# {
# pksim.pcs[pksim.pcs$DTXSID==this.id & pksim.pcs$Tissue==this.tissue,
# "HTTK.pred"] <-
# httk.PCs[regexpr(tolower(this.tissue),names(httk.PCs))!=-1][[1]]*
# httk.PCs[["Fup"]][[1]]
# }
# }
# }
# colnames(pksim.pcs)[colnames(pksim.pcs) ==
# "Tissue-to-plasma partition coefficient"] <- "PKSim.pred"
# colnames(pksim.pcs)[colnames(pksim.pcs) ==
# "Method for calculating tissue-to-plasma partition coefficients"] <-
# "Method"
# pksim.pcs[,"Method"] <- as.factor(pksim.pcs[,"Method"])
# pksim.pcs[,"Drug"] <- as.factor(pksim.pcs[,"Drug"])
# pksim.pcs[,"Tissue"] <- as.factor(pksim.pcs[,"Tissue"])
#
#
# pksim.pcs[,"PKSim.pred"] <- as.numeric(
# pksim.pcs[,"PKSim.pred"])
#
## ----compare_pksim_pcs, eval = execute.vignette-------------------------------
# pksim.pcs <- httk::pksim.pcs
#
# pksim.fit1 <- lm(data=pksim.pcs,
# log10(PKSim.pred)~log10(HTTK.pred))
# summary(pksim.fit1)
#
# pksim.fit2 <- lm(data=pksim.pcs,
# log10(PKSim.pred)~log10(HTTK.pred)+
# Drug+Tissue+Method)
# summary(pksim.fit2)
#
## ----wang2018_loaddata, eval = execute.vignette-------------------------------
# #Wangchems <- read_excel("Wang2018.xlsx",sheet=3,skip=2)
# Wangchems <- httk::wang2018
# maternal.list <- Wangchems$CASRN[Wangchems$CASRN %in%
# get_cheminfo(model="fetal_pbtk",
# suppress.messages=TRUE)]
# nonvols <- subset(chem.physical_and_invitro.data,logHenry < -4.5)$CAS
# nonfluoros <- chem.physical_and_invitro.data$CAS[
# regexpr("fluoro",tolower(chem.physical_and_invitro.data$Compound))==-1]
# maternal.list <- maternal.list[maternal.list %in% intersect(nonvols,nonfluoros)]
## ----wang2018_makepreds, eval = execute.vignette------------------------------
#
# pred.table1 <- subset(get_cheminfo(
# info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
# model="fetal_pbtk"),
# CAS %in% maternal.list,
# suppress.messages=TRUE)
# pred.table1$Compound <- gsub("\"","",pred.table1$Compound)
#
# for (this.cas in maternal.list)
# {
# Fup <- subset(get_cheminfo(
# info="all",
# suppress.messages=TRUE),
# CAS==this.cas)$Human.Funbound.plasma
# # Check if Fup is provided as a distribution:
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
# {
# skip <- TRUE
# } else skip <- FALSE
# Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
# Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
# # These results are actually from a Bayesian posterior, but we can approximate that
# # they are normally distributed about the median to estimate a standard deviation:
# Fup.sigma <- (Fup.upper - Fup)/1.96
# # If it's not a distribution use defaults (see Wambaugh et al. 2019)
# } else if (as.numeric(Fup)== 0)
# {
# skip <- TRUE
# } else {
# skip <- FALSE
# Fup <- as.numeric(Fup)
# Fup.sigma <- Fup*0.4
# Fup.upper <- min(1,(1+1.96*0.4)*Fup)
# }
#
# # Only run the simulation if we have the necessary parameters:
# if (!skip)
# {
# print(this.cas)
# out <- solve_fetal_pbtk(
# chem.cas=this.cas,
# dose=0,
# daily.dose=1,
# doses.per.day=3,
# fetal_fup_adjustenment=TRUE,
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# pred.table1[pred.table1$CAS==this.cas,"fup"] <- signif(Fup,3)
# pred.table1[pred.table1$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
# pred.table1[pred.table1$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
# pred.table1[pred.table1$CAS==this.cas,"Ratio.pred"] <-
# signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
# }
# }
## ----brain_uncertainty_prop_table, eval = execute.vignette--------------------
# PC.table <- pred.table1
# colnames(PC.table)[colnames(PC.table)=="Ratio.pred"] <- "R.plasma.FtoM"
# pred.table1$Uncertainty <- "Predicted F:M Plasma Ratio"
## ----Rfetmat_uncertainty, eval = execute.vignette-----------------------------
# pred.table3 <- pred.table1
# pred.table3$Uncertainty <- "Plasma Error (Fig. 4)"
# empirical.error <- RMSE(fit2sub)
# for (this.cas in maternal.list)
# {
#
# pred.table3[pred.table3$CAS==this.cas,"Ratio.pred"] <- signif(
# 1/((1/pred.table1[pred.table3$CAS==this.cas,"Ratio.pred"])*
# (1-1.96*empirical.error)), 3)
# }
# # Update final table for paper:
# PC.table$RMSE <- signif(RMSE(fit2sub),3)
# PC.table$R.plasma.FtoM.upper <- pred.table3$Ratio.pred
## ----Rfetbrain, eval = execute.vignette---------------------------------------
# pred.table4 <- pred.table1
# pred.table4$Uncertainty <- "Fetal Brain Partitioning"
# for (this.cas in maternal.list)
# {
# print(this.cas)
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=TRUE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
# fup <- pred.table4[pred.table4$CAS==this.cas,"fup"]
# pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"] <- signif(
# PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]*
# Kbrain2pu * fup, 3)
# PC.table[PC.table$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
# PC.table[PC.table$CAS==this.cas,"R.brain.FtoM"] <-
# pred.table4[pred.table4$CAS==this.cas,"Ratio.pred"]
# }
## ----Rfetbrain_uncertainty, eval = execute.vignette---------------------------
#
# pred.table5 <- pred.table1
# pred.table5$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
# for (this.cas in maternal.list)
# {
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=TRUE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
#
# fup <- pred.table5[pred.table5$CAS==this.cas,"fup"]
# sigma.fup <- pred.table5[pred.table5$CAS==this.cas,"fup.sigma"]
# Rmatfet <- 1/PC.table[PC.table$CAS==this.cas,"R.plasma.FtoM"]
# Rbrain2matblood <- Kbrain2pu * fup / Rmatfet
#
# # From Pearce et al. (2017) PC paper:
# sigma.logKbrain <- 0.647
# Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)
#
# error.Rmatfet <- PC.table[PC.table$CAS==this.cas,"RMSE"]
# sigma.Rbrain2matblood <- Rbrain2matblood *
# (log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2)
# Rbrain2matblood.upper <- Rbrain2matblood *
# (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2))
# pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"] <-
# signif(Rbrain2matblood.upper,3)
# PC.table[PC.table$CAS==this.cas,"sigma.logKbrain"] <-
# signif(sigma.logKbrain, 3)
# PC.table[PC.table$CAS==this.cas,"Kbrain2pu.upper"] <-
# signif(Kbrain2pu.upper, 3)
# PC.table[PC.table$CAS==this.cas,"sigma.Rbrain2matblood"] <-
# signif(sigma.Rbrain2matblood, 3)
# PC.table[PC.table$CAS==this.cas,"R.brain.FtoM.upper"] <-
# pred.table5[pred.table5$CAS==this.cas,"Ratio.pred"]
# }
## ----Wang_Uncertainty_Propagation, eval = execute.vignette--------------------
#
# pred.levels <- pred.table5$Compound[order(pred.table5$Ratio.pred)]
#
# pred.table <- rbind(
# pred.table1,
# # pred.table2,
# pred.table3,
# pred.table4,
# pred.table5)
# pred.table$Compound <- factor(pred.table$Compound,
# levels = pred.levels)
#
# pred.table$Uncertainty <- factor(pred.table$Uncertainty,
# levels = c(pred.table1[1,"Uncertainty"],
# # pred.table2[1,"Uncertainty"],
# pred.table3[1,"Uncertainty"],
# pred.table4[1,"Uncertainty"],
# pred.table5[1,"Uncertainty"]))
## ----wang2018_figure, eval = execute.vignette---------------------------------
# #Wang 2018 confirmed 6 chemical identities:
# confirmed.chemicals <- c(
# "2,4-Di-tert-butylphenol",
# "2,4-Dinitrophenol",
# "Pyrocatechol",
# "2'-Hydroxyacetophenone",
# "3,5-Di-tert-butylsalicylic acid",
# "4-Hydroxycoumarin"
# )
# confirmed.chemicals <- c(
# "96-76-4",
# "19715-19-6",
# "51-28-5",
# "120-80-9",
# "118-93-4",
# "1076-38-6")
#
#
# FigG <- ggplot(data=pred.table) +
# geom_point(aes(
# x=Ratio.pred,
# y=Compound,
# color = Uncertainty,
# shape = Uncertainty),
# size=3) +
# scale_shape_manual(values=c(15, 16,2, 23, 0, 1, 17, 5, 6))+
# scale_x_log10(limits=c(10^-2,10^3),label=scientific_10)+
# ylab(expression(paste(
# "Chemicals Found in Maternal Plasma by Wang ",italic("et al.")," (2018)"))) +
# xlab("Predicted Ratio to Maternal Plasma") +
# theme_bw() +
# # theme(legend.position="bottom")+
# theme( text = element_text(size=14))+
# theme(legend.text=element_text(size=10))#+
# # guides(color=guide_legend(nrow=4,byrow=TRUE))+
# #guides(shape=guide_legend(nrow=4,byrow=TRUE))
# #+
# # theme(legend.justification = c(0, 0), legend.position = c(0, 0))
#
# print(FigG)
## ----wang2018_details, eval = execute.vignette--------------------------------
#
# for (this.col in 7:14)
# PC.table[,this.col] <- signif(PC.table[,this.col],3)
#
# PC.table <- PC.table[order(PC.table$R.brain.FtoM.upper,decreasing=TRUE),]
#
# for (this.row in 1:dim(PC.table)[1])
# {
# out <- calc_ionization(
# pH=7.26,
# pKa_Donor=PC.table[this.row,"pKa_Donor"],
# pKa_Accept=PC.table[this.row,"pKa_Accept"])
# if (out$fraction_neutral>0.9) PC.table[this.row,"Charge_726"] <- "Neutral"
# else if (out$fraction_positive>0.1) PC.table[this.row,"Charge_726"] <-
# paste(signif(out$fraction_positive*100,2),"% Positive",sep="")
# else if (out$fraction_negative>0.1) PC.table[this.row,"Charge_726"] <-
# paste(signif(out$fraction_negative*100,2),"% Negative",sep="")
# else if (out$fraction_zwitter>0.1) PC.table[this.row,"Charge_726"] <-
# paste(signif(out$fraction_zwitter*100,2),"% Zwitterion",sep="")
# }
#
# PC.table <- PC.table[,c(
# "Compound",
# "CAS",
# "DTXSID",
# "logP",
# "Charge_726",
# "R.plasma.FtoM",
# "RMSE",
# "R.plasma.FtoM.upper",
# "Kbrain2pu",
# "fup",
# "R.brain.FtoM",
# "Kbrain2pu.upper",
# "R.brain.FtoM.upper")]
#
# write.csv(PC.table,
# file="WangTable.txt",
# row.names=FALSE)
## ----impact_of_model, eval = execute.vignette---------------------------------
# MFbrainfit <- lm(R.brain.FtoM.upper~Kbrain2pu.upper,data=PC.table)
# summary(MFbrainfit)
#
# cat(paste("As expected, the predicted fetal brain to maternal blood ratio was strongly correlated with the predicted brain partition coefficient (R2 = ",
# signif(summary(MFbrainfit)$adj.r.square,2),
# ", p-value ",
# signif(summary(MFbrainfit)$coefficients[2,4],2),
# "). However, the PBTK simulation impacted the plasma and tissue concentrations such that ",
# dim(subset(PC.table, rank(R.brain.FtoM.upper) < rank(Kbrain2pu.upper)))[1],
# " chemicals were ranked higher than they would have been using partitioning alone.",
# sep=""))
## ----wang2018_makepreds_nofup, eval = execute.vignette------------------------
#
# pred.table1.nofup <- subset(get_cheminfo(
# info=c("Compound","CAS","DTXSID","logP","pka_accept","pka_donor"),
# model="fetal_pbtk",
# suppress.messages=TRUE),
# CAS %in% maternal.list)
# pred.table1.nofup$Compound <- gsub("\"","",pred.table1.nofup$Compound)
#
# for (this.cas in maternal.list)
# {
# Fup <- subset(get_cheminfo(
# info="all",
# suppress.messages=TRUE),
# CAS==this.cas)$Human.Funbound.plasma
# # Check if Fup is provided as a distribution:
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.11))
# {
# skip <- TRUE
# } else skip <- FALSE
# Fup.upper <- as.numeric(strsplit(Fup,",")[[1]][3])
# Fup <- as.numeric(strsplit(Fup,",")[[1]][1])
# # These results are actually from a Bayesian posterior, but we can approximate that
# # they are normally distributed about the median to estimate a standard deviation:
# Fup.sigma <- (Fup.upper - Fup)/1.96
# # If it's not a distribution use defaults (see Wambaugh et al. 2019)
# } else if (as.numeric(Fup)== 0)
# {
# skip <- TRUE
# } else {
# skip <- FALSE
# Fup <- as.numeric(Fup)
# Fup.sigma <- Fup*0.4
# Fup.upper <- min(1,(1+1.96*0.4)*Fup)
# }
#
# # Only run the simulation if we have the necessary parameters:
# if (!skip)
# {
# out <- solve_fetal_pbtk(
# chem.cas=this.cas,
# dose=0,
# daily.dose=1,
# doses.per.day=3,
# fetal_fup_adjustenment=FALSE,
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup"] <- signif(Fup,3)
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.sigma"] <- signif(Fup.sigma,3)
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"fup.upper"] <- signif(Fup.upper,3)
# pred.table1.nofup[pred.table1.nofup$CAS==this.cas,"Ratio.pred"] <-
# signif(mean(out[last.row,"Cfplasma"])/mean(out[last.row,"Cplasma"]),3)
# }
# }
## ----brain_uncertainty_prop_table_nofup, eval = execute.vignette--------------
# PC.table.nofup <- pred.table1.nofup
# colnames(PC.table.nofup)[colnames(PC.table.nofup)=="Ratio.pred"] <- "R.plasma.FtoM"
# pred.table1.nofup$Uncertainty <- "Predicted F:M Plasma Ratio"
## ----Rfetmat_uncertainty_nofup, eval = execute.vignette-----------------------
# pred.table3.nofup <- pred.table1.nofup
# pred.table3.nofup$Uncertainty <- "Plasma Error (Fig. 4)"
# empirical.error <- RMSE(fit2sub)
# for (this.cas in maternal.list)
# {
#
# pred.table3.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
# 1/((1/pred.table1.nofup[pred.table3.nofup$CAS==this.cas,"Ratio.pred"])*
# (1-1.96*empirical.error)), 3)
# }
# # Update final table for paper:
# PC.table.nofup$RMSE <- signif(RMSE(fit2sub),3)
# PC.table.nofup$R.plasma.FtoM.upper <- pred.table3.nofup$Ratio.pred
## ----Rfetbrain_nofup, eval = execute.vignette---------------------------------
# pred.table4.nofup <- pred.table1.nofup
# pred.table4.nofup$Uncertainty <- "Fetal Brain Partitioning"
# for (this.cas in maternal.list)
# {
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=FALSE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
# fup <- pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"fup"]
# pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"] <- signif(
# PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]*
# Kbrain2pu * fup, 3)
# PC.table.nofup[PC.table.nofup$CAS==this.cas,"Kbrain2pu"] <- Kbrain2pu
# PC.table.nofup[PC.table.nofup$CAS==this.cas,"R.brain.FtoM"] <-
# pred.table4.nofup[pred.table4.nofup$CAS==this.cas,"Ratio.pred"]
# }
## ----Rfetbrain_uncertainty_nofup, eval = execute.vignette---------------------
#
# pred.table5.nofup <- pred.table1.nofup
# pred.table5.nofup$Uncertainty <- "Brain Partitioning Error (Pearce et al., 2017)"
# for (this.cas in maternal.list)
# {
# p <- parameterize_fetal_pbtk(chem.cas=this.cas,
# fetal_fup_adjustment=FALSE,
# suppress.messages=TRUE)
# Kbrain2pu <- p$Kfbrain2pu
#
# fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup"]
# sigma.fup <- pred.table5.nofup[pred.table5.nofup$CAS==this.cas,"fup.sigma"]
# Rmatfet <- 1/PC.table[PC.table.nofup$CAS==this.cas,"R.plasma.FtoM"]
# Rbrain2matblood <- Kbrain2pu * fup / Rmatfet
#
# # From Pearce et al. (2017) PC paper:
# sigma.logKbrain <- 0.647
# Kbrain2pu.upper <- signif(10^(log10(Kbrain2pu)+1.96*sigma.logKbrain),3)
#
# error.Rmatfet <- PC.table.nofup [PC.table.nofup $CAS==this.cas,"RMSE"]
# sigma.Rbrain2matblood <- Rbrain2matblood *
# (log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2)
# Rbrain2matblood.upper <- Rbrain2matblood *
# (1 + 1.96*(log(10)^2*sigma.logKbrain^2 +
# sigma.fup^2/fup^2 +
# error.Rmatfet^2/Rmatfet^2)^(1/2))
# pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"] <-
# signif(Rbrain2matblood.upper,3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.logKbrain"] <-
# signif(sigma.logKbrain, 3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"Kbrain2pu.upper"] <-
# signif(Kbrain2pu.upper, 3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"sigma.Rbrain2matblood"] <-
# signif(sigma.Rbrain2matblood, 3)
# PC.table.nofup [PC.table.nofup $CAS==this.cas,"R.brain.FtoM.upper"] <-
# pred.table5.nofup [pred.table5.nofup $CAS==this.cas,"Ratio.pred"]
# }
## ----compare_fup_correction_case_study,eval=FALSE-----------------------------
# case.study.fup.correct.table <- merge(PC.table,PC.table.nofup,by=c("Compound","CAS","DTXSID"))
#
# MFbrainfit.fup <- lm(R.brain.FtoM.upper.x~R.brain.FtoM.upper.y,
# data=case.study.fup.correct.table)
# summary(MFbrainfit.fup)
#
# cat(paste("The predictions for fetal brain to maternal blood ratio with or without the fetal plasma binding correction (Eqn. 1) were strongly correlated (R2 = ",
# signif(summary(MFbrainfit.fup)$adj.r.square,2),
# "). There were ",
# dim(subset(case.study.fup.correct.table,
# rank(R.brain.FtoM.upper.x) > rank(R.brain.FtoM.upper.y)))[1],
# " chemicals that were ranked higher with the correction than without, with an average increase of ",
# signif(mean(
# case.study.fup.correct.table$R.brain.FtoM.upper.y /
# case.study.fup.correct.table$R.brain.FtoM.upper.x),2),
# " times when the plasma binding change was included.\n",
# sep=""))
#
## ----fup_table, eval = execute.vignette---------------------------------------
# fup.table <- NULL
# all.chems <- get_cheminfo(
# model="fetal_pbtk",
# info="all",
# suppress.messages=TRUE)
# # Get rid of median fup 0:
# all.chems <- subset(all.chems,
# as.numeric(unlist(lapply(strsplit(
# all.chems$Human.Funbound.plasma,","),function(x) x[[1]])))!=0)
# for (this.chem in all.chems[,"CAS"])
# {
# temp <- parameterize_fetal_pbtk(chem.cas=this.chem,suppress.messages = TRUE)
# state <- calc_ionization(
# pH=7.26,
# pKa_Donor=temp$pKa_Donor,
# pKa_Accept=temp$pKa_Accept)
# if (state$fraction_positive > 0.5) this.charge <- "Positive"
# else if (state$fraction_negative > 0.5) this.charge <- "Negative"
# else this.charge <- "Neutral"
# this.row <- data.frame(DTXSID=all.chems[all.chems[,"CAS"]==this.chem,"DTXSID"],
# Compound=all.chems[all.chems[,"CAS"]==this.chem,"Compound"],
# CAS=this.chem,
# Fup.Mat.Pred = temp$Funbound.plasma,
# Fup.Neo.Pred = temp$Fraction_unbound_plasma_fetus,
# Charge = this.charge
# )
# fup.table <- rbind(fup.table,this.row)
# }
#
# fup.table[fup.table$Charge=="Positive","Charge"] <- paste("Positive (n=",
# sum(fup.table$Charge=="Positive"),
# ")",sep="")
# fup.table[fup.table$Charge=="Negative","Charge"] <- paste("Negative (n=",
# sum(fup.table$Charge=="Negative"),
# ")",sep="")
# fup.table[fup.table$Charge=="Neutral","Charge"] <- paste("Neutral (n=",
# sum(fup.table$Charge=="Neutral"),
# ")",sep="")
## ----fup_figure, eval = execute.vignette--------------------------------------
# FigA <- ggplot(data=fup.table) +
# geom_point(alpha=0.25, aes(
# x=Fup.Mat.Pred,
# y=Fup.Neo.Pred,
# shape=Charge,
# color=Charge),
# size=3) +
# geom_abline(slope=1, intercept=0) +
# ylab(expression(paste("Adjusted Neonate ",f[up]))) +
# xlab(expression(paste(italic("In vitro")," Measured Adult ",f[up]))) +
# scale_x_log10(label=scientific_10) +
# scale_y_log10(label=scientific_10) +
# theme_bw() +
# theme( text = element_text(size=14))
#
# print(FigA)
## ----MFratio_allhttk_preds, eval = execute.vignette---------------------------
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
#
# MFratio.pred <- NULL
# all.chems <- get_cheminfo(
# model="fetal_pbtk",
# info=c("Compound","DTXSID","CAS","Funbound.plasma"),
# suppress.messages=TRUE)
# for (this.cas in all.chems$CAS)
# if ((this.cas %in% nonvols) &
# !(this.cas %in% fluoros))
# {
# this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
# Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
# {
# skip <- TRUE
# } else skip <- FALSE
# } else if (Fup== 0)
# {
# skip <- TRUE
# } else skip <- FALSE
#
# if (!skip)
# {
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =TRUE,
# suppress.messages=TRUE)
# out <- solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# new.row <- data.frame(
# Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
# DTXSID = this.id,
# Mat.pred = mean(out[last.row,"Cplasma"]),
# Fet.pred = mean(out[last.row,"Cfplasma"]),
# MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# )
# MFratio.pred <- rbind(MFratio.pred,new.row)
# }
# }
## ----MFratio_allhttk_figure, eval = execute.vignette--------------------------
# FigD <- ggplot(data=MFratio.pred)+
# geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+
# xlab("Maternal:Fetal Plasma Concentration Ratio") +
# ylab("Number of chemicals")+
# theme_bw()+
# theme( text = element_text(size=14))
#
#
# print(FigD)
## ----MFratio_allhttk_stats, eval = execute.vignette---------------------------
# max.chem <- MFratio.pred[which(
# MFratio.pred$MFratio.pred==max(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
# min.chem <- MFratio.pred[which(
# MFratio.pred$MFratio.pred==min(MFratio.pred$MFratio.pred,na.rm=TRUE)),]
# cat(paste("In Figure X we examine the ratios predicted for the ",
# dim(MFratio.pred)[1],
# " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
# sep=""))
#
#
# cat(paste("We observe a median ratio of ",
# signif(median(MFratio.pred$MFratio.pred,na.rm=TRUE),3),
# ", with values ranging from ",
# signif(min.chem[,"MFratio.pred"],3),
# " for ",
# min.chem[,"DTXSID"],
# " to ",
# signif(max.chem[,"MFratio.pred"],3),
# " for ",
# max.chem[,"DTXSID"],
# ".\n",sep=""))
#
# # Check out phys-chem > 1.6, < 1:
# highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred>1.6)$DTXSID)
#
# cat(paste("There are",
# dim(highratio)[1],
# "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))
#
# # all highly bound
# highratio$Compound
# suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=TRUE)))
#
#
# lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred,MFratio.pred<0.9)$DTXSID)
# # No obvious pattern
## ----MFratio_allhttk_preds_nofup, eval = execute.vignette---------------------
# times <- sort(unique(c(seq(13 * 7, 40 * 7, 0.25),seq(278,280,.1))))
#
# MFratio.pred.nofup <- NULL
# all.chems <- get_cheminfo(
# model="fetal_pbtk",
# info=c("Compound","DTXSID","CAS","Funbound.plasma"),
# suppress.messages=TRUE)
# for (this.cas in all.chems$CAS)
# if ((this.cas %in% nonvols) &
# !(this.cas %in% fluoros))
# {
# this.id <- all.chems[all.chems$CAS==this.cas,"DTXSID"]
# Fup <- subset(all.chems,DTXSID==this.id)$Human.Funbound.plasma
# if (regexpr(",",Fup)!=-1)
# {
# if (as.numeric(strsplit(Fup,",")[[1]][1])==0 |
# (as.numeric(strsplit(Fup,",")[[1]][3])>0.9 &
# as.numeric(strsplit(Fup,",")[[1]][2])<0.1))
# {
# skip <- TRUE
# } else skip <- FALSE
# } else if (Fup== 0)
# {
# skip <- TRUE
# } else skip <- FALSE
#
# if (!skip)
# {
# p <- parameterize_fetal_pbtk(dtxsid=this.id,
# fetal_fup_adjustment =FALSE,
# suppress.messages=TRUE))
# out <- suppressWarnings(solve_fetal_pbtk(
# parameters=p,
# dose=0,
# times=times,
# daily.dose=1,
# doses.per.day=3,
# output.units = "uM",
# suppress.messages=TRUE)
# last.row <- which(out[,"time"]>279)
# last.row <- last.row[!duplicated(out[last.row,"time"])]
# new.row <- data.frame(
# Chemical = all.chems[all.chems$DTXSID==this.id,"Compound"],
# DTXSID = this.id,
# Mat.pred = mean(out[last.row,"Cplasma"]),
# Fet.pred = mean(out[last.row,"Cfplasma"]),
# MFratio.pred = mean(out[last.row,"Cplasma"])/mean(out[last.row,"Cfplasma"])
# )
# MFratio.pred.nofup <- rbind(MFratio.pred.nofup,new.row)
# }
# }
## ----MFratio_allhttk_figure_nofup, eval = execute.vignette--------------------
# FigSupD <- ggplot(data=MFratio.pred.nofup)+
# geom_histogram(binwidth = 0.05,fill="Red",aes(MFratio.pred))+
# xlab("Maternal:Fetal Plasma Concentration Ratio (No Fup Corr.)") +
# ylab("Number of chemicals")+
# theme_bw()+
# theme( text = element_text(size=14))
#
#
# print(FigSupD)
## ----MFratio_allhttk_stats_nofup, eval = execute.vignette---------------------
# max.chem <- MFratio.pred.nofup[which(
# MFratio.pred.nofup$MFratio.pred==max(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
# min.chem <- MFratio.pred.nofup[which(
# MFratio.pred.nofup$MFratio.pred==min(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE)),]
# cat(paste("In Figure X we examine the ratios predicted for the ",
# dim(MFratio.pred)[1],
# " appropriate (non-volatile or PFAS) chemicals with measured HTTK data.\n",
# sep=""))
#
#
# cat(paste("We observe a median ratio of ",
# signif(median(MFratio.pred.nofup$MFratio.pred,na.rm=TRUE),3),
# ", with values ranging from ",
# signif(min.chem[,"MFratio.pred"],3),
# " for ",
# min.chem[,"DTXSID"],
# " to ",
# signif(max.chem[,"MFratio.pred"],3),
# " for ",
# max.chem[,"DTXSID"],
# ".\n",sep=""))
#
# # Check out phys-chem > 1.6, < 1:
# highratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred>1.6)$DTXSID)
#
# cat(paste("There are",
# dim(highratio)[1],
# "chemicals with ratios greater than 1.6, indicating a tendency to have higher concentrations"))
#
# # all highly bound
# highratio$Compound
# suppressWarnings(apply(highratio,2,function(x) mean(as.numeric(x),na.rm=2)))
#
#
# lowratio <- subset(chem.physical_and_invitro.data,DTXSID%in%subset(MFratio.pred.nofup,MFratio.pred<0.9)$DTXSID)
# # No obvious pattern
## ----create_distributable_data, eval = execute.vignette-----------------------
# aylward2014 <-MFdata
# pregnonpregaucs <- TKstats
# fetalpcs <- Curley.pcs
# wang2018 <- Wangchems
#
# save(aylward2014,pregnonpregaucs,fetalpcs,wang2018,pksim.pcs,
# file="Kapraun2022Vignette.RData",version=2)
Try the httk package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.