Search
Question: Limma: Differential expression for micro RNA sequence data
0
gravatar for Michael Stone
3 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.

Thank you for your help.

 

Michael

ADD COMMENTlink modified 3 months ago by Gordon Smyth35k • written 3 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.

ADD REPLYlink written 3 months ago by Ryan C. Thompson7.1k
2
gravatar for Gordon Smyth
3 months ago by
Gordon Smyth35k
Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia
Gordon Smyth35k 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.

 

 

ADD COMMENTlink modified 3 months ago • written 3 months ago by Gordon Smyth35k

Thank you Gordon

The first solution works straight out of the box.

Much appreciated.

 

Michael

ADD REPLYlink written 3 months ago by Michael Stone0
Please log in to add an answer.

Help
Access

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.
Powered by Biostar version 2.2.0
Traffic: 195 users visited in the last hour