Limma contrast for variables disease, gender and age
1
1
Entering edit mode
miyuki_iii ▴ 10
@miyuki_iii-11085
Last seen 7.8 years ago

Hi,

I'm a student and I have doubts about how use Limma to contrast the differencial gene expression evaluating my three variables: disease (case/control), gender(male/female) and age.

My data is like this: (three columns with the gender, age and the disease)

> head(pData(gse))
                                                                                          gender age disease
GSM311511_jrs_U201u5_AR2_31_B112_R12558.CEL.gz female  88 control
GSM311528_jrs_U201u5_AR2_33_B112_R92737.CEL.gz    male  53 control
GSM311511_jrs_U201u5_AR2_31_B112_R97522.CEL.gz    male  94 control
GSM311530_jrs_U201u5_AR2_55_B013_R97210.CEL.gz female  40 control
GSM311511_jrs_U201u5_AR2_31_B110_R95417.CEL.gz     male 35 control
GSM311511_jrs_U201u5_AR2_31_B012_R97841.CEL.gz female  89 control

I know build my experiment over my variable disease:
        design= model.matrix(~ pData(gse)[,"disease"])
        colnames(design)= c("control", "case")

        fit= lmFit(gse, fesign)
        fit1= eBayes(fit)

        topTable(fit1, coef=2, adjust="BH", number= as.numeric(dim(gse)[1]))

But, How I can do a contrast Limma with my 3 variables?  Thank for your help!

age= pData(gse)[,"age"]
gender= pData(gse)[,"gender"]
disease= pData(gse)[,"disease"]

age.group= ifelse(age>40,"over40", "under40")

design = model.matrix(~0 + disease + gender + gender:age)

contr.matrix <- makeContrasts(

limma design matrix limma limma contrast matrix microarrays • 3.2k views
ADD COMMENT
6
Entering edit mode
Aaron Lun ★ 28k
@alun
Last seen 3 hours ago
The city by the bay

If you look at the column names of your design matrix, you should get something like this:

[1] "diseasecontrol"   "diseasesick"      "gendermale"       "genderfemale:age"
[5] "gendermale:age"  

(I've assumed that the other level of the disease factor is just called sick.) The first two coefficients represent average expression of in the control or sick condition; the third coefficient is the average log-fold change in males over females; and the last two coefficients represent the gender-specific effect of age.

Now, if you want to compare disease levels, you would do something like:

con <- makeContrasts(diseasesick - diseasecontrol, levels=design)

Then run this through contrasts.fit and use topTable with coef=1.

If you want to test for a gender effect, it's a bit more complicated. You not only need to compare the effect of gender directly, but also the interaction between gender and age. This is because if the age response is different between genders, then I would argue that gender does have an effect, even if the gendermale coefficient is equal to zero. To do this comparison:

con <- makeContrasts(gendermale, genderfemale_age - gendermale_age, levels=design)

Again, run this through contrasts.fit and use topTable with coef=NULL. This setting of coef will test both contrasts in con at once, because it doesn't make sense to consider the first contrast without also considering the second. Note that I've replaced the colons with underscores in the column names of design, otherwise makeContrasts won't be happy.

Finally, if you want to test for an age effect, you can do something like:

con <- makeContrasts(genderfemale_age, gendermale_age, levels=design)

This will test if there is a response to age in either male or female individuals. You can set coef appropriately in topTable to test either contrast separately (e.g., for gender-specific DE with age), or both contrasts at the same time (for any DE with age).

ADD COMMENT
0
Entering edit mode

Thanks a lot for for your answer and helping me!!!!

However, I tried to follow your instructions, but I don't know how to build it well :-S

age= pData(gse)[,"age"]
gender= pData(gse)[,"gender"]
disease= pData(gse)[,"disease"]
age.group= ifelse(age>40,"over40", "under40")

design = model.matrix(~0 + disease + gender + gender:age.group)
colnames(design) = c("diseasecontrol", "diseasecase","gendermale", "gendermale_age", "genderfemale_age")

fit.plus = lmFit(gse, design)

# 1. Compare by age:
contrast.matrix <- makeContrasts(genderfemale_age, gendermale_age, levels=design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
topTable(fit2, coef = 2, adjust = "BH")

(results= decideTests(fit2, adjust.method ="BH", p.value=0.005))
summary(results)

Error in contrasts.fit(fit, contrast.matrix) :
Number of rows of contrast matrix must match number of coefficients
More: Warning message:
In contrasts.fit(fit, contrast.matrix) :
row names of contrasts don't match col names of coefficients

# 2. Compare by gender:

contrast.matrix <- makeContrasts(gendermale, genderfemale_age - gendermale_age, levels=design)
> topTable(fit2, coef = NULL, adjust = "BH")

 

 

 

 

ADD REPLY
0
Entering edit mode

You're using a different design matrix from what you had in your original post, with a gender:age.group term rather than gender:age. This results in more coefficients in your design matrix, which probably causes the error and warning that you've observed. Why not just use the original design matrix? A threshold of 40 seems arbitrary anyway.

ADD REPLY
0
Entering edit mode

Again, many thanks for your useful help. I achieved results, and I understood and learned to do something else.

Problem solved!

ADD REPLY

Login before adding your answer.

Traffic: 443 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6