### R code from vignette source 'qusage.Rnw' ################################################### ### code chunk number 1: install (eval = FALSE) ################################################### ## if (!requireNamespace("BiocManager", quietly=TRUE)) ## install.packages("BiocManager") ## BiocManager::install("qusage") ################################################### ### code chunk number 2: setup ################################################### options(width=70) options(SweaveHooks=list(fig=function() par(mar=c(5.1, 4.1, 1.1, 2.1)))) library(qusage) library(nlme) data("fluExample") data("GeneSets") i = c(which(flu.meta$Hours==0), which(flu.meta$Hours==77)) eset = eset.full[,i] resp = flu.meta$Condition[i] ################################################### ### code chunk number 3: exprSimple ################################################### dim(eset) eset[1:5,1:5] ################################################### ### code chunk number 4: labelsSimple ################################################### labels = c(rep("t0",17),rep("t1",17)) labels ################################################### ### code chunk number 5: contrastSimple ################################################### contrast = "t1-t0" ################################################### ### code chunk number 6: geneSetSimple ################################################### ISG.geneSet[1:10] ################################################### ### code chunk number 7: qsSimple ################################################### qs.results = qusage(eset, labels, contrast, ISG.geneSet) ################################################### ### code chunk number 8: pValSimple ################################################### pdf.pVal(qs.results) ################################################### ### code chunk number 9: plotCodeSimple (eval = FALSE) ################################################### ## plot(qs.results, xlab="Gene Set Activation") ################################################### ### code chunk number 10: plotSimple ################################################### getOption("SweaveHooks")[["fig"]]() plot(qs.results, xlab="Gene Set Activation") ################################################### ### code chunk number 11: read.MSIG.geneSets (eval = FALSE) ################################################### ## MSIG.geneSets = read.gmt("c2.cp.kegg.v4.0.symbols.gmt") ################################################### ### code chunk number 12: MSIG.geneSets ################################################### summary(MSIG.geneSets[1:5]) MSIG.geneSets[2] qs.results.msig = qusage(eset, labels, contrast, MSIG.geneSets) numPathways(qs.results.msig) ################################################### ### code chunk number 13: multiHypCorr ################################################### p.vals = pdf.pVal(qs.results.msig) head(p.vals) q.vals = p.adjust(p.vals, method="fdr") head(q.vals) ################################################### ### code chunk number 14: qusage.Rnw:150-151 ################################################### options(width=100) ################################################### ### code chunk number 15: qsTable ################################################### qsTable(qs.results.msig, number=10) ################################################### ### code chunk number 16: qusage.Rnw:158-159 ################################################### options(width=70) ################################################### ### code chunk number 17: plotCIs ################################################### getOption("SweaveHooks")[["fig"]]() par(mar=c(10,4,1,2)) plot(qs.results.msig) ################################################### ### code chunk number 18: pairVector ################################################### pairs = c(1:17,1:17) pairs ################################################### ### code chunk number 19: twoWayLabels ################################################### resp twoWay.labels = paste(resp, labels, sep=".") twoWay.labels ################################################### ### code chunk number 20: twoWayRun ################################################### sx.results = qusage(eset, twoWay.labels, "sx.t1-sx.t0", MSIG.geneSets, pairVector=pairs) asx.results = qusage(eset, twoWay.labels, "asx.t1-asx.t0", MSIG.geneSets, pairVector=pairs) ################################################### ### code chunk number 21: twoCurvePval ################################################### p.vals = twoCurve.pVal(sx.results,asx.results) head(p.vals) ################################################### ### code chunk number 22: twoWayComplex ################################################### asx.vs.sx = qusage(eset, twoWay.labels, "(sx.t1-sx.t0) - (asx.t1-asx.t0)", MSIG.geneSets, pairVector=pairs) p.vals = pdf.pVal(asx.vs.sx) head(p.vals) ################################################### ### code chunk number 23: twoWayDensityCurvesCode ################################################### plotDensityCurves(asx.results,path.index=125,xlim=c(-0.2,0.5), col=3,main="CYTOSOLIC DNA SENSING") plotDensityCurves(sx.results,path.index=125,col=2,add=T) plotDensityCurves(asx.vs.sx, path.index=125,col=1,add=T, lwd=3) legend("topleft", legend=c("asx","sx","asx-sx"),lty=1,col=c(2:3,1), lwd=c(1,1,3)) ################################################### ### code chunk number 24: twoWayDensityCurves ################################################### getOption("SweaveHooks")[["fig"]]() plotDensityCurves(asx.results,path.index=125,xlim=c(-0.2,0.5), col=3,main="CYTOSOLIC DNA SENSING") plotDensityCurves(sx.results,path.index=125,col=2,add=T) plotDensityCurves(asx.vs.sx, path.index=125,col=1,add=T, lwd=3) legend("topleft", legend=c("asx","sx","asx-sx"),lty=1,col=c(2:3,1), lwd=c(1,1,3)) ################################################### ### code chunk number 25: plotGSD ################################################### getOption("SweaveHooks")[["fig"]]() plotGeneSetDistributions(sx.results,asx.results,path.index=125) ################################################### ### code chunk number 26: plotCIsGenes ################################################### getOption("SweaveHooks")[["fig"]]() par(mar=c(8,4,1,2)) plotCIsGenes(sx.results,path.index=125) ################################################### ### code chunk number 27: qgen_fluexpr ################################################### dim(eset.full) eset.full[1:5,1:5] ################################################### ### code chunk number 28: qgen_design ################################################### head(flu.meta) ################################################### ### code chunk number 29: qgen_fixedSimple ################################################### fixedeffects <- ~Hours+Condition+Condition*Hours+Age+Gender ################################################### ### code chunk number 30: qgen_contrast.factor ################################################### con.fac <- ~Condition*Hours ################################################### ### code chunk number 31: qgen_contrast2 ################################################### contrast = "sx77 - sx0" ################################################### ### code chunk number 32: qgen_geneSet ################################################### ISG.geneSet[1:10] ################################################### ### code chunk number 33: qgen_randomeffect ################################################### qusage.result.random<-qgen(eset.full,flu.meta, fixed=fixedeffects, contrast.factor=con.fac, contrast=contrast, geneSets=ISG.geneSet, random=~1|Subject, design.sampleid="SampleID") ################################################### ### code chunk number 34: qgen_randomResult ################################################### qsTable(qusage.result.random) ################################################### ### code chunk number 35: qgen_correlation (eval = FALSE) ################################################### ## corStruct = corExp(value=3,~Hours.Numeric|Subject) ## corStruct ## ## qusage.result.corstruct<-qgen(eset.full,flu.meta, ## fixed=fixedeffects, ## contrast.factor=con.fac, ## contrast=contrast, ## geneSets=ISG.geneSet, ## correlation=corStruct, ## design.sampleid="SampleID") ################################################### ### code chunk number 36: qgen_plot (eval = FALSE) ################################################### ## plotDensityCurves(qusage.result.random, main="ISG_geneSet", ## xlim=c(0,1.2), xlab="Day 77 - Day 0", col="red") ## plotDensityCurves(qusage.result.corstruct, col="blue",add=T) ################################################### ### code chunk number 37: plotQgen (eval = FALSE) ################################################### getOption("SweaveHooks")[["fig"]]() ## plotDensityCurves(qusage.result.random, main="ISG_geneSet", ## xlim=c(0,1.2), xlab="Day 77 - Day 0", col="red") ## plotDensityCurves(qusage.result.corstruct, col="blue",add=T) ################################################### ### code chunk number 38: LoadMetaSets ################################################### data("fluVaccine") ################################################### ### code chunk number 39: BTM (eval = FALSE) ################################################### ## read.gmt("BTM_for_GSEA_20131008.gmt") ################################################### ### code chunk number 40: GSE59635_desc ################################################### head(fluVaccine$phenoData[["GSE59635"]]) ################################################### ### code chunk number 41: GSE59654_desc ################################################### head(fluVaccine$phenoData[["GSE59654"]]) ################################################### ### code chunk number 42: metaAnalysis_Qusage ################################################### qs.results.list = lapply(c("GSE59635", "GSE59654"), function(n){ eset = fluVaccine$esets[[n]] phenoData = fluVaccine$phenoData[[n]] labels = paste(phenoData$responder, phenoData$bloodDrawDay, sep = ".") ##run QuSAGE qusage(eset, labels, "(high.d7-high.d0)-(low.d7-low.d0)", BTM.geneSets, pairVector=phenoData$subjectId, n.points=2^14) }) ################################################### ### code chunk number 43: metaAnalysis_combinePDFs ################################################### combined = combinePDFs(qs.results.list) ################################################### ### code chunk number 44: metaAnalysis_pval ################################################### head(pdf.pVal(combined)) ################################################### ### code chunk number 45: metaAnalysis_plot ################################################### plotCombinedPDF(combined, path.index=235) legend("topleft", legend=c("GSE59635", "GSE59654", "Meta-analysis"), lty=1, col=c("red", "blue","black")) ################################################### ### code chunk number 46: plotMeta ################################################### getOption("SweaveHooks")[["fig"]]() plotCombinedPDF(combined, path.index=235) legend("topleft", legend=c("GSE59635", "GSE59654", "Meta-analysis"), lty=1, col=c("red", "blue","black")) ################################################### ### code chunk number 47: prettyPDF ################################################### ##select the best 5 gene sets p.vals = pdf.pVal(qs.results.msig) topSets = order(p.vals)[1:5] ##plot plotDensityCurves(qs.results.msig, path.index=topSets, col=1:5, lwd=2, xlab=expression(log[2](Activity))) ################################################### ### code chunk number 48: plotPDF ################################################### getOption("SweaveHooks")[["fig"]]() ##select the best 5 gene sets p.vals = pdf.pVal(qs.results.msig) topSets = order(p.vals)[1:5] ##plot plotDensityCurves(qs.results.msig, path.index=topSets, col=1:5, lwd=2, xlab=expression(log[2](Activity))) ################################################### ### code chunk number 49: createCDF ################################################### ##get coordinates coord.list = plotDensityCurves(qs.results, plot=FALSE) coords = coord.list[[1]] x=coords$x ##Calculate the CDF as the integral of the PDF at each point. y = cumsum(coords$y)/sum(coords$y) plot(x,y, type="l", xlim=calcBayesCI(qs.results,low=0.001)[,1]) ################################################### ### code chunk number 50: plotCDF ################################################### getOption("SweaveHooks")[["fig"]]() ##get coordinates coord.list = plotDensityCurves(qs.results, plot=FALSE) coords = coord.list[[1]] x=coords$x ##Calculate the CDF as the integral of the PDF at each point. y = cumsum(coords$y)/sum(coords$y) plot(x,y, type="l", xlim=calcBayesCI(qs.results,low=0.001)[,1]) ################################################### ### code chunk number 51: plotISG_GSD ################################################### getOption("SweaveHooks")[["fig"]]() plotGeneSetDistributions(qs.results) ################################################### ### code chunk number 52: plotGSDHighlight ################################################### getOption("SweaveHooks")[["fig"]]() path = 125 ##the CYTOSOLIC_DNA_SENSING pathway path.gene.i = sx.results$pathways[[path]] ##the genes are stored as indexes path.genes = rownames(eset)[path.gene.i] cols = rep("#33333399",length(path.genes)) ##make the IFNA genes red, and the IFNB genes blue. cols[path.genes=="CXCL10"] = "#FF0000" plotGeneSetDistributions(sx.results,asx.results,path.index=125, groupLabel=c("Symptomatic","Asymptomatic"), colorScheme=cols, barcode.col=cols)