R/AllClasses.R

In mikelove/DESeq2: Differential gene expression analysis based on the negative binomial distribution

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

setValidity("DESeqDataSet", function(object) {
  if (! ("counts" %in% assayNames(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)
  # 'design' is either a formula or matrix
  stopifnot(is(design, "formula") | is(design, "matrix"))

  if (is(design, "formula")) {
    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"]
    # levels would duplicate after make.names()
    if (any(sapply(designFactors,function(v) {
      factor.lvls <- levels(colData(object)[[v]])
      factor.nms <- make.names(factor.lvls)
      any(duplicated(factor.nms))
    }))) {
      return("levels of factors in the design have non-unique level names after make.names() is applied.
  best to only use letters and numbers for levels of factors in the design")
    }
    # levels contain characters other than letters, numbers, and underscore
    if (any(sapply(designFactors,function(v) {
      factor.lvls <- levels(colData(object)[[v]])
      any(!grepl("^[A-Za-z0-9_.]+$",factor.lvls))
    }))) {
      # just a warning for now
      message("  Note: levels of factors in the design contain characters other than
  letters, numbers, '_' and '.'. It is recommended (but not required) to use
  only letters, numbers, and delimiters '_' or '.', as these are safe characters
  for column names in R. [This is a message, not a warning or an error]")
    }
  } else if (is(design, "matrix")) {
    # TODO add some more tests for if 'design' is matrix
    stopifnot(nrow(design) == ncol(object))
  }
  
  TRUE
})
  
#' DESeqDataSet object and constructors
#'
#' \code{DESeqDataSet} is a subclass of \code{RangedSummarizedExperiment},
#' 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 RangedSummarizedExperiment, a matrix, count files generated by
#' the python package HTSeq, or a list from the tximport function in the
#' tximport package.
#' See the vignette for examples of construction from different types.
#'
#' Note on the error message "assay colnames() must be NULL or equal colData rownames()":
#' this means that the colnames of countData are different than the rownames of colData.
#' Fix this with: \code{colnames(countData) <- NULL}
#'
#' @param se a \code{RangedSummarizedExperiment} with columns of variables
#' indicating sample information in \code{colData},
#' and the counts as the first element in the assays list, which will
#' be renamed "counts". A \code{RangedSummarizedExperiment} object can be
#' generated by the function \code{summarizeOverlaps} in the GenomicAlignments
#' package.
#' @param design a \code{formula} or \code{matrix}.
#' the \code{formula} expresses how the counts for each gene
#' depend on the variables in \code{colData}. Many R \code{formula} are valid,
#' including designs with multiple variables, e.g., \code{~ group + condition},
#' and designs with interactions, e.g., \code{~ genotype + treatment + genotype:treatment}.
#' See \code{\link{results}} for a variety of designs and how to extract results tables.
#' By default, the functions in this package will use 
#' the last variable in the formula for building results tables and plotting.
#' \code{~ 1} can be used for no design, although users need to remember
#' to switch to another design for differential testing.
#' @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 tidy for matrix input: whether the first column of countData is the rownames for the count matrix
#' @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 txi for tximport: the simple list output of the \code{tximport} function
#' @param directory for htseq-count: the directory relative to which the filenames are specified. defaults to current directory
#' @param ignoreRank use of this argument is reserved for DEXSeq developers only.
#' Users will immediately encounter an error upon trying to estimate dispersion
#' using a design with a model matrix which is not full rank.
#' @param ... arguments provided to \code{SummarizedExperiment} including rowRanges and metadata. Note that
#' for Bioconductor 3.1, rowRanges must be a GRanges or GRangesList, with potential metadata columns
#' as a DataFrame accessed and stored with \code{mcols}. If a user wants to store metadata columns
#' about the rows of the countData, but does not have GRanges or GRangesList information,
#' first construct the DESeqDataSet without rowRanges and then add the DataFrame with \code{mcols(dds)}.
#' 
#' @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:100,ncol=4)
#' condition <- factor(c("A","A","B","B"))
#' dds <- DESeqDataSetFromMatrix(countData, DataFrame(condition), ~ condition)
#'
#' @rdname DESeqDataSet
#' @importFrom utils packageVersion
#' @export
DESeqDataSet <- function(se, design, ignoreRank=FALSE) {
  if (!is(se, "RangedSummarizedExperiment")) {
    if (is(se, "SummarizedExperiment")) {
      se <- as(se, "RangedSummarizedExperiment")
    } else {
      stop("'se' must be a RangedSummarizedExperiment object")
    }
  }
  if (is.null(assayNames(se)) || assayNames(se)[1] != "counts") {
    # check if "counts" are in another assay spot:
    if (length(assayNames(se)) > 1) {
      if ("counts" %in% assayNames(se)[-1]) {
        stop("'counts' assay is not the first in assays list, re-arrange so that 'counts' is first")
      }
    }
    message("renaming the first element in assays to 'counts'")
    assayNames(se)[1] <- "counts"
  }
  # special tximeta processing
  if ("tximetaInfo" %in% names(metadata(se))) {
    se <- processTximeta(se)
  }  
  # before validity check, try to convert assay to integer mode
  if (any(is.na(assay(se))))
    stop("NA values are not allowed in the count matrix")
  if (any(assay(se) < 0)) {
    stop("some values in assay are negative")
  }
  if (!is.integer(assay(se))) {
    if (!is.numeric(assay(se))) {
      stop(paste("counts matrix should be numeric, currently it has mode:", mode(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)) {
    stop("all samples have 0 counts for all genes. check the counting script.")
  }
  
  if (all(MatrixGenerics::rowSums(assay(se) == assay(se)[,1]) == ncol(se))) {
    warning("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
  }

  if (is(design, "formula")) {
    designVars <- all.vars(design)
    if (!all(designVars %in% names(colData(se)))) {
      stop("all variables in design formula must be columns in colData")
    }
    for (v in designVars) {
      if (any(is.na(colData(se)[[v]])))
        stop(paste("variables in design formula cannot contain NA:", v))
    }
    
    designVarsClass <- sapply(designVars, function(v) class(colData(se)[[v]])[1])

    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]])
      }
    }
    
    if (length(designVars) == 1) {
      var <- colData(se)[[designVars]]
      if (all(var == var[1])) {
        stop("design has a single variable, with all samples having the same value.
  use instead a design of '~ 1'. estimateSizeFactors, rlog and the VST can then be used")
      }
    }

    designVarsNumeric <- sapply(designVars, function(v) is.numeric(colData(se)[[v]]))
    if (any(designVarsNumeric)) {
      msgIntVars <- FALSE
      msgCenterScale <- FALSE
      for (v in designVars[designVarsNumeric]) {
        if (all(colData(se)[[v]] == round(colData(se)[[v]]))) {
          msgIntVars <- TRUE
        }
        if (mean(colData(se)[[v]]) > 5 | sd(colData(se)[[v]]) > 5) {
          msgCenterScale <- TRUE
        }
      }
      if (msgIntVars) {
        message("  the design formula contains one or more numeric variables 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")
      }
      if (msgCenterScale) {
        message("  the design formula contains one or more numeric variables that have mean or
  standard deviation larger than 5 (an arbitrary threshold to trigger this message).
  Including numeric variables with large mean can induce collinearity with the intercept.
  Users should center and scale numeric variables in the design to improve GLM convergence.")
      }
    }

    if (any(designVarsClass == "ordered")) {
      stop("the design formula contains an ordered factor. The internal steps
do not work on ordered factors as a formula. Instead you should provide a matrix to
the 'design' slot or to the 'full' argument of DESeq(), constructed using model.matrix.")
    }
    
    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]])
      }
    }
    
    singleLevel <- sapply(designFactors, function(v) all(colData(se)[[v]] == colData(se)[[v]][1]))
    if (any(singleLevel)) {
      stop("design contains one or more variables with all samples having the same value,
  remove these variables from the design")
    }
    
    # if the last variable in the design formula is a
    # factor, and has a level 'control', check if it is
    # the reference 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 reference level. we recommend
  to use factor(...,levels=...) or relevel() to set this as the reference level
  before proceeding. for more information, please see the 'Note on factor levels'
  in vignette('DESeq2')."))
          }
        }
      }
    }

    modelMatrix <- stats::model.matrix.default(design, data=as.data.frame(colData(se)))
  } else if (is(design, "matrix")) {
    modelMatrix <- design
  } else {
    stop("'design' should be a formula or a matrix")
  }

  if (!ignoreRank) {
    checkFullRank(modelMatrix)
  }
  
  # 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)
  }
  
  object <- new("DESeqDataSet", se, design = design)
                                 
  # now we know we have at least an empty GRanges or GRangesList for rowRanges
  # 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(object))),
                         description=rep("",ncol(mcols(object))))
  mcols(mcols(object)) <- if (is.null(mcols(mcols(object)))) {
    mcolsRows
  } else if (all(names(mcols(mcols(object))) == c("type","description"))) {
    mcolsRows
  } else {
    cbind(mcols(mcols(object)), mcolsRows)
  }

  # stash the package version
  metadata(object)[["version"]] <- packageVersion("DESeq2")
  
  return(object)
}

#' @rdname DESeqDataSet
#' @export
DESeqDataSetFromMatrix <- function( countData, colData, design, tidy=FALSE, ignoreRank=FALSE, ... )
{
  
  if (tidy) {
    stopifnot(ncol(countData) > 1)
    rownms <- as.character(countData[,1])
    countData <- countData[,-1,drop=FALSE]
    rownames(countData) <- rownms
  }

  # check that these agree in number
  stopifnot(ncol(countData) == nrow(colData))

  if (is(countData, "data.frame")) {
    if (any(sapply(countData, is, "factor"))) {
      warning("\n\n  'countData' is a data.frame with one or more factor columns.
  Be aware that converting directly to numeric can lead to errors.
  Provide matrices of integers to avoid this error.")
    }
  }
  
  # we expect a matrix of counts, which are non-negative integers
  countData <- as.matrix( countData )

  if (is(colData,"data.frame"))
    colData <- as(colData, "DataFrame")
  
  # 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, ...)
  object <- DESeqDataSet(se, design = design, ignoreRank)

  return(object)
}

#' @rdname DESeqDataSet
#' @export
DESeqDataSetFromHTSeqCount <- function( sampleTable, directory=".", design, ignoreRank=FALSE, ...) 
{
  sampleTable <- as.data.frame(sampleTable) # in case tibble or data.table are provided
  if (missing(design)) stop("design is missing")
  l <- lapply( as.character( sampleTable[,2] ), function(fn) read.table( file.path( directory, fn ), fill=TRUE ) )
  if( ! all( sapply( l, function(a) all( a$V1 == l[[1]]$V1 ) ) ) )
    stop( "Gene IDs (first column) differ between files." )
  # select last column of 'a', works even if htseq was run with '--additional-attr'
  tbl <- sapply( l, function(a) a[,ncol(a)] )
  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 ]
  object <- DESeqDataSetFromMatrix(countData=tbl,colData=sampleTable[,-(1:2),drop=FALSE],design=design,ignoreRank, ...)
  return(object)
}   

#' @rdname DESeqDataSet
#' @export
DESeqDataSetFromTximport <- function(txi, colData, design, ...) 
{
  stopifnot(is(txi, "list"))
  counts <- round(txi$counts)
  mode(counts) <- "integer"
  object <- DESeqDataSetFromMatrix(countData=counts, colData=colData, design=design, ...)
  stopifnot(txi$countsFromAbundance %in% c("no","scaledTPM","lengthScaledTPM"))
  if (txi$countsFromAbundance %in% c("scaledTPM","lengthScaledTPM")) {
    message("using just counts from tximport")
  } else {
    message("using counts and average transcript lengths from tximport")
    lengths <- txi$length
    stopifnot(all(lengths > 0))
    dimnames(lengths) <- dimnames(object)
    assays(object)[["avgTxLength"]] <- lengths
  }
  return(object)
}   


#' @rdname DESeqResults
#' @export
setClass("DESeqResults",
         contains="DFrame",
         representation = representation( 
           priorInfo = "list")
         )

#' DESeqResults object and constructor
#'
#' This constructor function would not typically be used by "end users".
#' This simple class indirectly extends the DataFrame class defined in the
#' S4Vectors package to allow other packages to write methods for results
#' objects from the DESeq2 package. It is used by \code{\link{results}}
#' to wrap up the results table.
#'
#' @param DataFrame a DataFrame of results, standard column names are:
#' baseMean, log2FoldChange, lfcSE, stat, pvalue, padj.
#' @param priorInfo a list giving information on the log fold change prior
#'
#' @return a DESeqResults object
#' @docType class
#' @aliases DESeqResults-class
#' @rdname DESeqResults
#' @export
DESeqResults <- function(DataFrame, priorInfo=list()) {
  new("DESeqResults", DataFrame, priorInfo=priorInfo)
}

#' @rdname DESeqTransform
#' @export
setClass("DESeqTransform", contains="RangedSummarizedExperiment")

#' DESeqTransform object and constructor
#'
#' This constructor function would not typically be used by "end users".
#' This simple class extends the RangedSummarizedExperiment class of the
#' SummarizedExperiment package.
#' It is used by \code{\link{rlog}} and
#' \code{\link{varianceStabilizingTransformation}}
#' to wrap up the results into a class for downstream methods,
#' such as \code{\link{plotPCA}}.
#' 
#' @param SummarizedExperiment a RangedSummarizedExperiment
#'
#' @return a DESeqTransform object
#' @docType class
#' @aliases DESeqTransform-class
#' @rdname DESeqTransform
#' @export
DESeqTransform <- function(SummarizedExperiment) {
  se <- SummarizedExperiment
  if (!is(se, "RangedSummarizedExperiment")) {
    if (is(se, "SummarizedExperiment")) {
      se <- as(se, "RangedSummarizedExperiment")
    } else {
      stop("'SummarizedExperiment' must be a RangedSummarizedExperiment object")
    }
  }
  new("DESeqTransform", se)
}

# unexported helper function
processTximeta <- function(se) {
  assay(se) <- round(assay(se))
  mode(assay(se)) <- "integer"
  stopifnot("countsFromAbundance" %in% names(metadata(se)))
  cfa <- metadata(se)$countsFromAbundance
  stopifnot(cfa %in% c("no","scaledTPM","lengthScaledTPM"))
  lengthPresent <- "length" %in% assayNames(se)
  if (cfa %in% c("scaledTPM","lengthScaledTPM")) {
    message("using just counts from tximeta")
  } else {
    if (lengthPresent) {
      message("using counts and average transcript lengths from tximeta")
      stopifnot(all(assays(se)$length > 0))
      anms <- assayNames(se)
      assayNames(se)[anms == "length"] <- "avgTxLength"
    } else {
      message("length assay is missing, using original counts without length adjustment")
    }
  }
  return(se)
}