###################################################
### chunk number 1: loadPackage
###################################################
library(flowCore)
library(flowViz)


###################################################
### chunk number 2: Read.fcs
###################################################
myfile <- "data/a01"
fcs1 <- read.FCS(myfile)
fcs1


###################################################
### chunk number 3: Exploring eval=FALSE
###################################################
## exprs(fcs1)[1:10,]
## parameters(fcs1)
## pData(parameters(fcs1))
## summary(fcs1)


###################################################
### chunk number 4: SearchKeywords eval=FALSE
###################################################
## keyword(fcs1, c("$P1E", "$P2E", "$P3E", "$P4E"))


###################################################
### chunk number 5: Read.fcs2 eval=FALSE
###################################################
## fcs1 <- read.FCS(myfile, transformation = "linearize")
## summary(fcs1)
## fcs2 <- read.FCS(myfile, transformation = "scale")
## summary(fcs2)


###################################################
### chunk number 6: ReadFiles2 eval=FALSE
###################################################
## read.FCS(file.name, alter.names = TRUE) 


###################################################
### chunk number 7:  eval=FALSE
###################################################
## plot(fcs1)


###################################################
### chunk number 8: Plot1
###################################################
plot(fcs1)


###################################################
### chunk number 9: Plot2 eval=FALSE
###################################################
## plot(fcs1, c("FL3-H", "FL4-H"))


###################################################
### chunk number 10:  eval=FALSE
###################################################
## plot(fcs1, "SSC-H", breaks=256)


###################################################
### chunk number 11: Frames
###################################################
ffiles <- list.files("data", pattern = "^[a-h]")
fset <- read.flowSet(ffiles, path = "data")
fset


###################################################
### chunk number 12:  eval=FALSE
###################################################
## aliquot <- toupper(substr(ffiles, 1, 1))
## time <- c(-8, 0, 5, 12, 19, 27, 32, 39, 46, 53, 60, 67)[as.numeric(substr(ffiles, 2, 3))]
## pdata <- data.frame(aliquot, time)
## rownames(pdata) <- ffiles
## write.table(pdata, file = "data/phenodata.txt")


###################################################
### chunk number 13: 
###################################################
pdata <- read.table("data/phenodata.txt", header = TRUE)
pdata <- as(pdata, "AnnotatedDataFrame")
phenoData(fset) <- pdata
fset


###################################################
### chunk number 14: metaData
###################################################
varMetadata(phenoData(fset))[, "labelDescription"] <- 
    c("Aliquot", "Days Past Transplant")


###################################################
### chunk number 15: 
###################################################
fset


###################################################
### chunk number 16:  eval=FALSE
###################################################
## varMetadata(phenoData(fset))[, "labelDescription"] <- 
##     c("Aliquot", "Days Past Transplant")


###################################################
### chunk number 17: 
###################################################
fset$"a01"
fset[[1]]


###################################################
### chunk number 18: fsApply1 eval=FALSE
###################################################
## fsApply(fset, each_col, IQR)


###################################################
### chunk number 19: fsApply2 eval=FALSE
###################################################
## fsApply(fset, function(x) apply(x, 2, IQR), use.exprs = TRUE)


###################################################
### chunk number 20: Transfo1
###################################################
summary(transform(fset[[1]], log.FL1.H = log(`FL1-H`)))


###################################################
### chunk number 21: Transfo2
###################################################
summary(transform("FL1-H" = log) %on% fset[[1]])


###################################################
### chunk number 22:  eval=FALSE
###################################################
## splom(fset$h08[, 1:5])


###################################################
### chunk number 23:  eval=FALSE
###################################################
## splom(transform("FSC-H" = asinh, "SSC-H" = asinh,
##                 "FL1-H" = asinh, "FL2-H" = asinh,
##                 "FL3-H" = asinh) %on% fset$h08[, 1:5])


###################################################
### chunk number 24: 
###################################################
plot(splom(fset$h08[, 1:5], main = "Untransformed"), split = c(1, 1, 2, 1))
plot(splom(transform("FSC-H" = asinh, "SSC-H" = asinh,
                     "FL1-H" = asinh, "FL2-H" = asinh,
                     "FL3-H" = asinh) %on% fset$h08[, 1:5], , main = "asinh Transformed"),
     split = c(2, 1, 2, 1), newpage = FALSE)


###################################################
### chunk number 25: 
###################################################
fset.trans <- 
    transform("FSC-H" = asinh, "SSC-H" = asinh,
              "FL1-H" = asinh, "FL2-H" = asinh,
              "FL3-H" = asinh, "FL2-A" = asinh, 
              "FL4-H" = asinh) %on% fset  


###################################################
### chunk number 26: 
###################################################
xyplot(fset.trans[[3]])


###################################################
### chunk number 27: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 28: 
###################################################
qqmath(~`FSC-H` | factor(time), fset.trans, groups = aliquot, 
       f.value = ppoints(500), type = "l")


###################################################
### chunk number 29: 
###################################################
plot(trellis.last.object())


###################################################
### chunk number 30: 
###################################################
qqmath(~`FSC-H` | aliquot, fset.trans, f.value = ppoints(500), 
       type = c("l", "g"), subset = (time == 46))


###################################################
### chunk number 31: rectGate
###################################################
rectgate <- rectangleGate(filterId="Rectangle", 
                          "SSC-H" = c(2.5, 6.5), "FSC-H" = c(4, 5.75))


###################################################
### chunk number 32: 
###################################################
kmfilt <- kmeansFilter("kmfilt", "FSC-H" = c("Low", "High"))
bvnormfilt <- norm2Filter(filterId = "BVNorm", "FSC-H", "SSC-H", scale=2)
bvnormfilt


###################################################
### chunk number 33: 
###################################################

rectgate.results <- filter(fset.trans, rectgate)
kmfilt.results <- filter(fset.trans, kmfilt)
bvnormfilt.results <- filter(fset.trans, bvnormfilt)

summary(rectgate.results[[1]])



###################################################
### chunk number 34: 
###################################################
summary(kmfilt.results[[1]])
lapply(head(kmfilt.results, 4), summary)


###################################################
### chunk number 35: Norm2Filter
###################################################

smaller <- Subset(fset.trans, bvnormfilt)
fset.trans[[1]]
smaller[[1]]



###################################################
### chunk number 36: Split
###################################################
split(smaller[[1]], kmfilt)


###################################################
### chunk number 37: CombineFilter
###################################################

rectgate & bvnormfilt
rectgate | bvnormfilt
!bvnormfilt



###################################################
### chunk number 38:  eval=FALSE
###################################################
## xyplot(`FSC-H` ~ `SSC-H`, fset.trans, subset = (aliquot == "H" & time == 46),
##        filter = bvnormfilt)


###################################################
### chunk number 39: 
###################################################
xyplot(`FSC-H` ~ `SSC-H`, fset.trans, subset = (aliquot == "H" & time == 46),
       smooth = FALSE, pch = ".", cex = 2, filter = kmfilt)


###################################################
### chunk number 40: 
###################################################
plot(xyplot(`FSC-H` ~ `SSC-H`, fset.trans, subset = (aliquot == "H" & time == 46),
            filter = bvnormfilt, main = "a"), 
     split = c(1, 1, 2, 1))
plot(xyplot(`FSC-H` ~ `SSC-H`, fset.trans, subset = (aliquot == "H" & time == 46),
            smooth = FALSE, pch = ".", cex = 2, filter = kmfilt, main = "b"),
     split = c(2, 1, 2, 1), newpage = FALSE)


###################################################
### chunk number 41:  eval=FALSE
###################################################
## fset.high <- fsApply(fset.trans, function(x) split(x, kmfilt)$High)
## xyplot(`FSC-H` ~ `SSC-H` | factor(time), fset.high, subset = aliquot == "H",
##        filter = bvnormfilt)
## xyplot(`FSC-H` ~ `SSC-H` | factor(time), fset.trans, subset = aliquot == "H",
##        filter = bvnormfilt)


###################################################
### chunk number 42: 
###################################################
sessionInfo()


