segmentByMPCBS.RawGenomicSignals: Segment copy numbers using the multi-platform CBS (mpCBS)...

In aroma.core: Core Methods and Classes Used by 'aroma.*' Packages Part of the Aroma Framework

Segment copy numbers using the multi-platform CBS (mpCBS) method

Description

Segment copy numbers using the multi-platform CBS (mpCBS) method of the mpcbs package.

WARNING: The mpcbs package is an old package that is no longer maintained. It also has '_R_CHECK_LENGTH_1_CONDITION_' and '_R_CHECK_LENGTH_1_LOGIC2_' bugs, which give errors in R (>= 4.2.0). This means that 'segmentByMPCBS()' does not work in R (>= 4.2.0).

Usage

## S3 method for class 'RawGenomicSignals'
segmentByMPCBS(this, ..., cache=FALSE, force=FALSE, verbose=FALSE)

Arguments

...	Additional arguments passed to the segmentation function.
verbose	See Verbose.

Details

Internally mpcbs.mbic() of the mpcbs package is used for segmenting the signals. This segmentation method does not support weighted segmentation.

Value

Returns the fit object.

Author(s)

Henrik Bengtsson

See Also

For more information see RawGenomicSignals.

Examples

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Simulating copy-number data from multiple platforms
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Piecewise-constant copy-number state function
cnState <- function(x, ...) {
  n <- length(x)
  mu <- double(n)
  mu[20e6 <= x & x <= 30e6] <- +1
  mu[65e6 <= x & x <= 80e6] <- -1
  mu
} # cnState()

xMax <- 100e6

Js <- c(200, 400, 100)
bs <- c(1, 1.4, 0.5)
as <- c(0, +0.5, -0.5)
sds <- c(0.5, 0.3, 0.8)

cnList <- list()
for (kk in seq_along(Js)) {
  J <- Js[kk]
  a <- as[kk]
  b <- bs[kk]
  sd <- sds[kk]
  x <- sort(runif(J, max=xMax))
  mu <- a + b * cnState(x)
  eps <- rnorm(J, sd=sd)
  y <- mu + eps
  cn <- RawCopyNumbers(y, x)
  cnList[[kk]] <- cn
}


# - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Merge platform data (record their sources in 'id')
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
cn <- Reduce(append, cnList)
plot(cn, ylim=c(-3,3), col=cn$id)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Segment
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
legend <- c()

if (require("DNAcopy")) {
  fit <- segmentByCBS(cn)
  cnr <- extractCopyNumberRegions(fit)
  print(cnr)
  drawLevels(cnr, col="white", lwd=6)
  drawLevels(cnr, col="red", lwd=3)
  legend <- c(legend, red="CBS")
}


## WORKAROUND: There's a _R_CHECK_LENGTH_1_LOGIC2_ bug in
## mpcbs::mpcbs.mbic().  Until fixed, if ever, we cannot
## call segmentByMPCBS() here. /HB 2022-11-10
if (isTRUE(Sys.getenv("R_CHECK_FULL")) && require("mpcbs")) {
  fit <- segmentByMPCBS(cn)
  cnr <- extractCopyNumberRegions(fit)
  print(cnr)
  drawLevels(cnr, col="white", lwd=6)
  drawLevels(cnr, col="blue", lwd=3)
  legend <- c(legend, blue="MPCBS")
}


if (length(legend) > 0) {
  legend("topleft", pch=19, col=names(legend), legend, bty="n", horiz=TRUE)
}

aroma.core documentation built on June 25, 2024, 1:15 a.m.