Nothing
inst/doc/Vf_Wambaugh2019.R
In httk: High-Throughput Toxicokinetics
## ----configure_knitr, eval = TRUE---------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = '#>', fig.width=5, fig.height=3)
## ----clear_memory, eval = TRUE------------------------------------------------
rm(list=ls())
## ----runchunks, eval = TRUE---------------------------------------------------
# Set whether or not the following chunks will be executed (run):
execute.vignette <- FALSE
## ----setup, eval = execute.vignette-------------------------------------------
# library(data.table)
# library(httk)
# library(ggplot2)
# library(stringr)
# library(scales)
# library(grid)
# library(reshape)
## ----lmr2, eval = execute.vignette--------------------------------------------
# lm_R2 <- function(m,prefix=NULL){
# out <- substitute("RMSE = "~mse*","~~italic(R)^2~"="~r2,
# list(mse = signif(mean(m$residuals^2)^(1/2),3),
# r2 = format(summary(m)$adj.r.squared, digits = 2)))
# paste(prefix,as.character(as.expression(out)))
# }
## ----multi.plot, eval = execute.vignette--------------------------------------
# multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL, heights=NULL, widths=unit(rep_len(1, cols), "null")) {
#
# # Make a list from the ... arguments and plotlist
# plots <- c(list(...), plotlist)
#
# numPlots = length(plots)
#
# # If layout is NULL, then use 'cols' to determine layout
# if (is.null(layout)) {
# # Make the panel
# # ncol: Number of columns of plots
# # nrow: Number of rows needed, calculated from # of cols
# layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
# ncol = cols, nrow = ceiling(numPlots/cols))
# }
#
# if (numPlots==1) {
# print(plots[[1]])
#
# } else {
# # Set up the page
# grid.newpage()
# if (!is.null(heights))
# {
# pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout),heights=heights,widths=widths)))
# } else {
# pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout),widths=widths)))
# }
# # Make each plot, in the correct location
# for (i in 1:numPlots) {
# # Get the i,j matrix positions of the regions that contain this subplot
# matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
#
# print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
# layout.pos.col = matchidx$col))
# }
# }
# }
## ----scientific.notation, 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))
# }
## ----scientific.notation2, eval = execute.vignette----------------------------
# scientific_10_skip <- function(x) {
# out <- gsub("1e", "10^", scientific_format()(x))
# out <- gsub("\\+","",out)
# out <- gsub("10\\^01","10",out)
# out <- gsub("10\\^00","1",out)
# out <- gsub("10\\^-01",NA,out)
# return(parse(text=out))
# }
## ----cssfun.u, eval = execute.vignette----------------------------------------
# cssfun_u <- function(chemcas,
# MW,
# hepatic.bioavailability,
# Fgutabs=1,
# probs=0.95,
# minimum.Funbound.plasma=0.0001)
# {
# # Create a data table with uncertainty only:
# pop_u_httk <- as.data.table(parameterize_steadystate(chem.cas=chemcas,
# clint.pvalue.threshold=1,
# minimum.Funbound.plasma=0))
#
# # Get the physico-chemical properties:
# params.schmitt <-parameterize_schmitt(chem.cas=chemcas)
# psd <- params.schmitt$pKa_Donor
# psa <- params.schmitt$pKa_Accept
# params.schmitt <- params.schmitt[regexpr("pKa",names(params.schmitt))==-1]
# pop_u_httk[,names(params.schmitt):=params.schmitt]
# pop_u_httk[,pKa_Donor := paste(psd, collapse=",")]
# pop_u_httk[,pKa_Accept := paste(psa, collapse=",")]
# pop_u_httk <- pop_u_httk[rep(1,1000)]
#
# # Calculate the distribution coefficient:
# ions <- calc_ionization(pH=7.4,parameters=pop_u_httk)
# pop_u_httk[,Dow74 := Pow * (ions$fraction_neutral +
# 0.001 * ions$fraction_charged + ions$fraction_zwitter)]
#
# # Parse the Funbound.plasma.dist into quantiles:
# Funbound.plasma.dist <- pop_u_httk[1,Funbound.plasma.dist]
# if (nchar(Funbound.plasma.dist) -
# nchar(gsub(",","",Funbound.plasma.dist))!=2)
# {
# stop("Funbound.plasma distribution should be three values (median,low95th,high95th) separated by commas.")
# }
# temp <- strsplit(Funbound.plasma.dist,",")
# ppb.median <- as.numeric(temp[[1]][1])
# ppb.low95 <- as.numeric(temp[[1]][2])
# ppb.high95 <- as.numeric(temp[[1]][3])
#
# if (ppb.median>minimum.Funbound.plasma)
# {
# # Use optim to estimate alpha and beta for fup such that the median and 95%
# # credible interval approximate the estimate from MCMC:
# ppb.fit <- optim(c(2,(1-ppb.median)/ppb.median*2),
# function(x) (0.95-pbeta(ppb.high95,x[1],x[2])+
# pbeta(ppb.low95,x[1],x[2]))^2+
# (ppb.median-qbeta(0.5,x[1],x[2]))^2,
# method="BFGS")
# # We are drawing new values for the unadjusted Fup:
# pop_u_httk[, unadjusted.Funbound.plasma:=rbeta(1000,
# ppb.fit$par[1],
# ppb.fit$par[2])]
# } else if (ppb.high95>minimum.Funbound.plasma)
# {
# # Assume that since the median is zero but the u95 is not, that there is
# # an exponential distribution:
# # 97.5% of clearance values will be below Funbound.plasma.u95:
# pop_u_httk[,unadjusted.Funbound.plasma:=runif(1000,
# minimum.Funbound.plasma,
# min(1,minimum.Funbound.plasma+
# 2*(ppb.high95-minimum.Funbound.plasma)))]
# pop_u_httk[as.logical(rbinom(1000,1,.95)),
# unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# } else {
# pop_u_httk[, unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# }
#
# # then we need to adjust for in vitro binding:
# pop_u_httk[,Flipid:=subset(physiology.data,
# Parameter=='Plasma Effective Neutral Lipid Volume Fraction')[,
# which(colnames(physiology.data) == 'Human')]]
# pop_u_httk[,Funbound.plasma.adjustment:=1 / (Dow74 * Flipid +
# 1 / unadjusted.Funbound.plasma)/unadjusted.Funbound.plasma]
# pop_u_httk[,
# Funbound.plasma:=unadjusted.Funbound.plasma*Funbound.plasma.adjustment]
#
# # Enforce a minimum ppb:
# pop_u_httk[Funbound.plasma<minimum.Funbound.plasma,
# Funbound.plasma:=minimum.Funbound.plasma]
#
# # If Fup varies, then Rb2p and hepatic bioavailability also vary:
# Rblood2plasma <- get_rblood2plasma(chem.cas=chemcas)
# if (is.na(Rblood2plasma))
# {
# pop_u_httk[, Rblood2plasma:=calc_rblood2plasma(params=pop_u_httk)]
# # Right now we don't similate variability on a known blood:plasma ratio
# } else pop_u_httk[, Rblood2plasma:=Rblood2plasma]
#
# # Parse the clint.dist to retrieve the quantiles:
# Clint.dist <- pop_u_httk[1,Clint.dist]
# if (nchar(Clint.dist) -
# nchar(gsub(",","",Clint.dist))!=3)
# {
# stop("Clint distribution should be four values (median,low95th,high95th,pValue) separated by commas.")
# }
# temp <- strsplit(Clint.dist,",")
# Clint.median <- as.numeric(temp[[1]][1])
# Clint.low95 <- as.numeric(temp[[1]][2])
# Clint.high95 <- as.numeric(temp[[1]][3])
# Clint.pvalue <- as.numeric(temp[[1]][4])
#
# # Use optim to estimate mean and standard deviation of Clint such that the
# # median and 95% credible interval approximate the estimate from MCMC:
# if (Clint.high95>0)
# {
# if (Clint.median>0)
# {
# clint.fit <- optim(c(log(Clint.median),0.3), function(x)
# (0.95-plnorm(Clint.high95,x[1],x[2])+plnorm(Clint.low95,x[1],x[2]))^2+
# (Clint.median-qlnorm(0.5,x[1],x[2]))^2)
# pop_u_httk[,Clint:=rlnorm(1000,clint.fit$par[1],clint.fit$par[2])]
# } else {
# pop_u_httk[,Clint:=exp(runif(1000,log(1),
# (log(Clint.high95)-log(1))/0.975))]
# }
# } else pop_u_httk[,Clint:=0]
# pop_u_httk[as.logical(rbinom(1000,1,Clint.pvalue)),Clint:=0] # the Bayesian "p-value" here reflects how often there is no clearance
#
# pop_u_httk[, CLh:=calc_hepatic_clearance(chem.cas=chemcas,
# parameters=pop_u_httk,hepatic.model='unscaled',suppress.messages=TRUE,clint.pvalue.threshold=1)]
# pop_u_httk[, Qliver:=Qtotal.liverc*BW^0.75]
# pop_u_httk[, hepatic.bioavailability:= Qliver / (Qliver +
# Funbound.plasma * Clint / Rblood2plasma)]
# pop_u_httk[, (setdiff(names(pop_u_httk),
# names(httk::parameterize_steadystate(chem.cas=chemcas)))):=NULL]
#
#
# #compute Css for each individual
# # calc_analytic_css is vectorized so that this will work.
# css <- httk::calc_analytic_css(parameters=pop_u_httk,
# clint.pvalue.threshold=1,
# daily.dose=1,
# clint.pop.cv=NULL,
# fup.pop.cv=NULL,
# output.units="uM",
# model="3compartmentss",
# species="Human",
# suppress.messages=TRUE)#,
# # recalc.blood2plasma=TRUE)
# #get Css95
# cssquants <- quantile(css, probs=probs)
# return(cssquants)
# }
## ----hlfun.u, eval = execute.vignette-----------------------------------------
# hlfun_u <- function(chemcas,
# MW,
# probs=0.95,
# minimum.Funbound.plasma=0.0001)
# {
# # Create a data table with uncertainty only:
# print(chemcas)
# pop_u_httk <- as.data.table(parameterize_steadystate(chem.cas=chemcas,
# clint.pvalue.threshold=1,
# minimum.Funbound.plasma=0))
# params.schmitt <-parameterize_schmitt(chem.cas=chemcas)
# psd <- params.schmitt$pKa_Donor
# psa <- params.schmitt$pKa_Accept
# params.schmitt <- params.schmitt[regexpr("pKa",names(params.schmitt))==-1]
# pop_u_httk[,names(params.schmitt):=params.schmitt]
# pop_u_httk[,pKa_Donor := paste(psd, collapse=",")]
# pop_u_httk[,pKa_Accept := paste(psa, collapse=",")]
# pop_u_httk <- pop_u_httk[rep(1,1000)]
#
# ions <- calc_ionization(pH=7.4,parameters=pop_u_httk)
# pop_u_httk[,Dow74 := Pow * (ions$fraction_neutral + 0.001 * ions$fraction_charged + ions$fraction_zwitter)]
#
# # Parse the Funbound.plasma.dist into quantiles:
# Funbound.plasma.dist <- pop_u_httk[1,Funbound.plasma.dist]
# print(paste(chemcas,Funbound.plasma.dist))
# if (nchar(Funbound.plasma.dist) -
# nchar(gsub(",","",Funbound.plasma.dist))!=2)
# {
# stop("Funbound.plasma distribution should be three values (median,low95th,high95th) separated by commas.")
# }
# temp <- strsplit(Funbound.plasma.dist,",")
# ppb.median <- as.numeric(temp[[1]][1])
# ppb.low95 <- as.numeric(temp[[1]][2])
# ppb.high95 <- as.numeric(temp[[1]][3])
#
# if (ppb.median>minimum.Funbound.plasma)
# {
# # Use optim to estimate alpha and beta for fup such that the median and 95%
# # credible interval approximate the estimate from MCMC:
# ppb.fit <- optim(c(2,(1-ppb.median)/ppb.median*2),
# function(x) (0.95-pbeta(ppb.high95,x[1],x[2])+
# pbeta(ppb.low95,x[1],x[2]))^2+
# (ppb.median-qbeta(0.5,x[1],x[2]))^2,
# method="BFGS")
# # We are drawing new values for the unadjusted Fup:
# pop_u_httk[, unadjusted.Funbound.plasma:=rbeta(1000,
# ppb.fit$par[1],
# ppb.fit$par[2])]
# } else if (ppb.high95>minimum.Funbound.plasma)
# {
# # Assume that since the median is zero but the u95 is not, that there is
# # an exponential distribution:
# # 97.5% of clearance values will be below Funbound.plasma.u95:
# pop_u_httk[, unadjusted.Funbound.plasma:=runif(1000,
# minimum.Funbound.plasma,
# min(1,(ppb.high95-minimum.Funbound.plasma)/0.975))]
# pop_u_httk[as.logical(rbinom(1000,1,.975)),
# unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# } else {
# pop_u_httk[, unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# }
#
# # then we need to adjust for in vitro binding:
# pop_u_httk[,Flipid:=subset(physiology.data,
# Parameter=='Plasma Effective Neutral Lipid Volume Fraction')[,
# which(colnames(physiology.data) == 'Human')]]
# pop_u_httk[,Funbound.plasma.adjustment:=1 / (Dow74 * Flipid +
# 1 / unadjusted.Funbound.plasma)/unadjusted.Funbound.plasma]
# pop_u_httk[,
# Funbound.plasma:=unadjusted.Funbound.plasma*Funbound.plasma.adjustment]
#
# # Enforce a minimum ppb:
# pop_u_httk[Funbound.plasma<minimum.Funbound.plasma,
# Funbound.plasma:=minimum.Funbound.plasma]
#
# Rblood2plasma <- get_rblood2plasma(chem.cas=chemcas)
# if (is.na(Rblood2plasma))
# {
# pop_u_httk[, Rblood2plasma:=calc_rblood2plasma(params=pop_u_httk)]
# } else pop_u_httk[, Rblood2plasma:=Rblood2plasma] # Right now we don't similate variability on a known blood:plasma ratio
# # Parse the clint.dist to retrieve the quantiles:
# Clint.dist <- pop_u_httk[1,Clint.dist]
# if (nchar(Clint.dist) -
# nchar(gsub(",","",Clint.dist))!=3)
# {
# stop("Clint distribution should be four values (median,low95th,high95th,pValue) separated by commas.")
# }
# temp <- strsplit(Clint.dist,",")
# Clint.median <- as.numeric(temp[[1]][1])
# Clint.low95 <- as.numeric(temp[[1]][2])
# Clint.high95 <- as.numeric(temp[[1]][3])
# Clint.pvalue <- as.numeric(temp[[1]][4])
#
# # Use optim to estimate mean and standard deviation of Clint such that the
# # median and 95% credible interval approximate the estimate from MCMC:
# if (Clint.high95>0)
# {
# if (Clint.median>0)
# {
# clint.fit <- optim(c(log(Clint.median),0.3), function(x)
# (0.95-plnorm(Clint.high95,x[1],x[2])+plnorm(Clint.low95,x[1],x[2]))^2+
# (Clint.median-qlnorm(0.5,x[1],x[2]))^2)
# pop_u_httk[,Clint:=rlnorm(1000,clint.fit$par[1],clint.fit$par[2])]
# } else {
# pop_u_httk[,Clint:=exp(runif(1000,log(1),
# (log(Clint.high95)-log(1))/0.975))]
# }
# } else pop_u_httk[,Clint:=0]
# pop_u_httk[as.logical(rbinom(1000,1,Clint.pvalue)),Clint:=0] # the Bayesian "p-value" here reflects how often there is no clearance
#
# pop_u_httk[, CLh:=calc_hepatic_clearance(chem.cas=chemcas,parameters=pop_u_httk,hepatic.model='unscaled',suppress.messages=TRUE,clint.pvalue.threshold=1)]
# pop_u_httk[, Qliver:=Qtotal.liverc*BW^0.75]
# pop_u_httk[, hepatic.bioavailability:= Qliver / (Qliver + Funbound.plasma * Clint / Rblood2plasma)]
#
# PC.names <- names(predict_partitioning_schmitt(chem.cas="80-05-7"))
# pop_u_httk<-cbind(pop_u_httk,as.data.table(predict_partitioning_schmitt(parameters=pop_u_httk)))
# pop_u_httk[,hl:=calc_elimination_rate(parameters=pop_u_httk)]
#
# return(log(2)/quantile(unlist(pop_u_httk[,"hl",with=FALSE]),probs=probs))
# }
## ----Wambaugh2019.Fig1, eval = execute.vignette-------------------------------
# Fig1.table <- httk::wambaugh2019.raw
# Fig1.table$Error1 <- Fig1.table$Base.Fup.High - Fig1.table$Base.Fup.Low
# Fig1.table$Error2 <- Fig1.table$Affinity.Fup.High - Fig1.table$Affinity.Fup.Low
# Fig1.table$Sigma1 <- Fig1.table$Error1/(2*qnorm(0.975))
# Fig1.table$Mean1 <- (Fig1.table$Base.Fup.High + Fig1.table$Base.Fup.Low)/(qnorm(0.975))
# Fig1.table$CV1 <- Fig1.table$Sigma1/Fig1.table$Mean1
# Fig1.table$Sigma2 <- Fig1.table$Error2/(2*qnorm(0.975))
# Fig1.table$Mean2 <- (Fig1.table$Base.Fup.High + Fig1.table$Base.Fup.Low)/(qnorm(0.975))
# Fig1.table$CV2 <- Fig1.table$Sigma2/Fig1.table$Mean2
# Fig1.table$ErrorClint <- Fig1.table$CLint.1uM.High95th - Fig1.table$CLint.1uM.Low95th
# Fig1.table$SigmaClint <- Fig1.table$ErrorClint/(2*qnorm(0.975))
# Fig1.table$MeanClint <- (Fig1.table$CLint.1uM.High95th + Fig1.table$CLint.1uM.Low95th)/(qnorm(0.975))
# Fig1.table$CVClint <- Fig1.table$SigmaClint/Fig1.table$MeanClint
#
# Fupmeasured <- subset(Fig1.table,!is.na(Affinity.Fup.Med))
# Fupmeasured <- subset(Fupmeasured,!duplicated(CAS))
# CLintmeasured <- subset(Fig1.table,!is.na(CLint.1uM.Median))
# CLintmeasured <- subset(CLintmeasured,!duplicated(CAS))
#
#
#
# print(paste("New HTTK experimental measurements were made for",length(unique(Fig1.table$CAS)),"chemicals from the ToxCast HTS library."))
#
# print(paste("Fup was successfully measured for",length(unique(subset(Fig1.table,!is.na(Affinity.Fup.Med))$CAS)),"chemicals from the ToxCast HTS library."))
#
# print(paste("CLint was successfully measured for",length(unique(subset(Fig1.table,!is.na(CLint.1uM.Median))$CAS)),"chemicals from the ToxCast HTS library."))
#
#
#
#
# Fig1a <- ggplot(Fig1.table)+
# geom_freqpoly(binwidth = 0.5,lwd=1.5,color="Red",aes(Base.Fup.Med))+
# geom_freqpoly(binwidth = 0.5,lwd=1.5,lty=3,color="Blue",aes(Affinity.Fup.Med))+
# xlab(expression(paste("Measured ",F[up]))) +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(10^-6,1.05))+
# labs(title=paste(length(unique(subset(Fig1.table,!is.na(Affinity.Fup.Med))$CAS)),"Chemicals Measured"))+
# theme_bw()+
# theme( text = element_text(size=10)) +
# annotate("text", x=10^-6,y=80,size=8,label="a")
#
# print(paste("Median 100%:",signif(median(Fig1.table$Base.Fup.Med,na.rm=TRUE),2)))
# print(paste("Median Titration:",signif(median(Fig1.table$Affinity.Fup.Med,na.rm=TRUE),2)))
#
#
#
# Fig1b <- ggplot(Fig1.table)+
# geom_histogram(binwidth = 0.5,alpha=0.6,fill="green",aes(CLint.1uM.Median))+
# xlab(expression(paste("Measured ",Cl[int]," (uL/min/",10^6," hepatocytes)"))) +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(10^-1,5*10^3))+
# labs(title=paste(length(unique(subset(Fig1.table,!is.na(CLint.1uM.Median))$CAS)),"Chemicals Measured"))+
# theme_bw()+
# theme( text = element_text(size=10)) +
# annotate("text", x=2*10^-1,y=80,size=8,label="b")
#
# print(paste(sum(Fig1.table$CLint.1uM.Median==0,na.rm=TRUE),"Zero Measurments"))
# print(paste("Median Non-Zero:",signif(median(Fig1.table$CLint.1uM.Median,na.rm=TRUE),3)))
#
#
#
# Fig1c <- ggplot(Fig1.table)+
# geom_histogram(binwidth = 0.1,alpha=0.2,fill="Red",aes(Error1))+
# geom_histogram(binwidth = 0.1,alpha=0.2,fill="Blue",aes(Error2))+
# xlab("Width of Credible Interval") +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(5*10^-4,1.05))+
# labs(title=paste("CI:",signif(median(Fig1.table$Error1,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CV1,na.rm=TRUE),1),"/ CI:",signif(median(Fig1.table$Error2,na.rm=TRUE),1),"CV:",signif(median(Fig1.table$CV2,na.rm=TRUE),1)))+
# theme_bw()+
# theme( text = element_text(size=10))+
# annotate("text", x=5*10^-4,y=42,size=8,label="c")
#
# print(paste("Median 100% CI:",signif(median(Fig1.table$Error1,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CV1,na.rm=TRUE),2)))
# print(paste("Median Titr. CI:",signif(median(Fig1.table$Error2,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CV2,na.rm=TRUE),2)))
#
#
# Fig1d <- ggplot(Fig1.table)+
# geom_histogram(binwidth = 0.1,alpha=0.6,fill="green",aes(ErrorClint))+
# xlab("Width of Credible Interval") +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(10^0,10^3))+
# labs(title=paste("Median CI:",signif(median(Fig1.table$ErrorClint,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CVClint,na.rm=TRUE),2)))+
# theme_bw()+
# theme( text = element_text(size=10))+
# annotate("text", x=1.5,y=25,size=8,label="d")
#
# print(paste("Median CI:",signif(median(Fig1.table$ErrorClint,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CVClint,na.rm=TRUE),2)))
#
# multiplot(Fig1a,Fig1c,Fig1b,Fig1d,cols=2,widths=c(1.75,1.75))
## ----Wambaugh2019.Fig2, eval = execute.vignette-------------------------------
# Fig2.table <- httk::wambaugh2019.raw
# print(paste(sum(!is.na(Fig2.table$Fup.point)),"total successful Fup estimates using point estimation."))
# print(paste(sum(!is.na(Fig2.table$Base.Fup.Med)),"total successful Fup estimates using BAyesian analysis of top protein concentration."))
# Fig2.table$Error1 <- Fig2.table$Base.Fup.High - Fig2.table$Base.Fup.Low
# Fig2.table$Error2 <- Fig2.table$Affinity.Fup.High - Fig2.table$Affinity.Fup.Low
# Fig2.table$Sigma1 <- Fig2.table$Error1/(2*qnorm(0.975))
# Fig2.table$Mean1 <- (Fig2.table$Base.Fup.High + Fig2.table$Base.Fup.Low)/(qnorm(0.975))
# Fig2.table$CV1 <- Fig2.table$Sigma1/Fig2.table$Mean1
# Fig2.table$Sigma2 <- Fig2.table$Error2/(2*qnorm(0.975))
# Fig2.table$Mean2 <- (Fig2.table$Base.Fup.High + Fig2.table$Base.Fup.Low)/(qnorm(0.975))
# Fig2.table$CV2 <- Fig2.table$Sigma2/Fig2.table$Mean2
# Fig2.table$Point.Estimate <- "Good"
# Fig2.table[Fig2.table$Fup.point<0&!is.na(Fig2.table$Fup.point),"Point.Estimate"]<-"<0"
# Fig2.table[Fig2.table$Fup.point>1&!is.na(Fig2.table$Fup.point),"Point.Estimate"]<-">1"
# Fig2.table[is.na(Fig2.table$Fup.point),"Point.Estimate"]<-"No Value"
# Fig2.table[Fig2.table$Point.Estimate%in%c("<0","No Value"),"Fup.point"] <- Fig2.table[Fig2.table$Point.Estimate%in%c("<0","No Value"),"Base.Fup.Med"]
# Fig2.table$Uncertain<-"Yes"
# Fig2.table[is.na(Fig2.table$CV1),"CV1"] <- -999
# Fig2.table[Fig2.table$CV1<0.49,"Uncertain"] <- "No"
# Fig2.table[Fig2.table$CV1==-999,"CV1"] <- NA
#
# Fig2a <- ggplot(subset(Fig2.table,!is.na(Base.Fup.Med)),aes(x=Fup.point,y=Base.Fup.Med,shape=Point.Estimate,alpha=Uncertain,colour=Point.Estimate)) +
# geom_point(size=3)+
# # geom_errorbar(aes(ymax = Base.Fup.High, ymin=Base.Fup.Low))+
# scale_y_log10(label=scientific_10,limits=c(10^-6,2)) +
# scale_x_log10(label=scientific_10,limits=c(10^-6,2)) +
# ylab(expression(paste(F[up]," Bayesian Analysis (Single Conc.)"))) +
# xlab(expression(paste(F[up]," Point Estimate (Traditional Analysis)"))) +
# geom_abline(intercept = 0, slope = 1,linetype="dashed", colour="Blue") +
# # geom_vline(xintercept = 1, size=2,linetype="dashed", colour="red")+
# geom_vline(xintercept = 0.01, size=2,linetype="dotted", colour="red")+
# geom_hline(yintercept = 0.01, size=2,linetype="dotted", colour="red")+
# scale_colour_discrete(name="Point Estimate")+
# scale_shape_discrete(name="Point Estimate") +
# scale_alpha_manual(values=c(1,0.3))+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=12))+
# annotate("text", size=8,x = 10^-6, y = 2, label = "a")+
# guides(alpha=guide_legend(nrow=2,byrow=TRUE),colour=guide_legend(nrow=2,byrow=TRUE),shape=guide_legend(nrow=2,byrow=TRUE))
#
# #print(Fig2a)
#
# print(paste(sum(Fig2.table$Uncertain=="Yes"),"uncertain estimates using single conc."))
# print(paste(sum(Fig2.table$Point.Estimate=="<0"),"chemicals with negative point estimates."))
# print(paste(sum(Fig2.table$Point.Estimate==">1"),"chemicals with Fup > 1."))
#
#
#
# summary(lm(data=subset(Fig2.table,!is.na(Fup.point)&Fup.point>0&Base.Fup.Med>0),log(Fup.point)~log(Base.Fup.Med)))
#
#
#
#
# Fig2.table$Uncertain2<-"Yes"
# Fig2.table[is.na(Fig2.table$CV2),"CV2"] <- -999
# Fig2.table[Fig2.table$CV2<0.49,"Uncertain2"] <- "No"
# Fig2.table[Fig2.table$CV2==-999,"CV2"] <- NA
# Fig2.table$TopFail <- ""
# Fig2.table[!is.na(Fig2.table$Affinity.Fup.Med)&Fig2.table$Uncertain2=="No","TopFail"]<-"No"
# Fig2.table[Fig2.table$TopFail=="No"&!is.na(Fig2.table$Base.Fup.Med)&Fig2.table$Uncertain=="Yes","TopFail"]<-"Single Conc. Only"
# Fig2.table[!is.na(Fig2.table$Affinity.Fup.Med)&Fig2.table$Uncertain2=="Yes","TopFail"]<-"Yes"
# Fig2.table[!is.na(Fig2.table$Affinity.Fup.Med)&is.na(Fig2.table$Base.Fup.Med),"TopFail"] <- "Single Conc. Only"
# Fig2.table[is.na(Fig2.table$Base.Fup.Med),"Base.Fup.Med"] <- Fig2.table[is.na(Fig2.table$Base.Fup.Med),"Affinity.Fup.Med"]
#
#
# Fig2b <- ggplot(subset(Fig2.table,!is.na(Affinity.Fup.Med)),aes(x=Base.Fup.Med,y=Affinity.Fup.Med,shape=TopFail,alpha=TopFail,colour=TopFail)) +
# geom_point(size=3)+
# geom_point(data=subset(Fig2.table,TopFail=="No"),size=3)+
# # geom_errorbar(aes(ymax = Fup.High.x, ymin=Fup.Low.x))+
# # geom_errorbarh(aes(xmin=Fup.Low.y,xmax=Fup.High.y))+
# scale_y_log10(label=scientific_10,limits=c(10^-8,2)) +
# scale_x_log10(label=scientific_10,limits=c(10^-8,2)) +
# xlab(expression(paste(F[up]," Bayesian Analysis (Single Conc.)"))) +
# ylab(expression(paste(F[up]," Bayesian Analysis (Three Conc.)"))) +
# geom_abline(intercept = 0, slope = 1,linetype="dashed", colour="Blue")+
# scale_colour_discrete(name="Uncertain")+
# scale_shape_discrete(name="Uncertain") +
# # annotate("text",size=8, x = 10^-4, y = 1, label = "B")+
# scale_alpha_manual(name="Uncertain",values=c(1,0.6,0.3))+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=12))+
# annotate("text", size=8,x = 10^-8, y = 2, label = "b")+
# guides(alpha=guide_legend(nrow=2,byrow=TRUE),colour=guide_legend(nrow=2,byrow=TRUE),shape=guide_legend(nrow=2,byrow=TRUE))
#
# #print(Fig2b)
#
# # Need to reduce to unique chemicals.
# Fupbaseworked <- subset(Fig2.table,!is.na(Base.Fup.Med)&Uncertain=="No")
# Fupbasenot <- subset(Fig2.table,!is.na(Base.Fup.Med)&Uncertain=="Yes")
# Fupbaseworked <- subset(Fupbaseworked,!duplicated(CAS))
# Fupbasenot <- subset(Fupbasenot,!duplicated(CAS))
# # Also need to count a chemical as a success if it worked for either duplicate:
# Fupbasenot <- subset(Fupbasenot,!(CAS%in%Fupbaseworked$CAS))
#
# Fupaffinityworked <- subset(Fig2.table,!is.na(Affinity.Fup.Med)&Uncertain2=="No")
# Fupaffinitynot <- subset(Fig2.table,!is.na(Affinity.Fup.Med)&Uncertain2=="Yes")
# Fupaffinityworked <- subset(Fupaffinityworked,!duplicated(CAS))
# Fupaffinitynot <- subset(Fupaffinitynot,!duplicated(CAS))
# Fupaffinitynot <- subset(Fupaffinitynot,!(CAS%in%Fupaffinityworked$CAS))
#
#
# print(paste(dim(Fupbaseworked)[1],"total successful Fup estimates using single conc."))
# print(paste(signif(dim(Fupbasenot)[1]/(dim(Fupbaseworked)[1]+dim(Fupbasenot)[1])*100,3),"percent failure rate using single conc."))
# print(paste(dim(Fupaffinityworked)[1],"total successful Fup estimates using protein titration."))
# print(paste(signif(dim(Fupaffinitynot)[1]/(dim(Fupaffinityworked)[1]+dim(Fupaffinitynot)[1])*100,3),"percent failure rate using protein titration."))
# print(paste(dim(Fupbasenot)[1],"uncertain estimates using original protocol."))
# print(paste(dim(Fupaffinitynot)[1],"uncertain estimates using protein titration."))
# print(paste(dim(Fupaffinityworked)[1]-dim(Fupbaseworked)[1],"chemicals that could not be measured at 100% plasma protein."))
#
# Wetmore.data <- httk:::Wetmore.data
# W2015 <- subset(Wetmore.data,Reference=="Wetmore 2015")
#
# W2012 <- subset(Wetmore.data,Reference=="Wetmore 2012")
# print(paste("The Fup assay failed for ",signif(sum(W2012$Fub==0.005)/length(W2012$Fub)*100,3),"% of chemicals in {Wetmore, 2012} and ",signif(sum(W2015$Fub==0)/length(W2015$Fub)*100,3),"% of chemicals in {Wetmore, 2015}."))
#
#
# Fuplod <- subset(Fig2.table,TopFail=="Single Conc. Only")
# print(paste("Among the chemicals that were uncertain using the original protocol but could be measured using the new three concentration protocol, the median estimated Fup was ",signif(median(Fuplod$Affinity.Fup.Med),2),"with values as low as ",signif(min(Fuplod$Affinity.Fup.Med),2),"and as high as",signif(max(Fuplod$Affinity.Fup.Med),2),". Previously, a default of 0.005 has been assumed when Fup could not be measured {Rotroff, 2010}.}"))
#
#
# summary(lm(data=subset(Fig2.table,!is.na(Base.Fup.Med)&Affinity.Fup.Med>0&Base.Fup.Med>0),log(Base.Fup.Med)~log(Affinity.Fup.Med)))
#
# #dev.new(width=10,height=6)
# multiplot(Fig2a,Fig2b,cols=2,widths=c(1.75,1.75))
## ----Wambaugh2019.Fig3, eval = execute.vignette-------------------------------
# Fig1.table$Class <- "Other"
# Fig1.table[Fig1.table$DTXSID %in% pharma$DTXSID,"Class"]<-"Pharmaceutical"
# Fig1.table$Class <- as.factor(Fig1.table$Class)
#
# Fig3 <- ggplot(Fig1.table,aes(Affinity.Kd.Med,fill=Class))+
# geom_histogram(binwidth = 0.5)+
# xlab(expression(paste("Binding Affinity (", mu,"M)",sep=""))) +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10)+
# scale_fill_discrete(name="Chemical Class")+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=18))
#
# # dev.new()
# print(Fig3)
## ----Wambaugh2019.Fig4, eval = execute.vignette-------------------------------
# clearance.table <- subset(httk::wambaugh2019.raw,!is.na(CLint.1uM.Median))
#
# print(paste(sum(clearance.table$CLint.1uM.Median==0),"zero valued median Bayesian posteriors"))
# print(paste(sum(clearance.table$CLint.1uM.Point==0&clearance.table$CLint.1uM.Median>0,na.rm=TRUE),"zero valued point estimates where median posterior>0"))
# print(paste(sum(clearance.table$CLint.1uM.Point>0&clearance.table$CLint.1uM.Median==0,na.rm=TRUE),"zero valued median Bayesian posteriors where point estimate was non-zero"))
# print(paste(sum(is.na(clearance.table$CLint.1uM.Point)&clearance.table$CLint.1uM.Median==0,na.rm=TRUE),"zero valued median Bayesian posteriors where point estimate was NA"))
# print(paste(sum(is.na(clearance.table$CLint.1uM.Point)&clearance.table$CLint.1uM.Median>0,na.rm=TRUE),"non-zero valued median Bayesian posteriors where point estimate was NA"))
#
#
#
#
# clearance.table$Fit <- "Decreasing"
# clearance.table[clearance.table$CLint.1uM.Median == 0,"Fit"] <- "Median Zero"
# #clearance.table[is.na(clearance.table$CLint.1uM.Point == 0),"Fit"] <- "Point Est. Missing"
# for (i in 1:dim(clearance.table)[1])
# if (!is.na(clearance.table[i,"CLint.1uM.Point"]))
# {
# if (clearance.table[i,"CLint.1uM.Point"] == 0) clearance.table[i,"Fit"] <- "Point Est. Zero"
# } else clearance.table[i,"Fit"] <- "Point Est. Failed"
#
#
# set.seed(123456)
# clearance.table[is.na(clearance.table$CLint.1uM.Point),"CLint.1uM.Point"]<-runif(sum(is.na(clearance.table$CLint.1uM.Point)),0.1,0.2)
# clearance.table[clearance.table$CLint.1uM.Point==0,"CLint.1uM.Point"]<-runif(sum(clearance.table$CLint.1uM.Point==0),0.3,0.8)
# clearance.table[clearance.table$CLint.1uM.Low95th==0,"CLint.1uM.Low95th"]<-0.1
# clearance.table[clearance.table$CLint.1uM.Median==0,"CLint.1uM.Median"]<-0.1
# clearance.table[clearance.table$CLint.1uM.High95th==0,"CLint.1uM.High95th"]<-0.1
# clearance.table[clearance.table$CLint.1uM.High95th>1000,"CLint.1uM.High95th"]<-1000
# clearance.table[clearance.table$CLint.1uM.Low95th<0.1,"CLint.1uM.Low95th"]<-0.1
#
# clearance.fit <- lm(log10(CLint.1uM.Median)~log10(CLint.1uM.Point),data=subset(clearance.table,!is.na(CLint.1uM.Point)&CLint.1uM.Point>1))
# Fig4 <- ggplot(clearance.table, aes(x=CLint.1uM.Point,y=CLint.1uM.Median,colour=Fit)) +
# geom_segment(aes(x=CLint.1uM.Point,xend=CLint.1uM.Point,y=CLint.1uM.Low95th,yend=CLint.1uM.High95th),size=1,alpha=0.5)+
# geom_point(size=3) +
# scale_y_log10(label=scientific_10,limits=c(10^-1,1000)) +
# scale_x_log10(label=scientific_10_skip,limits=c(10^-1,1000)) +
# xlab(expression(paste(CL[int]," (",mu,"L/min/",10^{6}," Hep.) Point Estimate",sep="")))+
# ylab(expression(paste(CL[int]," (",mu,"L/min/",10^{6}," Hep.) Bayesian",sep="")))+
# geom_vline(xintercept = 0.25, size=2,colour="black")+
# geom_vline(xintercept = 0.9, size=2,colour="black")+
# geom_abline(intercept = 0, slope = 1,linetype="dashed")+
# scale_colour_discrete(name="Assay Result")+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=14))+
# guides(colour=guide_legend(nrow=2,byrow=TRUE))+
# annotate("text",x = 10^2, y = 10^0,size=6, label = lm_R2(clearance.fit,prefix=""),parse=TRUE)
#
# #dev.new()
# print(Fig4)
#
# summary(clearance.fit)
## ----Wambaugh2019.Fig5, eval = execute.vignette-------------------------------
# # This section takes a long time because calc_vdist is not yet data.table compatible:
# ceetox <- as.data.table(subset(httk::wambaugh2019,!is.na(Human.Funbound.plasma)&!is.na(Human.Clint)&!is.na(logP)))
# httk_cas <- ceetox$CAS[ceetox$CAS %in% get_cheminfo(model="3compartmentss")]
#
# # back up chem.physical_and_invtro.data:
# HTTK.data <- chem.physical_and_invitro.data
#
# # Use the point estimates:
# chem.physical_and_invitro.data <- add_chemtable(httk::wambaugh2019.raw,current.table=chem.physical_and_invitro.data,data.list=list(Compound="Name",CAS="CAS",Clint="CLint.1uM.Point",Funbound.plasma="Fup.point"),reference="Wambaugh 2019",species="Human",overwrite=TRUE)
#
#
#
# schmitt_cas <- httk_cas[httk_cas %in% get_cheminfo(model="schmitt")]
# ceetox[CAS%in%schmitt_cas,
# hlpoint:=log(2)/httk::calc_elimination_rate(chem.cas=CAS),
# by=.(CAS)]
#
# chem.physical_and_invitro.data <- HTTK.data
#
# ceetox[CAS%in%schmitt_cas,
# hl_975:=hlfun_u(chemcas=CAS,
# MW=MW,
# probs=0.975),
# by=.(CAS)]
#
# ceetox[CAS%in%schmitt_cas,
# hl_med:=hlfun_u(chemcas=CAS,
# MW=MW,
# probs=0.5),
# by=.(CAS)]
#
# ceetox[CAS%in%schmitt_cas,
# hl_25:=hlfun_u(chemcas=CAS,
# MW=MW,
# probs=0.025),
# by=.(CAS)]
#
# ceetox[,hlpointEstimated:=TRUE]
# ceetox[is.na(hlpoint),hlpointEstimated:=FALSE]
# ceetox[is.na(hlpoint),hlpoint:=hl_med]
#
# hlfit <- lm(data=subset(ceetox,hlpointEstimated),log10(hlpoint)~log10(hl_med))
# Fig5a <- ggplot(data=subset(ceetox,hlpointEstimated),aes(x=hlpoint,y=hl_med)) +
# geom_point(size=3,alpha=0.5) +
# geom_errorbar(aes(ymax = hl_975, ymin=hl_25),alpha=0.3)+
# scale_y_log10(label=scientific_10,limits=c(1,10^6)) +
# scale_x_log10(label=scientific_10,limits=c(1,10^6)) +
# xlab("half-life (h) Point Estimate")+
# ylab("Bayesian half-life (h)")+
# geom_abline(intercept = 0, slope = 1,linetype="dashed",color="blue")+
# annotate("text", size=8,x = 1, y = 10^6, label = "A")+
# theme_bw() +
# theme(text = element_text(size=18))+
# annotate("text",x = 3*10^3, y = 10^0,size=5, label = lm_R2(hlfit,prefix=""),parse=TRUE)
#
#
# cssfit <- lm(data=ceetox,log10(cssmed)~log10(cssmed_u))
# Fig5b <- ggplot(data=ceetox,aes(x=cssmed,y=cssmed_u)) +
# geom_point(size=3,alpha=0.5) +
# geom_errorbar(aes(ymin=css25_u,ymax=css975_u),alpha=0.3)+
# scale_y_log10(label=scientific_10,limits=c(10^-3,10^5)) +
# scale_x_log10(label=scientific_10,limits=c(10^-3,10^5)) +
# xlab(expression(paste(C[ss]," Point Estimate (uM)")))+
# ylab(expression(paste("Bayesian ",C[ss]," (uM)")))+
# annotate("text", size=8,x = 10^-3, y = 10^5, label = "B")+
# theme_bw() +
# theme(text = element_text(size=18)) +
# geom_abline(intercept = 0, slope = 1,linetype="dashed",color="blue")+
# annotate("text",x = 10^2, y = 10^-3,size=5, label = lm_R2(cssfit,prefix=""),parse=TRUE)
#
#
# #dev.new(width=8,height=5)
# multiplot(Fig5a,Fig5b,cols=2,widths=c(1.75,1.75))
## ----Wambaugh2019.Fig6, eval = execute.vignette-------------------------------
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css95_uv:=calc_mc_css(chem.cas=CAS,
# clint.pvalue.threshold=1, #Don't need to use p-values here (Bayesian)
# output.units="uM",
# model="3compartmentss",
# species="Human"),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css95_v:=calc_mc_css(chem.cas=CAS,
# clint.pvalue.threshold=1, #Don't need to use p-values here (Bayesian)
# output.units="uM",
# clint.meas.cv=NULL,
# fup.meas.cv=NULL,
# clint.pop.cv=0.3,
# fup.pop.cv=0.3,
# model="3compartmentss",
# species="Human"),
# by=.(CAS)]
#
#
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css95_u:=cssfun_u(chemcas=CAS,
# MW=MW),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# cssmed_u:=cssfun_u(chemcas=CAS,
# MW=MW,
# probs=0.5),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css25_u:=cssfun_u(chemcas=CAS,
# MW=MW,
# probs=0.025),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css975_u:=cssfun_u(chemcas=CAS,
# MW=MW,
# probs=0.975),
# by=.(CAS)]
#
# ceetox[CAS%in%httk_cas,
# cssmed:=httk::calc_analytic_css(chem.cas=CAS,
# clint.pvalue.threshold=1,
# output.units="uM",
# model="3compartmentss",
# species="Human",
# suppress.messages=TRUE),
# by=.(CAS)]
#
#
#
# # Go back to Bayes estimates:
# chem.physical_and_invitro.data <- HTTK.data
#
# dt_melt <- data.table::melt(ceetox[, .(Compound, CAS, css95_u, css95_v, css95_uv, cssmed)],
# id.vars=c("Compound","CAS", "cssmed"),
# variable.name="type",
# value.name="css95")
#
# dt_melt[, type:=gsub(x=type, pattern="css95_uv", replacement="Both", fixed=TRUE)]
# dt_melt[, type:=gsub(x=type, pattern="css95_u", replacement="Uncertainty", fixed=TRUE)]
# dt_melt[, type:=gsub(x=type, pattern="css95_v", replacement="Variability", fixed=TRUE)]
#
#
# #now set the factor levels explicitly
# dt_melt[, type:=factor(type,
# levels=c("Uncertainty", "Variability", "Both"))]
#
#
# dt_melt[, Norm:=css95/cssmed]
#
#
# #now set the factor levels explicitly
# dt_melt[, Compound:=factor(Compound,
# levels=dt_melt[type=="Both",Compound][order(dt_melt[type=="Both",Norm])])]
#
#
# Fig6 <- ggplot(data=dt_melt) +
# geom_point(size=2,alpha=0.7,aes(x=Compound,
# y=Norm,
# color=type,
# shape=type)) +
# scale_y_log10(label=scientific_10,limits=c(0.9,500)) +
# ylab(expression(paste("Ratio of ",C[ss]," ",95^th," Percentile to Median Estimate"))) +
# xlab("Chemicals")+
# scale_colour_discrete(name=expression(paste(C[ss]," Varied to Reflect")))+
# scale_shape_discrete(name=expression(paste(C[ss]," Varied to Reflect"))) +
# theme_bw() +
# theme(legend.position="bottom",axis.text.x = element_blank())+
# annotate("text", size=8,x = levels(dt_melt$Compound)[length(unique(dt_melt$Compound))/2], y = 400, label = paste("Median Ratio for Uncertainty:",signif(median(subset(dt_melt,type=="Uncertainty")$Norm),3)))+
# annotate("text", size=8,x = levels(dt_melt$Compound)[length(unique(dt_melt$Compound))/2], y = 200, label = paste("Median Ratio for Variability:",signif(median(subset(dt_melt,type=="Variability")$Norm),3)))+
# annotate("text", size=8,x = levels(dt_melt$Compound)[length(unique(dt_melt$Compound))/2], y = 100, label = paste("Median Ratio for Both:",signif(median(subset(dt_melt,type=="Both")$Norm),3)))
#
# #dev.new(width=10,height=6)
# print(Fig6)
## ----Wambaugh2019.Fig7, eval = execute.vignette------------------------------
# # Head to ftp://newftp.epa.gov/COMPTOX/High_Throughput_Screening_Data/Previous_Data/ToxCast_Data_Release_Oct_2015/ to get the ToxCast/Tox21 data:
# # These are the concentrations that caused activity in exess of the background (ACC) "hits".
# #
# # # Load the ACC table into an object Tox21.acc:
# # Tox21.acc <- read.csv("INVITRODB_V2_LEVEL5/EXPORT_LVL5&6_ASID7_TOX21_151020.csv",stringsAsFactors=FALSE)
# # # Subset to the chemicals of interest:
# # Tox21.acc <- subset(Tox21.acc,casn%in%Fig1.table$CAS)
# # # We only need hits:
# # Tox21.acc <- subset(Tox21.acc,hitc==1)
# #
# # # Pare this down to just the data we need:
# # Tox21.acc <- Tox21.acc[,c("code","chid","chnm","casn","aenm","modl_acc","flags")]
# #
# # # In this vignette we use the precalculated quantiles of the ACC's
# # # to help keep the HTTK package smaller:
# #
# # wambaugh2019.tox21 <- NULL
# # for (this.chem in sort(unique(Tox21.acc$casn)))
# # {
# # this.subset <- subset(Tox21.acc,casn==this.chem)
# # this.row <- data.frame(casn=this.chem,
# # med.conc =
# # median(10^(this.subset$modl_acc)),
# # q10.conc =
# # quantile(10^(this.subset$modl_acc),0.1),
# # low.conc =
# # min(10^(this.subset$modl_acc)),
# # q90.conc =
# # quantile(10^(this.subset$modl_acc),0.9),
# # high.conc =
# # max(10^(this.subset$modl_acc)),
# # stringsAsFactors = F)
# # wambaugh2019.tox21 <- rbind(wambaugh2019.tox21, this.row)
# # }
#
# chem.preds <- httk::wambaugh2019.seem3
# directnhanes.preds <- httk::wambaugh2019.nhanes
# Tox21.acc.numeric <- httk::wambaugh2019.tox21
#
# #Subset to those chemicals with HTS hits:
# human.hits <- subset(ceetox,CAS%in%Tox21.acc.numeric$casn)
#
# # Now calculate the oral equivalent dose for each chemical:
# # Must make sure that CSS is in units of uM!!
# for (this.chem in human.hits$CAS)
# {
# med.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"med.conc"]
# q10.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"q10.conc"]
# low.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"low.conc"]
# q90.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"q90.conc"]
# high.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"high.conc"]
# css95.v <- human.hits[human.hits$CAS==this.chem][["css95_v"]]
# css95.uv <- human.hits[human.hits$CAS==this.chem][["css95_uv"]]
# human.hits[human.hits$CAS==this.chem,"HTS.Median.ACC"] <- med.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Q10.ACC"] <- q10.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Q90.ACC"] <- q90.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Low.ACC"] <- low.conc
# human.hits[human.hits$CAS==this.chem,"HTS.High.ACC"] <- high.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Median.equivdose.v"] <- med.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Q10.equivdose.v"] <- q10.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Q90.equivdose.v"] <- q90.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Low.equivdose.v"] <- low.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.High.equivdose.v"] <- high.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Median.equivdose.uv"] <- med.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.Q10.equivdose.uv"] <- q10.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.Q90.equivdose.uv"] <- q90.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.Low.equivdose.uv"] <- low.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.High.equivdose.uv"] <- high.conc/css95.uv
# }
#
#
#
# #Add the exposure predictions:
# for (this.chem in human.hits$CAS)
# {
# print(this.chem)
# # If we have direct inferrences from NHANES, use those exposures:
# if (this.chem %in% directnhanes.preds$CASRN)
# {
# human.hits[human.hits$CAS==this.chem,"Exposure.median"] <- directnhanes.preds[directnhanes.preds$CASRN==this.chem,"lP"]
# human.hits[human.hits$CAS==this.chem,"Exposure.high"] <- directnhanes.preds[directnhanes.preds$CASRN==this.chem,"lP.max"]
# # Otherwise use the heuristics model (Wambaugh 2014):
# } else if (this.chem %in% chem.preds$CAS) {
# human.hits[human.hits$CAS==this.chem,"Exposure.median"] <- chem.preds[chem.preds$CAS==this.chem,"seem3"]
# human.hits[human.hits$CAS==this.chem,"Exposure.high"] <- chem.preds[chem.preds$CAS==this.chem,"seem3.u95"]
# }
# }
#
# # Drop chemicals without exposure predictions:
# human.hits<- subset(human.hits,!is.na(Exposure.median))
#
#
# # This code puts the chemicals into order by margin between exposure and activity:
# chem.names <- unique(human.hits$Compound)
# potency <- rep(Inf,times=length(chem.names))
# names(potency) <- chem.names
# potency.u <- potency
# for (this.chem in chem.names)
# {
# this.subset <- subset(human.hits,Compound==this.chem)
# potency[this.chem] <- this.subset$HTS.Q10.equivdose.v/this.subset$Exposure.high
# potency.u[this.chem] <- this.subset$HTS.Q10.equivdose.uv/this.subset$Exposure.high
# }
# potency <- sort(potency)
# human.hits$Compound <- factor(human.hits$Compound, levels=names(potency))
#
# potency.u <- potency.u[names(potency.u)[potency.u<1]]
# potency.u <- potency.u[names(potency)[potency>1]]
#
# new.overlaps <- names(potency.u[!is.na(potency.u)])
# first.change <- new.overlaps[1]
# last.change <- new.overlaps[length(new.overlaps)]
# window.start <- names(potency)[which(names(potency)==first.change)-5]
# window.end <- names(potency)[which(names(potency)==last.change)+5]
# window.names <- names(potency)[(which(names(potency)==first.change)-5):(which(names(potency)==last.change)+5)]
# # Initialize a new plot:
# Fig7a <- ggplot(data=human.hits) +
# geom_segment(aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.v,yend=HTS.Q90.equivdose.v),size=1,color="white",alpha=0.5) +
# geom_rect(mapping=aes(xmin=window.start,xmax=window.end,ymin=10^-9,ymax=10^3),color="grey",alpha=0.5)+
# geom_segment(aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.v,yend=HTS.Q90.equivdose.v),size=1,color="red",alpha=0.5) +
# geom_point(aes(x=Compound,y=HTS.Median.equivdose.v),shape=3) +
# geom_point(aes(x=Compound,y=HTS.Q10.equivdose.v)) +
# theme_bw() +
# theme(axis.text.x = element_blank()) +
# theme(axis.title.x = element_text(size=16)) +
# theme(axis.title.y = element_text(size=16)) +
# scale_y_log10(label=scientific_10,limits = c(10^-9,10^3)) +
# geom_segment(aes(x=Compound,xend=Compound,y=Exposure.high,yend=Exposure.median),size=1,color="blue",alpha=0.5) +
# geom_point(aes(x=Compound,y=Exposure.median),shape=2) +
# ylab("Bioactive Dose & Exposure\nmg/kg BW/day\n(Variability Only)") +
# xlab("Chemicals Ranked By Ratio Between Bioactive Dose and Exposure")+
# annotate("text", size=8,x = names(potency)[10], y = 100, label = "a")
#
# # Initialize a new plot:
# Fig7c <- ggplot() +
# geom_segment(data=human.hits,aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.v,yend=HTS.Q90.equivdose.v),size=3,color="grey") +
# geom_segment(data=human.hits,aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.uv,yend=HTS.Q90.equivdose.uv),size=3,color="red",alpha=0.5) +
# geom_point(data=human.hits,aes(x=Compound,y=HTS.Median.equivdose.uv),shape=3,size=2) +
# geom_point(data=human.hits,aes(x=Compound,y=HTS.Q10.equivdose.uv),size=2) +
# theme_bw() +
# theme(axis.text.x = element_text(angle = 60, hjust = 1,size=6)) +
# theme(axis.text.x = element_blank()) +
# theme(axis.title.x = element_text(size=16)) +
# theme(axis.title.y = element_text(size=16)) +
# xlim(window.names) +
# scale_y_log10(label=scientific_10,limits = c(10^-8,2*10^2)) +
# geom_segment(data=human.hits,aes(x=Compound,xend=Compound,y=Exposure.high,yend=Exposure.median),size=3,color="blue",alpha=0.5) +
# geom_segment(data=subset(human.hits,Compound%in%names(potency.u)),aes(x=Compound,xend=Compound,y=10^-8,yend=Exposure.median/1.5),size=2,arrow=arrow(length=unit(0.5,"cm")))+
# geom_point(data=human.hits,aes(x=Compound,y=Exposure.median),shape=2,size=2) +
# xlab("") +
# ylab("Bioactive Dose & Exposure\nmg/kg BW/day\n(Uncertainty and Variability)") +
# annotate("text", size=8,x = window.names[1], y = 100, label = "b")
#
#
# #print(Fig10b)
# #dev.new(width=10,height=6)
# multiplot(Fig7a,Fig7c,cols=1,heights=c(0.5,0.5))
#
# write.csv(subset(human.hits,Exposure.high>HTS.Q10.equivdose.uv)[,c("Compound","DSSTox_Substance_Id","Human.Clint","Human.Funbound.plasma","HTS.Q10.ACC","HTS.Q10.equivdose.v","HTS.Q10.equivdose.uv","Exposure.high")],file=paste("SupTable5-",Sys.Date(),".txt",sep=""),row.names=FALSE)
## ----calc.stats, eval = execute.vignette-------------------------------------
# Wetmore.chems <- subset(chem.physical_and_invitro.data,regexpr("Wetmore",Human.Clint.Reference)!=-1)[,c("Compound","CAS","DTXSID","Formula")]
# Wetmore.chems <- subset(Wetmore.chems,CAS%in%get_cheminfo())
#
# for (this.cas in Wetmore.chems$CAS)
# {
# Wetmore.chems[Wetmore.chems$CAS==this.cas,"Css.Lit"] <- get_lit_css(
# chem.cas=this.cas,
# output.units="uM")
# set.seed(42)
# Wetmore.chems[Wetmore.chems$CAS==this.cas,"Css.v192"] <- calc_mc_css(
# chem.cas=this.cas,
# fup.meas.cv=NULL,
# clint.meas.cv=NULL,
# fup.censored.dist=TRUE,
# output.units="uM")
# set.seed(42)
# Wetmore.chems[Wetmore.chems$CAS==this.cas,"Css.v110"] <- calc_mc_css(
# chem.cas=this.cas,
# output.units="uM")
# }
#
# median(Wetmore.chems$Css.v110/Wetmore.chems$Css.Lit,na.rm=TRUE)
#
# ceetox$Human.Clint2 <- sapply(ceetox$Human.Clint,function(x) as.numeric(strsplit(x,",")[[1]][1]))
# ceetox$Human.Fup2 <- sapply(ceetox$Human.Funbound.plasma,function(x) as.numeric(strsplit(x,",")[[1]][1]))
# write.csv(ceetox[,c(1,2,8,9,16,17,6,7,12,30,10,11,15,31,13,14,18:28)],row.names=FALSE,file=paste("Supplemental-Table4-",Sys.Date(),".txt",sep=""))
#
# print(paste("Complete HTTK data on",sum(get_cheminfo()%in%ceetox$CAS),"new chemicals."))
Try the httk package in your browser
Any scripts or data that you put into this service are public.
httk documentation built on March 7, 2023, 7:26 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 = '#>', fig.width=5, fig.height=3)
## ----clear_memory, eval = TRUE------------------------------------------------
rm(list=ls())
## ----runchunks, eval = TRUE---------------------------------------------------
# Set whether or not the following chunks will be executed (run):
execute.vignette <- FALSE
## ----setup, eval = execute.vignette-------------------------------------------
# library(data.table)
# library(httk)
# library(ggplot2)
# library(stringr)
# library(scales)
# library(grid)
# library(reshape)
## ----lmr2, eval = execute.vignette--------------------------------------------
# lm_R2 <- function(m,prefix=NULL){
# out <- substitute("RMSE = "~mse*","~~italic(R)^2~"="~r2,
# list(mse = signif(mean(m$residuals^2)^(1/2),3),
# r2 = format(summary(m)$adj.r.squared, digits = 2)))
# paste(prefix,as.character(as.expression(out)))
# }
## ----multi.plot, eval = execute.vignette--------------------------------------
# multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL, heights=NULL, widths=unit(rep_len(1, cols), "null")) {
#
# # Make a list from the ... arguments and plotlist
# plots <- c(list(...), plotlist)
#
# numPlots = length(plots)
#
# # If layout is NULL, then use 'cols' to determine layout
# if (is.null(layout)) {
# # Make the panel
# # ncol: Number of columns of plots
# # nrow: Number of rows needed, calculated from # of cols
# layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
# ncol = cols, nrow = ceiling(numPlots/cols))
# }
#
# if (numPlots==1) {
# print(plots[[1]])
#
# } else {
# # Set up the page
# grid.newpage()
# if (!is.null(heights))
# {
# pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout),heights=heights,widths=widths)))
# } else {
# pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout),widths=widths)))
# }
# # Make each plot, in the correct location
# for (i in 1:numPlots) {
# # Get the i,j matrix positions of the regions that contain this subplot
# matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
#
# print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
# layout.pos.col = matchidx$col))
# }
# }
# }
## ----scientific.notation, 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))
# }
## ----scientific.notation2, eval = execute.vignette----------------------------
# scientific_10_skip <- function(x) {
# out <- gsub("1e", "10^", scientific_format()(x))
# out <- gsub("\\+","",out)
# out <- gsub("10\\^01","10",out)
# out <- gsub("10\\^00","1",out)
# out <- gsub("10\\^-01",NA,out)
# return(parse(text=out))
# }
## ----cssfun.u, eval = execute.vignette----------------------------------------
# cssfun_u <- function(chemcas,
# MW,
# hepatic.bioavailability,
# Fgutabs=1,
# probs=0.95,
# minimum.Funbound.plasma=0.0001)
# {
# # Create a data table with uncertainty only:
# pop_u_httk <- as.data.table(parameterize_steadystate(chem.cas=chemcas,
# clint.pvalue.threshold=1,
# minimum.Funbound.plasma=0))
#
# # Get the physico-chemical properties:
# params.schmitt <-parameterize_schmitt(chem.cas=chemcas)
# psd <- params.schmitt$pKa_Donor
# psa <- params.schmitt$pKa_Accept
# params.schmitt <- params.schmitt[regexpr("pKa",names(params.schmitt))==-1]
# pop_u_httk[,names(params.schmitt):=params.schmitt]
# pop_u_httk[,pKa_Donor := paste(psd, collapse=",")]
# pop_u_httk[,pKa_Accept := paste(psa, collapse=",")]
# pop_u_httk <- pop_u_httk[rep(1,1000)]
#
# # Calculate the distribution coefficient:
# ions <- calc_ionization(pH=7.4,parameters=pop_u_httk)
# pop_u_httk[,Dow74 := Pow * (ions$fraction_neutral +
# 0.001 * ions$fraction_charged + ions$fraction_zwitter)]
#
# # Parse the Funbound.plasma.dist into quantiles:
# Funbound.plasma.dist <- pop_u_httk[1,Funbound.plasma.dist]
# if (nchar(Funbound.plasma.dist) -
# nchar(gsub(",","",Funbound.plasma.dist))!=2)
# {
# stop("Funbound.plasma distribution should be three values (median,low95th,high95th) separated by commas.")
# }
# temp <- strsplit(Funbound.plasma.dist,",")
# ppb.median <- as.numeric(temp[[1]][1])
# ppb.low95 <- as.numeric(temp[[1]][2])
# ppb.high95 <- as.numeric(temp[[1]][3])
#
# if (ppb.median>minimum.Funbound.plasma)
# {
# # Use optim to estimate alpha and beta for fup such that the median and 95%
# # credible interval approximate the estimate from MCMC:
# ppb.fit <- optim(c(2,(1-ppb.median)/ppb.median*2),
# function(x) (0.95-pbeta(ppb.high95,x[1],x[2])+
# pbeta(ppb.low95,x[1],x[2]))^2+
# (ppb.median-qbeta(0.5,x[1],x[2]))^2,
# method="BFGS")
# # We are drawing new values for the unadjusted Fup:
# pop_u_httk[, unadjusted.Funbound.plasma:=rbeta(1000,
# ppb.fit$par[1],
# ppb.fit$par[2])]
# } else if (ppb.high95>minimum.Funbound.plasma)
# {
# # Assume that since the median is zero but the u95 is not, that there is
# # an exponential distribution:
# # 97.5% of clearance values will be below Funbound.plasma.u95:
# pop_u_httk[,unadjusted.Funbound.plasma:=runif(1000,
# minimum.Funbound.plasma,
# min(1,minimum.Funbound.plasma+
# 2*(ppb.high95-minimum.Funbound.plasma)))]
# pop_u_httk[as.logical(rbinom(1000,1,.95)),
# unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# } else {
# pop_u_httk[, unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# }
#
# # then we need to adjust for in vitro binding:
# pop_u_httk[,Flipid:=subset(physiology.data,
# Parameter=='Plasma Effective Neutral Lipid Volume Fraction')[,
# which(colnames(physiology.data) == 'Human')]]
# pop_u_httk[,Funbound.plasma.adjustment:=1 / (Dow74 * Flipid +
# 1 / unadjusted.Funbound.plasma)/unadjusted.Funbound.plasma]
# pop_u_httk[,
# Funbound.plasma:=unadjusted.Funbound.plasma*Funbound.plasma.adjustment]
#
# # Enforce a minimum ppb:
# pop_u_httk[Funbound.plasma<minimum.Funbound.plasma,
# Funbound.plasma:=minimum.Funbound.plasma]
#
# # If Fup varies, then Rb2p and hepatic bioavailability also vary:
# Rblood2plasma <- get_rblood2plasma(chem.cas=chemcas)
# if (is.na(Rblood2plasma))
# {
# pop_u_httk[, Rblood2plasma:=calc_rblood2plasma(params=pop_u_httk)]
# # Right now we don't similate variability on a known blood:plasma ratio
# } else pop_u_httk[, Rblood2plasma:=Rblood2plasma]
#
# # Parse the clint.dist to retrieve the quantiles:
# Clint.dist <- pop_u_httk[1,Clint.dist]
# if (nchar(Clint.dist) -
# nchar(gsub(",","",Clint.dist))!=3)
# {
# stop("Clint distribution should be four values (median,low95th,high95th,pValue) separated by commas.")
# }
# temp <- strsplit(Clint.dist,",")
# Clint.median <- as.numeric(temp[[1]][1])
# Clint.low95 <- as.numeric(temp[[1]][2])
# Clint.high95 <- as.numeric(temp[[1]][3])
# Clint.pvalue <- as.numeric(temp[[1]][4])
#
# # Use optim to estimate mean and standard deviation of Clint such that the
# # median and 95% credible interval approximate the estimate from MCMC:
# if (Clint.high95>0)
# {
# if (Clint.median>0)
# {
# clint.fit <- optim(c(log(Clint.median),0.3), function(x)
# (0.95-plnorm(Clint.high95,x[1],x[2])+plnorm(Clint.low95,x[1],x[2]))^2+
# (Clint.median-qlnorm(0.5,x[1],x[2]))^2)
# pop_u_httk[,Clint:=rlnorm(1000,clint.fit$par[1],clint.fit$par[2])]
# } else {
# pop_u_httk[,Clint:=exp(runif(1000,log(1),
# (log(Clint.high95)-log(1))/0.975))]
# }
# } else pop_u_httk[,Clint:=0]
# pop_u_httk[as.logical(rbinom(1000,1,Clint.pvalue)),Clint:=0] # the Bayesian "p-value" here reflects how often there is no clearance
#
# pop_u_httk[, CLh:=calc_hepatic_clearance(chem.cas=chemcas,
# parameters=pop_u_httk,hepatic.model='unscaled',suppress.messages=TRUE,clint.pvalue.threshold=1)]
# pop_u_httk[, Qliver:=Qtotal.liverc*BW^0.75]
# pop_u_httk[, hepatic.bioavailability:= Qliver / (Qliver +
# Funbound.plasma * Clint / Rblood2plasma)]
# pop_u_httk[, (setdiff(names(pop_u_httk),
# names(httk::parameterize_steadystate(chem.cas=chemcas)))):=NULL]
#
#
# #compute Css for each individual
# # calc_analytic_css is vectorized so that this will work.
# css <- httk::calc_analytic_css(parameters=pop_u_httk,
# clint.pvalue.threshold=1,
# daily.dose=1,
# clint.pop.cv=NULL,
# fup.pop.cv=NULL,
# output.units="uM",
# model="3compartmentss",
# species="Human",
# suppress.messages=TRUE)#,
# # recalc.blood2plasma=TRUE)
# #get Css95
# cssquants <- quantile(css, probs=probs)
# return(cssquants)
# }
## ----hlfun.u, eval = execute.vignette-----------------------------------------
# hlfun_u <- function(chemcas,
# MW,
# probs=0.95,
# minimum.Funbound.plasma=0.0001)
# {
# # Create a data table with uncertainty only:
# print(chemcas)
# pop_u_httk <- as.data.table(parameterize_steadystate(chem.cas=chemcas,
# clint.pvalue.threshold=1,
# minimum.Funbound.plasma=0))
# params.schmitt <-parameterize_schmitt(chem.cas=chemcas)
# psd <- params.schmitt$pKa_Donor
# psa <- params.schmitt$pKa_Accept
# params.schmitt <- params.schmitt[regexpr("pKa",names(params.schmitt))==-1]
# pop_u_httk[,names(params.schmitt):=params.schmitt]
# pop_u_httk[,pKa_Donor := paste(psd, collapse=",")]
# pop_u_httk[,pKa_Accept := paste(psa, collapse=",")]
# pop_u_httk <- pop_u_httk[rep(1,1000)]
#
# ions <- calc_ionization(pH=7.4,parameters=pop_u_httk)
# pop_u_httk[,Dow74 := Pow * (ions$fraction_neutral + 0.001 * ions$fraction_charged + ions$fraction_zwitter)]
#
# # Parse the Funbound.plasma.dist into quantiles:
# Funbound.plasma.dist <- pop_u_httk[1,Funbound.plasma.dist]
# print(paste(chemcas,Funbound.plasma.dist))
# if (nchar(Funbound.plasma.dist) -
# nchar(gsub(",","",Funbound.plasma.dist))!=2)
# {
# stop("Funbound.plasma distribution should be three values (median,low95th,high95th) separated by commas.")
# }
# temp <- strsplit(Funbound.plasma.dist,",")
# ppb.median <- as.numeric(temp[[1]][1])
# ppb.low95 <- as.numeric(temp[[1]][2])
# ppb.high95 <- as.numeric(temp[[1]][3])
#
# if (ppb.median>minimum.Funbound.plasma)
# {
# # Use optim to estimate alpha and beta for fup such that the median and 95%
# # credible interval approximate the estimate from MCMC:
# ppb.fit <- optim(c(2,(1-ppb.median)/ppb.median*2),
# function(x) (0.95-pbeta(ppb.high95,x[1],x[2])+
# pbeta(ppb.low95,x[1],x[2]))^2+
# (ppb.median-qbeta(0.5,x[1],x[2]))^2,
# method="BFGS")
# # We are drawing new values for the unadjusted Fup:
# pop_u_httk[, unadjusted.Funbound.plasma:=rbeta(1000,
# ppb.fit$par[1],
# ppb.fit$par[2])]
# } else if (ppb.high95>minimum.Funbound.plasma)
# {
# # Assume that since the median is zero but the u95 is not, that there is
# # an exponential distribution:
# # 97.5% of clearance values will be below Funbound.plasma.u95:
# pop_u_httk[, unadjusted.Funbound.plasma:=runif(1000,
# minimum.Funbound.plasma,
# min(1,(ppb.high95-minimum.Funbound.plasma)/0.975))]
# pop_u_httk[as.logical(rbinom(1000,1,.975)),
# unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# } else {
# pop_u_httk[, unadjusted.Funbound.plasma:=minimum.Funbound.plasma]
# }
#
# # then we need to adjust for in vitro binding:
# pop_u_httk[,Flipid:=subset(physiology.data,
# Parameter=='Plasma Effective Neutral Lipid Volume Fraction')[,
# which(colnames(physiology.data) == 'Human')]]
# pop_u_httk[,Funbound.plasma.adjustment:=1 / (Dow74 * Flipid +
# 1 / unadjusted.Funbound.plasma)/unadjusted.Funbound.plasma]
# pop_u_httk[,
# Funbound.plasma:=unadjusted.Funbound.plasma*Funbound.plasma.adjustment]
#
# # Enforce a minimum ppb:
# pop_u_httk[Funbound.plasma<minimum.Funbound.plasma,
# Funbound.plasma:=minimum.Funbound.plasma]
#
# Rblood2plasma <- get_rblood2plasma(chem.cas=chemcas)
# if (is.na(Rblood2plasma))
# {
# pop_u_httk[, Rblood2plasma:=calc_rblood2plasma(params=pop_u_httk)]
# } else pop_u_httk[, Rblood2plasma:=Rblood2plasma] # Right now we don't similate variability on a known blood:plasma ratio
# # Parse the clint.dist to retrieve the quantiles:
# Clint.dist <- pop_u_httk[1,Clint.dist]
# if (nchar(Clint.dist) -
# nchar(gsub(",","",Clint.dist))!=3)
# {
# stop("Clint distribution should be four values (median,low95th,high95th,pValue) separated by commas.")
# }
# temp <- strsplit(Clint.dist,",")
# Clint.median <- as.numeric(temp[[1]][1])
# Clint.low95 <- as.numeric(temp[[1]][2])
# Clint.high95 <- as.numeric(temp[[1]][3])
# Clint.pvalue <- as.numeric(temp[[1]][4])
#
# # Use optim to estimate mean and standard deviation of Clint such that the
# # median and 95% credible interval approximate the estimate from MCMC:
# if (Clint.high95>0)
# {
# if (Clint.median>0)
# {
# clint.fit <- optim(c(log(Clint.median),0.3), function(x)
# (0.95-plnorm(Clint.high95,x[1],x[2])+plnorm(Clint.low95,x[1],x[2]))^2+
# (Clint.median-qlnorm(0.5,x[1],x[2]))^2)
# pop_u_httk[,Clint:=rlnorm(1000,clint.fit$par[1],clint.fit$par[2])]
# } else {
# pop_u_httk[,Clint:=exp(runif(1000,log(1),
# (log(Clint.high95)-log(1))/0.975))]
# }
# } else pop_u_httk[,Clint:=0]
# pop_u_httk[as.logical(rbinom(1000,1,Clint.pvalue)),Clint:=0] # the Bayesian "p-value" here reflects how often there is no clearance
#
# pop_u_httk[, CLh:=calc_hepatic_clearance(chem.cas=chemcas,parameters=pop_u_httk,hepatic.model='unscaled',suppress.messages=TRUE,clint.pvalue.threshold=1)]
# pop_u_httk[, Qliver:=Qtotal.liverc*BW^0.75]
# pop_u_httk[, hepatic.bioavailability:= Qliver / (Qliver + Funbound.plasma * Clint / Rblood2plasma)]
#
# PC.names <- names(predict_partitioning_schmitt(chem.cas="80-05-7"))
# pop_u_httk<-cbind(pop_u_httk,as.data.table(predict_partitioning_schmitt(parameters=pop_u_httk)))
# pop_u_httk[,hl:=calc_elimination_rate(parameters=pop_u_httk)]
#
# return(log(2)/quantile(unlist(pop_u_httk[,"hl",with=FALSE]),probs=probs))
# }
## ----Wambaugh2019.Fig1, eval = execute.vignette-------------------------------
# Fig1.table <- httk::wambaugh2019.raw
# Fig1.table$Error1 <- Fig1.table$Base.Fup.High - Fig1.table$Base.Fup.Low
# Fig1.table$Error2 <- Fig1.table$Affinity.Fup.High - Fig1.table$Affinity.Fup.Low
# Fig1.table$Sigma1 <- Fig1.table$Error1/(2*qnorm(0.975))
# Fig1.table$Mean1 <- (Fig1.table$Base.Fup.High + Fig1.table$Base.Fup.Low)/(qnorm(0.975))
# Fig1.table$CV1 <- Fig1.table$Sigma1/Fig1.table$Mean1
# Fig1.table$Sigma2 <- Fig1.table$Error2/(2*qnorm(0.975))
# Fig1.table$Mean2 <- (Fig1.table$Base.Fup.High + Fig1.table$Base.Fup.Low)/(qnorm(0.975))
# Fig1.table$CV2 <- Fig1.table$Sigma2/Fig1.table$Mean2
# Fig1.table$ErrorClint <- Fig1.table$CLint.1uM.High95th - Fig1.table$CLint.1uM.Low95th
# Fig1.table$SigmaClint <- Fig1.table$ErrorClint/(2*qnorm(0.975))
# Fig1.table$MeanClint <- (Fig1.table$CLint.1uM.High95th + Fig1.table$CLint.1uM.Low95th)/(qnorm(0.975))
# Fig1.table$CVClint <- Fig1.table$SigmaClint/Fig1.table$MeanClint
#
# Fupmeasured <- subset(Fig1.table,!is.na(Affinity.Fup.Med))
# Fupmeasured <- subset(Fupmeasured,!duplicated(CAS))
# CLintmeasured <- subset(Fig1.table,!is.na(CLint.1uM.Median))
# CLintmeasured <- subset(CLintmeasured,!duplicated(CAS))
#
#
#
# print(paste("New HTTK experimental measurements were made for",length(unique(Fig1.table$CAS)),"chemicals from the ToxCast HTS library."))
#
# print(paste("Fup was successfully measured for",length(unique(subset(Fig1.table,!is.na(Affinity.Fup.Med))$CAS)),"chemicals from the ToxCast HTS library."))
#
# print(paste("CLint was successfully measured for",length(unique(subset(Fig1.table,!is.na(CLint.1uM.Median))$CAS)),"chemicals from the ToxCast HTS library."))
#
#
#
#
# Fig1a <- ggplot(Fig1.table)+
# geom_freqpoly(binwidth = 0.5,lwd=1.5,color="Red",aes(Base.Fup.Med))+
# geom_freqpoly(binwidth = 0.5,lwd=1.5,lty=3,color="Blue",aes(Affinity.Fup.Med))+
# xlab(expression(paste("Measured ",F[up]))) +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(10^-6,1.05))+
# labs(title=paste(length(unique(subset(Fig1.table,!is.na(Affinity.Fup.Med))$CAS)),"Chemicals Measured"))+
# theme_bw()+
# theme( text = element_text(size=10)) +
# annotate("text", x=10^-6,y=80,size=8,label="a")
#
# print(paste("Median 100%:",signif(median(Fig1.table$Base.Fup.Med,na.rm=TRUE),2)))
# print(paste("Median Titration:",signif(median(Fig1.table$Affinity.Fup.Med,na.rm=TRUE),2)))
#
#
#
# Fig1b <- ggplot(Fig1.table)+
# geom_histogram(binwidth = 0.5,alpha=0.6,fill="green",aes(CLint.1uM.Median))+
# xlab(expression(paste("Measured ",Cl[int]," (uL/min/",10^6," hepatocytes)"))) +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(10^-1,5*10^3))+
# labs(title=paste(length(unique(subset(Fig1.table,!is.na(CLint.1uM.Median))$CAS)),"Chemicals Measured"))+
# theme_bw()+
# theme( text = element_text(size=10)) +
# annotate("text", x=2*10^-1,y=80,size=8,label="b")
#
# print(paste(sum(Fig1.table$CLint.1uM.Median==0,na.rm=TRUE),"Zero Measurments"))
# print(paste("Median Non-Zero:",signif(median(Fig1.table$CLint.1uM.Median,na.rm=TRUE),3)))
#
#
#
# Fig1c <- ggplot(Fig1.table)+
# geom_histogram(binwidth = 0.1,alpha=0.2,fill="Red",aes(Error1))+
# geom_histogram(binwidth = 0.1,alpha=0.2,fill="Blue",aes(Error2))+
# xlab("Width of Credible Interval") +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(5*10^-4,1.05))+
# labs(title=paste("CI:",signif(median(Fig1.table$Error1,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CV1,na.rm=TRUE),1),"/ CI:",signif(median(Fig1.table$Error2,na.rm=TRUE),1),"CV:",signif(median(Fig1.table$CV2,na.rm=TRUE),1)))+
# theme_bw()+
# theme( text = element_text(size=10))+
# annotate("text", x=5*10^-4,y=42,size=8,label="c")
#
# print(paste("Median 100% CI:",signif(median(Fig1.table$Error1,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CV1,na.rm=TRUE),2)))
# print(paste("Median Titr. CI:",signif(median(Fig1.table$Error2,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CV2,na.rm=TRUE),2)))
#
#
# Fig1d <- ggplot(Fig1.table)+
# geom_histogram(binwidth = 0.1,alpha=0.6,fill="green",aes(ErrorClint))+
# xlab("Width of Credible Interval") +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10,limits=c(10^0,10^3))+
# labs(title=paste("Median CI:",signif(median(Fig1.table$ErrorClint,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CVClint,na.rm=TRUE),2)))+
# theme_bw()+
# theme( text = element_text(size=10))+
# annotate("text", x=1.5,y=25,size=8,label="d")
#
# print(paste("Median CI:",signif(median(Fig1.table$ErrorClint,na.rm=TRUE),2),"CV:",signif(median(Fig1.table$CVClint,na.rm=TRUE),2)))
#
# multiplot(Fig1a,Fig1c,Fig1b,Fig1d,cols=2,widths=c(1.75,1.75))
## ----Wambaugh2019.Fig2, eval = execute.vignette-------------------------------
# Fig2.table <- httk::wambaugh2019.raw
# print(paste(sum(!is.na(Fig2.table$Fup.point)),"total successful Fup estimates using point estimation."))
# print(paste(sum(!is.na(Fig2.table$Base.Fup.Med)),"total successful Fup estimates using BAyesian analysis of top protein concentration."))
# Fig2.table$Error1 <- Fig2.table$Base.Fup.High - Fig2.table$Base.Fup.Low
# Fig2.table$Error2 <- Fig2.table$Affinity.Fup.High - Fig2.table$Affinity.Fup.Low
# Fig2.table$Sigma1 <- Fig2.table$Error1/(2*qnorm(0.975))
# Fig2.table$Mean1 <- (Fig2.table$Base.Fup.High + Fig2.table$Base.Fup.Low)/(qnorm(0.975))
# Fig2.table$CV1 <- Fig2.table$Sigma1/Fig2.table$Mean1
# Fig2.table$Sigma2 <- Fig2.table$Error2/(2*qnorm(0.975))
# Fig2.table$Mean2 <- (Fig2.table$Base.Fup.High + Fig2.table$Base.Fup.Low)/(qnorm(0.975))
# Fig2.table$CV2 <- Fig2.table$Sigma2/Fig2.table$Mean2
# Fig2.table$Point.Estimate <- "Good"
# Fig2.table[Fig2.table$Fup.point<0&!is.na(Fig2.table$Fup.point),"Point.Estimate"]<-"<0"
# Fig2.table[Fig2.table$Fup.point>1&!is.na(Fig2.table$Fup.point),"Point.Estimate"]<-">1"
# Fig2.table[is.na(Fig2.table$Fup.point),"Point.Estimate"]<-"No Value"
# Fig2.table[Fig2.table$Point.Estimate%in%c("<0","No Value"),"Fup.point"] <- Fig2.table[Fig2.table$Point.Estimate%in%c("<0","No Value"),"Base.Fup.Med"]
# Fig2.table$Uncertain<-"Yes"
# Fig2.table[is.na(Fig2.table$CV1),"CV1"] <- -999
# Fig2.table[Fig2.table$CV1<0.49,"Uncertain"] <- "No"
# Fig2.table[Fig2.table$CV1==-999,"CV1"] <- NA
#
# Fig2a <- ggplot(subset(Fig2.table,!is.na(Base.Fup.Med)),aes(x=Fup.point,y=Base.Fup.Med,shape=Point.Estimate,alpha=Uncertain,colour=Point.Estimate)) +
# geom_point(size=3)+
# # geom_errorbar(aes(ymax = Base.Fup.High, ymin=Base.Fup.Low))+
# scale_y_log10(label=scientific_10,limits=c(10^-6,2)) +
# scale_x_log10(label=scientific_10,limits=c(10^-6,2)) +
# ylab(expression(paste(F[up]," Bayesian Analysis (Single Conc.)"))) +
# xlab(expression(paste(F[up]," Point Estimate (Traditional Analysis)"))) +
# geom_abline(intercept = 0, slope = 1,linetype="dashed", colour="Blue") +
# # geom_vline(xintercept = 1, size=2,linetype="dashed", colour="red")+
# geom_vline(xintercept = 0.01, size=2,linetype="dotted", colour="red")+
# geom_hline(yintercept = 0.01, size=2,linetype="dotted", colour="red")+
# scale_colour_discrete(name="Point Estimate")+
# scale_shape_discrete(name="Point Estimate") +
# scale_alpha_manual(values=c(1,0.3))+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=12))+
# annotate("text", size=8,x = 10^-6, y = 2, label = "a")+
# guides(alpha=guide_legend(nrow=2,byrow=TRUE),colour=guide_legend(nrow=2,byrow=TRUE),shape=guide_legend(nrow=2,byrow=TRUE))
#
# #print(Fig2a)
#
# print(paste(sum(Fig2.table$Uncertain=="Yes"),"uncertain estimates using single conc."))
# print(paste(sum(Fig2.table$Point.Estimate=="<0"),"chemicals with negative point estimates."))
# print(paste(sum(Fig2.table$Point.Estimate==">1"),"chemicals with Fup > 1."))
#
#
#
# summary(lm(data=subset(Fig2.table,!is.na(Fup.point)&Fup.point>0&Base.Fup.Med>0),log(Fup.point)~log(Base.Fup.Med)))
#
#
#
#
# Fig2.table$Uncertain2<-"Yes"
# Fig2.table[is.na(Fig2.table$CV2),"CV2"] <- -999
# Fig2.table[Fig2.table$CV2<0.49,"Uncertain2"] <- "No"
# Fig2.table[Fig2.table$CV2==-999,"CV2"] <- NA
# Fig2.table$TopFail <- ""
# Fig2.table[!is.na(Fig2.table$Affinity.Fup.Med)&Fig2.table$Uncertain2=="No","TopFail"]<-"No"
# Fig2.table[Fig2.table$TopFail=="No"&!is.na(Fig2.table$Base.Fup.Med)&Fig2.table$Uncertain=="Yes","TopFail"]<-"Single Conc. Only"
# Fig2.table[!is.na(Fig2.table$Affinity.Fup.Med)&Fig2.table$Uncertain2=="Yes","TopFail"]<-"Yes"
# Fig2.table[!is.na(Fig2.table$Affinity.Fup.Med)&is.na(Fig2.table$Base.Fup.Med),"TopFail"] <- "Single Conc. Only"
# Fig2.table[is.na(Fig2.table$Base.Fup.Med),"Base.Fup.Med"] <- Fig2.table[is.na(Fig2.table$Base.Fup.Med),"Affinity.Fup.Med"]
#
#
# Fig2b <- ggplot(subset(Fig2.table,!is.na(Affinity.Fup.Med)),aes(x=Base.Fup.Med,y=Affinity.Fup.Med,shape=TopFail,alpha=TopFail,colour=TopFail)) +
# geom_point(size=3)+
# geom_point(data=subset(Fig2.table,TopFail=="No"),size=3)+
# # geom_errorbar(aes(ymax = Fup.High.x, ymin=Fup.Low.x))+
# # geom_errorbarh(aes(xmin=Fup.Low.y,xmax=Fup.High.y))+
# scale_y_log10(label=scientific_10,limits=c(10^-8,2)) +
# scale_x_log10(label=scientific_10,limits=c(10^-8,2)) +
# xlab(expression(paste(F[up]," Bayesian Analysis (Single Conc.)"))) +
# ylab(expression(paste(F[up]," Bayesian Analysis (Three Conc.)"))) +
# geom_abline(intercept = 0, slope = 1,linetype="dashed", colour="Blue")+
# scale_colour_discrete(name="Uncertain")+
# scale_shape_discrete(name="Uncertain") +
# # annotate("text",size=8, x = 10^-4, y = 1, label = "B")+
# scale_alpha_manual(name="Uncertain",values=c(1,0.6,0.3))+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=12))+
# annotate("text", size=8,x = 10^-8, y = 2, label = "b")+
# guides(alpha=guide_legend(nrow=2,byrow=TRUE),colour=guide_legend(nrow=2,byrow=TRUE),shape=guide_legend(nrow=2,byrow=TRUE))
#
# #print(Fig2b)
#
# # Need to reduce to unique chemicals.
# Fupbaseworked <- subset(Fig2.table,!is.na(Base.Fup.Med)&Uncertain=="No")
# Fupbasenot <- subset(Fig2.table,!is.na(Base.Fup.Med)&Uncertain=="Yes")
# Fupbaseworked <- subset(Fupbaseworked,!duplicated(CAS))
# Fupbasenot <- subset(Fupbasenot,!duplicated(CAS))
# # Also need to count a chemical as a success if it worked for either duplicate:
# Fupbasenot <- subset(Fupbasenot,!(CAS%in%Fupbaseworked$CAS))
#
# Fupaffinityworked <- subset(Fig2.table,!is.na(Affinity.Fup.Med)&Uncertain2=="No")
# Fupaffinitynot <- subset(Fig2.table,!is.na(Affinity.Fup.Med)&Uncertain2=="Yes")
# Fupaffinityworked <- subset(Fupaffinityworked,!duplicated(CAS))
# Fupaffinitynot <- subset(Fupaffinitynot,!duplicated(CAS))
# Fupaffinitynot <- subset(Fupaffinitynot,!(CAS%in%Fupaffinityworked$CAS))
#
#
# print(paste(dim(Fupbaseworked)[1],"total successful Fup estimates using single conc."))
# print(paste(signif(dim(Fupbasenot)[1]/(dim(Fupbaseworked)[1]+dim(Fupbasenot)[1])*100,3),"percent failure rate using single conc."))
# print(paste(dim(Fupaffinityworked)[1],"total successful Fup estimates using protein titration."))
# print(paste(signif(dim(Fupaffinitynot)[1]/(dim(Fupaffinityworked)[1]+dim(Fupaffinitynot)[1])*100,3),"percent failure rate using protein titration."))
# print(paste(dim(Fupbasenot)[1],"uncertain estimates using original protocol."))
# print(paste(dim(Fupaffinitynot)[1],"uncertain estimates using protein titration."))
# print(paste(dim(Fupaffinityworked)[1]-dim(Fupbaseworked)[1],"chemicals that could not be measured at 100% plasma protein."))
#
# Wetmore.data <- httk:::Wetmore.data
# W2015 <- subset(Wetmore.data,Reference=="Wetmore 2015")
#
# W2012 <- subset(Wetmore.data,Reference=="Wetmore 2012")
# print(paste("The Fup assay failed for ",signif(sum(W2012$Fub==0.005)/length(W2012$Fub)*100,3),"% of chemicals in {Wetmore, 2012} and ",signif(sum(W2015$Fub==0)/length(W2015$Fub)*100,3),"% of chemicals in {Wetmore, 2015}."))
#
#
# Fuplod <- subset(Fig2.table,TopFail=="Single Conc. Only")
# print(paste("Among the chemicals that were uncertain using the original protocol but could be measured using the new three concentration protocol, the median estimated Fup was ",signif(median(Fuplod$Affinity.Fup.Med),2),"with values as low as ",signif(min(Fuplod$Affinity.Fup.Med),2),"and as high as",signif(max(Fuplod$Affinity.Fup.Med),2),". Previously, a default of 0.005 has been assumed when Fup could not be measured {Rotroff, 2010}.}"))
#
#
# summary(lm(data=subset(Fig2.table,!is.na(Base.Fup.Med)&Affinity.Fup.Med>0&Base.Fup.Med>0),log(Base.Fup.Med)~log(Affinity.Fup.Med)))
#
# #dev.new(width=10,height=6)
# multiplot(Fig2a,Fig2b,cols=2,widths=c(1.75,1.75))
## ----Wambaugh2019.Fig3, eval = execute.vignette-------------------------------
# Fig1.table$Class <- "Other"
# Fig1.table[Fig1.table$DTXSID %in% pharma$DTXSID,"Class"]<-"Pharmaceutical"
# Fig1.table$Class <- as.factor(Fig1.table$Class)
#
# Fig3 <- ggplot(Fig1.table,aes(Affinity.Kd.Med,fill=Class))+
# geom_histogram(binwidth = 0.5)+
# xlab(expression(paste("Binding Affinity (", mu,"M)",sep=""))) +
# ylab("Number of chemicals")+
# scale_x_log10(label=scientific_10)+
# scale_fill_discrete(name="Chemical Class")+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=18))
#
# # dev.new()
# print(Fig3)
## ----Wambaugh2019.Fig4, eval = execute.vignette-------------------------------
# clearance.table <- subset(httk::wambaugh2019.raw,!is.na(CLint.1uM.Median))
#
# print(paste(sum(clearance.table$CLint.1uM.Median==0),"zero valued median Bayesian posteriors"))
# print(paste(sum(clearance.table$CLint.1uM.Point==0&clearance.table$CLint.1uM.Median>0,na.rm=TRUE),"zero valued point estimates where median posterior>0"))
# print(paste(sum(clearance.table$CLint.1uM.Point>0&clearance.table$CLint.1uM.Median==0,na.rm=TRUE),"zero valued median Bayesian posteriors where point estimate was non-zero"))
# print(paste(sum(is.na(clearance.table$CLint.1uM.Point)&clearance.table$CLint.1uM.Median==0,na.rm=TRUE),"zero valued median Bayesian posteriors where point estimate was NA"))
# print(paste(sum(is.na(clearance.table$CLint.1uM.Point)&clearance.table$CLint.1uM.Median>0,na.rm=TRUE),"non-zero valued median Bayesian posteriors where point estimate was NA"))
#
#
#
#
# clearance.table$Fit <- "Decreasing"
# clearance.table[clearance.table$CLint.1uM.Median == 0,"Fit"] <- "Median Zero"
# #clearance.table[is.na(clearance.table$CLint.1uM.Point == 0),"Fit"] <- "Point Est. Missing"
# for (i in 1:dim(clearance.table)[1])
# if (!is.na(clearance.table[i,"CLint.1uM.Point"]))
# {
# if (clearance.table[i,"CLint.1uM.Point"] == 0) clearance.table[i,"Fit"] <- "Point Est. Zero"
# } else clearance.table[i,"Fit"] <- "Point Est. Failed"
#
#
# set.seed(123456)
# clearance.table[is.na(clearance.table$CLint.1uM.Point),"CLint.1uM.Point"]<-runif(sum(is.na(clearance.table$CLint.1uM.Point)),0.1,0.2)
# clearance.table[clearance.table$CLint.1uM.Point==0,"CLint.1uM.Point"]<-runif(sum(clearance.table$CLint.1uM.Point==0),0.3,0.8)
# clearance.table[clearance.table$CLint.1uM.Low95th==0,"CLint.1uM.Low95th"]<-0.1
# clearance.table[clearance.table$CLint.1uM.Median==0,"CLint.1uM.Median"]<-0.1
# clearance.table[clearance.table$CLint.1uM.High95th==0,"CLint.1uM.High95th"]<-0.1
# clearance.table[clearance.table$CLint.1uM.High95th>1000,"CLint.1uM.High95th"]<-1000
# clearance.table[clearance.table$CLint.1uM.Low95th<0.1,"CLint.1uM.Low95th"]<-0.1
#
# clearance.fit <- lm(log10(CLint.1uM.Median)~log10(CLint.1uM.Point),data=subset(clearance.table,!is.na(CLint.1uM.Point)&CLint.1uM.Point>1))
# Fig4 <- ggplot(clearance.table, aes(x=CLint.1uM.Point,y=CLint.1uM.Median,colour=Fit)) +
# geom_segment(aes(x=CLint.1uM.Point,xend=CLint.1uM.Point,y=CLint.1uM.Low95th,yend=CLint.1uM.High95th),size=1,alpha=0.5)+
# geom_point(size=3) +
# scale_y_log10(label=scientific_10,limits=c(10^-1,1000)) +
# scale_x_log10(label=scientific_10_skip,limits=c(10^-1,1000)) +
# xlab(expression(paste(CL[int]," (",mu,"L/min/",10^{6}," Hep.) Point Estimate",sep="")))+
# ylab(expression(paste(CL[int]," (",mu,"L/min/",10^{6}," Hep.) Bayesian",sep="")))+
# geom_vline(xintercept = 0.25, size=2,colour="black")+
# geom_vline(xintercept = 0.9, size=2,colour="black")+
# geom_abline(intercept = 0, slope = 1,linetype="dashed")+
# scale_colour_discrete(name="Assay Result")+
# theme_bw() +
# theme(legend.position="bottom", text = element_text(size=14))+
# guides(colour=guide_legend(nrow=2,byrow=TRUE))+
# annotate("text",x = 10^2, y = 10^0,size=6, label = lm_R2(clearance.fit,prefix=""),parse=TRUE)
#
# #dev.new()
# print(Fig4)
#
# summary(clearance.fit)
## ----Wambaugh2019.Fig5, eval = execute.vignette-------------------------------
# # This section takes a long time because calc_vdist is not yet data.table compatible:
# ceetox <- as.data.table(subset(httk::wambaugh2019,!is.na(Human.Funbound.plasma)&!is.na(Human.Clint)&!is.na(logP)))
# httk_cas <- ceetox$CAS[ceetox$CAS %in% get_cheminfo(model="3compartmentss")]
#
# # back up chem.physical_and_invtro.data:
# HTTK.data <- chem.physical_and_invitro.data
#
# # Use the point estimates:
# chem.physical_and_invitro.data <- add_chemtable(httk::wambaugh2019.raw,current.table=chem.physical_and_invitro.data,data.list=list(Compound="Name",CAS="CAS",Clint="CLint.1uM.Point",Funbound.plasma="Fup.point"),reference="Wambaugh 2019",species="Human",overwrite=TRUE)
#
#
#
# schmitt_cas <- httk_cas[httk_cas %in% get_cheminfo(model="schmitt")]
# ceetox[CAS%in%schmitt_cas,
# hlpoint:=log(2)/httk::calc_elimination_rate(chem.cas=CAS),
# by=.(CAS)]
#
# chem.physical_and_invitro.data <- HTTK.data
#
# ceetox[CAS%in%schmitt_cas,
# hl_975:=hlfun_u(chemcas=CAS,
# MW=MW,
# probs=0.975),
# by=.(CAS)]
#
# ceetox[CAS%in%schmitt_cas,
# hl_med:=hlfun_u(chemcas=CAS,
# MW=MW,
# probs=0.5),
# by=.(CAS)]
#
# ceetox[CAS%in%schmitt_cas,
# hl_25:=hlfun_u(chemcas=CAS,
# MW=MW,
# probs=0.025),
# by=.(CAS)]
#
# ceetox[,hlpointEstimated:=TRUE]
# ceetox[is.na(hlpoint),hlpointEstimated:=FALSE]
# ceetox[is.na(hlpoint),hlpoint:=hl_med]
#
# hlfit <- lm(data=subset(ceetox,hlpointEstimated),log10(hlpoint)~log10(hl_med))
# Fig5a <- ggplot(data=subset(ceetox,hlpointEstimated),aes(x=hlpoint,y=hl_med)) +
# geom_point(size=3,alpha=0.5) +
# geom_errorbar(aes(ymax = hl_975, ymin=hl_25),alpha=0.3)+
# scale_y_log10(label=scientific_10,limits=c(1,10^6)) +
# scale_x_log10(label=scientific_10,limits=c(1,10^6)) +
# xlab("half-life (h) Point Estimate")+
# ylab("Bayesian half-life (h)")+
# geom_abline(intercept = 0, slope = 1,linetype="dashed",color="blue")+
# annotate("text", size=8,x = 1, y = 10^6, label = "A")+
# theme_bw() +
# theme(text = element_text(size=18))+
# annotate("text",x = 3*10^3, y = 10^0,size=5, label = lm_R2(hlfit,prefix=""),parse=TRUE)
#
#
# cssfit <- lm(data=ceetox,log10(cssmed)~log10(cssmed_u))
# Fig5b <- ggplot(data=ceetox,aes(x=cssmed,y=cssmed_u)) +
# geom_point(size=3,alpha=0.5) +
# geom_errorbar(aes(ymin=css25_u,ymax=css975_u),alpha=0.3)+
# scale_y_log10(label=scientific_10,limits=c(10^-3,10^5)) +
# scale_x_log10(label=scientific_10,limits=c(10^-3,10^5)) +
# xlab(expression(paste(C[ss]," Point Estimate (uM)")))+
# ylab(expression(paste("Bayesian ",C[ss]," (uM)")))+
# annotate("text", size=8,x = 10^-3, y = 10^5, label = "B")+
# theme_bw() +
# theme(text = element_text(size=18)) +
# geom_abline(intercept = 0, slope = 1,linetype="dashed",color="blue")+
# annotate("text",x = 10^2, y = 10^-3,size=5, label = lm_R2(cssfit,prefix=""),parse=TRUE)
#
#
# #dev.new(width=8,height=5)
# multiplot(Fig5a,Fig5b,cols=2,widths=c(1.75,1.75))
## ----Wambaugh2019.Fig6, eval = execute.vignette-------------------------------
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css95_uv:=calc_mc_css(chem.cas=CAS,
# clint.pvalue.threshold=1, #Don't need to use p-values here (Bayesian)
# output.units="uM",
# model="3compartmentss",
# species="Human"),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css95_v:=calc_mc_css(chem.cas=CAS,
# clint.pvalue.threshold=1, #Don't need to use p-values here (Bayesian)
# output.units="uM",
# clint.meas.cv=NULL,
# fup.meas.cv=NULL,
# clint.pop.cv=0.3,
# fup.pop.cv=0.3,
# model="3compartmentss",
# species="Human"),
# by=.(CAS)]
#
#
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css95_u:=cssfun_u(chemcas=CAS,
# MW=MW),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# cssmed_u:=cssfun_u(chemcas=CAS,
# MW=MW,
# probs=0.5),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css25_u:=cssfun_u(chemcas=CAS,
# MW=MW,
# probs=0.025),
# by=.(CAS)]
#
# set.seed(42)
# ceetox[CAS%in%httk_cas,
# css975_u:=cssfun_u(chemcas=CAS,
# MW=MW,
# probs=0.975),
# by=.(CAS)]
#
# ceetox[CAS%in%httk_cas,
# cssmed:=httk::calc_analytic_css(chem.cas=CAS,
# clint.pvalue.threshold=1,
# output.units="uM",
# model="3compartmentss",
# species="Human",
# suppress.messages=TRUE),
# by=.(CAS)]
#
#
#
# # Go back to Bayes estimates:
# chem.physical_and_invitro.data <- HTTK.data
#
# dt_melt <- data.table::melt(ceetox[, .(Compound, CAS, css95_u, css95_v, css95_uv, cssmed)],
# id.vars=c("Compound","CAS", "cssmed"),
# variable.name="type",
# value.name="css95")
#
# dt_melt[, type:=gsub(x=type, pattern="css95_uv", replacement="Both", fixed=TRUE)]
# dt_melt[, type:=gsub(x=type, pattern="css95_u", replacement="Uncertainty", fixed=TRUE)]
# dt_melt[, type:=gsub(x=type, pattern="css95_v", replacement="Variability", fixed=TRUE)]
#
#
# #now set the factor levels explicitly
# dt_melt[, type:=factor(type,
# levels=c("Uncertainty", "Variability", "Both"))]
#
#
# dt_melt[, Norm:=css95/cssmed]
#
#
# #now set the factor levels explicitly
# dt_melt[, Compound:=factor(Compound,
# levels=dt_melt[type=="Both",Compound][order(dt_melt[type=="Both",Norm])])]
#
#
# Fig6 <- ggplot(data=dt_melt) +
# geom_point(size=2,alpha=0.7,aes(x=Compound,
# y=Norm,
# color=type,
# shape=type)) +
# scale_y_log10(label=scientific_10,limits=c(0.9,500)) +
# ylab(expression(paste("Ratio of ",C[ss]," ",95^th," Percentile to Median Estimate"))) +
# xlab("Chemicals")+
# scale_colour_discrete(name=expression(paste(C[ss]," Varied to Reflect")))+
# scale_shape_discrete(name=expression(paste(C[ss]," Varied to Reflect"))) +
# theme_bw() +
# theme(legend.position="bottom",axis.text.x = element_blank())+
# annotate("text", size=8,x = levels(dt_melt$Compound)[length(unique(dt_melt$Compound))/2], y = 400, label = paste("Median Ratio for Uncertainty:",signif(median(subset(dt_melt,type=="Uncertainty")$Norm),3)))+
# annotate("text", size=8,x = levels(dt_melt$Compound)[length(unique(dt_melt$Compound))/2], y = 200, label = paste("Median Ratio for Variability:",signif(median(subset(dt_melt,type=="Variability")$Norm),3)))+
# annotate("text", size=8,x = levels(dt_melt$Compound)[length(unique(dt_melt$Compound))/2], y = 100, label = paste("Median Ratio for Both:",signif(median(subset(dt_melt,type=="Both")$Norm),3)))
#
# #dev.new(width=10,height=6)
# print(Fig6)
## ----Wambaugh2019.Fig7, eval = execute.vignette------------------------------
# # Head to ftp://newftp.epa.gov/COMPTOX/High_Throughput_Screening_Data/Previous_Data/ToxCast_Data_Release_Oct_2015/ to get the ToxCast/Tox21 data:
# # These are the concentrations that caused activity in exess of the background (ACC) "hits".
# #
# # # Load the ACC table into an object Tox21.acc:
# # Tox21.acc <- read.csv("INVITRODB_V2_LEVEL5/EXPORT_LVL5&6_ASID7_TOX21_151020.csv",stringsAsFactors=FALSE)
# # # Subset to the chemicals of interest:
# # Tox21.acc <- subset(Tox21.acc,casn%in%Fig1.table$CAS)
# # # We only need hits:
# # Tox21.acc <- subset(Tox21.acc,hitc==1)
# #
# # # Pare this down to just the data we need:
# # Tox21.acc <- Tox21.acc[,c("code","chid","chnm","casn","aenm","modl_acc","flags")]
# #
# # # In this vignette we use the precalculated quantiles of the ACC's
# # # to help keep the HTTK package smaller:
# #
# # wambaugh2019.tox21 <- NULL
# # for (this.chem in sort(unique(Tox21.acc$casn)))
# # {
# # this.subset <- subset(Tox21.acc,casn==this.chem)
# # this.row <- data.frame(casn=this.chem,
# # med.conc =
# # median(10^(this.subset$modl_acc)),
# # q10.conc =
# # quantile(10^(this.subset$modl_acc),0.1),
# # low.conc =
# # min(10^(this.subset$modl_acc)),
# # q90.conc =
# # quantile(10^(this.subset$modl_acc),0.9),
# # high.conc =
# # max(10^(this.subset$modl_acc)),
# # stringsAsFactors = F)
# # wambaugh2019.tox21 <- rbind(wambaugh2019.tox21, this.row)
# # }
#
# chem.preds <- httk::wambaugh2019.seem3
# directnhanes.preds <- httk::wambaugh2019.nhanes
# Tox21.acc.numeric <- httk::wambaugh2019.tox21
#
# #Subset to those chemicals with HTS hits:
# human.hits <- subset(ceetox,CAS%in%Tox21.acc.numeric$casn)
#
# # Now calculate the oral equivalent dose for each chemical:
# # Must make sure that CSS is in units of uM!!
# for (this.chem in human.hits$CAS)
# {
# med.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"med.conc"]
# q10.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"q10.conc"]
# low.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"low.conc"]
# q90.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"q90.conc"]
# high.conc <-
# Tox21.acc.numeric[Tox21.acc.numeric$casn==this.chem,"high.conc"]
# css95.v <- human.hits[human.hits$CAS==this.chem][["css95_v"]]
# css95.uv <- human.hits[human.hits$CAS==this.chem][["css95_uv"]]
# human.hits[human.hits$CAS==this.chem,"HTS.Median.ACC"] <- med.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Q10.ACC"] <- q10.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Q90.ACC"] <- q90.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Low.ACC"] <- low.conc
# human.hits[human.hits$CAS==this.chem,"HTS.High.ACC"] <- high.conc
# human.hits[human.hits$CAS==this.chem,"HTS.Median.equivdose.v"] <- med.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Q10.equivdose.v"] <- q10.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Q90.equivdose.v"] <- q90.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Low.equivdose.v"] <- low.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.High.equivdose.v"] <- high.conc/css95.v
# human.hits[human.hits$CAS==this.chem,"HTS.Median.equivdose.uv"] <- med.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.Q10.equivdose.uv"] <- q10.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.Q90.equivdose.uv"] <- q90.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.Low.equivdose.uv"] <- low.conc/css95.uv
# human.hits[human.hits$CAS==this.chem,"HTS.High.equivdose.uv"] <- high.conc/css95.uv
# }
#
#
#
# #Add the exposure predictions:
# for (this.chem in human.hits$CAS)
# {
# print(this.chem)
# # If we have direct inferrences from NHANES, use those exposures:
# if (this.chem %in% directnhanes.preds$CASRN)
# {
# human.hits[human.hits$CAS==this.chem,"Exposure.median"] <- directnhanes.preds[directnhanes.preds$CASRN==this.chem,"lP"]
# human.hits[human.hits$CAS==this.chem,"Exposure.high"] <- directnhanes.preds[directnhanes.preds$CASRN==this.chem,"lP.max"]
# # Otherwise use the heuristics model (Wambaugh 2014):
# } else if (this.chem %in% chem.preds$CAS) {
# human.hits[human.hits$CAS==this.chem,"Exposure.median"] <- chem.preds[chem.preds$CAS==this.chem,"seem3"]
# human.hits[human.hits$CAS==this.chem,"Exposure.high"] <- chem.preds[chem.preds$CAS==this.chem,"seem3.u95"]
# }
# }
#
# # Drop chemicals without exposure predictions:
# human.hits<- subset(human.hits,!is.na(Exposure.median))
#
#
# # This code puts the chemicals into order by margin between exposure and activity:
# chem.names <- unique(human.hits$Compound)
# potency <- rep(Inf,times=length(chem.names))
# names(potency) <- chem.names
# potency.u <- potency
# for (this.chem in chem.names)
# {
# this.subset <- subset(human.hits,Compound==this.chem)
# potency[this.chem] <- this.subset$HTS.Q10.equivdose.v/this.subset$Exposure.high
# potency.u[this.chem] <- this.subset$HTS.Q10.equivdose.uv/this.subset$Exposure.high
# }
# potency <- sort(potency)
# human.hits$Compound <- factor(human.hits$Compound, levels=names(potency))
#
# potency.u <- potency.u[names(potency.u)[potency.u<1]]
# potency.u <- potency.u[names(potency)[potency>1]]
#
# new.overlaps <- names(potency.u[!is.na(potency.u)])
# first.change <- new.overlaps[1]
# last.change <- new.overlaps[length(new.overlaps)]
# window.start <- names(potency)[which(names(potency)==first.change)-5]
# window.end <- names(potency)[which(names(potency)==last.change)+5]
# window.names <- names(potency)[(which(names(potency)==first.change)-5):(which(names(potency)==last.change)+5)]
# # Initialize a new plot:
# Fig7a <- ggplot(data=human.hits) +
# geom_segment(aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.v,yend=HTS.Q90.equivdose.v),size=1,color="white",alpha=0.5) +
# geom_rect(mapping=aes(xmin=window.start,xmax=window.end,ymin=10^-9,ymax=10^3),color="grey",alpha=0.5)+
# geom_segment(aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.v,yend=HTS.Q90.equivdose.v),size=1,color="red",alpha=0.5) +
# geom_point(aes(x=Compound,y=HTS.Median.equivdose.v),shape=3) +
# geom_point(aes(x=Compound,y=HTS.Q10.equivdose.v)) +
# theme_bw() +
# theme(axis.text.x = element_blank()) +
# theme(axis.title.x = element_text(size=16)) +
# theme(axis.title.y = element_text(size=16)) +
# scale_y_log10(label=scientific_10,limits = c(10^-9,10^3)) +
# geom_segment(aes(x=Compound,xend=Compound,y=Exposure.high,yend=Exposure.median),size=1,color="blue",alpha=0.5) +
# geom_point(aes(x=Compound,y=Exposure.median),shape=2) +
# ylab("Bioactive Dose & Exposure\nmg/kg BW/day\n(Variability Only)") +
# xlab("Chemicals Ranked By Ratio Between Bioactive Dose and Exposure")+
# annotate("text", size=8,x = names(potency)[10], y = 100, label = "a")
#
# # Initialize a new plot:
# Fig7c <- ggplot() +
# geom_segment(data=human.hits,aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.v,yend=HTS.Q90.equivdose.v),size=3,color="grey") +
# geom_segment(data=human.hits,aes(x=Compound,xend=Compound,y=HTS.Q10.equivdose.uv,yend=HTS.Q90.equivdose.uv),size=3,color="red",alpha=0.5) +
# geom_point(data=human.hits,aes(x=Compound,y=HTS.Median.equivdose.uv),shape=3,size=2) +
# geom_point(data=human.hits,aes(x=Compound,y=HTS.Q10.equivdose.uv),size=2) +
# theme_bw() +
# theme(axis.text.x = element_text(angle = 60, hjust = 1,size=6)) +
# theme(axis.text.x = element_blank()) +
# theme(axis.title.x = element_text(size=16)) +
# theme(axis.title.y = element_text(size=16)) +
# xlim(window.names) +
# scale_y_log10(label=scientific_10,limits = c(10^-8,2*10^2)) +
# geom_segment(data=human.hits,aes(x=Compound,xend=Compound,y=Exposure.high,yend=Exposure.median),size=3,color="blue",alpha=0.5) +
# geom_segment(data=subset(human.hits,Compound%in%names(potency.u)),aes(x=Compound,xend=Compound,y=10^-8,yend=Exposure.median/1.5),size=2,arrow=arrow(length=unit(0.5,"cm")))+
# geom_point(data=human.hits,aes(x=Compound,y=Exposure.median),shape=2,size=2) +
# xlab("") +
# ylab("Bioactive Dose & Exposure\nmg/kg BW/day\n(Uncertainty and Variability)") +
# annotate("text", size=8,x = window.names[1], y = 100, label = "b")
#
#
# #print(Fig10b)
# #dev.new(width=10,height=6)
# multiplot(Fig7a,Fig7c,cols=1,heights=c(0.5,0.5))
#
# write.csv(subset(human.hits,Exposure.high>HTS.Q10.equivdose.uv)[,c("Compound","DSSTox_Substance_Id","Human.Clint","Human.Funbound.plasma","HTS.Q10.ACC","HTS.Q10.equivdose.v","HTS.Q10.equivdose.uv","Exposure.high")],file=paste("SupTable5-",Sys.Date(),".txt",sep=""),row.names=FALSE)
## ----calc.stats, eval = execute.vignette-------------------------------------
# Wetmore.chems <- subset(chem.physical_and_invitro.data,regexpr("Wetmore",Human.Clint.Reference)!=-1)[,c("Compound","CAS","DTXSID","Formula")]
# Wetmore.chems <- subset(Wetmore.chems,CAS%in%get_cheminfo())
#
# for (this.cas in Wetmore.chems$CAS)
# {
# Wetmore.chems[Wetmore.chems$CAS==this.cas,"Css.Lit"] <- get_lit_css(
# chem.cas=this.cas,
# output.units="uM")
# set.seed(42)
# Wetmore.chems[Wetmore.chems$CAS==this.cas,"Css.v192"] <- calc_mc_css(
# chem.cas=this.cas,
# fup.meas.cv=NULL,
# clint.meas.cv=NULL,
# fup.censored.dist=TRUE,
# output.units="uM")
# set.seed(42)
# Wetmore.chems[Wetmore.chems$CAS==this.cas,"Css.v110"] <- calc_mc_css(
# chem.cas=this.cas,
# output.units="uM")
# }
#
# median(Wetmore.chems$Css.v110/Wetmore.chems$Css.Lit,na.rm=TRUE)
#
# ceetox$Human.Clint2 <- sapply(ceetox$Human.Clint,function(x) as.numeric(strsplit(x,",")[[1]][1]))
# ceetox$Human.Fup2 <- sapply(ceetox$Human.Funbound.plasma,function(x) as.numeric(strsplit(x,",")[[1]][1]))
# write.csv(ceetox[,c(1,2,8,9,16,17,6,7,12,30,10,11,15,31,13,14,18:28)],row.names=FALSE,file=paste("Supplemental-Table4-",Sys.Date(),".txt",sep=""))
#
# print(paste("Complete HTTK data on",sum(get_cheminfo()%in%ceetox$CAS),"new chemicals."))
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.