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"))
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/")