Shiva Retrospective Analysis

# Read in SHIVA raw data
merged_df <- readxl::read_excel("version_6.0/Supplementary_file/SHIVA_authorTable.xlsx"
                                  , sheet = 1
                                  , col_names = TRUE
 #                                , col_types = rep("text", times=24)
                                 )

# Import CNA data for the 496 profiled patients

cna_data <- readxl::read_excel("version_6.0/Supplementary_file/SHIVA_cna_data.xls"
                                , sheet = 1
                                , col_names = TRUE)

# Patient filtering
# ------------------------------------------------------------------------------
# All patients have CNA data, so we do not include any filter here. 
df_full_profile <- merged_df %>%
  # filter by patients that do not have a Complete profile: CNA are OK, no patients filtered there.
  # Filter those patient without a complete NGS profile (for SNV analysis)
  dplyr::filter(!NGS == "Not done") %>%
  # filter those patients with locations incompatible with the TCGA 
  dplyr::filter(!TCGA_tumor_ID == "") %>%
  # Remove patients that were not randomized because of problems
  # dplyr::filter(Reason.for.no.randomization == "ERROR")
  # Remove patients with tumor types  acyc, sclc, es, becuase have no CNA data in TCGA population
  dplyr::filter(!TCGA_tumor_ID %in% c("acyc", "sclc", "es"))
  
# Define tumor type list in the population
# -------------------------------------------------------------------------------
tumor_types <- as.character(unique(df_full_profile$TCGA_tumor_ID))
# calculate weights
tumor_weights <- table(as.character(df_full_profile$TCGA_tumor_ID))
# calculate frequency
tumor_freq <- prop.table(tumor_weights)

# Preview
# --------------------------------------------------------------------------------
datatable_withbuttons(df_full_profile, caption = "Patients with a complete profile")

data.frame(tumor_weights) %>% 
  dplyr::rename(Tumor = Var1, Count = Freq) %>%
  dplyr::mutate(Tumor = factor(Tumor
                               , level = Tumor[order(Count)])) %>%
  ggplot(aes(Tumor, Count)) +
    geom_bar(stat="identity", position=position_dodge() ) + 
    ggtitle("Number of patient for each tumor type in the SHIVA population") +
    coord_flip()

# DESIGN
SHIVA_population <- newCancerPanel(SHIVA_design)
#DOWNLOAD alternative
SHIVA_population <- getAlterations(SHIVA_population
                              , tumor_type = tumor_types
                              , repos = SHIVA_repo
                              )
#SIMULATE
SHIVA_population <- subsetAlterations(SHIVA_population)

#save a backup copy to speed up future analysis
saveRDS(SHIVA_population, file = "version_6.0/Temp/SHIVA_population.rds")

# check pop size from CNA
mya <- purrr::map_int(SHIVA_population@dataFull$copynumber$Samples, length) %>% 
  data.frame(count = .) %>%
  tibble::rownames_to_column("tumor_id") %>%
  dplyr::mutate(origin="CNA")

# check pop size from Mutations
myb <- purrr::map_int(SHIVA_population@dataFull$mutations$Samples, length) %>% 
  data.frame(count = .) %>%
  tibble::rownames_to_column("tumor_id") %>%
  dplyr::mutate(origin="mutations")

# calculate adjusting weight
# ------------------------------------------------------------------------------
# tumor_weights will be the vector to adjust the TCGA population
tumor_weights <- names(SHIVA_population@dataFull$copynumber$Samples) %>%
  purrr::map(function(tumor_name){
    
    mut = SHIVA_population@dataFull$mutations$Samples[[tumor_name]]
    cna = SHIVA_population@dataFull$copynumber$Samples[[tumor_name]]
    
    # find intersection of patients
    output <- length(mut %in% cna)
    names(output) = tumor_name
    
    return(output)
    
  }) %>% flatten_int

# check pop size used internally for adjusting TCGA
myc <- tumor_weights %>% 
  data.frame(count = .) %>%
  tibble::rownames_to_column("tumor_id") %>%
  dplyr::mutate(origin="new adjusting weight")

# generate a plot
rbind(mya, myb, myc) %>%
  ggplot(aes(tumor_id, count, fill=origin)) +
  geom_bar(stat="identity", position= position_dodge()) +
    coord_flip() + 
    ggtitle("number of patients per alteration type from SHIVA.")

# save the number of entries in the TCGA beofere any filtering
TCGA_original_number <- nrow(unique(TCGA_population@dataFull$mutations$data))

# Clean TCGA popualation running a cosmic filter
cosmic <- readr::read_csv("version_6.0/external_resources/cosmic_pathogenic_SHIVApanel")

cosmic_df <- data.frame(gene_symbol = cosmic$gene
                        , amino_acid_change = sub("^p." , "" , cosmic$AA.Mutation)
                        , stringsAsFactors = FALSE)
TCGA_population <- filterMutations(TCGA_population , filtered = cosmic_df , mode = "keep")

# SETTINGS
# ==============================================================================
# DEFINE GENES: define genes we are interested for the CNA analysis
CNA_genes <- TCGA_population@arguments$panel %>%
  dplyr::filter(alteration == "CNA") %>%
  dplyr::select(gene_symbol) %>%
  .[[1]]

# PATH: set global variable with the path to.rds files with CNA raw data
# PATH_2CNA_FILES <- "/Users/ale/IEO/projects/IEO_research/R_package_panel_and_projection/functionfactory/cna_firehose_rds/"
# Set a relative path for reproducibility
PATH_2CNA_FILES <- "version_6.0/Temp/cna_firehose_rds/"

# helper function that reads in a firehose RDS file, filters by gene, and apply the threshold configured
# ==============================================================================
get_cna_long <- function(tumor_type, PATH_2rds, genes_chr, del_thr = 1.1, amp_thr){
  
  
  # CHECKS
  # --------------------------------------------------------------------------------
  # check that PATH actually leads to a dir
  if(!dir.exists(PATH_2rds)){ stop("path to a dir that cannot be found")}

  # Check that all the files are found for each tumor type required in the analysis
  input_lgl <- purrr::map_lgl(tumor_types, ~ file.exists(paste0(PATH_2rds, toupper(.) , "_PTD.rds")))
  
  if(!all(input_lgl)){
    stop(paste("Files for some of the tumor types where not found:", tumor_types[-which(input_lgl)]))
  }
 
  # Checks
  if(!is.character(tumor_types) | is.null(tumor_types) | length(tumor_types) ==0){
    stop(paste("Tumor type is not confirm", tumor_types))
  }


  # Define an internal function to do the job
  # -------------------------------------------------------------------------------
  get_CNA_long_onetumor <- function(one_tumor, PATH_2rds, genes_chr, del_thr = 1.1, amp_thr){
    
    # read in CNA file
    # myfile <- readRDS(paste0(PATH_2rds, toupper(one_tumor) , "_PTD.rds"))
    # if(!file.exists(myfile)){
    #   message(myfile)
    #   return(NULL)
    # }
    # mat <- tryCatch(readRDS(myfile) , error = function(e) stop(myfile))
    
    mat <- readRDS(paste0(PATH_2rds, toupper(one_tumor) , "_PTD.rds"))
    
    # extract Patients
    pats <- colnames(mat)

    # Filter for gene of interest
    mat <- mat[ genes_chr , ]

    # Melt and adjust
    output <- reshape2::melt(mat , value.name = "CNAvalue")
    colnames(output) <- c("gene_symbol" , "case_id" , "CNAvalue")
      output$tumor_type <- one_tumor
    # Choose cutoff values
    # -----------------------------------------------------------------------------
    # With default threshold data are now like this: 
        # [0 , 1.1) "deletion"
        # 1.1 "shallow deletion"
        # 2 "normal"
        # 2.9 "shallow amplification"
        # (2.9 , Inf] "amplification"
    output$CNA <- ifelse( output$CNAvalue <del_thr 
                    , "deletion" 
                    , ifelse(output$CNAvalue > amp_thr , "amplification" , "normal"))
    # select
    output <- output %>%
      dplyr::select(gene_symbol , CNA , case_id , tumor_type, CNAvalue) %>%
      dplyr::mutate(gene_symbol = as.character(gene_symbol)) %>%
      dplyr::mutate(case_id = as.character(case_id)) 
    
    return(output)
  }

  # Main body of the function
  # -----------------------------------------------------------------------------
  cna_long_df <- purrr::map(tumor_types
                             , ~ get_CNA_long_onetumor(., PATH_2rds, genes_chr, del_thr, amp_thr)
                             ) %>%
    dplyr::bind_rows()
    
  return(cna_long_df)  

}
  

# RUN MAIN FUNCTION
# ==============================================================================
cna_long_df <- get_cna_long(tumor_types, PATH_2CNA_FILES, CNA_genes, del_thr = 1.1, amp_thr = 4)


# Create Sample
# ----------------------
# get list of tumors
tumor_chr <-  cna_long_df$tumor_type %>% 
  unique %>% 
  as.character

# prepare funtion to extract sample names for each tumor
sampleXtumor <- function(x, df){
  df %>% 
    dplyr::select(case_id, tumor_type) %>%
    dplyr::filter(tumor_type %in% x) %>% .[["case_id"]] %>%
    as.character() %>%
    unique()
}
# Extract samples for each tumors
Samples <- purrr::map(tumor_chr, ~ sampleXtumor(., df=cna_long_df))
names(Samples) <- tumor_chr

# Check if there are problems. 
testthat::expect_equal(length(unlist(Samples)), length(unique(cna_long_df$case_id)))

# Create list element for CNA for the DataFull argument of the panel (panel@DataFull)
# ==============================================================================
copynumber <- list(  data = data.frame(cna_long_df)
                   , Samples = Samples)


# Preview changes before committing
#TCGA_population@dataFull$copynumber <- copynumber
TCGA_cna_pop_comparison <- purrr::map2( 
    .x = list( TCGA_population@dataFull$copynumber$Samples, copynumber$Samples )
  , .y = list( "TCGA", "Firehose" )
  , .f = ~ purrr::map_int(.x, length) %>%
      data.frame(n_patients = .) %>%
      tibble::rownames_to_column("tumor_type") %>%
      dplyr::mutate(source = .y)
  ) %>% 
  dplyr::bind_rows() %>%
  ggplot(aes(tumor_type, n_patients, fill=source)) + 
    geom_bar(stat="identity", position= position_dodge()) +
    coord_flip() + 
    ggtitle("number of patients in the CNA datasets")

TCGA_cna_pop_comparison

# Apply changes
TCGA_population@dataFull$copynumber <- copynumber
TCGA_population <- subsetAlterations(TCGA_population)

purrr::map2( 
    .x = list( TCGA_population@dataFull$copynumber$Samples, SHIVA_population@dataFull$copynumber$Samples )
  , .y = list( "TCGA", "SHIVA" )
  , .f = ~ purrr::map_int(.x, length) %>%
      data.frame(n_patients = .) %>%
      tibble::rownames_to_column("tumor_type") %>%
      dplyr::mutate(source = .y)
  ) %>% 
  dplyr::bind_rows() %>%
  ggplot(aes(reorder(tumor_type, n_patients), n_patients, fill=source)) + 
    geom_bar(stat="identity", position= position_dodge()) +
    scale_fill_manual(values=MYPALETTE(10)[c(10,2)]) + 
    coord_flip() + 
    ggtitle("Number of samples available from TCGA and SHIVA group by tumor types")

# Run simulation
# ------------------------------------------------------------------------------
N_SIMULATIONS = 200

simul <- mclapply(1:N_SIMULATIONS , function(x) {
            cov <- coveragePlot(TCGA_population 
                    , alterationType=c("mutations" , "copynumber") 
                    , collapseByGene=TRUE
                    , tumor.weights = tumor_weights
                    , maxNumAlt = 7
                    , noPlot=TRUE
                     )
            return(cov$plottedTable/cov$Samples[1])
        } , mc.cores=4)

## Prepare function for plotting detected fraction with Confidenti intervals
# ------------------------------------------------------------------------------

# Convert to dataframe
simul <- data.frame(matrix(unlist(simul)
                        , nrow=length(simul)
                        , byrow=T )
                 , stringsAsFactors = FALSE)

# GET STATS
# ------------------------------------------------------------------------------
# Calculate for each n° of alteration, the mean
detectionPower <- colMeans(simul)

# Calculate lower CI for bootstrap
leftCI <- apply(simul,2,function(x){left <- quantile(x, 0.025)})

# Calculate upper CI for bootstrap
rightCI <- apply(simul,2,function(x){right <- quantile(x, 0.975)})

# Preapre for column with # alterations
n_alterations <- 1:length(detectionPower)

# Put them all together in a data frame
dtpow <- data.frame(n_alterations, detectionPower, leftCI, rightCI)

# Generate the plot
TCGA_detection_image <- dtpow %>%
  ggplot(aes(x=n_alterations, y=detectionPower)) + 
  geom_bar(stat="identity" , position=position_dodge(), fill=MYPALETTE(10)[2]) + 
  scale_y_continuous(limits = c(0,1)) +
  scale_x_continuous(breaks = n_alterations, trans = "reverse") +
  geom_errorbar(aes(ymax=rightCI
                    , ymin=leftCI)
                , position=position_dodge(0.9)
  ) + 
  ggtitle("Overall detected fraction in TCGA Population") +
  annotate("text"
           , x=n_alterations
           , y=detectionPower+0.15
           , label= sapply(detectionPower*100
                           , function(x) paste(c(round(x, digits=2), "%")
                                               , collapse=" ")
           )
  ) +
  coord_flip()

# Show image
TCGA_detection_image

# Preview table
datatable_withbuttons(dtpow, caption = "detected fraction TCGA")

# calculate shiva detected fraction
shiva_detect <-suppressWarnings(
                  coveragePlot(SHIVA_population
                   , alterationType = c("mutations", "copynumber")
                   , collapseByGene = TRUE
                   , maxNumAlt = 7
                   , noPlot= TRUE
                  )
              )

# calculate frequencies
shiva_freq <- as.numeric(shiva_detect$plottedTable/shiva_detect$Samples[1])

# Preapre for column with # alterations
n_alterations <- 1:length(shiva_freq)
  
# Put them all together in a data frame
dtpow_shiva <- data.frame(n_alterations, shiva_freq= shiva_freq)
  
# left CI and right CI
dtpow_shiva$leftCI <- purrr::map_dbl(shiva_freq, ~ leftCIprop(., sum(tumor_weights)))
dtpow_shiva$rightCI <- purrr::map_dbl(shiva_freq, ~ rightCIprop(., sum(tumor_weights)))

# Generate the plot
SHIVA_detection_image <- dtpow_shiva %>%
  ggplot(aes(x=n_alterations, y=shiva_freq)) + 
    geom_bar(stat="identity" , position=position_dodge(), fill=MYPALETTE(10)[10]) + 
    scale_y_continuous(limits = c(0,1)) +
    scale_x_continuous(breaks = n_alterations, trans = "reverse") +
    geom_errorbar(aes(   ymax= leftCI
                       , ymin= rightCI)
                   , position=position_dodge(0.9)
    ) + 
    ggtitle("Detected fraction SHIVA population") +
    annotate("text"
             , x=n_alterations
             , y=shiva_freq+0.15
             , label= sapply(shiva_freq*100
                             , function(x) paste(c(round(x, digits=2), "%")
                                                 , collapse=" ")
             )
    ) +
    coord_flip()
  
#Export as a svg
#img_file = "../Figures/fig4B.svg"
#ggsave(filename=img_file, plot=SHIVA_detection_image, device = "svg")
  
# Show image
SHIVA_detection_image

# Preview table
datatable_withbuttons(dtpow_shiva)

fig4a <- gridExtra::grid.arrange(TCGA_detection_image, SHIVA_detection_image)

#Export as a svg
#img_file = "../Figures/fig4A.svg"
#ggsave(filename=img_file, plot=fig4a, device = "svg", width=10, height = 8)


# Preview image
fig4a

###############################################################################
# Fetch SNV data from TCGA population 
# -----------------------------------------------------------------------------
tcga_wide <- purrr::map(1:N_SIMULATIONS
           , ~ cpFreq(TCGA_population, alteration="mutations", tumor.weights = tumor_weights) %>% 
                tidyr::spread(gene_symbol, freq)
           ) %>% bind_rows

# GET STATS
# ------------------------------------------------------------------------------
# Calculate for each n° of alteration, the mean
freq <- colMeans(tcga_wide)

# Calculate lower CI for bootstrap
leftCI <- apply(tcga_wide,2,function(x){left <- quantile(x, 0.025)})

# Calculate upper CI for bootstrap
rightCI <- apply(tcga_wide,2,function(x){right <- quantile(x, 0.975)})


# Put them all together in a data frame
tcga_snv <- data.frame(freq, leftCI, rightCI, population="TCGA") %>%
  tibble::rownames_to_column("gene") 

###############################################################################
# Fetch CNA data from SHIVA population 
# -----------------------------------------------------------------------------
temp <- cpFreq(SHIVA_population, alteration="mutations")
temp <- temp[which(temp[["gene_symbol"]] %in% SHIVA_population@arguments$panel$gene_symbol),]

# get CI
leftCI <- purrr::map_dbl(temp$freq, ~ leftCIprop(., sum(tumor_weights)))
rightCI <- purrr::map_dbl(temp$freq, ~ rightCIprop(., sum(tumor_weights)))

# get DF
shiva_snv <- temp %>%
  dplyr::rename(gene = gene_symbol) %>%
  dplyr::mutate(leftCI=leftCI
                , rightCI=rightCI
                ) %>%
  dplyr::mutate(population = "SHIVA")

###############################################################################
# put the plots together
# -----------------------------------------------------------------------------

# Put them together
merged_snv <- rbind(tcga_snv, shiva_snv)

# bind SHIVA and TCGA
fig4b <- merged_snv %>%
  dplyr::mutate(gene= factor(gene, levels = unique(gene[order(freq)]))) %>%
  ggplot(aes(x=gene, y=freq, fill=population)) + 
  geom_bar(stat="identity", position=position_dodge()) +
  scale_fill_manual(values=MYPALETTE(10)[c(2,10)]) +
  #scale_y_continuous(limits = c(0,.7)) +
  geom_errorbar(aes(ymax=rightCI
                    , ymin=leftCI)
                , position=position_dodge(0.9)
  ) + 
 ggtitle("SNV frequency comparison in the two populations") +   
 coord_flip()


#img_file = "../Figures/fig4B.svg"
#ggsave(filename=img_file, plot=fig4b, device = "svg", width=10, height = 8)

# Preview
fig4b

# Fetch CNA data from TCGA population 
# ------------------------------------------------------------------------------
simul_amplifications <- lapply(1:50, function(x) {
            cpFreq(TCGA_population, alteration="copynumber", tumor.weights = tumor_weights) %>%
                dplyr::select(amplification) %>% 
                tibble::rownames_to_column("gene") %>%
                tidyr::spread(gene, amplification)
        }) 

simul_amplifications <- simul_amplifications %>%  bind_rows()
  
simul_deletions <- lapply(1:50, function(x) {
            cpFreq(TCGA_population, alteration="copynumber", tumor.weights = tumor_weights) %>%
                dplyr::select(deletion) %>% 
                tibble::rownames_to_column("gene") %>%
                tidyr::spread(gene, deletion)
        }) %>%  bind_rows()
  
# GET STATS
# ------------------------------------------------------------------------------
TCGA_amplifications_df <- get_bootstrap_stats(simul_amplifications) %>%
  tibble::rownames_to_column("gene") %>%
  dplyr::mutate(population = "TCGA")

TCGA_deletion_df <- get_bootstrap_stats(simul_deletions) %>%
  tibble::rownames_to_column("gene") %>%
  dplyr::mutate(population = "TCGA")


# Fetch CNA data from SHIVA population 
# ------------------------------------------------------------------------------

# Amplifications
# ------------
shiva_cna_amplification <- cpFreq(SHIVA_population, alteration="copynumber") %>%
  tibble::rownames_to_column("gene") %>%
  dplyr::select(-deletion, -normal) %>%
  dplyr::rename(freq = amplification) %>%
  dplyr::mutate(population = "SHIVA") 

# Calculate left CI and right CI
shiva_cna_amplification$leftCI <- purrr::map_dbl(shiva_cna_amplification$freq, ~ leftCIprop(., sum(tumor_weights)))
shiva_cna_amplification$rightCI <- purrr::map_dbl(shiva_cna_amplification$freq, ~ rightCIprop(., sum(tumor_weights)))

# last formattings 
shiva_cna_amplification <-  shiva_cna_amplification %>%
  dplyr::mutate(row_counter = 1:nrow(shiva_cna_amplification)) %>%
  dplyr::select(gene, row_counter, freq, leftCI, rightCI, population)


# Deletion
# ------------
shiva_cna_deletion <- cpFreq(SHIVA_population, alteration="copynumber") %>%
  tibble::rownames_to_column("gene") %>%
  dplyr::select(-amplification, -normal) %>%
  dplyr::rename(freq = deletion) %>%
  dplyr::mutate(population = "SHIVA")

# Calculate left CI and right CI
shiva_cna_deletion$leftCI <- purrr::map_dbl(shiva_cna_deletion$freq, ~ leftCIprop(., sum(tumor_weights)))
shiva_cna_deletion$rightCI <- purrr::map_dbl(shiva_cna_deletion$freq, ~ rightCIprop(., sum(tumor_weights)))

# last formattings 
shiva_cna_deletion <-  shiva_cna_deletion %>%
  dplyr::mutate(row_counter = 1:nrow(shiva_cna_deletion)) %>%
  dplyr::select(gene, row_counter, freq, leftCI, rightCI, population)

# Prepare filter
# Removed genes not considered in the panel for amplifications or deletions
# ------------
amplification_filter <- SHIVA_design %>%
  dplyr::filter(exact_alteration == "amplification") %>%
  dplyr::select(gene_symbol) %>%
  .[[1]] %>% 
  as.character()

deletion_filter <- SHIVA_design %>%
  dplyr::filter(exact_alteration == "deletion") %>%
  dplyr::select(gene_symbol) %>%
  .[[1]] %>% 
  as.character()


# Merge datasets
# ------------
merged_cna_amplifications <- rbind(TCGA_amplifications_df, shiva_cna_amplification) %>% 
  dplyr::filter(gene %in% amplification_filter) 
merged_cna_deletion <- rbind(TCGA_deletion_df, shiva_cna_deletion) %>% 
  dplyr::filter(gene %in% deletion_filter) 


# bind SHIVA and TCGA | Amplification
# -----
CNA_plot_amplification <- merged_cna_amplifications %>%
  dplyr::mutate(gene= factor(gene, levels = unique(gene[order(freq)]))) %>%
  ggplot(aes(x=gene, y=freq, fill=population)) + 
  geom_bar(stat="identity", position=position_dodge()) +
  scale_fill_manual(values=MYPALETTE(10)[c(10,2)]) +
  #scale_y_continuous(limits = c(0,.7)) +
  geom_errorbar(aes(ymax=rightCI
                    , ymin=leftCI)
                , position=position_dodge(0.9)
  ) + 
 ggtitle("CNA Amplification frequency comparison in the two populations") +   
 coord_flip()


CNA_plot_deletion <- merged_cna_deletion %>%
  dplyr::mutate(gene= factor(gene, levels = unique(gene[order(freq)]))) %>%
  ggplot(aes(x=gene, y=freq, fill=population)) + 
  geom_bar(stat="identity", position=position_dodge()) +
  scale_fill_manual(values=MYPALETTE(10)[c(10,2)]) +
  #scale_y_continuous(limits = c(0,.7)) +
  geom_errorbar(aes(ymax=rightCI
                    , ymin=leftCI)
                , position=position_dodge(0.9)
  ) + 
 ggtitle("CNA Deletion frequency comparison in the two populations") +   
 coord_flip()


# Make the plot
fig4c <- gridExtra::grid.arrange(CNA_plot_amplification, CNA_plot_deletion, ncol=1, heights = c(20,5))

#Export as a svg
#img_file = "../Figures/fig4C.svg"
#ggsave(filename=img_file, plot=fig4c, device = "svg", width=10, height = 8)

# preview plot
fig4c

# ALLOCATION FUNCTION with not predefined rule
# ------------------
# Create function to perform the following tasks:
# 1. extract a weighted (according to SHIVA population) sample from the population
# 2. Simulate patient allocation according the defined schema
allocate_patients_fn <- function(population, weights=FALSE){
    
  #extract patient data from the panel
  if(weights){
    patient_lt <- dataExtractor(population
                  , alterationType = c("copynumber", "mutations")
                  , collapseByGene = TRUE 
                  , tumor.weights = weights
                  )
  } else {
    patient_lt <- dataExtractor(population
                  , alterationType = c("copynumber", "mutations")
                  , collapseByGene = TRUE 
                  )
  }
   
  #select general dataframe
  patient_df <- patient_lt$data
                
  # replace / with _
  patient_df[,2] <- gsub("/", "_", patient_df[,2])
  
  # generate a frequency table that indicates, for each patient, elibibility to group allocation, as well as tumor_type
  # one patient per row.
  patient_allocation_df <- patient_df %>%
    dplyr::select(group, case_id, tumor_type) %>%
    dplyr::mutate(found = TRUE) %>%
    unique() %>%
    tidyr::spread(group, found) 
  
  # SIMULATE PATIENTS ALLOCATION ACCORDING TO A SPECIFIC ORDER
  # ----------------------------------------
  # select ONLY PIK3_AKT_mTOR compatible patients
  PIK3_AKT_mTOR_assigned <- patient_allocation_df %>%
    dplyr::filter(is.na(RAF_MEK)) %>%
    dplyr::filter(PIK3_AKT_mTOR == TRUE)

  # select ONLY RAF_MEK compatible patients
  RAF_MEK_assigned <- patient_allocation_df %>%
    dplyr::filter(is.na(PIK3_AKT_mTOR)) %>%
    dplyr::filter(RAF_MEK == TRUE)
  
  # Select patients that would be allocated to either
  eligible_both <- patient_allocation_df %>%
    dplyr::filter(RAF_MEK == TRUE) %>%
    dplyr::filter(PIK3_AKT_mTOR == TRUE)
  
  # freq of patient allocation in each treatment group
  df <- data.frame(PIK3_AKT_mTOR = nrow(PIK3_AKT_mTOR_assigned) / nrow(patient_allocation_df)
             , RAF_MEK = nrow(RAF_MEK_assigned) / nrow(patient_allocation_df)
             , eligible_both = nrow(eligible_both) / nrow(patient_allocation_df)
  )
  return(df)
}

# TCGA
# ------------------------------------------------------------------------------
bootstrapping_df <- replicate(500
                                , allocate_patients_fn(TCGA_population, tumor_weights) 
                                , simplify=FALSE
                                ) %>% do.call("rbind", .)

# calculae the mean, sd and CI for each treatment group
bootstrappingMeans <- colMeans(bootstrapping_df)
lowerCI <- apply( bootstrapping_df , 2 , function(x) quantile(x, 0.025))
upperCI <- apply( bootstrapping_df , 2 , function(x) quantile(x, 0.975))


# Prepare dataframe for plotting
# ----------------------------
# covert to dataframe
tcga_df <- data.frame(mean=bootstrappingMeans
           , lowerSD = NA
           , upperSD = NA
           , lowerCI = lowerCI
           , upperCI = upperCI
           ) %>%
  tibble::rownames_to_column("group_treatment") %>%
  dplyr::arrange(mean) %>%
  dplyr::mutate(dataset="TCGA")



# SHIVA
# ------------------------------------------------------------------------------
shiva_df <- allocate_patients_fn(SHIVA_population)  %>%
  tidyr::gather(group_treatment, mean) %>%
  dplyr::mutate(lowerSD = NA
                , upperSD = NA
                , lowerCI = leftCIprop(mean, sum(tumor_weights))
                , upperCI = rightCIprop(mean, sum(tumor_weights))
                , dataset = "SHIVA"
                )


# merge dataframes together
merge_df <- rbind(shiva_df, tcga_df) %>%
  dplyr::rename(freq = mean)

# ----------------------------
#  PLOT2 - Same side barplot
# barplot, with only the genes in the panel
fig4E <- merge_df %>%
  dplyr::mutate(origin=paste(as.character(group_treatment), as.character(dataset), sep=":")) %>%
  dplyr::mutate(row_counter = as.integer(factor(dataset, labels = c(0,1)))) %>%
  ggplot(aes(x=origin, y=freq, fill=dataset)) + 
  geom_bar(stat="identity") + 
  scale_fill_manual(values=MYPALETTE(10)[c(10,2)]) + 
  geom_errorbar(aes(min=lowerCI, max=upperCI)) + 
  facet_grid(. ~ group_treatment, space="free", scales="free") +
  theme(axis.text.x=element_blank()
        , axis.ticks.x= element_blank()
        , axis.title.x = element_blank()
        ) +
    labs(title="Comparison between expected allocation to molecular pathways")


fig4E <- fig4E + geom_text(aes(x=(row_counter)
                  , y=upperCI+.05
                  , label= sapply(freq*100
                             , function(x) paste(c(round(x, digits=2), "%")
                                                 , collapse=" ")
                            )
                  )
)

#Export as a svg
#img_file = "../Figures/fig4E.svg"
#ggsave(filename=img_file, plot=fig4E, device = "svg", width=10, height = 8)

# preview image
fig4E

merge_df %>%
  dplyr::select(-lowerSD, -upperSD) %>%
  datatable_withbuttons

################################################################################
# ADD plot by Drug
################################################################################

# get the numer of samples used for the next computation
n_samples <- coveragePlot(TCGA_population 
                                  , alterationType=c("mutations" , "copynumber") 
                                  , collapseByGene=TRUE
                                    , tumor.weights = tumor_weights
                                    , maxNumAlt = 1
                                    , noPlot=TRUE
                                    , grouping = "drug"
                                     )$Sample[1]

# TCGA
# ------------------------------------------------------------------------------
bootstrapping_df <- replicate(500
          , coveragePlot(TCGA_population 
                                  , alterationType=c("mutations" , "copynumber") 
                                  , collapseByGene=TRUE
                                    , tumor.weights = tumor_weights
                                    , maxNumAlt = 1
                                    , noPlot=TRUE
                                    , grouping = "drug"
                                     ) %>% 
            .$plottedTable %>%
            data.frame %>%
            tibble::rownames_to_column("drug") %>%
            dplyr::rename(count=X1) %>%
            dplyr::mutate(freq = count/n_samples) %>%
            dplyr::select(-count) %>%
            dplyr::mutate(drug= gsub( "drug: ", "", drug) ) %>%
            tidyr::spread(drug, freq) %>%
            dplyr::mutate(Sorafenib = if (exists('Sorafenib', where = .)) Sorafenib else 0) %>%
            dplyr::mutate(Dasatinib = if (exists('Dasatinib', where = .)) Dasatinib else 0)
          , simplify=FALSE
  ) %>% do.call("rbind", .)


# calculae the mean, sd and CI for each treatment group
bootstrappingMeans <- colMeans(bootstrapping_df)
lowerCI <- apply( bootstrapping_df , 2 , function(x) quantile(x, 0.025))
upperCI <- apply( bootstrapping_df , 2 , function(x) quantile(x, 0.975))

# Prepare dataframe for plotting
# ----------------------------
# covert to dataframe
tcga_df <- data.frame(mean=bootstrappingMeans
           , lowerSD = NA
           , upperSD = NA
           , lowerCI = lowerCI
           , upperCI = upperCI
           ) %>%
  tibble::rownames_to_column("group_treatment") %>%
  dplyr::arrange(mean) %>%
  dplyr::mutate(dataset="TCGA")



# SHIVA
# ------------------------------------------------------------------------------
shiva_df <- coveragePlot(SHIVA_population 
                    , alterationType=c("mutations" , "copynumber") 
                    , collapseByGene=TRUE
                    , tumor.weights = tumor_weights
                    , maxNumAlt = 1
                    , noPlot=TRUE
                    , grouping = "drug"
                     ) %>% .$plottedTable %>%
  data.frame %>%
  tibble::rownames_to_column("drug") %>%
  dplyr::rename(count=X1) %>%
  dplyr::mutate(freq = count/n_samples) %>%
  dplyr::select(-count) %>%
  dplyr::mutate(drug= gsub( "drug: ", "", drug) ) %>%
  tidyr::spread(drug, freq) %>%
  dplyr::mutate(Sorafenib = if (exists('Sorafenib', where = .)) Sorafenib else 0) %>%
  dplyr::mutate(Dasatinib = if (exists('Dasatinib', where = .)) Dasatinib else 0)   %>%
  tidyr::gather(group_treatment, mean) %>%
  dplyr::mutate(lowerSD = NA
                , upperSD = NA
                , lowerCI = map_dbl(mean, ~ leftCIprop(., sum(tumor_weights)))
                , upperCI = map_dbl(mean, ~ rightCIprop(., sum(tumor_weights)))
                , dataset = "SHIVA"
                )
# merge dataframes together
merge_df <- rbind(shiva_df, tcga_df) %>%
  dplyr::rename(freq = mean) %>%
  mutate_cond(group_treatment %in% "Lapatinib+Trastuzumab", group_treatment = "Lapatinib+\nTrastuzumab")


# ----------------------------
#  PLOT2 - Same side barplot
# barplot, with only the genes in the panel
fig4F <- merge_df %>%
  dplyr::mutate(origin=paste(as.character(group_treatment), as.character(dataset), sep=":")) %>%
  dplyr::mutate(row_counter = as.integer(factor(dataset, labels = c(0,1)))) %>%
  ggplot(aes(x=origin, y=freq, fill=dataset)) + 
  geom_bar(stat="identity") + 
  scale_fill_manual(values=MYPALETTE(10)[c(10,2)]) + 
  geom_errorbar(aes(min=lowerCI, max=upperCI)) + 
  facet_grid(. ~ group_treatment, space="free", scales="free") +
  theme(axis.text.x=element_blank()
        , axis.ticks.x= element_blank()
        , axis.title.x = element_blank()
        ) +
    labs(title="Comparison between expected allocation to molecular pathways") 
# +
#   coord_flip()


#Export as a svg
#img_file = "../Figures/fig4F.svg"
#ggsave(filename=img_file, plot=fig4F, device = "svg")

# Preview
fig4F

# preview table
merge_df %>%
  dplyr::select(-lowerSD, -upperSD) %>%
  datatable_withbuttons()