CytoR/functionsCyto.R

gs_export_pop<-function(gs, pop, trim.channels=T){
  pop.dt<-c()
  # Initializes the progress bar
  pb <- txtProgressBar(min = 0,      # Minimum value of the progress bar
                       max = length(gs), # Maximum value of the progress bar
                       style = 3,    # Progress bar style (also available style = 1 and style = 2)
                       width = 50,   # Progress bar width. Defaults to getOption("width")
                       char = "=")   # Character used to create the bar
  
  
  for (i in 1:length(gs)){
    samp<-sampleNames(gs)[i]
    data<-as.data.frame(flowCore::exprs(gh_pop_get_data(gs[[samp]], pop))) %>% add_column("Sample"=samp)
    if (length(pop.dt) > 0){
      pop.dt<-rbind(pop.dt, data)
    }else{
      pop.dt<-data
    }
    setTxtProgressBar(pb, i)
  }
  close(pb)
  if (trim.channels){
    colnames(pop.dt)<-gsub(" [A-Za-z0-9.-]*","",colnames(pop.dt))
    pop.dt<-pop.dt[colnames(pop.dt)[!colnames(pop.dt) %in% c("FSC-H","SSC-H","FSC-Width","FVS")]]
  }
  return(pop.dt)
}

get_cut<-function(x, markers, ref, thr){
  x<-x %>%
    dplyr::filter(Condition == ref) %>%
    gather(Marker, Value, all_of(markers)) %>% 
    group_by(`ID.SAR-EM`, Marker, Tiempo) %>% 
    summarise(cut=predict_cut(Value, thr))
}

get_positive<-function(x, cut=NULL, markers, ref, thr){
  
  if (is.null(cut)){
    cut<-get_cut(cd8, markers = markers, ref = ref, thr = thr)
  }
  
  x<-merge(
    x %>% gather(Marker, Value, all_of(unique(cut$Marker))),
    cut
  )
  
  x<-x %>%
    mutate(Marker_cut=Value >= cut)
  
  x<-x %>%
    select(`ID.SAR-EM`, Tiempo, Condition, CellID, Marker, Marker_cut) %>%
    mutate(Marker_pos=paste0(Marker,".",Marker_cut)) %>%
    group_by(`ID.SAR-EM`,Tiempo, Condition, CellID) %>% 
    summarise(Pheno=paste0(Marker_pos, collapse = " ")) %>%
    # mutate(Pheno=paste0("IFNg.",IFNg_p," ","TNFa.",TNFa_p," ","CD69.",CD69_p," ","CD107a.",CD107a_p)) %>% 
    mutate(Pheno=gsub("[A-Za-z0-9]*.FALSE |[A-Za-z0-9]*.FALSE$| [A-Za-z0-9]*.FALSE","",Pheno),
           Pheno=gsub(".TRUE","+",Pheno),
           Pheno=gsub("^$","Neg",Pheno))
  
  
  x.dt<-x %>% 
    group_by(`ID.SAR-EM`, Tiempo, Condition, Pheno) %>% 
    count() %>%
    group_by(`ID.SAR-EM`, Tiempo, Condition) %>%
    mutate(perc=n*100/sum(n))
  return(x.dt)
}

transform_gs<-function(gs, saveparams=NULL, trans_params_obj=NULL, index=1){
  
  if (is.null(trans_params_obj)){
    log_mark<-colnames(gs)[!grepl("FS|SS|Time", colnames(gs))]
    trans_params<-lapply(log_mark, function(x) list(channel=x, scale="biexp",maxvalue=250000,pos=4, widthBasis=-400, max=10^3.7, min=0))
    names(trans_params)<-log_mark
    log_mark_lin<-colnames(gs)[grepl("FS|SS|Time", colnames(gs))]
    trans_params_lin<-lapply(log_mark_lin, function(x) list(channel=x, scale="lin", 
                                                            max=max(flowCore::exprs(gh_pop_get_data(gs, "root"))[,x])*1.1, min=0))
    names(trans_params_lin)<-log_mark_lin
    trans_params<-c(trans_params_lin, trans_params)
  }else{
    log_mark<-colnames(gs)[!grepl("FSC|SSC|Time", colnames(gs))]
    log_mark_lin<-colnames(gs)[grepl("FSC|SSC|Time", colnames(gs))]
  }
  
  shiny::runApp(shiny::shinyApp(
    ui = fluidPage(
      sidebarLayout(
        sidebarPanel(selectInput("channel", "Channel", colnames(gs), multiple = F),
                     uiOutput("maxvalue"),
                     uiOutput("pos"),
                     uiOutput("widthBasis"),
                     uiOutput("max"),
                     uiOutput("min"),
                     actionButton("applyBut", "Apply"),
                     actionButton("stopBut", "Quit")),
        mainPanel(plotOutput("hist"))
      )
    ),
    server = function(input, output) {
      output$hist <- renderPlot({
        gs.trans<-gs_clone(gs[index])
        max<-input$max
        if (trans_params[[input$channel]]$scale == "biexp"){
          trans.obj<-flowjo_biexp_trans(maxValue = input$maxvalue, pos=input$pos, 
                                        widthBasis=input$widthBasis)
          transList <- transformerList(input$channel, trans.obj)
          gs.trans <- transform(gs.trans, transList)
          max<-10^input$max
        }
        
        ggcyto(gs.trans, aes(x = .data[[input$channel]]), subset = "root")+
          geom_density(fill = "blue", alpha= 0.5)+
          ggcyto_par_set(limits = list(x = c(input$min,max)))
      })
      output$maxvalue<-renderUI({if(trans_params[[input$channel]]$scale == "biexp"){
        sliderInput("maxvalue", "maxValue", min=0, max=1000000, value=trans_params[[input$channel]]$maxvalue)}})
      output$pos<-renderUI({if(trans_params[[input$channel]]$scale == "biexp"){
        sliderInput("pos", "pos", min=0, max=10, value=trans_params[[input$channel]]$pos, step = 0.5)}})
      output$widthBasis<-renderUI({if(trans_params[[input$channel]]$scale == "biexp"){
        sliderInput("widthBasis", "widthBasis", min=-1000, max=0, value=trans_params[[input$channel]]$widthBasis)}})
      output$max<-renderUI({
        if(trans_params[[input$channel]]$scale == "biexp"){
          sliderInput("max", "Max X", min=0, max=6, value=log10(trans_params[[input$channel]]$max), step=0.1)
        }else{
          sliderInput("max", "Max X", min=0, max=trans_params[[input$channel]]$max, 
                      value=trans_params[[input$channel]]$max)
        }
      })
      output$min<-renderUI({
        if(trans_params[[input$channel]]$scale == "biexp"){
          sliderInput("min", "Min X", min=-1000, max=100000, value=trans_params[[input$channel]]$min)
        }else{
          sliderInput("min", "Min X", min=0, max=trans_params[[input$channel]]$max, value=trans_params[[input$channel]]$min)
        }
      })
      
      observeEvent(input$applyBut, {
        print("Transformation and Scaling Updated!")
        if (isolate(input$channel) %in% log_mark){
          trans_params[[isolate(input$channel)]]$maxvalue<<-isolate(input$maxvalue)
          trans_params[[isolate(input$channel)]]$pos<<-isolate(input$pos)
          trans_params[[isolate(input$channel)]]$widthBasis<<-isolate(input$widthBasis)
          trans_params[[isolate(input$channel)]]$max<<-10^isolate(input$max)
          trans_params[[isolate(input$channel)]]$min<<-isolate(input$min)
        }
        if (isolate(input$channel) %in% log_mark_lin){
          trans_params[[isolate(input$channel)]]$max<<-isolate(input$max)
          trans_params[[isolate(input$channel)]]$min<<-isolate(input$min)
        }
      })
      observeEvent(input$stopBut, {
        stopApp()
      })
    }
  ))
  
  trans_apply(gs, trans_params)
  
  print("Transformation Applied")
  if (!is.null(saveparams)){
    saveRDS(trans_params, file=saveparams)
  }
  return(trans_params)
}

trans_apply<-function(gs, trans_params){
  for (i in 1:length(trans_params)){
    if (trans_params[[i]]$scale == "biexp"){
      print(trans_params[[i]]$channel)
      trans.obj<-flowjo_biexp_trans(maxValue = trans_params[[i]]$maxvalue, 
                                    pos=trans_params[[i]]$pos, 
                                    widthBasis=trans_params[[i]]$widthBasis)
      transList <- transformerList(trans_params[[i]]$channel, trans.obj)
      gs <- transform(gs, transList) 
    }
  }
}

ggcyto_trans<-function(gs, x, y, trans_params_obj=trans_params, subset, gates=T, stats=T, ...){
  g<-ggcyto(gs, aes(.data[[x]], .data[[y]]), subset=subset, ...)+
    ggcyto_par_set(limits = list(x = c(trans_params_obj[[x]]$min,trans_params_obj[[x]]$max), 
                                 y=c(trans_params_obj[[y]]$min,trans_params_obj[[y]]$max)))
  
  if(gates==T){
    if (subset == "root"){
      gates<-gs_get_pop_paths(gs)[grepl(paste0("^/","[^/]*$"), gs_get_pop_paths(gs))]
    }else{
      gates<-gs_get_pop_paths(gs)[grepl(paste0(subset,"/","[^/]*$"), gs_get_pop_paths(gs))]  
    }
    g<-g+geom_point(size=0.1, alpha=0.1)+
      geom_gate(gates)
    if (stats== T){
      g<-g+geom_stats(gates)
    }
  }
  return(g)
}

ggcyto_trans_all<-function(gs, index=1, trans_params_obj=trans_params, stats=T,...){
  subsets<-unlist(sapply(gs_get_pop_paths(gs), function(x) strsplit(x, "/")[[1]][length(strsplit(x, "/")[[1]])-1]))
  subsets<-c("root",unique(subsets[subsets != ""]))
  
  g.list<-list()
  
  for (subset in subsets){
    if (subset == "root"){
      gates<-gs_get_pop_paths(gs)[grepl(paste0("^/","[^/]*$"), gs_get_pop_paths(gs))]
    }else{
      gates<-gs_get_pop_paths(gs)[grepl(paste0(subset,"/","[^/]*$"), gs_get_pop_paths(gs))]  
    }
    subset.x<-names(gs_pop_get_gate(gs, gates[1])[[1]]@parameters)[1]
    subset.y<-names(gs_pop_get_gate(gs, gates[1])[[1]]@parameters)[2]
    g.list[[subset]]<-ggcyto::as.ggplot(ggcyto_trans(gs[[index]], subset.x, subset.y, 
                                                     trans_params_obj, subset, stats=stats))
  }
  g<-do.call(ggpubr::ggarrange, c(g.list, ...))
  g<-ggpubr::annotate_figure(g, top = ggpubr::text_grob(sampleNames(gs)[index], 
                                                        color = "black", face = "bold", size = 14))
  return(g)
}

LMD2FCS<-function(files, output.dir=NULL){
  if(is.null(output.dir)){
    route<-paste0(gsub("/[^/]*$","",files[1]),"/")
  }
  for (lmd in files){
    lmd<-gsub(".*/","",lmd)
    fcs<-read.FCS(paste0(route,lmd), dataset = 2)
    keyword(fcs)['$FIL']<-paste0(gsub(".LMD","",lmd), ".fcs")
    write.FCS(fcs, paste0(route, gsub(".LMD","",lmd), ".fcs"))
  }
  print("Conversión completada")
}

gates_save<-function(gs, file="gates.rds", save=T, include=NULL){
  gates.list<-list()
  pop_paths<-gs_get_pop_paths(gs)[gs_get_pop_paths(gs) != "root"]
  if(!is.null(include)){pop_paths<-pop_paths[pop_paths %in% include]}
  pop_paths<-lapply(pop_paths, function(x) rev(rev(strsplit(x, "/")[[1]])[1:2]))
  for (gate in pop_paths){
    parent<-if (gate[1] == ""){"root"}else{gate[1]}
    gated_pop<-gate[2]
    gates.list[[paste0(parent,"/",gated_pop)]]<-gs_pop_get_gate(gs, gated_pop)
  }
  if(save){saveRDS(gates.list, file)}
  return(gates.list)
}

gates_apply<-function(gs, gates, exact=T){
  if (exact){
    for (gate in names(gates)){
      parent<-strsplit(gate, "/")[[1]][1]
      gated_pop<-strsplit(gate, "/")[[1]][2]
      gs_pop_add(gs, gates[[gate]], parent=parent, name=gated_pop)
      print(paste0("Gate ",gated_pop," applied!"))
    }
  }else{
    for (gate in names(gates)){
      parent<-strsplit(gate, "/")[[1]][1]
      gated_pop<-strsplit(gate, "/")[[1]][2]
      gs_pop_add(gs, gates[[gate]][[1]], parent=parent, name=gated_pop)
      print(paste0("Gate ",gated_pop," applied!"))
    }
  }
  recompute(gs)
  return(gs)
}

gs_pop_get_children_recursive<-function(gs, pop){
  childrens<-c()
  n<-0
  childrens<-c(childrens, gs_pop_get_children(gs, pop))
  while(n != length(childrens)){
    n<-length(childrens)
    for (i in 1:length(childrens)){
      childrens<-c(childrens, gs_pop_get_children(gs, childrens[i]))
    }
    childrens<-unique(childrens)
  }
  return(childrens)
}

gs_gate_interactive_regate2<-function(gs, filterId, sample=1, subset="root", ...){
  dims<-list(names(gs_pop_get_gate(gs, filterId)[[1]]@parameters)[1],
              names(gs_pop_get_gate(gs, filterId)[[1]]@parameters)[2])
  childrens<-gs_pop_get_children_recursive(gs, filterId)
  childrens_gate<-gates_save(gs, save=F, include=childrens)
  
  gs_pop_add(gs, gs_pop_get_gate(gs, filterId), parent=subset, name="duplicated")
  
  gs_gate_interactive(gs,
                      subset = subset,
                      filterId = filterId,
                      sample=sample,
                      dims = dims, regate=T, overlayGates = "duplicated")
  
  gs_pop_remove(gs, "duplicated")

  gates_apply(gs, childrens_gate)
}