## ----knitrPRrep, include=FALSE, eval=TRUE-------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = '#>', fig.width=5, fig.height=3)

## ----setup--------------------------------------------------------------------
# We love to give warning messages whenever assumptions are used by HTTK,
# but they will overwhelm the output of this vignette so we turn them
# off:
options(warn = -1)
library(data.table)
library(httk)
library(ggplot2)
library(gdata)
library(stringr)
library(scales)
library(grid)
library(reshape)

## ----lmr2---------------------------------------------------------------------
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---------------------------------------------------------------
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------------------------------------------------------
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-----------------------------------------------------
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-----------------------------------------------------------------
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=T,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------------------------------------------------------------------
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=T,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=F]),probs=probs))
}

## ----Wambaugh2019.Fig1--------------------------------------------------------
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=T),2)))
print(paste("Median Titration:",signif(median(Fig1.table$Affinity.Fup.Med,na.rm=T),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=T),"Zero Measurments"))
print(paste("Median Non-Zero:",signif(median(Fig1.table$CLint.1uM.Median,na.rm=T),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=T),2),"CV:",signif(median(Fig1.table$CV1,na.rm=T),1),"/ CI:",signif(median(Fig1.table$Error2,na.rm=T),1),"CV:",signif(median(Fig1.table$CV2,na.rm=T),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=T),2),"CV:",signif(median(Fig1.table$CV1,na.rm=T),2)))    
print(paste("Median Titr. CI:",signif(median(Fig1.table$Error2,na.rm=T),2),"CV:",signif(median(Fig1.table$CV2,na.rm=T),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=T),2),"CV:",signif(median(Fig1.table$CVClint,na.rm=T),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=T),2),"CV:",signif(median(Fig1.table$CVClint,na.rm=T),2)))

multiplot(Fig1a,Fig1c,Fig1b,Fig1d,cols=2,widths=c(1.75,1.75))

## ----Wambaugh2019.Fig2--------------------------------------------------------
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--------------------------------------------------------
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--------------------------------------------------------
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=T),"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=T),"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=T),"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=T),"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=T)

#dev.new()
print(Fig4)  

summary(clearance.fit)

## ----Wambaugh2019.Fig5, eval=F------------------------------------------------
#  # 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=T)
#  
#  
#  
#  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:=T]
#  ceetox[is.na(hlpoint),hlpointEstimated:=F]
#  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=T)
#  
#  
#  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=T)
#  
#  
#  #dev.new(width=8,height=5)
#  multiplot(Fig5a,Fig5b,cols=2,widths=c(1.75,1.75))

## ----Wambaugh2019.Fig6, eval=F------------------------------------------------
#  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=F------------------------------------------------
#  # 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=F)
#  # # 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 margine 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=F)

## ----calc.stats, eval=F-------------------------------------------------------
#  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=T,
#      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=T)
#  
#  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=F,file=paste("Supplemental-Table4-",Sys.Date(),".txt",sep=""))
#  
#  print(paste("Complete HTTK data on",sum(get_cheminfo()%in%ceetox$CAS),"new chemicals."))

