[BioC] Sweave & limma mystery

[Matthias Kohl]

Hi Philipp,

to generate your plot on different devices, the code in your document 
has to be executed more than once. Hence, the modification of your 
object MA works as expected in the first run, but also has an influence 
on further runs. This is not a specific problem which occurs only with 
plotMA.

I would recommend to introduce

MA0 <- MA
MA0$weights[MA0$genes$Status != 'gene', ] <- 0.0

and put this in the code chunk before the chunk which includes the 
plotting (and of course, use MA0 for the third and fourth plot).

Best,
Matthias


Philipp Pagel wrote:
> 	Hi!
>
> I have a pretty odd (and lengthly to describe) problem involving plotMA
> (from package limma) and Sweave: Inside a Sweave document I want to plot
> 4 MA-plots in one plot device. After the first two plots I set the
> weight of all control spots to zero in order remove them from the plot. 
> To my surprise, the changed weights appear to have a magic effect on the
> two first plots as well when judged by the resulting pdf file. BUT: the
> plots Sweave generates in the X11 device in parallel (why does it do
> that in the first place?) look fine.
>
> OK - that was quite abstract and hard to understand so here is the
> smallest example I was able to put together:
>
> I use the Bob.RData file found here:
> http://bioinf.wehi.edu.au/limma/data/Bob.RData
>
> SpotTypes.txt looks like this:
>
> ----------------------snip----------------------------------
> SpotType        ID      Library Color
> gene    *       *       black
> control *       Control red
> ----------------------snip----------------------------------
>
> And here is my report.Snw file:
>
> ----------------------snip----------------------------------
> \documentclass[11pt]{scrartcl}
> \begin{document}
> <<echo=F>>=
> 	require(limma)
> 	load('Bob.RData')
> 	spottypes = readSpotTypes()
> 	RG$genes$Status = controlStatus(spottypes, RG)
> 	RG$weights <- matrix( rep(1, dim(RG)[1]*dim(RG)[2]), ncol=dim(RG)[2] )
> @
>
> Some text here\dots
>
> <<>>=
> 	MA <- normalizeWithinArrays(RG, bc.method='normexp', method="robustspline", offset=50)
> 	MA <- normalizeBetweenArrays(MA, method="scale")
> @
>
> Even more text here\dots
>
> <<fig=T, echo=F>>=
> 	oldpar <- par(no.readonly = TRUE)
> 	par(las=1, cex=0.7, mfrow=c(2,2), oma=c(0,0,0,0))
>
> 	plotMA(MA, array=1, legend=F)
> 	plotMA(MA, array=2, legend=F)
> 	# hide all control spots
> 	MA$weights[MA$genes$Status != 'gene', ] <- 0.0
> 	plotMA(MA, array=1, legend=F)
> 	plotMA(MA, array=2, legend=F)
> 	par(oldpar)
> @
> \end {document}
> ----------------------snip----------------------------------
>
> in R:
>
>   
>> Sweave('report.Snw')
>>     
>
> and then
>
> pdflatex report.tex
>
> My R version/platform:
>
> platform       x86_64-pc-linux-gnu
> arch           x86_64
> os             linux-gnu
> system         x86_64, linux-gnu
> status
> major          2
> minor          7.0
> year           2008
> month          04
> day            22
> svn rev        45424
> language       R
> version.string R version 2.7.0 (2008-04-22)
>
> limma Version: 2.14.2 Date: 2008/05/22
>
> To see what I mean compare the plot in report.pdf to the screen output:
> On the screen the first row has red control spots and the second doesn't
> -- in the report no control spots are plotted in any of the 4 plots.
>
> Can anyone reproduce this? If so - any ideas how/why this happens?
>
> cu
> 	Philipp
>
>   

-- 
Dr. Matthias Kohl
www.stamats.de