Nothing
library("aroma.affymetrix")
log <- Arguments$getVerbose(-4, timestamp=TRUE)
dataSet <- "HapMap,CEU,testset"
chipTypes <- c("Mapping50K_Hind240", "Mapping50K_Xba240")
# Expected sample names
sampleNames <- c("NA06985", "NA06991", "NA06993",
"NA06994", "NA07000", "NA07019")
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Tests for setting up CEL sets and locating the CDF file
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
csRList <- list()
for (chipType in chipTypes) {
cs <- AffymetrixCelSet$byName(dataSet, chipType=chipType)
print(cs)
stopifnot(identical(getNames(cs), sampleNames))
csRList[[chipType]] <- cs
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Probe-level modelling test #1
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Do it on the raw data so that we can validate the results
csList <- csRList;
plmList <- list()
for (chipType in names(csList)) {
cs <- csList[[chipType]]
plm <- RmaPlm(cs)
print(plm)
fit(plm, verbose=log)
ces <- getChipEffectSet(plm)
print(ces)
stopifnot(identical(getNames(ces), getNames(cs)))
# Assert a few results
if (chipType == "Mapping250K_Nsp") {
unit <- 100
theta <- readUnits(ces, units=unit, verbose=log)
data <- unlist(theta, use.names=FALSE)
data <- sum(data, na.rm=TRUE)
if (!isZero(data - 57448.85473633, eps=.Machine$float.eps)) {
print(theta)
throw("Non-reproducible chip-effect estimates detected for unit ", unit)
}
}
plmList[[chipType]] <- plm
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# PLM residual test
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
rsList <- list()
for (chipType in names(plmList)) {
plm <- plmList[[chipType]]
print(plm)
rs <- calculateResidualSet(plm, verbose=log)
print(rs)
stopifnot(identical(getNames(rs), getNames(cs)))
# Assert a few results
if (chipType == "Mapping50K_Hind240") {
cs <- getDataSet(plm)
ces <- getChipEffectSet(plm)
paf <- getProbeAffinityFile(plm)
unit <- 500
y <- readUnits(plm, units=unit, verbose=log)
theta <- readUnits(ces, units=unit, verbose=log)
phi <- readUnits(paf, units=unit, verbose=log)
resids <- readUnits(rs, units=unit, verbose=log)
str(resids)
}
rsList[[chipType]] <- rs
}
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# Spatial residual plots test
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
ae <- ArrayExplorer(rsList)
setColorMaps(ae, c("log2,log2neg,rainbow", "log2,log2pos,rainbow"))
print(ae)
stopifnot(identical(unname(getArrays(ae)), getNames(cs)))
process(ae, arrays=1:2, interleaved="auto", verbose=log)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
# PLM weights (from residuals) test
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
for (chipType in names(plmList)[1]) {
plm <- plmList[[chipType]]
print(plm)
rs <- getResidualSet(plm, verbose=log)
print(rs)
stopifnot(identical(getNames(rs), getNames(cs)))
ws <- calculateWeights(plm, verbose=log)
print(rs)
rsList[[chipType]] <- rs
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.