[BioC] Designing a model with blocking and other interactions

[Gordon K Smyth]

On Thu, 3 Apr 2014, Eleanor Su wrote:

> Hi Gordon,
>
> When I enter the design that you've suggested,
>
> design1 <- model.matrix(~Family)
>  design2 <- model.matrix(~mitoHap*Treatment)
>  design <- cbind(design1,design2[,3:4])
>
> and test for the last coefficient, I see that I get DE for the interaction
> between Treatment:mitoHap (which is what I wanted to look at). As I look
> through the other columns in the design matrix, I see that I have data for
> Treatment (coef=6) but not for mitoHap.

No, naturally you can't have a mitoHap column because that factor is 
confounded with Family.

> If I use an equivalent formula for design2
>
> design2<-model.matrix(~mitoHap+Treatment+mitoHap:Treatment)
>
> would this allow me to see both factors (treatment and mitoHap)
> independently in other columns of the design matrix AND the interaction
> between the two in the last coefficient? I'd like to be able to look at
> differential expression in each factor independently and the interaction
> between the two.

Well, now you are ignoring Family, which previously you felt it was 
important to account for.

You are also asking for things that are impossible.  It isn't meaningful 
to test for factors independently of their interaction.

You might find it helpful to read the sections on "multi level designs" in 
the edgeR and limma User's Guides.

Gordon

> If so, how would the last "design" formula change?
>
> design<-cbind(design1,design2[?])
>
> Thanks for you help.
>
> Best,
> Eleanor
>
>
> On Wed, Apr 2, 2014 at 8:25 PM, Gordon K Smyth <smyth at wehi.edu.au> wrote:
>
>> Dear Eleanor,
>>
>>   design1 <- model.matrix(~Family)
>>   design2 <- model.matrix(~mitoHap*Treatment)
>>   design <- cbind(design1,design2[,3:4])
>>
>> Then test for the last coefficient.
>>
>> Best wishes
>> Gordon
>>
>>  Date: Tue, 1 Apr 2014 11:24:52 -0700
>>> From: Eleanor Su <eleanorjinsu at gmail.com>
>>> To: "bioconductor at stat.math.ethz.ch" <bioconductor at stat.math.ethz.ch>
>>> Subject: [BioC] Designing a model with blocking and other interactions
>>>
>>> Hi All,
>>>
>>> I'm trying to set up a model matrix where I can look at the interaction
>>> between Treatment and mitochondrial haplotypes in my paired samples. These
>>> are the preliminary commands that I've set up:
>>>
>>>  rawdata<-read.delim("piRNAtotalcount<10.txt", check.names=FALSE,
>>>>
>>> stringsAsFactors=FALSE)
>>>
>>>> y <- DGEList(counts=rawdata[,2:11], genes=rawdata[,1])
>>>> Family<-factor(c(6,6,9,9,11,11,26,26,28,28))
>>>> Treatment<-factor(c("C","H","C","H","C","H","C","H","C","H"))
>>>> mitoHap<-factor(c("S","S","S","S","S","S","D","D","D","D"))
>>>> data.frame(Sample=colnames(y),Family,Treatment,mitoHap)
>>>>
>>>    Sample Family Treatment mitoHap
>>> 1   6C (S)      6         C       S
>>> 2   6H (S)      6         H       S
>>> 3   9C (S)      9         C       S
>>> 4   9H (S)      9         H       S
>>> 5  11C (S)     11         C       S
>>> 6  11H (S)     11         H       S
>>> 7  26C (D)     26         C       D
>>> 8  26H (D)     26         H       D
>>> 9  28C (D)     28         C       D
>>> 10 28H (D)     28         H       D
>>>
>>>  design<-model.matrix(?)
>>>>
>>>
>>> I have 10 sequencing samples from 5 different families (a treatment and
>>> control sample from each family) and two different types of mitochondrial
>>> haplotypes. How do I set up a design where I can look at the interaction
>>> between the Treatments and mitoHap while still accounting for Family?
>>>
>>> Any help would be greatly appreciated. Thank you for your time.
>>>
>>> Best,
>>> Eleanor

______________________________________________________________________
The information in this email is confidential and intend...{{dropped:4}}