[BioC] limma: cannot repoduce older analysis

[Philipp Pagel]

	Dear list,

About 3 months ago I analyzed a simple two-color array experiment and got
results that looked quite reasonable and biologically sound. For some reason I
wanted to repeat the analysis and add a few plots that I had not included
before.

When I got VERY different results in my toptable, I assumed I must have
changed something in my approach so I simply ran my original analysis
script again and found I was unable to reproduce the original toptable.
I have spent quite some time trying to debug the problem and have to say
that I am stuck. I have the original data files and the original
R-script. The normalization is 100% reproducible - i.e. the normalized
MALists seem to be identical. Yet when searching for differential
expression I get totally different results.

The only difference between the two runs lies in updates to R and limma in
the meantime. Unfortunately, I did not record which version of R, limma etc. I
had used, originally. My current environment is this:

	> sessionInfo()
	R version 2.7.1 (2008-06-23)
	x86_64-pc-linux-gnu

	locale:
	LC_CTYPE=en_US.utf8;LC_NUMERIC=C;LC_TIME=en_US.utf8;LC_COLLATE=en_US.utf8;LC_MONETARY=C;LC_MESSAGES=en_US.utf8;LC_PAPER=en_US.utf8;LC_NAME=C;LC_ADDRESS=C;LC_TELEPHONE=C;LC_MEASUREMENT=en_US.utf8;LC_IDENTIFICATION=C

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

	other attached packages:
	[1] statmod_1.3.6   MASS_7.2-42     xtable_1.5-2    limma_2.14.2    lattice_0.17-10
	[6] cairoDevice_2.8

	loaded via a namespace (and not attached):
	[1] grid_2.7.1  tools_2.7.1


My search for differential expression seems pretty standard to me:

	MA$design <- modelMatrix(targets, ref="control")
	# flag out controls etc.
	MA$weights[MA$genes$Status != 'miRNA', ] = 0.0
	# sort spots by ID to put replicates next to each other
	MA2 <- MA[order(MA$genes$ID), ]

	dupfit <- duplicateCorrelation(MA2, ndups=4)
	fit <- lmFit(MA2, ndups=4, correlation=dupfit$consensus)
	fit <- eBayes(fit)
	tt <- topTable(fit, number=100)

I have siftet through the changelog of limma hoping to find a hint about
some changed default or behaviour in lmFit or eBayes but saw nothing
that seemed to expain my problem.

Any hints apprechiated.

cu
	Philipp

-- 
Dr. Philipp Pagel
Lehrstuhl für Genomorientierte Bioinformatik
Technische Universität München
Wissenschaftszentrum Weihenstephan
85350 Freising, Germany
http://mips.gsf.de/staff/pagel