R/calcVIFArm.R

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

#' A function to calculate the Variance Inflation Factor (VIF) for each of the gene sets input in geneSets This function does not depend on which method you used to calculate the variance originally.  If you assumed pooled variance, then use qusage calcVIF function for the variance will be calculated using LIMMA's "interGeneCorrelation" method (see Wu and Smyth, Nucleic Acids Res. 2012). This method will calculate the VIF separately for each group and return the average of each group's vif not using limma, not equal variances is assumed.
#' @param eset a matrix of log2(expression values). This must be the same dataset that was used to create geneResults
#' @param geneResults A QSarray object, as generated by either makeComparison or aggregateGeneSet
#' @param useAllData Boolean parameter determining whether to use all data in eset to calculated the VIF, or whether to only use data from the groups being contrasted. Only used if useCAMERA==FALSE
#' @import limma 
#' @export 
#' @return variance inflation factor 
calcVIFArm = function(eset,         ##a matrix of log2(expression values). This must be the same dataset that was used to create geneResults
                   geneResults,  ##A QSarray object, as generated by either makeComparison or aggregateGeneSet
                  useAllData = TRUE ##Boolean parameter determining whether to use all data in eset to calculated the VIF, or whether to only use data from the groups being contrasted. Only used if useCAMERA==FALSE
                ){
    if(is.null(geneResults$pathways)){stop("Pathway Information not found. Please provide a list of gene sets.")}
    geneSets = geneResults$pathways #numeric vector

  if(is(eset, "ExpressionSet")){eset = exprs(eset)}
  if(class(geneResults) != "QSarray"){stop("geneResults must be a QSarray object, as created by makeComparison")}
  ##run VIF calculation on each gene set
  vif = sapply(names(geneSets),function(i){
    GNames<-names(geneResults$mean)[geneSets[[i]]]
    gs.i = which(rownames(eset)%in%GNames)
    if(length(gs.i)<2){warning("GeneSet '",i,"' contains one or zero overlapping genes. NAs produced.");return(NA)}
     if(!is.null(geneResults$sd.alpha)){
       #call calcVIFarm.cpp with sdalpha return the list, and set attrbts
       t<-calcVIFarmalt(names(geneResults$mean),gs.i, geneResults$pathways[[i]],rownames(eset),eset,levels(geneResults$labels), geneResults$sd.alpha) #dispatches for each i in geneResults$pathways list
         }
         if(is.null(geneResults$sd.alpha)){
         #call calcVIFarm.cpp without sdalpha  return the list , set attrbts
        t<-calcVIFarm_nosdalphaalt(names(geneResults$mean),gs.i, geneResults$pathways[[i]],rownames(eset),eset,levels(geneResults$labels))
        }
  return(as.vector(t$vif))
   })
  geneResults$vif = vif
  if(!is.null(geneResults$path.PDF)){ ##if defined, rescale the pdf with the new vif values
    geneResults$path.PDF = t(t(geneResults$path.PDF) / pdfScaleFactor(geneResults))
  }
  return(geneResults)
}