Help with complex limma design
1
0
Entering edit mode
@michael-imbeault-3593
Last seen 9.6 years ago
Hello, I want to analyze a recent Affymetrix Exon array dataset with the limma package; I imported values from APT and created and ExpressionSet with it, but the experiment's design is not the simplest and the user guide is not really helping - no real explanations of the 'how' to get there with the various presented use cases. It's 27 arrays, 3 time points, 3 sample sources, 3 treatments (see data frame below) - I'm interested in comparing 'mock' vs 'pos' and 'mock' vs 'neg' at all time points (24vs24, 48vs48, etc). The samples are paired; it's from the same blood donor, so i'd like to capture this effect in the data to improve p-values. I did a limma analysis using onechannelGUI (with great results) but it doesn't allow for complex designs (to capture the paired samples effect and improve p-values, increase power), so I have to do this manually. So far I tried design <- model.matrix(~Treatment*Time*Sample, data=targets) as I seen on the mailing list, where Treatment <- factor(targets$Treatment), etc. Data frame: Sample ID Treatment Time Sample 1 Mock 24h A 2 Neg 24h A 3 Pos 24h A 4 Mock 24h B 5 Neg 24h B 6 Pos 24h B 7 Mock 24h C 8 Neg 24h C 9 Pos 24h C 10 Mock 48h A 11 Neg 48h A 12 Pos 48h A 13 Mock 48h B 14 Neg 48h B 15 Pos 48h B 16 Mock 48h C 17 Neg 48h C 18 Pos 48h C 19 Mock 72h A 20 Neg 72h A 21 Pos 72h A 22 Mock 72h B 23 Neg 72h B 24 Pos 72h B 25 Mock 72h C 26 Neg 72h C 27 Pos 72h C Any help would be greatly appreciated, Michael [[alternative HTML version deleted]]
limma limma • 1.0k views
ADD COMMENT
0
Entering edit mode
@sunny-srivastava-3793
Last seen 9.6 years ago
Hi Michael, I am pretty sure the senior members of this group will have more to say on the design - but this is what I think is a solution to your problem (I may be wrong though :-( ) If I understand your design correctly: It looks like a split plot design with repeated measures. Mock, Neg and Pos are whole plot treatments with repeated measurements at times = 24h, 48h and 72h. The samples A, B, C are crossed across all treatments? If I am correct then I have the following question Do you think that there is no correlation between the measurements at the three different times?. If no, then I would try the something similar to section 8.2 or 11.4 (but in this case correlation would be between times) If yes, then this is what I would do (similar to the treatment in limmaUsersGuide() section 11.4 -there is just an added the extra block of sample - A,B,C) you can treat your design as a simple Randomized Block Design - with A, B and C as your blocks. The treatments are - (mock,24), (mock,48), (mock, 72) .. similarly 6 more treatments for the remaining neg and pos. In situations where I am interested in contrasts like yours, I prefer to a cell means model matrix rather than treatment effect model matrix. (but you can take your call - you just need to remember what are your coeffecients representing) ## I may be wrong here - but this is how I think the treatments are assigned in the d.f. treatments = factor(c(rep(c(1,4,7),3),rep(c(2,5,8),3),rep(c(3,6,9),3)),labels=c("m2 4","n24","p24","m48","n48","p48","m72","n72","p72")) ## construct the design matrix for the cell means anova model design <- model.matrix(~0+treatments) ## get the correlation with in a block ie a sample - this models the correlations structure in a particular block rho <- duplicateCorrelation(eset, design, block=pData(eset)$sample)$consensus fit1<- lmFit(eset, design, block=pData(eset)$sample, cor=rho) contrast.matrix <- makeContrasts(m24-n24, m24-p24, levels = design) ## constrats for the two comparsions that you mentioned fit2 <- contrasts.fit(fit1, contrast.matrix) fit3 <- eBayes(fit2) Hope this helps! Please correct me if I am wrong. This forum is a good way to re-learn experimental design and analysis. Thanks and Best regards, S. On Thu, Dec 17, 2009 at 2:38 AM, Michael Imbeault < michael.imbeault@sympatico.ca> wrote: > Hello, > > I want to analyze a recent Affymetrix Exon array dataset with the limma > package; I imported values from APT and created and ExpressionSet with > it, but the experiment's design is not the simplest and the user guide > is not really helping - no real explanations of the 'how' to get there > with the various presented use cases. > > It's 27 arrays, 3 time points, 3 sample sources, 3 treatments (see data > frame below) - I'm interested in comparing 'mock' vs 'pos' and 'mock' vs > 'neg' at all time points (24vs24, 48vs48, etc). The samples are paired; > it's from the same blood donor, so i'd like to capture this effect in > the data to improve p-values. I did a limma analysis using onechannelGUI > (with great results) but it doesn't allow for complex designs (to > capture the paired samples effect and improve p-values, increase power), > so I have to do this manually. > > So far I tried design <- model.matrix(~Treatment*Time*Sample, > data=targets) as I seen on the mailing list, where Treatment <- > factor(targets$Treatment), etc. > > Data frame: > > Sample ID Treatment Time Sample > 1 Mock 24h A > 2 Neg 24h A > 3 Pos 24h A > 4 Mock 24h B > 5 Neg 24h B > 6 Pos 24h B > 7 Mock 24h C > 8 Neg 24h C > 9 Pos 24h C > 10 Mock 48h A > 11 Neg 48h A > 12 Pos 48h A > 13 Mock 48h B > 14 Neg 48h B > 15 Pos 48h B > 16 Mock 48h C > 17 Neg 48h C > 18 Pos 48h C > 19 Mock 72h A > 20 Neg 72h A > 21 Pos 72h A > 22 Mock 72h B > 23 Neg 72h B > 24 Pos 72h B > 25 Mock 72h C > 26 Neg 72h C > 27 Pos 72h C > > > Any help would be greatly appreciated, > Michael > > [[alternative HTML version deleted]] > > _______________________________________________ > Bioconductor mailing list > Bioconductor@stat.math.ethz.ch > 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

Login before adding your answer.

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