technical + biological replicate design for multiple categories
2
0
Entering edit mode
James Platt ▴ 20
@james-platt-6025
Last seen 9.6 years ago
I am using Limma for an experimental design with three categories, two biological replicates each and two technical replicates of each biological replicate for a total of 12 arrays. Since this goes just a little bit beyond the Limma user guide, I want to make sure I am setting this up correctly. Also, figure out why I am getting the error message below. The categories are: Cyt1 Cyt2 Mock (which is control for both of the other groups) The biological replicates are indicated with A and B. The technical replicates are indicated with .1 and .2. The variable expr has the array data and the sample columns are: MockA.1 MockA.2 Cyt1A.1 Cyt1A.2 Cyt2A.1 Cyt2A.2 MockB.1 MockB.2 Cyt1B.1 Cyt1B.2 Cyt2B.1 Cyt2B.2 > treatment [1] MockA MockA Cyt1A Cyt1A Cyt2A Cyt2A MockB MockB Cyt1B Cyt1B Cyt2B Cyt2B Levels: Cyt1A Cyt1B Cyt2A Cyt2B MockA MockB > design = model.matrix(~ 0 + treatment) > design Cyt1A Cyt1B Cyt2A Cyt2B MockA MockB 1 0 0 0 0 1 0 2 0 0 0 0 1 0 3 1 0 0 0 0 0 4 1 0 0 0 0 0 5 0 0 1 0 0 0 6 0 0 1 0 0 0 7 0 0 0 0 0 1 8 0 0 0 0 0 1 9 0 1 0 0 0 0 10 0 1 0 0 0 0 11 0 0 0 1 0 0 12 0 0 0 1 0 0 attr(,"assign") [1] 1 1 1 1 1 1 attr(,"contrasts") attr(,"contrasts")$treatment [1] "contr.treatment" > block [1] "MockA" "MockA" "Cyt1A" "Cyt1A" "Cyt2A" "Cyt2A" "MockB" "MockB" "Cyt1B" "Cyt1B" "Cyt2B" "Cyt2B" > dupcor = duplicateCorrelation(expr, design, block=block) > fit.raw = lmFit(expr, design, block=block, + correlation=dupcor$consensus) Error in chol.default(V) : the leading minor of order 2 is not positive definite [[alternative HTML version deleted]]
limma limma • 1.1k views
ADD COMMENT
0
Entering edit mode
Naomi Altman ★ 6.0k
@naomi-altman-380
Last seen 2.9 years ago
United States
You need to have a variable "biorep" with biorep number e.g. c(1,1,2,2,3,3,...6,6) This is a blocking variable. The treatments are just Cyt1, Cyt2 and Mock. You will need to follow the instructions for blocking variables. --Naomi At 05:15 PM 7/2/2013, James Platt wrote: >I am using Limma for an experimental design with three categories, >two biological replicates each and two technical replicates of each >biological replicate for a total of 12 arrays. Since this goes >just a little bit beyond the Limma user guide, I want to make sure I >am setting this up correctly. Also, figure out why I am getting the >error message below. > >The categories are: >Cyt1 >Cyt2 >Mock (which is control for both of the other groups) > >The biological replicates are indicated with A and B. The technical >replicates are indicated with .1 and .2. > >The variable expr has the array data and the sample columns are: >MockA.1 MockA.2 Cyt1A.1 Cyt1A.2 Cyt2A.1 Cyt2A.2 MockB.1 MockB.2 >Cyt1B.1 Cyt1B.2 Cyt2B.1 Cyt2B.2 > > > treatment > [1] MockA MockA Cyt1A Cyt1A Cyt2A Cyt2A MockB MockB Cyt1B Cyt1B Cyt2B Cyt2B >Levels: Cyt1A Cyt1B Cyt2A Cyt2B MockA MockB > > design = model.matrix(~ 0 + treatment) > > design > Cyt1A Cyt1B Cyt2A Cyt2B MockA MockB >1 0 0 0 0 1 0 >2 0 0 0 0 1 0 >3 1 0 0 0 0 0 >4 1 0 0 0 0 0 >5 0 0 1 0 0 0 >6 0 0 1 0 0 0 >7 0 0 0 0 0 1 >8 0 0 0 0 0 1 >9 0 1 0 0 0 0 >10 0 1 0 0 0 0 >11 0 0 0 1 0 0 >12 0 0 0 1 0 0 >attr(,"assign") >[1] 1 1 1 1 1 1 >attr(,"contrasts") >attr(,"contrasts")$treatment >[1] "contr.treatment" > > block > [1] "MockA" "MockA" "Cyt1A" "Cyt1A" "Cyt2A" "Cyt2A" "MockB" > "MockB" "Cyt1B" "Cyt1B" "Cyt2B" "Cyt2B" > > dupcor = duplicateCorrelation(expr, design, block=block) > > fit.raw = lmFit(expr, design, block=block, >+ correlation=dupcor$consensus) >Error in chol.default(V) : > the leading minor of order 2 is not positive definite > [[alternative HTML version deleted]] > >_______________________________________________ >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
ADD COMMENT
0
Entering edit mode
James Platt ▴ 20
@james-platt-6025
Last seen 9.6 years ago
On Jul 3, 2013, at 1:32 PM, Naomi Altman wrote: > You need to have a variable "biorep" with biorep number e.g. c(1,1,2,2,3,3,...6,6) > > This is a blocking variable. So I have to use numbers for the blocking variable? I can't use text? Otherwise, it's the same as the variable which I called "block." > The treatments are just Cyt1, Cyt2 and Mock. My design matrix should have just those three columns? Then, my variable "treatment" should actually be the following: treatment <- c("Mock", "Mock", "Cyt1", "Cyt1", "Cyt2", "Cyt2", "Mock", "Mock", "Cyt1", "Cyt1", "Cyt2", "Cyt2") > You will need to follow the instructions for blocking variables. If I understand the user guide correctly then it's the same as what I did before except for having the corrected "biorep" and "treatment" variables: design = model.matrix(~ 0 + treatment) dupcor = duplicateCorrelation(expr, design, block=biorep) fit.raw = lmFit(expr, design, block=biorep, correlation=dupcor$consensus) fit = eBayes(fit.raw) This runs without error but I just want to make sure that it's correct. The results do make sense. In other words, the top ranked genes include many that were expected. Thanks, James
ADD COMMENT
0
Entering edit mode
The blocking variable can have text. Your treatment variable is now correct and as far as I can see, your analysis should now be fine. --Naomi At 05:34 PM 7/3/2013, you wrote: >On Jul 3, 2013, at 1:32 PM, Naomi Altman wrote: > > > You need to have a variable "biorep" with biorep number e.g. > c(1,1,2,2,3,3,...6,6) > > > > This is a blocking variable. > >So I have to use numbers for the blocking variable? I can't use >text? Otherwise, it's the same as the variable which I called "block." > > > The treatments are just Cyt1, Cyt2 and Mock. > >My design matrix should have just those three columns? Then, my >variable "treatment" should actually be the following: > >treatment <- c("Mock", "Mock", "Cyt1", "Cyt1", "Cyt2", "Cyt2", >"Mock", "Mock", "Cyt1", "Cyt1", "Cyt2", "Cyt2") > > > > You will need to follow the instructions for blocking variables. > >If I understand the user guide correctly then it's the same as what >I did before except for having the corrected "biorep" and >"treatment" variables: > >design = model.matrix(~ 0 + treatment) >dupcor = duplicateCorrelation(expr, design, block=biorep) >fit.raw = lmFit(expr, design, block=biorep, correlation=dupcor$consensus) >fit = eBayes(fit.raw) > >This runs without error but I just want to make sure that it's >correct. The results do make sense. In other words, the top ranked >genes include many that were expected. > >Thanks, >James >_______________________________________________ >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
ADD REPLY

Login before adding your answer.

Traffic: 736 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