R/itemselect.R

In aursiber/DRomics: Dose Response for Omics

### select significantly responsive items 
itemselect <- function(omicdata, select.method = c("quadratic", "linear", "ANOVA"), FDR = 0.05, 
                       max.ties.prop = 0.2)
{
  # Checks
  if (!(inherits(omicdata, "microarraydata") | 
        inherits(omicdata, "RNAseqdata") |
        inherits(omicdata, "continuousomicdata") |
        inherits(omicdata, "continuousanchoringdata")))
    stop("Use only with 'microarraydata', 'RNAseqdata', 'continuousomicdata' or 
    'continuousanchoringdata' objects, respectively
    created with functions 'microarraydata()', 'RNAseqdata()', 'metabolomicdata()'
    or 'continuousomicdata()'
    and 'continuousanchoringdata()'.")
  select.method <- match.arg(select.method, c("quadratic", "linear", "ANOVA"))
  if (select.method == "ANOVA")
  {
    # Check that there are replicates for at least half of the doses
    nbdoses <- length(omicdata$design)
    nbdoseswithoutrep <- sum(omicdata$design == 1)
    if (nbdoseswithoutrep > nbdoses / 2) 
      stop("The ANOVA method should not be used when there are more than half of the doses 
           without replicates. Choose the quadratic or linear trend test for such a design.")
  }
  if (!is.numeric(FDR))
    stop("FDR, the false discovery rate, must a number in ]0; 1[ (generally under 0.1).")
  if ((FDR <=0) | (FDR >=1))
    stop("FDR, the false discovery rate, must in ]0; 1[.")
  if ((max.ties.prop <=0) | (max.ties.prop >0.5))
    stop("max.ties.prop, the maximal tolerated proportion of tied values per item, must in ]0; 0.5].")
  
  item <- omicdata$item
  dose <- omicdata$dose
  fdose <- as.factor(dose)
  doseranks <- as.numeric(as.factor(dose))
  doseranks2 <- doseranks*doseranks
  irow <- seq_len(length(item))
  
  if (inherits(omicdata,"microarraydata") | 
      inherits(omicdata,"continuousomicdata"))
  {
    data <- omicdata$data
    if (select.method == "quadratic")
    {
      design4lmFit <- stats::model.matrix(~ doseranks + doseranks2)
    } else
      if (select.method == "linear")
      {
        design4lmFit <- stats::model.matrix(~ doseranks)
      } else
        if (select.method == "ANOVA")
        {
          design4lmFit <- stats::model.matrix(~ fdose)
        } 
    
    # Selection using limma    
    fit <- lmFit(data, design4lmFit)
    fitBayes <- eBayes(fit)
    # all adjusted pvalues without sorting
    res <- topTable(fitBayes, adjust.method = "BH", number = nrow(data), sort.by = "none") 
    
    selectindexnonsorted <- irow[res$adj.P.Val < FDR]
    adjpvaluenonsorted <- res$adj.P.Val[selectindexnonsorted]
    irowsortedbypvalue <- order(adjpvaluenonsorted, decreasing = FALSE)
    selectindex <- selectindexnonsorted[irowsortedbypvalue]
    adjpvalue <- adjpvaluenonsorted[irowsortedbypvalue]
  } else
    
  if (inherits(omicdata,"continuousanchoringdata"))
  {
    data <- omicdata$data
    nitem <- nrow(data)
    pvalue <- numeric(length = nitem)
    # Selection using lm    
    if (select.method == "quadratic")
    {
      for (i in 1:nitem) # to write in a sapply in case there are a lot of endpoints
      {
        lmFit <- stats::lm(data[i, ] ~ doseranks + doseranks2)
        lmFitconst <- stats::lm(data[i, ] ~ 1)
        a <- stats::anova(lmFit, lmFitconst)
        pvalue[i] <- a[["Pr(>F)"]][2]
      }
    } else
    if (select.method == "linear")
    {
      for (i in 1:nitem) # to write in a sapply in case there are a lot of endpoints
      {
        lmFit <- stats::lm(data[i, ] ~ doseranks)
        lmFitconst <- stats::lm(data[i, ] ~ 1)
        a <- stats::anova(lmFit, lmFitconst)
        pvalue[i] <- a[["Pr(>F)"]][2]
      }
    } else
    if (select.method == "ANOVA")
    {
      for (i in 1:nitem) # to write in a sapply in case there are a lot of endpoints
      {
        lmFit <- stats::lm(data[i, ] ~ fdose)
        lmFitconst <- stats::lm(data[i, ] ~ 1)
        a <- stats::anova(lmFit, lmFitconst)
        pvalue[i] <- a[["Pr(>F)"]][2]
      }
    } 
      
    # all adjusted pvalues without sorting after Benjamini Hochberg procedure
    wholeadjpvalue <- stats::p.adjust(pvalue, method = "BH") 

    selectindexnonsorted <- irow[wholeadjpvalue < FDR]
    adjpvaluenonsorted <- wholeadjpvalue[selectindexnonsorted]
    irowsortedbypvalue <- order(adjpvaluenonsorted, decreasing = FALSE)
    selectindex <- selectindexnonsorted[irowsortedbypvalue]
    adjpvalue <- adjpvaluenonsorted[irowsortedbypvalue]
    
  } else
      
    
  if (inherits(omicdata,"RNAseqdata"))
  {
    data <- omicdata$raw.counts
    coldata <- data.frame(fdose = fdose, doseranks = doseranks, doseranks2 = doseranks2)
    rownames(coldata) <- colnames(data)
    
    if (select.method == "quadratic")
    {
      dds <- DESeqDataSetFromMatrix(
        countData = data,
        colData = coldata,
        design = ~ doseranks + doseranks2
      )
    } else
      if (select.method == "linear")
      {
        dds <- DESeqDataSetFromMatrix(
          countData = data,
          colData = coldata,
          design = ~ doseranks
        )
      } else
        if (select.method == "ANOVA")
        {
          dds <- DESeqDataSetFromMatrix(
            countData = data,
            colData = coldata,
            design = ~ fdose
          )
        } 
    
    # Selection using DESeq2    
    dds <- DESeq(dds, test = "LRT", reduced = ~ 1)
    res <- results(dds)
    # all adjusted pvalues without sorting in res$padj

    selectindexnonsorted <- irow[(res$padj < FDR) & !is.na(res$padj)]
    adjpvaluenonsorted <- res$padj[selectindexnonsorted]
    irowsortedbypvalue <- order(adjpvaluenonsorted, 
                                decreasing = FALSE, na.last = TRUE)
    selectindex <- selectindexnonsorted[irowsortedbypvalue]
    adjpvalue <- adjpvaluenonsorted[irowsortedbypvalue]
  }
  
  # elimination of first selected items with a proportion of tied values
  # above max.tied.values.prop (assuming tied values correspond to the min 
  # value, at which non detection are imputed
  nsample <- length(dose)
  max4nties <- nsample * max.ties.prop
  if (length(selectindex) != 0) 
  {
    check.ties <- function(index)
    {
      datai <- data[index, ]
      mini <- min(datai)
      nbtiesi <- length(which(datai == mini))
      return(nbtiesi < max4nties)
    }
    tokeep <- sapply(selectindex, check.ties)
    selectindex <- selectindex[tokeep]
    adjpvalue <- adjpvalue[tokeep]
  } else
  {
    warning(strwrap(prefix = "\n", initial = "\n", 
                    "NO ITEM WAS SELECTED."))
  }
  
  reslist <- list(adjpvalue = adjpvalue, selectindex = selectindex, 
                  omicdata = omicdata, select.method = select.method, FDR = FDR)  
  
  return(structure(reslist, class = "itemselect"))
}

print.itemselect <- function(x, nfirstitems = 20, ...)
{
  if (!inherits(x, "itemselect"))
    stop("Use only with 'itemselect' objects.")
  
  if (x$select.method == "ANOVA")
  {
    cat("Number of selected items using an ANOVA type test with an FDR of ", x$FDR, ": ", length(x$selectindex), "\n", sep = "")
  } else
    if (x$select.method == "linear")
    {
      cat("Number of selected items using a linear trend test with an FDR of ", x$FDR, ": ", length(x$selectindex), "\n", sep = "")
    } else
      if (x$select.method == "quadratic")
      {
        cat("Number of selected items using a quadratic trend test with an FDR of ", x$FDR, ": ", length(x$selectindex), "\n", sep = "")
      } 
  
  if (length(x$selectindex) > nfirstitems) 
  {
    cat("Identifiers of the first ", nfirstitems, " most responsive items:\n", sep = "")
    print(x$omicdata$item[x$selectindex[1:nfirstitems]])
  } else
  {
    cat("Identifiers of the responsive items:\n")
    print(x$omicdata$item[x$selectindex])
  }
}