R/difftest.R

In TSCAN: Tools for Single-Cell Analysis

#' difftest
#' 
#' testing differentially expressed genes
#'
#' This function tests whether a gene is significantly expressed given pseudotime ordering. Likelihood ratio test is performed to compare a generalized additive model (GAM) with a constant fit to get the p-values. The p-values are adjusted for multiple testing by fdr.
#'
#' @param data The raw single_cell data, which is a numeric matrix or data.frame. Rows represent genes/features and columns represent single cells.
#' @param TSCANorder The TSCAN ordering generated by function \code{\link{TSCANorder}}.
#' @param df Numeric value specifying the degree of freedom used in the GAM model.
#' @return Data frame containing pvalues and qvalues of testing differentially expression.
#' @export
#' @import mgcv
#' @author Zhicheng Ji, Hongkai Ji <zji4@@zji4.edu>
#' @examples
#' data(lpsdata)
#' procdata <- preprocess(lpsdata)
#' lpsorder <- TSCANorder(exprmclust(procdata))
#' diffval <- difftest(procdata,lpsorder)
#' #Selected differentially expressed genes under qvlue cutoff of 0.05
#' row.names(diffval)[diffval$qval < 0.05]

difftest <- function(data, TSCANorder, df = 3) {   
      ptime <- 1:length(TSCANorder)
      pval <- apply(data[,TSCANorder], 1, function(x) {
            if (sum(x) == 0) {
                  1
            } else {
                  model <- mgcv::gam(x~s(ptime,k=3))
                  pchisq(model$null.deviance - model$deviance, model$df.null - model$df.residual,lower.tail = F)                  
            }             
      })      
      qval <- p.adjust(pval, method = "fdr")
      data.frame(pval = pval,qval = qval)
}