R/runDEA.R

In xlucpu/MOVICS: Multi-Omics integration and VIsualization in Cancer Subtyping

#' @name twoclassdeseq2
#' @title Run two class comparison
#' @description Two class differential expression analysis using DESeq2 algorithm.
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param countsTable A matrix of RNA-Seq raw count data with rows for genes and columns for samples.
#' @param prefix A string value to indicate the prefix of output file.
#' @param overwt A logic value to indicate if to overwrite existing results; FALSE by default.
#' @param sort.p A logic value to indicate if to sort adjusted p value for output table; TRUE by default.
#' @param verbose A logic value to indicate if to only output id, log2fc, pvalue, and padj; TRUE by default.
#' @param res.path A string value to indicate the path for saving the results.
#' @importFrom DESeq2 DESeqDataSetFromMatrix DESeq results
#' @references Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15(12):550-558.
#' @export
#' @return Several .txt files storing differential expression analysis results by DESeq2
#' @keywords internal
#' @examples # There is no example and please refer to vignette.
twoclassdeseq2 <- function(moic.res    = NULL,
                           countsTable = NULL,
                           prefix      = NULL,
                           overwt      = FALSE,
                           sort.p      = TRUE,
                           verbose     = TRUE,
                           res.path    = getwd()) {

  # Create comparison list for differential expression analysis between two classes.
  createList  <- function(moic.res = NULL) {

    mo.method <- moic.res$mo.method
    moic.res <- moic.res$clust.res
    tumorsam <- moic.res$samID
    sampleList = list()
    treatsamList =list()
    treatnameList <- c()
    ctrlnameList <- c()

    n.moic <- length(unique(moic.res$clust))

    for (i in 1:n.moic) {
      sampleList[[i]] <- tumorsam
      treatsamList[[i]] = intersect(tumorsam, moic.res[which(moic.res$clust == i), "samID"])
      treatnameList[i] <- paste0("CS",i)
      ctrlnameList[i] <- "Others"
    }

    return(list(sampleList, treatsamList, treatnameList, ctrlnameList, mo.method))
  }

  complist <- createList(moic.res = moic.res)

  sampleList <- complist[[1]]
  treatsamList <- complist[[2]]
  treatnameList <- complist[[3]]
  ctrlnameList <- complist[[4]]
  mo.method <- complist[[5]]
  allsamples <- colnames(countsTable)

  options(warn = 1)
  for (k in 1:length(sampleList)) {
    samples <- sampleList[[k]]
    treatsam <- treatsamList[[k]]
    treatname <- treatnameList[k]
    ctrlname <- ctrlnameList[k]

    compname <- paste(treatname, "_vs_", ctrlname, sep = "")
    tmp <- rep("others", times = length(allsamples))
    names(tmp) <- allsamples
    tmp[samples] <- "control"
    tmp[treatsam] <- "treatment"

    if(!is.null(prefix)) {
      outfile <- file.path(res.path, paste(mo.method, "_", prefix, "_deseq2_test_result.", compname, ".txt", sep = ""))
    } else {
      outfile <- file.path(res.path, paste(mo.method ,"_deseq2_test_result.", compname, ".txt", sep = ""))
    }
    if (file.exists(outfile) & (overwt == FALSE)) {
      cat(paste0("deseq2 of ",compname, " exists and skipped...\n"))
      next
    }

    saminfo <- data.frame("Type" = as.factor(tmp[samples]),
                          "SampleID" = samples,
                          stringsAsFactors = FALSE)

    cts <- countsTable[,samples]
    coldata <- saminfo[samples,]

    dds <- quiet(DESeq2::DESeqDataSetFromMatrix(countData = cts,
                                                colData = coldata,
                                                design = as.formula("~ Type")))

    dds$Type <- relevel(dds$Type,ref = "control")

    dds <- quiet(DESeq2::DESeq(dds))
    res <- DESeq2::results(dds, contrast = c("Type","treatment","control"))

    if(sort.p) {
      resData <- as.data.frame(res[order(res$padj),])
    } else {
      resData <- as.data.frame(res)
    }
    resData$id <- rownames(resData)
    resData <- resData[,c("id","baseMean","log2FoldChange","lfcSE","stat","pvalue","padj")]
    colnames(resData) <- c("id","baseMean","log2fc","lfcSE","stat","pvalue","padj")
    resData$fc <- 2^resData$log2fc
    if(verbose) {
      resData <- resData[,c("id","fc","log2fc","pvalue","padj")]
    } else {
      resData <- resData[,c("id","fc","log2fc","baseMean","lfcSE","stat","pvalue","padj")]
    }

    write.table(resData, file = outfile, row.names = FALSE, sep = "\t", quote = FALSE)
    cat(paste0("deseq2 of ",compname, " done...\n"))
  }
  options(warn = 0)
}

#' twoclassedger
#' @title Run two class comparison
#' @description Two class differential expression analysis using edgeR algorithm.
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param countsTable A matrix of RNA-Seq raw count data with rows for genes and columns for samples.
#' @param prefix A string value to indicate the prefix of output file.
#' @param overwt A logic value to indicate if to overwrite existing results; FALSE by default.
#' @param sort.p A logic value to indicate if to sort adjusted p value for output table; TRUE by default.
#' @param verbose A logic value to indicate if to only output id, log2fc, pvalue, and padj; TRUE by default.
#' @param res.path A string value to indicate the path for saving the results.
#' @references Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 26(1):139-140.
#'
#' McCarthy DJ, Chen Y, Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40(10):4288-4297.
#' @export
#' @return Several .txt files storing differential expression analysis results by edgeR
#' @keywords internal
#' @importFrom edgeR DGEList calcNormFactors estimateDisp glmFit glmLRT topTags
#' @examples # There is no example and please refer to vignette.
twoclassedger <- function(moic.res    = NULL,
                          countsTable = NULL,
                          prefix      = NULL,
                          overwt      = FALSE,
                          sort.p      = TRUE,
                          verbose     = TRUE,
                          res.path    = getwd()) {

  # Create comparison list for differential expression analysis between two classes.
  createList  <- function(moic.res = NULL) {

    mo.method <- moic.res$mo.method
    moic.res <- moic.res$clust.res
    tumorsam <- moic.res$samID
    sampleList <- list()
    treatsamList <-list()
    treatnameList <- c()
    ctrlnameList <- c()

    n.moic <- length(unique(moic.res$clust))

    for (i in 1:n.moic) {
      sampleList[[i]] <- tumorsam
      treatsamList[[i]] <- intersect(tumorsam, moic.res[which(moic.res$clust == i), "samID"])
      treatnameList[i] <- paste0("CS",i)
      ctrlnameList[i] <- "Others"
    }

    return(list(sampleList, treatsamList, treatnameList, ctrlnameList, mo.method))
  }

  complist <- createList(moic.res = moic.res)

  sampleList <- complist[[1]]
  treatsamList <- complist[[2]]
  treatnameList <- complist[[3]]
  ctrlnameList <- complist[[4]]
  mo.method <- complist[[5]]
  allsamples <- colnames(countsTable)

  options(warn = 1)
  for (k in 1:length(sampleList)) {
    samples <- sampleList[[k]]
    treatsam <- treatsamList[[k]]
    treatname <- treatnameList[k]
    ctrlname <- ctrlnameList[k]

    compname <- paste(treatname, "_vs_", ctrlname, sep="")
    tmp <- rep("others", times = length(allsamples))
    names(tmp) <- allsamples
    tmp[samples] <- "control"
    tmp[treatsam] <- "treatment"

    if(!is.null(prefix)) {
      outfile <- file.path(res.path, paste(mo.method, "_", prefix, "_edger_test_result.", compname, ".txt", sep = ""))
    } else {
      outfile <- file.path(res.path, paste(mo.method,"_edger_test_result.", compname, ".txt", sep = ""))
    }
    if (file.exists(outfile) & (overwt==FALSE)) {
      cat(paste0("edger of ",compname, " exists and skipped...\n"))
      next
    }

    saminfo <- data.frame("Type" = tmp[samples],
                          "SampleID" = samples,
                          stringsAsFactors = FALSE)

    group = factor(saminfo$Type,levels = c("control","treatment"))

    design <- model.matrix(~ group)
    rownames(design) <- samples

    y <- edgeR::DGEList(counts = countsTable[,samples],group = saminfo$Type)
    y <- edgeR::calcNormFactors(y)
    y <- edgeR::estimateDisp(y, design, robust = TRUE)
    fit <- edgeR::glmFit(y, design)
    lrt <- edgeR::glmLRT(fit)
    ordered_tags <- edgeR::topTags(lrt, n = 100000)
    allDiff <- ordered_tags$table
    allDiff <- allDiff[is.na(allDiff$FDR) == FALSE,]
    diff <- allDiff

    diff$id <- rownames(diff)
    resData <- diff[,c("id","logFC","logCPM","LR","PValue","FDR")]
    colnames(resData) <- c("id","log2fc","logCPM","LR","pvalue","padj")
    resData$fc <- 2^resData$log2fc

    if(sort.p) {
      resData <- as.data.frame(resData[order(resData$padj),])
    } else {
      resData <- as.data.frame(resData)
    }
    if(verbose) {
      resData <- resData[,c("id","fc","log2fc","pvalue","padj")]
    } else {
      resData <- resData[,c("id","fc","log2fc","logCPM","LR","pvalue","padj")]
    }
    write.table(resData, file = outfile, row.names = FALSE, sep = "\t", quote = FALSE)
    cat(paste0("edger of ",compname, " done...\n"))
  }
  options(warn = 0)
}

#' twoclasslimma
#' @title Run two class comparison
#' @description Two class differential expression analysis using limma algorithm.
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param norm.expr A matrix of normalized expression data with rows for genes and columns for samples; FPKM or TPM without log2 transformation is recommended.
#' @param prefix A string value to indicate the prefix of output file.
#' @param overwt A logic value to indicate if to overwrite existing results; FALSE by default.
#' @param sort.p A logic value to indicate if to sort adjusted p value for output table; TRUE by default.
#' @param verbose A logic value to indicate if to only output id, log2fc, pvalue, and padj; TRUE by default.
#' @param res.path A string value to indicate the path for saving the results.
#' @references Ritchie, ME, Phipson, B, Wu, D, Hu, Y, Law, CW, Shi, W, and Smyth, GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 43(7):e47.
#' @export
#' @return Several .txt files storing differential expression analysis results by limma
#' @keywords internal
#' @importFrom limma lmFit makeContrasts eBayes topTable
#' @examples # There is no example and please refer to vignette.
twoclasslimma <- function(moic.res  = NULL,
                          norm.expr = NULL,
                          prefix    = NULL,
                          overwt    = FALSE,
                          sort.p    = TRUE,
                          verbose   = TRUE,
                          res.path  = getwd()) {

  # Create comparison list for differential expression analysis between two classes.
  createList  <- function(moic.res = NULL) {

    mo.method <- moic.res$mo.method
    moic.res <- moic.res$clust.res
    tumorsam <- moic.res$samID
    sampleList <- list()
    treatsamList <- list()
    treatnameList <- c()
    ctrlnameList <- c()

    n.moic <- length(unique(moic.res$clust))

    for (i in 1:n.moic) {
      sampleList[[i]] <- tumorsam
      treatsamList[[i]] <- intersect(tumorsam, moic.res[which(moic.res$clust == i), "samID"])
      treatnameList[i] <- paste0("CS",i)
      ctrlnameList[i] <- "Others"
    }

    return(list(sampleList, treatsamList, treatnameList, ctrlnameList, mo.method))
  }

  complist <- createList(moic.res = moic.res)

  sampleList <- complist[[1]]
  treatsamList <- complist[[2]]
  treatnameList <- complist[[3]]
  ctrlnameList <- complist[[4]]
  mo.method <- complist[[5]]
  allsamples <- colnames(norm.expr)

  # log transformation
  if(max(norm.expr) < 25 | (max(norm.expr) >= 25 & min(norm.expr) < 0)) {
    message("--expression profile seems to have been standardised (z-score or log transformation), no more action will be performed.")
    gset <- norm.expr
  }
  if(max(norm.expr) >= 25 & min(norm.expr) >= 0){
    message("--log2 transformation done for expression data.")
    gset <- log2(norm.expr + 1)
  }

  options(warn = 1)
  for (k in 1:length(sampleList)) {
    samples <- sampleList[[k]]
    treatsam <- treatsamList[[k]]
    treatname <- treatnameList[k]
    ctrlname <- ctrlnameList[k]

    compname <- paste(treatname, "_vs_", ctrlname, sep="")
    tmp <- rep("others", times = length(allsamples))
    names(tmp) <- allsamples
    tmp[samples] <- "control"
    tmp[treatsam] <- "treatment"

    if(!is.null(prefix)) {
      outfile <- file.path(res.path, paste(mo.method, "_", prefix, "_limma_test_result.", compname, ".txt", sep = ""))
    } else {
      outfile <- file.path(res.path, paste(mo.method, "_limma_test_result.", compname, ".txt", sep = ""))
    }

    if (file.exists(outfile) & (overwt == FALSE)) {
      cat(paste0("limma of ",compname, " exists and skipped...\n"))
      next
    }

    pd <- data.frame(Samples = names(tmp),
                     Group = as.character(tmp),
                     stringsAsFactors = FALSE)

    design <-model.matrix(~ -1 + factor(pd$Group, levels = c("treatment","control")))
    colnames(design) <- c("treatment","control")

    fit <- limma::lmFit(gset, design = design);
    contrastsMatrix <- limma::makeContrasts(treatment - control, levels = c("treatment", "control"))
    fit2 <- limma::contrasts.fit(fit, contrasts = contrastsMatrix)
    fit2 <- limma::eBayes(fit2, 0.01)
    resData <- limma::topTable(fit2, adjust = "fdr", sort.by = "B", number = 100000)
    resData <- as.data.frame(subset(resData, select=c("logFC","t","B","P.Value","adj.P.Val")))
    resData$id <- rownames(resData)
    colnames(resData) <- c("log2fc","t","B","pvalue","padj","id")
    resData$fc <- 2^resData$log2fc

    if(sort.p) {
      resData <- resData[order(resData$padj),]
    } else {
      resData <- as.data.frame(resData)
    }
    if(verbose) {
      resData <- resData[,c("id","fc","log2fc","pvalue","padj")]
    } else {
      resData <- resData[,c("id","fc","log2fc","t","B","pvalue","padj")]
    }
    write.table(resData, file = outfile, row.names = FALSE, sep = "\t", quote = FALSE)
    cat(paste0("limma of ",compname, " done...\n"))
  }
  options(warn = 0)
}

#' runDEA
#' @title Run differential expression analysis
#' @description Using choosen algorithm to run differential expression analysis between two classes identified by multi-omics clustering process.
#' @param dea.method A string value to indicate the algorithm for differential expression analysis. Allowed value contains c('deseq2', 'edger', 'limma'). The former two require RNA-Seq raw count data and the last one requires normalized expression data (FPKM or TPM without log2 transformation is recommended).
#' @param expr A matrix of expression data.
#' @param moic.res An object returned by `getMOIC()` with one specified algorithm or `get\%algorithm_name\%` or `getConsensusMOIC()` with a list of multiple algorithms.
#' @param prefix A string value to indicate the prefix of output file.
#' @param overwt A logic value to indicate if to overwrite existing results; FALSE by default.
#' @param sort.p A logic value to indicate if to sort adjusted p value for output table; TRUE by default.
#' @param verbose A logic value to indicate if to only output id, log2fc, pvalue, and padj; TRUE by default.
#' @param res.path A string value to indicate the path for saving the results.
#' @export
#' @return Several .txt files storing differential expression analysis results by specified algorithm
#' @importFrom limma lmFit makeContrasts eBayes topTable
#' @importFrom DESeq2 DESeqDataSetFromMatrix DESeq results
#' @importFrom edgeR DGEList calcNormFactors estimateDisp glmFit glmLRT topTags
#' @references Love, M.I., Huber, W., Anders, S. (2014) Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol, 15(12):550-558.
#'
#' Robinson MD, McCarthy DJ and Smyth GK (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26(1):139-140.
#'
#' McCarthy DJ, Chen Y, Smyth GK (2012). Differential expression analysis of multifactor RNA-Seq experiments with respect to biological variation. Nucleic Acids Res. 40(10):4288-4297.
#'
#' Ritchie, ME, Phipson, B, Wu, D, Hu, Y, Law, CW, Shi, W, and Smyth, GK (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res, 43(7):e47.
#'
#' @examples # There is no example and please refer to vignette.
runDEA <- function(dea.method = c("deseq2", "edger", "limma"),
                   expr       = NULL,
                   moic.res   = NULL,
                   prefix     = NULL,
                   overwt     = FALSE,
                   sort.p     = TRUE,
                   verbose    = TRUE,
                   res.path   = getwd()) {

  if(!is.element(dea.method, c("deseq2", "edger", "limma"))) {
    stop("unsupported algorithm: dea.method should be one of 'deseq2', 'edger', or 'limma'.")
  }

  method <- dea.method[1] # deseq2 by default

  comsam <- intersect(moic.res$clust.res$samID, colnames(expr))
  # check data
  if(length(comsam) == nrow(moic.res$clust.res)) {
    message("--all samples matched.")
  } else {
    message(paste0("--",(nrow(moic.res$clust.res)-length(comsam))," samples mismatched from current subtypes."))
  }

  moic.res$clust.res <- moic.res$clust.res[comsam,,drop = FALSE]
  expr <- expr[,comsam]

  rundea <- switch(method,
                   "deseq2" = twoclassdeseq2,
                   "edger"  = twoclassedger,
                   "limma"  = twoclasslimma)
  if(method %in% c("deseq2", "edger")) {
    message(paste0("--you choose ", method," and please make sure an RNA-Seq count data was provided."))
    rundea(moic.res    = moic.res,
           countsTable = expr,
           prefix      = prefix,
           overwt      = overwt,
           sort.p      = sort.p,
           verbose     = verbose,
           res.path    = res.path)
  } else {
    message(paste0("--you choose ", method," and please make sure a microarray profile or a normalized expression data [FPKM or TPM without log2 transformation is recommended] was provided."))
    rundea(moic.res    = moic.res,
           norm.expr   = expr,
           prefix      = prefix,
           overwt      = overwt,
           sort.p      = sort.p,
           verbose     = verbose,
           res.path    = res.path)
  }
}