R/normalize_counts_extra.R

In leandroroser/RNASeqFunctions: Several functions for rna-seq analysis

#' Normalization for counts matrix x
#' @description Normalization for counts matrix
#' @param x matrix to normalize
#' @param method method used for normalization
#' @param features_length length of features
#' @author Leandro Roser \email{leandroroser@@gmail.com}
#' @export

normalize_counts_extra <- function(x,
                             method = c("cpm", "rpkm", "quantile"),
                             features_length = NULL) {


  cpm_ <- function(x) {
    library_size <- colSums(x)
    t(apply(x, 1, function(y) 1E6 * y / library_size))
  }

  rpkm_ <- function(x) {
    library_size <- colSums(x)
    t(apply(x, 1, function(y) 1E9 * y  / (library_size*feature_length)))
  }


  quantile_ <- function(x) {

    # function to obtain a matrix with indices for ordering x by-column
    ord <- function(y, dec) apply(y, 2, function(z) order(z, decreasing = dec))

    x <- as.matrix(x)
    xcols <- ncol(x)
    xrows <- nrow(x)
    orden <- ord(x, dec = TRUE)
    x.ord <- sapply(seq_len(xcols), function(y)x[, y][orden[, y]])
    x.ord.mean <- apply(x.ord, 1, mean, na.rm = TRUE)
    out <- matrix(rep(x.ord.mean, xcols), nrow = xrows)

    # reorder the indices matrix with the original positions
    orden.inv <- ord(orden, dec = FALSE)
    x.ord <- sapply(seq_len(xcols), function(y)out[, y][orden.inv[, y]])
    x.ord
  }

  if(method == "cpm") {
    out <- cpm_(x)
  } else if(method == "rpkm") {
    out <- rpkm(x, feature_length)
  } else if(method == "quantile") {
    out <- quantile_(x)
  }
  rownames(out) <- rownames(x)
  colnames(out) <- colnames(x)
  out
}