[BioC] LIMMA on RT-PCR data

Gordon K Smyth smyth at wehi.EDU.AU
Wed Nov 13 00:46:17 CET 2013


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:4}}



More information about the Bioconductor mailing list