In sgosline/mpnstXenoModeling: MPNST Xenograft Modeling

knitr::opts_chunk$set(echo = TRUE)
if(!require(mpnstXenoModeling)){
  remotes::install_github('sgosline/mpnstXenoModeling')
  library(mpnstXenoModeling)
}
if(!require('dplyr')){
  BiocManager::install('dplyr')
  library(dplyr)
}
if(!require('tidyr')){
  BiocManager::install('tidyr')
  library(tidyr)
}
if(!require('drc')){
  BiocManager::install('drc')
  library(drc)
}
if(!require('ggplot2')){
  BiocManager::install('ggplot2')
  library(ggplot2)
}
if(!require('data.table')){
  BiocManager::install('data.table')
  library(data.table)
}
if(!require('pROC')){
  BiocManager::install('pROC')
  library(pROC)
}

Load all data from Synapse

Let's load all the data using the mpnstXenoModeling package function. This will pull the data from synapse and format it accordingly.

#this function simply loads all the data into memory
loadPDXData()
##let's clean up drug Dat

Load microtissue data and compare

The first thing we want to do is to load the microtissue data and carry out the bayesian correction. This is pretty straightforward.

filter_byDrug <- function(syn, mtd, drugID) {
  mtd <- subset(mtd, experimentalCondition == drugID)
  res = mpnstXenoModeling::getMicroTissueDrugData(syn,mtd)
  out = res[order(res$Conc),]
  # Assumes log(M) concentration
  return(out)
}


TryFit <- function(dose_response, fixed = c(NA, NA, NA, NA), names = c(NA, NA, NA, NA), nan.handle = c("LL4", "L4"), curve_type){

    if (var(dose_response$Viabilities) == 0) {
        dose_response$Viabilties[nrow(dose_response)] <- dose_response$Viabilities[nrow(dose_response)] + 10^-10
    }
    #dose_response[dose_response == 0] <- 10^-10
    nan.handle <- match.arg(nan.handle)
    #, logDose = exp(10)
    drug.model <- tryCatch({
        drcmod(dose_response, LL.4(fixed = fixed, names = names), curve_type)
    }, warning = function(w) {
        if(nan.handle == "L4"){
            drcmod(dose_response, L.4(fixed = fixed, names = names), curve_type)
    } else {
        drcmod(dose_response, LL.4(fixed = fixed, names = names), curve_type)
    }
    }, error = function(e) {
        drcmod(dose_response, L.4(fixed = fixed, names = names), curve_type)
    })
    return(drug.model = drug.model)
}


drcmod <- function(dose_resp, fctval, curve_type){
  drm(formula   = Viabilities ~ Conc
      , curveid   = CellLine
      , data      = dose_resp
      , fct       = fctval
      , na.action = na.omit
      , control   = drmc(errorm = FALSE)
  )
}


generate_DR_plots <- function(res, drugID) {
  require(data.table)
  dr.df <- res %>% dplyr::select(Conc, Viabilities, CellLine) #%>% setDT()
 # dr.df <- tidyr::unnest(dr.df, Viabilities)
  dr.df$Viabilities<- na_if(dr.df$Viabilities, "NA")
  dr.df <- dr.df[complete.cases(dr.df), ]
  dr.dt <- data.table(dr.df)
  dt2 <- data.table()
  dat.dt <- data.table(Drug=character(),CellLine=character(),Hill=numeric(),ec50=numeric(),
                       MinViability=numeric(),MaxViability=numeric(),ic50=numeric(),auc=numeric())
  # factor by CellLine
  dr.dt[,CellLine:=as.factor(CellLine)]
  scale.num <- nlevels(dr.dt$CellLine)
  # iterate over CellLine: Conc, Viabilities
  for (i in 1:nlevels(dr.dt$CellLine)) {
    # subset data.table by CellLine level
    dt <- dr.dt[CellLine %in% levels(CellLine)[i]]
    dt[, CellLine := as.character(CellLine)]
    # TryFit of subsetted data.table
    #dft <- dt[,.(Viabilities = unlist(Viabilities)), by = setdiff(names(dt), 'Viabilities')]
    df <- as.data.frame(dt)
    min_value <- min(df$Viabilities)
    max_value <- max(df$Viabilities)
    fit.LL4 <- TryFit(df, fixed = c(NA, min_value, max_value, NA), names = c("hill", "min_value", "max_value", "ec_50"), nan.handle = "L4", "CellLine")
    # hill from fit.LL4
    # ic50 code from https://rstudio-pubs-static.s3.amazonaws.com/378543_5b5bda32bf0541a485ccc49efbea680a.html
    coefs <- setNames(
      c(fit.LL4$coefficients, min_value, max_value),
      c("hill", "ec_50", "min_value","max_value"
      ))
    ic_50 <- with(as.list(coefs),
      exp(
        log(ec_50) + (1 / hill) * log(max_value / (max_value - 2 * min_value)) #log(ec_50)
      )
    )
    #cbind(auc,ED(fit.LL4,50))
    #rbindlist(dt.dat,auc)
    dt[, Pred := predict(object=fit.LL4)]
    sd.lst = dt[, list(SD=sd(Viabilities)/4),by=Conc]
    df <- as.data.frame(dt)
    cells <- unique(df$CellLine)
    # pROC::auc() https://cran.r-project.org/web/packages/pROC/pROC.pdf
    print(head(df$Conc))
    print(head(df$Viabilities))
    auc <- auc(df$Viabilities,df$Conc)
    meta2 <- data.table(Drug=drugID,CellLine=cells,Hill=coefs['hill'],ec50=coefs['ec_50'],
                       MinViability=min_value, MaxViability=max_value,ic50=ic_50,auc=as.numeric(auc))
    dat.dt <- rbind(dat.dt,meta2)
    dt = merge(dt,sd.lst)
    dt2 <- rbind(dt2, dt)
  }
  temp.plot = ggplot(dt2, aes(x=Conc, y=Viabilities, group=CellLine, color=CellLine, shape=CellLine)) +
            geom_ribbon(aes(y=Pred, ymin=Pred-SD, ymax=Pred+SD, fill=CellLine), alpha=0.2, color=NA) +
            geom_line(aes(y = Pred)) +
            scale_shape_manual(values=seq(0, scale.num)) +
            labs(x = "Conc log(M)", y = "Viabilities %", shape="Temp", color="Temp") +
            theme_bw() +
            ggtitle(paste("Dose-response curves for Drug:", drugID))
  return(list('plot'=temp.plot,'table'=dat.dt))
}


process_MicroTissueDrugData <- function(syn, master.table, upload=FALSE, path='.', parentID=NULL){
  drugIDs <- unique(master.table$experimentalCondition)
  meta.dt <- data.table(Drug=character(),CellLine=character(),Hill=numeric(),ec50=numeric(),MinViability=numeric(),MaxViability=numeric(),ic50=numeric(),auc=numeric())
  plt.list <- list()
  for (drug in drugIDs) {
    res <- filter_byDrug(syn, master.table, drug)%>%unnest(cols=c(`total cell count`, `live cell count`, Viabilities))##added unnest here
    plt = generate_DR_plots(res,drug)
    myplt <- plt$plot
    mytab <- plt$table
    meta.dt <- rbind(meta.dt, plt$table)
    plt.list[[drug]] = plt$plot
  }
  meta.dt<-meta.dt%>%tidyr::replace_na(list(ic50=0.0))%>%as.data.frame()

  fwrite(meta.dt,file="doseResponse_table.csv")
  pdf("doseResponsePlots.pdf", onefile = TRUE)
  for (drug in drugIDs) {print(plt.list[[drug]])}
  dev.off()
  if (isTRUE(upload)) {
    synapseStore(file.path(path, "doseResponsePlots.pdf"),parentId=parentID)
    synapseStore(file.path(path, 'doseResponse_table.csv'),parentId=parentID)
     #ile.path(path, "doseResponse_table.csv"),parentId=parentID)

  }
  return(meta.dt)
}


norm_MTData<-function(syn,mt.meta){
   drugIDs <- unique(mt.meta$experimentalCondition)
  meta.dt <- data.table(Drug=character(),CellLine=character(),Hill=numeric(),ec50=numeric(),MinViability=numeric(),MaxViability=numeric(),ic50=numeric(),auc=numeric())
  plt.list <- list()
  drug.tab<-do.call(rbind,lapply(drugIDs,function(x) unnest(filter_byDrug(syn,mt.meta,x),
                                                            cols = c(`total cell count`, `live cell count`, Viabilities))))

   ##now update viabilities based on altered rate of cell count
  drug.tab%>%
    mutate(prob_cell_count=ecdf(`total cell count`))%>%
    mutate(prob_live_count=ecdf(`live cell count`))

  prob_cell_count=c()
  prob_live_count=c()
  prob_cell_given_live=c()

  ##we want to estimate prob_live_give_cell= prob_cell_count_given_live x prob live / prob cell count

#  for (drug in drugIDs) {
#    res <- filter_byDrug(syn,mt.meta, drug)
#    plt = generate_DR_plots(res,drug)
#    myplt <- plt$plot
#    mytab <- plt$table
#    meta.dt <- rbind(meta.dt, plt$table)
#    plt.list[[drug]] = plt$plot
#  }
# meta.dt<-meta.dt%>%tidyr::replace_na(list(ic50=0.0))%>%as.data.frame()

}

Reprocess MT data

library(reticulate)
library(dplyr)
sync<-reticulate::import('synapseclient')
syn=sync$login()

meta.dt<-process_MicroTissueDrugData(syn, mt.meta, upload=TRUE, parentID='syn26066928')

new.dt <- norm_MTData(syn,met.meta)
res  = syn$tableQuery('select * from syn26136282')
syn$delete(res)
#res  = syn$tableQuery('select * from syn26136282')
#synTableStore(tabname='Miccrotissue Response Valsue',tab=meta.dt,parentId='syn21984813')

pd <- sync$build_table('Normalized Microtissue Response Values','syn21984813','doseResponse_table.csv')
syn$store(pd)