In arcolombo/junk: qusage: Quantitative Set Analysis for Gene Expression

SpeedSage Intro

qusage is published software that is slow for large runs, SpeedSage corrects for speed and efficiency at large orders. there is Bottlenecking of Functions Qusage can improve the speed of its algorithm by minimizing the cost of computaiton.

changes Armadillo C++

trading NA flexibility slows down qusage runs, but having the user input no NAs enforcing good input, this speeds up calcIndividualExpressions, as well as using C++ libraries.

calcVif Function

This test the local version which enforces no NA in Baseline or PostTreatment object, this reduces the flexibility. this test data is from the vignette where postTreatment was modified to be Baseline+20.4, a simple training set from the QuSAGE vignette.

library(inline)
library(microbenchmark)
library(Rcpp)
library(parallel)
library(speedSage)
library(qusage)
library(ggplot2)
eset<-system.file("extdata","eset.RData",package="speedSage")
load(eset)
labels<-c(rep("t0",134),rep("t1",134))
contrast<-"t1-t0"
colnames(eset)<-c(rep("t0",134),rep("t1",134))
fileISG<-system.file("extdata","c2.cgp.v5.1.symbols.gmt",package="speedSage")
ISG.geneSet<-read.gmt(fileISG)
geneSets<-ISG.geneSet[grepl("DER_IFN_GAMMA_RESPONSE_UP",names(ISG.geneSet))]
Baseline<-eset
PostTreatment<-eset+20.4
ncol(Baseline) #not splitting up eset
i<-1  #for testing 1 gene set 
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/sigmaArm.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/sigmaSingle.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/bayesEstimation.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/notbayesEstimation.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/calcVIFarm.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/calcVIFarmalt.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/calcVIFarm_nosdalphaalt.cpp")

pairVector<-NULL
var.equal<-FALSE
filter.genes<-FALSE
n.points<-2^12


#setting up calcVif call objects
results = makeComparisonArm(eset, labels, contrast, pairVector=pairVector,var.equal=var.equal)
nu = floor(min(results$dof,na.rm=T))
geneResults = aggregateGeneSet(results, geneSets, silent=F, n.points=n.points)
#eset, and results parameters for vif

useAllData<-TRUE
useCAMERA<-FALSE
##### qusage VIF calc
ogVIF<-calcVIF(eset,geneResults)

#for local code use  
ogGeneSets<-geneResults$pathways
GNames<-names(geneResults$mean)[ogGeneSets[[1]]]
gs.i = which(rownames(eset)%in%GNames)
###now to compare

library(speedSage)
useCAMERA<-FALSE
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/calcVIFarmalt.cpp")
sourceCpp(file="/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/calcVIFarm_nosdalphaalt.cpp")


t2<-calcVIFarmalt(names(geneResults$mean),gs.i, geneResults$pathways[[i]],rownames(eset),eset,levels(geneResults$labels), geneResults$sd.alpha) 



library(ggplot2)
speedUp<-microbenchmark(
calcVIFarmalt(names(geneResults$mean),gs.i, geneResults$pathways[[1]],rownames(eset),eset,levels(geneResults$labels), geneResults$sd.alpha),
calcVIF(eset,geneResults) )
speedUp

source("/home/anthonycolombo/Documents/qusage/qusage_repos/qusage_speed/R/calcVIFArm.R")
myVIF<-calcVIFArm(eset,geneResults)



### TO DO ensure that calcVIFarm matches qusage::calcVIF.R

identical(names(myVIF),names(ogVIF))
identical(myVIF$labels,ogVIF$labels)
 identical(myVIF$contrast,ogVIF$contrast)
identical(myVIF$n.samples,ogVIF$n.samples)
identical(myVIF$mean,ogVIF$mean)
 identical(myVIF$sd.alpha,ogVIF$sd.alpha)
 identical(myVIF$dof,ogVIF$dof)
identical(myVIF$var.method,ogVIF$var.method)
identical(myVIF$pathways,ogVIF$pathways)
identical(myVIF$path.mean,ogVIF$path.mean)
 identical(myVIF$path.size,ogVIF$path.size)
identical(myVIF$ranges,ogVIF$ranges)
 identical(myVIF$n.points,ogVIF$n.points)
all.equal(myVIF$path.PDF,ogVIF$path.PDF)
all.equal(myVIF$vif,ogVIF$vif)


qplot( as.vector(abs(myVIF$path.PDF-ogVIF$path.PDF)), xlab="Path PDF ") 

qplot( as.vector(abs(myVIF$vif-ogVIF$vif)), xlab="VIF ")

library(ggplot2)
mb<-microbenchmark(
myVIF<-calcVIFArm(eset,geneResults),
ogVIF<-calcVIF(eset,geneResults),times=1000 )

autoplot(mb)

mb