From c12d2ea726c07c910e8d729b09ecd4b298a15e7b Mon Sep 17 00:00:00 2001
From: Marcel Costa <m.costag@idibell.cat>
Date: Fri, 5 Apr 2024 14:19:03 +0000
Subject: [PATCH] first commit

---
 Analysis.R      | 104 +++++++++++++++++++
 README.md       |   0
 functionsCyto.R | 267 ++++++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 371 insertions(+)
 create mode 100644 Analysis.R
 create mode 100644 README.md
 create mode 100644 functionsCyto.R

diff --git a/Analysis.R b/Analysis.R
new file mode 100644
index 0000000..c8775b6
--- /dev/null
+++ b/Analysis.R
@@ -0,0 +1,104 @@
+library(tidyverse)
+library(flowWorkspace)
+library(Biobase)
+library(flowGate)
+source("functionsCyto.R")
+
+files<-list.files("Files", pattern = ".LMD", full.names = T)
+
+LMD2FCS(files)
+
+fs<-read.ncdfFlowSet(files=list.files("Files",".fcs", full.names = T), readonly = F)
+gs <- GatingSet(fs)
+
+comp<-as.matrix(read.csv("CompMatrix.csv", check.names = F, row.names = 1))
+gs<-compensate(gs, comp)
+
+trans_params<-transform_gs(gs)
+saveRDS(trans_params, "trans_params.rds")
+trans_params<-readRDS("trans_params.rds")
+trans_apply(gs, trans_params = trans_params)
+
+
+colnames(gs)[c(4,6,11)]<-c("CD16","CD56","CD107a")
+sampleNames(gs)<-gsub("\\s[0-9]*.fcs","",sampleNames(gs))
+sampleNames(gs)<-gsub(".*\\s","",sampleNames(gs))
+pData(gs)$name<-rownames(pData(gs))
+
+gates<-list()
+gs_gate_interactive(gs,
+  filterId = "Leukocytes",
+  dims = list("FS-A", "SS-A"))
+gates[["Leukocytes"]]<-gs_pop_get_gate(gs, "Leukocytes")
+
+gs_gate_interactive(gs,
+  subset = "Leukocytes",
+  filterId = "NKdim",
+  dims = list("CD16", "CD56"))
+gates[["NKdim"]]<-gs_pop_get_gate(gs, "NKdim")
+
+gs_gate_interactive(gs[["NKs-K562"]],
+  subset = "NKdim",
+  filterId = "CD107a+",
+  dims = list("CD107a", "CD56"), overlayGates = "CD107a+")
+gates[["CD107a+"]]<-gs_pop_get_gate(gs, "CD107a+")
+
+gates<-gates_save(gs, file = "gates.rds")
+gates<-readRDS("gates.rds")
+gs<-gates_apply(gs, gates)
+
+autoplot(gs, "CD107a+")
+ggcyto_trans(gs, "CD107a", "CD56", subset="NKdim")
+autoplot(gs[[5]], bins=0)
+ggcyto_trans_all(gs, index = 5, ncol=1)
+
+stats<-gs_pop_get_stats(gs, nodes="CD107a+", type="perc")
+
+ggpubr::ggarrange(
+  ggcyto_trans_all(gs, nrow=1),
+  ggcyto::as.ggplot(ggcyto_trans(gs, "CD107a", "CD56", subset="NKdim")+
+    facet_grid(.~name)),
+  ggplot(stats, aes(sample, percent))+
+    geom_bar(stat="identity", color="black", fill="grey70"),
+  ncol=1)
+ggsave("Analysis.jpg")
+
+save_gs(gs, "gs_analysis")
+gs<-load_gs("gs_analysis/")
+
+## For a new experiment with same settings
+
+files<-list.files("Files", pattern = ".LMD", full.names = T)
+LMD2FCS(files)
+
+fs<-read.ncdfFlowSet(files=list.files("Files",".fcs", full.names = T), readonly = F)
+gs <- GatingSet(fs)
+
+comp<-as.matrix(read.csv("CompMatrix.csv", check.names = F, row.names = 1))
+gs<-compensate(gs, comp)
+
+trans_params<-readRDS("trans_params.rds")
+trans_apply(gs, trans_params = trans_params)
+
+colnames(gs)[c(4,6,11)]<-c("CD16","CD56","CD107a")
+sampleNames(gs)<-gsub("\\s[0-9]*.fcs","",sampleNames(gs))
+sampleNames(gs)<-gsub(".*\\s","",sampleNames(gs))
+pData(gs)$name<-rownames(pData(gs))
+
+gates<-readRDS("gates.rds")
+gs<-gates_apply(gs, gates) #if not the same gates, apply a general with "exact=F"
+
+
+stats<-gs_pop_get_stats(gs, nodes="CD107a+", type="perc")
+
+ggpubr::ggarrange(
+  ggcyto_trans_all(gs, nrow=1),
+  ggcyto::as.ggplot(ggcyto_trans(gs, "CD107a", "CD56", subset="NKdim")+
+                      facet_grid(.~name)),
+  ggplot(stats, aes(sample, percent))+
+    geom_bar(stat="identity", color="black", fill="grey70"),
+  ncol=1)
+ggsave("Analysis.jpg")
+
+save_gs(gs, "gs_analysis")
+gs<-load_gs("gs_analysis/")
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..e69de29
diff --git a/functionsCyto.R b/functionsCyto.R
new file mode 100644
index 0000000..40954db
--- /dev/null
+++ b/functionsCyto.R
@@ -0,0 +1,267 @@
+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=5, 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){
+  gates.list<-list()
+  pop_paths<-gs_get_pop_paths(gs)[gs_get_pop_paths(gs) != "root"]
+  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)
+}