## ----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, # MFratiothis.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)