| normalizeFragmentLength {aroma.light} | R Documentation |
Normalizes signals for PCR fragment-length effects. Some or all signals are used to estimated the normalization function. All signals are normalized.
## Default S3 method: normalizeFragmentLength(y, fragmentLengths, targetFcns=NULL, subsetToFit=NULL, .isLogged=TRUE, ..., .returnFit=FALSE)
y |
A numeric vector of length K of signals to be normalized
across E enzymes. |
fragmentLengths |
An integer KxE matrix of fragment lengths. |
targetFcns |
A list of E functions - one per enzyme. |
subsetToFit |
The subset of data points used to fit the
normalization function.
If NULL, all data points are considered. |
.isLogged |
A logical. |
... |
Additional arguments passed to lowess. |
.returnFit |
A logical. |
Returns a numeric vector of the normalized signals.
It is assumed that the fragment-length effects from multiple enzymes added (with equal weights) on the intensity scale. The fragment-length effects are fitted for each enzyme separately based on units that are exclusively for that enzyme. If there are no or very such units for an enzyme, the assumptions of the model are not met and the fit will fail with an error. Then, from the above single-enzyme fits the average effect across enzymes is the calculated for each unit that is on multiple enzymes.
Henrik Bengtsson (http://www.braju.com/R/)
[1] H. Bengtsson, R. Irizarry, B. Carvalho, and T.P. Speed. Estimation and assessment of raw copy numbers at the single locus level, Bioinformatics, 2008.
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Example 1: Single-enzyme fragment-length normalization of 6 arrays
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Number samples
I <- 9;
# Number of loci
J <- 1000;
# Fragment lengths
fl <- seq(from=100, to=1000, length.out=J);
# Simulate data
y <- matrix(0, nrow=J, ncol=I);
maxY <- 12;
for (kk in 1:I) {
k <- runif(n=1, min=3, max=5);
mu <- function(fl) {
maxY - fl^{1/k};
}
eps <- rnorm(J, mean=0, sd=1);
y[,kk] <- mu(fl) + eps;
}
# Normalize data (to a zero baseline)
yN <- apply(y, MARGIN=2, FUN=function(y) {
normalizeFragmentLength(y, fragmentLengths=fl);
})
# Plot raw data
layout(matrix(1:9, ncol=3));
xlim <- c(0,max(fl));
ylim <- c(0,max(y));
xlab <- "Fragment length";
ylab <- expression(log2(theta));
for (kk in 1:I) {
plot(fl, y[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab);
lines(lowess(fl, y[,kk]), col="red", lwd=2);
}
# Plot normalized data
layout(matrix(1:9, ncol=3));
ylim <- c(-1,1)*max(y)/2;
for (kk in 1:I) {
plot(fl, yN[,kk], xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab);
lines(lowess(fl, yN[,kk]), col="blue", lwd=2);
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Example 2: Two-enzyme fragment-length normalization of 6 arrays
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
set.seed(0xbeef);
# Number samples
I <- 5;
# Number of loci
J <- 2000;
# Fragment lengths (two enzymes)
fl <- matrix(0, nrow=J, ncol=2);
fl[,1] <- seq(from=100, to=1000, length.out=J);
fl[,2] <- seq(from=1000, to=100, length.out=J);
# Let 1/2 of the units be on both enzymes
fl[seq(from=1, to=J, by=4),1] <- NA;
fl[seq(from=2, to=J, by=4),2] <- NA;
# Sty/Nsp mixing proportions:
rho <- rep(1, I);
rho[1] <- 1/3; # Less Sty in 1st sample
rho[3] <- 3/2; # More Sty in 3rd sample
# Simulate data
z <- array(0, dim=c(J,2,I));
maxLog2Theta <- 12;
for (ii in 1:I) {
# Common effect for both enzymes
mu <- function(fl) {
k <- runif(n=1, min=3, max=5);
maxLog2Theta - fl^{1/k};
}
# Calculate the effect for each data point
for (ee in 1:2) {
z[,ee,ii] <- mu(fl[,ee]);
}
# Update the Sty/Nsp mixing proportions
ee <- 2;
z[,ee,ii] <- rho[ii]*z[,ee,ii];
# Add random errors
for (ee in 1:2) {
eps <- rnorm(J, mean=0, sd=1/sqrt(2));
z[,ee,ii] <- z[,ee,ii] + eps;
}
}
hasFl <- is.finite(fl);
unitSets <- list(
nsp = which( hasFl[,1] & !hasFl[,2]),
sty = which(!hasFl[,1] & hasFl[,2]),
both = which( hasFl[,1] & hasFl[,2])
)
# The observed data is a mix of two enzymes
theta <- matrix(NA, nrow=J, ncol=I);
# Single-enzyme units
for (ee in 1:2) {
uu <- unitSets[[ee]];
theta[uu,] <- 2^z[uu,ee,];
}
# Both-enzyme units (sum on intensity scale)
uu <- unitSets$both;
theta[uu,] <- (2^z[uu,1,]+2^z[uu,2,])/2;
# Calculate target array
thetaT <- rowMeans(theta, na.rm=TRUE);
targetFcns <- list();
for (ee in 1:2) {
uu <- unitSets[[ee]];
fit <- lowess(fl[uu,ee], log2(thetaT[uu]));
class(fit) <- "lowess";
targetFcns[[ee]] <- function(fl, ...) {
predict(fit, newdata=fl);
}
}
# Normalize data (to a target baseline)
thetaN <- matrix(NA, nrow=J, ncol=I);
fits <- vector("list", I);
for (ii in 1:I) {
lthetaNi <- normalizeFragmentLength(log2(theta[,ii]), targetFcns=targetFcns, fragmentLengths=fl, .returnFit=TRUE);
fits[[ii]] <- attr(lthetaNi, "modelFit");
thetaN[,ii] <- 2^lthetaNi;
}
# Plot raw data
xlim <- c(0, max(fl, na.rm=TRUE));
ylim <- c(0, max(log2(theta), na.rm=TRUE));
Mlim <- c(-1,1)*4;
xlab <- "Fragment length";
ylab <- expression(log2(theta));
Mlab <- expression(M==log[2](theta/theta[R]));
layout(matrix(1:(3*I), ncol=I, byrow=TRUE));
for (ii in 1:I) {
plot(NA, xlim=xlim, ylim=ylim, xlab=xlab, ylab=ylab, main="raw");
# Single-enzyme units
for (ee in 1:2) {
# The raw data
uu <- unitSets[[ee]];
points(fl[uu,ee], log2(theta[uu,ii]), col=ee+1);
}
# Both-enzyme units (use fragment-length for enzyme #1)
uu <- unitSets$both;
points(fl[uu,1], log2(theta[uu,ii]), col=3+1);
for (ee in 1:2) {
# The true effects
uu <- unitSets[[ee]];
lines(lowess(fl[uu,ee], log2(theta[uu,ii])), col="black", lwd=4, lty=3);
# The estimated effects
fit <- fits[[ii]][[ee]]$fit;
lines(fit, col="orange", lwd=2);
muT <- targetFcns[[ee]](fl[uu,ee]);
lines(fl[uu,ee], muT, col="purple", lwd=1);
}
}
# Calculate log-ratios
thetaR <- rowMeans(thetaN, na.rm=TRUE);
M <- log2(thetaN/thetaR);
# Plot normalized data
for (ii in 1:I) {
plot(NA, xlim=xlim, ylim=Mlim, xlab=xlab, ylab=Mlab, main="normalized");
# Single-enzyme units
for (ee in 1:2) {
# The normalized data
uu <- unitSets[[ee]];
points(fl[uu,ee], M[uu,ii], col=ee+1);
}
# Both-enzyme units (use fragment-length for enzyme #1)
uu <- unitSets$both;
points(fl[uu,1], M[uu,ii], col=3+1);
}
ylim <- c(0,1.5);
for (ii in 1:I) {
data <- list();
for (ee in 1:2) {
# The normalized data
uu <- unitSets[[ee]];
data[[ee]] <- M[uu,ii];
}
uu <- unitSets$both;
if (length(uu) > 0)
data[[3]] <- M[uu,ii];
plotDensity(data, col=1:3+1, xlim=Mlim, xlab=Mlab, main="normalized");
abline(v=0, lty=2);
}