[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