R Notebook

####################################################################
i = "ID19_16_048_ACC2583_Tube_013.fcs"
i= "ID18_02_104_ACC2673_Tube_013.fcs"
file_full = paste0("D://TCyto/",i)
FCS = read.FCS(file_full)
colnames(FCS)[which(colnames(FCS)=="SSC-A")] <- "SSCA"
colnames(FCS)[which(colnames(FCS)=="FSC-A")] <- "FSCA"
colnames(FCS)[which(colnames(FCS)=="Pacific Orange-A")] <- "CD8"
colnames(FCS)[which(colnames(FCS)=="Pacific Blue-A")] <- "CD4"
colnames(FCS)[which(colnames(FCS)=="PE-Texas Red-A")] <- "CD45RA"
colnames(FCS)[which(colnames(FCS)== "APC-Cy7-A")] <- "CD27"
colnames(FCS)[which(colnames(FCS)=="PE-Cy7-A")] <- "CCR7"
colnames(FCS)[which(colnames(FCS)=="PE-Cy5-A")] <- "CD3"
colnames(FCS)[which(colnames(FCS)=="PE-A")] <- "GranzA"
colnames(FCS)[which(colnames(FCS)=="FITC-A")] <- "Perf"
colnames(FCS)[which(colnames(FCS)=="APC-A")] <- "CD28"
colnames(FCS)[which(colnames(FCS)=="Alexa Fluor 700-A")] <- "GranB"
library(ggcyto)
biexp  <- biexponentialTransform("myTransform",w=10)
myTrans <-  transformList(c('CD8','CD4',"SSCA","FSCA","CD3","Perf",
                            "CD28","CD27","GranzA","CD45RA","CCR7","GranB"), biexp)
FCS1 = transform(FCS, myTrans)
# write FCS file in Flow 
outFile <- file.path("D://FCS2.fcs")
## now write out into a file
write.FCS(FCS1, outFile)

[1] "D://FCS2.fcs"

library(ncdfFlow)
file_full1 = "D://FCS2.fcs"
fs1<-read.flowSet(file_full1)
p4 <- ggcyto(fs1, aes(x =CD4, y = CD8)) + geom_hex(bins = 128)
rect.g <- rectangleGate(list("CD4" =  c(17,20), "CD8" = c(12,20)))
rect.gates <- sapply(sampleNames(fs1), function(sn)rect.g)
p4 + geom_gate(rect.gates) + geom_stats() + annotate(geom = "text", x = 18, y = 14, label  = "CD4+", color = "red")

T2 <- (FCS1@exprs)
M = T2[,c("SSCA","FSCA","CD3","CD4","CD8")]
M1 = T2[,c("SSCA","FSCA","CD4","CD8")]
cellevents <- M1
b_tree = CytomeTree(M1, force_first_markers = c("FSCA","SSCA","CD4","CD8"), minleaf = 20)

  |                                                                                        
  |                                                                                  |   0%mc_uni_df 2 mc_uni_loglik  -181377.3 
mc_mix_df 4 mc_mix_loglik  -114064.8 
n 224129 
mc_uni_df 2 mc_uni_loglik  -181377.3 
mc_mix_df 4 mc_mix_loglik  -114064.8 
n 224129 

  |                                                                                        
  |====================                                                              |  25%mc_uni_df 2 mc_uni_loglik  -1365.216 
mc_mix_df 4 mc_mix_loglik  33.64744 
n 57038 
mc_uni_df 2 mc_uni_loglik  -1365.216 
mc_mix_df 4 mc_mix_loglik  33.64744 
n 57038 
mc_uni_df 2 mc_uni_loglik  12030.84 
mc_mix_df 4 mc_mix_loglik  17767.38 
n 167091 
mc_uni_df 2 mc_uni_loglik  12030.84 
mc_mix_df 4 mc_mix_loglik  17767.38 
n 167091 

  |                                                                                        
  |=========================================                                         |  50%mc_uni_df 2 mc_uni_loglik  -4238.662 
mc_mix_df 4 mc_mix_loglik  -3701.914 
n 1566 
mc_uni_df 2 mc_uni_loglik  -4238.662 
mc_mix_df 4 mc_mix_loglik  -3701.914 
n 1566 
mc_uni_df 2 mc_uni_loglik  -130580.7 
mc_mix_df 4 mc_mix_loglik  -115339.8 
n 55472 
mc_uni_df 2 mc_uni_loglik  -130580.7 
mc_mix_df 4 mc_mix_loglik  -115339.8 
n 55472 
mc_uni_df 2 mc_uni_loglik  -21312.24 
mc_mix_df 4 mc_mix_loglik  -21312.85 
n 20717 

Unable to force split on CD4 for some node at level3

mc_uni_df 2 mc_uni_loglik  -105660.2 
mc_mix_df 4 mc_mix_loglik  -105666.3 
n 146374 
mc_uni_df 2 mc_uni_loglik  -105660.2 
mc_mix_df 4 mc_mix_loglik  -105666.3 
n 146374 

  |                                                                                        
  |==============================================================                    |  75%mc_uni_df 2 mc_uni_loglik  -748.2666 
mc_mix_df 4 mc_mix_loglik  -519.4429 
n 319 
mc_uni_df 2 mc_uni_loglik  -748.2666 
mc_mix_df 4 mc_mix_loglik  -519.4429 
n 319 
mc_uni_df 2 mc_uni_loglik  -2472.853 
mc_mix_df 4 mc_mix_loglik  -1918.037 
n 1247 
mc_uni_df 2 mc_uni_loglik  -2472.853 
mc_mix_df 4 mc_mix_loglik  -1918.037 
n 1247 
mc_uni_df 2 mc_uni_loglik  -6441.649 
mc_mix_df 4 mc_mix_loglik  -6441.541 
n 3203 
mc_uni_df 2 mc_uni_loglik  -6441.649 
mc_mix_df 4 mc_mix_loglik  -6441.541 
n 3203 
mc_uni_df 2 mc_uni_loglik  -88231.23 
mc_mix_df 4 mc_mix_loglik  -82805 
n 52269 
mc_uni_df 2 mc_uni_loglik  -88231.23 
mc_mix_df 4 mc_mix_loglik  -82805 
n 52269 
mc_uni_df 2 mc_uni_loglik  -52.05245 
mc_mix_df 4 mc_mix_loglik  -52.05768 
n 115 
mc_uni_df 2 mc_uni_loglik  -52.05245 
mc_mix_df 4 mc_mix_loglik  -52.05768 
n 115 
mc_uni_df 2 mc_uni_loglik  -95787.72 
mc_mix_df 4 mc_mix_loglik  -95793.47 
n 146259 
mc_uni_df 2 mc_uni_loglik  -95787.72 
mc_mix_df 4 mc_mix_loglik  -95793.47 
n 146259 

  |                                                                                        
  |==================================================================================| 100%

## ------------------------------------------------------------------------
# getting the only the cell events with the 3 markers measured:
#load("C:/Users/tranh/Documents/GitHub/Best.Rdata")
Tree_Partition <- b_tree$labels
Annot = Annotation(b_tree,K2markers = c("SSCA","FSCA","CD4","CD8"))

#Annot = Annbtree_annot
phenotypes <- list()
AA <- rbind(c("SSCA", 0), 
                 c("FSCA", 0))
                 
phenotypes[[1]] <- AA
BB <- rbind(c("SSCA", 1), 
                  c("FSCA", 1))
                  
phenotypes[[2]] <- BB
AB <- rbind(c("SSCA", 0), 
                 c("FSCA", 1))
                 
phenotypes[[3]] <- AB
BA <- rbind(c("SSCA", 1), 
                        c("FSCA", 0))
                        
            
phenotypes[[4]] <- BA 
methodd <- Populations <- X <- Y <- c()
#PhenoInfos <- RetrievePops(Annot, phenotypes)
PhenoInfos <- RetrievePops_p(b_tree, phenotypes)
MergedPops <- PhenoInfos$phenotypesinfo
AA <- MergedPops[[1]][length(MergedPops[[1]])]
BB <- MergedPops[[2]][length(MergedPops[[2]])]
AB <- MergedPops[[3]][length(MergedPops[[3]])]
BA <- MergedPops[[4]][length(MergedPops[[4]])]
MergedPartition <- PhenoInfos$Mergedleaves
ind_ <- MergedPartition%in%c(AA,BB,AB,BA) 
label2 <- MergedPartition[ind_]
LABEL2 <- rep(NA, length(label2))
LABEL2[label2 ==  AA] <- "SSCA-/FSCA-"
LABEL2[label2 ==  BB] <- "SSCA+/FSCA+"
LABEL2[label2 ==  AB] <- "SSCA-/FSCA+"
LABEL2[label2 ==  BA] <- "SSCA+/FSCA-"
mydf = data.frame(M1[,c("SSCA","FSCA")],LABEL2)
q <- ggplot(mydf, aes(x = SSCA, y =  FSCA,
                      color = LABEL2)) +
  geom_point(alpha = 0.2,cex = 1)+
  #Dscale_colour_manual(values = c("green" , "cyan","blue","red") )+
   theme_bw()+ geom_hex(bins=120) + scale_fill_viridis()+
  
  
  theme(legend.position="bottom")+
  ylab("SSCA")+
  xlab("FSCA")+
  guides(colour = guide_legend(override.aes = list(alpha=1, size = 2)))
q 

# 
# q <- ggplot(mydf, aes(x = SSCA, y =  FSCA
#                       )) +
#   geom_point(alpha = 0.2,cex = 1)+
#   #Dscale_colour_manual(values = c("green" , "cyan","blue","red") )+
#   theme_bw()+ geom_hex(bins=120) + scale_fill_viridis()+
#   
#   
#   theme(legend.position="bottom")+
#   ylab("SSCA")+
#   xlab("FSCA")+
#   guides(colour = guide_legend(override.aes = list(alpha=1, size = 2)))
# q 
#