backgroundBiasEstimation: Background bias estimation

In benja0x40/MRA.TA: Multi Resolution Analysis of Tiling Arrays

Background bias estimation

Description

Estimate a background bias within A, M values from a tiling array.

Usage

backgroundBiasEstimation(
  x,
  y,
  AM.scale.compensation = T,
  smoothness = 0.08,
  epsilon = 0.001,
  nsteps = 11,
  plots = F,
  xlab = "A",
  ylab = "M"
)

Arguments

x	numeric vector of A values, e.g. the average of log2 of the specific and the reference signals.
y	numeric vector of M values, e.g. the log2 ratio of the specific over the reference signal.
AM.scale.compensation	logical determining if A and M values should be rescaled prior to bias estimation. Scale compensation is required and activated by default.
smoothness	numeric coefficient controlling the meanshift radius in the rescaled (A,M) plane.
epsilon	convergence limit corresponding to the minimum displacement required per meanshift iteration in the rescaled (A,M) plane.
nsteps	minimum number of meanshift iterations.
plots	logical used to activate control plots (default = F, no plots)
xlab	label of the x axis for control plots.
ylab	label of the y axis for control plots.

Value

numeric value of the estimated bias angle in radians

Author(s)

Benjamin Leblanc

References

Leblanc B., Comet I., Bantignies F., and Cavalli G., Chromosome Conformation Capture on Chip (4C): data processing. Book chapter in Polycomb Group Proteins: Methods and Protocols. Lanzuolo C., Bodega B. editors, Methods in Molecular Biology (2016). http://dx.doi.org/10.1007/978-1-4939-6380-5_21

See Also

backgroundBiasCorrection, normalizeArrayData

Examples

# A simple test of the accuracy of bias estimations
ntst <- 10 # Increase to over 1000 for a reliable test of accuracy

rotate <- function(x, y, theta) {
  u <- x * cos(theta) - y * sin(theta)
  v <- x * sin(theta) + y * cos(theta)
  list(x=u, y=v)
}

tst <- c()
plots <- T
pb <- txtProgressBar(min=1, max=ntst, char="|", style=3)
for(i in 1:ntst) {
  n.bg <- sample(c(70000, 80000, 90000), 1) # Background level data
  n.sp <- 100000 - n.bg                     # Enriched level data
  a <- c(rnorm(n.bg, 10, 2),  rnorm(n.sp, 10, 1))
  m <- c(rnorm(n.bg, 0,  0.5), rnorm(n.sp,  1.5, 1))
  alpha <- runif(1, -pi/4, pi/4)            # Pick a random angle
  r <- rotate(a, m, theta = alpha)          # Simulate bias
  xtm <- Sys.time()
  theta <- backgroundBiasEstimation(        # Estimate bias
    r$x, r$y, plots=plots, AM.scale.compensation = F, smoothness = 0.07
  )
  xtm <- Sys.time() - xtm
  alpha <- 180/pi * alpha
  theta <- 180/pi * theta
  tst <- rbind(
    tst, c(n.bg=n.bg, n.sp=n.sp, time=xtm, simulated=alpha, estimated=theta)
  )
  setTxtProgressBar(pb, i)
  plots <- F
}
close(pb)
write.table(
  tst, "testBackgroundBiasEstimation.txt",
  quote = F, sep = '\t', row.names=F, col.names=T
)

tst <- read.delim("testBackgroundBiasEstimation.txt", stringsAsFactors=F)
boxplot(tst$estimated - tst$simulated, range=0, ylab="estimated - simulated (degrees)")
legend("topright", paste("N =", ntst), bty="n")

# Result with ntst=1000, delta <2.5 degrees in 95% of the tests:
# delta <- abs(tst$estimated - tst$simulated)
# q <- quantile(delta, probs=c(0.5, 0.75, 0.9, 0.95, 0.99), na.rm=T)

benja0x40/MRA.TA documentation built on March 13, 2023, 5:15 a.m.