[BioC] vsn and scaling

[Wolfgang Huber]

Dear Roslin,

I have now had a look at your data (which you sent me off-list). The 
quick bottomline for hurried listreaders is: consider the parameter 
"lts.quantile" and try different choices from the default, especially if 
you are looking at quite variable data!


Now in detail. There are two issues:

1. your two datasets not only differ in the order of the arrays 
(columns), but also in the set of features (rows of the data matrix): 
the first one has 3884, the second has 3746. Out of these, 3588 are in 
common,  296 are only in the first and 158 are only in the second.

2. the data appear relatively noisy, and hence the precision of the 
estimated normalization parameters is not great.

For the purpose of comparison, (1) can be addressed by fitting vsn only 
to the common features (you can then apply the fit to all using the 
predict method).

(2) can often be addressed by varying the parameter "lts.quantile". Its 
defafult value 0.9 in vsn2 is probably too optimistic for your data (it 
assumes less variability between arrays), and in the example code below 
I set it to 0.5. With that, the results for R10001 and R10002 are very 
similar (see commonfeatures.png).

The output of vsn2 can depend sensitively on the input (e.g. the order 
of the arrays or the addition/omission of a few rows) BUT if that 
happens, it indicates a problem, either because the vsn model is not 
appropriate for the data or the numerical optimisation algorithm failed 
to find the optimal fit. In such cases, tuning of the parameters (in 
particular, lts.quantile and pstart) can help a lot.

Then, in some cases (and I think this applies to yours), the results can 
seemingly differ, but they could both be within the range of the 
plausible. Think of them as two different estimates for the same 
underlying parameter, which are both consistent with the available data 
and the model. This can happen if the likelihood landscape is very 
shallow / the model is overparameterized. In particular, this will 
happen for vsn2 when the relative size of additive noise is small, 
"glog" is overkill, and "log" would do just as well.

The resulting plots ans script are also here: 
http://www.ebi.ac.uk/~huber/pub/roslin

I hope this helps, please let me know.

Best wishes
	Wolfgang

---------------------------------------------------------------
library("vsn")

x1 = as.matrix(read.table("R10001.vsn.filtered.raw.txt", sep="\t", 
quote="", header=TRUE, row.names=1L))
x2 = as.matrix(read.table("R10002.vsn.filtered.raw.txt", sep="\t", 
quote="", header=TRUE, row.names=1L))

common = intersect(rownames(x1), rownames(x2))
onlyx1 = setdiff(rownames(x1), rownames(x2))
onlyx2 = setdiff(rownames(x2), rownames(x1))

cat(length(common), length(onlyx1), length(onlyx2), "\n")
## 3588 296 158

myPlot = function(a, b, name) {
   png(paste(name, "png", sep="."), width=700, height=500)
   par(mfrow=c(2,3))
   for(cn in colnames(a)){
     plot(a[common, cn], b[common, cn], pch=".", main=cn)
     abline(a=0, b=1, col="blue")
   }
   dev.off()
}


myPlot(log2(x1), log2(x2), "beforevsn")

## with all features
v1 = vsn2(x1, lts.quantile=0.5)
v2 = vsn2(x2, lts.quantile=0.5)
myPlot(exprs(v1), exprs(v2), "withallfeatures")

## with common features only
w1 = vsn2(x1[common,], lts.quantile=0.5)
w2 = vsn2(x2[common,], lts.quantile=0.5)
myPlot(exprs(w1), exprs(w2), "commonfeatures")

-------------------------------------------------------------------
 > sessionInfo()
R version 2.6.0 Patched (2007-10-03 r43075)
i686-pc-linux-gnu

locale:
LC_CTYPE=en_GB.UTF-8;LC_NUMERIC=C;LC_TIME=en_GB.UTF-8;LC_COLLATE=en_GB.UTF-8;LC_MONETARY=en_GB.UTF-8;LC_MESSAGES=en_GB.UTF-8;LC_PAPER=en_GB.UTF-8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_GB.UTF-8;LC_IDENTIFICATION=C

attached base packages:
[1] tools     stats     graphics  grDevices utils     datasets  methods
[8] base

other attached packages:
[1] vsn_3.3.1              limma_2.11.14          affy_1.15.12
[4] preprocessCore_0.99.22 affyio_1.5.11          Biobase_1.15.36
[7] fortunes_1.3-3

loaded via a namespace (and not attached):
[1] grid_2.6.0      lattice_0.16-5  rcompgen_0.1-15
-------------------------------------------------------------------



Wolfgang Huber ha scritto:
> Dear Roslin,
> 
>> Am I correct in saying that the default setting of vsn2 is to scale 
>> everything to the first array?
> 
> This is not the case. The model that vsn fits is invariant under 
> permutation of the arrays (you can look at the vignette "Likelihood 
> Calculations for vsn" if you are interested in details).
> 
> However, the fitting does includes an iterative numeric optimisation; if 
> that fails, and runs into a local optimum, such a behaviour can depend 
> sensitively on the data (not just the order of the arrays).
> 
>> Could this be potentially problematic? For example, we have observed 
>> strange results when we change the order of loading of the data into 
>> vsn2. The resulting normalised data looks different when vieved in an 
>> MA plot.
> 
> Can you send me an example that I can reproduce? Such examples would be 
> helpful for me to improve the diagnostics that can be done by vsn2, 
> either automatically, or in a subsequent quality control step.
> 
> Best wishes
>   Wolfgang
> 
> ------------------------------------------------------------------
> Wolfgang Huber  EBI/EMBL  Cambridge UK  http://www.ebi.ac.uk/huber
>