Search
Question: Design matrix for continuous independent variable in limma
1
2.4 years ago by
Linda Chaba10
Kenya
Linda Chaba10 wrote:

I'm using limma for the first time. I have a problem in creating a design matrix. I'm interested in finding genes that are differentially expressed given only one continuous outcome for example age. I have 35 samples. How will my design matrix looks like?

Thank you.

modified 2.4 years ago • written 2.4 years ago by Linda Chaba10
3
2.4 years ago by
Aaron Lun19k
Cambridge, United Kingdom
Aaron Lun19k wrote:

I would probably do something like this:

age <- runif(35, 20, 50) # dummy variable
require(splines)
X <- ns(age, df=5)
design <- model.matrix(~X)

This will fit a spline to the expression of each gene, using the age of each sample as a covariate. The use of a spline is more flexible than just putting age directly into the design matrix, as it will allow for non-linear trends. Splines with more df provide more flexible fits at the cost of burning up the residual d.f. available for variance estimation. In general, 3-5 d.f. work well as it's rare to get a trend with lots of peaks/troughs in real data. In your case, you've got plenty of samples and age is the only covariate, so you using 5 d.f. won't do too much damage to your variance estimates.

For testing, you can't interpret the spline coefficients individually. Rather, you'll have to drop them all at once:

# Assume we get a fit object out of lmFit/eBayes.
results <- topTable(fit, coef=2:ncol(design), n=Inf, sort.by="none")

This will test for any response of expression to age, be it an increasing/decreasing/squiggly trend. You'll have to look at the results on a gene-by-gene basis to determine what the trend is, because that's not easily summarized into a single statistic. For purposes of seeing whether it's generally up or down with age, you may consider fitting a simpler model using age directly as a covariate in design:

design2 <- model.matrix(~age)
# ... after processing with lmFit/eBayes to get fit2 ...
results2 <- topTable(fit2, coef=2, n=Inf, sort.by="none")

... and reporting the log-fold change from results2 with the p-values in results. Positive log-fold changes correspond to a positive gradient with respect to age (i.e., expression goes up with age). I wouldn't use the p-values from results2, as it assumes linearity.

Thank you Aaron for your answer, it was very helpful. However, I'm getting an error (Error: could not find function "topTags") when using "topTags" to get my results. It could be because I'm missing a package. What difference would it have if I use "topTable" to get my results?  I am using limma package.

Thank you

2

Aaron got his wires crossed a bit, is all. topTags is the analogous function to topTable in edgeR. You should go forward with using topTable

Thank you Steve.

Thanks Steve, fixed above.

Most people like to see which genes are more significant, so allow topTable() to sort the results for you:

topTable(fit, n=20)

say.

The reason why I set sort.by="none" was to allow for easier cross-referencing between results and results2, if the gradient of the age response was of interest. Otherwise, we'd have to use match or merge.

Dear Aaron,

One more question.

I don't have raw data set.I have  processed (normalized and filtered) data set consisting of log2 rations of 23,360 genes on 35 samples in an excel format. I imported the data the usual way and used it to get a fit object out of lmFit/eBayes. My question is, is it okay to use processed data since I see most examples using raw data set. Secondly, with the results2 you mentioned, can I use adjusted p-values to get differentially expressed genes?  I'm adjusting using fdr approach.

Thank you

1
1. This depends on how the data is normalized and filtered, but generally speaking, yes. lmFit expects log-transformed expression values after normalization and filtering - see Chapters 6 and 7 of the limma user's guide.
2. Use the adjusted p-values from results, not results2. The only purpose of results2 is to give you a single log-fold change statistic that you can interpret as the gradient with respect to age. The fit with the splines will provide more appropriate (adjusted) p-values in results as it doesn't assume linearity in expression with respect to age.