R/adjustExprLevel.R

In scottzijiezhang/m6Amonster: Analyze m6A seq data

#' @title adjustExprLevel
#' @param x The data list generated by normlizeLibrary() function.
#' @param adjustBy By default, adjust post-IP count by INPUT geneSum. Can also choose "pos" to use current position count to adjust for expression level.
#' @param mode Running "bin" mode or "peak" mode.
#' @return  x The data list now with IP-count adjusted for expression level added
#' @export
adjustExprLevel <- function(x,adjustBy = "geneSum",mode){

  if(mode == "bin"){
    gene.size <- t( apply(x$geneSum,1,function(x){x/mean(x)}) )
    gene.size.factor <- gene.size[x$geneBins$gene,]

    if(adjustBy == "geneSum"){
      ip_norm_geneSum <- x$norm.ip/gene.size.factor
      ip_norm_geneSum <- round(ip_norm_geneSum)
      na.flag <- apply(ip_norm_geneSum,1, function(x){any(is.na(x)) | any(is.infinite(x))})
      ip_norm_geneSum <-  ip_norm_geneSum[!na.flag,]
      x <- c(x,list('ip_adjExpr'=ip_norm_geneSum))
    }else if(adjustBy == "pos"){
      pos.size <-  t( apply(x$norm.input,1,function(x){x/mean(x)}) )
      ip_norm_pos <- x$norm.ip/pos.size
      ip_norm_pos <- round(ip_norm_pos)
      na.flag <- apply(ip_norm_pos,1, function(x){any(is.na(x)) | any(is.infinite(x))})
      ip_norm_pos <-  ip_norm_pos[!na.flag,]
      x <- c(x,list('ip_adjExpr'=ip_norm_pos))
    }else{
      stop("Must adjust by \"geneSum\" or by \"pos\"...")
    }
    return(x)
  }else if(mode == "peak"){

    if(adjustBy == "geneSum"){
      geneSum <- x$geneSum
      geneSum <- t(apply(geneSum,1,.noZero))
      gene.size <- t( apply(geneSum,1,function(x){x/mean(x)}) )
      gene.size.factor <- gene.size[x$geneBins[rownames(x$norm.jointPeak_ip),"gene"],]
      ip_norm_geneSum <- x$norm.jointPeak_ip/gene.size.factor
      ip_norm_geneSum <- round(ip_norm_geneSum)
      x <- c(x,list('jointPeak_adjExpr'=ip_norm_geneSum))
    }else if(adjustBy == "pos"){
      norm.jointPeak_input <-t( t(x$jointPeak_input) / x$sizeFactor$input )
      norm.jointPeak_input <- t(apply(norm.jointPeak_input,1,.noZero))
      pos.size <-  t( apply(norm.jointPeak_input,1,function(x){x/mean(x)}) )
      ip_norm_pos <- x$norm.jointPeak_ip/pos.size
      ip_norm_pos <- round(ip_norm_pos)
      x <- c(x,list('jointPeak_adjExpr'=ip_norm_pos))
    }else{
      stop("Must adjust by \"geneSum\" or by \"pos\"...")
    }
    return(x)
  }else{
    stop("Must choose the mode, either bin mode or peak mode.")
  }

}


## a helper function to remove 0 from vector
.noZero <- function(x){sapply(x,max,1)}