## ----framework_figure, echo=FALSE,fig.cap='The OpenCyto Gating Framework is a collection of R/BioConductor packages for easily building reproducible flow data analysis pipelines.',message=FALSE,fig.align='center',fig.width=8,cache=TRUE---- require(png) require(grid) img<-readPNG("resources/Framework.png") grid.raster(img) ## ----eval=FALSE---------------------------------------------------------- ## require(BiocInstaller) ## biocLite("openCyto") ## ----eval=FALSE---------------------------------------------------------- ## require(devtools) ## packages<-c("RGLab/flowStats","RGLab/flowCore","RGLab/flowViz","RGLab/ncdfFlow","RGLab/flowWorkspace","RGLab/openCyto") ## install_github(packages,quick=TRUE) ## ----libs,echo=FALSE,message=FALSE,warning=FALSE------------------------- require(openCyto) require(data.table) require(reshape2) require(clue) require(ggplot2) ## ----echo=1,results='markup',results='hide'------------------------------ ws<-openWorkspace("data/workspace/080 batch 0882.xml") ws ## ----parseWorkspace_echo,echo=TRUE,eval=FALSE---------------------------- ## gating_set<-parseWorkspace(ws,name="0882 Samples",path="data/FCS/",isNcdf=TRUE) ## ----parseWorkspace,echo=FALSE,results='hide',eval=TRUE------------------ if(!file.exists("data/manual_gating/")){ gating_set<-parseWorkspace(ws,name="0882 Samples",path="data/FCS/",isNcdf=TRUE,overwrite=TRUE) save_gs(gating_set,path="data/manual_gating") }else{ gating_set<-load_gs("data/manual_gating/") } ## ----eval=FALSE,echo=TRUE,results='hide'--------------------------------- ## save_gs(gating_set,path="data/manual_gating") ## ----clean_manual,echo=FALSE,results='hide',warning=FALSE,message=FALSE---- Rm("Not 4+",gating_set) ## ----vis_manual_gates,echo=TRUE,eval=TRUE,results='hide',cache=FALSE,message=FALSE,warning=FALSE,fig.align='center',fig.cap='Layout of manual gates'---- plotGate(gating_set[[1]],xbin=16,gpar=list(ncol=5)) # Binning for faster plotting ## ----plot_manual_tree,echo=2,results="hide",message=FALSE,warning=FALSE,fig.cap=''---- plot(gating_set) ## ----keyword,echo=c(2,3,4,8,9,10),eval=TRUE------------------------------ getKeywords<-function(gs,kv){ r<-as.data.frame(do.call(cbind,lapply(kv,function(k){ keyword(gs,k)[1] }))) data.table::setnames(r,"$FIL","name") r } keyword_vars<-c("$FIL","Stim","Sample Order","EXPERIMENT NAME") #relevant keywords pd<-data.table(getKeywords(gating_set,keyword_vars)) #extract relevant keywords to a data table annotations<-data.table:::fread("data/workspace/pd_submit.csv") # read the annotations from flowrepository setnames(annotations,"File Name","name") setkey(annotations,"name") setkey(pd,"name") pd<-data.frame(annotations[pd]) #data.table style merge setnames(pd,c("Timepoint","Individual"),c("VISITNO","PTID")) pData(gating_set)<-pd #annotate gating_subset<-gating_set[subset(pd,!is.na(Condition))$name] #subset only the GAG and negctrl ## ----keywords_tbl,eval=TRUE,echo=FALSE,results='asis',cache=FALSE-------- knitr::kable(head(subset(na.omit(pd)[,c(1:5)],PTID%in%"080-17"),4),row.names = FALSE) ## ----copy_and_save,eval=FALSE,message=FALSE,warining=FALSE,results='hide'---- ## auto_gating<-clone(gating_subset) ## Rm("S",auto_gating) ## save_gs(auto_gating,path="data/autogating",overwrite=TRUE) ## ----load_empty,eval=TRUE,echo=FALSE,results='hide',message=FALSE,warning=FALSE---- if(!file.exists("data/autogating")){ auto_gating<-clone(gating_subset) try(Rm("S",auto_gating)) save_gs(auto_gating,path="data/autogating") }else{ auto_gating<-load_gs("data/autogating/") } ## ----echo=TRUE,eval=TRUE------------------------------------------------- list.files("data/autogating") ## ----read_template,echo=FALSE,results='asis',cache=FALSE----------------- template<-read.csv("data//template//gt_080.csv",as.is=FALSE,allowEscapes=TRUE,stringsAsFactors=FALSE) template$collapseDataForGating[is.na(template$collapseDataForGating)]<-" " Kmisc::htmlTable(data.frame(template)[1:8,c(1:6,8,9)],attr="width=100%") ## ----read_plot_template,echo=TRUE,message=FALSE,warning=FALSE,results='hide',fig.align="center",fig.width=8,fig.height=5,out.width=600,dpi=200,fig.cap=''---- gt<-gatingTemplate("data/template/gt_080.csv") plot(gt) ## ----gate_subset,echo=TRUE,message=FALSE,warning=FALSE,results='hide',eval=FALSE,fig.cap=''---- ## gating(x = gt, y = auto_gating) ## ## Some output.. ## plot(auto_gating) ## ----plot_autogate_tree,echo=FALSE,message=FALSE,warning=FALSE,fig.align='center',out.width=600,results='hide'---- if(length(getNodes(auto_gating))<1){ gating(x = gt, y = auto_gating) save_gs(auto_gating,path="data//autogating",overwrite = TRUE) } setNode(auto_gating,"nonNeutro",FALSE) setNode(auto_gating,"DebrisGate",FALSE) setNode(auto_gating,"cd4/cd57-/gzb_gate",FALSE) setNode(auto_gating,"cd4/cd57-/prf_gate",FALSE) setNode(auto_gating,"cd8/cd57-/gzb_gate",FALSE) setNode(auto_gating,"cd8/cd57-",FALSE) setNode(auto_gating,"cd4/cd57-",FALSE) setNode(auto_gating,"cd8/cd57-/prf_gate",FALSE) setNode(auto_gating,"cd4neg",FALSE) setNode(auto_gating,"cd4pos",FALSE) setNode(auto_gating,"cd8+",FALSE) setNode(auto_gating,"cd8-",FALSE) plot(auto_gating) ## ----compare_to_manual,echo=TRUE,fig.cap='Automated and Manual Gates for CD4/CD8',cache=FALSE,fig.align='center',fig.width=5,fig.height=2.5,out.width=600---- p1<-plotGate(auto_gating[[1]],c("cd8","cd4"),arrange=FALSE, projections=list("cd4"=c(x="CD8",y="CD4"),"cd8"=c(x="CD8",y="CD4")), main="Automated Gate",path=2)[[1]] p2<-plotGate(gating_subset[[1]],c("8+","4+"), projections=list("4+"=c(x="CD8",y="CD4"),"8+"=c(x="CD8",y="CD4")), arrange=FALSE,main="Manual Gate",path=2)[[1]] grid.arrange(arrangeGrob(p1,p2,ncol=2)) ## ----extract_stats,echo=FALSE,message=FALSE, warning=FALSE,error=FALSE,results='hide',fig.cap='Comparison of Manual vs Automated Gating Cell Subset Counts',fig.width=12,fig.height=5,cache=FALSE,fig.align='center',out.width=800---- auto_stats<-getPopStats(auto_gating,statistic="count") manual_stats<-getPopStats(gating_subset,statistic="count") gates_to_remove<-rownames(manual_stats)[which(sapply(basename(rownames(manual_stats)),function(x)length(which(strsplit(x,"")[[1]]=="+"))>1))] gates_to_remove<-c(gates_to_remove,rownames(manual_stats)[which(sapply(basename(rownames(manual_stats)),function(x)length(which(strsplit(x,"")[[1]]=="-"))>1))],"/S/Lv/L/Not 4+","/S/Lv/L/3+/8+/Granzyme B-","/S/Lv/L/3+/8+/IFN\\IL2","/S/Lv/L/3+/4+/IFN\\IL2","/S/Lv/L/3+/8+/IFN+IL2-","/S/Lv/L/3+/4+/IFN+IL2-","/S/Lv/L/3+/8+/IFN-IL2+", "/S/Lv/L/3+/4+/Granzyme B-","4+/IFN-IL2+","/S/Lv/L/3+/8+/57-", "/S/Lv/L/3+/4+/57-" ) for(i in gates_to_remove){ try(Rm(i,gating_subset),silent=TRUE) } .getCounts<-function(stat="count"){ auto_stats<-getPopStats(auto_gating,statistic=stat) manual_stats<-getPopStats(gating_subset,statistic=stat) #combine melted_auto<-melt(auto_stats) melted_manual<-melt(manual_stats) setnames(melted_manual,c("population","file","counts")) setnames(melted_auto,c("population","file","counts")) melted_auto<-subset(melted_auto,!population%in%c("boundary","root")) melted_auto$population<-factor(melted_auto$population) melted_manual<-subset(melted_manual,!population%in%c("4+/57-","4+/Granzyme B-","4+/IFN+IL2-","4+/IFN-IL2+","4+/IFN\\IL2","8+/57-","8+/Granzyme B-","8+/IFN+IL2-", "8+/IFN-IL2+" ,"8+/IFN\\IL2","root")) melted_manual$population<-factor(melted_manual$population) D<-adist((levels(melted_auto$population)),(levels(melted_manual$population)),ignore.case=TRUE,partial=FALSE) D<-solve_LSAP(t(D)) data.frame(levels(melted_auto$population)[D],levels(melted_manual$population)) levels(melted_auto$population)[D]<-levels(melted_manual$population) merged<-rbind(cbind(melted_auto,method="auto"),cbind(melted_manual,method="manual")) merged } merged_counts<-.getCounts(stat="count") merged_props<-.getCounts(stat="freq") corr.coef<-cor(dcast(merged_counts,file+population~method,value.var = "counts")[,c("auto","manual")],use="complete") corr.coef.freq<-cor(dcast(merged_props,file+population~method,value.var = "counts")[,c("auto","manual")],use="complete") et<-element_text(size=20) et2<-element_text(size=12) p1<-ggplot(dcast(merged_counts,file+population~method,value.var = "counts"))+geom_point(aes(x=auto,y=manual,color=basename(as.character(population))))+scale_y_log10(Kmisc::wrap("Manual Gating Population Count",30))+scale_x_log10(Kmisc::wrap("Autogating Population Count",30))+theme_bw()+geom_abline(lty=3)+theme(axis.text.x=et,axis.text.y=et,axis.title.x=et,axis.title.y=et,legend.text=et2,legend.position="none",plot.title=et)+scale_color_discrete("Population")+geom_text(aes(x=10,y=100000,label=sprintf("rho=%s",signif(corr.coef[1,2],4))),data=data.frame(corr.coef))+ggtitle("Counts") p2<-ggplot(dcast(merged_props,file+population~method,value.var = "counts"))+geom_point(aes(x=auto,y=manual,color=basename(as.character(population))))+scale_y_log10(Kmisc::wrap("Manual Gating Population Proportion",30))+scale_x_log10(Kmisc::wrap("Autogating Population Proportion",30))+theme_bw()+geom_abline(lty=3)+theme(axis.text.x=et,axis.text.y=et,axis.title.x=et,axis.title.y=et,legend.text=et2,legend.position="left",plot.title=et)+scale_color_discrete("Population")+geom_text(aes(x=0.01,y=0.75,label=sprintf("rho=%s",signif(corr.coef.freq[1,2],4))),data=data.frame(corr.coef.freq))+ggtitle("Proportions") grid.draw(cbind(ggplotGrob(p1), ggplotGrob(p2), size="last")) ## ----extract_stats_echo,echo=TRUE,eval=FALSE----------------------------- ## #Extract stats ## auto_stats<-getPopStats(auto_gating,statistic="count") ## manual_stats<-getPopStats(gating_subset,statistic="count") ## ----example_Perforin,echo=FALSE,eval=TRUE,fig.cap='Automated and Manual Gating of Perforin/CD8',cache=FALSE,fig.align='center',fig.width=7.5,fig.height=3.5,out.width=500---- p1<-plotGate(auto_gating[[sampleNames(gating_subset)[7]]],"cd8/Prf",default.y="",xbin=256,margin=FALSE,path=2,arrange=FALSE,main="Automated",xlab=list(cex=1.5),ylab=list(cex=1.5),scales=list(cex=1.5),par.strip.text=list(cex=1.5))[[1]] p2<-plotGate(gating_subset[[sampleNames(gating_subset)[7]]],"8+/Perforin+",default.y="",xbin=256,margin=FALSE,path=2,xlab=list(cex=1.5),ylab=list(cex=1.5),scales=list(cex=1.5),par.strip.text=list(cex=1.5),arrange=FALSE,main="Manual")[[1]] grid.arrange(p1,p2,ncol=2) ## ----example_derivative_gate,echo=FALSE,eval=TRUE,fig.align='center',out.width=500,fig.width=7.5,fig.height=3.5---- fr<-getData(auto_gating[[7]],"cd8/cd57-") cut<-tailgate(fr,channel="",filter_id="perforin_gate") dens<-density(exprs(fr[,""]),adjust=2) second_deriv <- diff(sign(diff(dens$y))) which_maxima <- which(second_deriv == -2) + 1 which_maxima <- which_maxima[findInterval(dens$x[which_maxima], range(exprs(fr[,""]))) == 1] which_maxima <- which_maxima[order(dens$y[which_maxima], decreasing = TRUE)] peaks <- dens$x[which_maxima] deriv_out <- openCyto:::.deriv_density(x = exprs(fr[,""]), adjust = 2, deriv = 1) par(mfrow=c(1,2)) plot(dens,main="Density") abline(v=peaks) abline(v=cut@min,col="red") plot(deriv_out,type="l",main="First Derivative",xlab="X",ylab="Density") abline(v=cut@min,col="red") abline(v=peaks) ## ----example_TNFa,echo=FALSE,eval=TRUE,cache=FALSE,fig.cap="Automated and manual gating of TFNa",fig.align='center'---- p1<-plotGate(auto_gating[[sampleNames(gating_subset)[7]]],c("cd8/TNFa","cd4/TNFa"),default.y="",xbin=256,margin=FALSE,path=2,arrange=FALSE,main="Automated",xlab=list(cex=1.1),ylab=list(cex=1.1),scales=list(cex=1.1),par.strip.text=list(cex=1.1)) p2<-plotGate(gating_subset[[sampleNames(gating_subset)[7]]],c("8+/TNFa+","4+/TNFa+"),default.y="",xbin=256,margin=FALSE,path=2,xlab=list(cex=1.1),ylab=list(cex=1.1),scales=list(cex=1.1),par.strip.text=list(cex=1.1),arrange=FALSE,main="Manual") do.call(grid.arrange,c(p1,p2,ncol=2)) ## ----getData,echo=TRUE,eval=FALSE---------------------------------------- ## cd3_population<-getData(auto_gating,"cd3") ## ----plotGate,echo=TRUE,eval=FALSE--------------------------------------- ## plotGate(auto_gating,"cd3") ## ----subset,echo=TRUE,eval=FALSE----------------------------------------- ## first_ten_fcs_files<-auto_gating[1:10] ## ----listgtMethods,echo=TRUE,eval=FALSE---------------------------------- ## listgtMethods() ## ----register_plugin,echo=TRUE,eval=FALSE-------------------------------- ## registerPlugins(myfunction,methodName,dependencies, "preprocessing"|"gating") ## ----gen_template,echo=TRUE,eval=FALSE,results='asis'-------------------- ## templateGen(gating_subset[[1]]) ## ----gen_template_noecho,echo=FALSE,eval=TRUE,results='asis'------------- hd<-head(data.frame(templateGen(gating_subset[[1]])),7) hd[is.na(hd)]<-"" Kmisc::htmlTable(hd) ## ----echo=TRUE,eval=FALSE------------------------------------------------ ## # Work with a subset of the data ## auto_subset<-auto_gating[1:20] # first 20 samples ## # Make changes in the template. SAVE it, then re-load it in R. ## gt<-gatingTemplate("data/template/gt_080.csv") ## # Remove the old gate from auto_gating ## Rm("viable",auto_subset) ## # Gate with the new template, stop after the viable gate (which is nonNeutro) ## gating(gt,auto_subset,stop.at="nonNeutro") ## # View the new result ## plotGate(auto_subset,"viable")