<%@meta title="RMA reproducibility: justRMA() versus affy::justRMA()"%>
<%@meta author="Henrik Bengtsson, Mark Robinson"%>

# <%@meta name="title"%>

by <%@meta name="author"%>

This test verifies that `aroma.affymetrix::justRMA()` of aroma.affymetrix can reproduce the RMA chip-effect estimates as estimated by `affy::justRMA()`.

<%
library("R.utils")
use("R.devices")
verbose <- Arguments$getVerbose(-8, timestamp=TRUE)
options(digits=4L)
plotDensity <- aroma.light::plotDensity
pairTag <- "justRMA_vs_justRMA"
%>

## Setup
```r
<%=withCapture({
use("aroma.affymetrix")
use("affy")
use("hgu133plus2cdf")
})%>
```

## Data
```r
<%=withCapture({
dataset <- "GSE9890"
chiptype <- "HG-U133_Plus_2"
cdf <- AffymetrixCdfFile$byChipType(chiptype)
cdf

csR <- AffymetrixCelSet$byName(dataset, cdf=cdf)
csR
})%>
```

## RMA by aroma.affymetrix
```r
<%=withCapture({
eset <- justRMA(csR, verbose=verbose)
print(eset)
})%>
```

## RMA by affy
```r
<%=withCapture({
eset0 <- justRMA(filenames=getPathnames(csR))
print(eset0)
})%>
```

<%---------------------------------------------------------------
 Assertions
 ---------------------------------------------------------------%>
<%
# Sample names
sampleNames0 <- sampleNames(eset0)
sampleNames <- sampleNames(eset)
stopifnot(identical(sampleNames, tools::file_path_sans_ext(sampleNames0)))

# Feature names
stopifnot(identical(featureNames(eset), featureNames(eset0)))

# Annotation
stopifnot(identical(annotation(eset), annotation(eset0))) 
%>


## Comparisons

Overall correlation of chip effects per each sample:
```r
<%=withCapture({
theta0 <- exprs(eset0)
theta <- exprs(eset)
rho <- diag(cor(theta, theta0))
rho
range(rho)
})%>
```

Chip-effect differences:
```r
<%=withCapture({
d <- (theta - theta0)
summary(d)
})%>
```

![](<%=toPNG(dataset, tags=c(chiptype, pairTag, "pairwise"), aspectRatio=5/2, {
  par(mar=c(5,5,4,2)+0.1, cex.main=2, cex.lab=2, cex.axis=1.5)
  layout(matrix(seq_along(csR), ncol=2L, byrow=TRUE))

  xlab <- expression(log[2](theta[affyPLM]))
  ylab <- expression(log[2](theta[aroma.affymetrix]))
  for (kk in seq_len(ncol(theta))) {
    main <- colnames(theta)[kk]
    plot(theta0[,kk], theta[,kk], pch=".", xlab=xlab, ylab=ylab, main=main)
    abline(0,1, col="blue")
    stext(side=3, pos=0, line=-1.1, cex=1.2, substitute(rho==x, list(x=rho[kk])))
  }
})%>)

_Figure: Pairwise chip-effect estimates (theta) from `aroma.affymetrix::justRMA()` and `affy::justRMA()` across <%=length(csR)%> <%=chiptype%> arrays._



![](<%=toPNG(dataset, tags=c(chiptype, pairTag, "density_of_differences"), aspectRatio=2/3, {
  par(mar=c(5,5,4,2)+0.1, cex.main=2, cex.lab=2, cex.axis=1.5)
  xlab <- expression(log[2](theta[aroma.affymetrix]/theta[affyPLM]))
  plotDensity(d, xlab=xlab)
})%>)

_Figure: Empirical density of differences in chip-effect estimates (theta) between `aroma.affymetrix::justRMA()` and `affy::justRMA()` across <%=length(csR)%> <%=chiptype%> arrays._


<%----
<%
# (b) Assert correlations
print(rho)
print(range(rho))
stopifnot(all(rho > 0.99995))

# (c) Assert differences
stopifnot(mean(as.vector(d^2)) < 1e-3)
stopifnot(sd(as.vector(d^2)) < 1e-3)
stopifnot(quantile(abs(d), 0.99) < 0.05)
stopifnot(max(abs(d)) < 0.100)

verbose && print(verbose, sessionInfo())
%>

----%>

## Appendix
```r
<%=withCapture({
sessionInfo()
})%>
```

--------------------------------------------------------------------
Last generated: <%=format(Sys.time(), "%Y-%m-%d %H:%M:%S UTC%z")%>.
Powered by [RSP](http://cran.r-project.org/package=R.rsp).
