R/topCounts.R

In baySeq: Empirical Bayesian analysis of patterns of differential expression in count data

# modification on git from copied files
selectTop <- function(cD, group, ordering, orderings = TRUE, decreasing = TRUE, number = 10, likelihood, FDR, FWER, posteriors) {
  if(missing(likelihood)) likelihood <- NULL
  if(missing(FDR)) FDR <- NULL
  if(missing(FWER)) FWER <- NULL
  if(missing(ordering)) ordering <- NULL

  if(missing(posteriors)) posteriors <- NULL
  
  if(!missing(group)) {
      st <- .selectTags(cD, group = group, ordering = ordering, decreasing = decreasing, number = number, likelihood = likelihood, FDR = FDR, FWER = FWER, posteriors = posteriors)
      if(!is.null(st)) return(cD[st,]) else return(NULL)
  }

  if(length(cD@nullPosts) > 0) {
    st <- .selectTags(cD, NULL, ordering = NULL, decreasing = decreasing, number = number, likelihood, FDR, FWER, posteriors = posteriors)
    if(!is.null(st)) nullCD <- cD[st,] else nullCD <- NULL
  } else nullCD <- NULL
  
  
  
  if(!orderings) {
    selOrd <- lapply(1:length(cD@groups), function(ii) cD[.selectTags(cD, ii, ordering = NULL, decreasing = decreasing, number = number, likelihood, FDR, FWER, posteriors = posteriors),])                     
    names(selOrd) <- names(cD@groups)
    } else {
      selOrd <- do.call("c", lapply(1:ncol(cD@orderings),function(ii) {
        selord <- lapply(levels(cD@orderings[,ii]), function(ord) {
          st <- .selectTags(cD, ii, ordering = ord, decreasing = decreasing, number = number, likelihood, FDR, FWER, posteriors = posteriors)
          if(!is.null(st)) cD[st,] else NULL
        })
        names(selord) <- paste(colnames(cD@orderings)[ii], ":", levels(cD@orderings[,ii]), sep = "")
        names(selord) <- gsub(":$", "", names(selord))
        selord
      }))
    }
  if(length(cD@nullPosts) > 0) return(c(list(null = nullCD), selOrd)) else return(selOrd)
}
  
#.selectTop <- function(cD, group, ordering, decreasing = TRUE, number = 10, likelihood, FDR, FWER) #{
#  if(missing(likelihood)) likelihood <- NULL
#  if(missing(FDR)) FDR <- NULL
#  if(missing(FWER)) FWER <- NULL
#  if(missing(ordering)) ordering <- NULL  
  
#  selTags <- .selectTags(cD, group, ordering, decreasing = decreasing, number = number, likelihood, FDR, FWER)
#  cD[selTags,]
#}

.selectTags <- function(cD, group, ordering, decreasing = TRUE, number = 10, likelihood, FDR, FWER, posteriors) {
    if(!inherits(cD, what = "countData"))
        stop("variable 'cD' must be of or descend from class 'countData'")
    if(nrow(cD@posteriors) == 0)
        stop("The '@posteriors' slot of cD is empty!")
  
    if(is.character(group))
        group <- pmatch(group, names(cD@groups))
    if(!is.null(group) && is.na(group)) stop("Can't match this group name.")
  
    if(!is.null(ordering)) {
        ordCD <- which(cD@orderings[,group] == ordering)
        cD <- cD[ordCD,]
        if(!is.null(posteriors)) posteriors <- posteriors[ordCD]
    }
    
    if(is.null(posteriors)) {
        if(is.null(group)) {
            if(length(cD@nullPosts) == 0)
                stop("The '@nullPosts' slot of cD is empty - you can't use 'group = NULL'.")
            likes <- cD@nullPosts
            neglikes <- cD@posteriors
        } else {
            likes <- cD@posteriors[,group,drop = FALSE]    
            if(length(cD@nullPosts) > 0) neglikes <- cbind(cD@posteriors[,-group,drop=FALSE], cD@nullPosts) else neglikes <- cD@posteriors[,-group, drop = FALSE]
        }
    } else {
        likes <- matrix(posteriors, ncol = 1)
        neglikes <- matrix(log(1 - exp(likes)), ncol = 1)
    }
  
    ordgroups <- order(.logRowSum(neglikes), decreasing = !decreasing)

  cutNumber <- c()
  if(!is.null(likelihood))
      cutNumber <- c(cutNumber, sum(likes > log(likelihood), na.rm = TRUE))
  if (!is.null(FDR))
      cutNumber <- c(cutNumber, sum(cumsum(1 - exp(likes[ordgroups[1:sum(!is.na(likes[,1]))],1])) / 1:sum(!is.na(likes[,1])) < FDR, na.rm = TRUE))
  if (!is.null(FWER))
      cutNumber <- c(cutNumber, sum(1 - cumprod(exp(likes[ordgroups,1])) < FWER, na.rm = TRUE))
  
  if(!is.null(likelihood) | !is.null(FDR) | !is.null(FWER)) {
      number <- min(cutNumber)
      if(number == 0) warning("No features were found using the cutoffs for likelihood, FDR and FWER specified")
   }
  number <- min(number, nrow(likes))

  if(number > 0) {
    selTags <- ordgroups[1:number]
    if(!is.null(ordering)) selTags <- ordCD[selTags]
    return(selTags)
  } else return(NULL)
}
  
`topCounts` <-
function(cD, group, ordering, decreasing = TRUE, number = 10, likelihood, FDR, FWER, normaliseData = FALSE, posteriors)
  {
      if(missing(likelihood)) likelihood <- NULL
      if(missing(FDR)) FDR <- NULL
      if(missing(FWER)) FWER <- NULL
      if(missing(ordering)) ordering <- NULL
      
      if(is.character(group))
          group <- pmatch(group, names(cD@groups))
      if(!is.null(group) && is.na(group)) stop("Can't match this group name.")
      
      if(missing(posteriors)) posteriors <- NULL
      
      selTags <- .selectTags(cD, group, ordering, decreasing = decreasing, number = number, likelihood, FDR, FWER, posteriors)
      if(all(is.null(selTags))) return(NULL)
      
      if(is.null(posteriors)) {
          if(!is.null(group))
              likes <- cD@posteriors[selTags,group]
          else likes <- cD@nullPosts[selTags,1]
      } else {
          likes <- posteriors[selTags,]
      }

      topCD <- cD[selTags,]
      topCD@posteriors <- matrix(likes, nrow = nrow(topCD))
      if(length(cD@orderings) > 0 && !is.null(group)) topCD@orderings <- topCD@orderings[, group, drop = FALSE] else topCD@orderings <- data.frame()
      if(nrow(topCD@annotation) == 0) topCD@annotation <- data.frame(rowID = selTags)
      
      topTags <- flatten(topCD)
      topTags <- cbind(topTags,
                       FDR = cumsum(1 - exp(likes)) / 1:length(selTags),
                       FWER = 1 - cumprod(exp(likes)))
      names(topTags)[names(topTags) == "FDR"] <- paste("FDR", names(cD@groups)[group[1]], sep = ".")
      names(topTags)[names(topTags) == "FWER"] <- paste("FWER", names(cD@groups)[group[1]], sep = ".")      
      topTags
  }