Question: Limma: Differential expression for micro RNA sequence data
0
8 months ago by
Perth
Michael Stone0 wrote:

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

Sample_id Age Sex Group Batch
mouse_1 6 F group1 batch1
mouse_2 10 M group1 batch2
mouse_3 25 F group2 batch1
mouse_4 36 M group2 batch2

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;

Symbol logFC AveExp t P.Value adj.P.Val B
RNA1 RNA1 0.76 8.6 4.6 1.4e-05 0.006 2.9
RNA2 RNA2 0.73 5.23 4.44 2.5e-05 0.006 2.3

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.

Michael

limma mirna-seq • 215 views
modified 8 months ago by Gordon Smyth37k • written 8 months ago by Michael Stone0

Did you modify the column names of your design matrix? The coefficient names in your call to makeContrasts don't look like what I would expect model.matrix to produce.

Answer: Limma: Differential expression for micro RNA sequence data
2
8 months ago by
Gordon Smyth37k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth37k wrote:

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.

Thank you Gordon

The first solution works straight out of the box.

Much appreciated.

Michael