We tested whether we could exploit the cell-derived Boolean models to identify novel personalized combinatorial treatments. To these aims, each patient’s mutational profile (Fig. S6A and D, Table S5) was turned into the initial activity of a personalized Boolean model (Fig. 5A, panel c).
Next, for each patient, we performed a simulation of the following conditions in silico:
Importantly, our approach enabled us to obtain patient-specific predictive Boolean models able to describe the drug-induced signaling rewiring (Fig. 5F and Fig. S7) and to quantify apoptosis inhibition and proliferation levels (Fig. 5A, panel d and e, Fig. 5E and Fig. S6E).
# Read the mutation file of the patients
mutations <- read_tsv('./input/FLT3_ITD_patients_mutational_profile.tsv')
mutations[mutations == 0] <- 3 #Set NA mutation to a non Boolean value
mutations[mutations == 1] <- 0 #Set LOF mutation to Boolean value 0
mutations[mutations == 2] <- 1 #Set GOF mutation to Boolean value 1
# Create an empty list for each patient' plots
patients_plot_list <- list()
# Counter for plot list
plot_idx = 1
# For each patient
for(p_idx in 1:nrow(mutations)){
# Patient name and mutation vector
patient <- mutations$Patient[p_idx]
mutations_vector <- unlist(mutations[p_idx,2:6])
# Get the two cell-derived models subfamily
JMD_session <- load(paste0('./results/JMD_subfamily_threshold.RData'))
JMD_best <- as.vector(unlist(best_100_models[1,]))
JMD_model <- modelT0
JMDsimList <- simList
TKD_session <- load(paste0('./results/TKD_subfamily_threshold.RData'))
TKD_best <- as.vector(unlist(best_100_models[1,]))
TKD_model <- modelT0
TKDsimList <- simList
# set the conditions we want to simulate
# 1) patients' mutational profile + active TNF IGF1 receptors
# 2) condition 1 + FLT3i
# 3) condition1 + FLT3i + druggable node
# To create the collection of nodes to inhibit
# we start from the node inhibited in the perturbation experiment
inhibited <- colnames(CNOlist@inhibitors)
# We expand the 'inhibited' condition with
# the 'mutations' and four additive kinases 'KRAS', 'ERK1/2' and 'PDPK1'
conditions <- c('mutations', inhibited,
'ERK1_2', 'KRAS', 'PDPK1')
# For each condition
for (condition in conditions){
# SET THE INITIAL STATE OF EACH CONDITION
if(condition == 'mutations'){
# set stimuli conditions as a named vector of 1
# and assign it to CNOlist@stimuli
# --> TNF and IGF1 receptors
d <- c(1,1)
names(d) <- c('TNF', 'IGF1R')
# --> GOF mutations
s <- mutations_vector[mutations_vector == 1]
sp <- c(s,d)
sf <- rep(1, length(sp))
names(sf) <- names(sp)
CNOlist@stimuli <- t(as.matrix(sf))
# Set inhibitory conditions as a named vector of 1
# and assign it to CNOlist@inhibitors
# --> LOF mutations
ini <- mutations_vector[mutations_vector == 0]
i_new <- rep(1, length(ini))
names(i_new) <- names(ini)
# Modify the built-in object
CNOlist@inhibitors <- t(as.matrix(i_new))
}else if(condition == 'FLT3'){ # set stimuli conditions
# --> TNF and IGF1 receptors
d <- c(1,1)
names(d) <- c('TNF', 'IGF1R')
# --> all the GOF
# BUT we exclude FLT3 since we inhibit it!
s <- mutations_vector[mutations_vector == 1]
sp <- c(s,d)
sf <- rep(1, length(sp))
names(sf) <- names(sp)
sf <- sf[names(sf) != 'FLT3'] # exclude FLT3
# Modify the built-in object
CNOlist@stimuli <- t(as.matrix(sf))
# Set inhibitory condition
# --> FLT3i
in_val <- c(1)
names(in_val) <- condition
# --> LOF mutations
i <- mutations_vector[mutations_vector == 0]
i_new <- rep(1, length(i))
names(i_new) <- names(i)
# combine the two
i_trans <- c(i_new, in_val)
names(i_trans) <- c(names(i_new), names(in_val))
# modify the built-in object
CNOlist@inhibitors <- t(as.matrix(i_trans))}else{
## MUTATIONS + FLT3i + NODES ##
# set stimuli conditions
# --> TNF and IGF1 receptor
d <- c(1,1)
names(d) <- c('TNF', 'IGF1R')
# --> GOF mutations excluding FLT3
s <- mutations_vector[mutations_vector == 1]
sp <- c(s,d)
sf <- rep(1, length(sp))
names(sf) <- names(sp)
sf <- sf[names(sf) != 'FLT3'] #exclude FLT3
CNOlist@stimuli <- t(as.matrix(sf)) # modify the built-in object
# set inhibitory conditions
# --> LOFs mutations
i <- mutations_vector[mutations_vector == 0]
i_new <- rep(1, length(i))
names(i_new) <- names(i)
# --> FLT3i and new druggable node specified in condition
in_val <- c(1, 1)
names(in_val) <- c('FLT3', condition)
i_trans <- c(i_new, in_val) #combine them
names(i_trans) <- c(names(i_new), names(in_val))
CNOlist@inhibitors <- t(as.matrix(i_trans))
}
# RUN THE SIMULATION
matrix_sim_q <- list()
# Select the cellular model corresponding to the patient mutation type
if(!grepl('^TKD', patient)){
#print('selecting JMD model')
bString <- JMD_best
modelT0 <- JMD_model
simList <- JMDsimList
}else{
#print('Selecting TKD model')
bString <- TKD_best
modelT0 <- TKD_model
simList <- TKDsimList
}
modelCut <- cutModel(modelT0, list(bString)[[1]]) # cutting expanded gates to optimal bString
newSimList <- cutSimList(simList, list(bString)[[1]]) # creating a new simList from cutted model
indexList <- indexFinder(CNOlist, modelT0) # creating index list
# T0 simulation (before treatment)
simResT0 <- simulatorT1(CNOlist=CNOlist, model=modelCut, simList=newSimList,
indexList=indexList, mode = 0)
colnames(simResT0) <- modelT0$namesSpecies
# T1 simulation (after treatment)
simResT1 <- simulatorT1(CNOlist=CNOlist, model=modelCut, simList=newSimList,
indexList=indexList, mode = 1)
colnames(simResT1) <- modelT0$namesSpecies
simRes <- list(T0 = simResT0, T1 = simResT1) # store results in a list
# If there is some oscillatory node, set to 50
simRes$T0[is.na(simRes$T0)] <- 50
simRes$T1[is.na(simRes$T1)] <- 50
# Create a matrix of the simulation result
sim_matrix <- tibble(nodes = colnames(simRes$T0), T0 = simRes$T0[1,], T1 = simRes$T1[1,])
colnames(sim_matrix) <- c('nodes', 'T0', paste0('FLT3i+',condition,'i'))
sim_matrix <- sim_matrix %>%
column_to_rownames('nodes') %>% as.matrix
if(condition == 'mutations'){
sim_matrix_new <- sim_matrix
conditions_names <- c('T0', 'FLT3-ITD+TNF+IGF')
colnames(sim_matrix_new) <- conditions_names
}else{
conditions_names <- colnames(sim_matrix_new)
sim_matrix_new <- cbind(sim_matrix_new, sim_matrix[,2])
if(condition == 'FLT3'){
colnames(sim_matrix_new) <- c(colnames(sim_matrix_new)[1:length(colnames(sim_matrix_new))-1],
paste0(condition,'i'))
}else{
colnames(sim_matrix_new) <- c(colnames(sim_matrix_new)[1:length(colnames(sim_matrix_new))-1],
paste0('FLT3i+',condition,'i'))
}
}
}
# Remove the T0 column (not important)
sim_matrix_new <- sim_matrix_new[,-1]
# Store simulation result in a tibble
sim_tibble_p <- tibble(rownames_to_column(data.frame(sim_matrix_new), 'node'))
# Remove combinatorial treatments with oscillatory outcome
sim_tibble_p <- sim_tibble_p[,colSums(sim_tibble_p == 50) == 0]
# Create a file with model nodes' states for each patient
write_tsv(sim_tibble_p,
paste0('./results/patients_simulations/', patient, '_steady_states.tsv'))
# PHENOTYPES INFERENCE ACTIVITY
# Infer proliferation and apoptosis activity in all conditions
phens <- c('APOPTOSIS', 'PROLIFERATION')
for(phen_i in 1:length(phens)){
phen <- phens[phen_i]
# read ProxPath output
apoptosis_annotation <- read_tsv('./results/phenotypes_analysis/paths_to_apoptosis_proliferation.txt') %>%
select(node = QueryNode, phenotype = EndNode, apoptosis = Final_Effect)
# select the phenotype and remove duplicated nodes
apoptosis_annotation <- apoptosis_annotation %>%
filter(phenotype == phen) %>%
select(node, apoptosis) %>%
distinct()
# Join nodes activity in simulation with the regulatory role
# on the phenotype
if(grepl('JMD', patient)){
annotated <- left_join(sim_tibble_p, apoptosis_annotation %>%
select(node, apoptosis), by = 'node')
# keep only terminal nodes in the signal cascade
if(phen == 'APOPTOSIS'){
annotated <- annotated %>%
filter(node %in% c('MTOR', 'STAT5A'))
}else{
annotated <- annotated %>%
filter(node %in% c('PI3K', 'MAPK14', 'STAT5A', 'MTOR'))
}
}else{
annotated <- left_join(sim_tibble_p, apoptosis_annotation %>%
select(node, apoptosis), by = 'node')
if(phen == 'APOPTOSIS'){
annotated <- annotated %>%
filter(node %in% c('STAT5A', 'AKT'))
}else{
annotated <- annotated %>%
filter(node %in% c('PI3K', 'KRAS', 'CREB1', 'STAT5A'))
}
}
# modify the table for the phenotype barplot
annotated <- annotated %>% column_to_rownames('node')
annotated <- annotated*annotated$apoptosis
annotated$apoptosis <- NULL
pheno_df <- tibble(condition = str_replace_all(names(colSums(annotated)), '\\.', '+'),
value = colSums(annotated))
#revert the APOPTOSIS sign to create 'apoptosis inhibition'
if(phen == 'APOPTOSIS'){
pheno_df$value <- pheno_df$value * -1
}
barplot_apoptosis <- ggplot(pheno_df, aes(x = condition, y = value, fill)) +
geom_bar(stat = 'identity', position = 'dodge', width = 0.8) +
#labs(title = str_replace(conditions_sim, '\\.', '+'))+
xlab('') + ylab(phen)+
theme_classic() +
ylim(0,5)+
ggtitle(paste0(patient, ' ', phen))+
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1, size = 10),
axis.text.y = element_text(size = 10),
axis.title = element_text(size = 12))
patients_plot_list[[plot_idx]] <- barplot_apoptosis
plot_idx = plot_idx + 1
print(plot_idx)
# create a table for proliferation and apoptosis phenotype for heatmap
if(p_idx == 1){
if(phen == 'APOPTOSIS'){
apoptosis_df <- pheno_df
colnames(apoptosis_df) <- c('condition', patient)
}else{
proliferation_df <- pheno_df
colnames(proliferation_df) <- c('condition', patient)
}
}else{
if(phen == 'APOPTOSIS'){
colnames(pheno_df) <- c('condition', patient)
apoptosis_df <- full_join(apoptosis_df, pheno_df, by = 'condition')
}else{
colnames(pheno_df) <- c('condition', patient)
proliferation_df <- full_join(proliferation_df, pheno_df, by = 'condition')
}
}
}
patients_plot_list[[plot_idx]] <- plot_list
}
# save df for apoptosis and proliferation levels
# arrange all the plots in a grid
combined_plots <- plot_grid(plotlist = patients_plot_list, labels = "AUTO", ncol=8,nrow = 4)
write_rds(combined_plots, './results/patients_simulations/barplot_combinations.rds')
write_tsv(apoptosis_df, './results/patients_simulations/apoptosis_combination_patients.tsv')
write_tsv(proliferation_df, './results/patients_simulations/proliferation_combination_patients.tsv')
combined_plots <- readRDS('./results/patients_simulations/barplot_combinations.rds')
combined_plots
Figure S6E Each couple of bar plots show the in-silico apoptosis inhibition (right) and proliferation activation (left) levels in control and FLT3 inhibition conditions in combination with the knock-out of each of 10 kinases in FLT3 ITD-JMD and FLT3 ITD-JMD + TKD and FLT3 ITD-TKD patients.
i_plot = 1
plot_list <- list()
for(phenotype in c('apoptosis', 'proliferation')){
pheno_tb <- read_tsv(paste0('./results/patients_simulations/', phenotype, '_combination_patients.tsv'))
pheno_matrix <- pheno_tb %>% column_to_rownames('condition')
if(phenotype == 'apoptosis'){
myColor <- c('#557887', '#FFE2E2', 'red3')
}else{
myColor <- colorRampPalette(c("#FFE2E2", "red3"))(4)
}
p <- pheatmap(as.matrix(pheno_matrix),
color = myColor,
cellwidth = 20,
cellheight = 20,
cluster_rows = F,
cluster_cols = T,
border_color = 'white',
fontsize = 13,
width = 5,
heigth = 30,
silent = TRUE,
main = paste0('in silico ', phenotype,
ifelse(phenotype == 'apoptosis',
' inhibition ', ' '),
'levels'),
angle_col = c('90'))
plot_list[[i_plot]] <- p[[4]]
i_plot = i_plot + 1
}
plot_grid(plotlist = plot_list, ncol = 2)
Figure 5E Heatmap representing patient-specific in silico apoptosis inhibition (left panel) and proliferation levels (right panel) upon each simulation condition.