R/ebseq.R

In Coraline66/RNASeqAnalysis: Perform RNASeq Data Analysis

#' Function "ebseq"
#'
#' This function Perform Differential Gene Expression Analysis using EBSeq.
#' @import EBSeq
#' @param type a character type variable indicating whether to perform "Gene-level analysis" or "Isoform-level analysis".
#' @param data either an unnormalized counts matrix for "Gene-level analysis", or a list for "Isoform-level analysis".
#' @param condition a length two variable containing two conditions.
#' @param n1 an integer indicating the number of samples in the first condition.
#' @param n2 an integer indicating the number of samples in the second condition.
#' @param N an integer indicating the number of iterations in EB Test.
#' @param fdr a numeric variable indicating the target FDR in EB Test.
#' @param method "robust" or "classic". Refer to GetDEResults() in EBSeq. Default is "robust".
#' @param fdrmethod "hard" or "soft". Refer to GetDEResults() in EBSeq. Default is "hard".
#' @param fc Threshold for the fold change statistics.Default is 0.7.
#' @param Qcut Transcrips with 100 th quantile <= Qcut will be removed before testing. Default is 0.
#' @param tt The number of uncertainty groups the user wish to define. Default is 3.
#' @return a list consisting of the output and results from EB Test.
#' @export
#'

ebseq <- function(type, data, condition, n1, n2, N=5, fdr=0.05, method="robust", fdrmethod="hard", fc=0.7, Qcut=0, tt=3) {
  if (type == "Gene") {
    ## get library size factor
    sizes <- MedianNorm(data)
    groups <- factor(c(rep(condition[1], n1), rep(condition[2], n2)), levels = condition)
    
    ## EB test
    EBout <- EBTest(Data = data, Conditions = groups, sizeFactors = sizes, maxround = N, QtrmCut = Qcut)
    
    ## calculate the fold change
    GeneFC <- PostFC(EBout)
    
    ## extract result of EB test
    EBresult <- GetDEResults(EBout, FDR = fdr, Method = method, FDRMethod = fdrmethod, Threshold_FC = fc)
    
    ## return results
    return(list("FC" = GeneFC, "output" = EBout, "result" = EBresult))
  }
  
  if (type =="Isoform") {
    groups <- factor(c(rep(condition[1], n1), rep(condition[2], n2)), levels = condition)
    IsoMat <- data$IsoMat
    IsoNames <- data$IsoNames
    IsoGeneNames <- data$IsoGeneNames
    
    ## get library size factor
    IsoSizes <- MedianNorm(IsoMat)
    
    ## get the I_g vector
    NgList <- GetNg(IsoformName = IsoNames, GeneName = IsoGeneNames, TrunThre = tt)
    
    ## EB test
    IsoEBout <- EBTest(Data = IsoMat, NgVector = NgList$IsoformNgTrun, Conditions = groups, sizeFactors = IsoSizes, maxround = N, QtrmCut = Qcut)
    
    ## calculate fold change
    IsoFC <- PostFC(IsoEBout)
    
    ## extract result of EB test
    IsoEBresult <- GetDEResults(IsoEBout, FDR = fdr, Method = method, FDRMethod = fdrmethod, Threshold_FC = fc)
    
    ## return results
    return(list("Nglist" = NgList, "FC" = IsoFC, "output" = IsoEBout, "result" = IsoEBresult))
  }
}