Dashboards/msa/global.R

# ***************************
# **   I M P O R T A N T   **
# ***************************

# DO NOT USE OR MODIFY THIS.
# USE THE ONE IN THE 'Dashboards'
# REPO





library(shinycssloaders)
library(DT)
library(NGLVieweR)
library(hash)

# FIXME This is slow and should happen *once only*
#source("~/git/LSHTM_analysis/scripts/Header_TT.R")

# FIXME: these are needed but slow to load every time
#  source("~/git/LSHTM_analysis/config/alr.R")
#  source("~/git/LSHTM_analysis/config/gid.R")

#  source("~/git/LSHTM_analysis/config/pnca.R")
#  source("~/git/LSHTM_analysis/config/rpob.R")
#  source("~/git/LSHTM_analysis/config/katg.R")

# TODO: this is TEMPORARY and will shortly get replaced with a target picker
# that'll reload everything when changing targets. the lapply() is *much* quicker
# than previous approaches

# system.time({
source("~/git/LSHTM_analysis/scripts/Header_TT.R")

load_target_globals=function(target){
  cat(paste0("Reloading Target: ", target))
  source(paste0("~/git/LSHTM_analysis/config/", target, ".R")) # load per-target config file

  invisible(assign(paste0(target, "_merged_df3"), read.csv(paste0("~/git/Misc/shiny_dashboard/data/",target,"-merged_df3.csv")), envir = .GlobalEnv))
  invisible(assign(paste0(target, "_merged_df2"), read.csv(paste0("~/git/Misc/shiny_dashboard/data/",target,"-merged_df2.csv")), envir = .GlobalEnv))
  invisible(assign(paste0(target, "_corr_df_m3_f"), read.csv(paste0("~/git/Misc/shiny_dashboard/data/",target,"-corr_df_m3_f.csv")), envir = .GlobalEnv))
  invisible(assign(paste0(target, "_lin_lf"), read.csv(paste0("~/git/Misc/shiny_dashboard/data/",target,"-lin_lf.csv")), envir = .GlobalEnv))
  invisible(assign(paste0(target, "_lin_wf"), read.csv(paste0("~/git/Misc/shiny_dashboard/data/",target,"-lin_wf.csv")), envir = .GlobalEnv))
  lapply(
    c(
      "duet",
      "mcsm_lig",
      "foldx",
      "deepddg",
      "dynamut2",
      "consurf",
      "snap2",
      "provean",
      "dist_gen",
      "mcsm_ppi2"#,
      #"mcsm_na"
    ), function(x){
      wf_filename=paste0("~/git/Misc/shiny_dashboard/data/", tolower(gene), "-wf_", x ,".csv")
      wf_var=paste0("wf_",x)
      if (file.exists(wf_filename)){
        invisible(assign(wf_var,read.csv(wf_filename), envir = .GlobalEnv)) # FILTH
      }
      lf_filename=paste0("~/git/Misc/shiny_dashboard/data/", tolower(gene), "-lf_", x ,".csv")
      lf_var=paste0(target, "_lf_",x)
      if (file.exists(lf_filename)){
        invisible(assign(lf_var,read.csv(lf_filename), envir = .GlobalEnv)) # FILTH
      }
    }
  )
}
# populate target-specific *_unified_msa vars
load_msa_global=function(gene){
  drug=target_map[[gene]]
  in_filename_msa = paste0(tolower(gene), "_msa.csv")
  infile_msa = paste0("~/git/Data/", drug, "/output/", in_filename_msa)
  print(infile_msa)
  msa1 = read.csv(infile_msa, header = F)
  msa_seq = msa1$V1

  infile_fasta = paste0("~/git/Data/", drug, "/input/", tolower(gene), "2_f2.fasta")
  print(infile_fasta)
  msa2 = read.csv(infile_fasta, header = F)
  wt_seq = msa2$V1

  target_name=paste0(gene, '_unified_msa')
  #print(target_name)
  invisible(assign(target_name, list(msa_seq = msa_seq, wt_seq = wt_seq), envir = .GlobalEnv))
}

#### Local Functions ####


# FIXME: passing in the per-plot params is broken
lin_sc=function(x, all_lineages = F, ...){
  lf_var = get(paste0(x,"_lin_lf"))
  wf_var = get(paste0(x,"_lin_wf"))
  cowplot::plot_grid(lin_count_bp_diversity(wf_var, all_lineages, ...), lin_count_bp(lf_var, all_lineages, ...))
}


options(shiny.port = 8000)
options(shiny.host = '0.0.0.0') # This means "listen to all addresses on all interfaces"
options(shiny.launch.browser = FALSE)
options(width=120)
options(DT.options = list(scrollX = TRUE))


################ STATIC GLOBALS ONLY ##############
# never quite sure where "outdir" gets set :-|

# using dataframes instead of lists lets us avoid use of map()
plot_functions_df=data.frame(
  tab_name=c(
    "LogoP SNP",
    "Lineage Sample Count",
    "Site SNP count",
    "Stability SNP by site",
    "DM OM Plots",
    "Correlation",
    "Lineage Distribution",
    "Consurf",
    "LogoP OR",
    "LogoP ED"
  ),
  plot_function=c(
    "LogoPlotSnps",
    "lin_sc",
    "site_snp_count_bp",
    "bp_stability_hmap",
    "lf_bp2",
    "my_corr_pairs",
    "lineage_distP",
    "wideP_consurf3",
    "LogoPlotCustomH",
    "LogoPlotMSA"
  ),
  plot_df=c(
    "mutable_df3"  ,
    "lin_lf",
    "mutable_df3",
    "merged_df3" ,
    "lf_duet" ,
    "corr_df_m3_f",
    "merged_df2",
    "merged_df3",
    "merged_df2",
    "unified_msa"
  )
)

stability_boxes_df=data.frame(
  outcome_colname=c("duet_outcome",
                    "foldx_outcome",
                    "deepddg_outcome",
                    "ddg_dynamut2_outcome",
                    "mcsm_na_outcome",
                    "mcsm_ppi2_outcome",
                    "snap2_outcome",
                    "consurf_outcome",
                    "avg_stability_outcome"),
  stability_type=c(
    "DUET",
    "FoldX",
    "DeepDDG",
    "Dynamut2",
    "mCSM-NA",
    "mCSM-ppi2",
    "SNAP2",
    "Consurf",
    "Average"
  ),
  stability_colname=c(
    "duet_scaled",
    "foldx_scaled",
    "deepddg_scaled",
    "ddg_dynamut2_scaled",
    "mcsm_na_scaled",
    "mcsm_ppi2_scaled",
    "snap2_scaled",
    "consurf_scaled",
    "avg_stability_scaled"
  )
)

table_columns = c(
  "position",
  "mutationinformation",
  "sensitivity",
  "ligand_distance",
  "avg_lig_affinity",
  "avg_lig_affinity_outcome",
  "avg_stability",
  "avg_stability_outcome",
  "or_mychisq",
  "maf",
  "snap2_outcome",
  "consurf_outcome",
  "provean_outcome",
  "rsa",
  "kd_values" ,
  "rd_values"
)

logoPlotSchemes <- list("chemistry"
                        , "taylor"
                        , "hydrophobicity"
                        , "clustalx")
dm_om_methods = c("DUET ΔΔG"
                  , "Consurf"
                  , "Deepddg ΔΔG"
                  , "Dynamut2 ΔΔG"
                  , "FoldX ΔΔG"
                  , "Ligand affinity (log fold change)"
                  , "mCSM-NA affinity ΔΔG"
                  , "SNAP2")
#### target_map: handy gene/drug mapping hash ####
target_map = hash(
  c(
    "alr",
    "gid",
    "embb",
    "pnca",
    "rpob",
    "katg"),
  c(
    "cycloserine",
    "streptomycin",
    "ethambutol",
    "pyrazinamide",
    "rifampicin",
    "isoniazid")
)

#load E V E R Y T H I N G
lapply(c(
  "alr",
  "embb",
  "gid",
  "katg",
  "pnca",
  "rpob"
),function(x){
  invisible(load_target_globals(x))
  invisible(load_msa_global(x))
}
)

#### Shiny UI #####
if(interactive()){
ui <- dashboardPage(
  #dashboardHeader(title = paste0(gene, "/", drug)),
  dashboardHeader(title = "Sequence Alignment"),
  
  dashboardSidebar(
    sidebarMenu( id = "sidebar",
                 selectInput(
                   "switch_target",
                   label="Target",
                   choices = c(
                     "alr",
                     "embb",
                     "gid",
                     "katg",
                     "pnca",
                     "rpob"
                   ),
                   selected="embb"),
                 menuItem("LogoP ED", tabName="LogoP ED"),
                 
                 sliderInput(
                   "display_position_full_range"
                   , "Display Positions"
                   , min=1, max=150, value=c(1,150)
                 ),
                 
                 conditionalPanel(
                   condition="
                                  input.sidebar == 'LogoP SNP' ||
                                  input.sidebar == 'LogoP OR' ||
                                  input.sidebar == 'LogoP ED'", 
                   selectInput(
                     "logoplot_colour_scheme", 
                     label="Logo Plot Colour Scheme", 
                     choices = logoPlotSchemes, 
                     selected="chemistry"
                   )
                 )
    )
  ),
  body <- dashboardBody(
    
    tabItems(
      tabItem(tabName = "dashboard",
              h2("Dashboard tab content")
      ),
      
      tabItem(tabName = "widgets",
              h2("Widgets tab content")
      )
    ),
    # creates a 'Conditional Panel' containing a plot object from each of our 
    # ggplot plot functions (and its associated data frame)
    fluidRow(
      column(
        width=12,
        plotOutput("LogoPlotMSA",
                   click = "plot_click") %>% withSpinner(color="#0dc5c1")
      )
    )
  )
)
#### Shiny Backend Server #####
server <- function(input, output, session) {
  observeEvent(
    {
      input$display_position_full_range #special-purpose for MSA
      input$logoplot_colour_scheme
      input$switch_target
    },
    {
      # C O M P A T I B I L I T Y
      #gene=input$switch_target
      #drug=target_map[[gene]]
      merged_df3=cbind(get(paste0(input$switch_target, '_merged_df3')))
      
      position_max=max(merged_df3[['position']])
      position_min=min(merged_df3[['position']])
      min_ligand_distance=min(merged_df3$ligand_distance)
      max_ligand_distance=max(merged_df3$ligand_distance)
      # FIXME: these are IMPORTANT
      # # add "pos_count" position count column
      # merged_df3=merged_df3 %>% dplyr::add_count(position)
      # merged_df3$pos_count=merged_df3$n
      # merged_df3$n=NULL
      # 
      mutable_df3 = cbind(merged_df3)
      unified_msa = get(paste0(input$switch_target, '_unified_msa'))
      # 
      # # re-sort the dataframe according to position count
      sorted_df = cbind(merged_df3)
      sorted_df = sorted_df %>% arrange(pos_count)
      
      # 
      outdir = paste0("~/git/Data/", drug, '/output/')
      indir  = paste0("~/git/Data/", drug , "/input/")
      # 
      # source("~/git/LSHTM_analysis/scripts/plotting/logo_data_msa.R") # probably unnecessary...
      
      #  source("~/git/LSHTM_analysis/scripts/plotting/get_plotting_dfs.R")
      
      #### nasty special-purpose merged_df3 variants ####
      # FIXME: SLOW
      # corr_plotdf = corr_data_extract(
      #   merged_df3
      #   , gene = gene
      #   , drug = drug
      #   , extract_scaled_cols = F
      # )
      
      #input$stability_snp_param
      
      updateSliderInput(
        session, 
        "display_position_range",
        min = position_min,
        max = position_max
        #, value = c(position_min, position_min+150)
      )
      updateSliderInput(
        session,
        "display_position_full_range",
        min = 1, 
        max = position_max #, 
        # value = c(position_min, position_min+150)
      )
      updateNumericInput(session, "selected_logop_snp_position", min = position_min, max = position_max, value = position_min)
      updateNumericInput(session, "selected_logop_ed_position", min = position_min, max = position_max, value = position_min)
      updateNumericInput(session, "corr_lig_dist", min = min_ligand_distance, max = max_ligand_distance, value = min_ligand_distance)
      
      stability_colname      = stability_boxes_df[stability_boxes_df$stability_type==input$stability_snp_param,"stability_colname"]
      outcome_colname        = stability_boxes_df[stability_boxes_df$stability_type==input$stability_snp_param,"outcome_colname"]
      
      display_position_range = input$display_position_range
      plot_min=display_position_range[1]
      plot_max=display_position_range[2]
      
      display_position_full_range = input$display_position_full_range
      full_range_min=display_position_full_range[1]
      full_range_max=display_position_full_range[2]
      
      logoplot_colour_scheme = input$logoplot_colour_scheme
      omit_snp_count = input$omit_snp_count
      
      print(paste0('Plotting positions between: ', plot_min, ' and ', plot_max))
      
      subset_mutable_df3=mutable_df3[(mutable_df3$position>=plot_min & mutable_df3$position <=plot_max),]
      
      subset_mutable_df3=mutable_df3[(mutable_df3$position>=plot_min & mutable_df3$position <=plot_max),]
      subset_sorted_df=sorted_df[(sorted_df$position>=plot_min & sorted_df$position <=plot_max),]
      
      
      
      #### LogoPlotMSA/Logo Plot ED ####
      output$LogoPlotMSA = renderPlot(
        LogoPlotMSA(target=input$switch_target,
                    plot_positions=full_range_min:full_range_max,
                    my_logo_col = logoplot_colour_scheme,
                    aa_pos_drug = aa_pos_drug,
                    active_aa_pos = active_aa_pos,
                    aa_pos_lig1 = aa_pos_lig1,
                    aa_pos_lig2 = aa_pos_lig2,
                    aa_pos_lig3 = aa_pos_lig3
        )
      )
    }
  )
  
  #### EOF Shiny Server ####
  
  
  
}
################ Running Server ##############
app <- shinyApp(ui, server)


runApp(app)
}