In sgosline/mpnstXenoModeling: MPNST Xenograft Modeling

knitr::opts_chunk$set(echo = TRUE)
require(remotes)

if(!require(mpnstXenoModeling)){
  remotes::install_github('sgosline/mpnstXenoModeling')
  library(mpnstXenoModeling)
}
#remotes::install_github('https://github.com/chapmandu2/IncucyteDRC')

if(!require('leapr')){
  remotes::install_github('biodataganache/leapR')
  library(leapr)
}
if(!require('Xeva')){
  BiocManager::install('Xeva')
  library(Xeva)
}
if(!require('foreach')){
    install.packages('foreach')
    library(foreach)
}
library(data.table)
library(dplyr)
library(repr)
library(pheatmap)
library(ggplot2)
library(drc)
library(pROC)
library(tidyr)

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.

mpnstXenoModeling::loadPDXData()

icyteData<-tidyr::unite(icyteData, "compound_name", c(experimentalCondition,dosage,synid,replicate),na.rm=TRUE,remove=FALSE)

treatTab <- icyteData[icyteData$experimentalCondition != 'DMSO',]
contTab <- icyteData[icyteData$experimentalCondition == 'DMSO',]

### function taken from xenoResponseMetrics, calculating AUC with Xeva package
computeAUC<-function(treatedTab,contTab){
  #https://link.springer.com/article/10.1208/s12248-018-0284-8
  tauc=treatedTab%>%mutate(volume=as.numeric(response),time=experimental_time_point)%>%
        group_by(compound_name)%>%
        group_map(~ unlist(Xeva::AUC(.x$time,.x$volume))[['value']],.keep=TRUE)
  tauc.list <- as.numeric(tauc)
  names(tauc.list) <- unique(treatedTab$compound_name)
  cauc=contTab%>%mutate(volume=as.numeric(response),time=experimental_time_point)%>%
    group_by(compound_name)%>%
    dplyr::select(compound_name,time,volume)%>%
    group_map(~ unlist(Xeva::AUC(.x$time,.x$volume))[['value']],.keep=TRUE)
  caucval <- mean(as.numeric(unlist(cauc)),na.rm=T)
  aucfunc = function(x){1-((caucval-x)/caucval)}
  ret <- lapply(tauc.list,aucfunc)
  return(ret)
}

auc.all <- computeAUC(treatTab,contTab)

filter_byDrug <- function(icd, drugID) {
  icd <- subset(icd, compound_name == drugID)
  icd <- icd[order(icd$experimental_time_point),]
  #icd$dosage <- log((icd$dosage)/1000000000)
  # Assumes log(M) concentration
  return(icd)
}

TryFit <- function(time_response, fixed = c(NA, NA, NA, NA), names = c(NA, NA, NA, NA), nan.handle = c("LL4", "L4"), curve_type){
    if (var(time_response$response) == 0) {
        time_response$response[nrow(time_response)] <- time_response$response[nrow(time_response)] + 10^-10
    }
    #dose_response[dose_response == 0] <- 10^-10
    nan.handle <- match.arg(nan.handle)
    #, logDose = exp(10)
    drug.model <- tryCatch({
        drcmod(time_response, LL.4(fixed = fixed, names = names), curve_type)
    }, warning = function(w) {
        if(nan.handle == "L4"){
            drcmod(time_response, L.4(fixed = fixed, names = names), curve_type)
    } else {
        drcmod(time_response, LL.4(fixed = fixed, names = names), curve_type)
    }
    }, error = function(e) {
        drcmod(time_response, L.4(fixed = fixed, names = names), curve_type)
    })
    return(drug.model = drug.model)
}
drcmod <- function(time_resp, fctval, curve_type){
  temp <- drm(formula   = response ~ experimental_time_point
      , curveid   = model_system_name
      , data      = time_resp
      , fct       = fctval
      , na.action = na.omit
      , control   = drmc(errorm = FALSE)
  )
  summary(temp)
  return(temp)
}
generate_TR_plots <- function(res, drugID) {
  require(data.table)
  tr.df <- res %>% dplyr::select(experimental_time_point, response, model_system_name, dosage) #%>% setDT()
  tr.df <- tidyr::unnest(tr.df, response)
  tr.df$response<- na_if(tr.df$response, "NA")
  tr.df <- tr.df[complete.cases(tr.df), ]
  tr.dt <- data.table(tr.df)
  dt2 <- data.table()
  dat.dt <- data.table(compound_name=character(),model_system_name=character(),dosage=numeric(),Hill=numeric(),ec50=numeric(),
                       MinViability=numeric(),MaxViability=numeric(),ic50=numeric())
  # factor by CellLine
  tr.dt[,model_system_name:=as.factor(model_system_name)]
  scale.num <- nlevels(tr.dt$model_system_name)
  # iterate over CellLine: Conc, Viabilities
  for (i in 1:nlevels(tr.dt$model_system_name)) {
    # subset data.table by CellLine level
    dt <- tr.dt[model_system_name %in% levels(model_system_name)[i]]
    dt[, model_system_name := as.character(model_system_name)]
    # TryFit of subsetted data.table
    #dft <- dt[,.(Viabilities = unlist(Viabilities)), by = setdiff(names(dt), 'Viabilities')]
    df <- as.data.frame(dt)
    min_value <- min(df$response)
    max_value <- max(df$response)
    fit.LL4 <- TryFit(df, fixed = c(NA, min_value, max_value, NA), names = c("hill", "min_value", "max_value", "ec_50"), nan.handle = "L4", "model_system_name")
    # 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(response)/4),by=experimental_time_point]
    df <- as.data.frame(dt)
    cells <- unique(df$model_system_name)
    dose <- unique(df$dosage)
    # pROC::auc() https://cran.r-project.org/web/packages/pROC/pROC.pdf
    #auc <- auc(df$response,df$experimental_time_point)
    meta2 <- data.table(compound_name=drugID,model_system_name=cells,dosage=dose,Hill=coefs['hill'],ec50=coefs['ec_50'],
                       MinViability=min_value, MaxViability=max_value,ic50=ic_50)
    dat.dt <- rbind(dat.dt,meta2)
    dt = merge(dt,sd.lst)
    dt2 <- rbind(dt2, dt)
  }
  return(list('plot'=dt2,'table'=dat.dt))
}

tr_plot<-function(dt2) {
    drug=unique(dt2$drug)
    dose=unique(dt2$dosage)
    plot=ggplot(dt2, aes(x=experimental_time_point, y=response, group=model_system_name, color=model_system_name)) + #, shape=model_system_name
            geom_ribbon(aes(y=Pred, ymin=Pred-SD, ymax=Pred+SD, fill=model_system_name), alpha=0.2, color=NA) +
            geom_line(aes(y = Pred)) +
            scale_shape_manual(values=seq(0, nlevels(dt2$model_system_name))) +
            labs(x = "Time", y = "Confluency %", shape="Temp", color="Temp") +
            theme_bw() +
            ggtitle(paste("Time-response curves for", drug, "at dose", dose, "nM"))
    return(plot)
    }

process_IncucyteDrugData <- function(treatTab, upload=FALSE, path='.', parentID=NULL){
  test=treatTab%>%dplyr::select(compound_name,dosage, synid)%>%distinct()
  drugIDs <- test$compound_name
  synIDs <- unique(test$synid)
  meta.dt <- data.table(compound_name=character(),model_system_name=character(),dosage=numeric(),Hill=numeric(),ec50=numeric(),MinViability=numeric(),MaxViability=numeric(),ic50=numeric())
  plt.list <- list()
  pi=1
  for (drug in drugIDs) {
    dres <- subset(treatTab, compound_name == drug)
    dres <- split(dres,dres$dosage)
    for (res in dres){
      dose <- unique(res$dosage)
      plt = generate_TR_plots(res,drug)
      meta.dt <- rbind(meta.dt, plt$table)
      plt.list[[pi]] <- cbind(plt$plot,drug)
      pi = pi+1

    }
  }
  meta.dt<-meta.dt%>%tidyr::replace_na(list(ic50=0.0))%>%as.data.frame()
  #fwrite(meta.dt,file="incucyteDoseResponse_table.csv")
  pdf("incucyteDoseResponsePlots.pdf")#, onefile = TRUE
  master.plot<-lapply(plt.list, tr_plot)
  #invisible(lapply(master.plot, print))
  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)
}

meta.dt<-process_IncucyteDrugData(treatTab,contTab, upload=FALSE)

auc.df <- as.data.frame(t(as.data.frame(do.call(cbind, auc.all))))%>%tibble::rownames_to_column()
colnames(auc.df) <- c("compound_name", "auc")
meta.dt <- merge(meta.dt,auc.df,by='compound_name')
out.dt <- merge(meta.dt,icyteData%>%
               dplyr::select(model_system_name, experimentalCondition, compound_name, dosage)%>%
               distinct(),
               by=c('compound_name','model_system_name', 'dosage'))

#fwrite(out.dt,file="incucyteDoseResponse_table.csv")

res  = syn$tableQuery('select * from syn26608026')
syn$delete(res)
#pd <- sync$build_table('Incucyte Response Value','syn21984813','incucyteDoseResponse_table.csv')
pd <- sync$build_table('Incucyte Response Value','syn21984813',out.dt)
syn$store(pd)

Analysis of IncuCyte % confluence

icStats0 <- treatTab%>%
            mutate(maxKilling = (response/100))
minDrugs <- icStats0%>%
            mutate(hasRNASeq=(model_system_name%in%rnaSeq$Sample))%>%
            group_by(compound_name)%>%
            summarise(samples=n_distinct(model_system_name),samplesWithRNA=count(hasRNASeq==TRUE))%>%
            subset(samples>2)
#DT::datatable(minDrugs)
icStats0<-tidyr::unite(icStats0,'Drug',c(experimentalCondition,dosage))

spi.mat<-icStats0%>%
  dplyr::select(Drug,model_system_name,maxKilling)%>%
  group_by(Drug,model_system_name)%>%
  summarise(maxKilling=max(maxKilling))%>%
  tidyr::pivot_wider(values_from=maxKilling,names_from=model_system_name)%>%#,values_fill=list(maxKilling=0))%>%
  tibble::column_to_rownames('Drug')%>%
  as.matrix()
#pheatmap(spi.mat,cellwidth = 10,cellheight=10,main='max killing across Incucyte data')
pheatmap(spi.mat,cellwidth = 10,cellheight=10,main='max killing across Incucyte data',filename='MK_ic.png')

##let's plot the stats for each drug
icStats<-merge(meta.dt,icyteData%>%
               dplyr::select(model_system_name, experimentalCondition, compound_name, dosage)%>%
               distinct(),
               by=c('compound_name','model_system_name', 'dosage'))
icStats<-tidyr::unite(icStats,'Drug',c(experimentalCondition,dosage))

auc.mat<-icStats%>%
  dplyr::select(Drug,model_system_name,auc)%>%
  aggregate(auc~Drug+model_system_name,data=.,mean)%>%
  tidyr::pivot_wider(values_from=auc,names_from=model_system_name)%>%#,values_fill=list(auc=0))#%>%
  tibble::column_to_rownames('Drug')%>%
  as.matrix()

pheatmap(auc.mat, cellwidth = 10,cellheight=10,main='AUC across Incucyte data')

tgi.mat<-icStats%>%
  dplyr::select(Drug,model_system_name,ic50)%>%
  aggregate(ic50~Drug+model_system_name,data=.,mean)%>%
  tidyr::pivot_wider(values_from=ic50,names_from=model_system_name)%>%#,values_fill=list(ic50=0))#%>%
  tibble::column_to_rownames('Drug')%>%
  as.matrix()
skip <- which(rowMeans(tgi.mat)%in%c(1,0))
if(length(skip)>0)
  tgi.mat <-tgi.mat[-skip,]
tgi.mat
pheatmap(tgi.mat,cellwidth = 10,cellheight=10,main='IC50 across Incucyte data',filename='ic50_ic.png')
#pheatmap(tgi.mat,cellwidth = 10,cellheight=10,main='IC50 across Incucyte data')