## ----,eval=TRUE,echo=TRUE------------------------------------------------ # Common files FCSPATH="/data/OpenCyto/data/FCS/" WORKSPACE="/data/OpenCyto/data/workspace/080 batch 0882.xml" OUTPATH="~" MANUAL<-file.path(OUTPATH,"data","manual") AUTO<-file.path(OUTPATH,"data","auto") TEMPLATE="/data/OpenCyto/data/template/gt_080.csv" ANNOTATIONS="/data/OpenCyto/data/workspace/pd_submit.csv" # Move the template into your home directory if(!file.exists("~/data")){ dir.create("~/data") } file.copy(TEMPLATE,file.path("~/data",basename(TEMPLATE)),overwrite = TRUE) file.copy(ANNOTATIONS,file.path("~/data",basename(ANNOTATIONS)),overwrite = TRUE) TEMPLATE<-"~/data/gt_080.csv" ANNOTATIONS<-"~/data/pd_submit.csv" ## ----eval=TRUE----------------------------------------------------------- require(openCyto) require(ggplot2) require(plyr) require(reshape2) require(data.table) ## ----eval=TRUE, echo=TRUE------------------------------------------------ ws<-openWorkspace(WORKSPACE) show(ws) ## ----eval=TRUE, echo=TRUE,message=FALSE,results='hide',warning=FALSE,error=FALSE---- # Parse the worspace only if we haven't already done so, and save it. if(!file.exists(MANUAL)){ G<-parseWorkspace(ws,isNcdf=TRUE,name="0882 Samples",path=FCSPATH) save_gs(G,MANUAL) }else{ G<-load_gs(MANUAL) } ## ----eval=TRUE,echo=TRUE------------------------------------------------- show(G) ## ----echo=TRUE, eval=TRUE------------------------------------------------ G[1:10] # Subset by numeric index.. not the preferred way of doing things. sampleNames(G) # Get the sample names G[sampleNames(G)[1:10]] # Subset by sample name.. preferred. ## ----echo=TRUE,eval=TRUE------------------------------------------------- raw<-seq(1,2.5e5,by=1000) #Input data range on the raw scale transformed<-lapply(getTransformations(G[[1]]),function(x)matrix(x(raw))) # transform with each transformation transformed<-do.call(cbind,transformed) # group output colnames(transformed)<-gsub("^ ","",names(getTransformations(G[[1]]))) transformed<-melt(transformed) transformed<-(cbind(transformed,raw)) setnames(transformed,c("index","parameter","transformed","raw")) ggplot(data.frame(raw,transformed))+geom_line(aes(x=raw,y=transformed))+facet_wrap(~parameter)+theme_minimal() ## ----eval=TRUE, echo=TRUE------------------------------------------------ getCompensationMatrices(G[[1]]) #Compensation for the first GatingHierarchy ggplot(melt(getCompensationMatrices(G[[1]])@spillover,value.name = "Coefficient"))+geom_tile(aes(x=Var1,y=Var2,fill=Coefficient))+scale_fill_continuous(guide="colourbar")+theme(axis.text.x=element_text(angle=45,hjust=1)) ## ----eval=TRUE, echo=TRUE------------------------------------------------ getNodes(G) ## ----path_auto----------------------------------------------------------- getNodes(G,path="auto") ## ----getPopStats--------------------------------------------------------- freq_gs<-getPopStats(G[1:5],statistic="freq") head(freq_gs) counts_gh<-getPopStats(G[1:5],statistic="count") head(counts_gh) ## ----plot---------------------------------------------------------------- # The gating tree for this gating set plot(G[[1]],boolean=FALSE) ## ----hiding_nodes,echo=FALSE,results="hide"------------------------------ # setNode doesn't have a "GatingSet" interface, so we need to loop over all samples lapply(G,function(x)setNode(x,y = "Not 4+/IFNg+", FALSE)) #set the node as hidden plot(G) ## ----plotGate------------------------------------------------------------ #plot one specific gate from one sample p1<-plotGate(G[[2]],"4+/TNFa+",arrange=FALSE) #plot one gate from two samples - automatic faceting p2<-plotGate(G[2:3],"4+/TNFa+",arrange=FALSE) #oops! population names must be unique try(plotGate(G[1:2],"TNFa+")) #binning using xbin (or not) p4<-plotGate(G[[2]], "4+/TNFa+",xbin=128,arrange=FALSE) #smaller bins (default is 32), more detail # use gridExtra for layout require(gridExtra) grid.arrange(p1[[1]],p2,p4[[1]]) # Single plot is a list with element 1 being the plot object. Multiple plots are trellis objects directly. #plot an entire layout (boolean gates are skipped automatically) p3<-plotGate(G[[2]]) # a little slower since we read all the event-level data. Arrange is TRUE by default ## ----getData------------------------------------------------------------- cd3 <- getData(G[[1]],"3+") cd3 ## ----getData_flowSet----------------------------------------------------- cd3_set<-getData(G[1:10],"3+") cd3_set ## ----ncdfFlowSet--------------------------------------------------------- cd3_set@file #where the NetCDF file lives flowData(G) #Return the entire ncdfFlowSet associated with a GatingSet flowData(G[1:5]) #Subsettig works, and still points to the same file. ## ------------------------------------------------------------------------ # Get cells that express IL2 or IFNg and are CD8+. # map tells us how the node names map to markers on channels. sce<-getSingleCellExpression(G[1:5],c("8+/IL2+","8+/IFNg+"),map=list("8+/IL2+"="IL2","8+/IFNg+"="IFNg")) str(sce) # list of 5 n x 2 matrices. ## ----keywords------------------------------------------------------------ require(flowIncubator) # the `getKeywords` function is in the `flowIncubator` package at http://www.github.com/RGLab/flowIncubator names(keyword(G[[1]])) # valid keywords keyword_vars<-c("$FIL","Stim","Sample Order","EXPERIMENT NAME") #relevant keywords pd<-data.table(getKeywords(G,keyword_vars)) #extract relevant keywords to a data table annotations<-data.table:::fread(ANNOTATIONS) # read the annotations from flowrepository setnames(annotations,"File Name","name") setkey(annotations,"name") setkey(pd,"name") head(pd) head(annotations) # We match on the name column, which matches the sampleNames(G). pd<-data.frame(annotations[pd]) #data.table style merge setnames(pd,c("Timepoint","Individual"),c("VISITNO","PTID")) #Rename Timepoint and Indivudual to VISITNO and PTID which are used in the template files pData(G)<-pd #annotate ## ----clone,message=FALSE------------------------------------------------- automated<-clone(G) ## ----remove_nodes-------------------------------------------------------- Rm("S",automated) #Rm is the API to remove a node and all nodes beneath it. ## ----load_gating_template------------------------------------------------ gating_template<-gatingTemplate(TEMPLATE) ## ----plot_template------------------------------------------------------- plot(gating_template) ## ----boundary,warning=FALSE---------------------------------------------- options("warn"=-1) head(fread(TEMPLATE)[1,]) ## ----available_methods--------------------------------------------------- listgtMethods() ## ----flowframe----------------------------------------------------------- getData(automated[[1]],use.exprs=FALSE) #use.exprs=FALSE don't retrieve events.. faster. ## ----template_boundary--------------------------------------------------- gating_template@nodeData@data[["/boundary"]] # node named "/boundary" gating_template@edgeData@data[["root|/boundary"]] #edge from root to the boundary node str(gating_template@edgeData@data[["root|/boundary"]]$gtMethod) ## ------------------------------------------------------------------------ head(fread(TEMPLATE)[3:6,]) # Thre four entries in the csv template. ## ----lymph_csv----------------------------------------------------------- head(fread(TEMPLATE)[7,]) # The lymph population in the csv template. ## ----cd8_csv------------------------------------------------------------- head(fread(TEMPLATE)[9,]) # The four entries in the csv template. ## ----expanded_template_cd8----------------------------------------------- # Two nodes for the single input row "cd8+/-" gating_template@nodeData@data[["/boundary/singlet/viable/nonDebris/lymph/cd3/cd8+"]] gating_template@nodeData@data[["/boundary/singlet/viable/nonDebris/lymph/cd3/cd8-"]] ## ----run_gating_to_tnfa-------------------------------------------------- #gate the first 5 samples nodes(gating_template) gating(gating_template,automated[1:10]) ## ----plot_cd4------------------------------------------------------------ plotGate(automated[[1]],c("cd8-","cd4pos","cd4"),gpar=list(nrow=1)) #gpar adjusts graphical parameters so we can layout by row. ## ----echo=FALSE,results='hide'------------------------------------------- plot(gating_template) ## ------------------------------------------------------------------------ fread(TEMPLATE)[14,] ## ------------------------------------------------------------------------ fread(TEMPLATE)[20,] ## ------------------------------------------------------------------------ args(registerPlugins) ?registerPlugins ## ----plugin_preprocessing------------------------------------------------ # The gate method will subset the flowFrame on the "control" events and gate those, returning a "negative control" gate. # This could work equally well with FMOs. .negGate <- function(fr, pp_res, xChannel=NA, yChannel=NA, filterId="ppgate", ...){ my_gate <- tailgate(fr[pp_res,],channel=yChannel, filter_id=filterId, ...) return(my_gate) } registerPlugins(fun=.negGate,methodName='negGate',dep=NA) # Identifies which events belong to the negative control sample vs the treatment sample. # groupBy includes PTID # Treatment information comes from the pData of the flowSet # passes a vector of 0/1 for treamtent / control events .ppnegGate <- function(fs, gs, gm, xChannel, yChannel, groupBy, isCollapse, ...) { d <- c() for(i in c(1:length(fs))) { d = c(d,rep.int(pData(fs[i])$control,nrow(exprs(fs[[i]])))) } return(as.logical(d)) } registerPlugins(fun=.ppnegGate, methodName='ppnegGate', dep=NA, "preprocessing")