Limma design code with technical replicates, paired samples and two factors
1
0
Entering edit mode
@otto-benjamin-6552
Last seen 9.6 years ago
Dear list members, though I am aware of each of the single design models, I'm not quite sure if I got the design for the following setting correct: The experiment consists of 24 single-color samples in total obtained from 6 animals. From each animal two time points (prior to treatment and post treatment) exist. In addition each sample was obtained in a technical duplicate. 3 of the animals were treated with a reagent, the other three with a control. The matrix is provided below. Now what I would like to do is a classic calculation of time, treatment and time x treatment effects ... under consideration of 1) sample pairing and 2) technical replication. My first (simple) approach neglecting treatment is as follows: ----------------------------- Line 1 > biolrep <- pData$Replicate Line 2 > timep <- factor(pData$Timepoint, levels=c("Pre", "Post")) Line 3 > treat <- factor(pData$Treatment, levels=c("reagent","control")) Line 4 > sibship <- factor(pData$Animal) Line 5 > design <- model.matrix(~sibship+timep) # I will first skip treatment here for simplicity Line 6 > corfit<-duplicateCorrelation(eset,design, ndups=1, block=biolrep) # account for tech-reps. Line 7 > fit <- lmFit(eset, design, block=biolrep, cor=corfit$consensus) #include it in fit Line 8 > fit <- lmFit(eset, design) Line 9 > fit <- eBayes(fit) Line 10 > topTable(fit, adjust = "BH", coef="timepPost") -------------------------------- Once I incorporate the treatment instead of the time Line 5 > design <- model.matrix(~sibship+treat) ... I get an error: >> Partial NA coefficients for 96 probe(s) Note: Our array has 96 probes in total so actually all probe coefficients are NA. If I remove the animal source info and swap line 5 by: Line 5 > design <- model.matrix(~treat) ... then everything is fine. Probably because of the overlap of the groups (all samples of each animal are either part of the treatment or the non treatment group). So my question is: Is there some design model where I can incorporate time AND treatment under consideration of sample pairing and technical replication at the same time? Thanks for the help Benjamin Sample Animal Replicate Timepoint Treatment 1 Animal 1 1 Pre control 2 Animal 1 1 Pre control 3 Animal 1 2 Post control 4 Animal 1 2 Post control 5 Animal 2 3 Pre reagent 6 Animal 2 3 Pre reagent 7 Animal 2 4 Post reagent 8 Animal 2 4 Post reagent 9 Animal 3 5 Pre control 10 Animal 3 5 Pre control 11 Animal 3 6 Post control 12 Animal 3 6 Post control 13 Animal 4 7 Pre reagent 14 Animal 4 7 Pre reagent 15 Animal 4 8 Post reagent 16 Animal 4 8 Post reagent 17 Animal 5 9 Pre control 18 Animal 5 9 Pre control 19 Animal 5 10 Post control 20 Animal 5 10 Post control 21 Animal 6 11 Pre reagent 22 Animal 6 11 Pre reagent 23 Animal 6 12 Post reagent 24 Animal 6 12 Post reagent -- DANKE FÜR 125 JAHRE ENGAGEMENT UND VERTRAUEN. www.uke.de/125 _____________________________________________________________________ Besuchen Sie uns auf: www.uke.de _____________________________________________________________________ Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg Vorstandsmitglieder: Prof. Dr. Christian Gerloff (Vertreter des Vorsitzenden), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Prölß, Rainer Schoppik _____________________________________________________________________ SAVE PAPER - THINK BEFORE PRINTING [[alternative HTML version deleted]]
• 4.1k views
ADD COMMENT
0
Entering edit mode
@ryan-c-thompson-5618
Last seen 7 months ago
Scripps Research, La Jolla, CA
Hi Benjamin, Your fundamental issue is that your sibship and treatment variables are confounded, since each animal belongs to one and only one treatment. So, when you put both animal and treatment into the same design matrix, the result is not full rank. I believe that you must treat the animal factor not as a single factor with 6 levels, but rather as two separate factors with 3 levels each, within the two treatments. I'm not aware of a way to do this automatically, but you can accomplish it by manually constructing an appropriate contrasts matrix for the animal factor, which requires a bit of extra work. I've prepared some example code here: I am reasonably sure that the design matrix constructed above accurately represents your experimental design, but you may want to check with a statistician at your local institution to check my logic. Also, note that, when considering only biological replicates, you have 12 samples with a design matrix of rank 8, leaving only 4 residual degrees of freedom. Because of this, you might consider doing the analysis without using the animal factor at all, which would leave 8 residual df but lose the advantage of having a paired design for the timepoint test. -Ryan On 5/14/14, 7:47 AM, Otto, Benjamin wrote: > Dear list members, > > though I am aware of each of the single design models, I'm not quite sure if I got the design for the following setting correct: > > The experiment consists of 24 single-color samples in total obtained from 6 animals. From each animal two time points (prior to treatment and post treatment) exist. In addition each sample was obtained in a technical duplicate. 3 of the animals were treated with a reagent, the other three with a control. The matrix is provided below. > > Now what I would like to do is a classic calculation of time, treatment and time x treatment effects ... under consideration of 1) sample pairing and 2) technical replication. > > My first (simple) approach neglecting treatment is as follows: > ----------------------------- > Line 1 > biolrep <- pData$Replicate > Line 2 > timep <- factor(pData$Timepoint, levels=c("Pre", "Post")) > Line 3 > treat <- factor(pData$Treatment, levels=c("reagent","control")) > Line 4 > sibship <- factor(pData$Animal) > Line 5 > design <- model.matrix(~sibship+timep) # I will first skip treatment here for simplicity > > Line 6 > corfit<-duplicateCorrelation(eset,design, ndups=1, block=biolrep) # account for tech-reps. > Line 7 > fit <- lmFit(eset, design, block=biolrep, cor=corfit$consensus) #include it in fit > > Line 8 > fit <- lmFit(eset, design) > Line 9 > fit <- eBayes(fit) > Line 10 > topTable(fit, adjust = "BH", coef="timepPost") > -------------------------------- > > Once I incorporate the treatment instead of the time > > Line 5 > design <- model.matrix(~sibship+treat) > > ... I get an error: > >>> Partial NA coefficients for 96 probe(s) > Note: Our array has 96 probes in total so actually all probe coefficients are NA. > > > If I remove the animal source info and swap line 5 by: > > Line 5 > design <- model.matrix(~treat) > > ... then everything is fine. Probably because of the overlap of the groups (all samples of each animal are either part of the treatment or the non treatment group). > > > So my question is: Is there some design model where I can incorporate time AND treatment under consideration of sample pairing and technical replication at the same time? > > > Thanks for the help > > Benjamin > > > > > > > Sample Animal Replicate Timepoint Treatment > 1 Animal 1 1 Pre control > 2 Animal 1 1 Pre control > 3 Animal 1 2 Post control > 4 Animal 1 2 Post control > 5 Animal 2 3 Pre reagent > 6 Animal 2 3 Pre reagent > 7 Animal 2 4 Post reagent > 8 Animal 2 4 Post reagent > 9 Animal 3 5 Pre control > 10 Animal 3 5 Pre control > 11 Animal 3 6 Post control > 12 Animal 3 6 Post control > 13 Animal 4 7 Pre reagent > 14 Animal 4 7 Pre reagent > 15 Animal 4 8 Post reagent > 16 Animal 4 8 Post reagent > 17 Animal 5 9 Pre control > 18 Animal 5 9 Pre control > 19 Animal 5 10 Post control > 20 Animal 5 10 Post control > 21 Animal 6 11 Pre reagent > 22 Animal 6 11 Pre reagent > 23 Animal 6 12 Post reagent > 24 Animal 6 12 Post reagent > -- > > DANKE FÜR 125 JAHRE ENGAGEMENT UND VERTRAUEN. > www.uke.de/125 > _____________________________________________________________________ > > Besuchen Sie uns auf: www.uke.de > _____________________________________________________________________ > > Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg > Vorstandsmitglieder: Prof. Dr. Christian Gerloff (Vertreter des Vorsitzenden), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Prölß, Rainer Schoppik > _____________________________________________________________________ > > SAVE PAPER - THINK BEFORE PRINTING > > [[alternative HTML version deleted]] > > > > _______________________________________________ > Bioconductor mailing list > Bioconductor@r-project.org > https://stat.ethz.ch/mailman/listinfo/bioconductor > Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor [[alternative HTML version deleted]]
ADD COMMENT
0
Entering edit mode
Hi Ryan, thanks for super quick response and the marvelous elaborate code. regards Benjamin ________________________________ Von: Ryan [rct@thompsonclan.org] Gesendet: Donnerstag, 15. Mai 2014 00:49 An: Otto, Benjamin; bioconductor@r-project.org Betreff: Re: [BioC] Limma design code with technical replicates, paired samples and two factors Hi Benjamin, Your fundamental issue is that your sibship and treatment variables are confounded, since each animal belongs to one and only one treatment. So, when you put both animal and treatment into the same design matrix, the result is not full rank. I believe that you must treat the animal factor not as a single factor with 6 levels, but rather as two separate factors with 3 levels each, within the two treatments. I'm not aware of a way to do this automatically, but you can accomplish it by manually constructing an appropriate contrasts matrix for the animal factor, which requires a bit of extra work. I've prepared some example code here: I am reasonably sure that the design matrix constructed above accurately represents your experimental design, but you may want to check with a statistician at your local institution to check my logic. Also, note that, when considering only biological replicates, you have 12 samples with a design matrix of rank 8, leaving only 4 residual degrees of freedom. Because of this, you might consider doing the analysis without using the animal factor at all, which would leave 8 residual df but lose the advantage of having a paired design for the timepoint test. -Ryan On 5/14/14, 7:47 AM, Otto, Benjamin wrote: Dear list members, though I am aware of each of the single design models, I'm not quite sure if I got the design for the following setting correct: The experiment consists of 24 single-color samples in total obtained from 6 animals. From each animal two time points (prior to treatment and post treatment) exist. In addition each sample was obtained in a technical duplicate. 3 of the animals were treated with a reagent, the other three with a control. The matrix is provided below. Now what I would like to do is a classic calculation of time, treatment and time x treatment effects ... under consideration of 1) sample pairing and 2) technical replication. My first (simple) approach neglecting treatment is as follows: ----------------------------- Line 1 > biolrep <- pData$Replicate Line 2 > timep <- factor(pData$Timepoint, levels=c("Pre", "Post")) Line 3 > treat <- factor(pData$Treatment, levels=c("reagent","control")) Line 4 > sibship <- factor(pData$Animal) Line 5 > design <- model.matrix(~sibship+timep) # I will first skip treatment here for simplicity Line 6 > corfit<-duplicateCorrelation(eset,design, ndups=1, block=biolrep) # account for tech-reps. Line 7 > fit <- lmFit(eset, design, block=biolrep, cor=corfit$consensus) #include it in fit Line 8 > fit <- lmFit(eset, design) Line 9 > fit <- eBayes(fit) Line 10 > topTable(fit, adjust = "BH", coef="timepPost") -------------------------------- Once I incorporate the treatment instead of the time Line 5 > design <- model.matrix(~sibship+treat) ... I get an error: Partial NA coefficients for 96 probe(s) Note: Our array has 96 probes in total so actually all probe coefficients are NA. If I remove the animal source info and swap line 5 by: Line 5 > design <- model.matrix(~treat) ... then everything is fine. Probably because of the overlap of the groups (all samples of each animal are either part of the treatment or the non treatment group). So my question is: Is there some design model where I can incorporate time AND treatment under consideration of sample pairing and technical replication at the same time? Thanks for the help Benjamin Sample Animal Replicate Timepoint Treatment 1 Animal 1 1 Pre control 2 Animal 1 1 Pre control 3 Animal 1 2 Post control 4 Animal 1 2 Post control 5 Animal 2 3 Pre reagent 6 Animal 2 3 Pre reagent 7 Animal 2 4 Post reagent 8 Animal 2 4 Post reagent 9 Animal 3 5 Pre control 10 Animal 3 5 Pre control 11 Animal 3 6 Post control 12 Animal 3 6 Post control 13 Animal 4 7 Pre reagent 14 Animal 4 7 Pre reagent 15 Animal 4 8 Post reagent 16 Animal 4 8 Post reagent 17 Animal 5 9 Pre control 18 Animal 5 9 Pre control 19 Animal 5 10 Post control 20 Animal 5 10 Post control 21 Animal 6 11 Pre reagent 22 Animal 6 11 Pre reagent 23 Animal 6 12 Post reagent 24 Animal 6 12 Post reagent -- DANKE FÜR 125 JAHRE ENGAGEMENT UND VERTRAUEN. www.uke.de/125<http: www.uke.de="" 125=""> _____________________________________________________________________ Besuchen Sie uns auf: www.uke.de<http: www.uke.de=""> _____________________________________________________________________ Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg Vorstandsmitglieder: Prof. Dr. Christian Gerloff (Vertreter des Vorsitzenden), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Prölß, Rainer Schoppik _____________________________________________________________________ SAVE PAPER - THINK BEFORE PRINTING [[alternative HTML version deleted]] _______________________________________________ Bioconductor mailing list Bioconductor@r-project.org<mailto:bioconductor@r-project.org> https://stat.ethz.ch/mailman/listinfo/bioconductor Search the archives: http://news.gmane.org/gmane.science.biology.informatics.conductor -- DANKE FÜR 125 JAHRE ENGAGEMENT UND VERTRAUEN. www.uke.de/125 _____________________________________________________________________ Besuchen Sie uns auf: www.uke.de _____________________________________________________________________ Universitätsklinikum Hamburg-Eppendorf; Körperschaft des öffentlichen Rechts; Gerichtsstand: Hamburg Vorstandsmitglieder: Prof. Dr. Christian Gerloff (Vertreter des Vorsitzenden), Prof. Dr. Dr. Uwe Koch-Gromus, Joachim Prölß, Rainer Schoppik _____________________________________________________________________ SAVE PAPER - THINK BEFORE PRINTING [[alternative HTML version deleted]]
ADD REPLY

Login before adding your answer.

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