troubleshoot_biexp.R

library(ggcyto)

## Loading required package: ggplot2

## Loading required package: flowCore

## Loading required package: ncdfFlow

## Loading required package: RcppArmadillo

## Loading required package: BH

## Loading required package: flowWorkspace

# #-----------
# #parse diva workspace
# path <- "/loc/no-backup/mike/shared/diva/2001-D-g001/"
# ws <- openDiva(file.path(path, "/2001-D-g001.xml"))
# gs <- parseWorkspace(ws, name = 2, subset = 6)
# save_gs(gs, file.path(path, "gs"))
# autoplot(gs[[1]])
# autoplot(gs, "P2") + ggcyto_par_set(limits = "data")
#
# g <- getGate(gs[[1]], "P2")
# g
# gs <- load_gs(file.path(path, "gs"))
# ### extract and save the data of interest
# fr <- getData(gs, "L")[[1, c("FSC-A", "G660-A")]]
# saveRDS(fr, file = "/fh/fast/gottardo_r/mike_working/diva_parse_troubleshoot/fr.rds")
# trans <- getTransformations(gs[[1]], only = F)[["G660-A"]]
# saveRDS(trans, file = "/fh/fast/gottardo_r/mike_working/diva_parse_troubleshoot/trans.rds")
# saveRDS(g, "/fh/fast/gottardo_r/mike_working/diva_parse_troubleshoot/gate.rds")


#load the archived the data
fr_saved <- readRDS("/fh/fast/gottardo_r/mike_working/diva_parse_troubleshoot/fr.rds")
g <- readRDS("/fh/fast/gottardo_r/mike_working/diva_parse_troubleshoot/gate.rds")
trans_saved <- readRDS("/fh/fast/gottardo_r/mike_working/diva_parse_troubleshoot/trans.rds")
##subsample data to speed up plot
set.seed(1)
fr <- fr_saved[sample.int(nrow(fr_saved), nrow(fr_saved) *0.01),]

#now we try different ways of constructing biexp to see if we can get closer to DIVA results

##inverse transform the data to restore the raw data
exprs(fr)[, "G660-A"] <- trans_saved$inverse(exprs(fr)[, "G660-A"])

##copy the original gate coordinates from diva xml
g@min[[2]] <- 1.4
g@max[[2]] <- 2.86

##try to fiddle with various parameters recorded in xml
xml.min <- 0
xml.max <- 5.418
xml.biexp_scale <- 805
xml.comp_biexp_scale <- 10941
xml.manual_biexp_scale <- 8392

#use the equation suggested by BD engineer last year
generate_trans <- function(maxValue, pos, r)
{
  if(r == 0)
    r <- maxValue/10^pos
  w <- (pos - log10(maxValue/r))/2
  if(w < 0)
    w <- 0
  logicleTransform(w=w, t = maxValue, m = pos) #
}

#construct the new trans
trans <- generate_trans(maxValue = 262144, pos = 4.5, r = xml.comp_biexp_scale)
as.list(environment(trans))[3:6]

## $w
## [1] 1.560259
## 
## $t
## [1] 262144
## 
## $m
## [1] 4.5
## 
## $a
## [1] 0

#transform the gate
g@min[[2]] <- trans(10^1.4) #inverse log10 coordinates to raw scale before transform it to logicle
g@max[[2]] <- trans(10^2.86)
g

## Rectangular gate 'P2' with dimensions:
##   FSC-A: (74770.453125,209924.421875)
##   G660-A: (1.56461719640015,1.68571960007558)

#transform the data
exprs(fr)[, "G660-A"] <- trans(exprs(fr)[, "G660-A"])
#show results
autoplot(fr, "FSC-A", "G660-A") + geom_gate(g) + geom_stats()

## Warning: Removed 5 rows containing missing values (geom_hex).

#didn't work, try different one
set.seed(1)
fr <- fr_saved[sample.int(nrow(fr_saved), nrow(fr_saved) *0.01),]
exprs(fr)[, "G660-A"] <- trans_saved$inverse(exprs(fr)[, "G660-A"])

trans <- generate_trans(maxValue = 262144, pos = 4.5, r = xml.manual_biexp_scale)
as.list(environment(trans))[3:6]

## $w
## [1] 1.502663
## 
## $t
## [1] 262144
## 
## $m
## [1] 4.5
## 
## $a
## [1] 0

g@min[[2]] <- trans(10^1.4) #inverse log10 coordinates to raw scale before transform it to logicle
g@max[[2]] <- trans(10^2.86)
g

## Rectangular gate 'P2' with dimensions:
##   FSC-A: (74770.453125,209924.421875)
##   G660-A: (1.50787710470241,1.6525980706805)

exprs(fr)[, "G660-A"] <- trans(exprs(fr)[, "G660-A"])
autoplot(fr, "FSC-A", "G660-A") + geom_gate(g) + geom_stats()

## Warning: Removed 5 rows containing missing values (geom_hex).

##nop! try something new
trans <- generate_trans(maxValue = 10^xml.max, pos = xml.max, r = xml.comp_biexp_scale)
as.list(environment(trans))[3:6]

## $w
## [1] 2.019529
## 
## $t
## [1] 261818.3
## 
## $m
## [1] 5.418
## 
## $a
## [1] 0

g@min[[2]] <- trans(10^1.4) #inverse log10 coordinates to raw scale before transform it to logicle
g@max[[2]] <- trans(10^2.86)
g

## Rectangular gate 'P2' with dimensions:
##   FSC-A: (74770.453125,209924.421875)
##   G660-A: (2.02756293507503,2.25024613264607)

set.seed(1)
fr <- fr_saved[sample.int(nrow(fr_saved), nrow(fr_saved) *0.01),]
exprs(fr)[, "G660-A"] <- trans_saved$inverse(exprs(fr)[, "G660-A"])

exprs(fr)[, "G660-A"] <- trans(exprs(fr)[, "G660-A"])
autoplot(fr, "FSC-A", "G660-A") + geom_gate(g) + geom_stats()

## Warning: Removed 6 rows containing missing values (geom_hex).

##try to manually fill the parameters for logicle
# trans <- logicleTransform(w = 2, t = 262144, m = 6.4) #
# g@min[[2]] <- trans(10^1.4) #inverse log10 coordinates to raw scale before transform it to logicle
# g@max[[2]] <- trans(10^2.86)
# g
# set.seed(1)
# fr <- fr_saved[sample.int(nrow(fr_saved), nrow(fr_saved) *0.01),]
# exprs(fr)[, "G660-A"] <- trans_saved$inverse(exprs(fr)[, "G660-A"])
#
# exprs(fr)[, "G660-A"] <- trans(exprs(fr)[, "G660-A"])
# autoplot(fr, "FSC-A", "G660-A") + geom_gate(g) + geom_stats(digits = 4)