R/AllClasses.R

In aghozlane/DESeq2shaman: Differential gene expression analysis based on the negative binomial distribution

#' @rdname DESeqDataSet
#' @export
#setClass("SummarizedExperiment")
setClass("DESeqDataSet",
         contains = "SummarizedExperiment",
         representation = representation( 
           design = "formula",
           dispersionFunction = "function"))

setValidity( "DESeqDataSet", function( object ) {
  if (! ("counts" %in% names(assays(object))))
    return( "the assays slot must contain a matrix named 'counts'" )
  if ( !is.numeric( counts(object) ) )
    return( "the count data is not numeric" )
  if ( any( is.na( counts(object) ) ) )
    return( "NA values are not allowed in the count matrix" )
  if ( !is.integer( counts(object) ) )
    return( "the count data is not in integer mode" )
  if ( any( counts(object) < 0 ) )
    return( "the count data contains negative values" )
  design <- design(object)
  designVars <- all.vars(design)
  if (!all(designVars %in% names(colData(object)))) {
    return("all variables in design formula must be columns in colData")
  }
  designVarsClass <- sapply(designVars, function(v) class(colData(object)[[v]]))
  if (any(designVarsClass == "character")) {
    return("variables in design formula are character vectors.
  convert these columns of colData(object) to factors before including in the design formula")
  }
  designFactors <- designVars[designVarsClass == "factor"]
  if (any(sapply(designFactors,function(v) any(table(colData(object)[[v]]) == 0)))) {
    return("factors in design formula must have samples for each level.
  this error can arise when subsetting a DESeqDataSet, in which
  all the samples for one or more levels of a factor in the design were removed.
  if this was intentional, use droplevels() to remove these levels, e.g.:

  dds$condition <- droplevels(dds$condition)
")
  }
  TRUE
} )

#' DESeqDataSet object and constructors
#'
#' The \code{DESeqDataSet} is a subclass of \code{SummarizedExperiment},
#' used to store the input values, intermediate calculations and results of an
#' analysis of differential expression.  The \code{DESeqDataSet} class
#' enforces non-negative integer values in the "counts" matrix stored as
#' the first element in the assay list.
#' In addition, a formula which specifies the design of the experiment must be provided.
#' The constructor functions create a DESeqDataSet object
#' from various types of input:
#' a SummarizedExperiment, a matrix, or count files generated by
#' the python package HTSeq.  See the vignette for examples of construction
#' from all three input types.
#'
#' @param se a \code{SummarizedExperiment} with at least one column in colData,
#' and the counts as the first element in the assays list, which will be renamed
#' "counts".  A \code{SummarizedExperiment} object can be generated by the
#' function \code{summarizeOverlaps} in the GenomicRanges package.
#' @param design a formula which specifies the design of the experiment, taking the form
#' \code{formula(~ x + y + z)}.  By default, the functions in this package will use
#' the last variable in the formula (e.g. z) for presenting results (fold changes, etc.) and plotting.
#' @param countData for matrix input: a matrix of non-negative integers
#' @param colData for matrix input: a \code{DataFrame} or \code{data.frame} with at least a single column.
#' Rows of colData correspond to columns of countData.
#' @param sampleTable for htseq-count: a \code{data.frame} with three or more columns. Each row
#' describes one sample. The first column is the sample name, the second column
#' the file name of the count file generated by htseq-count, and the remaining
#' columns are sample metadata which will be stored in \code{colData}
#' @param directory for htseq-count: the directory relative to which the filenames are specified
#' @param ignoreRank for advanced use only, allows creation of a
#' DESeqDataSet which is not of full rank
#' @param ... arguments provided to \code{SummarizedExperiment} including rowData and exptData
#' 
#' @return A DESeqDataSet object.
#' 
#' @aliases DESeqDataSet DESeqDataSet-class DESeqDataSetFromMatrix DESeqDataSetFromHTSeqCount
#'
#' @references See \url{http://www-huber.embl.de/users/anders/HTSeq} for htseq-count
#'
#' @docType class
#'
#' @examples
#'
#' countData <- matrix(1:4,ncol=2)
#' colData <- data.frame(condition=factor(c("a","b")))
#' dds <- DESeqDataSetFromMatrix(countData, colData, formula(~ condition))
#'
#' @rdname DESeqDataSet
#' @export
DESeqDataSet <- function(se, design, ignoreRank=FALSE) {
  if (is.null(names(assays(se))) || names(assays(se))[1] != "counts") {
    message("renaming the first element in assays to 'counts'")
    names(assays(se, withDimnames=FALSE))[1] <- "counts"
  }
  # before validity check, try to convert assay to integer mode
  if (any(assay(se) < 0)) {
    stop("some values in assay are negative")
  }
  if (!is.integer(assay(se))) {
    if (any(round(assay(se)) != assay(se))) {
      stop("some values in assay are not integers")
    }
    message("converting counts to integer mode")
    mode(assay(se)) <- "integer"
  }

  if (all(assay(se) == 0)) {
    message("all samples have 0 counts for all genes. check the counting script.")
  }
  
  if (all(rowSums(assay(se) == assay(se)[,1]) == ncol(se))) {
    message("all genes have equal values for all samples. will not be able to perform differential analysis")
  }

  if (any(duplicated(rownames(se)))) {
    warning(sum(duplicated(rownames(se)))," duplicate rownames were renamed by adding numbers")
    rnms <- rownames(se)
    dups <- unique(rnms[duplicated(rnms)])
    for (rn in dups) {
      idx <- which(rnms == rn)
      rnms[idx[-1]] <- paste(rnms[idx[-1]], c(seq_len(length(idx) - 1)), sep=".")
    }
    rownames(se) <- rnms
  }
  
  designVars <- all.vars(design)
  if (!all(designVars %in% names(colData(se)))) {
    stop("all variables in design formula must be columns in colData")
  }

  designVarsClass <- sapply(designVars, function(v) class(colData(se)[[v]]))
  if (any(designVarsClass == "character")) {
    warning("some variables in design formula are characters, converting to factors")
    for (v in designVars[designVarsClass == "character"]) {
      colData(se)[[v]] <- factor(colData(se)[[v]])
    }
  }

  designVarsNumeric <- sapply(designVars, function(v) is.numeric(colData(se)[[v]]))
  if (any(designVarsNumeric)) {
    warnIntVars <- FALSE
    for (v in designVars[designVarsNumeric]) {
      if (all(colData(se)[[v]] == round(colData(se)[[v]]))) {
        warnIntVars <- TRUE
      }
    }
    if (warnIntVars) {
      message(paste0("the design formula contains a numeric variable with integer values,
  specifying a model with increasing fold change for higher values.
  did you mean for this to be a factor? if so, first convert
  this variable to a factor using the factor() function"))
    }
  }

  designFactors <- designVars[designVarsClass == "factor"]
  missingLevels <- sapply(designFactors,function(v) any(table(colData(se)[[v]]) == 0))
  if (any(missingLevels)) {
    message("factor levels were dropped which had no samples")
    for (v in designFactors[missingLevels]) {
      colData(se)[[v]] <- droplevels(colData(se)[[v]])
    }
  }
  
  modelMatrix <- model.matrix(design, data=as.data.frame(colData(se)))
  if (!ignoreRank) {
    if (qr(modelMatrix)$rank < ncol(modelMatrix)) {
      stop("the model matrix is not full rank, so the model cannot be fit as specified.
  one or more variables or interaction terms in the design formula
  are linear combinations of the others and must be removed")
    }
  }

  # if the last variable in the design formula is a
  # factor, and has a level 'control', check if it is
  # the base level and if not print a message
  lastDV <- length(designVars)
  if (length(designVars) > 0 && designVarsClass[lastDV] == "factor") {
    lastDVLvls <- levels(colData(se)[[designVars[lastDV]]])
    controlSynonyms <- c("control","Control","CONTROL")
    for (cSyn in controlSynonyms) {
      if (cSyn %in% lastDVLvls) {
        if (cSyn != lastDVLvls[1]) {
          message(paste0("it appears that the last variable in the design formula, '",designVars[lastDV],"',
  has a factor level, '",cSyn,"', which is not the base level. we recommend
  to use factor(...,levels=...) or relevel() to set this as the base level
  before proceeding. for more information, please see the 'Note on factor levels'
  in vignette('DESeq2')."))
        }
      }
    }
  }
  
  # Add columns on the columns
  mcolsCols <- DataFrame(type=rep("input",ncol(colData(se))),
                         description=rep("",ncol(colData(se))))
  mcols(colData(se)) <- if (is.null(mcols(colData(se)))) {
    mcolsCols
  } else if (all(names(mcols(colData(se))) == c("type","description"))) {
    mcolsCols
  } else {
    cbind(mcols(colData(se)), mcolsCols)
  }
  dds <- new("DESeqDataSet", se, design = design)
                                 
  # now we know we have at least an empty GRanges or GRangesList for rowData
  # so we can create a metadata column 'type' for the mcols
  # and we label any incoming columns as 'input'

  # this is metadata columns on the rows
  mcolsRows <- DataFrame(type=rep("input",ncol(mcols(dds))),
                         description=rep("",ncol(mcols(dds))))
  mcols(mcols(dds)) <- if (is.null(mcols(mcols(dds)))) {
    mcolsRows
  } else if (all(names(mcols(mcols(dds))) == c("type","description"))) {
    mcolsRows
  } else {
    cbind(mcols(mcols(dds)), mcolsRows)
  }
  
  return(dds)
}

#' @rdname DESeqDataSet
#' @export
DESeqDataSetFromMatrix <- function( countData, colData, design, ignoreRank=FALSE, ... )
{
  # we expect a matrix of counts, which are non-negative integers
  countData <- as.matrix( countData )

  if (is(colData,"data.frame")) colData <- DataFrame(colData, row.names=rownames(colData))

  # check if the rownames of colData are simply in different order
  # than the colnames of the countData, if so throw an error
  # as the user probably should investigate what's wrong
  if (!is.null(rownames(colData)) & !is.null(colnames(countData))) {
    if (all(sort(rownames(colData)) == sort(colnames(countData)))) {
      if (!all(rownames(colData) == colnames(countData))) {
        stop(paste("rownames of the colData:
  ",paste(rownames(colData),collapse=","),"
  are not in the same order as the colnames of the countData:
  ",paste(colnames(countData),collapse=",")))
      }
    }
  }
  if (is.null(rownames(colData)) & !is.null(colnames(countData))) {
    rownames(colData) <- colnames(countData)
  }
  
  se <- SummarizedExperiment(assays = SimpleList(counts=countData), colData = colData, ...)
  dds <- DESeqDataSet(se, design = design, ignoreRank)
  return(dds)
}

#' @rdname DESeqDataSet
#' @export
DESeqDataSetFromHTSeqCount <- function( sampleTable, directory="", design, ignoreRank=FALSE, ...) 
{
  if (missing(design)) {
    stop("design is missing")
  }
  l <- lapply( as.character( sampleTable[,2] ), function(fn) 
              read.table( file.path( directory, fn ) ) )
  if( ! all( sapply( l, function(a) all( a$V1 == l[[1]]$V1 ) ) ) )
    stop( "Gene IDs (first column) differ between files." )
  tbl <- sapply( l, function(a) a$V2 )
  colnames(tbl) <- sampleTable[,1]
  rownames(tbl) <- l[[1]]$V1
  rownames(sampleTable) <- sampleTable[,1]
  oldSpecialNames <- c( "no_feature", "ambiguous",
                       "too_low_aQual", "not_aligned",
                       "alignment_not_unique" )
  # either starts with two underscores
  # or is one of the old special names (htseq-count backward compatability)
  specialRows <- (substr(rownames(tbl),1,1) == "_") | rownames(tbl) %in% oldSpecialNames
  tbl <- tbl[ !specialRows, , drop=FALSE ]
  dds <- DESeqDataSetFromMatrix(countData = tbl,
                                colData = sampleTable[,-(1:2),drop=FALSE],
                                design = design,
                                ignoreRank, ...)
  return(dds)
}   

#' @rdname DESeqResults
#' @export
setClass("DESeqResults", contains="DataFrame")

#' DESeqResults object and constructor
#'
#' This class extends the DataFrame class of the IRanges package
#' simply to allow other packages to write methods for results
#' objects from the DESeq2 package.
#'
#' @param DataFrame a DataFrame of results, standard column names are:
#' baseMean, log2FoldChange, lfcSE, stat, pvalue, padj.
#'
#' @return a DESeqResults object
#'
#' @docType class
#'
#' @aliases DESeqResults-class
#' 
#' @rdname DESeqResults
#' @export
DESeqResults <- function(DataFrame) {
  new("DESeqResults", DataFrame)
}