# [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..
>>