examples/example.R

library(QUALIFIER)


#save_db(db, path = "/home/wjiang2/rglab/workspace/QUALIFIER/output/preprocessedData", overwrite = T, cdf = "copy")
db <- load_db(path = "/home/wjiang2/rglab/workspace/QUALIFIER/output/preprocessedData")
qaTask.list <- db$qaTaskList



###############################################################################
#1.parse gating template
###############################################################################
ws<-openWorkspace("/shared/silo_researcher/Gottardo_R/mike_working/ITN029ST/QA_template.xml")
GT<-parseWorkspace(ws
					,name=2
					,execute=F
					,subset=1
					)
gh_template<-GT[[1]]					
gs_pop_get_count_fast(gh_template)[,2:3]

#test the matchNode

###############################################################################
#2.apply gating template to new data
###############################################################################
			
datapath<-"/shared/silo_researcher/Gottardo_R/mike_working/ITN029ST/"
#newSamples<-getSample(gh_template)
newSamples<-list.files(datapath)[1:200]
length(newSamples)
G<-GatingSet(gh_template
			,newSamples[1:100]
			,path=datapath
			,isNcdf=T
#			,dMode=4
			)
gs_pop_get_count_fast(G[[1]])[,2:3]
plotGate(G[[1]],merge=F)
################################################################################  
#3.extract stats (6 min for 500 samples, 12 min for 1k samples)
###############################################################################
#library(parallel)
db<-new.env()
initDB(db)
metaFile="/shared/silo_researcher/Gottardo_R/mike_working/ITN029ST/FCS_File_mapping.csv"

qaPreprocess(db=db,gs=G
		,metaFile=metaFile
		,fcs.colname="FCS_Files"
		,date.colname=c("RecdDt","AnalysisDt")
		,nslave=1
#		,type="SOCK"
		,isMFI=T,isSpike=T, isChannel = T
)
getQAStats(G[[1]],isMFI=T,isSpike=T)
#colnames(pData(db$gs[[1]]))
#pData(G)
#head(db$stats)
#getQAStats(G[1])

#saveToDB(db=db,gs=G
#		,metaFile=metaFile
#		,fcs.colname="FCS_Files"
#		,date.colname=c("RecdDt","AnalysisDt")
#	)
################################################################################  
#4.load QA check list
###############################################################################
checkListFile<-file.path(system.file("data",package="QUALIFIER"),"qaCheckList.csv.gz")
qaTask.list <- read.qaTask(db,checkListFile=checkListFile)


pData(db$gs[[1]])$RecdDt<-as.Date(pData(db$gs[[1]])$RecdDt)
QUALIFIER:::.parseTubeID(db)#parse TubeID from FCS filenames

#read pre-determined events number for tubes from csv file
##pannel name should be in place of tube name since the entire package is using pannel name 
##to represent the tube

tubesEvents<-read.csv(file.path(system.file("data",package="QUALIFIER"),"tubesevents.csv.gz"),row.names=1)
tubesEventsOrig<-QUALIFIER:::.TubeNameMapping(db,tubesEvents=tubesEvents[,3,drop=F])
tubesEvents20090825<-QUALIFIER:::.TubeNameMapping(db,tubesEvents=tubesEvents[,2,drop=F])
tubesEvents20090622<-QUALIFIER:::.TubeNameMapping(db,tubesEvents=tubesEvents[,1,drop=F])


###80% of the pre-defined the value for each pannel
qaCheck(qaTask.list[["NumberOfEvents"]]
		,formula=count ~ RecdDt | Tube
		,outlierfunc=list(func=outlier.cutoff,args=list(lBound=0.8*tubesEvents20090825))
		,subset=as.Date(RecdDt,"%m/%d/%y")>='2009-08-25'
)

qaCheck(qaTask.list[["NumberOfEvents"]]
		,formula=count ~ RecdDt | Tube
		,outlierfunc=list(func=outlier.cutoff,args=list(lBound=0.8*tubesEvents20090622))
		,subset=as.Date(RecdDt,"%m/%d/%y")<'2009-08-25'&as.Date(RecdDt,"%m/%d/%y")>='2009-06-22'
)
qaCheck(qaTask.list[["NumberOfEvents"]]
		,formula=count ~ RecdDt | Tube
		,outlierfunc=list(func=outlier.cutoff,args=list(lBound=0.8*tubesEventsOrig))
		,subset=as.Date(RecdDt,"%m/%d/%y")<'2009-06-22'
)

CairoX11()

plot(qaTask.list[["NumberOfEvents"]]
#		,count ~ RecdDt
		,subset=Tube=='CD8/CD25/CD4/CD3/CD62L'
#,dest="image"
#		,scales=list(x=list(rot=45
#							,cex=0.5
#							))
#		,pch=19
)

#qaTask.list[["NumberOfEvents"]]@formula<-count ~ `RecdDt` | `t/Tube`:`t/PTID`

plot(qaTask.list[["NumberOfEvents"]]
#		,subset=id=='245'
#		,scatterPlot=TRUE
#		,scatterPar=list(stat=T)
)
clearCheck(qaTask.list[["NumberOfEvents"]])


##add new aggregated stats
QUALIFIER:::.addStats(qaTask.list[["BoundaryEvents"]]
		,definition=sum(proportion)~RecdDt|fileid+gsid
		,pop="/root/MNC/margin"
		,statName="sum.prop"
        )

head(subset(
				queryStats(qaTask.list[["BoundaryEvents"]]
							,y=sum.prop ~ RecdDt 
							,pop="margin"
#							,subset=value>0&id==806
							)
				,outlier==TRUE)
)
#check on the new stats
qaCheck(qaTask.list[["BoundaryEvents"]]
		,sum.prop ~ RecdDt 
		,outlierfunc= list(func = outlier.cutoff,args=list(uBound=0.0003))
		)

plot(qaTask.list[["BoundaryEvents"]]
		,sum.prop ~ RecdDt 
#		,subset=channel=="PE-A"
#				&id==806
#		,ylab="sum(percent)"
#		,scatterPlot=T
#		,scatterPar=list(
#						xlog=T
#						,stat=T
#						)
		,scales=list(format="%m/%d/%y")
#		,plotAll="none"
#		,dest="image"
)

#plot(qaTask.list[["BoundaryEvents"]]
#		,proportion ~ RecdDt|channel 
##		,subset=channel=="PE-A"
##				&id==806
#		,ylab="percent"
##		,scatterPlot=T
##		,scatterPar=list(
##						xlog=T
##						,stat=T
##						)
#		,scales=list(format="%m/%d/%y")
###		,plotAll=F
##		,dest="image"
#)

clearCheck(qaTask.list[["BoundaryEvents"]])
		
		
		
		
qaCheck(qaTask.list[["MFIOverTime"]]
#		,outlierfunc=outlier.norm
		,rFunc=rlm
#		,Subset=channel%in%c('PE-Cy7-A')
#		,z.cutoff=10
)

head(subset(
        queryStats(qaTask.list[["MFIOverTime"]]
#            ,y=sum.prop ~ RecdDt 
            ,pop="MFI"
#							,subset=value>0&id==806
        )
        ,outlier==TRUE)
)
plot(qaTask.list[["MFIOverTime"]]
#		,y=MFI~RecdDt|stain
		,subset=channel%in%c('APC-A')
				&stain=="CD62L APC-A"
##				&id==806
#		,rFunc=rlm
#		,scales=list(format="%m/%d/%y")
##		,scatterPlot=TRUE
#		,scatterPar=list(xlog=F)
#		,dest="image"

)

clearCheck(qaTask.list[["MFIOverTime"]])


qaCheck(qaTask.list[["RBCLysis"]]
		,outlierfunc=list(func=outlier.cutoff
                        ,args=list(lBound=0.9))
)

subset(
		queryStats(qaTask.list[["RBCLysis"]]
					,subset=Tube=='CD8/CD25/CD4/CD3/CD62L')
		,outlier==TRUE)

plot(qaTask.list[["RBCLysis"]]
#		,subset=Tube=='CD8/CD25/CD4/CD3/CD62L'
#				&id%in%c(270)
#		, RecdDt~proportion | Tube
#		,scales=list(format="%m/%d/%y")
#		,ylab="percent"
#		,scatterPlot=T
#		,scatterPar=list(stat=T
#						,xbin=128)
#		,horiz=T
#		,dest="image"
#		,highlight="coresampleid"
#		,plotAll="none"
#		,width=27,height=13
)	

clearCheck(qaTask.list[["RBCLysis"]])


qaCheck(qaTask.list[["spike"]]
		,outlierfunc = list(func = outlier.t
                  		    ,args = list(alpha=0.001
                                        ,isLower=FALSE)
                          )
		
)
plot(qaTask.list[["spike"]]
		,y=spike~RecdDt|channel
		,subset=Tube=='CD11c/CD80/DUMP/HLADr/CD123'
#	,dest="image"
#	,plotAll=T
)

plot(qaTask.list[["spike"]],y=spike~RecdDt|channel
		,subset=channel=='FITC-A'
#					&id%in%c(245,119)
#		,scatterPlot=TRUE
		,scatterPar=list(ylog=T
						,xlim=c(0,100)
#						,xbin=128
						)
#		,dest="image"
#		,plotAll="none"
)


qaCheck(qaTask.list[["MNC"]]
#		,Subset=coresampleid%in%c(11730,8780)
#		,z.cutoff=0.1
)

plot(qaTask.list[["MNC"]]
		,proportion~`coresampleid`
#		, factor(coresampleid)~proportion
#		,par=list(horiz=TRUE)
)

#scatter plot for a sample group	
plot(qaTask.list[["MNC"]]
		,proportion~factor(coresampleid)
#		,par=list(xlab="coresampleid")
#		, coresampleid ~proportion
#		,par=list(horiz=TRUE)
#		,subset=coresampleid%in%c(
##									11730
#									8780
#		)
#		,scatterPlot=TRUE
		,scatterPar=list(xbin=128
						,stat=T)
		,dest="~/rglab/workspace/QUALIFIER/output/image"
		,plotAll=F
)
#scatter for one sample
plot(qaTask.list[["MNC"]]
		,scatterPlot=TRUE
		,subset=coresampleid==11730&id==245)

clearCheck(qaTask.list[["MNC"]])

qaCheck(qaTask.list[["RedundantStain"]]
#			,gOutlierfunc=qoutlier
#			,outlierfunc=qoutlier
#			,alpha=1.5
#			,z.cutoff=2
)

##example of passing lattice arguments		
plot(qaTask.list[["RedundantStain"]]
		,subset=channel=='APC-A'
#				&stain%in%c('CD123','CD3')
#				&coresampleid==11730
		,y=proportion~factor(coresampleid)|channel:stain
#		,scatterPlot=T
		,scatterPar=list(xlog=F
						,stat=T
						)
		,scales=list(x=list(relation="free"))
#		,layout=c(2,NA,1)
#		,dest="image"
)
################################################################################  
#4.qa report in html format 
#set plotAll=TRUE to generate the scatter plots for all the individual FCS files 
#otherwise only plots for outliers are generated.
###############################################################################

##customerize some of the task before pass them to report method
htmlReport(qaTask.list[["MFIOverTime"]])<-TRUE
rFunc(qaTask.list[["MFIOverTime"]])<-rlm

highlight(qaTask.list[["BoundaryEvents"]])<-"coresampleid"
#scatterPar(qaTask.list[["MFIOverTime"]])<-list(xlog=TRUE)
#scatterPar(qaTask.list[["BoundaryEvents"]])<-list(xlog=TRUE)
#scatterPar(qaTask.list[["RedundantStain"]])<-list(xlog=TRUE)
qpar(qaTask.list[["RedundantStain"]])<-list(scales=list(x=list(relation="free")))

#modify functions within package namespace
funcToinsert <- ".postProcessSVG" 
funcSym <- as.symbol(funcToinsert)
eval(substitute(environment(ff) <- getNamespace("QUALIFIER"), list(ff = funcSym)))
assignInNamespace(funcToinsert, eval(funcSym), ns = "QUALIFIER")


db <- load_db(path = "/home/wjiang2/rglab/workspace/QUALIFIER/output/preprocessedData")
qaTask.list <- db$qaTaskList
db$stats[sid == 3, value:= 0.1]
qaCheck(qaTask.list[["MNC"]] 
        , gOutlierfunc=list(func=outlier.norm
                            ,args=list(alpha=0.1))
        )
         

plot(qaTask.list[["MNC"]]
#    , dest = "rglab/workspace/QUALIFIER/output/image"
#    , plotAll = "none"
#    , scatterPlot = T
#    , subset = participantid == "002005"
#    , scatterPar = list(xbin = 32)
    )
  
    
qaReport(qaTask.list, outDir = "rglab/workspace/QUALIFIER/output/"
          , plotAll = T
            , subset = participantid == "002005"
          )

Try the QUALIFIER package in your browser

Any scripts or data that you put into this service are public.

QUALIFIER documentation built on Oct. 31, 2019, 3:24 a.m.