R/ZIQML.fit.R

In qianli10000/lncDIFF: lncRNA Normalized Counts Differential Expression Analysis

#' Group effect estimate on lncRNA normalized counts
#'
#' ZIQML.fit estimates the group effect on genes using zero-inflated exponential quasi likelihood.  
#' @param edata Normalized counts matrix with genes in rows and samples in columns.
#' @param design.matrix Design matrix for groups, generated by model.matrix().
#' @param link Link function for the generalized linear model and likelihood function,either 'log' or 'identity'. The default is 'log'. 
#' @return 
#'    \item{Estimates}{Estimated group effect on input genes by zero-inflated exponential quasi maximum likelihood (ZIQML) estimator.}
#'    \item{logLikelihood}{The value of zero-inflated quasi likelihood at estimated group effect.} 
#'    \item{edata}{The input gene expression matrix.} 
#'    \item{design.matrix}{The input design matrix.} 
#'    \item{link}{The specified link function.} 
#' @examples 
#' 
#' data('tcga.hnsc.match.edata','design') 
#' # 'tcga.hnsc.match.edata' contains RPKM of 1132 lncRNA genes and 80 samples 
#' # 'design' is the design matrix for tumor vs normal tissue.
#'                                        
#' fit.log=ZIQML.fit(edata=tcga.hnsc.match.edata,design.matrix=design,link='log') 
#' # Fit GLM by ZIQML with logarithmic link function 
#' 
#' fit.identity=ZIQML.fit(edata=tcga.hnsc.match.edata,design.matrix=design,link='identity') 
#' # Fit GLM by ZIQML with identity link function
#' @export 
#' @importFrom stats optim p.adjust pchisq 

ZIQML.fit=function(edata,design.matrix,link='log'){
  
  if(any(is.na(edata))) stop('Invalid values in expression data: NA')
  if(any(edata<0)) stop('Invalid values in expression data: Negative')
  if(any(is.na(design.matrix))) stop('Design matrix contains NA')
  if(class(design.matrix) %in% c('numeric','factor','character')) stop('Design matrix not in correct form')
  if(nrow(design.matrix)!=ncol(edata)) stop('Dimensions of input matrices do not match')
  if(sum(design.matrix[,1]==0)>0)stop('Design matrix must contain intercept')
  if(!link %in% c('log','identity'))stop('Link function is not valid')
  
  est=exact.llh.value=NULL
  edata=as.matrix(edata)
  design.matrix=as.matrix(design.matrix)
  for(i in 1:nrow(edata)){
    ini=c(mean(edata[i,]),rep(0,(ncol(design.matrix)-1)))#,sum(edata[i,]!=0)/length(edata[i,]))
    
    if(link=='identity'){
    mfit=optim(par = ini,
               fn = function(par)llh.identity(c(par,sum(edata[i,]!=0)/length(edata[i,])),edata[i,],design.matrix))
    }else{ 
      mfit=optim(par = c(log(ini[1]),ini[-1]),
                 fn = function(par)llh.log(c(par,sum(edata[i,]!=0)/length(edata[i,])),edata[i,],design.matrix))
    }
  est=rbind(est,mfit$par)
  exact.llh.value=c(exact.llh.value,-mfit$value)
  }
  rownames(est)=rownames(edata)
  colnames(est)=colnames(design.matrix)
  
  output=list(Estimates=est,logLikelihood=exact.llh.value,edata=edata,design.matrix=design.matrix,link=link)
  return(output)  
}