Hi, 

I'd like to use limma to look for probes that correlation strongly and with a large slope (positive or negative), based on a linear variable. I've scanned the Limma users guide but couldn't find such an example. 

design              <- model.matrix(~0 + df$AGE)
fit                 <- lmFit(exprs(lumi.norm), design)
fit                 <- eBayes(fit)
topTable(fit, adjust.method="BH", number=10, p.value=0.01, sort.by="M")

The above code shows what I've done so far, which is provide the design of the model with a linear variable. I was wondering if this was the correct way to do it? If so, how would i filter for strength of correlation of size of slope?

Is there a better way to do it?

Thanks!

If you are going to fit a conventional linear model, you should not constrain the intercept to zero (which is what you are doing). If you were to do such a thing, you are in effect saying that the expression for all of your genes at birth should be equal to zero, which is only true for zombie babies.

Instead you should do

design <- model.matrix(~AGE, df)

and

topTable(fit, 2, <other args you like>)

And then your logFC column will be the slope of the line, correlating age and gene expression. The p-value will test the hypotheses Ho: slope == 0 vs Ha: slope != 0.

But do note that this linear model assumes that the expression is a linear function of age, which may not be true, especially if you have a wide range of ages. In which case you might want a quadratic term as well, but that can make interpretation tricky.

As for 'strength of correlation of size of slope', I am not sure that's a thing. You can test that the slope is not equal to zero, or you could use treat() and require that the slope be greater than some value that you think is the limit of biological meaningfulness (which seems like it shouldn't be a word, yet is). But in the microarray context, there isn't much else you can do.

I suppose you could look at individual genes using lm() and associated model diagnostics, but I would normally wait for the validation stage for that sort of thing.

The main issue with your code is that you eliminated the necessary intercept term by putting the zero in your model. That's essentially forcing all your lines to go through the origin (i.e. zero RMA value at age zero), which makes no biological sense. You should probably just use a design of model.matrix(~1+AGE, data=df). Secondly, if you want to enrich your results for slopes that are both large and significant, then you should experiment with the treat function, which combines a fold change (i.e. slope) cutoff with a significance cutoff, instead of topTable, which simply tests whether the slope is significantly different from zero.