R/pruneCNA.R

In aroma.cn: Copy-Number Analysis of Large Microarray Data Sets

# \section{Required method implementations}{
#  In order for this method to work, the following methods need to be
#  implemented for the class of argument \code{fit}:
#  \itemize{
#   \item \code{findAtomicAberrations()}
#   \item \code{mergeTwoSegments()}
#  }
# }
setMethodS3("pruneCNA", "AbstractCBS", function(fit, ..., maxGeneration=Inf, onAtomicIsland=NULL, verbose=FALSE) {
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Validate arguments
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Argument 'maxGeneration':
  maxGeneration <- Arguments$getDouble(maxGeneration, range=c(0,Inf))

  # Argument 'verbose':
  verbose <- Arguments$getVerbose(verbose)
  if (verbose) {
    pushState(verbose)
    on.exit(popState(verbose))
  }


  fitT <- fit

  fitList <- list()
  gg <- 1
  while (gg <= maxGeneration) {
    verbose && enter(verbose, sprintf("Generation #%d", gg))

    fitList[[gg]] <- fitT

    nbrOfSegments <- nbrOfSegments(fitT)
    maxH <- nbrOfSegments-2L

    # Done
    if (maxH < 0) {
      break
    }

    hasChanged <- FALSE
    for (hh in 0:maxH) {
#    for (hh in 0) {
      verbose && enter(verbose, sprintf("Block size H=%d of %d", hh, maxH))

      res <- findAtomicAberrations(fitT, H=hh, ..., verbose=verbose)
      verbose && str(verbose, res)

      # (i) Atomic islands?
      if (hh == 0) {
        atomicIslands <- res$atomicRegions
      } else {
        atomicIslands <- res$atomicIslands
      }

      if (length(atomicIslands) > 0) {
        verbose && printf(verbose, "Atomic islands found (H=%d):\n", hh)
        verbose && print(verbose, atomicIslands)

        # Overlapping atomic islands?
        if (length(atomicIslands) > 1) {
          regionsHH <- matrix(c(atomicIslands, atomicIslands+hh), ncol=2L, byrow=FALSE)
          colnames(regionsHH) <- c("from", "to")
          rownames(regionsHH) <- sprintf("Atomic island #%d", seq_along(atomicIslands))
          verbose && print(verbose, regionsHH)

          froms <- regionsHH[-1,"from"]
          tos <- regionsHH[-length(atomicIslands),"to"]
          isOverlapping <- (froms <= tos)
          verbose && printf(verbose, "Overlapping: %s\n", any(isOverlapping))

          if (any(isOverlapping)) {
            verbose && cat(verbose, "Overlapping atomic islands. Dropping only the first.")
            atomicIslands <- atomicIslands[1]
          }
        }

        # Drop atomic islands and merge flanking segments
        dropList <- list()
        atomicIslands <- sort(atomicIslands, decreasing=TRUE)
        for (kk in seq_along(atomicIslands)) {
          atomicIsland <- atomicIslands[kk]
          if (hh == 0) {
            atomicIslandTag <- sprintf("change point #%d", atomicIsland)
          } else if (hh == 1) {
            atomicIslandTag <- sprintf("segment #%d", atomicIsland)
          } else {
            atomicIslandTag <- sprintf("segments #%d-#%d", atomicIsland, atomicIsland+(hh-1L))
          }
          verbose && enter(verbose, sprintf("Atomic island #%d ('%s') of %d", kk, atomicIslandTag, length(atomicIslands)))
          n0 <- nbrOfSegments(fitT)
          verbose && cat(verbose, "Number of segments before: ", n0)
          fitTT <- fitT

          if (hh > 0) {
            verbose && enter(verbose, sprintf("Dropping %d segments (%s)", hh, atomicIslandTag))
            fitTT <- dropRegions(fitTT, regions=atomicIsland, H=hh)
            nT <- nbrOfSegments(fitTT)
            verbose && exit(verbose)
            verbose && cat(verbose, "Number of segments left: ", nT)
            # Sanity check
            .stop_if_not(n0-nT == hh)
          }

          fitDrop <- fitTT$dropped
          fitTT$dropped <- NULL
          dropList[[kk]] <- fitDrop

          verbose && enter(verbose, sprintf("Merging segments (#%d and #%d)", atomicIsland-1L, atomicIsland))
          fitTT <- mergeTwoSegments(fitTT, left=atomicIsland-1L)

          if (is.function(onAtomicIsland)) {
            onAtomicIsland(fit0=fitT, fit1=fitTT, fitD=fitDrop,
                           atomicIsland=atomicIsland, H=hh)
          }
          n1 <- nbrOfSegments(fitTT)
          verbose && exit(verbose)
          verbose && cat(verbose, "Number of segments left: ", n1)
          # Sanity check
          .stop_if_not(n0-n1 == hh+1)

          fitT <- fitTT

          verbose && exit(verbose)
        } # for (kk ...)

        fitT$dropped <- dropList
        fitT$atomicIslands <- rev(atomicIslands)
        fitT$H <- hh

        hasChanged <- TRUE

        verbose && exit(verbose)

        # Go to next generation
        break
      } # if (length(atomicIslands) > 0)

      verbose && exit(verbose)
    } # for (hh ...)

    if (!hasChanged) {
      break
    }

    # Next generation
    gg <- gg + 1L

    verbose && exit(verbose)
  } # while(gg <= maxGeneration)

  class(fitList) <- c("PruneCNA", class(fitList))

  fitList
}) # prunceCNA()