GettingStartedWithFlowStats.R
In flowStats: Statistical methods for the analysis of flow cytometry data

### R code from vignette source 'GettingStartedWithFlowStats.Rnw'

###################################################
### code chunk number 1: loadGvHD
###################################################
library(flowCore)
library(flowStats)
library(flowWorkspace)
library(ggcyto)
data(ITN)


###################################################
### code chunk number 2: transform
###################################################
require(scales)
gs <- GatingSet(ITN)
trans.func <- asinh
inv.func <- sinh
trans.obj <- trans_new("myAsinh", trans.func, inv.func)
transList <- transformerList(colnames(ITN)[3:7], trans.obj)
gs <- transform(gs, transList)


###################################################
### code chunk number 3: lymphGate
###################################################
lg <- lymphGate(gs_cyto_data(gs), channels=c("CD3", "SSC"),
                preselection="CD4", filterId="TCells",
                scale=2.5)
gs_pop_add(gs, lg)
recompute(gs)


###################################################
### code chunk number 4: lymphGatePlot
###################################################
ggcyto(gs, aes(x = CD3, y = SSC)) + geom_hex(bins = 32) + geom_gate("TCells")


###################################################
### code chunk number 5: variation
###################################################
require(ggridges)
require(gridExtra)
pars <- colnames(gs)[c(3,4,5,7)]

data <- gs_pop_get_data(gs, "TCells")
plots <- lapply(pars, function(par)
  as.ggplot(ggcyto(data, aes(x = !!par, fill = GroupID)) +
    geom_density_ridges(mapping = aes(y = name), alpha = 0.4) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank()) +
    facet_null()
  ))
do.call("grid.arrange", c(plots, nrow=1))


###################################################
### code chunk number 6: quad_gate
###################################################
qgate <- quadrantGate(gs_pop_get_data(gs, "TCells"), stains=c("CD4", "CD8"),
                      filterId="CD4CD8", sd=3)
gs_pop_add(gs, qgate, parent = "TCells")
recompute(gs)


###################################################
### code chunk number 7: quad_gate_norm
###################################################
gs <- normalize(gs, populations=c("CD4+CD8+", "CD4+CD8-", "CD4-CD8+", "CD4-CD8-"), 
                     dims=c("CD4", "CD8"), minCountThreshold = 50)


###################################################
### code chunk number 8: variation_norm
###################################################
data <- gs_pop_get_data(gs, "TCells")
plots <- lapply(pars, function(par)
  as.ggplot(ggcyto(data, aes(x = !!par, fill = GroupID)) +
    geom_density_ridges(mapping = aes(y = name), alpha = 0.4) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank()) +
    facet_null()
  ))
do.call("grid.arrange", c(plots, nrow=1))


###################################################
### code chunk number 9: quad_gate_norm_plot
###################################################
ggcyto(gs_pop_get_data(gs, "TCells"), aes(x=CD4, y=CD8)) +
  geom_hex(bins=32) +
  geom_gate(gs_pop_get_gate(gs, "CD4-CD8-")) +
  geom_gate(gs_pop_get_gate(gs, "CD4-CD8+")) +
  geom_gate(gs_pop_get_gate(gs, "CD4+CD8-")) +
  geom_gate(gs_pop_get_gate(gs, "CD4+CD8+"))


###################################################
### code chunk number 10: rangeGate
###################################################
CD69rg <- rangeGate(gs_cyto_data(gs), stain="CD69",
                    alpha=0.75, filterId="CD4+CD8-CD69", sd=2.5)
gs_pop_add(gs, CD69rg, parent="CD4+CD8-")
# gs_pop_add(gs, CD69rg, parent="CD4+CD8-", name = "CD4+CD8-CD69-", negated=TRUE)
recompute(gs)


###################################################
### code chunk number 11: rangeGatePlot
###################################################
ggcyto(gs_pop_get_data(gs, "CD4+CD8-"), aes(x=CD69, fill = GroupID)) +
    geom_density_ridges(mapping = aes(y = name), alpha = 0.4) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank()) +
    geom_vline(xintercept = CD69rg@min) +
    facet_null() +
    ggtitle("CD4+")


###################################################
### code chunk number 12: rangeGate_norm
###################################################
gs <- normalize(gs, populations=c("CD4+CD8-CD69"), 
                     dims=c("CD69"), minCountThreshold = 50)


###################################################
### code chunk number 13: rangeGate_norm_plot
###################################################
ggcyto(gs_pop_get_data(gs, "CD4+CD8-"), aes(x=CD69, fill = GroupID)) +
    geom_density_ridges(mapping = aes(y = name), alpha = 0.4) +
    theme(axis.title.y = element_blank(),
          axis.text.y = element_blank(),
          axis.ticks.y = element_blank()) +
    geom_vline(xintercept = CD69rg@min) +
    facet_null() +
    ggtitle("CD4+")


###################################################
### code chunk number 14: createData
###################################################
dat <- gs_cyto_data(gs)



###################################################
### code chunk number 15: rawData
###################################################
autoplot(dat, "CD4", "CD8") +
  ggtitle("Experimental Dataset")


###################################################
### code chunk number 16: createControlData
###################################################
datComb <- as(dat,"flowFrame")
subCount <- nrow(exprs(datComb))/length(dat)
	sf <- sampleFilter(filterId="mySampleFilter", size=subCount)
	fres <- filter(datComb, sf)
	ctrlData <- Subset(datComb, fres)
	ctrlData <- ctrlData[,-ncol(ctrlData)] ##remove the  column name "original"
	


###################################################
### code chunk number 17: BinControlData
###################################################
minRow=subCount*0.05
refBins<-proBin(ctrlData,minRow,channels=c("CD4","CD8"))
	


###################################################
### code chunk number 18: controlBinsPlot
###################################################
plotBins(refBins,ctrlData,channels=c("CD4","CD8"),title="Control Data")



###################################################
### code chunk number 19: binSampleData
###################################################
sampBins <- fsApply(dat,function(x){
		   binByRef(refBins,x)
		   })


###################################################
### code chunk number 20: pearsonStat
###################################################
pearsonStat <- lapply(sampBins,function(x){
		      calcPearsonChi(refBins,x)
                     })


###################################################
### code chunk number 21: Roderers PBin metric
###################################################
sCount <- fsApply(dat,nrow)
pBStat <-lapply(seq_along(sampBins),function(x){
		calcPBChiSquare(refBins,sampBins[[x]],subCount,sCount[x])
		})


###################################################
### code chunk number 22: plotBinsresiduals
###################################################
par(mfrow=c(4,4),mar=c(1.5,1.5,1.5,1.5))

plotBins(refBins,ctrlData,channels=c("CD4","CD8"),title="Control Data")
patNames <-sampleNames(dat)
tm<-lapply(seq_len(length(dat)),function(x){
		plotBins(refBins,dat[[x]],channels=c("CD4","CD8"),
			title=patNames[x],
			residuals=pearsonStat[[x]]$residuals[2,],
			shadeFactor=0.7)
		
		}
      )


###################################################
### code chunk number 23: chiSqstatisticvalues
###################################################
library(xtable)
chi_Square_Statistic <- unlist(lapply(pearsonStat,function(x){
		x$statistic
		}))

pBin_Statistic <-unlist(lapply(pBStat,function(x){
                x$pbStat
						                        })) 
	
frame <- data.frame(chi_Square_Statistic, pBin_Statistic)
rownames(frame) <- patNames
 	


###################################################
### code chunk number 24: GettingStartedWithFlowStats.Rnw:395-396
###################################################
print(xtable(frame))


###################################################
### code chunk number 25: GettingStartedWithFlowStats.Rnw:402-403
###################################################
toLatex(sessionInfo())