R/PairedPSCBS.backgroundCorrect.R

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

setMethodS3("translateC1C2", "PairedPSCBS", function(fit, dC1=0, dC2=0, sC1=1, sC2=1, ..., verbose=FALSE) {
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Validate arguments
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  dC1 <- Arguments$getNumeric(dC1)
  dC2 <- Arguments$getNumeric(dC2)
  sC1 <- Arguments$getNumeric(sC1)
  sC2 <- Arguments$getNumeric(sC2)

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


  verbose && enter(verbose, "Scaling and translating by (C1,C2)")
  verbose && cat(verbose, "sC1: ", sC1)
  verbose && cat(verbose, "sC2: ", sC2)
  verbose && cat(verbose, "dC1: ", dC1)
  verbose && cat(verbose, "dC2: ", dC2)

  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Extract data
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  segs <- getSegments(fit, splitters=FALSE)
  .stop_if_not(!is.null(segs))

  nbrOfSegments <- nrow(segs)
  verbose && cat(verbose, "Number of segments: ", nbrOfSegments)

  # (C1,C2,...)
  X <- extractC1C2(fit, splitters=FALSE)

  # (C1,C2)
  C1C2 <- X[,1:2, drop=FALSE]

  C1C2[,1] <- dC1 + sC1*C1C2[,1]
  C1C2[,2] <- dC2 + sC2*C1C2[,2]

  verbose && enter(verbose, "(C1,C2) to (TCN,DH)")
  # (C1,C2) -> (TCN,DH)
  gamma <- rowSums(C1C2, na.rm=TRUE)
  dh <- 2*(C1C2[,2]/gamma - 1/2)
  verbose && exit(verbose)

  # Update segmentation means
  segs[,"tcnMean"] <- gamma
  segs[,"dhMean"] <- dh
  segs[,"c1Mean"] <- C1C2[,1]
  segs[,"c2Mean"] <- C1C2[,2]

  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Return results
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  fitO <- fit
  fitO$output <- segs

  verbose && exit(verbose)

  fitO
})



setMethodS3("transformC1C2", "PairedPSCBS", function(fit, fcn, ..., verbose=FALSE) {
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Validate arguments
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Argument 'fcn':
  .stop_if_not(is.function(fcn))

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


  verbose && enter(verbose, "Transform (C1,C2) by a function")

  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Extract data
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  segs <- getSegments(fit, splitters=FALSE)
  .stop_if_not(!is.null(segs))

  nbrOfSegments <- nrow(segs)
  verbose && cat(verbose, "Number of segments: ", nbrOfSegments)

  # (C1,C2,...)
  X <- extractC1C2(fit, splitters=FALSE)

  # (C1,C2)
  C1C2 <- X[,1:2, drop=FALSE]

  C1C2 <- fcn(C1C2, ...)

  verbose && enter(verbose, "(C1,C2) to (TCN,DH)")
  # (C1,C2) -> (TCN,DH)
  gamma <- rowSums(C1C2, na.rm=TRUE)
  dh <- 2*(C1C2[,2]/gamma - 1/2)
  verbose && exit(verbose)

  # Update segmentation means
  segs[,"tcnMean"] <- gamma
  segs[,"dhMean"] <- dh
  segs[,"c1Mean"] <- C1C2[,1]
  segs[,"c2Mean"] <- C1C2[,2]

  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Return results
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  fitO <- fit
  fitO$output <- segs

  verbose && exit(verbose)

  fitO
})



setMethodS3("estimateC2Bias", "PairedPSCBS", function(fit, ...) {
  # Identify region in allelic balance
  segs <- getSegments(fit, splitters=FALSE)
  abCall <- segs$abCall
  if (is.null(abCall)) {
    throw("Allelic balance has not been called.")
  }
  idxs <- which(segs$abCall)
  segs <- segs[idxs,,drop=FALSE]

  # Extract (TCN,DH)
  gamma <- segs[, "tcnMean"]
  rho <- segs[, "dhMean"]

  # Calculate (C1,C2)
  C1 <- 1/2 * (1 - rho) * gamma
  C2 <- gamma - C1

  # Calculate bias in C2
  dC2 <- C2 - C1

  # Calculate weighted average of all C2 biases
  w <- segs[,"dhNbrOfLoci"]
  w <- w / sum(w, na.rm=TRUE)
  dC2 <- weightedMedian(dC2, w=w)

  dC2
}) # estimateC2Bias()



setMethodS3("backgroundCorrect", "PairedPSCBS", function(fit, targetMedian=c("same", "none"), ...) {
  # Argument 'targetMedian':
  if (is.character(targetMedian)) {
    targetMedian <- match.arg(targetMedian)
  } else {
    targetMedian <- Arguments$getDouble(targetMedian, range=c(0,Inf))
  }

  modelFit <- list()

  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # (1) Estimate background (e.g. normal contamination and more)
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  kappa <- estimateKappa(fit)
  modelFit$kappa <- kappa
  fitBG <- translateC1C2(fit, dC1=-kappa, dC2=-kappa)


  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # (2) Rescale
  # - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
  # Calculate current median
  modelFit$targetMedian <- targetMedian
  if (!identical(targetMedian, "none")) {
    segsBG <- getSegments(fitBG)
    count <- segsBG[["tcnNbrOfLoci"]]
    # chrs <- 1:22
    # if (!is.null(chrs)) segsBG <- subset(segsBG, chromosome %in% chrs)
    CTBG <- segsBG[["tcnMean"]]
    muBG <- weightedMedian(CTBG, w=count)
    muBG <- Arguments$getDouble(muBG, range=c(0,Inf))
    modelFit$muBG <- muBG

    if (identical(targetMedian, "same")) {
      segs <- getSegments(fit)
      # if (!is.null(chrs)) segs <- subset(segs, chromosome %in% chrs)
      CT <- segs[["tcnMean"]]
      mu <- weightedMedian(CT, w=count)
      mu <- Arguments$getDouble(mu, range=c(0,Inf))
    } else if (is.numeric(targetMedian)) {
      mu <- targetMedian
    }
    modelFit$mu <- mu

    scale <- mu/muBG
    modelFit$scale <- scale

    fitBG <- translateC1C2(fitBG, sC1=scale, sC2=scale)
  } # if (!identical(targetMedian, "none"))


  # Store model fit
  postMethods <- fit$postMethods
  if (is.null(postMethods)) postMethods <- list()
  postMethods$backgroundCorrection <- list(modelFit=modelFit)
  fitBG$postMethods <- postMethods

  fitBG
}) # backgroundCorrect()



##############################################################################
# HISTORY
# 2013-08-21
# o Retired this version of deShearC1C2() by renaming it.
# 2012-09-22
# o Added backgroundCorrect() for PairedPSCBS.
# o Now fitDeltaC1C2ShearModel() "handles" also single elements.
# o BUG FIX: deShearC1C2() would incorrectly set the TCN level to zero
#   if (C1,C2) = (NA,NA).
# o Now deShearC1C2() stores the modelFit in postMethods$deShearC1C2$modelFit.
# o CORRECTION: The peak caller of fitDeltaXYShearModel() swapped the labels
#   of the horizontal and vertical peaks. Labels are only used for display.
# 2012-09-21 [HB]
# o BUG FIX: fitDeltaC1C2ShearModel(... dropABChangePoints=TRUE) could
#   generate an incorrect number change-point weights if dropping.
# o Added argument 'onError' to fitDeltaXYShearModel() for matrix.
# o Now attribute 'modelFit' fitDeltaXYShearModel() also contains
#   debug$cpAngleDensity.
# 2012-09-21 [PN]
# o Now fitDeltaC1C2ShearModel() will apply a fast naive AB caller,
#   iff AB calls are not available.
# o BUG FIX: horizontal/vertical axes were swapped for direction '|_'
# 2012-09-19 [HB]
# o MEMORY/"BUG FIX": fitDeltaXYShearModel() would return "H" deshearing
#   functions that each would hold the very large fitDeltaXYShearModel()
#   call environment.  Now we make sure each of those objects only
#   contains the few scalar parameter estimates needed.
# 2012-09-19 [PN]
# o Added direction '|_': similar to '|-' but enforces equal amount of
#   de-shearing in both dimensions and does not change TCN.
# 2011-10-17 [HB]
# o Added argument 'dropABChangePoints' to fitDeltaC1C2ShearModel().
# 2011-10-16 [HB]
# o Now deShearC1C2(), translateC1C2() and transformC1C2() also update
#   C1 and C2 mean levels.
# 2011-07-10 [HB]
# o Updated code to work with the new column names in PSCBS v0.11.0.
# 2010-10-20 [HB]
# o Now fitDeltaXYShearModel() uses both -pi/4 and +pi/4 to estimate
#   the diagonal parameters.
# 2010-10-08 [HB]
# o Added estimateC2Bias().
# o Added fitDeltaXYShearModel().
# o Now deShearC1C2() uses fitDeltaC1C2ShearModel().
# o Added fitDeltaC1C2ShearModel().
# o Now deShearC1C2() returns the 'modelFit'.
# o Now deShearC1C2() calls peaks using callPeaks() for PeaksAndValleys.
# 2010-09-26 [HB]
# o Added argument 'adjust' to deShearC1C2() with new default.
# o Added sanity checks to deShearC1C2().
# o Now normalizeBAFsByRegions() for PairedPSCBS handles multiple chromosomes.
# 2010-09-22 [PN]
# o Added deShearC1C2() for PairedPSCBS.
# 2010-09-19 [HB+PN]
# o Added orthogonalizeC1C2() for PairedPSCBS.
# 2010-09-15 [HB]
# o Added Rdocs for callCopyNeutralRegions().
# 2010-09-09 [HB]
# o Added callCopyNeutralRegions() for PairedPSCBS.
# 2010-09-08 [HB]
# o Added subsetBySegments() for PairedPSCBS.
# o Added Rdocs with an example.
# o Added normalizeBAFsByRegions() for PairedPCSBS.
# o Created.
##############################################################################