I'm using limma in BioConductor to make a simple paired comparison in gene expression.  Now I'd like to run an auxiliary analysis to see if there are individual differences which correlate with changes in gene expression.  I'm having difficulty understanding how to fold this additional research question in.

First - the simple paired comparison.  I have animals which learn a task on only 1 side of the body.  I then harvest CNS tissue from the trained and untrained sides.  Gene expression due to learning is then identified by comparing the trained and untrained samples on a two-color array (each array is within-subjects with one dye for the trained sample, the other for the untrained sample; dyes are counterbalanced across condition; 16 biological replicates (i.e. 16 different microarrays, each with trained/untrained samples from 1 animal)).  The essential part of the R script I'm using is:

design <- modelMatrix(targets, ref="Control")
fit <- lmFit(MA2.avg, design)
regulated <- treat(fit, lfc=0.137503524)

Even though all the animals have learned to criterion, there are large individual differences in their ability to retain the information.  So each animal has a 'memory score' that measures the degree to which the learned behavior has decayed between training and harvesting the samples. 

So: I'd like to find changes in gene expression that relate to memory scores.  To be clear, I'm interested in relating memory scores to the log-fold-change scores (the comparisons of trained/untrained), though I suppose it could also be possible to look for correlations with expression overall (e.g. the average expression across both sides).

I know that limma can handle quantitative conditions.  But none of the examples/papers I've been able to find so far seem to quite fit the design I have.  Hope I'm not asking an obvious our previously answered question. 

Any help or pointers much appreciated,

Bob

Just to be clear, here's the targets array for the paired comparison.

SampleNumber    FileName    Cy3    Cy5
1   Animal1.txt    Control    Trained
...
16    Animal16.txt    Trained    Control

So I could extend it with the memory scores:

SampleNumber    FileName    Cy3    Cy5    MScore
1    Animal1.txt    Control    Trained    2.1
...
16    Animal168.txt    Trained    Control    1.3

but then I'm pretty lost as to how to make a proper contrast to find transcripts which have a LFC correlated with the memory scores.

First, I strongly recommend you include a probe-specific dye effect in the design matrix:

design <- cbind(Dye=1, design)

This exploits your dye-swap design and often improves the analysis markedly for Agilent arrays.

For the MScore, you need to make the sign of this variable match the sign of the M-values:

i <- targets$Cy5=="Control"
MScore[i] <- -MScore[i]

Then just add the MScore to the design matrix:

design <- cbind(design, MScore=Mscore)

Then test for that coefficient:

fit <- lmFit(MA2.avg, design)
fit <- treat(fit, lfc=log2(1.1))
topTreat(fit, coef="MScore")