library(flowClust)
library(flowWorkspace)
library(openCyto)

load the gating set

gs <- load_gs("lyoplate_out/gated_data/auto/gslist-tcell/D54tFo7RPl/")

## loading R object...
## loading tree object...
## Done

grab CD3+ tcell population from one sample

parent <- getData(gs[[1]], "CD3")

## Error in getData(gs[[1]], "CD3"): unused argument ("CD3")

plot it on cd8 vs cd4 projection

xyplot(`CD8`~`CD4`, parent, smooth = FALSE, xbin = 32)

Using flowClust with priors (openCyto::quadGate.tmix is a wrapper function based on flowClust)

prior=list(Mu0 = matrix(c(1000,3500
                          ,1000,1000
                          ,2800,2000
                          )
                        , nrow = 3
                        , ncol = 2
                        ,byrow=T
                        )
           ,w0=c(10,10,10)
           ,nu0 = c(4,4,4)
           ,Lambda0=array(c(rep(1e5,3)
                            ,rep(0,3)
                            ,rep(0,3)
                            ,rep(1e5,3)
                            )
                          ,dim=c(3,2,2))
           ,Omega0=array(c(rep(5e4,3)
                           ,rep(0,3)
                           ,rep(0,3)
                           ,rep(5e4,3)
                            )
                         ,dim=c(3,2,2)
                         )
           )

gates <- quadGate.tmix(parent, c("CD4", "CD8"), K=3, prior = prior, quantile3=0.8,quantile1=0.8)

## Using the parallel (multicore) version of flowClust with 16 cores

xyplot(`CD8`~`CD4`, parent, smooth = FALSE, xbin = 32, filter = gates, stats = TRUE)

Using flowClust without priors (tmixGate is a wrapper function based on flowClust, it constructs flowCore gates from tmixFilterResult)

for(K in 2:4){
#   res <- flowClust(x = parent, varNames = c("CD4", "CD8"), K = k)
#   plot(res, data = parent, main = paste0("K = ", k), xlim = c(-1000, 4600), ylim =c(-1000, 4600))
    gates <- tmixGate(fr = parent, channels = c("CD4", "CD8"), K = K)  
    print(xyplot(`CD8`~`CD4`, parent, smooth = FALSE, xbin = 32, filter = gates, stats = TRUE, main = paste0("K = ", K)))
}

## The prior specification has no effect when usePrior=no
## Using the parallel (multicore) version of flowClust with 16 cores
## The prior specification has no effect when usePrior=no
## Using the parallel (multicore) version of flowClust with 16 cores

## The prior specification has no effect when usePrior=no
## Using the parallel (multicore) version of flowClust with 16 cores

Using flowMerge

library(flowMerge)
flowClust.res <- flowClust(parent, varNames=c("CD4", "CD8"), K = 1: 8,trans = 0, nu.est=1,randomStart=20);

## Using the parallel (multicore) version of flowClust with 16 cores

plot(flowMerge:::BIC(flowClust.res), main="BIC for 1 through 8 cluster flowClust solutions",xlab="K",ylab="BIC",type="o");

flowClust.maxBIC<-flowClust.res[[which.max(BIC(flowClust.res))]];
flowClust.flowobj<-flowObj(flowClust.maxBIC,parent);
flowClust.merge<-merge(flowClust.flowobj,metric="entropy");

## Merged to 6 clusters
## Merged to 5 clusters
## Merged to 4 clusters
## Merged to 3 clusters
## Merged to 2 clusters
## Merged to 1 clusters
## Updating model statistics

## Rule of identifying outliers: 90% quantile

## Updated model 1

## Rule of identifying outliers: 90% quantile

## Updated model 2

## Rule of identifying outliers: 90% quantile

## Updated model 3

## Rule of identifying outliers: 90% quantile

## Updated model 4

## Rule of identifying outliers: 90% quantile

## Updated model 5

## Rule of identifying outliers: 90% quantile

## Updated model 6

## Rule of identifying outliers: 90% quantile

## Updated model 7

i<-fitPiecewiseLinreg(flowClust.merge)
par(mfrow=c(2,2));
flowClust.mergeopt<-flowClust.merge[[i]];
plot(flowClust.res[[4]],data=parent,main="Max BIC solution", xlim = c(-1000, 4600), ylim =c(-1000, 4600));

## Rule of identifying outliers: 90% quantile

plot(flowClust.res[[which.max(flowMerge:::ICL(flowClust.res))]],data=parent,main="Max ICL solution", xlim = c(-1000, 4600), ylim =c(-1000, 4600));

## Rule of identifying outliers: 90% quantile

plot(flowClust.mergeopt,level=0.75,pch=20,main="Merged Solution", xlim = c(-1000, 4600), ylim =c(-1000, 4600))

## Rule of identifying outliers: 75% quantile

# pop <- which(apply(apply(getEstimates(flowClust.mergeopt)$locations,2,function(x)order(x,decreasing=TRUE)==3),1,all));
# lymphocytes<-split(flowClust.mergeopt,population=list("lymphocytes"=pop))$lymphocytes;
# lymphocytes<-lymphocytes[,c(3,5)];
# l.flowC<-flowClust(lymphocytes,varNames=c("FL1.H","FL3.H"),K=1:8,B=1000,B.init=100,tol=1e-5)
# 
# par(mfrow=c(2,2));
# l.flowO<-flowObj(l.flowC[[which.max(flowMerge:::BIC(l.flowC))]],lymphocytes);
# plot(l.flowO,main="max BIC solution",new.window=F,pch=20,level=0.9);
# 
# plot(flowObj(l.flowC[[which.max(flowMerge:::ICL(l.flowC))]],lymphocytes),main="max ICL solution",new.window=F,pch=20,level=0.9);
#