Hi Gordon and Ryan, Thank you so much for all your help. Relieved after knowing that I can use limma on this gene set, because for microarray I used limma. This is validation experiment with RT PCR. So I wanted to maintain uniformity (if that makes sense) My experimental design is spread across 85 samples like this- 4 batches- A, B, C and D In each Batch, 3 treatments - T1, T2, T3 Same patients were taken for treatments of a batch. I don't have comparisons across the batches, but I have comparisons within batches SampleName Batch Trt Patient S1 A 1 1 S2 A 2 1 S3 A 3 1 S4 A 1 2 S5 A 2 2 S6 A 3 2 .................... ..................... S7 D 1 30 S8 D 2 30 S9 D 3 30 Observation are missing on some samples. TD = factor(paste(Batch, Trt, sep="_")) Patient <- factor(Patient) design <- model.matrix(~0+TD+Patient) fit <- lmFit(dCT, design) cont.matrix <- makeContrasts(TDA_2 - TDA_1, levels = design) fit1 = contrasts.fit(fit, cont.matrix) fit2 = eBayes(fit1) Is this correct? I hope I don't need to use duplicateCorrelation as I don't want to compare Batch A with Batch B. Please advise. Thanks, Sandhya -----Original Message----- From: Gordon K Smyth [mailto:smyth@wehi.EDU.AU] Sent: Wednesday, November 13, 2013 5:16 AM To: Ryan Cc: Bioconductor mailing list; Sandhya Pemmasani Kiran Subject: Re: [BioC] LIMMA on RT-PCR data Hi Ryan, Here's a simulated example, just to illustrate: Simulate 50 patients with before and after treatment: > library(limma) > y <- matrix(rnorm(4*100),4,100) > patient <- gl(50,2) > treatment <- gl(2,1,100) Drop 15 observations as missing: > keep <- sample(1:100,85) > y <- y[,keep] > patient <- factor(patient[keep]) > treatment <- treatment[keep] > design <- model.matrix(~treatment) > dupcor <- duplicateCorrelation(y2,design,block=patient) > dupcor$consensus [1] 0.0166576 Theoretically the correlation should be zero. > fit <- lmFit(y,design,block=patient,correlation=dupcor$consensus) > fit <- eBayes(fit) > fit$df.total [1] 115.5603 115.5603 115.5603 115.5603 > fit$df.prior [1] 32.56032 > fit$df.residual [1] 83 83 83 83 The perfect solution here is to pool the variances which would give df.total = 4*83 = 332. The estimated value is slighty more conservative. limma takes care never to use df.total greater than the pooled value. > topTable(fit,coef=2) logFC AveExpr t P.Value adj.P.Val B 4 0.47806544 0.07233048 1.8941941 0.06069901 0.2427960 -3.954093 1 -0.12812909 0.04446857 -0.5713735 0.56885611 0.6946831 -5.112549 3 -0.09334881 0.08623914 -0.4705573 0.63884385 0.6946831 -5.150553 2 0.09463560 -0.00908128 0.3934894 0.69468313 0.6946831 -5.174667 No DE, which is correct. Cheers Gordon On Tue, 12 Nov 2013, Ryan wrote: > Wow, I wasn't aware that limma could be beneficial for such a small > number of genes. That's good to know. > > On Tue Nov 12 14:50:21 2013, Gordon K Smyth wrote: >> Hi Sandhya, >> >> Yes, limma should work fine on this data, although you are on the >> lower boundary in terms of number of genes. Theoretically, the >> minimum number of genes for the empirical Bayes procedure to be >> beneficial is three. Four genes is probably the minimum from a >> practical point of view. >> >> You may already know how to use duplicateCorrelation. If you have a >> simple before vs after paired comparison with some unmatched samples, >> then you could proceed like this: >> >> design <- model.matrix(~treatment) >> dupcor <- duplicateCorrelation(y, design, block=patient) fit <- >> lmFit(y, design, block=patient, correlation=dupcor$consensus) fit <- >> eBayes(fit) >> topTable(fit,coef=2) >> >> Be sure to check that dupcor$consensus is greater than zero. >> >> We used this strategy to compare tumour vs normal tissue in the >> presence of unmatched samples in >> >> http://www.ncbi.nlm.nih.gov/pubmed/17236199 >> >> although that was microarray data rather than PCR. >> >> Best wishes >> Gordon >> >> >> ---------- original message ---------- BioC] LIMMA on RT-PCR data >> Sandhya Pemmasani Kiran sandhya.p at ocimumbio.com Tue Nov 12 >> 12:16:25 CET 2013 >> >> Dear list, >> >> I have RT-PCR data on 4 genes and 85 samples. >> Can I use 'limma' on this small set of genes. >> >> I want to use limma rather than usual paired t-test because I have >> missing values and I don't want to miss the information available on >> paired samples.. >> >> Please advise. >> >> >> Thanks, >> Sandhya ______________________________________________________________________ The information in this email is confidential and intend...{{dropped:25}}