R/qsmooth.R

In stephaniehicks/qsmooth: Smooth quantile normalization

#' @title qsmooth
#' 
#' @description This function applies a generalization of 
#' quantile normalization called smoothed quantile 
#' normalization. This function defines the qsmooth class 
#' and constructor. 
#'
#' @param object an object which is a \code{matrix} or 
#' \code{data.frame} with observations (e.g. probes or genes) on 
#' the rows and samples as the columns. Alternatively, 
#' a user can provide a \code{SummarizedExperiment} object
#' and the \code{assay(object, "counts")} will be used as input 
#' for the qsmooth normalization. 
#' @param group_factor a group level continuous or categorial 
#' covariate associated with each sample or column in the 
#' \code{object}. The order of the \code{group_factor} must 
#' match the order of the columns in \code{object}.
#' @param batch (Optional) batch covariate (multiple 
#' batches are not allowed). 
#' If batch covariate is provided, \code{Combat()} from 
#' \code{sva} is used prior to qsmooth normalization to 
#' remove batch effects. See \code{Combat()} for more details. 
#' @param norm_factors optional normalization scaling factors.
#' @param window window size for running median which is 
#' a fraction of the number of rows in \code{object}. 
#' Default is 0.05. 
#' @export 
#' 
#' @return A object of the class \code{qsmooth} that 
#' contains a numeric vector of the qsmooth weights in 
#' the \code{qsmoothWeights} slot and a matrix of normalized  
#' values after applying smoothed quantile normalization in 
#' the \code{qsmoothData} slot. 
#' 
#' @details 
#' Quantile normalization is one of the most widely used 
#' normalization tools for data analysis in genomics. Although it 
#' was originally developed for gene expression microarrays it is
#' now used across many different high-throughput applications 
#' including RNAseq and ChIPseq. The methodology 
#' relies on the assumption that observed changes in the empirical 
#' distribution of samples are due to unwanted variability. 
#' Because the data is transformed to remove these differences 
#' it has the potential to remove interesting biologically driven
#' global variation. Therefore, applying quantile normalization, 
#' or other global normalization methods that rely on similar 
#' assumptions, may not be an appropriate depending on the type 
#' and source of variation. 
#' 
#' This function computes a weight at every 
#' quantile that compares the variability between groups 
#' relative to within groups. In one extreme quantile 
#' normalization is applied and in the other extreme quantile 
#' normalization within each biological condition is applied. 
#' The weight shrinks the group-level quantile normalized data 
#' towards the overall reference quantiles if variability 
#' between groups is sufficiently smaller than the variability 
#' within groups. 
# 
#' See the vignette for more details. 
#' 
#' @aliases qsmooth
#' 
#' @docType methods
#' @name qsmooth
#' @importFrom sva ComBat
#' @importFrom SummarizedExperiment assays
#' @importFrom SummarizedExperiment assay
#' @importFrom stats ave
#' 
#' @examples
#' dat <- cbind(matrix(rnorm(1000), nrow=100, ncol=10), 
#'              matrix(rnorm(1000, .1, .7), nrow=100, ncol=10))
#' dat_qs <- qsmooth(object = dat, 
#'                   group_factor = rep(c(0,1), each=10))
#' 
#' @rdname qsmooth
#' @export
qsmooth <- function(object, group_factor,
                    batch = NULL, norm_factors = NULL, 
                    window = 0.05)
{

  if(!any(is(object, "matrix") | is(object, "data.frame") |
     is(object, "SummarizedExperiment"))){
    stop("The class of the object must be a matrix, 
         data.frame or SummarizedExperiment")
  }
  
  if(is.data.frame(object)){ object <- as.matrix(object) }
  
  if(is(object, "SummarizedExperiment")){ 
    if("counts" %in% names(assays(object))){ 
      object <- assay(object, "counts")
    } else {
      stop("There is no assay slot named 'counts' inside 
           the object. Please check the names of the 
           assay slots using names(assays(object)).")
    }
  }
  
  if(is.null(group_factor)){  
    stop("Must provide group_factor to specify the group 
        level information associated with each sample or 
             or column in object.")
  }
  
  if(ncol(object) != length(group_factor)){
    stop("Number of columns in object does not match 
             length of group_factor.")
  }
  
  if(is.factor(group_factor) & length(levels(group_factor)) < 2){
    stop(paste0("group_factor is a factor and number of levels in 
                  group_factor is less than 2 (levels(group_factor): ",
                levels(group_factor), "). Must provide a factor with 
                  2 or more levels to use qsmooth."))
  }
  
  if(any(is.na(object))){ stop("Object must not contains NAs.") }    
  
  # Scale normalization step
  if(!is.null(norm_factors)) { 
    object <- sweep(object, 2, norm_factors, FUN = "/") 
  }
  
  # If batch is provided, run sva::ComBat to remove batch effects
  if(!is.null(batch)){
    object <- ComBat(object, factor(batch))
  }
  
  # Compute quantile statistics
  qs <- qstats(object=object, group_factor=group_factor,
               window=window)
  Qref <- qs$Qref 
  Qhat <- qs$Qhat  
  w <- qs$smoothWeights
  
  # Weighted quantiles
  objectNorm = w * Qref + (1 - w) * Qhat
  
  # Re-order objectNorm by rank of object (columnwise)
  for (i in seq_len(ncol(objectNorm))) {
    # Grab ref. i
    ref = objectNorm[,i]
    
    # Grab object column i
    x = object[,i]
    
    # Grab ranks of x (using min rank for ties)
    rmin = rank(x, ties.method="min")
    
    # If x has rank ties then average the values of ref at 
    # those ranks
    dups = duplicated(rmin)
    
    if (any(dups)) {
      # Grab ranks of x (using random ranks for ties) 
      # (needed to uniquely identify the indices of tied ranks)
      rrand = rank(x, ties.method="random")
      
      # Grab tied ranks
      tied.ranks = unique(rmin[dups])
      
      for (k in tied.ranks) {
        # Select the indices of tied ranks 
        sel = rrand[rmin == k] 
        ref[sel] = ave(ref[sel])
      }
    }
    
    # Re-order ref and replace in objectNorm
    objectNorm[,i] = ref[rmin]
  }
  
  results <- new("qsmooth")
  results@qsmoothWeights <- w
  
  rownames(objectNorm) = rownames(object)
  colnames(objectNorm) = colnames(object)
  results@qsmoothData <- objectNorm
  
  return(results)
}