Hi all, 

I have run a differential expression analysis using limma looking for DE between group1 and group2, and male versus female

I have RNAseq data in a single tissue type for 80 individuals, with age groups being group1 = 1:20wks old and group2 being 21:40 weeks old. Groups 1 and 2 are a mix of male and female individuals. There are two seperate sequencing runs with associated batch effects. No individual is sampled more than once eg

Up to now I have used a design and contrast matrix as follows;

design <- model.matrix(~ Group + Sex + Batch)
cont.mat_Age <- c(0,1,0,0) # when looking for DE between age groups
cont.mat_Sex <- c(0,0,1,0) # when looking for DE between sexes

followed by voomWithQualityWeights(), lmfit(), contrasts.fit() and eBays().

What I would like to do is to test for DE between age group 1 and age group 2, treating male and female as seperate groups, but as I understand it I shouldn't just subset the counts table and run seperate regression, but rather control for sex in the model. I found a post (Limma: Continuous variable model designs) from 3 yrs ago and adapted the code to;

design_independent <- model.matrix(~0+Sex+Sex:Group + Batch)

contrast.independent <- makeContrasts(Sex_F_Group2 - Sex_M_Group2, levels=design_independent)

followed by voomWithQualityWeights(), lmfit(), contrasts.fit() and eBays().

with a (mock) topTable as follows;

My code runs, but I'm unsure if I'm getting what I'm asking for.

Questions:

1. Is the design and contrast correct?

2. If I want to make sure the batch effect is accounted for, how do I include it in the contrast.independent?

3. I'm unsure how to interpret the results: Does this give me the difference for females only between age group 1 and age group 2.

Thank you for your help.

Michael

No, the contrast you are making isn't correct. The correct code is somewhat simpler.

You can either use a factorial model formula, in which case you don't need to make contrasts at all, or you can treat your experiment as a oneway layout and make the obvious contrast. Either way, there's no need for any trickery.

Option 1: Factorial formula and no contrast

A nested factorial model allows you to compare groups within each sex. This achieves your aim of treating the two genders separately:

design <- model.matrix(~ Sex + Sex:Group + Batch)
fit <- lmFit(y, design)
fit <- eBayes(fit)
topTable(fit, coef=3) # group 2 vs group 1 for females
topTable(fit, coef=4) # group 2 vs group 1 for males

Option 2: One big factor and contrasts 

Alternatively you can combine your two factors of interest into one combined factor and make the contrasts you want explicitly:

SexGroup <- factor(paste(Sex,Group,sep="."))
design <- model.matrix(~0 + SexGroup + Batch)
colnames(design) <- c(levels(SexGroup), "Batch")
cont.mat <- makeContrasts(
                Group.Female = F.group2 - F.group1,
                Group.Male =   M.group2 - M.group1,
            levels=design)
fit <- lmFit(y, design)
fit2 <- contrasts.fit(fit, cont.mat)
fit2 <- eBayes(fit2)
topTable(fit2, coef="Group.Female")
topTable(fit2, coef="Group.Male")

I recommend the first approach, although the two are in principle equivalent.