R/deGenesANOVA.R

In maigesPack: Functions to handle cDNA microarray data, including several methods of data analysis

## Function to fit linear models using limma functions
##
## Parameters: data    -> object of maigesANOVA class
##             eBayes  -> logical indicating the use of empirical Bayes
##                        statistics from limma package
##             retOrig -> return original objects generated by limma. Defaults
##                        to FALSE
##             retF    -> Must return the results associated with the F test
##                        (TRUE) or with the individual contrasts (FALSE -
##                        default)?
##             doClust -> Logical. The object generated from this analysis will
##                        be used for cluster analysis? Defaults to TRUE
##             ...     -> additional parameters for functions lmFit (from
##                        limma)
##
## Gustavo Esteves
## 14/05/07
##
##


deGenesANOVA <- function(data=NULL, eBayes=FALSE, retOrig=FALSE, retF=FALSE,
doClust=TRUE, ...) {
    
    
    ## Making some basic tests
    if(is.null(data))
        stop("The data object must be specified (class maiges)!")
    
    
    ## Picking R and packages version information
    tmp <- sessionInfo()
    vInfo <- list()
    vInfo$R.version <- tmp$R.version$version.string
    vInfo$BasePacks <- tmp$basePkgs
    tmp1 <- NULL
    for (i in 1:length(tmp$otherPkgs))
        tmp1 <- c(tmp1, paste(tmp$otherPkgs[[i]]$Package, "version",
        tmp$otherPkgs[[i]]$Version))
    
    vInfo$AddPacks <- tmp1
    
    
    ## Defining the table to be used to do the tests
    wTable <- calcW(data)
    
    
    ## Converting the object to MAList class
    toFit <- as(data, "MAList")
    
    
    
    ## indexing samples where Ref was labelled with 'red'
    idxMult <- tolower(getLabels(data, "Ref")) == "red"
    if (sum(idxMult) > 0) 
        toFit$M[, idxMult] <- (-1)*toFit$M[, idxMult]
    
    
    ##
    ## Doing the tests
    
    ## Adjusting the linear model by limma functions
    fitTmp <- lmFit(toFit, data@Dmatrix, weights=matrix(1, dim(data)[1],
    dim(data)[2]), ...)
    
    
    ## Doing the tests...
    if(eBayes) {
        fit <- eBayes(contrasts.fit(fitTmp, data@Cmatrix))
        if(retF)
            test <- "ANOVA - moderated F test"
        else
            test <- "ANOVA - moderated t tests"
    }
    else {
        fit <- contrastsFitM(fitTmp, data@Cmatrix)
        if(retF)
            test <- "ANOVA - F test"
        else
            test <- "ANOVA - t tests"
    }
    
    if(retOrig)
        return(fit)
    else if(retF) {
        if(doClust)
            result <- new("maigesDEcluster")
        else
            result <- new("maigesDE")
        result@GeneInfo <- data@Glabels
        result@SampleInfo <- data@Slabels
        if(doClust)
            result@W <- wTable
        result@stat <- as.matrix(fit$F)
        result@p.value <- as.matrix(fit$F.p.value)
        result@factors <- colnames(data@Dmatrix)
        result@test <- test
        result@Date <- date()
        result@V.info <- vInfo
        
        for(i in 1:length(result@GeneInfo))
            result@GeneInfo[[i]][data@BadSpots] <-
            paste(result@GeneInfo[[i]][data@BadSpots], "(*)")
        
        return(result)
    }
    else {
        if(doClust)
            result <- new("maigesDEcluster", GeneInfo=data@Glabels,
            SampleInfo=data@Slabels, W=wTable, fold=fit$coefficients,
            stat=as.matrix(fit$t), p.value=as.matrix(fit$p.value),
            factors=colnames(data@Dmatrix), test=test, Date=date(),
            V.info=vInfo)
        else
            result <- new("maigesDE", GeneInfo=data@Glabels,
            SampleInfo=data@Slabels, fold=fit$coefficients,
            stat=as.matrix(fit$t), p.value=as.matrix(fit$p.value),
            factors=colnames(data@Dmatrix), test=test, Date=date(),
            V.info=vInfo)
        
        for(i in 1:length(result@GeneInfo))
            result@GeneInfo[[i]][data@BadSpots] =
            paste(result@GeneInfo[[i]][data@BadSpots], "(*)")
        
        return(result)
        
    }
    
}