Question

Limma; design matrix for treatments and sibship

0

Entering edit mode

John ▴ 10

@john-7115

Last seen 6.2 years ago

United Kingdom

Hello.

I have a targets file as follows;

treat sibship
file1 GeneInsert 1
file2 GeneInsert 2
file4 EmptyVector 1
file5 EmptyVector 2
file7 no_treatment 1
file8 no_treatment 2

Treat = treatment by tranfsection of an interesting gene, or with an empty vector transfected, and finally cells not treated at all. They are also paired by sibship.

How can I make differential gene expression of GeneInsert vs. no_treatment (but taking into account both the Empty vector and sibship)?

I have read the pertinent sections, 9.4 and 9.5 of the Limma manual, but I am not sure which is, if any, are applicable?

I will appreciate any help, thank you.

John.

Note; for an unpaired analysis I did;

contrast.matrix <- makeContrasts('InsertVscontrol'=GeneInsert-((no_treatment+EmptyVector)/2), levels=design)

Looks to be working as plots of expression values of the signif. diff. expressed genes look as they should (but I would like to incorporate the sibship as well to be precise).

limma design and contrast matrix • 1.5k views

ADD COMMENT • link 9.1 years ago John ▴ 10

score 1 · Answer 1 · 2015-10-20

1

Entering edit mode

James W. MacDonald 67k

@james-w-macdonald-5106

Last seen 1 day ago

United States

I don't know what you mean by 'taking into account...Empty vector'. That's a completely different treatment, that has nothing to do with the other two groups, so you can't really take it into account. You can however block on sibship to account for a shift in expression due to the different litters.

> treat <- factor(rep(c("GeneInsert","EmptyVector","NoTreatment"), each = 2))
> sibship <- factor(rep(1:2, each = 3))
> treat
[1] GeneInsert  GeneInsert  EmptyVector EmptyVector NoTreatment NoTreatment
Levels: EmptyVector GeneInsert NoTreatment
> treat <- relevel(treat, "NoTreatment")
> design <- model.matrix(~treat+sibship)
> design
  (Intercept) treatEmptyVector treatGeneInsert sibship2
1           1                0               1        0
2           1                0               1        0
3           1                1               0        0
4           1                1               0        1
5           1                0               0        1
6           1                0               0        1

The second coefficient compares the EmptyVector to NoTreatment, and the third compares GeneInsert to NoTreatment. If you want to compare GeneInsert and EmptyVector, you have to use a contrast.

ADD COMMENT • link 9.1 years ago James W. MacDonald 67k

1

Entering edit mode

James, I think your sibship definition might be wrong; it should be factor(rep(1:2, 3)).

Anyway, to elaborate on James' answer, the comparison between the inserted and empty vectors can be done with:

con <- makeContrasts(treatGeneInsert - treatEmptyVector, levels=design)

There's no point or need to consider the no-treatment group (besides for residual d.f. considerations, as Gavin has noted). The above comparison is equivalent to asking if the effect of the gene insert relative to no-treatment is different to the effect of the empty vector relative to no-treatment. This is because the no-treatment expression is represented by the intercept:

(treatGeneInsert - intercept) - (treatEmptyVector - intercept)

... which just cancels out when you expand past the parentheses.

ADD REPLY • link 9.1 years ago Aaron Lun ★ 28k

0

Entering edit mode

Thanks for the answer and code. I understand this is similar to section 9.4.1 of the Limma manual.

I thought the empty vector is related to the geneInsert samples, as that same vector is used in both sample groups; 1) to insert the gene and 2) the vector alone without the gene insert. Some genes will change due to the vector insertion in both sample groups.

Am I wrong in this thinking?

John.

ADD REPLY • link 9.1 years ago John ▴ 10

0

Entering edit mode

You could look at genes that were significant according to the third coefficient (treatGeneInsert) but not the second (treatEmptyVector). I would not recommended this, however, as you're going to decrease overall power due to the double round of testing. As the GeneInsert is effectively NoTreatment+InsertEffect+TreatmentEffect, and your EmptyVector is NoTreatment+InsertEffect, then you can estimate TreatmentEffect as GeneInsert - EmptyVector, so

treat <- relevel(treat, "EmptyVector")

and refitting will give a coefficient treatGeneInsert that estimates the Insert effect over and above the vector effect. The noTreatment group is not necessary to estimate effect size (and strictly needn't have been assayed), but the presence of these extra samples will increase power by providing better estimates of within-group variability. Sibship will have been dealt with correctly by the inclusion of it as a covariate.

ADD REPLY • link 9.1 years ago Gavin Kelly ▴ 690

0

Entering edit mode

Thanks Gavin for the clear logic in your answer and a solution, that helps a lot.

John.

ADD REPLY • link 9.1 years ago John ▴ 10

score 0 · Answer 2 · 2015-10-20

0

Entering edit mode

John ▴ 10

@john-7115

Last seen 6.2 years ago

United Kingdom

Added as a comment.

ADD COMMENT • link 9.1 years ago John ▴ 10