###################################################
### chunk number 1: Load Libraries
###################################################
library(flowCore)
library(flowQ)
library(flowViz)
library(flowStats)


###################################################
### chunk number 2: ReadFlowSet
###################################################

flowData <- read.flowSet(path="." , phenoData="annotation.txt",transformation=FALSE)
sampleNames(flowData) <-  as.character(pData(flowData)[,"PatientID"])



###################################################
### chunk number 3: ShowName/description fields
###################################################
pData(parameters(flowData[[1]]))


###################################################
### chunk number 4: UpdateDescriptionFields
###################################################
for(i in seq_len(length(flowData))){
pData(parameters(flowData[[i]]))[,"desc"] <- c("NA","NA", "CD8", "CD69",
                                               "CD4","CD3", "HLADr","NA")
}
flowData[[1]]


###################################################
### chunk number 5: Plot Forward/Side scatter
###################################################
print(xyplot(`SSC-A` ~ `FSC-A`,data=flowData))



###################################################
### chunk number 6: PlotUntransformedData
###################################################

print(splom(flowData[[1]]))



###################################################
### chunk number 7: transformData
###################################################
tData <- transform(flowData, transformList(colnames(flowData)[3:7], asinh))


###################################################
### chunk number 8: PlotTransformedData
###################################################
print(splom(tData[[1]]))



###################################################
### chunk number 9: QABoundaryEvents
###################################################

dest <- "QAResults"
qp1 <- qaProcess.marginevents(tData, channels=c("FSC-A","SSC-A"), outdir=dest,
                              pdf=TRUE,sum.dimensions=c(5,3))



###################################################
### chunk number 10: QAProcessCellNumber
###################################################
qp2 <- qaProcess.cellnumber(tData, outdir=dest, cFactor=2, pdf=TRUE)



###################################################
### chunk number 11: QAProcessTimeLine
###################################################
qp3 <- qaProcess.timeline(tData, channels="FSC-A", outdir=dest, cutoff=1,
                          pdf=TRUE)

	

###################################################
### chunk number 12: WriteHTMLResports
###################################################
 url <- writeQAReport(tData, processes=list(qp1, qp2,qp3), outdir=dest,
                       pdf=TRUE)



###################################################
### chunk number 13: start the work flow
###################################################
colnames(flowData) <- c("FSC", "SSC","CD8","CD69", "CD4","CD3","HLADr", "Time")
wf <- workFlow(flowData, name="biocExample")
wf


###################################################
### chunk number 14: Transform the data
###################################################
tf <- transformList(colnames(Data(wf[["base view"]]))[3:7], asinh,
                    transformationId="asinh")
add(wf, tf)
wf


###################################################
### chunk number 15: Boundary filter FSC SSC channels
###################################################
boundFilt <-boundaryFilter(filterId="boundFilt", x=c("FSC","SSC"))
add(wf,boundFilt,parent="asinh")


###################################################
### chunk number 16: lymphGate
###################################################
lg <- lymphGate(Data(wf[["boundFilt+"]]), channels=c("FSC", "SSC"),
                preselection="CD3", filterId="TCells", eval=FALSE,
                scale=2.5)
add(wf, lg$n2gate, parent="boundFilt+")


###################################################
### chunk number 17: lymphGatePlot
###################################################

xyplot(SSC ~ FSC| PatientID, wf[["TCells+"]],
      par.settings=list(gate=list(col="red",
         fill="red", alpha=0.3)))

###################################################
### chunk number 18:  Plot flowFrame before & after Lymph gate 
###################################################

xyplot(SSC ~ FSC ,Data(wf[["TCells+"]])[[1]],smooth=F,pch=1.3,col="red")
dev.new()
xyplot(SSC ~ FSC ,Data(wf[["TCells-"]])[[1]],smooth=F,pch=1.3,col="red")

		 
###################################################
### chunk number 19: Density plots for CD3, CD4, CD8
###################################################

print(densityplot(PatientID~., Data(wf[["TCells+"]]), channels=c("CD3","CD4","CD8") ,
                  groups=GroupID, scales=list(y=list(draw=F)),
                  filter=lapply(c("CD3","CD4","CD8"), curv1Filter), layout=c(3,1)))


###################################################
### chunk number 20: Data after normalization
###################################################
pars <- colnames(Data(wf[["base view"]]))[c(3,4,5,6)]
norm <- normalization(normFun=function(x, parameters, ...)
                      warpSet(x, parameters,...),
                      parameters=pars,
                      normalizationId="Warping")
add(wf, norm, parent="TCells+")


###################################################
### chunk number 21: PlotNormalized CD3 ,CD4, CD8
###################################################
print(densityplot(PatientID~., Data(wf[["Warping"]]),
                  channels=c("CD3","CD4","CD8"), groups=GroupID,
                  scales=list(y=list(draw=F)),
                  filter=lapply(c("CD3","CD4","CD8"), curv1Filter),
                  layout=c(3,1)))


###################################################
### chunk number 22: QuadrantgateCD3+CD4+  (Helper/inducer cells T cells)
###################################################
qgate <- quadrantGate(Data(wf[["Warping"]]), stains=c("CD3", "CD4"),
                      plot=FALSE, filterId="CD3CD4")
add(wf, qgate, parent="Warping")


###################################################
### chunk number 23: PlotQuadrantGateCD3+CD4+
###################################################
xyplot(CD3 ~ CD4 | PatientID, wf[["CD3+CD4+"]] )


###################################################
### chunk number 24: QuadrantGateCD3+CD8+ (Cytotoxic/suppressor T cells)
###################################################
qgate <- quadrantGate(Data(wf[["Warping"]]), stains=c("CD3", "CD8"),
                      plot=FALSE, filterId="CD3CD8")
add(wf, qgate, parent="Warping")


###################################################
### chunk number 25: PlotQuadrantGateCD3+CD4+
###################################################
xyplot(CD3 ~ CD8 | PatientID, wf[["CD3+CD8+"]])


###################################################
### chunk number 26: HLADrActivatedCD3+CD4+
###################################################
HLADr1 <- rangeGate(Data(wf[["CD3+CD4+"]]), stain="HLADr", plot=FALSE,
                    alpha=0.75, filterId="CD3+CD4+HLAct")
add(wf, HLADr1, parent="CD3+CD4+")


###################################################
### chunk number 27: PlotHLADrActivatedCD3+CD4+
###################################################
densityplot(PatientID ~ HLADr, Data(wf[["CD3+CD4+"]]), refline=HLADr1@min)


###################################################
### chunk number 28: HLADrActivatedCD3+CD8+
###################################################
HLADr2 <- rangeGate(Data(wf[["CD3+CD8+"]]), stain="HLADr", plot=FALSE,
                    alpha=0.75, filterId="CD3+CD8+HLAct")
add(wf, HLADr2, parent="CD3+CD8+")


###################################################
### chunk number 29: PlotHLADrActivatedCD3+CD8+
###################################################
densityplot(PatientID ~ HLADr, Data(wf[["CD3+CD8+"]]), refline=HLADr2@min)


###################################################
### chunk number 30: Percent of HLADrActivatedCD3+CD4+ to original population
###################################################
pr <-fsApply(Data(wf[["CD3+CD4+HLAct-"]]),nrow)*100 /fsApply(Data(wf[["boundFilt+"]]),nrow)
res <- data.frame(pr, pData(Data(wf[["CD3+CD4+"]]))[c("GroupID","PatientID")])


barchart(reorder(PatientID,as.numeric(factor(GroupID))) ~ pr,data=res,
                         groups=GroupID,stack=TRUE,
                         auto.key = list(points = FALSE,
                                         rectangles = TRUE,
                                         space = "right"),
                          main="HlADr Activated CD3+CD4+ T cells",
                          xlab="Percentage of activated cells"
         )


###################################################
### chunk number 31: Percentage of HLADRrActivatedCD3+CD8+ to orignal population
###################################################
pr <- fsApply(Data(wf[["CD3+CD8+HLAct-"]]),nrow)*100 /fsApply(Data(wf[["boundFilt+"]]),nrow)
res <- data.frame(pr,pData(Data(wf[["CD3+CD8+"]]))[c("GroupID","PatientID")])

barchart(reorder(PatientID,as.numeric(factor(GroupID))) ~ pr,data=res,
                         groups=GroupID,stack=TRUE,
                         auto.key = list(points = FALSE,
                                         rectangles = TRUE,
                                         space = "right"),
                          main="Activated  CD3+CD8+ T cells",
                          xlab=" Percentage of activated cells "
         )


###################################################
### chunk number 32: PlotWorkFlow
###################################################
plot(wf)