R/center_data.R

In hyochoi/SCISSOR: Shape changes in selecting sample outliers in RNA-seq

#' Centering data
#'
#' This function finds a robust mean vector along the gene length using a trimmed mean. Also, it provides
#' a new data object
#'
#' @export
center_data = function(inputData, Ranges, average="mean",trim=0.1,adjval=NULL, ...) {

  d = dim(inputData)[1]
  nexons = nrow(Ranges$lRanges)
  mean.vec = apply(inputData,1,FUN=function(x){adj.center(x,average=average,trim=trim,adjval=adjval)})

  exonbaseMat = matrix(0,ncol=nexons,nrow=d)
  colnames(exonbaseMat) = rownames(Ranges$lRanges)
  zeromean = c()
  for (ie in 1:nexons) {
    x = mean.vec[Ranges$lRanges[ie,2]:Ranges$lRanges[ie,3]]
    if (sum(x^2)>1e-05) {
      exonbaseMat[Ranges$lRanges[ie,2]:Ranges$lRanges[ie,3],ie] = x/sum(x^2)
    } else {
      zeromean = c(zeromean,ie)
    }
  }
  msf.exons = t(exonbaseMat)%*%inputData
  keepexons = which(! c(1:nrow(msf.exons)) %in% zeromean)
  if (length(keepexons)>1) {
    msf = apply(msf.exons[which(! c(1:nrow(msf.exons)) %in% zeromean),],2,median)
  } else if (length(keepexons)==1) {
    msf = msf.exons[which(! c(1:nrow(msf.exons)) %in% zeromean),]
  } else {
    msf = rep(0,dim(inputData)[2])
  }
  NewX = inputData - mean.vec%*%t(msf)
  # msf2=rep(1,length(msf))
  # msf2[which(msf>1)]=1/msf[which(msf>1)]
  # NewX = sweep(NewX,2,msf2,"*")
  return(list(outputData=NewX,msf=msf,data.center=mean.vec))
}

center_data_v1 = function(inputData, robustLM=FALSE, ...) {
  msfres = scale.factor(X=inputData, robustLM=robustLM, ...)
  data.center = msfres$mean.vec
  msf=msfres$msf
  msf2=rep(1,length(msf))
  msf2[which(msf>1)]=1/msf[which(msf>1)]
  outdata = sweep(msfres$NewX,2,msf2,"*")
  return(list(outputData=outdata,msf=msf,data.center=data.center))
}