Designing a model with blocking and other interactions
2
0
Entering edit mode
Eleanor Su ▴ 110
@eleanor-su-6460
Last seen 9.6 years ago
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 [[alternative HTML version deleted]]
Sequencing Sequencing • 1.6k views
ADD COMMENT
0
Entering edit mode
@gordon-smyth
Last seen 8 hours ago
WEHI, Melbourne, Australia
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}}
ADD COMMENT
0
Entering edit mode
Hi Gordon, Could you elaborate on how to test for the last coefficient? Sorry if that's an elementary question. I'm a master's student working with VERY limited resources in my department, and I appreciate all the help I can get. Thanks for all your help! Best, Eleanor On 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 inte...{{dropped:6}}
ADD REPLY
0
Entering edit mode
I mean when you get to the glmLRT() step of the analysis pipeline that you set coef=ncol(design). The design matrix has 7 columns, so lrt <- glmLRT(fit,coef=7) etc. One of your earlier emails said you had completed a differential expression analysis in edgeR, so I assume that you know what the above means. If you don't, then follow one of the case studies in the edgeR User's Guide. Gordon On Wed, 2 Apr 2014, Eleanor wrote: > Hi Gordon, > > Could you elaborate on how to test for the last coefficient? Sorry if > that's an elementary question. I'm a master's student working with VERY > limited resources in my department, and I appreciate all the help I can > get. Thanks for all your help! > > Best, > Eleanor > > On 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 intended solely for the addressee. >> You must not disclose, forward, print or use it without the permission of the sender. >> ______________________________________________________________________ > > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}
ADD REPLY
0
Entering edit mode
Hi Gordon, I understand now. Thank you so much for your help. Best, Eleanor On Wed, Apr 2, 2014 at 11:17 PM, Gordon K Smyth <smyth@wehi.edu.au> wrote: > I mean when you get to the glmLRT() step of the analysis pipeline that you > set coef=ncol(design). The design matrix has 7 columns, so > > lrt <- glmLRT(fit,coef=7) > > etc. > > One of your earlier emails said you had completed a differential > expression analysis in edgeR, so I assume that you know what the above > means. If you don't, then follow one of the case studies in the edgeR > User's Guide. > > Gordon > > > On Wed, 2 Apr 2014, Eleanor wrote: > > Hi Gordon, >> >> Could you elaborate on how to test for the last coefficient? Sorry if >> that's an elementary question. I'm a master's student working with VERY >> limited resources in my department, and I appreciate all the help I can >> get. Thanks for all your help! >> >> Best, >> Eleanor >> >> On Apr 2, 2014, at 8:25 PM, Gordon K Smyth <smyth@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@gmail.com> >>>> To: "bioconductor@stat.math.ethz.ch" <bioconductor@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 intended solely for >>> the addressee. >>> You must not disclose, forward, print or use it without the permission >>> of the sender. >>> ______________________________________________________________________ >>> >> >> >> > ______________________________________________________________________ > The information in this email is confidential and inte...{{dropped:10}}
ADD REPLY
0
Entering edit mode
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. 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. 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@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@gmail.com> >> To: "bioconductor@stat.math.ethz.ch" <bioconductor@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 inte...{{dropped:10}}
ADD REPLY
0
Entering edit mode
@gordon-smyth
Last seen 8 hours ago
WEHI, Melbourne, Australia
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}}
ADD COMMENT
0
Entering edit mode
Is this not a case for using duplicateCorrelation across the Family units, which would account for the correlation within a family, while still allowing for a mitoHap main effect? -Aaron On Thu, Apr 3, 2014 at 7:04 PM, Gordon K Smyth <smyth@wehi.edu.au> wrote: > 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@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@gmail.com> >>>> To: "bioconductor@stat.math.ethz.ch" <bioconductor@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 inte...{{dropped:13}}
ADD REPLY
0
Entering edit mode
It is, but that would require to a shift from edgeR to voom and limma, and one still can't test for a main effect in the presence of an interaction. I have refered Eleanor to a section in the limma User's Guide which explains the use of duplicateCorrelation for this purpose. Gordon On Thu, 3 Apr 2014, Aaron Mackey wrote: > Is this not a case for using duplicateCorrelation across the Family units, > which would account for the correlation within a family, while still > allowing for a mitoHap main effect? > > -Aaron > > > On Thu, Apr 3, 2014 at 7:04 PM, Gordon K Smyth <smyth at="" wehi.edu.au=""> wrote: > >> 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}} >> >> _______________________________________________ >> Bioconductor mailing list >> Bioconductor at r-project.org >> https://stat.ethz.ch/mailman/listinfo/bioconductor >> Search the archives: http://news.gmane.org/gmane. >> science.biology.informatics.conductor >> > ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:4}}
ADD REPLY

Login before adding your answer.

Traffic: 847 users visited in the last hour
Help About
FAQ
Access RSS
API
Stats

Use of this site constitutes acceptance of our User Agreement and Privacy Policy.

Powered by the version 2.3.6