test.R

In Thibaudpmx/VirtualTumor: What the Package Does (One Line, Title Case)

library(VirtualTumor)
library(peccary)
library(RxODE)

# Step 1 create or load project ---------------------------------------------------
create_VT_Project("D:/these/Second_project/QSP/modeling_work/Lind_eq_VTpckg")



bags <- VT_extract_bags(VT)
saveRDS(bags, "D:/these/Second_project/QSP/modeling_work/Lind_eq_VTpckg/3_virtual_tumors/SUDHL4-Venetoclax_bags.RDS")



bags %>%
  group_by(bag) %>%
  tally %>%
  arrange(desc(n)) %>%
  tail


# Tentative number 1: minimal sum BIM PUMA NOA

bags %>%
  mutate(BIM0s = BIM0 / mean(BIM0),
         PUMA0s = PUMA0 / mean(PUMA0),
         NOXA0s = NOXA0 / mean(NOXA0)) %>%
  mutate(sumpro = BIM0s + PUMA0s + NOXA0s) %>%
  arrange(sumpro) %>%
  group_by(bag) %>%
  slice(1) %>%
  ungroup %>%
  gather("Prot", "value", BAK0, BAXc0, Bcl20, Bclxl0, Mcl10, BIM0, PUMA0, NOXA0) %>%
  ggplot()+
  geom_histogram(aes(value))+
  facet_wrap(~Prot, scales = "free")



# Analyse VT -----------------------------------------------

Vt <- load_VT()

bags <- Vt$extract_bags()


desc <- Vt$description(returnAna = T)

desc %>%
  distinct() %>%
  mutate(nsimilar = map_dbl(id, ~ nrow(bags %>% filter(cellid_in_harvest == .x)))) %>%
  mutate(ninharvest = map_dbl(id, ~ sum(Vt$harvest == .x)))  -> analyse_complete


analyse_complete %>%
  gather("key", "value", drug1Level, drug2Level) %>%
  ggplot()+
  geom_boxplot(aes(as.factor(value), ninharvest))+
  facet_wrap(~key, scales = "free")

analyse_complete %>%
  # gather("key", "value", drug1Level, drug2Level) %>%
  ggplot()+
  geom_point(aes(ninharvest, total))+
  theme_bw()+
  labs(x = "Pct of that cell in the VT", y = "Number of concentrations leading to death")


analyse_complete %>%
  # filter(nsimilar < 10000) %>%
  # gather("key", "value", drug1Level, drug2Level) %>%
  ggplot()+
  geom_point(aes(nsimilar, total))+
  theme_bw()+
  labs(x = "Number of cells with same profil", y = "Number of concentrations leading to death")+
  scale_x_log10()



analyse_complete %>%
  # filter(nsimilar < 10000) %>%
  # gather("key", "value", drug1Level, drug2Level) %>%
  ggplot()+
  geom_point(aes(nsimilar, ninharvest))+
  theme_bw()+
  labs(x = "Number of cells with same profil", y = "Number of concentrations leading to death")+
  scale_x_log10()



# GOF ---------------------------------------------------------------------

Vt$plot(drug = c(1,4), accesRes = T) %>%
  mutate(group = c("Venetoclax")) %>%
  bind_rows(

    Vt$plot(drug = c(2,4), accesRes = T) %>%
      mutate(group = c("A-1155463"))

  ) %>%
  spread(key = Reconst, value = res) %>%
  mutate(Value = if_else(Value > 100, 100, Value)) -> reconst

reconst %>%
ggplot()+
  geom_point(aes(Value, res, col = group))+
  geom_abline(slope = 1, intercept = 0)+
  # scale_y_log10()+
  # scale_x_log10()+
  theme_bw()+
  labs(x = "Observation (pct viability)", y = "Reconstitution (pct viability)", title = "Pred vs Obs") +
  theme(plot.title = element_text(hjust = 0.5)) -> gof1;gof1

reconst %>%
  mutate(residuals = res-Value) %>%
  mutate(conc = if_else(conc1 == 0, conc2, conc1)) %>%
  ggplot()+
  geom_point(aes(conc, residuals, col = group))+
  scale_x_log10()+
  theme_bw()+
  geom_hline(yintercept = 0)+
  labs(x = "Concentration (µM)", y = "Residuals (Prev - obs)", title = "Residuals") +
  theme(plot.title = element_text(hjust = 0.5)) -> gof2;gof2

  plot_grid(gof1, gof2)


# Analysis 3 drugs --------------------------------------------------------

desc <- Vt$description()
  desc
self <- Vt

  celltheque <- self$celltheque

  harvest <- self$harvest


  ## now let's compute the analysis


  tibble(id = unique(harvest))  %>%
    mutate(analyse = map_chr(id, function(id){
      # print(id)
      #empty tibble to fill for futur unnest
      # results <- tibble(a = NA)

    veneto <-  matrix_cell(celltheque = celltheque %>% filter(group == "Drug1_4"), id, drug = c(1,4), plot = F)
    a12 <-  matrix_cell(celltheque = celltheque %>% filter(group == "Drug1_4"), id, drug = c(4,1), plot = F)
    a11 <-  matrix_cell(celltheque = celltheque %>% filter(group == "Drug2_4"), id, drug = c(2,4), plot = F)

    veneto_single <- if_else(sum(veneto$`0`) > 0, T , F )
    a12_single <- if_else(sum(a12$`0`) > 0, T , F )
    a11_single <- if_else(sum(a11$`0`) > 0, T , F )

    veneto_a12 <- if_else(sum(as.matrix(veneto[ , -1])) > 0 , T, F)
    a11_a12 <-  if_else(sum(as.matrix(a11[ , -1])) > 0 , T, F)

    case_when(
      veneto_single == T & a12_single == T & a11_single == T ~ "any_of_three",
      veneto_single == T & a12_single == F & a11_single == F ~ "Veneto-only",
      veneto_single == F & a12_single == T & a11_single == F ~ "a12-only",
      veneto_single == F & a12_single == F & a11_single == T ~ "a11-only",
      veneto_single == T & a12_single == T & a11_single == F ~ "Veneto_or_a12",
      veneto_single == T & a12_single == F & a11_single == T ~ "Veneto_or_a11",
      veneto_single == F & a12_single == T & a11_single == T ~ "a12_or_a11",
      veneto_a12 == T & a11_a12 == F ~ "veneto_and_a12",
      veneto_a12 == F & a11_a12 == T ~ "a11_and_a12",
      veneto_a12 == T & a11_a12 == T ~ "any_combination",
      T ~ "todetermine"
      )

    })) -> toana

tibble(id = Vt$harvest) %>%
  left_join(toana) %>%
  group_by(analyse) %>%
  tally %>%
  arrange(desc(n))


# Recompute protein bags after equilibrium --------------------------------


Vt <- load_VT()
celltheque <- Vt$celltheque
harvest <- Vt$harvest

matrix_cell(celltheque = celltheque, cellID = 833, drug = c(1,4))

bags <- Vt$extract_bags()

bags %>%
  group_by(bag) %>%
  tally

bags1 <- bags %>%
  filter(bag == 1)


bags1


compute_at_equilibrium <- function( df = bags1){


  df %>%
    select(BAK0, BAXc0, Bcl20, Mcl10, BIM0, PUMA0, NOXA0,Bclxl0 ) %>%
    rowid_to_column() %>%
    mutate(new = map(rowid, function(rowid){

    line <- df %>%
      slice(rowid)


    add_events <- tibble(cmt = c("Bclxl_I",  "Mcl1_I", "Bcl2_I")) %>%  # Administration sampling, with "concX" being replaced by concentration
      mutate(time = 700, amt = c(0, 0, 0))

    res <- simulations(ind_param = line, add_events = add_events, returnSim = T)

    res %>%
      select(time, Bcl2, BIM, PUMA, NOXA, BAK, BAXc, Bclxl, Mcl1) %>%
      filter(time == 700) %>%
      select(-time)-> temp

    names(temp) <- paste0(names(temp) , "_eq")

    temp
    })) -> test


  test %>%
    group_by(BAK0, BAXc0) %>%
    tally

  test %>%
    filter(BAK0 == 1000, BAXc0 == 0) %>%
    unnest() %>%
    gather("prot", "value", Bcl2_eq, BIM_eq, PUMA_eq, NOXA_eq, BAK_eq, BAXc_eq, Bclxl_eq, Mcl1_eq) %>%
    ggplot()+
    geom_histogram(aes(value))+
    facet_wrap(~prot, scales = "free")+
    ggtitle("At equilibrium")-> eq

  test %>%
    filter(BAK0 == 1000, BAXc0 == 0) %>%
    unnest() %>%
    gather("prot", "value", Bcl20, BIM0, PUMA0, NOXA0, BAK0, BAXc0, Bclxl0, Mcl10) %>%
    ggplot()+
    geom_histogram(aes(value))+
    facet_wrap(~prot, scales = "free")+
    ggtitle("At t0")-> noteq

  plot_grid(noteq ,eq )

  test %>%
    unnest() %>%
    filter(Mcl1_eq  > Mcl10 * 0.8)
}



# Find missing  profiles --------------------------------------------------

library(VirtualTumor)
library(peccary)
library(RxODE)

# Step 1 create or load project ---------------------------------------------------
create_VT_Project("D:/these/Second_project/QSP/modeling_work/Lind_eq_VTpckg")

# All we have right now

cellt <- cellthequeLoad(drug = c(1,4), return = 1:2)

allth <- celltheque_theo_full()
# saveRDS(cellt, "D:/these/Second_project/QSP/modeling_work/Lind_eq_VTpckg/temp.RDS")

allth %>%
  group_by(cellid) %>%
  slice(1)
# get all the theorical with the number of first Trye
allth %>%
  mutate(group = paste0("conc1=", conc1, "conc4=", conc2)) %>%
  select(cellid, group, res, A, B,C, D) %>%
  spread(group, res) %>%
  select(-cellid) -> alltheospread

# missing <-
# alltheospread %>%
#   left_join( cellt[[2]] %>% mutate(test = T)) %>%
#   filter(is.na(test)) %>%
#   select(A, B, C, D) %>%
#   mutate(sumSurvival = A + B + C+D)
#
#
# cellt

#
# Try to have the cells in celltheque with same number system
cellt[[1]] %>%
  # filter(cellid == 17 & conc4 == 0)
  group_by(cellid, conc4) %>%
  filter(res == T) %>%
  slice(1) %>%
  arrange(cellid) %>%
  select(cellid, conc1, conc4) %>%
  rename(drug1 = conc1) %>%
  mutate(n = map_dbl(drug1, function(x){

    which(conc1 == x) - 1

  })) %>%
  select(-drug1) %>%
spread(conc4, n) %>%
  map_df(function(x){


    x[is.na(x)] <- 10
    x
  }) -> celltfornow

names(celltfornow)  <- c("cellid", "A2", "B2", "C2", "D2")

matrix_cell(celltheque =cellt[[1]], cellID = 241, drug = c(1,4) )

# missing ones
missing <- alltheospread %>%
  left_join( cellt[[2]] %>% mutate(test = T)) %>%
  filter(is.na(test)) #%>%
  # select(A, B, C, D)

crossing(celltfornow, missing %>% select(A, B, C, D)) %>%
  mutate(difA = A2 - A, difB = B2 - B, difC = C2 - C, difD = D2 - D ) %>%
  mutate(ngroupdif = pmap_dbl(list(difA, difB, difC, difD), function(difA, difB, difC, difD){

sum(c(difA, difB, difC, difD) != 0 )


  })) %>%
  filter(ngroupdif == 1) %>%
  mutate(nmax = pmap_dbl(list(difA, difB, difC, difD), function(difA, difB, difC, difD){

    max(abs(c(difA, difB, difC, difD )))


  })) %>%
  arrange(nmax) %>%
  filter(nmax == 1) -> un_pret#%>%
  # group_by(A, B, C, D ) %>%
  # tally %>% arrange(desc(n))

un_pret %>%
group_by(A, B, C, D ) %>%
tally %>% arrange(desc(n))

# les cellids à faire en priorité !!!
un_pret %>%
  group_by(cellid) %>%
  tally %>% arrange(desc(n))

un_pret %>% filter(cellid == 61819)

matrix_cell(cellt[[1]], cellID = 61819, drug = c(1,4))

cellID
line <- cellt[[2]] %>% filter(cellid == 61819)

toadd <- crossing(Bcl20 = line$Bcl20 * seq(0.8,1.2,0.1),
                  Mcl10 = line$Mcl10 * seq(0.8,1.2,0.1),
                  Bclxl0 = line$Bclxl0 * seq(0.8,1.2,0.1),
                  BIM0 = line$BIM0 * seq(0.8,1.2,0.1),
                  PUMA0 = line$PUMA0 * seq(0.8,1.2,0.1),
                  NOXA0 = line$NOXA0 * seq(0.8,1.2,0.1),
                  BAXc0 = line$BAXc0,
                  BAK0 = line$BAK0)



add_events <- tibble(cmt = c("Bclxl_I",  "Mcl1_I", "Bcl2_I")) %>%  # Administration sampling, with "concX" being replaced by concentration
  mutate(time = 700, amt = c("0", "conc4", "conc1"))


celltheque_produc(file.name = "temp.RDS", drug = list(c(1,4)), toadd = toadd, add_events = add_events, update_at_end = F, saven = 200,time_compteur = F)

celltheque_verification(test, add_events = add_events, drug = c(1,4))

test <- readRDS("D:/these/Second_project/QSP/modeling_work/Lind_eq_VTpckg/2_celltheques/celltheques/temp.RDS")

candidates <- celltheque_spread(drug = c(1,4), celltheque = test, returnOnIDperGroup = T)


alltheospread %>%
  left_join( cellt[[2]] %>% mutate(test = T)) -> allt

allt[ , grepl("(conc)|(test)", names(allt))] %>%
  select(test, everything()) %>%
  mutate(test=if_else(is.na(test), "missing", "already have")) -> allt2

candidates %>%
  left_join(

    allt2
  ) %>%
  filter(test == "missing") %>%
  pull(cellid) %>%
  map(~ matrix_cell(celltheque = test, cellID = .x, drug = c(1,4))) %>%
  invoke(.fn = plot_grid)





newlines <- readRDS("D:/these/Second_project/QSP/modeling_work/Lind_eq_VTpckg/newprofiles.RDS") %>%
  reduce(bind_rows)


newlines %>%
  select(-n, - cellid, -Bcl20, - Mcl10, - Bclxl0, - BIM0, -PUMA0, - NOXA0, -BAXc0, - BAK0,-test) %>%
  distinct()