R/SBL.R

In isglobal-brge/R-GADA: Genome Alteration Detection Algorithm (GADA)

`SBL` <-
  function (x, sigma2, aAlpha=0.2, estim.sigma2 = FALSE, maxit = 10000, 
            tol = 1e-08, verbose = FALSE, saveInfo = TRUE) 
{
  if (!inherits(x, "setupGADA")) 
    stop("object must be of class 'setupGADA' ")
  if (missing(sigma2) & !estim.sigma2) 
    stop("parameter 'sigma2' is missing. \n It can be estimated using estim.sigma2=TRUE")
  if (estim.sigma2) 
    sigma2 <- -1
  gen.info <- attr(x, "gen.info")
  
  debug <- ifelse(verbose, 1, 0)

  x$log.ratio[is.na(x$log.ratio)] <- 0  # RPR we finaly use this to deal with the NA?

  if (is.null(gen.info)) {
    out <- list()
    y <- x$log.ratio
    M <- length(y)
    if (sigma2 < 0) {
      sigma2 <- 0.5 * mean(diff(y)^2, trim = 0.01, na.rm = TRUE)
      cat("    The estimated sigma2 =", sigma2, "\n")
    }
    out <- SBL.fit(0, y, NULL, NULL, NULL, sigma2, aAlpha, 
                   maxit, tol, debug)
  }
  else {
    # chr <- unique(gen.info$chr) # Modified JRG 08/2018!!! important!
    # chr <- chr[!is.na(chr)]
    chr <- levels(gen.info$chr)
    n.chr <- length(chr)
    control <- 1:nrow(gen.info)
    if (sigma2 < 0) {
      sigma2 <- 0
      totalM <- 0
 
      for (i in 1:n.chr) {

# Changed to control missing data, JRG
#        selec <- control[gen.info$chr == chr[i]]

        tt <- gen.info$chr==chr[i] & !is.na(gen.info$chr)
        selec <- control[tt]

        y <- x$log.ratio[selec]
        M <- length(y)
        totalM <- totalM + (M - 1)
        sigma2 <- sigma2 + (M - 1) * 0.5 * mean(diff(y)^2, 
                                                trim = 0.01, na.rm=TRUE)
      }
      sigma2 <- sigma2/totalM
      cat("    The estimated sigma2 =", sigma2, "\n")
    }
    
    out <- lapply(1:n.chr, SBL.fit, x = x$log.ratio, 
                  control = control, gen.info = gen.info, chr = chr, 
                  sigma2 = sigma2, aAlpha = aAlpha, maxit = maxit, 
                  tol = tol, debug = debug)
    attr(out, "chr") <- chr
  }

  # RPR Potential BUG!!! Do we also account for the removed NA probes in the plots?
  
  attr(out, "data") <- x$log.ratio
  if (saveInfo) {
    attr(out, "gen.info") <- attr(x, "gen.info")
  }
  class(out) <- "SBL"
  out
}